An Accelerated Uzawa Method for Application to Frictionless Contact Problem

TL;DR

This paper introduces an accelerated version of the Uzawa method, improving convergence speed for solving frictionless contact problems in computational mechanics by applying an accelerated projected gradient approach.

Contribution

The paper proposes a novel accelerated Uzawa method that enhances convergence speed for contact mechanics problems compared to the traditional approach.

Findings

Faster convergence demonstrated in preliminary experiments

Application of accelerated projected gradient to dual problem

Potential for improved computational efficiency in contact mechanics

Abstract

The Uzawa method is a method for solving constrained optimization problems, and is often used in computational contact mechanics. The simplicity of this method is an advantage, but its convergence is slow. This paper presents an accelerated variant of the Uzawa method. The proposed method can be viewed as application of an accelerated projected gradient method to the Lagrangian dual problem. Preliminary numerical experiments suggest that the convergence of the proposed method is much faster than the original Uzawa method.

Equations

$$\min \{\mathbf{x},\mathbf{y}\}=(\min \{x_1,y_1\},\dots ,\min \{x_n,y_n\}{)}^{\top }.$$

$${\Pi }_{-}(\mathbf{y})=\mathrm{argmin}\{\Vert \mathbf{x}-\mathbf{y}\Vert \mid \mathbf{x}\in {\mathbb{R}}_{-}^n\}=\min \{\mathbf{y},0\}.$$

$$g_i(\mathbf{u})\ge 0,i=1,\dots ,m.$$

$$\mathrm{Minimize}$$

$$\mathrm{subject}\mathrm{to}$$

$$\displaystyle L(\boldsymbol{u},\boldsymbol{r})=\begin{dcases*}\pi(\boldsymbol{u})+\sum_{i=1}^{m}r_{i}g_{i}(\boldsymbol{u})&if $\boldsymbol{r}\leq\boldsymbol{0}$,\\ -\infty&otherwise.\end{dcases*}$$

$$\psi (\mathbf{r})=\inf \{L(\mathbf{u},\mathbf{r})\mid \mathbf{u}\in {\mathbb{R}}^d\},$$

$$\mathrm{Maximize}$$

$$\mathrm{subject}\mathrm{to}$$

$$\mathbf{r}={\Pi }_{-}(\mathbf{r}+\alpha \nabla \psi (\mathbf{r})).$$

$${\mathbf{r}}^{(k+1)}:={\Pi }_{-}({\mathbf{r}}^{(k)}+\alpha \nabla \psi ({\mathbf{r}}^{(k)})).$$

$${\mathbf{u}}_{\bar{\mathbf{r}}}=\mathrm{argmin}\{L(\mathbf{u},\bar{\mathbf{r}})\mid \mathbf{u}\in {\mathbb{R}}^d\}.$$

$$\frac{\partial }{\partial r_i}\psi (\bar{\mathbf{r}})=g_i({\mathbf{u}}_{\bar{\mathbf{r}}}),i=1,\dots ,m$$

$$\Vert {\mathbf{r}}^{(k)}-{\mathbf{r}}^{(k+1)}\Vert \le ϵ$$

$$\pi (\mathbf{u})=\frac{1}{2}{\mathbf{u}}^{\top }K\mathbf{u}-{\mathbf{p}}^{\top }\mathbf{u}.$$

$$g_i(\mathbf{u})=h_i-{\mathbf{n}}_i^{\top }\mathbf{u},i=1,\dots ,m.$$

$$\mathbf{g}(\mathbf{u})=\mathbf{h}-N\mathbf{u},$$

$$\mathrm{Minimize}$$

$$\mathrm{subject}\mathrm{to}$$

