- Email: [email protected]
Chapter 14 Numerical & Biological Shape Optimization
Unification of Finite Element Methods H. Kardestuncer (Editor) 0 Elsevier Science Publishers B.V. (North-Holland), 1984
CHAPTER 14 NUMERICAL & BIOLOGICAL SHAPE OPTIMIZATION 4. Philpott & G. Strung
This paper describes a shape optimization code which moves the joints of the structure toward optimal positions and handles the possibilfty of degeneracy in the underlying linear program. The code i s applied to a design problem arising in biology, attempting to analyze the framework of internal fibers in the patella and eventually to predict their pattern of remodeling after fracture. The underlying problems also appear in engineering design applications, involving optimal bounds on the effective properties of composite materials and nonconvex functionals in the calculus of variations. INTRODUCTION We report here on an optimization problem which might have arisen in engineering design, but actually came from biology. under given constraints, to find It is a classical problem a structure of minimum weight. We take the structure to be a plane truss subject to forces along its boundary. The forces will roughly approximate those to which the human patella is actually subjected, and we reduce to two dimensions (plane strain) by considering a vertical cross-section taken from front to back in the knee. Figure 1 indicates the tensile forces on the upper and lower edges of the patella, from the quadriceps and the ligaments, and the equilibrating compressive force from the support provided by the femur.
-
When these forces are concentrated at single points as in Figure 2, Michell [l] found the solution. It is a truss in which the members under compression are orthogonal t o those
32 1
322
A . Pliilpott & G. Strutig
i n tension, as i s always the case i n an optimal design. "he precise geometry depends on the points a t which the forces a c t , but t h e design i t s e l f i s straightforward: t h e compress i o n members a r e d i s t r i b u t e d r a d i a l l y and t h e tension is concentrated on a s i n g l e circumferential member connecting t h e two forces. Certainly t h a t is not an accurate description of t h e fibers i n s i d e t h e p a t e l l a . On t h e other hand it is not absurd; both s t e r e o l o g i c a l d a t a and f i n i t e element analyses [2-41 confirm t h i s general p a t t e r n . Photographs of t h e trabeculae show a framework too s t r i k i n g t o be an accident. Already i n t h e l a s t century c l i n i c a l observations led t o the formulation of Wolff I s law, t h a t t h e f i b e r s a r e i n the d i r e c t i o n s of p r i n c i p a l s t r e s s . A more exact statement of t h a t law, and of the mechanisms thatgovern t h e laying down and resorbing of bone, i s badly needed. It i s a challenge t h a t s e v e r a l orthopaedic research groups have accepted, and we thank W. C. Hayes and h i s colleagues a t Beth I s r a e l Hospital f o r p a t i e n t encouragement t o explore t h e r e l a t i o n s h i p t o optimal design. A f i r s t s t e p (our immediate goal) I s t o r e f i n e Figure 2 by admitting a more r e a l i s t i c s e t of loads. This paper discusses t h e numerical optimization and not t h e biology. Even with t h i s l i l n i t a t i o n i t is incomplete: we could not pursue a l l possible approaches. Some t h a t we d i d pursue were not successful o r r a t h e r , t h e i r c o s t increased so rapidly t h a t a r e a l i s t i c three-dimensional p i c t u r e (with t h e bone allowed t o remodel a f t e r f r a c t u r e ) was inaccessible. Other algorithms a r e more promlslng. Perhaps we can l i s t here some of the a l t e r n a t i v e s , with preliminary comments:
-
1. Fixed Geometrg. Fix t h e nodes of t h e t r u s s and allow a specified subset of edges as candidates i n t h e design optimization. This y i e l d s a standard problem i n l i n e a r programming, with a sparse c o e f f i c i e n t matrix described below. The c o s t i s e f f e c t i v e l y t h e sum C a 1s I of edge lengths multiplied &Ill by member f o r c e s . The.simplex method w i l l y i e l d an optimal (and s t a t i c a l l y determinate) t r u s s .
Numerical & Biological Shape Optimizatiotz
323
It was hoped that the simplicity of t h i s model would allow a comparatively large number of nodes and edges, g i v i n g the design enough degrees of freedom t o imitate curved f i b e r s . We now believe t h a t t h i s approach i s not e f f i c i e n t . The approximation of a nonlinear design problem by l i n e a r programming can be very slowly convergent. In our experience the l i n e a r programs themselves were a l s o highly degenerate, w i t h multiple pivot steps before the cost was reduced.
Variable Geometry. The positions of the nodes a r e optimized a s well as the weight of the t r u s s t h a t connects them. Normally these steps are separate; a l i n e a r program computes the t r u s s f o r fixed nodes, and a nonlinear calculat i o n moves the nodes. O f the codes we know about (and there must be many othera) we mention three: 2.
( I ) The CONMIN code developed by (3. Vanderplaats [ S ] : This was generously made available t o us and was the basis f o r the conference paper by Snyder e t a 1 [ 6 ] on trabecular orientation i n the p a t e l l a . The code admits a much wider range of problem i n 2 and 3 dimensions, w i t h constraints on displacements and buckling s t i f f n e s s and frequency response. (11) The Michell structure code developed by McConnel [7]: This uses a d i r e c t search method of Powell t o optimize the nodal positions (with fixed node-edge incidence matrix)
.
Experimental r e s u l t s on problems of bridge design studied by Hemp [see 81 look very s a t i s f a c t o r y and we hope t o t e s t the code on biological structures. (iii) A newly developed code:
We are conpletlng a s p e c i a l purpose code t o t e s t other algorithms f o r this c l a s s of optimization problems. This paper describes the use of the dual variables from the fixed-node program t o choose the directions of movement of the nodes. It seems c l e a r t h a t t h i s variation of geometry i s valuable and probably e s s e n t i a l .
3. Michell Displacement Field. It was the great discovery of Michell that under the presribed load, the Optimsl design
324
A . Philpott & G. Strang
has p r i n c i p a l s t r a i n s of constant magnitude. ( I n a region where no members are required i n one o r both directions, t h e corresponding s t r a i n s a r e uniformly below t h a t constant value.) Lagache [ 9 ] has found a b e a u t i f u l existence proof f o r Michell t r u s s e s based on a convergent sequence of approximations t o t h i s s t r a i n f i e l d . Although those p a r t i c u l a r approximations nay be i n e f f i c i e n t i n computations they a r e based on fixed nodes, as i n 1. above t h e r e seems t o be no reason why a program t o compute f i n i t e element dlsplacements f o r Michell's s t r a i n f i e l d should not succeed.
-
-
Michell s t r e s s f i e l d . I n the s t r e s s formulation of t h e problem, t h e c o n s t r a i n t d i v cf = 0 can be removed by deriving t h e sOress u from Airy's function t(x,y): 4.
same
noted i n [lo], t h i s y i e l d s an ordinary (but nonquadratlc) problem in t h e calculus of v a r i a t i o n s . Again a combination of f i n i t e elements and nonlinear optimization would give a d i r e c t approximation of the Michell d i r e c t i o n s in t h e optimal design. The d e n s i t i e s and o r i e n t a t i o n s then come d i r e c t l y from 6 r a t h e r than i n d i r e c t l y from a d i s c r e t e t r u s s . As
This paper describes t h e new v a r i a b l e geometry code mentioned above. It was prepared i n l a t e 1983 and t e s t e d on t h e
-
b i o l o g i c a l applications although it applies equally t o engineering design. The next s e c t i o n of t h e paper r e c a l l s t h e s t e p s of the simplex method, modified only t o account f o r the absolute values i n the c o s t function; members can be i n tension or compression. Section 3 applies the simplex method t o c o s t ( o r volume) minimization, f o r a given placement of the nodes. Section 4 develops an algorithm f o r moving t h e nodes toward an improved design. There a r e numerical d i f f i c u l t i e s associated w i t h degeneracy, handled by a proj e c t i o n technique. In t h e l a s t s e c t i o n we b r i e f l y discuss the implementation and the applications.
Niitvierical& Biological Shape Optimization
325
THE FIXED GEOMETRY DESIGN PROBLEM
number of authors,startimg w i t h Dorn, Greenberg, and Gomory [U],have considered the application of l i n e a r programming t o the minimum weight problem f o r a plane t r u s s subjected t o forces along i t s boundary. "he problem is t o s e l e c t a s e t of members connecting a given collection of pin j o i n t s which a r e fixed i n number and position. Each member i n the optimal t r u s s is i n tension or compression a t i t s corresponding l i m i t s t r e s s . A
The formulation as a l i n e a r program is as follows. We consider m pin j o i n t s I n the plane, and l e t XiDyi, i = 1, .,m be the coordinates (with respect t o some o r i g i n ) of the l t h j o i n t . We l e t [ i , k ] denote the member connecting j o i n t s 1 and k . To avoid unnecessary subscripts we associa t e with each [ i , k ] a unique index j E (1,2,. .,n), and as the force i n member j . Here sj is positive define s j i f j is in tension and negative i f j is i n compression. j in "he number of members n can be a t most - , practice we r e s t r i c t our choice of members t o a r e l a t i v e l y small set, usually those having length below some chosen bound.
--
The forces must be such t h a t each j o i n t is in s t a t i c "j equilibrium. These 2 m equality constraints may be formulated by resolving each member force in the two coordinate directions and summing the forces i n each direction a t each node. The right hand sides of these constrainta are the x and y components of the applied forces. For j o i n t I, these a r e denoted by fli and fgi respectively. For many j o i n t s , including those i n t e r n a l t o the structure, the applied forces are zero. If member j = [ i , k ] makes en angle 0 w i t h the positive j x-axis a t j o i n t i, we may write the constraints f o r node i as
A . Plrilpotr & G. Strarig
326
where Ii is the set of members connected to the ith joint. Since each member connects exactly two joints, these equilibrium constraints in matrix form become
where 8 is the vector of n member forces, fl and f2 are the x and y components of the boundary forces, and A1 and A2 are m X n matrices with zero entries in column j = [ i , k ] except in row8 i and k where
-
(A21ij = -sin
6
3 ’
(A21kj = s i n 0j
I;:[
-
In what fOllOW8 we shall refer to the matrix and the right hand side vector
as
f.
A n example of a pinjolnted trusla and its constraint matrix
is given in Figure 3.
If
aj
is the length of member j, minimizing the volume of a planar truss is equivalent to minimizing
n over those vectors 8 which satisfy the equilibrlun condition As = f. The derivation of this result requires that the l i m i t stresses in tension (a=) be the same for
327
Nu 17 i er icul & Biological Sliap e Opt irri izat io n
each member, and t h a t the l i m i t s t r e s s e s under compression (udn) be the same f o r each member. Since W is independent of the values of uand %in , we may assume without 108s of generality that u- udn = u. When each member j i s f u l l y stressed i t s cross-sectional area is proportional t o
-
9:
t!i
aJ
= u
, for
j
i n tension,
for
j
i n compression.
= -6,
1s I measures the thickness required in member 3 j j the problem of minimizing the weight of a fixed j o i n t planar t r u s s may thus be formulated as
In this way
Problem P:
n
C djls,l j=l
Minimize
subject t o As = f
.
(1)
As posed, t h i s objective function is nonlinear. However it is easy t o reformulate P as a l i n e a r program by defining
variables
whence W become
C 1 t
J j A(t
+
Z
-
1u
U)
j j =
is l i n e a r and the constraints
f, t
2
0,
U
2
0
.
We s h a l l write the problem in the formulation P, but make considerable use of the l i n e a r programming formulation and the large body of theory which t h i s allows us. W e can formulate the problem dual t o D: Maximize
vTf l
+
P
as follows:
wTf2 subject t o I(vTA,+wTA2),I j =
1 , e . m
an
5 dj, (21
A. Philpott & G. Strang
328
Multiplying both s i d e s by I S il duclng f l = A1s and f 2 = A2s easy half of d u a l i t y theory.
I j
on jj and i n t r 0 leads t o weak d u a l i t y the
-
SUmming
Theorem 1: If v and w a r e f e a s i b l e f o r t h e dual and s i s f e a s i b l e f o r the primal then n T T (3) v f l + w f2 5 c i j l S j l
J-1
The complementary slackness conditions, which must hold a t the minimizing 6 and maximizing v and w, a r e as follows : if
sj
>
0, then
if
8
<
0, then
j
( V TA
1
+ wTA p ) j
aj (5)
The dual problem i s o f t e n given the following physical l n t e r pretation, with E a s Young's modulus. A s e t of dual variables v i J w i J corresponds t o a v i r t u a l deformation of the t r u s s I n which each j o i n t 1 i s displaced by v i ( i ) i n the y d i r e c t i o n . The elongation e j of each member j under such a deformation i s given by
The dual c o n s t r a i n t s t h e r e f o r e require of each elongation e j that
I n other words t h e l i n e a r s t r a i n i n each member caused by t h e v i r t u a l deformation must be a t most t h e value which corresponds t o the member being a t i t s limit s t r e s s . The dual objective function vTf l
+
wTf 2 can be i n t e r p r e t e d
as E/a times t h e work done by e x t e r n a l f o r c e s t o produce t h i s v i r t u a l deformation. Since t h e optimal s o l u t i o n s w i l l s a t i s f y (4) and ( 5 ) J t h e s t r a i n in a member with sj o produced by t h e optimal displacements v and w must be
+
Numerical & Biological Shape Optimization
329
the strain which would occur with j at its limit stress (a if sj > 0 and - u if sj < 0). Members with s = 0 J and they are below do not appear in the optimal structure the limit stress when subjected to the optimal v and w.
-
AN ALUORITHM FOR THE FIXED GEOMETRY PROBUM
Since the fixed geometry design problem P may be formulated as a linear program, the simplex algorithm will solve it. This method constructs a sequence of basic feasible solutions with decreasing weight C l? I s I . The optimum does occur at J 3 a basic feasible solution (defined below) and the sequence terminates at an optimal solution after a finite number of steps. For a planar truss with n joints we have 2m equilibrium constraints. However the rows of A are not linearly independent, and for equilibrium the boundary loads must sum to zero in both coordinate directions. Furthermore the sum of moments about any point in the truss must be zero. The number of independent equilibrium constraints is therefore 2m-3, and three constraints must be removed before the simplex algorithm can be applied.
In practice these degrees of freedom are removed by fixing two or more boundary Joints. Then we remove from A1 and the rows which correspond to the anchored joints, and we remove from the vectors v and w the displacements which are no longer possible. We shall assume in what follows that the constraint matrix A has undergone this process, eqaal to the and take m to be the number of rows of A number of displacement unknowns. This removes at the 8ame t h e any restrictions on the boundary loads; they need not be self-equilibrating since any additional forces necessary for equilibrium will be taken up by the foundation. A2
-
The basic fersible solutions for problem P are found by choosing 2m linearly independent columns of A and solving (1) uniquely for the member forces sj corresponding
A. Philpott & G. Strang
330
t o these columns. A l l other member forces a r e zero. Thus each basic f e a s i b l e s o l u t i o n corresponds t o a e t a t i c a l l y deteranate truss. W e s h a l l c a l l any such truss a basic f e a s i b l e t r u s s . S t r i c t l y speaking, t h e basic f e a s i b l e s o l u t i o n t o t h e l i n e a r program is given by t h e nonnegative u : variables s a t i s f y i n g s = t
-
and With this understanding we apply t h e term "basic f e a s i b l e solution" t o the member forces 8 and also t o t h e variables
t
and u.
If t h e choeen columns of
a r e not l i n e a r l y independent, t h e corresponding truss i s no longer r i g i d under an a r b i t r a r y set of boundary loads, and it becomes a mechanism. On the other hand, i f more than 2m columns a r e chosen (with 2 m l i n e a r l y independent) then t h e corresponding t r u s s i s s t a t i c a l l y Indeterminate, and we cannot solve uniquely f o r the member f o r c e s . A t h i r d p o s s i b i l i t y is t h a t t h e 2m columns a r e l i n e a r l y independent and y e t some of t h e member forces in t h e basic f e a s i b l e truss a r e zero. I n t h i s case t h e basic f e a s i b l e s o l u t i o n i s degenerate. A
The algorithm f o r the f i x e d geometry optimal design follows t h e s t e p s of t h e slmplex method, modified to work w i t h t h e formulation P. If we assume that degenerate basic f e a s i b l e s o l u t i o n s do not occur, t h e s t e p s a r e as follows:
[i;]
1. Determine a s e t of 2 m members which form a basic f e a s i b l e t r u s s . Denote by B = t h e 2m by 2m
"basis matrix" formed by t h e corresponding columns of 2.
A.
Solve Bs = f .
3. For t h e s o l u t i o n S, construct a complementary s l a c k dual s o l u t i o n v#w. I n o t h e r words, solve t h e l i n e a r system (vTB1
+
wTB )
2 j
=
-+ 4j
(choosing t h e s i g n of
Sj)
33 1
Numerical & Biological Shape Optimization
4.
Choose the value of
'J 5. If
2
=
aJ
-
j (say
j = q)
I(vTAl + wTA )
2J
I
which minimizes
.
2 0 then stop3 the current t r u s s is optimal. Q ra < 0 then member q enters the basic f e a s i b l e t r u s s in t'ension o r compression depending on the sign of T T = ( V A1 + w If
r
9
6.
Determine the member t o leave the basic f e a s i b l e truss, over given by the Index 1 = p which minimizes l q i l indices 1 i n I. Here P
si (B'lA
iq
I = ( I ; (B-lA)iq
0
and sign
=
sign zs)
.
7.
Replace member p i n the basic f e a s i b l e truss (and the column p in the matrix B) w i t h member q . Return t o Step 2.
It is easy t o see t h a t t h i s algorithm follows exactly the same steps as the simplex method. Step 4 f i n d s a m i n i m u m reduced cost, and ti i n Step 6 i s the r a t i o t e s t t h a t determines which basic variable becomes non-basic a f t e r a pivot operation. I n t h i s step, i f no quotients of approp r i a t e sign e x i s t then the problem has an unbounded solut ion.
Since the constraint matrix A i s very sparse, the implementation of the algorithm uses an LU f a c t o r i z a t i o n of the basis matrix B and updates the triangular f a c t o r s L and U a t each i t e r a t i o n . We employ the method of Forrest and Tomlin [12]. To take f u r t h e r advantage of sparsity, the columns of A corresponding t o nonbasic variables a r e recomputed a t Step 4 i n each i t e r a t i o n . The steps described above are followed w i t h no e s s e n t i a l changes when degeneracy is encountered. If the s t a r t i n g
332
A. Philpott & G. Strarig
basic feasible t r u s s i s degenerate we a r b i t r a r i l y assign positive or negative signs t o the zero forces. Since by Step 5 every member entering the basis does so e i t h e r under tension or compression, f o r every subsequent basic feasible t r u s s we can assign positive and negative signs t o the (possibly zero) member forces. This convention gives a unique complementary slack dual solution i n Step 3. The choice of member p t o leave the basic feasible t r u s s is however not uniquely determined, since any 1 w i t h si 0 w i l l render l$il a minimum i n Step 6. In the code we choose the variable with lowest index, and although t h i s rule does not preclude cycling, no instances of cycling have been observed. (Dantzig has observed only one i n long experience with simplex calculations. ) Of course a more expensive pivot rule could ensure t h a t cycling never occurs.
-
VARIABLE GEOMETRY
For fixed j o i n t positions the above algorithm gives an optimal solution it finds the minimum weight of a plane t r u s s which supports given boundary loads. Being e s s e n t i a l l y the simplex method, the algorithm terminates a t the optimum a f t e r a f i n i t e number of I t e r a t i o n s . Unfortunately, f o r large n and m (many j o i n t s and many possible members) the time taken t o reach the optimum is prohibitive. If members are permitted between any two j o i n t s then the problem s i z e Increases rapidly with the number of j o i n t s . To avoid t h i s problem we may r e s t r i c t the allowed members t o those whose length l i e s below some fixed bound. However, even a small problem w i t h 100 j o i n t s and 800 members can take hours of CPU time
-
.
One reason f o r t h i s poor convergence behavior i s the high degree of degeneracy inherent i n the l i n e a r programs associated w i t h pin jointed structures. Often a large number of j o i n t s have no external forces, and play no p a r t i n the designs the members incident upon these j o i n t s carry zero force. Even If these j o i n t s and the members connected t o them are deleted,
333
Numerical & Biological Shape Optimization
it is often the case i n a basic f e a s i b l e truss that other zero-force members e x i s t . As observed i n [ l l ] , these cannot be deleted since the resulting truss would no longer be r i g i d . The presence of degenerate basic f e a s i b l e solutions causes poor performance in the simplex algorithm, since many pivot operations occur w i t h no change i n the objective f unc t ion. To approximate accurately the trabecular architecture of the human patella, a large number of i n t e r n a l j o i n t s would be
-
i f t h e i r positions a r e fixed. needed i n the design process This p r o l i f e r a t i o n of j o i n t s subject t o zero loads exacerbates the degeneracy problem, and the optimization becomes Idq)ract i c a b l e . A n a l t e r n a t i v e is t o admit a much smaller number of j o i n t s and then alter t h e i r positions i n such a way a8 t o improve the design. We proceed t o describe a method f o r computing the gradient of the cost W in terms of the j o i n t positions, thus giving a s e t of directions of steepest descent; s h i f t i n g the j o i n t s i n these directions produces a truss w i t h l e s s weight.
To f i n d the gradient, we investigate the e f f e c t on the objective function W of a change i n position of the j o i n t s from (Xisyi) t o ( x i + 8Xj.s y i + 8yi). A t optimal solution of problem P# the cost I s n w = jc-1 aj l a3 I .
The sum is taken over members or" the optimal basic f e a s i b l e t r u s s . L e t this truss correspond t o the basis matrix B =
LB,'] 1 rB,
where
B1
and
B2
are
m x 2 m matrices.
We
assume-for the present that this truss is not degenerate. Changing the positions of the j o i n t s w i l l change the length8 a j and a l s o a f f e c t the e n t r i e s of B implying that the member forces will a l s o change. Suppose that the forces become s j + 8 s j 8 and the changes in angles a l t e r B t o B + bB. Then equation (1) gives
-
A. Philpott & G. Strang
334
(B
+
+
bB)(s
6s) = f ,
whence t o f i r s t order
+
B6s
(6B)s = 0 .
The f i r s t order change i n W = C
a 3 1s3
+
C 0aj
a
8s
J J
given by Sjl
+
2 Jj61Sj(.
sj'o
The last term i n (7) is absent by the assumption of nondegeneracy. It then follows from the definition of the complementary slack dual variables (the displacements v w ) that
and thus from (6) we obtain
r
i
L
J
(7)
and
We now express &4? and 6B i n terms of the changes in j o i n t positions, Ox and by. If j = [ i , k ] , then the is (xk + (yk Y , ) ~ . Since length
$
2
-
-
=
COB
ej,
aa 4 --
COB 0
j '
& aa = - s i n e j , 4, Yi
2 = s i n e
it follows from the d e f i n i t i o n of 6a
-
T B16x
T + B26y.
B1
and
Be
that (9)
The matrix 8B is the f i r s t order change i n the matrix B and a r i s e s from the changes i n COB e and sin 8 f o r J 3 each member a f t e r the positions o f some j o i n t s are changed. It is easy t o ehow that if j = [ i , k ] ,
335
Numerical & Biological Shape Optimization
7 -
cos e
yi
s i n ej ,
a
cos eijsin ej
cos dyk e
j
Therefore the changes i n cosine and sine a r e given t o f i r s t order by
The change ( 6 B ) j i n column j of B I s zero except i n , k, i + m, and k + mj i n these four rows we have rows I
w i t h the same formulas and opposite signs f o r ( 6B)k+m, j . The term i n brackets is j u s t
If we l e t Sjj
=
2, j
8W
S
be the
then
6W
=I4T[ BTl .TI[(.]- 8 x
2m x 2m
(FJB),~
diagonal matrix w i t h
I n ( 8 ) can be written as follows:
[vT wT 1[82] s[B; OB1
We simplify t h i s expression by defining
-B3
E];
and
336
A. Philpott & G. Strang
leaving 6W = gT b x + hT b y .
(13)
For given loads and a given nondegenerate optimal truss, the vector [If] defined by (11) and (12) i s a direction of steepest descent f o r W when expressed as a function of the nodal positions. Therefore it i s natural t o extend the algorithm t o include a s h i f t of positions i n t h i s descent direction. This variable-geometry algorithm proceeds as follows : Step 1. For a given s e t of j o i n t coordinates the fixed geometry problem P
-x,y,-
solve
Step 2.
For the optimal basic feasible solution t o calculate g and h from (11) and (12)
Step 3.
For the given boundary loads, find the value of Xh) the stepsize 1 which minimizes W(z Xg,Y
Step 4.
If the change I n W i s l e s s than some tolerance then terminate; otherwise update 'ji and and return t o S t e p 1.
-
Pa
-
The c r u c i a l step is the l i n e search i n Step 31 It is here t h a t the Improvement In geometry I s made. There are a number of ways t o carry out t h i s search. One p o s s i b i l i t y I s t o r e t a i n the current basic feasible truss and recalculate W f o r this truss a t each point i n the l i n e search (determining t h e member forces from equilibrium Ba = f ) . Alternatively, we may return t o the simplex algorithm and find the optimal basic feasible truss a t each point of the search. We tend t o favor the former approach) re-solving the l i n e a r
331
Numerical & Biological Shape Optimization
programming problem a t each point requires a greater comput a t i o n a l e f f o r t , a s i t u a t i o n which is not improved by the presence of degeneracy. One might hope t h a t t h i s algorithm gives a reasonably robust method of improving the geometry. Unfortunately t h i s is not the case.
In practice, l i n e searches often terminate a t values of A f o r which the resulting basic f e a s i b l e t r u s s is degenerate. This seem reasonable since the cross-sectional area of the degenerate member has a l o c a l minimum exactly a t the value of X which renders it degenerate. This contradicts our assumpt i o n ofnondegeneratebasic f e a s i b l e t r u s s e s j even i f we begin the l i n e search a t a well-behaved truss, we a r e l i k e l y t o meet degeneracy. possible remedy is t o continue as i f our assumption were valid, and hope that the sum X 1 61s I over members with J J sj = 0 (which was neglected i n ( 7 ) ) is small. However it soon becomes obvious t h a t this approach is untenable. In practice, l i n e searches from current degenerate basic f e a s i b l e trusses often terminate a t a steplength A which is zero t o working accuracy, and the algorithm is j m e d a t a suboptimal solution. Clearly a t some point, t h i s neglected fourth term i n (7) has a considerable e f f e c t . A
The solution i s t o compute a more e f f e c t i v e search direction. Since the undesirable behavior ie caused by changes i n 1s I J f o r degenerate members, It makes sense t o seek a direction which minimizes these changes. We therefore search in a direction which, a t l e a s t locally, keeps degenerate members degenerate. Then one expects t h e i r contribution t o 8W t o be small. The adjusted search direction is computed as follows. we have 8s = -B-l(8B)s
,
From (6)
A. Philpott & G. Strang
338
whence, using (lo),
Let Id be the set of indices of columns of B corresponding to degenerate members, and let Bd be the matrix formed by the rows of Bol with these indices. Then if
It;] -
D = -Bd[:il]
S [BE -Bl], T any vector
_
I
in the nullspace of
D will render (6s) = 0 for j E Id. Thus to obtain a descent direction which keeps the changes in degenerate onto the members as smell as possible, we project nullspace of D :
[E]
[f] =
[g] - DT(DDT)-l
D[E]
.
Since DDT is often singular, we use a Gram-Schmidt procedure t o construct a matrix Q with orthonormal columns spanning the column space of DT, whence
Even with this projected search direction, we are not for large step sizes guaranteed a decrease in W(%XE,T-Ali) A. Indeed, if a step is so large that some member forces change sign, the current choice of direction may no longer give descent. For this reason, we compute from (14) the values of X for which a member force changes sign, and evaluate the objection function W at these values in succession. (The equilibrium equation is solved for 8.) When an increase in W is detected, this gives an interval in X over which a golden section line search minimizes W. C ONC LUS ION
The algorithm described above was implemented in FORTRAN on a VAX 780, and applied to a model of the human patella. The loads correspond to a 60’ flexion of the leg (as in
Numerical & Biological Shape Optimization
339
climbing s t a i r s ) . We worked with 40 j o i n t s , and allowed a l l edges w i t h length not exceeding one centimeter. An i n i t i a l basic feasible t r u s s i s shown i n Figure 4a. The 14 fixed nodes along the bottom provide reactions t o the 5 applied loads. Solid l i n e s denote members i n tension and dashed l i n e s indicate compression. The thickness of the members i s represented by the number of l i n e s i n the diagram, and bars denoted by single l i n e s i n Figure 4a are degenerate. They have zero thickness, as i n the two members which meet a t the unloaded node a t the top right hand side; they a r e present only t o avoid a mechanism. The weight of the i n i t i a l t r u s s i s W = 2.064 grams. Figure 4b shows the output from the fixed geometry design algorithm, s t a r t i n g from the truss of Figure 4a. Degenerate members ,are now omitted. The optimal t r u s s has a weight of W = 1.829 gram and was obtained a f t e r 25 pivot operations i n the simplex method. For a problem of t h i s size, each pivot step takes approximately .5 seconds, a value which increases substantially f o r l a r g e r problems. Figure 4c shows the output from the variable geometry algorithm. It i s the o p t i m l basic f e a s i b l e t r u s s a f t e r a succession of changes i n nodal positions. The weight i s reduced t o W = 1.708 gram a f t e r t e n l i n e searches, a f t e r each of which the fixed geometry problem was resolved and new directions of movement were computed (with projection t o control degeneracy). Subsequent l i n e searches f a i l e d t o give a significant impuovement i n the weight of the trusa and the program terminated. ACKNOWLEDGEMENT We thank the Army Research Office and the National Science Foundation for their support under contracts We are also grateful to the National Institute for Health for the opportunity, under
DAAG 29-83K-0025 and 81-02371.
contract AM30875, to work in the Orthopaedics Research Laboratory at Beth Israel Hospital.
3 40
A . Plzilpott & G. Strang
REFERENCES [l] Mlchell, A.G.M., The limit of economy of material in frame structures, Phil. Mag-8 (1904) 589-597. [2] Hayes, W.C., Snyder, B., Levine, B.M. and Ramaswamy, S., Stress-morphology relationships in trabecular bone of the patella,in: R. H. Gallagher et al., eds., Finite Elements in Biomechanics ( J o h n Wiley, New York, 1983). [3] Huberti, H.H., Hayes, W.C., Stone, J . L . and Shybut, G.T., Force ratios In the quadriceps tendon and ligamentum patellae, J. Biomechanics, to appear. [ 41
Stone, J . L . , Three-dimensional stress-morphology relationships In trabecular bone, M.Sc. Thesis, M. I . T . (1983).
-
, Numerical Optimization Techniques for Engineering Design With Applications (McGraw-Hill, New York, 1984).
[ 5 1 Vanderplaats, G.
[61 Snyder, B., Strang, Go, Hayes, W.C. and Norris, G., Application of structural geometry optimiration techniques in the microstructural remodeling of trabecular bone, ASME Annual Meeting, Boston, 1983.
[TI MCCOMel, ROE., Least-weight frameworks for load across Spar
ASCE J.
H e w W.5.r 19731 191
m.MeCh.
DIv. 100 (1974) 885-901.
O;ptlmU Structures (Clarendon Press, Oxford,
Lagache, J.-M., Treillis do volume minimal dans une region donnh!, J. Mdcanique 20 (1981) 415-448.
[lo] Strang, G. and Kohn, R.V., Hencky-Prandtl nets and constrained Michell trusses, Comp. Methods in Appl. Mech. E W . 36 (1983) 207-222. 1111 Dorn, W.S., Gomory, R.E. and Greenb,erg, H . J . , Automatic design of optimal structures, J. Mecanique 3 (1964) 25-52.
1121 Forrest, J . J . H . and Tomlin, J.A., Updated triangular factors of the basis to maintain sparsit in the product form simplex method, Math. Prog. 2 (19727 263-278.
Numerical & Biological Shape Optimization
FIGURE 1. HUMAN PATELLA: ROTATED CROSS SECTION, WITH APPLIED TENSIONS AND SUPPORT REGION.
F' FIGURE 2. MICHELL TRUSS, 3-POINT LOAD.
341
A . Philpott & G. Strang
342
3
FIGURE 3. A SIMPLE PINJOINTED TRUSS.
scale:
1 cm.
FIGURE 4A. INITIAL TRUSS.
Numerical & Biological Shape Optimization
F I G U R E 4 B . OPTIMAL TRUSS WITH F I X E D NODES.
FIGURE 4C.
OPTIMAL TRUSS AFTER 10 GEOMETRY CHANGES.
3 43
CHAPTER 14 NUMERICAL & BIOLOGICAL SHAPE OPTIMIZATION 4. Philpott & G. Strung
This paper describes a shape optimization code which moves the joints of the structure toward optimal positions and handles the possibilfty of degeneracy in the underlying linear program. The code i s applied to a design problem arising in biology, attempting to analyze the framework of internal fibers in the patella and eventually to predict their pattern of remodeling after fracture. The underlying problems also appear in engineering design applications, involving optimal bounds on the effective properties of composite materials and nonconvex functionals in the calculus of variations. INTRODUCTION We report here on an optimization problem which might have arisen in engineering design, but actually came from biology. under given constraints, to find It is a classical problem a structure of minimum weight. We take the structure to be a plane truss subject to forces along its boundary. The forces will roughly approximate those to which the human patella is actually subjected, and we reduce to two dimensions (plane strain) by considering a vertical cross-section taken from front to back in the knee. Figure 1 indicates the tensile forces on the upper and lower edges of the patella, from the quadriceps and the ligaments, and the equilibrating compressive force from the support provided by the femur.
-
When these forces are concentrated at single points as in Figure 2, Michell [l] found the solution. It is a truss in which the members under compression are orthogonal t o those
32 1
322
A . Pliilpott & G. Strutig
i n tension, as i s always the case i n an optimal design. "he precise geometry depends on the points a t which the forces a c t , but t h e design i t s e l f i s straightforward: t h e compress i o n members a r e d i s t r i b u t e d r a d i a l l y and t h e tension is concentrated on a s i n g l e circumferential member connecting t h e two forces. Certainly t h a t is not an accurate description of t h e fibers i n s i d e t h e p a t e l l a . On t h e other hand it is not absurd; both s t e r e o l o g i c a l d a t a and f i n i t e element analyses [2-41 confirm t h i s general p a t t e r n . Photographs of t h e trabeculae show a framework too s t r i k i n g t o be an accident. Already i n t h e l a s t century c l i n i c a l observations led t o the formulation of Wolff I s law, t h a t t h e f i b e r s a r e i n the d i r e c t i o n s of p r i n c i p a l s t r e s s . A more exact statement of t h a t law, and of the mechanisms thatgovern t h e laying down and resorbing of bone, i s badly needed. It i s a challenge t h a t s e v e r a l orthopaedic research groups have accepted, and we thank W. C. Hayes and h i s colleagues a t Beth I s r a e l Hospital f o r p a t i e n t encouragement t o explore t h e r e l a t i o n s h i p t o optimal design. A f i r s t s t e p (our immediate goal) I s t o r e f i n e Figure 2 by admitting a more r e a l i s t i c s e t of loads. This paper discusses t h e numerical optimization and not t h e biology. Even with t h i s l i l n i t a t i o n i t is incomplete: we could not pursue a l l possible approaches. Some t h a t we d i d pursue were not successful o r r a t h e r , t h e i r c o s t increased so rapidly t h a t a r e a l i s t i c three-dimensional p i c t u r e (with t h e bone allowed t o remodel a f t e r f r a c t u r e ) was inaccessible. Other algorithms a r e more promlslng. Perhaps we can l i s t here some of the a l t e r n a t i v e s , with preliminary comments:
-
1. Fixed Geometrg. Fix t h e nodes of t h e t r u s s and allow a specified subset of edges as candidates i n t h e design optimization. This y i e l d s a standard problem i n l i n e a r programming, with a sparse c o e f f i c i e n t matrix described below. The c o s t i s e f f e c t i v e l y t h e sum C a 1s I of edge lengths multiplied &Ill by member f o r c e s . The.simplex method w i l l y i e l d an optimal (and s t a t i c a l l y determinate) t r u s s .
Numerical & Biological Shape Optimizatiotz
323
It was hoped that the simplicity of t h i s model would allow a comparatively large number of nodes and edges, g i v i n g the design enough degrees of freedom t o imitate curved f i b e r s . We now believe t h a t t h i s approach i s not e f f i c i e n t . The approximation of a nonlinear design problem by l i n e a r programming can be very slowly convergent. In our experience the l i n e a r programs themselves were a l s o highly degenerate, w i t h multiple pivot steps before the cost was reduced.
Variable Geometry. The positions of the nodes a r e optimized a s well as the weight of the t r u s s t h a t connects them. Normally these steps are separate; a l i n e a r program computes the t r u s s f o r fixed nodes, and a nonlinear calculat i o n moves the nodes. O f the codes we know about (and there must be many othera) we mention three: 2.
( I ) The CONMIN code developed by (3. Vanderplaats [ S ] : This was generously made available t o us and was the basis f o r the conference paper by Snyder e t a 1 [ 6 ] on trabecular orientation i n the p a t e l l a . The code admits a much wider range of problem i n 2 and 3 dimensions, w i t h constraints on displacements and buckling s t i f f n e s s and frequency response. (11) The Michell structure code developed by McConnel [7]: This uses a d i r e c t search method of Powell t o optimize the nodal positions (with fixed node-edge incidence matrix)
.
Experimental r e s u l t s on problems of bridge design studied by Hemp [see 81 look very s a t i s f a c t o r y and we hope t o t e s t the code on biological structures. (iii) A newly developed code:
We are conpletlng a s p e c i a l purpose code t o t e s t other algorithms f o r this c l a s s of optimization problems. This paper describes the use of the dual variables from the fixed-node program t o choose the directions of movement of the nodes. It seems c l e a r t h a t t h i s variation of geometry i s valuable and probably e s s e n t i a l .
3. Michell Displacement Field. It was the great discovery of Michell that under the presribed load, the Optimsl design
324
A . Philpott & G. Strang
has p r i n c i p a l s t r a i n s of constant magnitude. ( I n a region where no members are required i n one o r both directions, t h e corresponding s t r a i n s a r e uniformly below t h a t constant value.) Lagache [ 9 ] has found a b e a u t i f u l existence proof f o r Michell t r u s s e s based on a convergent sequence of approximations t o t h i s s t r a i n f i e l d . Although those p a r t i c u l a r approximations nay be i n e f f i c i e n t i n computations they a r e based on fixed nodes, as i n 1. above t h e r e seems t o be no reason why a program t o compute f i n i t e element dlsplacements f o r Michell's s t r a i n f i e l d should not succeed.
-
-
Michell s t r e s s f i e l d . I n the s t r e s s formulation of t h e problem, t h e c o n s t r a i n t d i v cf = 0 can be removed by deriving t h e sOress u from Airy's function t(x,y): 4.
same
noted i n [lo], t h i s y i e l d s an ordinary (but nonquadratlc) problem in t h e calculus of v a r i a t i o n s . Again a combination of f i n i t e elements and nonlinear optimization would give a d i r e c t approximation of the Michell d i r e c t i o n s in t h e optimal design. The d e n s i t i e s and o r i e n t a t i o n s then come d i r e c t l y from 6 r a t h e r than i n d i r e c t l y from a d i s c r e t e t r u s s . As
This paper describes t h e new v a r i a b l e geometry code mentioned above. It was prepared i n l a t e 1983 and t e s t e d on t h e
-
b i o l o g i c a l applications although it applies equally t o engineering design. The next s e c t i o n of t h e paper r e c a l l s t h e s t e p s of the simplex method, modified only t o account f o r the absolute values i n the c o s t function; members can be i n tension or compression. Section 3 applies the simplex method t o c o s t ( o r volume) minimization, f o r a given placement of the nodes. Section 4 develops an algorithm f o r moving t h e nodes toward an improved design. There a r e numerical d i f f i c u l t i e s associated w i t h degeneracy, handled by a proj e c t i o n technique. In t h e l a s t s e c t i o n we b r i e f l y discuss the implementation and the applications.
Niitvierical& Biological Shape Optimization
325
THE FIXED GEOMETRY DESIGN PROBLEM
number of authors,startimg w i t h Dorn, Greenberg, and Gomory [U],have considered the application of l i n e a r programming t o the minimum weight problem f o r a plane t r u s s subjected t o forces along i t s boundary. "he problem is t o s e l e c t a s e t of members connecting a given collection of pin j o i n t s which a r e fixed i n number and position. Each member i n the optimal t r u s s is i n tension or compression a t i t s corresponding l i m i t s t r e s s . A
The formulation as a l i n e a r program is as follows. We consider m pin j o i n t s I n the plane, and l e t XiDyi, i = 1, .,m be the coordinates (with respect t o some o r i g i n ) of the l t h j o i n t . We l e t [ i , k ] denote the member connecting j o i n t s 1 and k . To avoid unnecessary subscripts we associa t e with each [ i , k ] a unique index j E (1,2,. .,n), and as the force i n member j . Here sj is positive define s j i f j is in tension and negative i f j is i n compression. j in "he number of members n can be a t most - , practice we r e s t r i c t our choice of members t o a r e l a t i v e l y small set, usually those having length below some chosen bound.
--
The forces must be such t h a t each j o i n t is in s t a t i c "j equilibrium. These 2 m equality constraints may be formulated by resolving each member force in the two coordinate directions and summing the forces i n each direction a t each node. The right hand sides of these constrainta are the x and y components of the applied forces. For j o i n t I, these a r e denoted by fli and fgi respectively. For many j o i n t s , including those i n t e r n a l t o the structure, the applied forces are zero. If member j = [ i , k ] makes en angle 0 w i t h the positive j x-axis a t j o i n t i, we may write the constraints f o r node i as
A . Plrilpotr & G. Strarig
326
where Ii is the set of members connected to the ith joint. Since each member connects exactly two joints, these equilibrium constraints in matrix form become
where 8 is the vector of n member forces, fl and f2 are the x and y components of the boundary forces, and A1 and A2 are m X n matrices with zero entries in column j = [ i , k ] except in row8 i and k where
-
(A21ij = -sin
6
3 ’
(A21kj = s i n 0j
I;:[
-
In what fOllOW8 we shall refer to the matrix and the right hand side vector
as
f.
A n example of a pinjolnted trusla and its constraint matrix
is given in Figure 3.
If
aj
is the length of member j, minimizing the volume of a planar truss is equivalent to minimizing
n over those vectors 8 which satisfy the equilibrlun condition As = f. The derivation of this result requires that the l i m i t stresses in tension (a=) be the same for
327
Nu 17 i er icul & Biological Sliap e Opt irri izat io n
each member, and t h a t the l i m i t s t r e s s e s under compression (udn) be the same f o r each member. Since W is independent of the values of uand %in , we may assume without 108s of generality that u- udn = u. When each member j i s f u l l y stressed i t s cross-sectional area is proportional t o
-
9:
t!i
aJ
= u
, for
j
i n tension,
for
j
i n compression.
= -6,
1s I measures the thickness required in member 3 j j the problem of minimizing the weight of a fixed j o i n t planar t r u s s may thus be formulated as
In this way
Problem P:
n
C djls,l j=l
Minimize
subject t o As = f
.
(1)
As posed, t h i s objective function is nonlinear. However it is easy t o reformulate P as a l i n e a r program by defining
variables
whence W become
C 1 t
J j A(t
+
Z
-
1u
U)
j j =
is l i n e a r and the constraints
f, t
2
0,
U
2
0
.
We s h a l l write the problem in the formulation P, but make considerable use of the l i n e a r programming formulation and the large body of theory which t h i s allows us. W e can formulate the problem dual t o D: Maximize
vTf l
+
P
as follows:
wTf2 subject t o I(vTA,+wTA2),I j =
1 , e . m
an
5 dj, (21
A. Philpott & G. Strang
328
Multiplying both s i d e s by I S il duclng f l = A1s and f 2 = A2s easy half of d u a l i t y theory.
I j
on jj and i n t r 0 leads t o weak d u a l i t y the
-
SUmming
Theorem 1: If v and w a r e f e a s i b l e f o r t h e dual and s i s f e a s i b l e f o r the primal then n T T (3) v f l + w f2 5 c i j l S j l
J-1
The complementary slackness conditions, which must hold a t the minimizing 6 and maximizing v and w, a r e as follows : if
sj
>
0, then
if
8
<
0, then
j
( V TA
1
+ wTA p ) j
aj (5)
The dual problem i s o f t e n given the following physical l n t e r pretation, with E a s Young's modulus. A s e t of dual variables v i J w i J corresponds t o a v i r t u a l deformation of the t r u s s I n which each j o i n t 1 i s displaced by v i ( i ) i n the y d i r e c t i o n . The elongation e j of each member j under such a deformation i s given by
The dual c o n s t r a i n t s t h e r e f o r e require of each elongation e j that
I n other words t h e l i n e a r s t r a i n i n each member caused by t h e v i r t u a l deformation must be a t most t h e value which corresponds t o the member being a t i t s limit s t r e s s . The dual objective function vTf l
+
wTf 2 can be i n t e r p r e t e d
as E/a times t h e work done by e x t e r n a l f o r c e s t o produce t h i s v i r t u a l deformation. Since t h e optimal s o l u t i o n s w i l l s a t i s f y (4) and ( 5 ) J t h e s t r a i n in a member with sj o produced by t h e optimal displacements v and w must be
+
Numerical & Biological Shape Optimization
329
the strain which would occur with j at its limit stress (a if sj > 0 and - u if sj < 0). Members with s = 0 J and they are below do not appear in the optimal structure the limit stress when subjected to the optimal v and w.
-
AN ALUORITHM FOR THE FIXED GEOMETRY PROBUM
Since the fixed geometry design problem P may be formulated as a linear program, the simplex algorithm will solve it. This method constructs a sequence of basic feasible solutions with decreasing weight C l? I s I . The optimum does occur at J 3 a basic feasible solution (defined below) and the sequence terminates at an optimal solution after a finite number of steps. For a planar truss with n joints we have 2m equilibrium constraints. However the rows of A are not linearly independent, and for equilibrium the boundary loads must sum to zero in both coordinate directions. Furthermore the sum of moments about any point in the truss must be zero. The number of independent equilibrium constraints is therefore 2m-3, and three constraints must be removed before the simplex algorithm can be applied.
In practice these degrees of freedom are removed by fixing two or more boundary Joints. Then we remove from A1 and the rows which correspond to the anchored joints, and we remove from the vectors v and w the displacements which are no longer possible. We shall assume in what follows that the constraint matrix A has undergone this process, eqaal to the and take m to be the number of rows of A number of displacement unknowns. This removes at the 8ame t h e any restrictions on the boundary loads; they need not be self-equilibrating since any additional forces necessary for equilibrium will be taken up by the foundation. A2
-
The basic fersible solutions for problem P are found by choosing 2m linearly independent columns of A and solving (1) uniquely for the member forces sj corresponding
A. Philpott & G. Strang
330
t o these columns. A l l other member forces a r e zero. Thus each basic f e a s i b l e s o l u t i o n corresponds t o a e t a t i c a l l y deteranate truss. W e s h a l l c a l l any such truss a basic f e a s i b l e t r u s s . S t r i c t l y speaking, t h e basic f e a s i b l e s o l u t i o n t o t h e l i n e a r program is given by t h e nonnegative u : variables s a t i s f y i n g s = t
-
and With this understanding we apply t h e term "basic f e a s i b l e solution" t o the member forces 8 and also t o t h e variables
t
and u.
If t h e choeen columns of
a r e not l i n e a r l y independent, t h e corresponding truss i s no longer r i g i d under an a r b i t r a r y set of boundary loads, and it becomes a mechanism. On the other hand, i f more than 2m columns a r e chosen (with 2 m l i n e a r l y independent) then t h e corresponding t r u s s i s s t a t i c a l l y Indeterminate, and we cannot solve uniquely f o r the member f o r c e s . A t h i r d p o s s i b i l i t y is t h a t t h e 2m columns a r e l i n e a r l y independent and y e t some of t h e member forces in t h e basic f e a s i b l e truss a r e zero. I n t h i s case t h e basic f e a s i b l e s o l u t i o n i s degenerate. A
The algorithm f o r the f i x e d geometry optimal design follows t h e s t e p s of t h e slmplex method, modified to work w i t h t h e formulation P. If we assume that degenerate basic f e a s i b l e s o l u t i o n s do not occur, t h e s t e p s a r e as follows:
[i;]
1. Determine a s e t of 2 m members which form a basic f e a s i b l e t r u s s . Denote by B = t h e 2m by 2m
"basis matrix" formed by t h e corresponding columns of 2.
A.
Solve Bs = f .
3. For t h e s o l u t i o n S, construct a complementary s l a c k dual s o l u t i o n v#w. I n o t h e r words, solve t h e l i n e a r system (vTB1
+
wTB )
2 j
=
-+ 4j
(choosing t h e s i g n of
Sj)
33 1
Numerical & Biological Shape Optimization
4.
Choose the value of
'J 5. If
2
=
aJ
-
j (say
j = q)
I(vTAl + wTA )
2J
I
which minimizes
.
2 0 then stop3 the current t r u s s is optimal. Q ra < 0 then member q enters the basic f e a s i b l e t r u s s in t'ension o r compression depending on the sign of T T = ( V A1 + w If
r
9
6.
Determine the member t o leave the basic f e a s i b l e truss, over given by the Index 1 = p which minimizes l q i l indices 1 i n I. Here P
si (B'lA
iq
I = ( I ; (B-lA)iq
0
and sign
=
sign zs)
.
7.
Replace member p i n the basic f e a s i b l e truss (and the column p in the matrix B) w i t h member q . Return t o Step 2.
It is easy t o see t h a t t h i s algorithm follows exactly the same steps as the simplex method. Step 4 f i n d s a m i n i m u m reduced cost, and ti i n Step 6 i s the r a t i o t e s t t h a t determines which basic variable becomes non-basic a f t e r a pivot operation. I n t h i s step, i f no quotients of approp r i a t e sign e x i s t then the problem has an unbounded solut ion.
Since the constraint matrix A i s very sparse, the implementation of the algorithm uses an LU f a c t o r i z a t i o n of the basis matrix B and updates the triangular f a c t o r s L and U a t each i t e r a t i o n . We employ the method of Forrest and Tomlin [12]. To take f u r t h e r advantage of sparsity, the columns of A corresponding t o nonbasic variables a r e recomputed a t Step 4 i n each i t e r a t i o n . The steps described above are followed w i t h no e s s e n t i a l changes when degeneracy is encountered. If the s t a r t i n g
332
A. Philpott & G. Strarig
basic feasible t r u s s i s degenerate we a r b i t r a r i l y assign positive or negative signs t o the zero forces. Since by Step 5 every member entering the basis does so e i t h e r under tension or compression, f o r every subsequent basic feasible t r u s s we can assign positive and negative signs t o the (possibly zero) member forces. This convention gives a unique complementary slack dual solution i n Step 3. The choice of member p t o leave the basic feasible t r u s s is however not uniquely determined, since any 1 w i t h si 0 w i l l render l$il a minimum i n Step 6. In the code we choose the variable with lowest index, and although t h i s rule does not preclude cycling, no instances of cycling have been observed. (Dantzig has observed only one i n long experience with simplex calculations. ) Of course a more expensive pivot rule could ensure t h a t cycling never occurs.
-
VARIABLE GEOMETRY
For fixed j o i n t positions the above algorithm gives an optimal solution it finds the minimum weight of a plane t r u s s which supports given boundary loads. Being e s s e n t i a l l y the simplex method, the algorithm terminates a t the optimum a f t e r a f i n i t e number of I t e r a t i o n s . Unfortunately, f o r large n and m (many j o i n t s and many possible members) the time taken t o reach the optimum is prohibitive. If members are permitted between any two j o i n t s then the problem s i z e Increases rapidly with the number of j o i n t s . To avoid t h i s problem we may r e s t r i c t the allowed members t o those whose length l i e s below some fixed bound. However, even a small problem w i t h 100 j o i n t s and 800 members can take hours of CPU time
-
.
One reason f o r t h i s poor convergence behavior i s the high degree of degeneracy inherent i n the l i n e a r programs associated w i t h pin jointed structures. Often a large number of j o i n t s have no external forces, and play no p a r t i n the designs the members incident upon these j o i n t s carry zero force. Even If these j o i n t s and the members connected t o them are deleted,
333
Numerical & Biological Shape Optimization
it is often the case i n a basic f e a s i b l e truss that other zero-force members e x i s t . As observed i n [ l l ] , these cannot be deleted since the resulting truss would no longer be r i g i d . The presence of degenerate basic f e a s i b l e solutions causes poor performance in the simplex algorithm, since many pivot operations occur w i t h no change i n the objective f unc t ion. To approximate accurately the trabecular architecture of the human patella, a large number of i n t e r n a l j o i n t s would be
-
i f t h e i r positions a r e fixed. needed i n the design process This p r o l i f e r a t i o n of j o i n t s subject t o zero loads exacerbates the degeneracy problem, and the optimization becomes Idq)ract i c a b l e . A n a l t e r n a t i v e is t o admit a much smaller number of j o i n t s and then alter t h e i r positions i n such a way a8 t o improve the design. We proceed t o describe a method f o r computing the gradient of the cost W in terms of the j o i n t positions, thus giving a s e t of directions of steepest descent; s h i f t i n g the j o i n t s i n these directions produces a truss w i t h l e s s weight.
To f i n d the gradient, we investigate the e f f e c t on the objective function W of a change i n position of the j o i n t s from (Xisyi) t o ( x i + 8Xj.s y i + 8yi). A t optimal solution of problem P# the cost I s n w = jc-1 aj l a3 I .
The sum is taken over members or" the optimal basic f e a s i b l e t r u s s . L e t this truss correspond t o the basis matrix B =
LB,'] 1 rB,
where
B1
and
B2
are
m x 2 m matrices.
We
assume-for the present that this truss is not degenerate. Changing the positions of the j o i n t s w i l l change the length8 a j and a l s o a f f e c t the e n t r i e s of B implying that the member forces will a l s o change. Suppose that the forces become s j + 8 s j 8 and the changes in angles a l t e r B t o B + bB. Then equation (1) gives
-
A. Philpott & G. Strang
334
(B
+
+
bB)(s
6s) = f ,
whence t o f i r s t order
+
B6s
(6B)s = 0 .
The f i r s t order change i n W = C
a 3 1s3
+
C 0aj
a
8s
J J
given by Sjl
+
2 Jj61Sj(.
sj'o
The last term i n (7) is absent by the assumption of nondegeneracy. It then follows from the definition of the complementary slack dual variables (the displacements v w ) that
and thus from (6) we obtain
r
i
L
J
(7)
and
We now express &4? and 6B i n terms of the changes in j o i n t positions, Ox and by. If j = [ i , k ] , then the is (xk + (yk Y , ) ~ . Since length
$
2
-
-
=
COB
ej,
aa 4 --
COB 0
j '
& aa = - s i n e j , 4, Yi
2 = s i n e
it follows from the d e f i n i t i o n of 6a
-
T B16x
T + B26y.
B1
and
Be
that (9)
The matrix 8B is the f i r s t order change i n the matrix B and a r i s e s from the changes i n COB e and sin 8 f o r J 3 each member a f t e r the positions o f some j o i n t s are changed. It is easy t o ehow that if j = [ i , k ] ,
335
Numerical & Biological Shape Optimization
7 -
cos e
yi
s i n ej ,
a
cos eijsin ej
cos dyk e
j
Therefore the changes i n cosine and sine a r e given t o f i r s t order by
The change ( 6 B ) j i n column j of B I s zero except i n , k, i + m, and k + mj i n these four rows we have rows I
w i t h the same formulas and opposite signs f o r ( 6B)k+m, j . The term i n brackets is j u s t
If we l e t Sjj
=
2, j
8W
S
be the
then
6W
=I4T[ BTl .TI[(.]- 8 x
2m x 2m
(FJB),~
diagonal matrix w i t h
I n ( 8 ) can be written as follows:
[vT wT 1[82] s[B; OB1
We simplify t h i s expression by defining
-B3
E];
and
336
A. Philpott & G. Strang
leaving 6W = gT b x + hT b y .
(13)
For given loads and a given nondegenerate optimal truss, the vector [If] defined by (11) and (12) i s a direction of steepest descent f o r W when expressed as a function of the nodal positions. Therefore it i s natural t o extend the algorithm t o include a s h i f t of positions i n t h i s descent direction. This variable-geometry algorithm proceeds as follows : Step 1. For a given s e t of j o i n t coordinates the fixed geometry problem P
-x,y,-
solve
Step 2.
For the optimal basic feasible solution t o calculate g and h from (11) and (12)
Step 3.
For the given boundary loads, find the value of Xh) the stepsize 1 which minimizes W(z Xg,Y
Step 4.
If the change I n W i s l e s s than some tolerance then terminate; otherwise update 'ji and and return t o S t e p 1.
-
Pa
-
The c r u c i a l step is the l i n e search i n Step 31 It is here t h a t the Improvement In geometry I s made. There are a number of ways t o carry out t h i s search. One p o s s i b i l i t y I s t o r e t a i n the current basic feasible truss and recalculate W f o r this truss a t each point i n the l i n e search (determining t h e member forces from equilibrium Ba = f ) . Alternatively, we may return t o the simplex algorithm and find the optimal basic feasible truss a t each point of the search. We tend t o favor the former approach) re-solving the l i n e a r
331
Numerical & Biological Shape Optimization
programming problem a t each point requires a greater comput a t i o n a l e f f o r t , a s i t u a t i o n which is not improved by the presence of degeneracy. One might hope t h a t t h i s algorithm gives a reasonably robust method of improving the geometry. Unfortunately t h i s is not the case.
In practice, l i n e searches often terminate a t values of A f o r which the resulting basic f e a s i b l e t r u s s is degenerate. This seem reasonable since the cross-sectional area of the degenerate member has a l o c a l minimum exactly a t the value of X which renders it degenerate. This contradicts our assumpt i o n ofnondegeneratebasic f e a s i b l e t r u s s e s j even i f we begin the l i n e search a t a well-behaved truss, we a r e l i k e l y t o meet degeneracy. possible remedy is t o continue as i f our assumption were valid, and hope that the sum X 1 61s I over members with J J sj = 0 (which was neglected i n ( 7 ) ) is small. However it soon becomes obvious t h a t this approach is untenable. In practice, l i n e searches from current degenerate basic f e a s i b l e trusses often terminate a t a steplength A which is zero t o working accuracy, and the algorithm is j m e d a t a suboptimal solution. Clearly a t some point, t h i s neglected fourth term i n (7) has a considerable e f f e c t . A
The solution i s t o compute a more e f f e c t i v e search direction. Since the undesirable behavior ie caused by changes i n 1s I J f o r degenerate members, It makes sense t o seek a direction which minimizes these changes. We therefore search in a direction which, a t l e a s t locally, keeps degenerate members degenerate. Then one expects t h e i r contribution t o 8W t o be small. The adjusted search direction is computed as follows. we have 8s = -B-l(8B)s
,
From (6)
A. Philpott & G. Strang
338
whence, using (lo),
Let Id be the set of indices of columns of B corresponding to degenerate members, and let Bd be the matrix formed by the rows of Bol with these indices. Then if
It;] -
D = -Bd[:il]
S [BE -Bl], T any vector
_
I
in the nullspace of
D will render (6s) = 0 for j E Id. Thus to obtain a descent direction which keeps the changes in degenerate onto the members as smell as possible, we project nullspace of D :
[E]
[f] =
[g] - DT(DDT)-l
D[E]
.
Since DDT is often singular, we use a Gram-Schmidt procedure t o construct a matrix Q with orthonormal columns spanning the column space of DT, whence
Even with this projected search direction, we are not for large step sizes guaranteed a decrease in W(%XE,T-Ali) A. Indeed, if a step is so large that some member forces change sign, the current choice of direction may no longer give descent. For this reason, we compute from (14) the values of X for which a member force changes sign, and evaluate the objection function W at these values in succession. (The equilibrium equation is solved for 8.) When an increase in W is detected, this gives an interval in X over which a golden section line search minimizes W. C ONC LUS ION
The algorithm described above was implemented in FORTRAN on a VAX 780, and applied to a model of the human patella. The loads correspond to a 60’ flexion of the leg (as in
Numerical & Biological Shape Optimization
339
climbing s t a i r s ) . We worked with 40 j o i n t s , and allowed a l l edges w i t h length not exceeding one centimeter. An i n i t i a l basic feasible t r u s s i s shown i n Figure 4a. The 14 fixed nodes along the bottom provide reactions t o the 5 applied loads. Solid l i n e s denote members i n tension and dashed l i n e s indicate compression. The thickness of the members i s represented by the number of l i n e s i n the diagram, and bars denoted by single l i n e s i n Figure 4a are degenerate. They have zero thickness, as i n the two members which meet a t the unloaded node a t the top right hand side; they a r e present only t o avoid a mechanism. The weight of the i n i t i a l t r u s s i s W = 2.064 grams. Figure 4b shows the output from the fixed geometry design algorithm, s t a r t i n g from the truss of Figure 4a. Degenerate members ,are now omitted. The optimal t r u s s has a weight of W = 1.829 gram and was obtained a f t e r 25 pivot operations i n the simplex method. For a problem of t h i s size, each pivot step takes approximately .5 seconds, a value which increases substantially f o r l a r g e r problems. Figure 4c shows the output from the variable geometry algorithm. It i s the o p t i m l basic f e a s i b l e t r u s s a f t e r a succession of changes i n nodal positions. The weight i s reduced t o W = 1.708 gram a f t e r t e n l i n e searches, a f t e r each of which the fixed geometry problem was resolved and new directions of movement were computed (with projection t o control degeneracy). Subsequent l i n e searches f a i l e d t o give a significant impuovement i n the weight of the trusa and the program terminated. ACKNOWLEDGEMENT We thank the Army Research Office and the National Science Foundation for their support under contracts We are also grateful to the National Institute for Health for the opportunity, under
DAAG 29-83K-0025 and 81-02371.
contract AM30875, to work in the Orthopaedics Research Laboratory at Beth Israel Hospital.
3 40
A . Plzilpott & G. Strang
REFERENCES [l] Mlchell, A.G.M., The limit of economy of material in frame structures, Phil. Mag-8 (1904) 589-597. [2] Hayes, W.C., Snyder, B., Levine, B.M. and Ramaswamy, S., Stress-morphology relationships in trabecular bone of the patella,in: R. H. Gallagher et al., eds., Finite Elements in Biomechanics ( J o h n Wiley, New York, 1983). [3] Huberti, H.H., Hayes, W.C., Stone, J . L . and Shybut, G.T., Force ratios In the quadriceps tendon and ligamentum patellae, J. Biomechanics, to appear. [ 41
Stone, J . L . , Three-dimensional stress-morphology relationships In trabecular bone, M.Sc. Thesis, M. I . T . (1983).
-
, Numerical Optimization Techniques for Engineering Design With Applications (McGraw-Hill, New York, 1984).
[ 5 1 Vanderplaats, G.
[61 Snyder, B., Strang, Go, Hayes, W.C. and Norris, G., Application of structural geometry optimiration techniques in the microstructural remodeling of trabecular bone, ASME Annual Meeting, Boston, 1983.
[TI MCCOMel, ROE., Least-weight frameworks for load across Spar
ASCE J.
H e w W.5.r 19731 191
m.MeCh.
DIv. 100 (1974) 885-901.
O;ptlmU Structures (Clarendon Press, Oxford,
Lagache, J.-M., Treillis do volume minimal dans une region donnh!, J. Mdcanique 20 (1981) 415-448.
[lo] Strang, G. and Kohn, R.V., Hencky-Prandtl nets and constrained Michell trusses, Comp. Methods in Appl. Mech. E W . 36 (1983) 207-222. 1111 Dorn, W.S., Gomory, R.E. and Greenb,erg, H . J . , Automatic design of optimal structures, J. Mecanique 3 (1964) 25-52.
1121 Forrest, J . J . H . and Tomlin, J.A., Updated triangular factors of the basis to maintain sparsit in the product form simplex method, Math. Prog. 2 (19727 263-278.
Numerical & Biological Shape Optimization
FIGURE 1. HUMAN PATELLA: ROTATED CROSS SECTION, WITH APPLIED TENSIONS AND SUPPORT REGION.
F' FIGURE 2. MICHELL TRUSS, 3-POINT LOAD.
341
A . Philpott & G. Strang
342
3
FIGURE 3. A SIMPLE PINJOINTED TRUSS.
scale:
1 cm.
FIGURE 4A. INITIAL TRUSS.
Numerical & Biological Shape Optimization
F I G U R E 4 B . OPTIMAL TRUSS WITH F I X E D NODES.
FIGURE 4C.
OPTIMAL TRUSS AFTER 10 GEOMETRY CHANGES.
3 43