Minimization of maximum unilateral force

Computer methods in applied mechanics and engineering ELSEVIER

Comput. Methods Appl. Mech. Engrg. 177 (1999) 215-234 www.elsevier.com/locate/cma

Minimization of maximum unilateral force* D a n i e l H t l d l n" g " ", A n d e r s K l a r b r m g " a . . , J o n g - S h i o, , a n ~ b ~Department qf Mechanical Engineering. Division of Mechanics Linkiiping University, Linkiiping, S-581 83, Sweden hDepartment of Mathematical Sciences, The Johns Hopkins University Baltimore, MD 21218, USA Received 5 April 1998

Abstract The paper treats the structural optimization problem of minimizing forces arising from unilateral contact, and/or non-compressive constraints. We formulate this design problem as a mathematical program with complementarity constraints and establish some analytical results that pertain to the existence of an optimal design and the first-order necessary conditions for such a design. We apply three numerical algorithms and compare their performance on several examples. The computational results suggest that together these algorithms are well suited for solving this class of structural optimization problems. © 1999 Elsevier Science S.A. All rights reserved.

1. Introduction This paper deals with a large class o f structural optimization p r o b l e m i n v o l v i n g unilateral contact, a n d / o r n o n - c o m p r e s s i v e constraints. Specifically, we consider the theory, algorithms and numerical e x a m p l e s for the f o l l o w i n g unilateral optimal design p r o b l e m in the variables (s,u,p): minimize subject to

max Pi

(1)

s E b ° C R ''~'

(2)

K(s)u + CI, p = f

(3)

O<~g - C , , u l p > ~ O .

(4)

Eq. (3) and the c o m p l e m e n t a r i t y constraint (4) represent a unilateral structural m e c h a n i c s state problem, parametrized by the design vector s; n, is the n u m b e r o f design variables. The constraint (4) c o m p r i s e s of three conditions. The first, 0 <~g - C,,u, is a g e o m e t r i c condition that imposes a restriction on the n -dimensional d i s p l a c e m e n t - l i k e vector u; with C,, being a kinematic transformation matrix of order n, × n,, (where the subscript n in C,, refers to the " n o r m a l d i r e c t i o n " ) , the transformed displacement vector C,,u is restricted not to e x c e e d s o m e g i v e n distance vector g. The vector p is the force associated with the g e o m e t r i c condition; this force is restricted to be n o n n e g a t i v e as it pertains to the normal direction. Finally, the orthogonality condition

The research of Hilding and Klarbring was based on work supported by the Center tor Industrial Information Technology (CENIIT) and the Swedish Research Council for Engineering Sciences (TFR). The research of Pang was based on work supported by the United States National Science Foundation under grant CCR-9624018. *Corresponding author. Tel.: +46-13-281-117; fax: +46-13-281-101; e-mail: [email protected] 0045-7825/99/$ - see front matter © 1999 Elsevier Science S.A. All rights reserved. PlI: S0045-7825(98)00382-X

216

D. Hilding et al. / Comput. Methods Appl. Mech. Engrg. 177 (1999) 215-234

g - C,,u±p (the notation x l y

means that the two vectors x and y are perpendicular), along with the nonnegativity conditions, implies that a force component is active only if the geometric condition corresponding to that component is fulfilled as an equality. Eq. (3) represents the standard structural equilibrium equation for a linear elastic structure; K(s) is the stiffness matrix that depends on the design s, and f is the external load. The set 50 describes the constraints on the design. The cost function m a x ~ i ~ , , Pi represents the maximum of the force components; the minimization of this function is the objective of the overall structural design problem. The most natural example of the above problem is perhaps that of structural design under unilateral frictionless contact. In this case, the state problem (3) and (4) for a given design s may then be a set of conditions obtained after a finite element discretization of a continuous problem, or it may represent a naturally discrete problem such as that of a truss or frame structure with nodes that may come into unilateral contact. Another possible interpretation of the above design problem is that of designing a linear elastic structure supported by members such as ropes or cables that can only carry load in tension. Both the frictionless contact and the unilateral rope interpretations are used in this paper and these are the contexts within which the design problem (1)-(4) are derived. To minimize the peak or maximum contact force is a natural design objective in contact problems. This objective was used in some of the first papers on optimum shape design of such problems [3,4]. However, due to the difficulty in giving a mathematical meaning to this objective in the continuum case and due to the nondifferentiability of the maximum function, the emphasis was shifted towards minimizing equilibrium potential energy [1]. The rationale is that the objective of equilibrium potential energy can be shown to give a almost uniform contact pressure if the shape of the contact surface is the design variable [7]. In this paper we treat problems where the shape of the obstacle surface is not the design variable; for such problems minimizing equilibrium potential energy is an inappropriate objective. Indeed, the latter objective would imply minimizing stiffness; this would give rise to very flexible designs that are not particularly useful in practice. Therefore, we have chosen to go back to the objective of min max contact force. It should also be noted that for a priori discrete structures, the minimax objective causes no difficulty in terms of the mathematical well-posedness of the problem; this is in contrast to the continuum case. Furthermore, the problem applies naturally to the optimization of cable forces. The basic objective of the present paper is twofold: one, to present a rigorous mathematical treatment of several analytical aspects of the problem ( 1)-(4); and two, to present and compare several numerical algorithms for solving this optimal design problem. The presence of the complementarity condition (4) distinguishes this problem from a standard, smooth nonlinear program; indeed, the condition renders the computation of a global minimizer to the problem extremely difficult. The authors of the recent text [9] have undertaken a comprehensive study of the class of optimization problems with complementarity constraints of which our problem (1)-(4) is a special case. The term Mathematical Programs with Equilibrium Constraints (MPEC) was coined for such an optimization problem. In this text, detailed explanation is given to the computational difficulty of an MPEC in general; by using some advanced sensitivity results for piecewise smooth equations, first and second order optimality conditions are derived; and several solution approaches are described. In this paper, we shall rely on the results and algorithms obtained therein. The review paper [5] gives a fairly complete list of references concerning design optimization contact problems and also discusses in more detail the relation of these problems to MPECs. The theoretical difficulties of the optimal design problem are plenty; two major ones are: the nonconvexity of the feasible region and the nondifferentiability (and possible multi-valuedness) of the equilibrium state variables (u,p) satisfying (3) and (4) as a function of the design variable s. In spite of these difficulties, we are able to give a rigorous mathematical treatment to the problem. This treatment is presented in Sections 4 and 6. Specifically, we give a set of sufficient conditions for the existence of an optimal solution to this problem; we next present a set of first-order optimality conditions that are necessarily satisfied by a local minimizer of this problem. Following the analytical treatment, we present three numerical algorithms for solving the optimal design problem. Two of these, the implicit programming algorithm (IMPA) and the penalty interior point algorithm (PIPA) are adoptions of algorithms developed for solving general MPECs [9]. The third algorithm is the Method of Moving Asymptotes (MMA) [1 1]; this is a sequential programming method for differentiable nonlinear programs. The algorithm is known to work well for a large variety of structural optimization problems. Although MMA is only designed for differentiable optimization problems it is not difficult to apply it, with

D. Hilding et al. / Comput. Methods Appl. Mech. Engrg. 177 (1999) 215-234

217

