A fast first-order optimization approach to elastoplastic analysis of skeletal structures

Abstract

It is classical that, when the small deformation is assumed, the incremental analysis problem of an elastoplastic structure with a piecewise-linear yield condition and a linear strain hardening model can be formulated as a convex quadratic programming problem. Alternatively, this paper presents a different formulation, an unconstrained nonsmooth convex optimization problem, and proposes to solve it with an accelerated gradient-like method. Specifically, we adopt an accelerated proximal gradient method, that has been developed for a regularized least squares problem. Numerical experiments show that the presented algorithm is effective for large-scale elastoplastic analysis. Also, a simple warm-start strategy can speed up the algorithm when the path-dependent incremental analysis is carried out.

Appendices

Appendix 1: SOCP formulation of problem (12)

In this section, we explain how problem (12) is recast as a second-order cone programming (SOCP) problem.

The second-order cone in ${\mathbb {R}}^{n}$ is defined by

$${ \mathcal {L}}^{n} = \left\{ (x_{0},x_{1},\dots ,x_{n-1})^{\top } \in {\mathbb {R}}^{n} \Bigm | x_{0} \ge \sqrt{x_{1}^{2}+\dots +x_{n-1}^{2}} \right\} . $$

SOCP is a minimization (or maximization) of a linear objective function under some second-order cone constraints and affine constraints.

The inequality constraints in (12c) can be written as second-order cone constraints as

$$\begin{aligned} \begin{bmatrix} \Delta \gamma _{i} \\ \Delta c_{{\text{p}}i} \\ \end{bmatrix} \in{ \mathcal {L}}^{2} , \quad i=1,\dots ,m . \end{aligned}$$

The constraints in (12b) are affine (i.e., linear equality) constraints. To convert the objective function to a linear one, we introduce auxiliary variables, $\xi \in {\mathbb {R}}$ and $\zeta \in {\mathbb {R}}$, that serve as upper bounds for the quadratic terms in (12a). Namely, we consider the following constraints:

$$\begin{aligned} \xi&\ge \sum _{i=1}^{m} \frac{1}{2} k_{i} \Delta c_{{\rm e}i}^{2} , \end{aligned}$$

(61)

$$\begin{aligned} \zeta&\ge \sum _{i=1}^{m} \frac{1}{2} h_{i} \Delta \gamma _{i}^{2} . \end{aligned}$$

(62)

The convex quadratic inequality constraint in (61) can be rewritten equivalently as (Ben-Tal and Nemirovski 2001)

$$\xi + 1 \ge \left\| \left[ {\begin{array}{c} {\xi + 1} \\ \sqrt {2k_{1} } \Delta c_{\rm e1} \\ \vdots \\ \sqrt {2k_{m} } \Delta c_{{\rm e}m} \\ \end{array} } \right] \right\|. $$

This is a second-order cone constraint. Constraint (62) can be rewritten in the same manner.

The upshot is that problem (12) can be converted to the following SOCP problem:¹

$$\begin{aligned} \displaystyle {\text {Minimize}} {\quad } \displaystyle & \sum _{i=1}^{m} q^{(t)}_{i} \Delta c_{{\rm e}{i}} + \xi + \sum _{i=1}^{m} R^{(t)}_{i} \Delta \gamma _{i} + \zeta - {\varvec{f}}^{\top } \Delta {\varvec{u}} \\ {\mathrm {subject\;to}}\, \displaystyle & \quad \Delta c_{{\rm e}i} + \Delta c_{{\text{p}}i} = \varvec{b}_{i}^{\top } \Delta {\varvec{u}} , \quad i=1,\dots ,m, \\&\displaystyle \\&\begin{bmatrix} \Delta \gamma _{i} \\ \Delta c_{{\text{p}}i} \\ \end{bmatrix} \in{ \mathcal {L}}^{2} , \quad i=1,\dots ,m , \\ & \displaystyle \begin{bmatrix} \xi + 1 \\ \xi - 1 \\ \sqrt{2 k_{1}} \Delta c_{{\rm e}1} \\ \vdots \\ \sqrt{2 k_{m}} \Delta c_{{\rm e}m} \\ \end{bmatrix} \in{ \mathcal {L}}^{m+1} , \quad \begin{bmatrix} \zeta + 1 \\ \zeta - 1 \\ \sqrt{2 h_{1}} \Delta \gamma _{1} \\ \vdots \\ \sqrt{2 h_{m}} \Delta \gamma _{m} \\ \end{bmatrix} \in{ \mathcal {L}}^{m+1} . \end{aligned}$$

Here, $\Delta c_{{\rm e1}},\dots ,\Delta c_{{\rm e}m}$, $\xi$, $\Delta \gamma _{1},\dots ,\Delta \gamma _{m}$, $\zeta$, and $\Delta {\varvec{u}}$ are variables to be optimized.

Appendix 2: Equivalence of (21) and (23)

As one of fundamental properties of the proximal mapping, we can show, for any $\alpha > 0$, that ${\varvec{p}} \in {\mathbb {R}}^{m}$ satisfies

$$\begin{aligned} {\mathbf{0}}&\in \nabla _{{\varvec{p}}}g_{1}({\varvec{v}},{\varvec{p}}) + \partial g_{2}({\varvec{p}}) \end{aligned}$$

