{'fraction_non_alphanumeric': 0.047619047619047616, 'fraction_numerical': 0.016048602331788172, 'mean_word_length': 4.37200528202226, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 18, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Exploiting Lagrange duality for topology optimization with frictionless unilateral contact

This paper presents tractable reformulations of topology optimization problems of structures subject to frictionless unilateral contact conditions. Specifically, we consider stiffness maximization problems of trusses and continua. Based on the Lagrange duality theory, we derive formulations that do not involve complementarity constraints. It is often that a structural optimization problem with contact conditions is formulated as a mathematical programming problem with complementarity constraints (MPCC problem). However, MPCC usually requires special treatment for numerical solution, because it does not satisfy standard constraint qualifications. In contrast, to the formulation presented in this paper, we can apply standard optimization approaches. Numerical experiments of trusses and continua are performed to examine efficiency of the proposed approach.


Introduction
In this paper, we explore topology optimization of elastic structures (specifically, trusses and continua) subject to frictionless unilateral contact conditions. We assume that the structure can possibly make contact with the surface of a rigid obstacle. The obstacle is fixed in space, and the contact takes place without friction and adhesion. In this problem setting, we seek to find a structural design that has the maximal stiffness. This problem has quite a long history in the field of structural optimization; see a survey by Hilding et al. [14] for early contributions. This paper attempts to shed new light on this problem from a perspective of Lagrange duality theory.
An equilibrium state of a unilateral contact problem corresponds to a solution of a complementarity problem [16,35,47]. Therefore, it is often (and also natural) that a structural optimization problem with unilateral contact is formulated as a mathematical programming problem with complementarity constrains (MPCC problem); such formulations can be found in, e.g., [11,13,44,46]. Special treatment is usually required to solve an MPCC problem, because any feasible solution to an MPCC problem does not satisfy standard constraint qualifications [33]. A typical remedy is to apply reformulation and smoothing to the complementarity constraints so that the structural optimization problem can be handled with a conventional nonlinear programming (NLP) approach [11,46] or a gradient-based topology optimization approach [44]. Since the state variables (e.g., displacements, contact reactions, etc.) of a contact problem are in general nonsmooth with respect to design variables, another remedy is to apply a nonsmooth optimization method to the structural optimization problem. For example, Stavroulakis [41] used a bundle method, and Petersson and Patriksson [38] adopted a subgradient method. Alternatively, in spite of awareness of this nonsmoothness, a conventional gradient-based topology optimization approach (e.g., the method of moving asymptotes) were sometimes adopted [12,28]. Strömberg [42,43] performed sensitivity analysis by using the (augmented) Lagrangian for the contact problem. As another approach for dealing with complementarity constraints in a structural optimization problem, Hilding et al. [13] used a penalty interior-point method for MPCC.
Other than approaches based on complementarity problems, the unilateral contact conditions in a structural optimization problem have been treated (often in approximated manners) by using so-called penalty methods (e.g., with gap-elements) [9,25,26,32,34]. Lawry and Maute [31] adopted the so-called stabilized Lagrangian multiplier method in the extended FEM [10].
The approach presented in this paper is different from the ones in literature cited above. Namely, by using the Lagrange duality, we recast the topology optimization problem under consideration as a standard optimization problem (in the sense that, unlike MPCC problems, it satisfies standard constraint qualifications). Moreover, this approach does not resort to any approximation, which means that an optimal solution of the proposed formulation precisely satisfies unilateral contact conditions.
More specifically, we show that the stiffness maximization problem of trusses subject to unilateral contacts can be recast as a second-order cone programming (SOCP) problem. This is a convex optimization problem, and can be solved efficiently with a primal-dual interior-point method [2]. We next extend the presented formulation to continua. Here, we adopt the conventional solid isotropic material with penalization (SIMP) approach [4] with the density filter [5,6]. Due to the SIMP penalization, the formulation extended to continua is nonconvex. However, this formulation is suitable for application of a sequential SOCP approximation method, where the SIMP penalizations on densities are sequentially linearized. Thus, the proposed reformulation for continua can also be solved (in, in turn, a local sense) with a standard mathematical optimization approach (i.e., it does not require any special treatment for complementarity constraints).
From the perspective of convexity in truss problems, the approach presented in Kočvara et al. [30] is also of interest. Namely, they proposed to solve a convex optimization problem, which is obtained by eliminating the design variables (i.e., the member cross-sectional areas); similar approaches can be found also in [3,27]. In section 2.3, we establish clear relationship between the approach presented in this paper and the one in [30], by using the Legendre-Fenchel transform and the minimax theorem. This analysis provides us with a deeper understanding of the source of convexity in these two different approaches. From a practical point of view, one significant advantage of the approach in this paper over the one in [30] is that the former retains the member cross-sectional areas as optimization variables while the latter has eliminated them. Therefore, we can incorporate diverse constraints on truss design into the formulation in this paper. In the numerical experiments described in section 4.1, we demonstrate some concrete examples of such additional design constraints.
The paper is organized as follows. In section 2, we derive an SOCP formulation for trusses, by using the Lagrange duality theory. We also pursue investigation of relations between this formulation and existing other formulations [27,30,38]. Section 3 extends the formulation presented in section 2 to continua. In section 4, we perform numerical experiments on trusses and continua. In section 5, we draw some conclusions.
In our notation, denotes the transpose of a vector or a matrix. All vectors are column vectors. For notational simplicity, we often write (n + m)-dimensional column vector (x , y ) consisting of x ∈ R n and y ∈ R m as (x, y). For a vector x ∈ R n , the notation x designates its Euclidean norm, i.e., x = √ x x. We use 1 to denote an all-ones column vector. The notation S n designates the set of n × n real symmetric matrices. For a proper function f : R n → R ∪ {+∞}, its conjugate function is defined by where f → f * is called the Legendre-Fenchel transform (a.k.a. the Fenchel transform). We use R n + to denote the nonnegative orthant, i.e., R n + = {(x 1 , . . . , x n ) ∈ R n | x i ≥ 0 (i = 1, . . . , n)}. The n-dimensional rotated second-order cone, denoted K n , is defined by It is easy to verify that condition (x, y, z) ∈ K n is equivalent to namely, K n can be expressed with an n-dimensional second-order cone constraint.

Truss topology optimization
This section presents a convex formulation for truss topology optimization under the unilateral contact conditions.

Problem setting
Following the conventional ground structure method, consider an initial truss that consists of m members and has n degrees of freedom of the nodal displacements. Suppose that several nodes of the truss can possibly make contact with a rigid obstacle, and that the set of contact candidate nodes is specified in advance. Throughout the paper, we assume that the obstacle is fixed in space, and that contact between each contact candidate node and the obstacle surface is frictionless and adhesionless. Also, we assume linear elasticity and small deformation.
Let u ∈ R n denote the nodal displacement vector. We use c e ∈ R to denote the elongation of member e (e = 1, . . . , m). The compatibility relations can be written in the form c e = b e u, e = 1, . . . , m, where b e ∈ R n is a constant vector. Let l e and E denote the initial length of member e and the Young modulus, respectively. We use x 1 , . . . , x m to denote the member cross-sectional areas, which are considered design variables of the topology optimization problem. The stiffness matrix of the truss is given by Let c denote the number of contact candidate nodes. We use g j ≥ 0 (j = 1, . . . , c) to denote the initial gap between the jth contact candidate node and the obstacle surface. The non-penetration conditions for these nodes can be written as where C n ∈ R c×n is a constant matrix consisting of the unit inner normal vectors of the obstacle surface; see, e.g., [16,35,47] for fundamentals on kinematics in contact mechanics.
In structural optimization, the compliance is conventionally used as a measure of global flexibility of a structure. If all the prescribed nodal displacements are equal to zero, then the compliance can be defined as the external work done by the prescribed external nodal forces. Therefore, to maximize the stiffness of a structure, it is often that the external work is minimized. However, when there exists a non-zero prescribed displacement, it is known that the external work is not suitable as a measure of structural flexibility and the compliance should be defined by using the total potential energy (in the manner described below) [29,36,43]. Since the situation considered in this paper can possibly involves non-zero prescribed displacements at the contact candidate nodes (unless g = 0), we adopt the latter definition. Let f ∈ R n denote the prescribed external nodal load vector. For a given truss design x ∈ R m + , its compliance is defined by The truss topology optimization problem is then formulated as follows: Here, v > 0 is the specified upper bound for the structural volume. 4

Reformulation
In this section, we show that problem (5) can be recast as a convex optimization problem.
A key step for this is formulating the Lagrange dual problem of the maximization problem on the right side of (4), as performed below.
Theorem 1. For any x ∈ R m + , π(x) defined by (4) coincides with the optimal value of the following optimization problem: Here, w 1 , . . . , w m ∈ R, q 1 , . . . , q m ∈ R, and r ∈ R c are variables to be optimized, and the optimal value is defined to be +∞ if the problem is infeasible. 1 Proof. We first observe that, by substituting (1) and (2) into (4), π(x) coincides with the optimal value of the following optimization problem: subject to c e = b e u, e = 1, . . . , m, (6b) Here, u ∈ R n and c ∈ R m are variables to be optimized. For notational simplicity, let Associated with problem (6), the Lagrangian, L : where q ∈ R m and r ∈ R c are the Lagrange multipliers. Since all the constraints of problem (6) are affine, the Lagrange duality theory ensures the strong duality holds [7,8]: Here, the left side of (9) corresponds to problem (6), and the right side is its dual problem.
In the remainder, we show that this dual problem is equivalent to the minimization problem stated in this theorem. By direct calculations, for r ≤ 0 we obtain We can further reduce the first term of the last expression to where we have used x e ≥ 0, E > 0, and l e > 0. The proof is completed by substituting (11) into (10).
By using Theorem 1, we can recast problem (5) equivalently as follows: This problem is an SOCP problem; to see this explicitly, rewrite (12b) and (12e) as   l e /Eq e w e x e   ∈ K 3 , e = 1, . . . , m.
We can solve SOCP problems efficiently with a primal-dual interior-point method [2].
Remark 1. Constraint (12c) can be interpreted as the force-balance equation, where q e (e = 1, . . . , m) and r j (j = 1, . . . , c) correspond to the axial force of truss member e and the normal contact reaction at contact candidate node j, respectively. Constraint (12d) corresponds to the non-adhesion condition of the unilateral contact.
Remark 2. Extension of problem (12) to a multiple load case is obvious.
Remark 3. The form in (6) is crucial for deriving our convex reformulation in (12). For example, associate with the form on the right side of (4), we can define the Lagrangian asL for any r ≤ 0, the Lagrange dual problem is formulated as follows: Therefore, topology optimization problem (5) is equivalent also to the following problem in variables x ∈ R m , u ∈ R n , and r ∈ R c : This problem, however, has a nonconvex objective function and nonconvex equality constraints. Therefore, problem (12) is much preferable to problem (13). Alternatively, if K(ρ) is regular, then we can eliminate variable u from problem (13) by using (13b). This yields the following form: Here, > 0 is a small constant to avoid singularity of K(x). Since the objective function of this problem is nonconvex, problem (12) is again much preferable.

Relation with existing formulations
To obtain an optimal solution of problem (5), Kočvara et al. [30] proposed to solve a convex optimization problem that is different from problem (12). 2 Specifically, their approach solves an optimization problem with a linear objective function and some convex quadratic constraints. This problem can further be reduced to a quadratic programming (QP) problem, where a convex quadratic function is minimized under some linear inequality constraints. A similar formulation was presented also by Klarbring et al. [27]. 3 In this section, we establish the relation between this approach and our approach, i.e., problem (12). The analysis using the perturbation function for the problem under consideration attempts to provide deep understanding of the attribute of the problem. For notational simplicity, define X by We also use U and V defined by (7). For every x ∈ X, define Φ = η e (e = 1, . . . , m), It is easy to verify that Φ x is a closed proper convex function for any x ∈ X. By definition in (4), the compliance can be written as Thus, Φ is considered a perturbation function associated with the optimization problem on the right side of (4). The topology optimization problem (5) can be identified with finding x at which the value is attained. 4 Since −Φ x (u, c; 0, 0) is closed, convex with respect to x, and concave with respect to (u, c), the standard minimax theorem asserts [40] inf On the other hand, the duality theory asserts [8] sup from which problem (16) can be reduced to In a nutshell, the QP approach uses (17), as well as the facts that Φ x (u, c; 0, 0) is linear with respect to x and vertices of polyhedron X can be obtained explicitly. In contrast, our formulation in (12) is essentially based on (19). In the remainder, we give brief exposition of the two approaches. We begin with the QP approach. For notational simplicity, define U F by By direct calculations, the right side of (17) is reduced to From the standard theory of linear programming, we obtain where the vertices of X have been considered to evaluate the optimal value of a linear programming problem. Substitution of (21) into (20) yields the following optimization problem: This is minimization of a linear objective function under the convex quadratic constraints, and can be recast as an SOCP problem. Problem (22) can further be reduced to the following form: c e = b e u, e = 1, . . . , m, (23c) This is a QP problem (i.e., by converting maximization to minimization, we obtain a problem that minimizes a convex quadratic function under some linear constraints). 5 Then the right side of (17) can be written as sup u∈R n ψ(u).

Petersson and Patriksson [38] applied a subgradient method to this optimization problem.
We next see that the approach presented in this paper is essentially based on (18). For this purpose, we show that (9), used in the proof of Theorem 1, implies (18). From definition (14) of Φ x and definition (8) This relation is the key to showing which imply (18). Namely, for any x ∈ X and any (u, c) ∈ U , since Φ x (u, c; · , · ) is a closed proper convex function, its biconjugate function coincides with itself [8], i.e., which establishes (25). Moreover, application of the Legendre-Fenchel transform to Φ x yields which establishes (26). Thus, the difference between the QP approach in literature and the approach presented in this paper can be clearly captured in (17) and (19). One of advantages of the proposed formulation is that it retains the design variable x and hence we can incorporate various practicality constraints on truss design. 6 In contrast, it is very difficult for the QP approach to handle constraints concerning truss design, because in this approach the design variable x is eliminated; see (21).
Remark 5. It is worth of noting that (19) holds even if Φ x (u, c; η, λ) is not convex with respect to x. This allows extending the presented approach to continua, as seen in section 3.
Remark 6. For problem (16), approaches based on MPCC for problem can also be found in literature [11,13,44,46]. In these approaches, the KKT condition for the inner maximization problem of (16) is treated as constraint for the outer minimization problem.

Continuum topology optimization
This section extends the formulation in section 2 to the continuum-based topology optimization problem. We adopt the SIMP approach [4] with the density filter [5,6].
For simplicity, suppose that a planar design domain is discretized as regular mesh with square four-node quadrilateral (Q4) finite elements. Extensions to a three-dimensional case and other finite elements are straightforward.
Any notation unexplained in this section is the same as the one in section 2. Let u ∈ Rn and f ∈ Rn denote the nodal displacement vector and the external nodal load vector, respectively, wheren is twice the number of nodes. The displacement boundary conditions (i.e., the Dirichlet boundary conditions) are written in the form where D ∈ R d×n is a constant matrix. Note that the number of degrees of freedom of the nodal displacements isn − d. The non-penetration conditions of the contact candidate nodes are given as (3) with C n ∈ R c×n .
Let m denote the number of finite elements. We use ρ e ∈ [0, 1] (e = 1, . . . , m) to denote the density of element e, where ρ e = 1 means that element e exists and ρ e = 0 means that it is absent. The stiffness matrix, denoted K(ρ) ∈ Sn, is given in the form where p > 1 is the penalization power of SIMP, and K e ∈ Sn is a constant positive semidefinite matrix. Diagonalization of K e results in the form where κ ∈ S 5 is a positive definite constant matrix, and B e ∈ Rn ×5 is a constant matrix. Note that κ is identical for all finite elements because we suppose a regular mesh. The compliance is defined by The following theorem provides another expression of π(ρ), that is useful in topology optimization.
Theorem 2. For any ρ ∈ R m + , π(ρ) in (29) coincides with the optimal value of the following optimization problem: Here, if the problem is infeasible, then the optimal value is defined to be +∞.
Proof. By substituting (27) and (28) into (29), we see that π(ρ) is the optimal value of the following optimization problem: subject to c e = B e u, e = 1, . . . , m, (31b) Here, optimization variables are c 1 , . . . , c m ∈ R 5 and u ∈ Rn. Analogous to Theorem 1, the assertion of this theorem can be obtained as the Lagrange dual problem of (31).
Remark 7. Constraint (30c) is the force-balance equation, where s e ∈ R 5 is the generalized stress, t ∈ R d is the reaction vector stemming from the the Dirichlet boundary conditions, and r ∈ R c is the normal contact reaction vector. It is worth noting that s e is work-conjugate to the generalized strain c e defined by (31b). Constraint (30d) corresponds to the non-penetration condition of the unilateral contact law.
We are now in position to formulate the topology optimization problem. Application of the density filter can be written in the form where H ∈ R m×m is a constant matrix, x e ∈ [0, 1] is the original density, and ρ e is the so-called physical density. The topology optimization problem minimizes the compliance evaluated at ρ under the volume constraint: Application of Theorem 2 to problem (32) yields the following form: Minimize ρ, x, w, s, t, r m e=1 w e − 2g r (33a) subject to w e ρ p e ≥ s e κ −1 s e , e = 1, . . . , m, Here, ρ ∈ R m , x ∈ R m , w ∈ R m , s e ∈ R 5 (e = 1, . . . , m), t ∈ R d , and r ∈ R c are variables to be optimized. It is worth noting that an artificial small positive lower bound for x e , which is usually used in topology optimization to avoid singularity of the stiffness matrix, is unnecessary for this formulation.

Numerical examples
This section presents some numerical examples for the approaches proposed in section 2 and section 3. Numerical experiments were carried out on a 2.2 GHz Intel Core i5 processor with 8 GB RAM.

Example (I): Trusses
In this section, we solve the truss topology optimization problem in (12), with incorporating some additional design constraints. Associated with each truss member, we   introduce a binary variable to indicate whether the member exists or vanishes. Then, we can consider various design constraints, including lower bound constraints on the cross-sectional areas of existing members [23], an upper bound constraint on the number of nodes [20], limitation of the number of different values of cross-sectional areas [17,19],  and upper bound constraints on the degrees of nodes [21]. All of these constraints can be treated within the framework of mixed-integer second-order cone programming (MIS-OCP). SOCP and MISOCP problems were solved with CPLEX ver. 12.8.0 [15], where the problem data were prepared in the CPLEX LP file format by the code implemented in Matlab ver. 9.0.0.
As for MIP parameters of CPLEX, we set mip tolerances integrality (the amount by which each 0-1 variable can be different from an integer) to 10 −6 , mip tolerances mipgap (the relative tolerance on the gap between the best feasible objective value and the objective value of the best node remaining) to 10 −6 , emphasis numerical (the emphasis on numerical precision) to yes, emphasis mip to 4 (i.e., finding hidden feasible solutions is emphasized). Figure 1 depicts a problem setting. The nodes of the truss are aligned on a 1 m × 1 m grid. Any pair of two nodes is connected by members aligned in a straight manner, where overlapping of members is avoided by removing the longer member when two members overlap. Accordingly, the truss has m = 1361 members and n = 132 degrees of freedom of the nodal displacements. It has no node at which the nodal displacement is prescribed, i.e., d = 0. The number of contact candidate nodes is c = 20. The initial gaps between the contact candidate nodes and the obstacle surface are uniformly distributed, i.e., g 1 = · · · = g c . The Young modulus of the members is E = 20 GPa. The downward vertical external force of 100 kN is applied at the right bottom node. Figure 2 collects the obtained optimal solutions of problem (12) with various values of initial gaps. Here, the width of each member is proportional to its cross-sectional area, and filled circles indicate the contact candidate nodes that are in contact at the  Figure 3 is obtained. The computational results of these examples are listed in the first five rows of Table 1. Here, "obj." reports the optimal value, and "time" is the computational time spent by CPLEX. We next consider an upper bound constraint on the degree of each node [21]. 7 It is worth noting that, in every solution shown in Figure 2, the degree of the bottom center node is 7. We use δ max to denote the specified upper bound for the nodal degree. Also, we consider the lower and upper bounds, denoted x min and x max , respectively, for the member cross-sectional areas. Namely, we incorporate the constraints [17,23] x e ∈ {0} ∪ [x min , x max ], e = 1, . . . , m with x min = 500 mm 2 and x max = 5000 mm 2 . Figure 4 and Figure 5 collects the obtained optimal solutions with δ max = 5 and δ max = 4, respectively. The computational results are reported in Table 1, where "#BnB-node" is the number of enumeration nodes explored by CPLEX.
We focus attention on the solution in Figure 5d, which has a pair of mutually crossing members. By incorporating the constraints that prohibit presence of mutually crossing members [17], we obtain the optimal solution shown in Figure 6. The computational results are reported in the bottom row of Table 1.
Thus, we can solve optimization problems with various design constraints, because our formulation in (12) is convex and retains the member cross-sectional areas as explicit optimization variables.

Example (II): Continua
Consider a problem instance outlined in Figure 7. The rectangular elastic body is discretized as 80 × 32 mesh, i.e., m = 2560. There exists no node at which the nodal displacement is prescribed, i.e., d = 0 and n = 5346. The number of contact candidate nodes is c = 147. The initial gaps between the elastic body and the obstacle surface are uniformly distributed, i.e., g 1 = · · · = g c . We omit the units of quantities for simplicity, as often done in literature on continuum-based topology optimization. The Young modulus and the Poisson ratio are 1 and 0.3, respectively. The external force of 1 is applied at the right bottom corner. The sequential SOCP described in section 3 was implemented in Matlab ver. 9.0.0. The implementation of the SIMP approach was based on the Matlab 88 line code in [1]. The penalization power and the filter radius divided by the element size are p = 3 and 1.5, respectively. The specified upper bound for the volume fraction is 0.5. As for the initial point for the sequential SOCP, ρ (0) 1 , . . . , ρ (0) m are set to the volume fraction. We solved SOCP subproblems with SeDuMi ver. 1.3 [39,45], which implements a primal-dual interior-point method. Figure 8 collects the obtained solutions for various values of the initial gaps, when the downward external force is applied. Figure 9 collects the solutions obtained for the upward external force. Table 2 reports the computational results. Here, "obj." is the obtained objective value, "time" is the computational time, and "#iter." is the number of iterations of the sequential SOCP (i.e., the number of SOCP subproblems solved before convergence).
In Figure 8 and Figure 9, small open circles indicate the contact candidate nodes that are in contact at the equilibrium state; the remaining contact candidate nodes are free. We can observe that the larger the initial gaps, the wider the optimal configuration becomes.

Concluding remarks
This paper has presented new formulations for topology optimization of structures with frictionless unilateral contact. Specifically, for trusses, we have seen that the stiffness maximization problem can be recast as an SOCP (second-order cone programming) problem. Also, for continua, the stiffness maximization problem is formulated as a nonlinear SOCP problem. This nonlinear SOCP formulation is suitable for application of a simple successive linearization method, in which we solve an SOCP subproblem at each iteration. The key to deriving these formulations is the Lagrange duality in the optimization problem defining the compliance of a structural design. Unlike the convex optimization approach in Kočvara et al. [30] for trusses, the SOCP formulation presented in this paper retains the member cross-sectional areas as optimization variables. Therefore, various additional design constraints, e.g., limitation of the number of different cross-sectional areas [17] and limitation of the number of nodes [20], as well as the self-weight load [23], can be incorporated into this SOCP formulation. Problems with such constraints are formulated as MISOCP (mixed-integer second-order cone programming) problems, which can be solved globally with a standard mixedinteger programming solver. Extension of the presented approach to frame topology optimization with discrete design variables [18] is straightforward.