suitable modifications, to our present problem which has an implicit programming formulation with a nondifferentiable objective function. The modification amounts to the use of directional derivatives in place of standard derivatives at points where the objective function is not differentiable. Such an approach has previously been presented for a similar problem in [8]. There are several important benefits for testing different algorithms. Being a nonconvex constrained optimization problem, the optimal design problem is computationally very challenging; since none of the algorithms can guarantee to produce a globally optimal solution, by comparing their outputs, we gain confidence in the quality of the solutions obtained. Moreover, in case an algorithm fails to terminate successfully, thereby yielding no solution to the problem at all, hopefully a different approach will produce a reasonable solution. Finally, the experience that we gather herein for this class of optimal design problems will help us understand the practical performance of the algorithms; knowledge like this will in turn be very useful when we consider applying these algorithms to other optimal design problems.

2. The state problem We consider a discrete linear elastic structure that is unilaterally supported. Our model includes two important situations; namely there are planar rigid obstacles that induce frictionless contact with the structure, and/or there are elastic cables that cannot support compression. The structure consists of n, interconnected nodes, n,, fixed plane obstacles and n r unilateral cables (or ropes). Obstacle i E {1 . . . . . n } , with unit normal vector fii which is either a 3-dimensional vector (for 3D problems) or a 2-dimensional vector (for 2D problems), may come into frictionless contact with node m(i) E {1 . . . . . n,,} of the structure. In Fig. 1 the notation for the contact condition is explained. The initial distance from obstacle i to the node m(i) is gi. The displacement of node j E {1 . . . . . n,,} is denoted d J; this is a vector in R 3 for 3-dimensional problems and a vector in R 2 for 2-dimensional problems. The scalar p~ is the value of the contact force in the negative normal direction. The frictionless small displacement contact conditions for each obstacle i E {1 . . . . . n } can be written as fli'd"'(i)--gi<~O,

Pi~O,

(5)

(~i'd"'(i)-gi)Pi=O.

These conditions express non-penetration of the obstacles, non-adhesion of the contact force, and the fact that this force must be zero when there is no contact. The central dot denotes the scalar product of geometric vectors. When interpreting these vectors as column vectors we may write ~i. d "(~)= ( ~ ) ' d ''(~, where the superscript t denotes transposition. The condition (5) is known as Signorini's unilateral contact condition. In the other situation, cable a E {1 . . . . . nr} is connected to node n ( a ) ~ {1 . . . . . n,}. The elongation of the cable is denoted by e,~ and the unit normal vector in the direction of the cable towards the node is r~'~. The force in the cable is f ~ and ~ denotes the initial slack. The following conditions describes the constitutive behavior of cable ol ~ {1 . . . . . nr}:

ft"'d"(")-e~-g.<~O,

p~>~O,

(fzi'd"('~)-e~-g~)p,~=O.

(6)

It is convenient to work with a vector u E R ''J that contains all the elongations of the cables and all the displacements of the nodes in the structure, where n z =- 2n,, + n,. in the 2D case and nj =-- 3n,, + n r in the 3D

i

~