(63)

if and only if it satisfies

$$\begin{aligned} {\varvec{p}}&= {\mathbf {\mathsf{{prox}}}}_{\alpha g_{2}} ({\varvec{p}} - \alpha \nabla _{{\varvec{p}}} g_{1}({\varvec{v}},{\varvec{p}})) . \end{aligned}$$

(64)

See, e.g., Parikh and Boyd (2014). For the reader’s convenience, essentials of the proof are repeated here.

Suppose that ${\varvec{p}}$ satisfies (63). This is equivalent to

$$\begin{aligned} {\mathbf{0}}&\in \alpha \nabla _{{\varvec{p}}}g_{1}({\varvec{v}},{\varvec{p}}) + \alpha \partial g_{2}({\varvec{p}}) \nonumber \\&= \alpha \nabla _{{\varvec{p}}}g_{1}({\varvec{v}},{\varvec{p}}) - {\varvec{p}} + {\varvec{p}} + \alpha \partial g_{2}({\varvec{p}}) . \end{aligned}$$

(65)

Let ${\varvec{s}} := {\varvec{p}} - \alpha \nabla _{{\varvec{p}}}g_{1}({\varvec{v}},{\varvec{p}})$ for notational simplicity. Then (65) is rewritten as

$$\begin{aligned} {\mathbf{0}}&\in \alpha \partial g_{2}({\varvec{p}}) + ({\varvec{p}} - {\varvec{s}}) , \end{aligned}$$

which is equivalent to

$$\begin{aligned} {\varvec{p}} = \mathop{\mathrm{arg\,min}}\limits_{{\varvec{z}}} \left\{ \alpha g_{2}({\varvec{z}}) + \frac{1}{2} \Vert {\varvec{z}} - {\varvec{s}} \Vert ^{2} \right\} . \end{aligned}$$

(66)

By definition, (66) is equivalent to (64).

Appendix 3: Algorithm for piecewise-linear hardening

We begin with computation of the gradient of $g_{1}$ defined by (58). In a manner similar to Sect. 3.3, it is convenient to define ${\varvec{e}} \in {\mathbb {R}}^{m}$ by

$$ {\varvec{e}} = B {\varvec{v}} - {\varvec{p}} - {\varvec{s}} , $$

which corresponds to the vector of incremental elastic elongation, $\Delta \varvec{c}_{{\rm e}}$, in problem (55). Then the gradient of $g_{1}$ can be calculated as

$$\begin{aligned} \nabla _{{\varvec{v}}}g_{1}({\varvec{v}},{\varvec{p}},{\varvec{s}})&= B^{\top } ( {\text{diag}} (\varvec{k}) {\varvec{e}} + \varvec{q}^{(t)}) - {\varvec{f}} , \\ \nabla _{{\varvec{p}}}g_{1}({\varvec{v}},{\varvec{p}},{\varvec{s}})&= {\text{diag}} (\varvec{h}_{1}) {\varvec{p}} - {\text{diag}} (\varvec{k}) {\varvec{e}} - \varvec{q}^{(t)} , \\ \nabla _{{\varvec{s}}}g_{1}({\varvec{v}},{\varvec{p}},{\varvec{s}})&= {\text{diag}} (\varvec{\eta }) {\varvec{s}} - {\text{diag}} (\varvec{k}) {\varvec{e}} - \varvec{q}^{(t)} , \end{aligned}$$

where

$$\begin{aligned} \nabla _{{\varvec{v}}}g_{1} = \frac{\partial g_{1}}{\partial {\varvec{v}}} , \quad \nabla _{{\varvec{p}}}g_{1} = \frac{\partial g_{1}}{\partial {\varvec{p}}} , \quad \nabla _{{\varvec{s}}}g_{1} = \frac{\partial g_{1}}{\partial {\varvec{s}}} . \end{aligned}$$

Moreover, the Hessian matrix of $g_{1}$ is written as

(67)

Since $k_{i}>0$, $h_{i1}>0$, $\eta _{i}>0$ $(i=1,\dots ,m)$ and B is of row full rank for a kinematically determinate truss, $\nabla ^{2} g_{1}({\varvec{v}},{\varvec{p}},{\varvec{s}})$ is positive definite. In a manner similar to Sect. 3.3, the proximal mapping of $\alpha g_{2}$ with $\alpha > 0$ can be computed as

$$\begin{aligned} {\mathbf {\mathsf{{prox}}}}_{\alpha g_{2}}({\varvec{w}},{\varvec{z}}) = \left[ {\begin{array}{l} {\text{diag}} ( {\mathrm {sgn}} ({\varvec{w}})) \max \{ |{\varvec{w}}| - \alpha \varvec{R}^{(t)}, {\mathbf{0}} \} \\ {\text{diag}} ( {\mathrm {sgn}} ({\varvec{z}})) \max \{ |{\varvec{z}}| - \alpha \varvec{R}^{{\text{s}}}, {\mathbf{0}} \} \\ \end{array} } \right] \end{aligned}$$

We are now in position to describe an accelerated proximal gradient method for solving problem (57).

At step 1 of Algorithm 3, auxiliary variables $\varvec{\varepsilon }_{l}$, ${\varvec{w}}_{l}$, and ${\varvec{s}}_{l}$ are used for convenience of computation.