$$\alpha \in ]0,2{\lambda }_1(K)/{\sigma }_d(N)[,$$

$${\mathbf{r}}^{(k+1)}$$

$${\mathbf{\rho }}^{(k+1)}$$

$${\tau }_{k+1}$$

$${\omega }_{k+1}$$

$$\nabla \psi ({\mathbf{r}}^{(k)}{)}^{\top }({\mathbf{r}}^{(k+1)}-{\mathbf{r}}^{(k)})<0.$$

$$\Vert {\mathbf{\rho }}^{(k)}-{\mathbf{r}}^{(k+1)}\Vert \le ϵ$$

$$K\mathbf{u}=\mathbf{p}+N^{\top }{\mathbf{\rho }}^{(k)}.$$

$${\mathbf{e}}_1$$

$${\mathbf{e}}_2$$

$${\mathbf{e}}_3$$

$${\mathbf{e}}_4$$

Full text

Keywords

Fast first-order method, accelerated gradient scheme, contact mechanics, nonsmooth mechanics, Uzawa method, convex optimization.

1 Introduction

It has been recognized well that contact mechanics has close relation with optimization and variational inequalities [6, 27]. The static frictionless contact problem of a linear elastic body, also called Signorini’s problem, is one of the most fundamental problems in contact mechanics. This is a boundary value problem to find the equilibrium configuration of an elastic body, where some portion of the boundary of the body can possibly touch the surface of a rigid obstacle (or the surface of another elastic body). Positive distance between the elastic body and the obstacle surface (i.e., positive gap) implies zero contact pressure (i.e., zero reaction), while nonzero reaction implies zero gap. This disjunction nature can be described by using complementarity conditions. Moreover, the frictionless contact problem can be formulated as a continuous optimization problem under inequality constraints [27].

The Uzawa method is known as a classical method for solving constrained optimization problems [2, 5, 25]. Due to ease in implementation, the Uzawa method is often applied to contact problems [22, 18, 12, 10, 21, 26, 20, 19]. Major drawback of the Uzawa method is that its convergence is slow; it exhibits only linear convergence in general.

Recently, accelerated, or “optimal” [14, 15], first-order methods have received substantial attention, particularly for solving large-scale optimization problems; see, e.g., [3, 4, 16, 13]. Advantages of most of these methods include ease of implementation, cheap computation per each iteration, and fast local convergence. Application of an accelerated first-order method to computational mechanics can be found in [11]. In this paper, we apply the acceleration scheme in [3] to the Uzawa method.

The paper is organized as follows. Section 2 provides an overview of the necessary background of the frictionless contact problem and the Uzawa method. Section 3 presents an accelerated Uzawa method. Section 4 reports the results of preliminary numerical experiments. Some conclusions are drawn in section 5.

In our notation, ⊤ denotes the transpose of a vector or a matrix. For two vectors $\boldsymbol{x}=(x_{i})\in\mathbb{R}^{n}$ and $\boldsymbol{y}=(y_{i})\in\mathbb{R}^{n}$, we define

[TABLE]

We use $\mathbb{R}_{-}^{n}$ to denote the nonpositive orthant, i.e., $\mathbb{R}_{-}^{n}=\{\boldsymbol{x}=(x_{i})\in\mathbb{R}^{n}\mid x_{i}\leq 0\ (i=1,\dots,n)\}$. For $\boldsymbol{y}\in\mathbb{R}^{n}$, we use $\Pi_{-}(\boldsymbol{y})$ to denote the projection of $\boldsymbol{y}$ on $\mathbb{R}_{-}^{n}$, i.e.,

[TABLE]

We use $\mathop{\mathrm{diag}}\nolimits(\boldsymbol{x})$ to denote a diagonal matrix, the vector of diagonal components of which is $\boldsymbol{x}$.

2 Fundamentals of frictionless contact and Uzawa method

We briefly introduce the frictionless contact problem of an elastic body; see, e.g., [27] and [6] for fundamentals of contact mechanics.

Consider an elastic body subjected to a static load and a rigid obstacle fixed in space. Since the body cannot penetrate the surface of the obstacle, the deformation of the body is constrained from one side by the obstacle surface. We also assume the absence of friction and adhesion between the body surface and the obstacle surface. The set of these conditions is called the frictionless unilateral contact.

Suppose that the conventional finite element procedure is adopted for discretization of the elastic body. Let $\boldsymbol{u}\in\mathbb{R}^{d}$ denote the nodal displacement vector, where $d$ is the number of degrees of freedom of displacements. We use $\pi(\boldsymbol{u})$ to denote the total potential energy caused by $\boldsymbol{u}$. At the (unknown) equilibrium state, the nodes on a portion of the body surface can possibly make contact with the obstacle surface. Such nodes are called the contact candidate nodes, and the number of them is denoted by $m$. We use $g_{i}(\boldsymbol{u})$ $(i=1,\dots,m)$ to denote the gap between the $i$th contact candidate node and the obstacle surface. The non-penetration conditions are then formulated as

[TABLE]

The displacement vector at the equilibrium state minimizes the total potential energy under the constraints in (1). Namely, the frictionless contact problem is formulated as follows:

[TABLE]

The Uzawa method solving problem (2) is listed in Algorithm 1 [2, 5, 25, 1]. It is worth noting that step 3 of Algorithm 1 can be performed as the equilibrium analysis with a conventional finite element code. Also, step 4 is simple to implement. Due to such ease in implementation, the Uzawa method is widely used in computational contact mechanics.

As shown in Ciarlet [5, section 9.4], the Uzawa method can be viewed as a projected gradient method solving the Lagrange dual problem. Essentials of this observation are repeated here. The Lagrangian $L:\mathbb{R}^{d}\times\mathbb{R}^{m}\to\mathbb{R}\cup\{-\infty\}$ associated with problem (2) is given by

[TABLE]

Here, the Lagrange multipliers, $r_{1},\dots,r_{m}$, correspond to the reactions. The inequality constraints imposed on the reactions in (3) correspond to the non-adhesion conditions. Define $\psi:\mathbb{R}^{m}\to\mathbb{R}\cup\{-\infty\}$ by

[TABLE]

which is the Lagrange dual function. The Lagrange dual problem of (2) is then formulated as follows:

[TABLE]

Since $\psi$ is the pointwise infimum of a family of affine functions of $\boldsymbol{r}$, it is concave (even if problem (2) is not convex). Therefore, the dual problem (5) is convex.

Let $\alpha>0$ be an arbitrary constant. A point $\boldsymbol{r}\in\mathbb{R}^{m}$ is optimal for problem (5) if and only if it satisfies

[TABLE]

This fixed point relation yields the iteration of the projected gradient method as follows (see, e.g., [17]):

[TABLE]

Here, $\alpha$ plays a role of the step size. For given $\bar{\boldsymbol{r}}\in\mathbb{R}_{-}^{m}$, define $\boldsymbol{u}_{\bar{\boldsymbol{r}}}$ by

[TABLE]

It can be shown that the relations

[TABLE]

holds, when $\pi$ is convex and $g_{1},\dots,g_{m}$ are concave [5, Theorem 9.3-3]. Substitution of (9) into (7) results in step 4 of Algorithm 1. Thus, the Uzawa method for problem (2) is viewed as the projected gradient method applied to the Lagrange dual problem (5). It is worth noting that a reasonable stopping criterion,

[TABLE]

with threshold $\epsilon$, is derived from the optimality condition in (6).

If we assume the small deformation, the total potential energy is given by

[TABLE]

Here, $K\in\mathbb{R}^{d\times d}$ is the stiffness matrix, which is a constant positive definite symmetric matrix, and $\boldsymbol{p}\in\mathbb{R}^{d}$ is the external nodal force vector. Moreover, $g_{i}$ is linearized as

[TABLE]

For notational simplicity, we rewrite (10) as

[TABLE]

where $\boldsymbol{h}\in\mathbb{R}^{m}$ is a constant vector, and $N\in\mathbb{R}^{m\times d}$ is a constant matrix. The upshot is that, in the small deformation theory, the frictionless contact problem is reduced to the following convex quadratic programming (QP) problem:

[TABLE]

The step size, $\alpha$, of the Uzawa method for solving problem (11) is chosen as follows. Let $\lambda_{1}(K)$ and $\sigma_{d}(N)$ denote the minimum eigenvalue of $K$ and the maximum singular value of $\sigma_{d}(N)$. If $\alpha$ is chosen so that

[TABLE]

then it is shown that the Uzawa method converges [5, section 9.4].

3 Accelerated Uzawa method

In this section, we present an accelerated version of the Uzawa method.

The Uzawa method in Algorithm 1 is essentially viewed as the projected gradient method (7). This method is considered a special case of the proximal gradient method, because the proximal operator of the indicator function of a nonempty closed convex set is reduced to the projection onto the set [17]. Therefore, it is natural to apply the acceleration scheme for the proximal gradient method [3] to the projected gradient update (7). This results in the update

[TABLE]

Here, $\omega_{k}\in[0,1)$ is an extrapolation parameter which is to be determined so that the convergence acceleration is achieved. In [3], $\omega_{k}$ is chosen as

[TABLE]

with $\tau_{0}:=1$. The sequence of the dual objective values, $\{\psi(\boldsymbol{r}^{(k)})\}$, converges to the optimal value with rate $O(1/k^{2})$ [3].

The accelerated method in (12) and (13) is not guaranteed to be monotone in the dual objective value. Therefore, we follow O’Donoghue and Candès [16] in incorporating an adaptive restart technique of the acceleration scheme. Namely, we perform the restart procedure whenever the momentum term, $\boldsymbol{r}^{(k+1)}-\boldsymbol{r}^{(k)}$, and the gradient of the objective function, $\nabla\psi(\boldsymbol{r}^{(k)})$, make an obtuse angle, i.e.,

[TABLE]

As the upshot, the accelerated Uzawa method with adaptive restart is listed in Algorithm 2. One reasonable stopping criterion is

[TABLE]

at step 5.

In the case of the small deformation theory, steps 3 and 4 of Algorithm 2 are simplified as follows.

• •

Step 3: Let $\boldsymbol{u}^{(k)}$ be the solution to the system of linear equations

[TABLE]

Here, the coefficient matrix $K$ is common to all the iterations. Hence, we carry out the Cholesky factorization only at the first iteration; at the following iterations, we can solve (15) only with the back-substitutions.

• •

Step 4: $\boldsymbol{\gamma}^{(k)}:=\boldsymbol{h}-N\boldsymbol{u}^{(k)}$.

4 Preliminary numerical experiments

Consider an elastic body shown in Figure 1. The body is in the plane-stress state, with thickness $5\,\mathrm{mm}$, width $60\,\mathrm{mm}$, and height $20\,\mathrm{mm}$. It consists of an isotropic homogeneous material with Young’s modulus $200\,\mathrm{GPa}$ and Poisson’s ratio $0.3$. The bottom edge of the body is on the rigid obstacle, and the left edge is fixed by supports. The uniform downward traction of $50\,\mathrm{kPa}$ is applied to the top edge, and the uniform upward traction of $500\,\mathrm{kPa}$ is applied to the right edge. The body is discretized into $N_{X}\times N_{Y}$ four-node quadrilateral (Q4) finite elements, where $N_{X}$ $(=3N_{Y})$ is varied to generate problem instances with different sizes. The number of degrees of freedom of displacements is $d=2N_{X}(N_{Y}+1)$ and the number of contact candidate nodes is $m=N_{X}$. At the equilibrium state, about 73.3% of contact candidate nodes are in contact with nonzero reactions. We assume the small deformation, and solve QP problem (11). Computation was carried out on a $2.2\,\mathrm{GHz}$ Intel Core i5 processor with $8\,\mathrm{GB}$ RAM.

The proposed algorithm was implemented with MATLAB ver. 9.0.0. The threshold for the stopping criteria in (14) is set to $\epsilon=10^{-6}$. Comparison was performed with QUADPROG [23], SDPT3 ver. 4.0 [24] with CVX ver. 2.1 [9], and PATH ver. 4.7.03 [8] via the MATLAB Interface [7]. QUADPROG is a MATLAB built-in function for QP. We use an implementation of an interior-point method with setting the termination threshold, the parameter TolFun, to $10^{-8}$. SDPT3 implements a primal-dual interior-point method for solving conic programming problems. A QP problem is converted to the standard form of conic programming by CVX, where the parameter cvx_precision for controlling the solver precision is set to high. PATH is a nonsmooth Newton method to solve mixed complementarity problems, and hence can solve the KKT condition for QP.

Figure 2 reports the convergence history of the proposed algorithm (Algorithm 2). It also shows the result of the accelerated Uzawa method without restart scheme, and that of Algorithm 1 (i.e., the Uzawa method without acceleration). Figure 2 shows the convergence history of the dual objective function. It is observed that the acceleration and restart schemes drastically speed up the convergence. Particularly, with the proposed algorithm, the dual objective value seems to converge quadratically and monotonically. The residual in Figure 2 is defined by using the KKT condition. Namely, we define vectors $\boldsymbol{e}_{j}$ $(j=1,2,3,4)$ by

[TABLE]

and the value of $\|(\boldsymbol{e}_{1},\boldsymbol{e}_{2},\boldsymbol{e}_{3},\boldsymbol{e}_{4})\|$ is reported in Figure 2.

Figure 3 compares the computational results of the five methods. It is observed in Figure 3 that QUADPROG is fastest for almost all the problem instances. The computational time required by PATH is very small for small instances, but drastically increases as the instance size increases. On the other hand, it is observed in Figure 3 that PATH converges to the highest accuracy. For large instances, the proposed method is the second-fastest method among the five methods. The computational time of SDPT3 is larger than that of the proposed method. Also, SDPT3 converges to the lowest accuracy. The residuals of the solutions obtained by the proposed method, the unaccelerated Uzawa method, and QUADPROG are similar. The computational time required by the unaccelerated Uzawa method is very large.

The upshot is that, within the small deformation theory, the proposed method is quite efficient, but the interior-point method for QP (QUADPROG) is more efficient. If we consider large deformation, problem (2) is not convex, and hence QUADPROG and SDPT3 cannot be adopted. Efficiency of the proposed method, that solves the convex dual problem, remains to be studied.

5 Concluding remarks

This paper has presented an accelerated Uzawa method. The method can be recognized as an accelerated projected gradient method solving the Lagrangian dual problem of a constrained optimization problem. It has been shown in the numerical experiments that the acceleration and restart schemes drastically speed up the convergence of the Uzawa method. Besides, the proposed method is easy to implement.

Many possibilities of extensions could be considered. To update the primal variables, the Uzawa method solves a system of linear equations. For large-scale problems, this process might be replaced by a fast first-order minimization method of a convex quadratic function, e.g., [13]. Also, an extension to large deformation problems remains to be explored.

Acknowledgments

This work is partially supported by JSPS KAKENHI 26420545 and 17K06633.