m(O

F p~, I~ ~ ~e~, , , g~, rope %.

~(~)

nl Fig. 1. Illustration of an unilateral contact condition and a cable constraint.

218

D. Hilding et al.

Comput. Methods Appl. Mech. Engrg. 177 (1999) 2 1 5 - 2 3 4

case. Defining u = (d ~. . . . . d"",e I . . . . e,,,), we can then write the contact conditions (5) for all the n. obstacles and condition (6) for all n, cables jointly as ~ (7)

O ~ g - C,,u ± p >~ O.

Here, p = (p~ . . . . .

P,,,, f i l , . . . , f i , , , ) E R " < ,

where n =--n,, + n,., is the vector of contact and cable forces,

g = (gl . . . . . g,,, ff~ . . . . . if,,,) C R"< is the vector of initial gaps and slacks and

LC,,] where row i of the matrix C',, ~ R '''x'''' contains the components of r~i placed in the appropriate columns and row o~ of the matrix C,, ~ R '''×'''~ contains the components of r~" and - 1 placed in the appropriate columns. The contact force on node m ( i ) due to the obstacle i is -pir~ i and similarly for the cable force on node n(ce) due to cable ce. If the total contact and cable force on node i is f ' i , then

[f',, ..... f<,;,,,p, ....

,]

- Cl, p.

Let K be the stiffness matrix of the structure and cables. Since we are dealing with stable elastic material, we may assume that K is symmetric and at least positive semidefinite. Let f ~ - ( f ~ . . . . . f,,, 0 . . . . . 0) E R"" be a vector of external nodal loads. Force equilibrium together with the contact and cable constraints (7) now give the state conditions of a structure in frictionless contact with a rigid obstacle: Find (u, p) E R"" × R"' such that K u + C;,p = f O << g _ C,,uJ_p >~O '

(8)

where the first equation describes the force equilibrium. Problem (8) is a mixed linear complementarity problem which we denote by MLCP.

3. The optimal design problem In structural optimization the data of the state problem (8) are functions of the design, i.e. the stiffness matrix K, the matrix C, that reflects the geometry of the obstacles and cables, the load vector f and the gap vector g, may all depend on the design. In [6] the problem where the gap vector itself is the design variable was considered. Here, we will assume that the stiffness matrix depends on a design vector s while the other data is independent of the design; we denote such a design dependent stiffness matrix by K(s). This situation includes certain shape optimization problems where the shape is modified away from the contact surface as well as certain truss sizing problems where cross sectional areas of bars are optimized. In general, we need to impose constraints on the design variable s in order to formulate a reasonable design problem. One common constraint is of the box type: Sio w ~ S ~ Shig h.

(9)

For instance, if the structure to be optimized is a truss and s represents cross-sectional areas, Sh,w > 0 and Sh,gh introduce natural constraints on the size of each cross-sectional area. Another constraint that is reasonable to impose is that of a limit on the available amount of material to be distributed optimally in the structure. In shape optimization this constraint is generally nonlinear in s. However, for a truss or in the case of sizing of a sheet, the constraint becomes affine; in fact, it is of the form: tt s

c~is i <~ Maxvolume,

(10)

i=1

tWe use the following notation: if a, b E R" then a ~ b means that for all i = I. . . . . n: a ~ b. The meaning of a = b, a/> O, etc. is analogous to that of a ~ b.

D. Hilding et al. / Comput. Methods Appl. Mech. Engrg. 177 (1999) 215-234

219

where n, is the number of bars. Moreover, in this case, the stiffness matrix K(s) as a function of s also has a special form; specifically, tt s

I,:(s) =~s,I<, + R,

(ll)

i=l

where/~ E R '''~x'''~ is the stiffness of the bars and cables whose design is not to be altered and the K i ~ R '''~x'''~ are the stiffness matrices of the bars subject to change. Constraints (9) and (10) are restrictions on the design. Clearly, constraints involving the state, u or p or both, are also physically reasonable; for instance, constraints on the maximum allowed stress. However, since constraints of the latter type involve the state and design variables jointly, they impose considerable difficulty on the analytical and algorithmic aspects of the problem; for this reason, we shall not include such joint constraints in our treatment of the design problem. Summarizing the above description of the basic model, we obtain the formulation (1)-(4) for the optimal design problem of minimizing the maximum contact forces and/or cable forces of a linear elastic structure subject to frictionless contact, non-compressive constraints and design constraints. This mathematical program is a bilevel optimization problem with complementarity constraints, and as mentioned before, is a special case of an MPEC. The two conditions (3) and (4) define the state problem for a particular design s; this state problem is a parametric mixed linear complementarity problem which we denote by MLCP(s) in the subsequent sections.

4. E x i s t e n c e o f an o p t i m a l s o l u t i o n

The first analytical issue we will address is the existence of an optimal solution to design problem (1)-(4). For this purpose, we introduce the following assumptions: (A1) 5f is a nonempty compact subset of R"'; (A2) the obstacle constraint matrix C,, has linearly independent rows which we denote as C~,,, for i = 1. . . . .

(A3) (A4) (A5)

n;

there exists a feasible displacement vector u satisfying C,,u <~g; the stiffness matrix K(s) is symmetric positive semidefinite for all feasible designs s E ~; moreover, K(s) is a continuous function of s. For each s C

~/'(K(s)) f3 cg = {0}, where N(K(s)) denotes the null space of K(s) and ~ denotes the recession cone of the permissible displacements; i.e.

~ ~ {u ~ R"": C,,u ~<0}. Assumption (A1) clearly holds for the particular case of 5P discussed above; that is, when 5g is defined by (9) and (10). Assumption (A2) obviously holds for a contact problem involving at most one obstacle at each node as in the problem described in the previous section. Assumption (A3) is obviously needed for the design problem to be solvable. Assumption (A4) is true, for instance, in the case when K(s) is represented by (11). Assumption (A5) stipulates that no rigid body motions are allowed by the feasible designs; this assumption is trivially satisfied if the stiffness matrix K(s) is positive definite for s E 55. Notice that although the design set 5g is by assumption compact, this assumption alone is not sufficient to ensure the existence of an optimal solution to the problem. Another remark is that since the stiffness matrix K(s) is not assumed to be positive definite, there may exist multiple vectors u satisfying the state constraints defined by the complementarity system MLCP(s).

PROPOSITION 1. Under assumptions (A 1)-(A5), the optimal design problem (1)-(4) has a solution. PROOF. Assumptions (A4) and (A5) imply that for each s E ~ the state system MLCP(s) has a solution (u, p). Thus, the problem (1)-(4) has a feasible solution. Moreover, assumption (A2) implies that C,,C',, is a nonsingular matrix; thus,

220

D. Hilding et al. / Comput. Methods Appl. Mech. Engrg. 177 (1999) 215-234 p = ( C , , C i , ) - l ( f - K(s)u)

which implies that

[Jpl[ ~< ]l(c,,c;,)-'/l(][fll + ]l K(s)/I Ilu][).

(12)

Therefore, if we can show that u is bounded, then it follows that the problem ( 1 ) - ( 4 ) has a bounded feasible region and the proof will be completed. By way of contradiction, assume that there exists a sequence {(s k, u k, pk)} of feasible vectors to ( 1 ) - ( 4 ) such that lim

]]u' ]]--

k-->~

The bound (12) o f p in terms of u implies that the sequence {pk/[]u~]l} is bounded. Without loss of generality, let us assume that {(pk u%/llu~ll} converges to a pair (p~, u~); we may further assume, since ow is compact, that {s k} converges to a vector s ~ E ~ For each k, we have K(sk)u k + Ct,,p k = f O < ~ g - C , , u k ± pk>~O.

Normalizing by ]lukll and passing to the limit k--4 ~c, we deduce K(s~)u ~ + Ci, p ~ = 0 O>~C,,u ~ ± p ~ O .

Since u ~ # 0 (A5). []

and K(s ~) is symmetric positive semidefinite, the above conditions clearly violate assumption

5. The implicit programming reformulation It would be useful to present an equivalent formulation for the optimal design problem (1)-(4). For this purpose, we state a preliminary result which shows that the contact pressure vector p is a uniquely defined continuous function of the state variable s. This is in contrast to the displacement vector u which in general is a multivalued function of s. P R O P O S I T I O N 2. Under assumptions ( A I ) - ( A S ) , there exists f o r each s ~ S Y a unique vector p(s) satisfying MLCP(s); moreover, p(s) is continuous in s ~ 5~. PROOF. For any s E ow, the MLCP(s) is the K a r u s h - K u h n - T u c k e r optimality conditions of the convex quadratic program: ]

I

minimize

~u K(s)u - f t u

subject to

C,,u >1g.

By a result of Mangasarian, the gradient of the objective function, K(s)u - f is constant for all optimal solutions u. Thus so is CI, p. Since CI, has linearly independent columns, it follows that p is unique. Denoting this unique multiplier as p(s), we can easily establish that p(s) is continuous in s, by using this characterization of p(s). [] Employing the uniquely defined (but only implicitly known) function p(s), we may formulate the design problem ( 1 ) - ( 4 ) equivalently as a constrained (nonsmooth) optimization problem in the state variable s alone: minimize subject to

max pi(s) sC

(13)

This compact formulation is useful in two major ways. Algorithmically, it is the basis for the IMPA;

D. Hilding et al. / Comput. Methods Appl. Mech. Engrg. 177 (1999) 215-234

221

analytically, it is the basis for the derivation of the first-order optimality conditions for a minimizer of the optimal design problem. It should be noted that although the max operator is a convex function of its arguments, the objective function of (13) is not a convex function of the state variable s; the reason is that each pi(s) is in general not convex in s. Consequently, it is not realistic to expect that one can actually compute the global minimum of (13) in practice. The objective function maxj<_~<_,,,p~(s) in (13) is clearly nondifferentiable due to the max operator. This complication may be removed by introducing one extra variable and n, extra constraints. That is, we have the following equivalent formulation: minimize

p ....

subjectto

pi(s)--Pmax~O, s ~ bO CR",

i-1 . . . . . n,

(14)

6. First-order stationarity conditions We derive the first-order stationarity conditions of the problem (1)-(4) via the equivalent formulation (13). Since every local minimizer must be a stationary point, these conditions must necessarily be satisfied by every local minimizer of the design problem. We will focus on the case where the set 50 is polyhedral and of the form:

~=-{s~R",:As<~b} where A E R m×'',' and b E R'". Further, for simplicity, we will also take the stiffness matrix K(s) to be given by (11), where each K~ a n d / ( are symmetric positive semidefinite matrices. This setting, though simple, suffices to illustrate the main approach for obtaining the stationarity conditions in general. Let s* E 5e be a feasible design; assume that the matrix K(s*) is nonsingular, thus positive definite. (It is possible to relax this positive definiteness assumption somewhat; for simplicity, we choose not to do so in order to keep the discussion and notation simple.) Under this assumption and also by Proposition 2, there exists a unique pair (u*, p*) that solves the state problem MLCP(s) with s = s*. We define three fundamental index sets associated with the triple z*-= (s*, u*, p*):

o~(Z*) ~(z*) y(z*)

~ {i:(g-C,,u*)i=O O = p ~ * } .

For an obvious reason, we call fl(z*) the degenerate index set associated with z*. If fl(z*) is empty, then (u*, p*) is a nondegenerate solution of the state problem MLCP(s) corresponding to the state s*. The following result gives an important consequence of the positive definite assumption of the matrix K(s*). This result follows from the sensitivity theory of parameter dependent complementarity problems. It uses the notion of B-differentiable functions; a function is B-differentiable if it is Lipschitz continuous and directionally differentiable [9].

PROPOSITION 3. Let K(s) be given by (11). Suppose that K(s*) is positive definite. There exists a neighborhood 2¢" of s* and B-differentiable functions u: Y--~ R"" and p: ?¢'--+R"' such that for each s ~2(, (u(s), p(s)) is the unique pair of vectors satisfying MLCP(s); moreover, for every vector d ~ R"', the directional derivatives, u "(s* ; d) and p "(s* ; d), of the functions u and p at s* along d, respectively, are the unique pair of vectors (v, r)satisfying i1 s

diKiu* + K(s*)v + i--I

O= Ci, v,

~ i E ~ ( = * ) U 13(z*)

i ~ a(z*)

O>~Ci,,o±ri>~O, i ~ f l ( z * ) O = r i,

i~y(z*).

Ci,,r i = 0

D. Hilding et al. / Comput. Methods Appl. Mech. Engrg. 177 ( 1 9 9 9 ) 2 1 5 - 2 3 4

222

If the triple z* ~-(s*, u*, p * ) is a local minimizer of (1)-(4), then s* is a local minimizer of the implicit program (13). Thus, s* is a stationary point of the latter program. If K(s*) is positive definite, then each function pi(s) is B-differentiable at s*. Letting 5~* ~ - { j : p *

= max,_
and ~ / * be the binding indices at s*, we conclude that a necessary condition for z* to be a local minimizer of ( 1 ) - ( 4 ) is that d = 0 is the global minimizer of the following minimization problem in the variable d: r

#

minimize

max Pi(S ; d)

subject to

A id <~O,

i~l':

i C ag *.

where A i denotes the ith row of the matrix A. Based on this observation, we obtain the following result which gives a set of first-order necessary conditions for z* to be a locally optimal triple of the design problem. The remaining proof of the result is given in the Appendix A. PROPOSITION 4. Let s* ~ 5~ be such that K(s* ) is positive definite. Let (u*, p* ) be the unique pair satisfying MLCP(s) at s*. Then z* ~- (s*, u*. p* ) is a stationary point o f ( I ) - ( 4 ) ~f and only if for every pair (~1, 13_,) of index sets that partition ~(z*), there exists (A, ~, r r ) ~ R '''+''''+''' such that (Aj)tA+~:'~u*=0,

