Non-Monotone Projected Gradient Method in Linear Elasticity Contact Problems with Given Friction

: We are focusing on the algorithms for solving the large-scale convex optimization problem in linear elasticity contact problems discretized by Finite Element method (FEM). The unknowns of the problem are the displacements of the FEM nodes, the corresponding objective function is deﬁned as a convex quadratic function with symmetric positive deﬁnite stiffness matrix and additional non-linear term representing the friction in contact. The feasible set constraints the displacement subject to non-penetration conditions. The dual formulation of this optimization problem is well-known as a Quadratic Programming (QP) problem and can be considered as a most basic non-linear optimization problem. Understanding these problems and the development of efﬁcient algorithms for solving them play the crucial role in the large-scale problems in practical applications. We shortly review the theory and examine the behavior and the efﬁciency of Spectral Projected Gradient method modiﬁed for QP problems (SPG-QP) on the solution of a toy example in MATLAB environment.


Introduction
One of the key components of the structural design is the contact problem. This problem is crucial, for example, during the design process of steel components, namely welded and screwed contacts. From the mathematical point of view, we are dealing with variational equalities and variational inequalities [1]. Since the problem is analytically solvable only for simple geometries, the typical approach is to discretize the problem and using the Finite Element Method (FEM, [2,3]). Using the reformulation of the problem by the so-called Ritz-Galerkin method and under certain assumptions, one can mathematically prove that the solution of the discretized problem is the approximation of the contiguous one. Additionally, the error converges asymptotically to zero with the density increase in the used mesh [1,3]. The advantage of the method is the possibility of the formulation of the problem as the (constrained) energy functional minimization problem. In the case of the linear elasticity, this corresponding problem is the Quadratic programming (QP) problem [4].
QP problems belong among to the most important optimization problems in practical applications; for example, in the case of the least square methods for regression analysis with the Euclidean distance [5], the solution of the problems with symmetric positive definite system matrix [6,7], Support Vector Machine (SVM, [8,9]), the inner subproblem of more complex (non-linear) optimization problems [10,11], and the numerical solution of variational inequalities in the linear elasticity contact problems [1,4,12,13].
In our case, the Hessian matrix represents the stiffness matrix and the linear term of the objective function is the force density vector. The unknown is the vector of the displacements in the FEM nodes. In the case of contact problems, the non-penetration conditions of the bodies and obstacles are defining the feasible set of the corresponding optimization problem with the linear inequality constraints. Additionally, if we consider also a friction between the bodies and/or the obstacles, the original contiguous problem captures this property by the new term influencing the solution displacement. For example, in the case of the linear elasticity contact problem with the Tresca (given) friction, the discretized minimization problem of the energy functional with the additional non-linear term representing the response of displacement caused by friction and constraints representing the non-penetration, can be reformulated using the dualization into the QP problem with bound and separable quadratic inequality constraints [1,12,14].
Although the QP problems belong to the basic non-linear optimization problems necessary for understanding non-linear optimization theory in general, there does not exist the best algorithm for solving problems from different applications, especially when one introduces the Dirichlet boundary conditions by equality constraints and/or uses non-overlapping domain decomposition methods, such as the Finite Element Tearing and Interconnecting Method (FETI, [4,[15][16][17][18][19][20][21]). The contact problems of linear elasticity are characterized by the bounded spectrum independent of the used mesh [22]. The utilization of this property has been proven and practically applied/demonstrated only in the case of active-set algorithms, namely the MPRGP algorithm (Modified Proportioning with Reduced Gradient Projections, [4,8,12,23,24]). Another approach is to use Interior-Point (IP, [25]) method-we eliminate the inequality constraints using barrier functions and solve the sequence of resulting saddle-point problems. In this paper, we examine the projected gradient descent method [26][27][28], especially SPG-QP (Spectral Projected Gradients [29] for QP [5]).

Problem Definition
As a simple benchmark, we consider the block of homogeneous material Ω with fixed zero displacement on boundary Γ D and imposed traction forces F on Γ F . The part Γ C denotes the part of boundary that may get into contact with a rigid obstacle. The block is attracted to the obstacle as consequence of the gravity force F g , see Figure 1. We solve the discretized form of the problem using FEM (see [2,3,30]), which leads to the optimization problem u * := arg min where N ∈ N is a number of used FEM nodes, n = 3N is a number of primal variables, m c ≤ N is a number of FEM nodes in Γ C , the terms of the objective function are defined as and • u ∈ R n is a vector of unknown displacements in FEM nodes; • Ω C ⊂ R n is a set of feasible (non-penetrated) displacements; • K ∈ R n,n is a symmetric positive definite (SPD) stiffness matrix (the Dirichlet boundary condition is implemented by the modification of the stiffness matrix); • f ∈ R n is a vector of forces density resulting from the stresses imposed on the structure during a displacement (the discretized form of F and F g ); • j h : R n → R is a numerical integration of functional describing the friction forces in the weak formulation of the problem; . . , m c are basis vectors formed by appropriately placed multiples of the unit tangential vectors in such a way that the jump of tangential displacement in i-th FEM node is given by T i u; . . , m are slip bound coefficients associated with T i .
Since our problem has a simple geometry, see Figure   For every contact node (i-th node of Γ C ) T i ∈ R 2,n is given by a zero matrix with 1 in the first row on appropriate x-coordinate of i-th node and in the second row on appropriate y-coordinate of i-th node. Consequently, the matrix composed as is a full rank matrix. We can modify the non-differentiable term j h in (1) into equivalent form (see [1]) where τ i ∈ R 2 are regulation variables. We denote the function and the vector and we write constraints τ i 2 ≤ ψ i in (2) in the form After the substitution into (1), we get where we denoted Λ τ ⊂ R 2m c as a feasible set defined by constraints (4). We consider L(u, τ) as a Lagrange function and τ as a vector of the Lagrange multipliers (in notation (3)) and we use the duality theorem (see Dostál [4]) to reformulate problem (5) into We include condition u ∈ Ω C by creating new Lagrange multipliers where matrix B ∈ R m c ,n and vector c ∈ R m c are constructed in such way that Due to geometry in our problem we can construct B simply; B is zero matrix with −1 in every i-th row (which is corresponding to i-th node in Γ C ) on appropriate z-coordinate of i-th node (see the choice of normal vectors for nodes in Γ C in Figure 2). Using these values, the set (8) consists of the feasible displacements which do not penetrate the rigid obstacle.
Problem (6) is equivalent to the saddle point problem where λ := τ λ C ,B := T B ,ĉ := 0 c We derive the dual problem corresponding to (9). From the first KKT condition (since stiffness matrix is SPD, we can use standard inversion K −1 ) we get and substitute into Lagrange function (7). After the simplifications, we get We obtain the function of only one variable λ.
Since we want to find maximizer (see saddle-point problem (9)), we omit the constant term and change signs and λ * solves equivalent minimization problem where we denoted After solving the minimization problem (12), the corresponding solution u * of primal problem (1) can be evaluated using (11).
Obviously, A ∈ R 3m c ,3m c is SPD and problem (12) is the QP problem with separable quadratic constraints combinated with bound constraints. Such a problem always has a unique solution [4].

Numerical Solution
To solve the problem (12), we use the Spectral Projected Gradient method (SPG, [29]), which is a projection gradient method for solving the general minimization problems The method is based on the Barzilai-Borwein (BB, [31,32]) step-size where ∇ f (x) ∈ R n is a gradient of objective function f : R n → R and P Θ : R n → Θ is the projection onto closed convex feasible set Θ ⊂ R n . However, since this procedure does not necessarily provide the decreasing sequence of function values (see [33]), an additional step is introduced where d it ∈ R n is called spectral projected gradient. The appropriate step-size β it ∈ (0, ∞) is computed using the iterative Grippo-Lampariello-Lucidi line-search method (GLL, [34]) to satisfy so-called generalized Armijo condition Here, τ ∈ (0, 1) represents a safeguarding parameter and with predefined constant m ∈ N. Because of the enforcement of condition (14), the algorithm generates the sequence of approximations with a decrease in objective function in every m iterations and the algorithm is converging to the minimizer. Recently, [5] introduced the simplification of the method in the case of the quadratic optimization problem. If the objective function is quadratic, one can derive the analytical formula for the computation of step-size, which satisfies (14) β it ∈ σ 1 , min{σ 2 ,β} with safeguarding parameters 0 < σ 1 < σ 2 < 1 and where A = ∇ 2 f (x it ) ∈ R n,n is a Hessian matrix of quadratic objective function f : R n → R. The algorithm can be written with only single matrix-vector multiplication during one iteration, see Algorithm 1.
Return the approximation of solution x it and the number of performed iterations it.
We used the presented Spectral Projected Gradient method for QP for solving the dual problem (12) with variable λ ∈ R m c , objective function Θ : R 3m c → R, and feasible set (10) defined by separable quadratic constraints and bound constraints.

Results
We implement and solve the problem introduced in Section 2 using MATLAB. For the FEM discretization, we adopt the open-source library presented in [35]. Instead of evaluating and computing the inverse of the matrix in the dualization (13), we solve the system of linear equations with multiple right-hand side vectors. In MATLAB, we evaluate In our numerical experiment, we consider steel brick with E = 2 × 10 5 [MPa], µ = 0.33, ρ = 7.85 × 10 −2 [kg · m −3 ] of size 2 × 1 × 0. 25 [m] in contact with the rigid obstacle. The displacement and the friction are caused by the force F = [450, 0, 0] [MPa]. We consider given (Tresca) friction between the brick and the obstacle. To discretize the problem, we construct regular grid with the number of FEM nodes in axis equal to N x = 25, N y = 13, N z = 4. In total, we obtain N = 1300 FEM nodes, the primal problem of dimension n = 3900, the number of FEM nodes in contact m c = 325, and dual problem of 975 unknowns. The feasible set of dual problem is composed from 325 bound constraints and 325 separable quadratic constraints.
We implemented SPG-QP in MATLAB and set the parameters of the Armijo criterion to recommended values m = 10, τ = 0.9 and initial approximation equal to the (feasible) zero vector. The first BB step-size is equal to steepest descent (Cauchy) step-size. We solved the problem (see Figures 3 and 4) with the coefficient of friction (slip bound) ψ i = 100S i [MPa], where S i is the contact surface corresponding to i-th node. Figure 5 demonstrates the decrease in the Euclidean norm of scaled projected gradient whereᾱ ∈ 0, 1/λ A max (see [36]) and λ A max is an upper estimation of the largest eigenvalue of the Hessian matrix A, which we compute using Gershgorin circle theorem. The norm of the projected gradient is used as a stopping criterion of the algorithm-i.e., we stop the algorithm if In the SPG-QP algorithm, we decided to choose the largest possible β it to perform the steps in as large a manner as possible. This approach has been suggested by Fletcher [33].   during the solution process. The notation x * is used for the solution obtained in the last iteration. Figure 6 shows the dependency of the number of iterations on the value of friction coefficient (the size of quadratic constraints) and the size of the problem. In the first case, we solve the problem of dimension N = 1300 and absolute error ε = 10 −6 for the norm of scaled projected gradient. In the second case, we are solving the problem with friction coefficientψ = 100 [MPa] while varying the problem size using the same stopping criterion. To improve the performance of the solution process, we utilize the GPU capability of the MATLAB. In our case, we compute on GeForce RTX 2080 Ti (4352 cuda cores, 11GB GDDR6 memory) and compare with Intel Core i9-7900X in MATLAB R2017a. Since our code is highly vectorized, the only prerequisite to perform the computation on GPU is to transfer the data using MATLAB command gpuArray. The overloaded MATLAB operations perform the matrix and vector operations on GPU. To analyze the performance of this approach, we fix the number of iterations of the SPG-QP algorithm to be equal to 10 4 and measure the computational time. From the results presented in Figure 7, it is clear that the main computational bottleneck is not the solution of dual QP problem, but the dualization-i.e., the solution of the system of linear Equations (15). The results suggest the suitability of used gradient projected method on GPU; however, the problem has to be sufficiently large to efficiently use this architecture. The performance of SPG-QP is limited by the matrix-vector multiplication with dense dual Hessian matrix, which is the most time-consuming operation. Consequently, the time complexity of the method scales quadratically with the dimension of dual problem, see Figure 7.

Conclusions
We examined SPG-QP algorithm for solving the strictly convex QPQC in the solution of linear elasticity contact problem with given friction. The results demonstrate the suitability of the algorithm for this type of problems. However, with the increasing size of the primal problem, we observe the increase in the computational complexity of the corresponding Hessian matrix inverse (which has polynomial complexity, in general). We suggested using domain decomposition methods, for instance the Finite Element Tearing and Interconnecting method (FETI). This extension will be the topic of our upcoming research.