j=l

. . . . . n,

t

K(s)( + C,,rr = 0 Ci,, sc <~O,

i E fll

Ci,,~>~O,

i~J*

E

1

i C,] *

Ci,(=O,

i E c r ( z * ) ~ 5 ~*

O<~A k b - As* >t0 ~>-0,

iE~2

~. = 0,

i E y(Z* ).

The above stationarity conditions simplify significantly in the case where the degenerate index set /3(z*) is empty and the argmax set ..¢* is a singleton. We state this special case in the corollary below. C O R O L L A R Y 1. Let s* ~ 5~ be such that K(s*) is positive definite. Let (u*, p*) be the unique pair satisfying MLCP(s) at s*. Suppose that p* + (g - C,,u*) > 0 and that J * is a singleton, say 5~* = {i*}. Then z* -= (s*, u*, p * ) is a stationary point of ( 1 ) - ( 4 ) tf and only if there exists (A, sc, rr)E R ...... "+"' such that (a/)'a+sC'Kiu*=0,

j=

K(s)¢ + CI, rr = 0 Ci.,,~= 1 Ci,,( = 0,

i E ~(z*)\{i*}

O<~ a k b - As* >~O = O,

i E y(z*).

1. . . . . n

D. Hilding et al. / Comput. Methods Appl. Mech. Engrg. 177 (1999) 215-234

223

7. The implicit programming algorithm (IMPA) The implicit programming approach to the unilateral force optimization problem is based on the formulation (13). We assume throughout that K(s) has the additive form ( i 1) and K(s) is symmetric positive semidefinite for each s E A general iteration of the IMPA is as follows. Let s k E 5g be a given feasible design and pk =p(sk). Also, let u k be such that (pk, u k) solves the MLCP(sk). Without loss of generality, we may assume that pk ¢ 0 (otherwise, the design problem is solved with the optimal contact force identically equal to zero). Define several index sets: d k

[3k

--=

{i:p~=max

=-

{j:

p~}

=-- {]: g . i - C j . u ' = O = p ~ } =-

Since p k ¢ 0, we have d k C_ ak. Let Qk ~ R"" +"'~×~"' +"'~ be a given positive definite matrix. We compute a search direction ds k by solving the following convex minimization problem in the variables (ds, dp, du), which can easily be converted into a convex program: minimize

i/dsV /ds\ max @i +v-~@) Q@dp)

subject to

sk + ds E Se

i E:.~4k

~'~ds i K iu k + K(s k) du +

Z

1

c,,, dpj = 0

i=l

O=Ci, du,

jEce k

0 >>-Ci,, d u ] O<~dpi ~

jC~

0 = dpi,

jE~.

(]5)

The following lemma asserts that this quadratic program always has an optimal solution with a nonpositive objective value.

LEMMA 1. Let Qk be symmetric positive definite; assume that the matrix

K(s*) is nonsingular. The following two statements hold." (a) the directional subprogram (15) has an optimal solution with a nonpositive optimum objective value; and (b) if the optimum objective value of (15) is zero, then s k is a stationary point of (13). Let (ds k, dp k, du k) be an optimal solution of (15). Suppose that the optimum objective value of (15) is negative. This implies max i E .,4t

< o.

Under certain theoretical hypotheses (such as the linear independence constraint qualification), it can be shown that the so-obtained vector ds k is a direction of descent for the max contact force function max~i~,, ' pi(s) at s k. In practical implementation, we determine a stepsize T~ @ [~,1] as follows, where ~ > 0 is a prescribed lower bound on the stepsize, say e E (0, 0.1 ]. Set

D. Hildil g e t al. / (_'ompttt. Methods Appl. Mech. Engrg. 177 (1999) 215-Z74

224

s'~(l ") ~-- s 1' + "rds k

and

-~p(s~(r)).

if(r)

Notice that to obtain p~(r), we need to solve the MLCP(s~(r)). For given positive scalars p, ~rc~C (0, 1 ), c > 0, let r~ =- p'"' where m k is the smallest nonnegative integer m such that (a) 7- ~- p'" ~< 6-, or (b) the following Armijo inequality holds: /,

max pi(r) <- max p ~ + ~ r r

maxdp).

(16)

To determine m k, we start with m = 0 and check (a) and (b), increasing m by one when both of these conditions fail. If condition (b) holds, then we have computed a new design s k * ~ =- s ~ + r~ ds ~ with a lower maximum contact force than s ~. Notice that it is possible that (a) holds for a stepsize r~ before (b) does. In this case, we k+l accept the new design s anyway. The above iterative procedure is repeated until an iterate s'k is reached with the (negative) optimal objective value of the problem (15) not exceeding a prescribed positive tolerance. A summary of the IMPA is presented below.

The IMPA for program (13) Step 0 (hHtiali:ation). Select an initial design s ° E 5'] a symmetric positive definite matrix Qo, line search parameters p,y ~ (0,1 ) and set k = 0. Step 1 (State problem). Solve MLCP(s k) for (U(Sk),p(sk)). Step 2 (Direction generation). Let (ds k, dtt k, d f f ) be the solution of the quadratic subprogram (15). Step 3 (Armijo Linesearch). Let m k be the smallest integer m/> 0 such that the Armijo inequality (16) holds. Set the step length r~ = p"~ and new design s k+ ~= s~+ r~df. Step 4 (Termination check). Check the prescribed rule for termination, if it is satisfied then s k' ~ contains the final design else set k = k + l, create Qk and go to step 1.

8. The penalty interior point algorithm (PIPA) An alternative approach for solving the minimax optimal design problem is by the Penalty Interior Point Algorithm (PIPA) as described in [9]. This algorithm computes a B-stationary point of the problem based on the following formulation: minimize subject to

max p~ sEJ

K(s)u + C t p - f = 0 u+C,,u-g=O

(17)

diag(v)p = 0 (v, p)~>0, where diag(v) is the square diagonal matrix formed from the vector v. A general iteration of the algorithm is follows. An iterate (s k, u ~, i f , v k) that satisfies the following properties is given: (PI) (feasible design) s k ~ J'; (P2) (positivity of gap and pressure) (v k, i f ) > 0; and (P3) (centrality condition) diag(v ~) pk ~> p/z~l, where p E (0, 1) is a given constant, 1 is the vector of all ones and

( u~ )rpk ~-Let

D. Hilding et al. I Comput. Methods Appl. Mech. Engrg. 177 (1999) 215-Z74

F* =--K(sk)u * + C',,p k - f

225

r k =- v k + C,,u k - g ,

and

and define the index set:

=

s~k =- { i ' p ~

max

I -"---/~n

pkt. #j

A search direction is computed by solving the following convex quadratic program in the variables (ds, du, dp, dr): minimize

! ~

max dpi + ½ ds Qk ds i ~ .~4k

subject to

sk + ds E NdsH~< c(llFkll + Nrkll + (v~)'p ~) i

(18)

ds, K,u k + K(s k) du + C;, dp = 0

F k +~ i=l k

r +dv+C,,du=O diag(v k) dp + d i a g ( p k) dv + diag(vk)p ~ = o'dxk 1, where c and o-k are positive scalars with o-k < 1, and Qk ~ R"'×"' is a given positive definite matrix. In practical implementation, we can use the last three equations in the constraints to solve for the variables (du, de, dp) in terms of ds; the above quadratic subprogram thus reduces to one in the variable ds only. Define

uk(r)/ is kk'~ pk(r)l

I;kl

du

+ r

vk,,l \o 1

,

rE

[0,l];

t d:)

and let t M r ) =-

(Vk(T)) ~pk(r) n c

We seek a step size r k E (0, 1] so that the new iterate

(sk+l, Uk+l, pk+,, Vk+l)=--(s k (rk),u k (rk), P k (rk),v /, (~k)) satisfies the positivity and centrality copditions and a penalty function decreases in value. Notice that we have sk(r) E 5p for all r C [0, 1]. Specifically, first determine the scalar ~k--:sup

r E (0,1]:

min [vj(r)pi(r)-plxk(r)]>~O

.

(19)

This can be accomplished very easily because for each i, v~(r)pT(r) is a quadratic function in r. The following penalty update rule is then applied: with a k_ ~ > 0 being the current value of the penalty parameter (penalty update rule). Let ~ be the smallest integer ( / > 1 such that (

max d p , - 2ak_,(ilFkl] 2 + Llrkll 2 + '~(l _

\z o-~)tv

kxl k ) p

<

-

(Nfk][ 2 + II/112 +

(v%'p%.

i ~ .c'[k -Let cek = a k[k_ ~. Define the penalty function:

P.~ (s, u, p, v) =- max pi + cek(liK(s)u + CI, p _ fl]2 + I]v + C.u - g[]2 + v'p). l ~ i <~

Perform the

Armijo inexact line search: for a given backtracking factor 6 E (0, 1) and a constant y ~ (0, 1), let r k =--0.9999 8 ~~ ~

(20)

D. Hilding et al. I Comput. Methods AppI. Mech, Engrg. 177 (1999) 215-234

226

where t~£ is the smallest nonnegative integer ( such that with T --= 0.9999 6 ( ~, we have

P,,~(s k (z), u k (T), pk(z), vk('r)) - P,~(s k, u k , p k , v k) ~< yz[max•,~.~.¢,d p , - 2 a ;

,(IIF*II 2 + t1711: + #(1_ - a-~)(vk)'pk)].

(21)

This completes the description of iteration k of the PIPA. A summary of this algorithm is presented below.

The P1PA for program (17) Step 0 (Initialization). Set centrality parameter p ~ (0, 1), o-o E (0, 1), c > 0, a j > 1, linesearch parameters 8, y E (0, 1), a symmetric positive definite matrix Qo and set k = 0. Select an intial design s ° E ow and (u °, pO, v °) fulfilling (P1-3). Step 1 (Direction generation). Let (ds k, du k, dp k, du k) be the solution of the quadratic subprogram (18). Step 2 (Step reduction). Determine the maximum step length ~k maintaining (P2-3) from (19). Step 3 (Penalty update). Let ga/> 1 be the smallest integer such that the inequality (20) holds. Set the new penalty to a k = a~L j. Step 4 (Armijo line search). Let (~ be the smallest integer (~> 0 such that the Armijo inequality (21) holds for v = 0.9999 6 ( _zk. Set the step length z, = 0.9999 6 (t~k and the new iterate (s k+ 1 u k + l p k + l , vk+l ) = (s k u k pk, u k ) + "r~(ds*, du*, dp k, dvk). Step 5 (Termination check). Check the prescribed rule for termination, if it is satisfied then s k+~ contains the final design else set k = k + 1, create Qk, o-~ E (0,1) and go to step 1.

9. An alternative to PIPA and IMPA During the computations, both the PIPA and IMPA occasionally showed an annoying tendency to zig-zagging that leads to slow convergence. We can overcome this difficulty by the following simple modification. In the subprograms (18) and (15) replace max dp i i ~ "'~Ik

with max (Pi + dPi) i C .'.~k

where

and 0 ~< r/~< 1 is a constant (we used r / = 10 4). In the forthcoming sections we will refer to the modified (alternative) algorithms as PIPAA and IMPAA respectively.

10. The Method of Moving Asymptotes (MMA) For benchmarking of PIPA/PIPAA and IMPA/IMPAA we use the Method of Moving Asymptotes, MMA, Svanberg [11], which is a sequential programming method for differentiable nonlinear programs. The algorithm is known to work well for a large variety of structural optimization problems. We note, however, that MMA is not suited to solve our design problem ( 1 ) - ( 4 ) directly. It is better to apply it to the optimization formulation (14). MMA requires derivatives of constraints and objective. For non-differentiable points of P i ( ' ) one

D. Hilding et al. / Comput. Methods Appl. Mech. Engrg. 177 (1999)215-234

227

possibility is to use one of the directional derivatives instead. A detailed description of MMA is given in Appendix B.

11. Computational experiments As an evaluation of the algorithms we present below their performance when applied to a suite of test problems. As the algorithms are quite different it is difficult to find a completely fair measure of the performance and computational costs. We have chosen to report the following statistics: • # e q : Number of factorizations. Due to the implementation, see Section 12, the main work, besides solving the direction finding subprogram, in each algorithm is to factor a ( n + n ) × (n,, + n ) matrix. To factor this matrix is of similar complexity for all algorithms. • # subs: Number of direction finding subprograms that have been solved. For all algorithms this problem is of similar complexity, because the subprograms in PIPA/PIPAA and IPA/IPAA can all be reduced to subprograms of the same size as the one in MMA, with n, + 1 variables; see also Section 12. • maxpi: Value of the objective function after optimization. The following statistics, which do not apply to MMA, are also reported • # line: Number of iterations in the Armijo line search procedure. • rain ~': Minimum step length during the optimization. • avg T: Average step length during the optimization. The following stopping criteria was used for MMA and IMPA/IMPAA:

Hpk ' - p k l ] < l

× 10 -6"

(22)

For PIPA/PIPAA the stopping criterion was

lip ~-' - p k ] l < 1 × lO -6

(23)

II(Ku + Cl, p

(24)

and - F, min(g - C,,u,p))N < 1 × 10 4.

The condition (24) assures that when the computations are terminated with some design s* the state variables (u*, pk) fulfill the state equation MLCP(s*) with sufficient accuracy. In all figures representing mechanical structures solid lines are bars whereas dashed lines are elastic cables. Both cables and bars have a Young's modulus of 200 GPa. All elastic ropes have a cross section area of 0.02 m z and the letter F is used to represents a force with magnitude I MN. In all problems the bar cross-section areas are taken as design variables and a lower and upper limit of 0.001 m 2 and 0.1 m-', respectively, on these are used. 11.1. Test p r o b l e m 1

The structure illustrated in Fig. 2 is optimized for two different initial slacks in the rope, (i) 0.001 m and (ii) 2 HI

V: II

3 0.5 m

Fig. 2. Test problem 1, two bars and an elastic rope supporting the force F.

228

D. Hihling et al. / Comput. Methods Appl. Mech. El grg 177 (1999)215-Z74

Table I Statistics for test problem I with initial slack 0.001 m (degenerate case) Property

PIPA and PIPAA

IMPA and IMPAA

MMA

# eq # subs

36 36

6 5

8 7

0

0

-

0.294 0.969 1.23 x 10 "

I 1 3.93 X 10 "

8.51 X 10 "'

#

line

r avg r max p, (MN) rain

Table 2 Statistics lbr tesl problem I with initial slack 0.0005 m (nondegenerate case) Properly

PIPA and PIPAA

IMPA and IMPAA

MMA

# eq # subs

22 22

9 8

10 9

0

0

-

0.575 0.981 1.132281

I I 1.132281

1.132281

#

line

r avg r max Pi (MN) nfin

0.0005 m. The objective ft, nction max, p~ is in this case the m a x i m u m tension in the rope. The v o l u m e is constrained to be less than 0.02 m ~ and the initial design is s ° = (0.0015, 0.0025) m e. The value of the objective function is for the initial design in case (i) 1.6714 M N and in case (ii) 2.3531 MN. In both cases the problem has one global o p t i m u m , the difference is that in case (i) the o p t i m u m solution is degenerate whereas in case (ii) it is not. Such a solution is said to be degenerate if both the rope tension and the slack is zero in equilibrium. This test p r o b l e m is interesting because theoretically the algorithms may experience difficulties when applied to -k p r o b l e m s with degenerate solutions. The matrices Qk and Q were chosen to be the scaled identity matrix 1 x 10 7diag(1) lor all algorithms. The statistics from the different algorithms is shown in Tables 1 and 2. Note that for this problem there is no difference in p e r f o r m a n c e between PIPA and PIPAA or IMPA and IMPAA. 11.2.

Test p r o b l e m

2

The rope suspended truss structure is illustrated in Fig. 3. As in test problem 1, the objective function in this case is the m a x i m u m tension in the ropes. The initial slack of the ropes is zero, the m a x i m u m allowed v o l u m e of the truss is 1.5 times that of the initial design: .¢] = 0.02 m-,~ i = 1 , . , ., n.. The matrices Q~ and 0 a were chosen to 1 × 10 "~diag(1) for all algorithms. The execution statistics is s u m m a r i z e d in Table 3; the results clearly

, FI

IFI

; X××

x×x

t

. . . . .

,._x × ×

--,i.

gT

l-'ig.

3. Test

rope~,.

l*~roblenl

"~ a truss .;,

in

s u p p o r t e d by. a n u n l b e r

of

elastic

flz

ga

IF

i

NIl lm

.(]4 .(],5

Fig. 4. Test problem 3, a truss in contact with a number o f rigid obstacles (5 X 5 nodes case shown.)

D. Hihli# g e t al. I Comlmt. Methods Appl. Mech. El grt, 177 (1999) 215-Z~4

229

Table 3 Statistics for lest problem 2 Property

PIPA"

PIPAA

IMPA

IMPAA

MMA

# eq # subs

I01 101

37 37

1082 82

10 7

8 7

# line

()

0

965

rain 7avg 7max p, ( MN )

2.67 × 10 <' 0.301 1.820273

I 1 1.575103

3.81 × I0 +' 0.0393 1.575108

2

-

0.5 0.857 1.575103

1.575103

" Faihn'e, stopping criteria not fulfilled after 101 iterations tit which point the progress of the algorithm was extremely slow. Table 4 Statistics for test 9roblem 3, 5 × 5 nodes and n, = 72 Property

PIPA'

PIPAA

IMPA

IMPAA

MMA

# eq # subs # line inin r avg 7" max p, I MN )

101 101 l) 1.5X 10 " 0.0333

33 33 0 0.0315 0.353

7(t 7 50 7.4× 10 " 0.228

8 5 0 I I

7 4 -

2.500(X)0

2.5(X)000

1.697981

1.666666

2.50(!(100

:' Failure, stopping criteria not fulfilled after 101 iterations at which point the progress ~1" the algorithm was extremely slow. Table 5 Slatislics for lest problem 3. t0 × 10 nodes n = 342 Property

PIPAA

IMPAA

MMA

# eq stlbs .,~ line nlin 7:.lVg 7" lilax /~ (MN ~

39 39

I(I 4 I) I I

9 4 -

0 0.0209 0.256

2.5(X1000

5.000000

5.00001R)0

s h o w the s u p e r i o r i t y o f P I P A A a n d 1 M P A A o v e r P I P A a n d I M P A .

1 1.3. 7-est l,'ohh'/n .7 T h e s t r u c t u r e is i l h i s t r a t o d in Fig. 4 M

the c a s e 5 b y 5 n o d e s . T h e o b j e c t i v e is in this c a s e the n l a x i m u m

contact force b e t w e en truss and obstacle. The maxim um

a l l o w e d v o h i l n e o f the t r u s s is 1.5 t i m e s that o f the

initial d e s i g n : s~' = 0 . 0 2 in-', i i = 1. . . . .

n

1. . . . . J t a n d the g a p s b e t w e e n t r u s s and o b s t a c l e s are g~ = 0.001 l n ( i - 1) 15, ~l T h e m a t r i c e s Qk a n d Q w e r e c h o s e n to 1 × 10 ~"diagt 1) Ior all a l g o r i t h m s . E x e c u t i o n s t a t i s t i c s

is g i v e n f o r t w o d i f f e r e n t s i z e s o f the s t r u c t u r e in T a b l e s 4 a n d 5. T h e resulb; s h o w that P I P A a n d P I P A A yield b e l i e r o p t i m a ~hail I M P A . I M P A A a n d M M A .

12. I m p l e m e n t a t i o n

Proper scaling of fhc problem is itnportant for all the algoritlm>, t,~ a'~tfid mHnorical difficulties anti to obtain fast c o m c r g c n c e fl~r PII~,\/P[I~AA. It seems that P I P A / P I P A A i~, !v,orc scil~ilk.e to scaling than [ M P A / I M P A A and MKL,"~. By ex,~erimcnt ~,e d e d t i c e d ih~,t p r o p e r ~,calin.,_, n;c~,~, fi:m file ~,,riables u. p a n d J h a v e s i m i l a r

il/agnil;ld¢. !:or t!'lv c~m'-ide!cd i0st problenls we achieved il~ix i > u!,r',_, kiN; :> force ullit and nun its length unit in lhc ,:alctilatlll!!-. All h/il'llellleni;.itiolts where done in Miltl,;!b 5.1 rtlnniii~ ,,)i! ~.l~l 60 M H z He,elksialion eqtiipped with H P - U X IO.](I and 4','~ Mt~ R A M

D. Hilding et al. / Comput. Methods Appl. Mech. Engrg. 177 (1999)215-234

230

12.1. Solving the state equations In both IMPA/IPAA and MMA the state problem (8) must be solved to evaluate (u(.), p(. )) for a particular design s. MLCP(s) is reformulated as the following non-smooth equations:

H(s, u(s), p(s)) = 0,

_ F K(s)u + C:,P - F] H(s, u, p) - [ min(g - C,,u, p) J"

(25)

This system is then solved using the BN-algorithm, [2], which is a damped Newton method for B-differentiable system of equations. The main operation in the BN-algorithm is to solve the linear system

Vl,,.mH d = - H of n, + n equations, which is done by factoring the coefficient matrix V,,.mH. The BN-algorithm proved efficient when used together with IMPA/IMPAA and MMA. By using the values of (u(s~-~), p(s~-~)), from a previous design iterate, as starting point when solving (25) for (u(s~), p(s~)), the BN-algorithm converged typically after only 1 or 2 iterations. (This means that the V~,q,~H matrix was factored 1 respectively 2 times.)

12.2. hnplementation of MMA A MatLab code for MMA was provided by Krister Svanberg. What was needed to run this code for the min max unilateral problem was the sensitivity analysis; MMA needs gradients of objective function and constraints. The only nontrivial gradients are those of p, which were obtained using an adjoint method: If H is differentiable at s then by differentiating the equation (25) it is possible to obtain the following expression for V p,(s): .

,I-~,.(~)7

(s) = - a l l

0

J'

V,,,,,HA i =P,,

(26)

where P~ is an n,, + n zero column except for a one on row n,, + i. The evaluation of V p~(s) is very cheap as the factorization of V,,.mH already is available from the solution of the state problem. At points where H is not differentiable an arbitrary directional Jacobian of H is used instead of V,,.mH. Other approaches are possible for handling nondifferentiable points. However, H is non-differentiable in a very small set and it is thus unlikely that an iterate s ~ would be in this set. To conclude we believe that the handling of these points is of minor importance for the overall performance of MMA.

12.3. bnplementation of IMPA/IMPAA The performance of the algorithm seems not to be very sensitive to the choice of line search parameters and we set them to the following values: p = 0.5, o - = 0.2 and the minimum allowed step length c = 1 × 10 -~°. The main work in each each iteration, besides solving (25), is to solve the subprogram (15). In the implementation we used the equality constraints in the subprogram to solve for dp, thus producing a new reduced subprogram. If the set /3 is empty then at an iterate k dp~ = (V~p~(s*)) t ds,

i = 1. . . . . n ,

where V~p i is given by the expression (26). The reduced subprogram in ds only, was solved using MatLab's qp routine. If/3 is not empty then the relation between dpi and ds becomes more involved. It is, however, unlikely that/3 is non-empty and during the computational testing this never occurred. That case is therefore left without consideration.

12.4. Implementation of PIPA/PIPAA The performance of the PIPA/PIPAA proved sensitive mainly to: (i) the scaling of the problem; (ii) the update rule for the parameter o-h; (iii) the choice of the Q~ matrix; and (iv) the starting point for the variables (u, p, v). The performance seems less dependent on the line search parameters. We used the following values: ~ j = 1 with update rule trk =0.95o-x__~, initial penalty parameter ce_j = 1.01, centrality parameter p = 0 . 2 , c = 1 × 10 ~- and in the line search y = 0.5 and 6 = 0.95. The starting point for (u, p, v) was chosen to the vector 1 of all ones.

D. Hilding et al. / Comput. Methods Appl. Mech. Engrg. 177 (1999) 215-234

231

The main work in each iteration is to solve the subprogram (18). In the implementation we used the equality constraints in the subprogram to solve for (du, dp, dr) resulting in a reduced subprogram in ds only, which was solved using MatLab's qp routine. The reduced subprogram in ds and the subsequent calculation of (du, dp, dr) was, by an appropriate ordering of the calculations, obtained by factoring the (n, + n ) X (n, + n ) matrix

K(s) _ diag(pk)C,,

1

C;, diag(vk)j,

(27)

and no more than n + 1 forward and backward substitutions.

13. Conclusions We have extended two methods for solution of design problems with complementarity constraints. For comparison we have also used the MMA, which is well known and well tested in structural optimization. The extensions, as compared to the presentation in [9], is that a non-smooth objective function of max-type is used. This opens the possibility to effectively and correctly treat other problems of max-type such as those involving maximum effective stress. The numerical results show that the implicit method (IMPA) (in modified form, see below) and MMA were comparably effective and generally convergent in about the same number of iterations. Although MMA does not take into account the possibility of degenerate contact, this has not caused any problems during execution. Since MMA and IMPA are based on very similar ideas the overall comparable behavior may be expected. The PIPA generally required more iterations than the other two algorithms. However, in one case a significantly lower objective value was achieved by PIPA. The intuitive explanation for this we believe is that PIPA 'feels' nonactive unilateral constraints during iterations and therefore can take advantage of these constraints to create lower optima. The implicit methods use more local approximations and therefore have a tendency to get stuck in local minima. Finally, one of the examples showed the necessity to extend the basic formulations of IMPA and PIPA. The modification makes the subproblems less local. The MMA essentially has this modification in its original form.

Acknowledgements The authors thank Krister Svanberg for the use of the MMA MatLab code.

Appendix A PROOF OF PROPOSITION 3. Continuing the analysis in Section 6, we deduce that the following implication must hold: [Aid <~ 0

V i E ~¢*] ~

[max Pi' (s * ; d) >i 0]. Li~.~*

(A.1)

To avoid triviality, we may assume that p* ~ 0 ; thus 5¢*C a(z*). By the characterization of the directional derivative p'(s*; d) as described in Proposition 3, we may rephrase the implication (A.1) as saying that there exists no tuple (d, v, r) such that

Aid<~O,

i~s4[*

II s

diKiu* + K(s*)v + i= I

O= Ci,,v, O>r~

~ i~a(z*)Ufl(:*)

i ~ a(z*)

i E 5 t*

Cti, ri = 0

D. HiMing et al. I Compta. Methods Appl. Mech. Engt\~. 177 (1999) 215-_934

232

0 >1C~,,v ] O~
~

iEfi(z*)

riCi,,v : 0

O = r i,

i~3/(=*).

In turn, this is equivalent to saying htat for every complementary subsets [3~ and [32 of fi(z*), there exists no tuple (d, v, r) such that

Aid<~O, E i

iE~t*

diKiu*+K(s*)v+

E

I

O=Ci,,v, O>r,

C;,,ri=O

i ~ ~r(=;: ) U ,(3 I

iEa(z*)U~l

iE#*

O>~C,,,v, i E f i 2 O<~ri,

i~j

O = r i,

i~y(z*)Ufi2.

By dualizing the above system for all partitionings of the degenerate index set fi(z*), we obtain the desired system asserted in the proposition. []

Appendix B

Some details of The MMA. The Method of Moving Asymptotes is described in Svanberg I11,12]. Consider the lbllowing nonlinear program: minimize

f)(x)

subject to

f,(x)~
i= I...m,

(B.I)

Xlllin ~ .~ ~ ~._lll;IX

where each f : R" --->R is a differentiable function, x @ R" is the variable of the problem and .[lllill arid .It"1 . . . . . ~ R II are constants. (The variables introduced in this subsection are local and not related to those in the rest of the article.) B . / A COlll'#.~

apl~roxJmation

The heart of M M A

is the f o l l o w i n g

convex

L

approximation of ( B . I ) near a given poin! .~ .

iii

t

subject to

I

ff(.v) <~ O, I~ IllJ n

~

i = 1... m,

(B.2)

~ lllll\

v ~(), ~[,

~,,,here v C R'" are so called artificial variables, x E R", Mi some large numbers and lhe ,[', functions are con,,ex approximations of tile f defined as follows:

D. Hilding et al. I Comput. Methods Appl. Mech. Engrg. 177 (1999)215-234 .

=

i ~2_, i

(Ujk - x,k )2 . k

¢ .i~l

(i)

)

,k,2

(B.3)

( txi-Lj' \ (xi-Li)

XS

\

(<

" 't,t 7 - { 7 -s7- : s)

s

233

k-L

)'

where I(i) =/./:~x/(X ) k

and

=

~xi(X ) >

(B.4)

k

Here, U ] and L.! are so called asymptotes that determine how convex the program is and they must satisfy k , k Lk t Li, _< ~ X/ ~ U j . The closer U i, / are to Xi the more conservative is the approximation (B.2) of (B. 1). To prevent numerical ill-conditioning of the .~ approximation functions the bounds on x in (B.2) are set to be: . . i. ,O. Ixik + 0 . 9 L ~ ) , x....i i. =. max(x/ xi

=mln(xj

,0.1x + 0 . 9 U ),

(B.5)

j=l...n.

Due to the form of the program (B.2) it always has a solution 2 and can be solved efficiently by dual methods as described in Svanberg [ 10,11 ]. When it comes to the penalty factors M~, selecting too small M~ will lead to that (B.2) is not a good approximation of (B. 1) near x k, selecting Mi to large may give numerical difficulties. A good rule for selecting M~ is that M~ should be 'much' bigger than the Lagrange multiplier of the ith constraint of (B.2).

B.2. Description

of

the MMA

The M M A is an iterative method in which a sequence of subprograms (B.2) are solved: The M M A for program ( B . I ) Step 0 (hiitialization). Select a starting point x °, the penalty factors M i, asymptotes U~/, L~j) and set the iteration counter k = 0. Step 1 (New iterate). Let the new iterate x k+~ be the solution of the convex subprogram (30). Step 2 (Termination check). Check the prescribed rule for termination, if it is satisfied then x *+ ~ contains the final iterate else set k = k + 1, update the asymptotes to L~,U~ and go to step I. The asymptotes U ki and L jk are selected according to the following rules, see Svanberg [12]: k=0andk=l" L .ik = X ~ - C , ( X ~

Max

-Xi

Mill

)

(B.6)

Uk.i = X~ + Cl(xiIll - x ' i ~'i')

(B.7)

il x

k~>2:

c i =

,,

< - , - : ')(.,-: '-,-,

)<0

(.rj - .,,j*

)= 0

(B.8)

2/(1 +c2), L*=X~ - C ~ k-k_L k i ) i (.r, l i

:If we assume that for all i E {1 .., n]: .r:"" ~.v:'"'.

(B.9)

D. Hilding et al. / Comput. Methods Appl. Mech. Engrg. 177 (1999) 215-2,34

234 k

k

U i =X~ +

ck..

i(ui

k-I

-Xi

k-I

),

(B.10)

w h e r e w e u s e d C~ = 0.5 a n d C 2 = 0.7 in o u r n u m e r i c a l c a l c u l a t i o n s . S m a l l e r v a l u e s o f C~ a n d C : m a k e s t h e a l g o r i t h m m o r e r o b u s t b u t p o t e n t i a l l y a l s o s l o w e r . W e set all t h e p e n a l t y f a c t o r s M, = 1000 a n d this w o r k e d w e l l in o u r c a s e .

References [1] R.L. Benedict and J.E. Taylor, Optimal design for elastic bodies in contact, in: E.J. Haug and J. Cea, eds., Optimization of Distributed Parameter Structures (Sijthoff and Noordhoff, Alphen aan den Rijn, Holland, 1981) 917-930. [2] P.W. Christensen, A. Klarbring, J.S. Pang and N. Strtimberg. Formulation and comparison of algorithms for frictional contact problems, Int. J. Numer. Methods Engrg. 42 (1998) 145-173. [3] T.F. Conry and A. Seireg, A mathematical programming method for design of elastic bodies in contact, ASME J. Appl. Mech. 38(2) (1971) 387-392. [4] E.J. Haug and B.M. Kwak, Contact stress minimization by contour design. Int. J. Numer. Methods Engrg. 12 (1978) 917-930. [5] D. Hilding, A. Klarbring and J. Petersson. Optimization of structures in unilateral contact. Technical Report LiTH-IKP-R-996, Link6ping University, IKP, 1998, to appear in Appl. Mech. Rev. [6] A. Klarbring, On the problem of optimizing contact force distributions. J. Opt. Th. App. 74 (1992) 131-150. [7] A. Klarbring, J. Haslinger. On almost constant contact stress distributions by shape optimization. Struct. Opt. 5 (1993) 213-216. [8] A. Klarbring and M. R6nnqvist. Nested approach to structural optimization in nonsmooth mechanics. Struct. Opt. 10 (1995) 79-86. [9] Z.-Q. Luo, J.-S. Pang and D. Ralph, Mathematical Programs with Equilibrium Constraints (Cambridge University Press, 1997). [10] K. Svanberg, An algorithm for optimum structural design using duality. Math. Prog. Study. 20 ( 1981 ) 161-177. [11] K. Svanberg, The method of moving asymptotes - - a new method for structural optimization. Int. J. Numer. Methods Engrg. 24 (1987) 359-373. [12] K. Svanberg, A globally convergent version of MMA without linesearch, in: N. Olhoff and G.I.N. Rozvany, eds., Proc. of the 1995 1st World Congress of Structural and Multidisciplinary Optimization, ISSMO (Elsevier Science Ltd., Oxford, 1995) 9-16.