- Email: [email protected]
Random search method based on exact reanalysis for topology optimization of trusses with discrete cross-sectional areas
Computers and Structures 79 (2001) 673±679
www.elsevier.com/locate/compstruc
Random search method based on exact reanalysis for topology optimization of trusses with discrete cross-sectional areas Makoto Ohsaki * Department of Architecture and Architectural Systems, Kyoto University, Sakyo, Kyoto 606-8501, Japan Received 15 September 1999; accepted 18 April 2000
Abstract A simple formulation applicable to topology modi®cation is presented for exact reanalysis of static response of trusses based on general methods of inverting modi®ed matrices. A random search algorithm is also presented for topology optimization of trusses with discrete list of cross-sectional areas. It is shown that the exact reanalysis is very eective for the case where the cross-sectional area of only one member is modi®ed at each step for generating a neighboring solution. The diculty due to discontinuity in the stress constraints is successfully avoided by considering the cross-sectional areas as discrete variables. Ó 2001 Elsevier Science Ltd. All rights reserved. Keywords: Exact reanalysis; Truss; Topology optimization; Random search; Discrete variable
1. Introduction In the process of designing structures, response analysis for a given set of design loads is carried out many times before reaching an admissible design that is optimal in the designer's point of view. It is important especially for the optimum design process that the computational costs for reanalysis should be reduced. Therefore, extensive research has been carried out for seeking ecient methods of reanalysis. First order sensitivity analysis is the most convenient way of estimating the response after modi®cation of design variables. It is well known for skeletal structures that the use of reciprocal variables leads to more accurate estimate of response than the case where the crosssectional areas or stinesses are directly used as design variables [1±3]. Higher order sensitivity coecients may be incorporated for more accurate estimates [4,5]. Substantial computational eort is needed, however, for calculation of Hessian that is needed for computing the second order sensitivity coecients. Polynomial ®tting *
Tel.: +81-75-753-5733; fax: +81-75-753-5733. E-mail address: [email protected] (M. Ohsaki).
algorithms may also be used for obtaining approximate response [6]. Many gradient-based algorithms have been developed for optimum design of ®nite dimensional structures assuming that the cross-sectional areas of a skeletal structure can have arbitrary non-negative values. In the design practice, however, the cross-sectional areas are to be selected from a catalogue of available values. Although the optimum designs in the context of continuous variables can be found within a moderate computational cost, a search method that requires large computational cost for reanalysis should be carried out for ®nding the optimal solution with discrete variables. Therefore, ecient reanalysis method for the given changes in the design variables is strongly desired in the ®eld of structural optimization with discrete variables. Melosh and Luik [7] presented a simple approach for reanalysis. Argyris and Roy [8] presented more general approach allowing changes in the number of freedom of displacements and support conditions. Kirsch and Rubinstein [9] discussed convergence properties of the iterative process. Recently, Kirsch [10] presented an ecient reanalysis method combining the reduced basis method and the expansion method. His method has been shown to lead to exact response of trusses [11]. The
0045-7949/01/$ - see front matter Ó 2001 Elsevier Science Ltd. All rights reserved. PII: S 0 0 4 5 - 7 9 4 9 ( 0 0 ) 0 0 1 6 8 - 1
674
M. Ohsaki / Computers and Structures 79 (2001) 673±679
formulation is general and in a form applicable to any ®nite dimensional structures, and is also applicable to shape modi®cation and topological changes [12]. Calculation of inverse of a modi®ed matrix has been discussed in many ®elds of engineering and mathematics. Developments in general form of inverting matrices are summarized in the review article by Henderson and Searle [13]. It is critical in application to static analysis problem of structures that the inverse of the stiness matrix of the initial structure is not usually known; i.e. the matrix is only decomposed to triangular matrices. Kavlie et al. [14] developed a stiness based method for computing the displacements of the modi®ed design based on the Sherman±Morrison formula [15]. They showed that the displacements are found without computation of the inverse matrix. In this paper, the exact reanalysis method [7,11,14] is reformulated for trusses based on the Sherman±Morrison formula for single-rank modi®cation. A simple formulation without inverse of the stiness matrix is presented for modi®cation of several cross-sectional areas including addition of members. It is shown that the formulation similar to that by Melosh and Luik [7] can be derived from the formula by Woodbury [16] for multi-rank modi®cation. A random search algorithm is also presented for topology optimization of trusses with discrete list of cross-sectional areas. In view of computational cost for reanalysis, the cross-sectional area of only one member is changed at each step for generating a neighboring solution. The diculty due to discontinuity in the stress constraints is successfully avoided by considering the cross-sectional areas as discrete variables.
denoted by subscript i. Since the components of K are proportional to Ai , K is written as K
m X Ai Ki ;
2
i1
where Ki is the stiness matrix of the ith member. Note that member properties are de®ned by the matrices and vector of the size n for convenience of the matrix operation. Computation is carried out, however, by elementsize matrices and vectors. Consider a case where the cross-sectional area of the ith member is increased by DAi . Eq. (1) for the modi®ed truss is formulated as K DAi Ki u1 P:
3
The objective here is to ®nd the displacement vector u1 of the modi®ed truss for the given load vector P. In the following, the superscripts 0 and 1 denote the values of initial and modi®ed structures, respectively. There have been variety of studies for computing the inverse of the modi®ed matrix from the known inverse of the initial matrix. For our purpose, however, those formulas cannot be directly used, because the inverse of K is not usually known; i.e. K is only decomposed for computing u0 . Therefore, it is critical to derive a more explicit formula for computing u1 of the modi®ed truss. Let N 2 Rm be the vector of axial forces. The relation between Ni0 and u0 is written by using a vector ci as Ni0 Ai cTi u0 ;
4
where the superscript T denotes the transpose of a vector. The relation between Ni0 and the vector of equivalent nodal forces F0i 2 Rn of the ith member is written as F0i Ni0 bi :
2. Exact reanalysis of trusses There are three directions for approximate or exact reanalysis of responses of structures subjected to static loads; i.e. expansion method, reduced basis method, and inversion of modi®ed matrices. Kirsch [10] presented an ecient and general approach combining expansion and reduced basis methods, and showed that exact response can be found for trusses. In this section, a more explicit formulation of reanalysis of trusses is presented based on general formulas for inverting modi®ed matrices [13]. Consider a truss subjected to static nodal loads P 2 Rn , where n is the number of freedom of displacements. Let K 2 Rnn denote the stiness matrix. The nodal displacement vector u0 2 Rn is calculated from Ku0 P:
1 m
The matrix K is a function of the vector A 2 R of the cross-sectional areas, where m is the number of members. In the following, the ith element of a vector is
5
Let E and Li denote the elastic modulus and the length of the ith member. Then the relation ci
E bi Li
6
holds and Ki is de®ned by using ci and bi as Ki bi cTi :
7
Note that bi is the nodal load vector corresponding to unit axial force of the ith member. The inverse of the matrix K DAi Ki obtained from [13] K DAi Ki ÿ1 K DAi bi cTi ÿ1 Kÿ1 ÿ
DAi Kÿ1 bi cTi Kÿ1 ; 1 ÿ DAi cTi Kÿ1 bi
8
where the inverse Kÿ1 , which is not usually known, is needed for calculating the inverse of the modi®ed
M. Ohsaki / Computers and Structures 79 (2001) 673±679
matrix. By post-multiplying P to Eq. (8), the following relation is derived: DAi Kÿ1 bi cTi u0 u1 u0 ÿ ; 1 ÿ DAi cTi Kÿ1 b
9
DAi ui cTi u0 DAi Ee0i u0 ÿ ui ; T i 1 ÿ DAi ci u 1 ÿ DAi Eei i
10
where Eq. (4) has been used, and ei is the strain of the ith member. This way, u1 can be obtained without using Kÿ1 . The strain of the member with null cross-sectional area may be obtained from the displacements of the nodes connected to the member. Therefore, Eq. (10) can be used for addition of a member. Next, we consider a case where the cross-sectional areas of q members are modi®ed simultaneously. It is assumed without loss of generality that Ai of members 1; 2; . . . ; q are changed. The stiness matrix of the modi®ed truss is formulated as K
q q X X DAi E T DAi Ki K bi bi K BDBT L i i1 i1
11
where the ith column of B 2 Rnq is bi , and D 2 Rqq is a diagonal matrix de®ned as D diagfDAi E=Li g:
Dÿ1 Y x y;
17
where the i; j component of Y is equal to ei j Lj . Then u1 is computed from u1 u0 U x:
where Eqs. (1) and (3) have been used. Let superscript i denote the values of the initial truss subjected to the load vector bi . Then Eq. (9) is rewritten as u1 u0 ÿ
675
12
18
It may be easily seen that Eq. (10) for modi®cation of single variable is included in Eqs. (17) and (18) for multiple variable modi®cation. Note that the formulation derived for the multi-variable modi®cation is similar to that by Melosh and Luik [7]. As an illustrative example, reanalysis has been carried out for the plane truss as shown in Fig. 1. The length of the members in x- or y-direction are 2.0 m; and P 100:0 kN. The displacements in the y-direction of the nodes a±c for the case of Ai 10:0 cm2 for i 6 1 and A1 0 are as shown in the ®rst row of Table 1, where members 1±3 are as de®ned in Fig. 1. Reanalysis has been carried out for the following three cases: Case 1: Use Eq. (10) for estimating the displacements for A2 15:0 cm2 . Case 2: Increase A1 from 0 to 5.0 cm2 and use Eq. (10) for reanalysis. Case 3: Increase A2 and A3 to 15.0 cm2 simultaneously and use Eqs. (17) and (18). The results are as listed in Table 1. It has been con®rmed that the exact responses have to be found for all the cases including Case 2 corresponding to the addition of a member.
The inverse of the modi®ed stiness matrix is written as [13] K BDBT ÿ1 Kÿ1 Kÿ1 B Dÿ1 BT Kÿ1 Bÿ1 BT Kÿ1 : 13 By post-multiplying P to Eq. (13), the following relation is derived: u1 u0 Kÿ1 B Dÿ1 BT Kÿ1 Bÿ1 DBT u0 u0 U Dÿ1 BT Kÿ1 Bÿ1 y;
14 Fig. 1. A 2 2 plane grid.
where the ith component of the vector y 2 Rq is e0i Li , and the ith column of the matrix U 2 Rq is ui . A temporary variable x 2 Rq is introduced as x Dÿ1 BT Kÿ1 Bÿ1 y:
15
The vector x is found by solving Dÿ1 BT Kÿ1 Bx y which is rewritten as
16
Table 1 Result of reanalysis for displacements (cm) Initial Case 1 Case 2 Case 3
ua
ub
uc
ÿ0.223732 ÿ0.211225 ÿ0.220397 ÿ0.210492
ÿ0.243724 ÿ0.230694 ÿ0.238722 ÿ0.229644
ÿ0.305510 ÿ0.290569 ÿ0.302175 ÿ0.289815
676
M. Ohsaki / Computers and Structures 79 (2001) 673±679
3. Optimization of trusses from discrete list of crosssectional areas Consider inequality constraints for the response quantity g^j A such as nodal displacement and member stress as g^j A 6 gj
j 1; 2; . . . ; f ;
19
where gj is the upper bound for g^j A, and f is the number of constraints. Lower bound may be given similarly. Then the constraint is normalized assuming gj 6 0 as gj A
g^j A ÿ 1 6 0: gj
C A
subject to gj A 6 0
4.
20
The objective function C A is the total structural volume. The value of Ai has to be selected from the list fA0k g. Then consider the following optimum design problem: minimize
3.
21 j 1; 2; . . . ; f
Ai 2 fA0k g k 1; 2; . . . ; r;
22
where r is the number of members in the list. There have been many methods presented for optimization of structures with discrete variables [17]. The purpose of this section is to show eectiveness of the exact reanalysis method applied to a topology optimization of trusses. A random search algorithm is used in the following, where a new technique is presented for improving convergence to the optimal solution and for avoiding the diculty in de®ning appropriate value for the penalty terms. The most simple approach for using a random search for constrained optimization problem may be to add a penalty term of the constraints to the objective function to formulate a performance measure. In this case, however, it is very dicult to assign an appropriate value for the penalty coecient so that the magnitude of objective function is equivalent to that of the penalty term. In the proposed method, the penalty for violating the constraints is implicitly given by the threshold level for de®ning the increase or decrease of the design variables. The optimization algorithm is as follows: 1. Assign initial value of A from the given list fA0k g, and set parameters G; s and d. 2. Evaluate the objective function C A and the constraints gj A, and ®nd the maximum value G max gj . At this step, the matrix K is reconstructed j and decomposed. Note that G is negative if all the constraints are satis®ed. Assign the parameter
5. 6.
k 0:5 sG, where s moves the threshold level for increase or decrease of the variable. The parameter k is replaced with 0 or 1 if k < 0 or k > 1; respectively. Generate a uniform random number between 0 and 1 to select the member k for which cross-sectional area is to be modi®ed. Suppose there are r allowable values in the increasing order in the list fA0j g, and Ak A0i at the current step. Generate another random number p and calculate the smallest integer q that is not less than jk ÿ pj=d. Modify Ak to A0j if p < k, where j minfr ÿ i; qg, otherwise decrease Ak to A0j , where j maxfi ÿ q; 1g. Compute constraint functions for the modi®ed design, where the exact reanalysis method is to be used here. Accept the design and go to 2 if G < G and the number of the accepted designs is less than the speci®ed limit. Go to 3 if the number of trial at this stage is less than the speci®ed limit. Find the optimal design from the list of the accepted designs.
As an example, the proposed algorithm has been applied to topology optimization of plane square grids under constraints on member stress. The case of 2 2 and 5 5 are as shown in Figs. 1 and 2. The length of the members in x- or y-direction are 2.0 m, and P 200:0 kN. Stress constraints are given for all the members, where the upper bound for the absolute value is 40.0 MPa. The parameters are G 0:3, s 10:0 and d 0:1. Note that local buckling constraints are not considered here. The main diculty of topology optimization under stress constraints arises from the fact that the constraint
Fig. 2. A 5 5 plane grid.
M. Ohsaki / Computers and Structures 79 (2001) 673±679
677
Table 2 Result of optimization with 100 dierent random value parameters 22 33 44 55
u
sa
ta
Va
V min
V max
0.0025 0.005 0.007 0.01
469 288 564 1026
8278 4303 7933 14 297
0.00729698 0.130939 0.184428 0.256714
0.00692918 0.103967 0.138661 0.172274
0.0892801 0.179940 0.254913 0.364895
need not be satis®ed by the vanishing member [18]. Therefore, the optimal solution is sometimes an isolated point of the feasible region in the design variable space, and it is not useful to use the optimal solution with continuous variable [19] as the initial solution for the problem with discrete variable. The list of allowable cross-sectional areas (cm2 ) is given as fA0j g f0:01; 20; 40; 60; 80g. Note that the member with Ai 0:01 cm2 is to be removed and the stress constraint need not be satis®ed by such a member. An arti®cial upper bound u has been given for the nodal displacements to prevent the design such that the crosssectional areas of all the members are equal to 0.01 cm2 is regarded as the optimum design because obviously it has the minimum volume and no constraint should be satis®ed [20]. Optimization has been carried out for 2 2, 3 3, 4 4 and 5 5 grids, where 100 dierent initial values of the random value are generated for each case. The limit for number of modi®cation of the cross-sectional area is 5000, and the limit for the trial for generating an acceptable neighboring solution at each step is 2000. Let sa and ta denote the average numbers of steps and trial analysis before reaching the minimum value of the objective value within 5000 steps. Average, minimum and maximum values of the minimum objective function
within 5000 steps among 100 dierent random value parameters are denoted by V a , V min and V max respectively. Table 2 shows the values of u, sa , ta , V a , V min and max V for each case. Note that optimal solutions have not been reached for all the cases of dierent random value parameters. It has been con®rmed that an appropriate value of u leads to rapid convergence to nearly optimal solutions. The values of u is a little greater than the maximum displacement of the truss with optimal topology for each case. The optimum design of 5 5 grid is as shown in Fig. 3. The cross-sectional areas are 60.0 cm2 for members 6± 10 and 80.0 cm2 for the remaining existing members. The value of V min for each case in Table 2 corresponds to the objective value of the apparent optimal solution similar with the truss in Fig. 3. It may be observed from Table 2 that the average number sa of analysis needed to decompose the stiness matrix is less than 10% of the average number ta of trial analysis that can be done by the reanalysis formula. The ratio of the CPU time for obtaining the optimal topology of the 5 5 grid by using the reanalysis method to that without reanalysis is 0.2250. Therefore, computational cost may be dramatically reduced by using the proposed method of reanalysis. Next, optimal topology has been found for the rectangular 10 5 grid as shown in Fig. 4, where P 100:0 kN. Other parameters are same as those for the square grids. Optimization has been carried out for 100 dierent values of random value parameter. The results for
Fig. 3. The optimal solution for the 5 5 plane grid.
Fig. 4. A 10 5 plane grid.
678
M. Ohsaki / Computers and Structures 79 (2001) 673±679
Fig. 5. The optimal solution for the 10 5 plane grid.
the arti®cial upper-bound displacement u 0:02 cm are sa 3732, ta 40 690, V a 0:443299, V min 0:253592 and V max 0:651226. The solution corresponding to V min that coincides with the obvious solution is as shown in Fig. 5. Although the ratio of V a to V min is not suciently small for this case, it is important to obtain the optimal topology with simple algorithm that does not involve diculty in assignment of the parameters.
4. Conclusion The exact reanalysis method has been reformulated based on the general formula of inverting a modi®ed matrix. A simple expression that leads to exact estimate of the responses has been derived for modi®cation of cross-sectional areas of a truss. In this method, only the responses of the initial truss is computed for the equivalent load corresponding to the unit axial force of the member of which the cross-sectional areas should be modi®ed, and complicated matrix operation that are necessary for estimating the inverse of the stiness matrix of the modi®ed structure are not needed. It has been shown that the formulation similar to that by Melosh and Luik [7] can be derived from the formula by Woodbury [16] for multi-rank modi®cation. Note that the formulation based on the strains should be used for estimating the responses after addition of a member, because the stress of the member with null cross-sectional area cannot be calculated while the strain of such a member can be obtained from the displacements of the nodes connected to the member. A random search algorithm has been presented for trusses with discrete list of cross-sectional areas. The penalty for the violation of the constraints is incorporated as the threshold level for deciding increase or decrease of the cross-sectional areas for a given random number. Therefore, there is no diculty in de®ning the penalty coecients considering the balance of the ob-
jective function and the penalty for the constraint violation. It has been shown that the exact reanalysis is very useful for the case where the cross-sectional area of only one member is modi®ed at each step for generating a neighboring solution. Eectiveness of the proposed method has been demonstrated through examples of square and rectangular square grids. The diculty due to discontinuity in the stress constraints has been successfully avoided by considering the cross-sectional areas as discrete variables. It has been shown that an appropriate value of arti®cial upper bound displacement leads to a rapid convergence to the optimal topology under stress constraints.
References [1] Noor AK, Lowder HE. Structural reanalysis via a mixed method. Comput Struct 1975;5:9±12. [2] Storaasli OO, Sbieszczanski J. On the accuracy of the taylor approximation for structural resizing. AIAA J 1983;21:1571±80. [3] Haftka RT, G urdal Z, Kamat MP. Elements of structural optimization. London: Kluwer Academic Publishers, Dordrecht, 1990. [4] Fleury C. Ecient approximation concepts using second order information. Int J Num Meth Engng 1989;28:2041± 58. [5] Fuchs MB. Nth-order stiness sensitivities in structural analysis. Struct Optimizat 1993;5:207±12. [6] Haftka RT, Nachlas JA, Watson LT, Rizzo T, Desai R. Two-point constraint approximation in structural optimization. Comp Meth Appl Mech Engng 1987;60:289± 301. [7] Melosh RJ, Luik R. Multiple con®guration analysis of structures. Proc ASCE 1968;94(ST11):2581±96. [8] Argyris JH, Roy R. General treatment of structural modi®cation. Proc ASCE 1972;98(ST2):465±92. [9] Kirsch U, Rubinstein MF. Structural reanalysis by iteration. Comput Struct 1972;2:497±510. [10] Kirsch U. Ecient sensitivity analysis for structural optimization. Comp Meth Appl Mech Engng 1994;117: 143±56. [11] Kirsch U, Liu S. Exact structural reanalysis by a ®rst-order reduced basis approach. Struct Optimizat 1995;10:153±8. [12] Kirsch U, Liu S. Structural reanalysis for general layout modi®cation. Proc 6th AIAA/NASA/ISSMO Symposium on Multidisciplinary Analysis and Optimization. AIAA 1996. p. 1142±9. [13] Henderson HV, Searle R. On deriving the inverse of a sum of matrices. SIAM Rev 1981;23:53±60. [14] Kavlie D, Graham H, Powell GH. Ecient reanalysis of modi®ed structures. Proc ASCE 1971;97(ST1):377±92. [15] Sherman J, Morrison WJ. Adjustment of an inverse matrix corresponding to changes in one element of a given matrix. Ann Math Statist 1950;21:124±7. [16] Woodbury MA. Inverting modi®ed matrices. Memorandum Report 42, Princeton: NJ, Statistical Research Group, 1950.
M. Ohsaki / Computers and Structures 79 (2001) 673±679 [17] Arora JS, Elwakeil OA, Chahande AI, Hsier CC. Global optimization methods for engineering applications. Struct Optimizat 1995;9:137±59. [18] Sved G, Ginos Z. Structural optimization under multiple loading. Int J Mech Sci 1968;10:803±5.
679
[19] Cheng G, Guo X. e-relaxed approach in structural topology optimization. Struct Optimizat 1997;13:258± 66. [20] Ohsaki M. Genetic algorithm for topology optimization of trusses. Comput Struct 1995;57(2):219±25.
www.elsevier.com/locate/compstruc
Random search method based on exact reanalysis for topology optimization of trusses with discrete cross-sectional areas Makoto Ohsaki * Department of Architecture and Architectural Systems, Kyoto University, Sakyo, Kyoto 606-8501, Japan Received 15 September 1999; accepted 18 April 2000
Abstract A simple formulation applicable to topology modi®cation is presented for exact reanalysis of static response of trusses based on general methods of inverting modi®ed matrices. A random search algorithm is also presented for topology optimization of trusses with discrete list of cross-sectional areas. It is shown that the exact reanalysis is very eective for the case where the cross-sectional area of only one member is modi®ed at each step for generating a neighboring solution. The diculty due to discontinuity in the stress constraints is successfully avoided by considering the cross-sectional areas as discrete variables. Ó 2001 Elsevier Science Ltd. All rights reserved. Keywords: Exact reanalysis; Truss; Topology optimization; Random search; Discrete variable
1. Introduction In the process of designing structures, response analysis for a given set of design loads is carried out many times before reaching an admissible design that is optimal in the designer's point of view. It is important especially for the optimum design process that the computational costs for reanalysis should be reduced. Therefore, extensive research has been carried out for seeking ecient methods of reanalysis. First order sensitivity analysis is the most convenient way of estimating the response after modi®cation of design variables. It is well known for skeletal structures that the use of reciprocal variables leads to more accurate estimate of response than the case where the crosssectional areas or stinesses are directly used as design variables [1±3]. Higher order sensitivity coecients may be incorporated for more accurate estimates [4,5]. Substantial computational eort is needed, however, for calculation of Hessian that is needed for computing the second order sensitivity coecients. Polynomial ®tting *
Tel.: +81-75-753-5733; fax: +81-75-753-5733. E-mail address: [email protected] (M. Ohsaki).
algorithms may also be used for obtaining approximate response [6]. Many gradient-based algorithms have been developed for optimum design of ®nite dimensional structures assuming that the cross-sectional areas of a skeletal structure can have arbitrary non-negative values. In the design practice, however, the cross-sectional areas are to be selected from a catalogue of available values. Although the optimum designs in the context of continuous variables can be found within a moderate computational cost, a search method that requires large computational cost for reanalysis should be carried out for ®nding the optimal solution with discrete variables. Therefore, ecient reanalysis method for the given changes in the design variables is strongly desired in the ®eld of structural optimization with discrete variables. Melosh and Luik [7] presented a simple approach for reanalysis. Argyris and Roy [8] presented more general approach allowing changes in the number of freedom of displacements and support conditions. Kirsch and Rubinstein [9] discussed convergence properties of the iterative process. Recently, Kirsch [10] presented an ecient reanalysis method combining the reduced basis method and the expansion method. His method has been shown to lead to exact response of trusses [11]. The
0045-7949/01/$ - see front matter Ó 2001 Elsevier Science Ltd. All rights reserved. PII: S 0 0 4 5 - 7 9 4 9 ( 0 0 ) 0 0 1 6 8 - 1
674
M. Ohsaki / Computers and Structures 79 (2001) 673±679
formulation is general and in a form applicable to any ®nite dimensional structures, and is also applicable to shape modi®cation and topological changes [12]. Calculation of inverse of a modi®ed matrix has been discussed in many ®elds of engineering and mathematics. Developments in general form of inverting matrices are summarized in the review article by Henderson and Searle [13]. It is critical in application to static analysis problem of structures that the inverse of the stiness matrix of the initial structure is not usually known; i.e. the matrix is only decomposed to triangular matrices. Kavlie et al. [14] developed a stiness based method for computing the displacements of the modi®ed design based on the Sherman±Morrison formula [15]. They showed that the displacements are found without computation of the inverse matrix. In this paper, the exact reanalysis method [7,11,14] is reformulated for trusses based on the Sherman±Morrison formula for single-rank modi®cation. A simple formulation without inverse of the stiness matrix is presented for modi®cation of several cross-sectional areas including addition of members. It is shown that the formulation similar to that by Melosh and Luik [7] can be derived from the formula by Woodbury [16] for multi-rank modi®cation. A random search algorithm is also presented for topology optimization of trusses with discrete list of cross-sectional areas. In view of computational cost for reanalysis, the cross-sectional area of only one member is changed at each step for generating a neighboring solution. The diculty due to discontinuity in the stress constraints is successfully avoided by considering the cross-sectional areas as discrete variables.
denoted by subscript i. Since the components of K are proportional to Ai , K is written as K
m X Ai Ki ;
2
i1
where Ki is the stiness matrix of the ith member. Note that member properties are de®ned by the matrices and vector of the size n for convenience of the matrix operation. Computation is carried out, however, by elementsize matrices and vectors. Consider a case where the cross-sectional area of the ith member is increased by DAi . Eq. (1) for the modi®ed truss is formulated as K DAi Ki u1 P:
3
The objective here is to ®nd the displacement vector u1 of the modi®ed truss for the given load vector P. In the following, the superscripts 0 and 1 denote the values of initial and modi®ed structures, respectively. There have been variety of studies for computing the inverse of the modi®ed matrix from the known inverse of the initial matrix. For our purpose, however, those formulas cannot be directly used, because the inverse of K is not usually known; i.e. K is only decomposed for computing u0 . Therefore, it is critical to derive a more explicit formula for computing u1 of the modi®ed truss. Let N 2 Rm be the vector of axial forces. The relation between Ni0 and u0 is written by using a vector ci as Ni0 Ai cTi u0 ;
4
where the superscript T denotes the transpose of a vector. The relation between Ni0 and the vector of equivalent nodal forces F0i 2 Rn of the ith member is written as F0i Ni0 bi :
2. Exact reanalysis of trusses There are three directions for approximate or exact reanalysis of responses of structures subjected to static loads; i.e. expansion method, reduced basis method, and inversion of modi®ed matrices. Kirsch [10] presented an ecient and general approach combining expansion and reduced basis methods, and showed that exact response can be found for trusses. In this section, a more explicit formulation of reanalysis of trusses is presented based on general formulas for inverting modi®ed matrices [13]. Consider a truss subjected to static nodal loads P 2 Rn , where n is the number of freedom of displacements. Let K 2 Rnn denote the stiness matrix. The nodal displacement vector u0 2 Rn is calculated from Ku0 P:
1 m
The matrix K is a function of the vector A 2 R of the cross-sectional areas, where m is the number of members. In the following, the ith element of a vector is
5
Let E and Li denote the elastic modulus and the length of the ith member. Then the relation ci
E bi Li
6
holds and Ki is de®ned by using ci and bi as Ki bi cTi :
7
Note that bi is the nodal load vector corresponding to unit axial force of the ith member. The inverse of the matrix K DAi Ki obtained from [13] K DAi Ki ÿ1 K DAi bi cTi ÿ1 Kÿ1 ÿ
DAi Kÿ1 bi cTi Kÿ1 ; 1 ÿ DAi cTi Kÿ1 bi
8
where the inverse Kÿ1 , which is not usually known, is needed for calculating the inverse of the modi®ed
M. Ohsaki / Computers and Structures 79 (2001) 673±679
matrix. By post-multiplying P to Eq. (8), the following relation is derived: DAi Kÿ1 bi cTi u0 u1 u0 ÿ ; 1 ÿ DAi cTi Kÿ1 b
9
DAi ui cTi u0 DAi Ee0i u0 ÿ ui ; T i 1 ÿ DAi ci u 1 ÿ DAi Eei i
10
where Eq. (4) has been used, and ei is the strain of the ith member. This way, u1 can be obtained without using Kÿ1 . The strain of the member with null cross-sectional area may be obtained from the displacements of the nodes connected to the member. Therefore, Eq. (10) can be used for addition of a member. Next, we consider a case where the cross-sectional areas of q members are modi®ed simultaneously. It is assumed without loss of generality that Ai of members 1; 2; . . . ; q are changed. The stiness matrix of the modi®ed truss is formulated as K
q q X X DAi E T DAi Ki K bi bi K BDBT L i i1 i1
11
where the ith column of B 2 Rnq is bi , and D 2 Rqq is a diagonal matrix de®ned as D diagfDAi E=Li g:
Dÿ1 Y x y;
17
where the i; j component of Y is equal to ei j Lj . Then u1 is computed from u1 u0 U x:
where Eqs. (1) and (3) have been used. Let superscript i denote the values of the initial truss subjected to the load vector bi . Then Eq. (9) is rewritten as u1 u0 ÿ
675
12
18
It may be easily seen that Eq. (10) for modi®cation of single variable is included in Eqs. (17) and (18) for multiple variable modi®cation. Note that the formulation derived for the multi-variable modi®cation is similar to that by Melosh and Luik [7]. As an illustrative example, reanalysis has been carried out for the plane truss as shown in Fig. 1. The length of the members in x- or y-direction are 2.0 m; and P 100:0 kN. The displacements in the y-direction of the nodes a±c for the case of Ai 10:0 cm2 for i 6 1 and A1 0 are as shown in the ®rst row of Table 1, where members 1±3 are as de®ned in Fig. 1. Reanalysis has been carried out for the following three cases: Case 1: Use Eq. (10) for estimating the displacements for A2 15:0 cm2 . Case 2: Increase A1 from 0 to 5.0 cm2 and use Eq. (10) for reanalysis. Case 3: Increase A2 and A3 to 15.0 cm2 simultaneously and use Eqs. (17) and (18). The results are as listed in Table 1. It has been con®rmed that the exact responses have to be found for all the cases including Case 2 corresponding to the addition of a member.
The inverse of the modi®ed stiness matrix is written as [13] K BDBT ÿ1 Kÿ1 Kÿ1 B Dÿ1 BT Kÿ1 Bÿ1 BT Kÿ1 : 13 By post-multiplying P to Eq. (13), the following relation is derived: u1 u0 Kÿ1 B Dÿ1 BT Kÿ1 Bÿ1 DBT u0 u0 U Dÿ1 BT Kÿ1 Bÿ1 y;
14 Fig. 1. A 2 2 plane grid.
where the ith component of the vector y 2 Rq is e0i Li , and the ith column of the matrix U 2 Rq is ui . A temporary variable x 2 Rq is introduced as x Dÿ1 BT Kÿ1 Bÿ1 y:
15
The vector x is found by solving Dÿ1 BT Kÿ1 Bx y which is rewritten as
16
Table 1 Result of reanalysis for displacements (cm) Initial Case 1 Case 2 Case 3
ua
ub
uc
ÿ0.223732 ÿ0.211225 ÿ0.220397 ÿ0.210492
ÿ0.243724 ÿ0.230694 ÿ0.238722 ÿ0.229644
ÿ0.305510 ÿ0.290569 ÿ0.302175 ÿ0.289815
676
M. Ohsaki / Computers and Structures 79 (2001) 673±679
3. Optimization of trusses from discrete list of crosssectional areas Consider inequality constraints for the response quantity g^j A such as nodal displacement and member stress as g^j A 6 gj
j 1; 2; . . . ; f ;
19
where gj is the upper bound for g^j A, and f is the number of constraints. Lower bound may be given similarly. Then the constraint is normalized assuming gj 6 0 as gj A
g^j A ÿ 1 6 0: gj
C A
subject to gj A 6 0
4.
20
The objective function C A is the total structural volume. The value of Ai has to be selected from the list fA0k g. Then consider the following optimum design problem: minimize
3.
21 j 1; 2; . . . ; f
Ai 2 fA0k g k 1; 2; . . . ; r;
22
where r is the number of members in the list. There have been many methods presented for optimization of structures with discrete variables [17]. The purpose of this section is to show eectiveness of the exact reanalysis method applied to a topology optimization of trusses. A random search algorithm is used in the following, where a new technique is presented for improving convergence to the optimal solution and for avoiding the diculty in de®ning appropriate value for the penalty terms. The most simple approach for using a random search for constrained optimization problem may be to add a penalty term of the constraints to the objective function to formulate a performance measure. In this case, however, it is very dicult to assign an appropriate value for the penalty coecient so that the magnitude of objective function is equivalent to that of the penalty term. In the proposed method, the penalty for violating the constraints is implicitly given by the threshold level for de®ning the increase or decrease of the design variables. The optimization algorithm is as follows: 1. Assign initial value of A from the given list fA0k g, and set parameters G; s and d. 2. Evaluate the objective function C A and the constraints gj A, and ®nd the maximum value G max gj . At this step, the matrix K is reconstructed j and decomposed. Note that G is negative if all the constraints are satis®ed. Assign the parameter
5. 6.
k 0:5 sG, where s moves the threshold level for increase or decrease of the variable. The parameter k is replaced with 0 or 1 if k < 0 or k > 1; respectively. Generate a uniform random number between 0 and 1 to select the member k for which cross-sectional area is to be modi®ed. Suppose there are r allowable values in the increasing order in the list fA0j g, and Ak A0i at the current step. Generate another random number p and calculate the smallest integer q that is not less than jk ÿ pj=d. Modify Ak to A0j if p < k, where j minfr ÿ i; qg, otherwise decrease Ak to A0j , where j maxfi ÿ q; 1g. Compute constraint functions for the modi®ed design, where the exact reanalysis method is to be used here. Accept the design and go to 2 if G < G and the number of the accepted designs is less than the speci®ed limit. Go to 3 if the number of trial at this stage is less than the speci®ed limit. Find the optimal design from the list of the accepted designs.
As an example, the proposed algorithm has been applied to topology optimization of plane square grids under constraints on member stress. The case of 2 2 and 5 5 are as shown in Figs. 1 and 2. The length of the members in x- or y-direction are 2.0 m, and P 200:0 kN. Stress constraints are given for all the members, where the upper bound for the absolute value is 40.0 MPa. The parameters are G 0:3, s 10:0 and d 0:1. Note that local buckling constraints are not considered here. The main diculty of topology optimization under stress constraints arises from the fact that the constraint
Fig. 2. A 5 5 plane grid.
M. Ohsaki / Computers and Structures 79 (2001) 673±679
677
Table 2 Result of optimization with 100 dierent random value parameters 22 33 44 55
u
sa
ta
Va
V min
V max
0.0025 0.005 0.007 0.01
469 288 564 1026
8278 4303 7933 14 297
0.00729698 0.130939 0.184428 0.256714
0.00692918 0.103967 0.138661 0.172274
0.0892801 0.179940 0.254913 0.364895
need not be satis®ed by the vanishing member [18]. Therefore, the optimal solution is sometimes an isolated point of the feasible region in the design variable space, and it is not useful to use the optimal solution with continuous variable [19] as the initial solution for the problem with discrete variable. The list of allowable cross-sectional areas (cm2 ) is given as fA0j g f0:01; 20; 40; 60; 80g. Note that the member with Ai 0:01 cm2 is to be removed and the stress constraint need not be satis®ed by such a member. An arti®cial upper bound u has been given for the nodal displacements to prevent the design such that the crosssectional areas of all the members are equal to 0.01 cm2 is regarded as the optimum design because obviously it has the minimum volume and no constraint should be satis®ed [20]. Optimization has been carried out for 2 2, 3 3, 4 4 and 5 5 grids, where 100 dierent initial values of the random value are generated for each case. The limit for number of modi®cation of the cross-sectional area is 5000, and the limit for the trial for generating an acceptable neighboring solution at each step is 2000. Let sa and ta denote the average numbers of steps and trial analysis before reaching the minimum value of the objective value within 5000 steps. Average, minimum and maximum values of the minimum objective function
within 5000 steps among 100 dierent random value parameters are denoted by V a , V min and V max respectively. Table 2 shows the values of u, sa , ta , V a , V min and max V for each case. Note that optimal solutions have not been reached for all the cases of dierent random value parameters. It has been con®rmed that an appropriate value of u leads to rapid convergence to nearly optimal solutions. The values of u is a little greater than the maximum displacement of the truss with optimal topology for each case. The optimum design of 5 5 grid is as shown in Fig. 3. The cross-sectional areas are 60.0 cm2 for members 6± 10 and 80.0 cm2 for the remaining existing members. The value of V min for each case in Table 2 corresponds to the objective value of the apparent optimal solution similar with the truss in Fig. 3. It may be observed from Table 2 that the average number sa of analysis needed to decompose the stiness matrix is less than 10% of the average number ta of trial analysis that can be done by the reanalysis formula. The ratio of the CPU time for obtaining the optimal topology of the 5 5 grid by using the reanalysis method to that without reanalysis is 0.2250. Therefore, computational cost may be dramatically reduced by using the proposed method of reanalysis. Next, optimal topology has been found for the rectangular 10 5 grid as shown in Fig. 4, where P 100:0 kN. Other parameters are same as those for the square grids. Optimization has been carried out for 100 dierent values of random value parameter. The results for
Fig. 3. The optimal solution for the 5 5 plane grid.
Fig. 4. A 10 5 plane grid.
678
M. Ohsaki / Computers and Structures 79 (2001) 673±679
Fig. 5. The optimal solution for the 10 5 plane grid.
the arti®cial upper-bound displacement u 0:02 cm are sa 3732, ta 40 690, V a 0:443299, V min 0:253592 and V max 0:651226. The solution corresponding to V min that coincides with the obvious solution is as shown in Fig. 5. Although the ratio of V a to V min is not suciently small for this case, it is important to obtain the optimal topology with simple algorithm that does not involve diculty in assignment of the parameters.
4. Conclusion The exact reanalysis method has been reformulated based on the general formula of inverting a modi®ed matrix. A simple expression that leads to exact estimate of the responses has been derived for modi®cation of cross-sectional areas of a truss. In this method, only the responses of the initial truss is computed for the equivalent load corresponding to the unit axial force of the member of which the cross-sectional areas should be modi®ed, and complicated matrix operation that are necessary for estimating the inverse of the stiness matrix of the modi®ed structure are not needed. It has been shown that the formulation similar to that by Melosh and Luik [7] can be derived from the formula by Woodbury [16] for multi-rank modi®cation. Note that the formulation based on the strains should be used for estimating the responses after addition of a member, because the stress of the member with null cross-sectional area cannot be calculated while the strain of such a member can be obtained from the displacements of the nodes connected to the member. A random search algorithm has been presented for trusses with discrete list of cross-sectional areas. The penalty for the violation of the constraints is incorporated as the threshold level for deciding increase or decrease of the cross-sectional areas for a given random number. Therefore, there is no diculty in de®ning the penalty coecients considering the balance of the ob-
jective function and the penalty for the constraint violation. It has been shown that the exact reanalysis is very useful for the case where the cross-sectional area of only one member is modi®ed at each step for generating a neighboring solution. Eectiveness of the proposed method has been demonstrated through examples of square and rectangular square grids. The diculty due to discontinuity in the stress constraints has been successfully avoided by considering the cross-sectional areas as discrete variables. It has been shown that an appropriate value of arti®cial upper bound displacement leads to a rapid convergence to the optimal topology under stress constraints.
References [1] Noor AK, Lowder HE. Structural reanalysis via a mixed method. Comput Struct 1975;5:9±12. [2] Storaasli OO, Sbieszczanski J. On the accuracy of the taylor approximation for structural resizing. AIAA J 1983;21:1571±80. [3] Haftka RT, G urdal Z, Kamat MP. Elements of structural optimization. London: Kluwer Academic Publishers, Dordrecht, 1990. [4] Fleury C. Ecient approximation concepts using second order information. Int J Num Meth Engng 1989;28:2041± 58. [5] Fuchs MB. Nth-order stiness sensitivities in structural analysis. Struct Optimizat 1993;5:207±12. [6] Haftka RT, Nachlas JA, Watson LT, Rizzo T, Desai R. Two-point constraint approximation in structural optimization. Comp Meth Appl Mech Engng 1987;60:289± 301. [7] Melosh RJ, Luik R. Multiple con®guration analysis of structures. Proc ASCE 1968;94(ST11):2581±96. [8] Argyris JH, Roy R. General treatment of structural modi®cation. Proc ASCE 1972;98(ST2):465±92. [9] Kirsch U, Rubinstein MF. Structural reanalysis by iteration. Comput Struct 1972;2:497±510. [10] Kirsch U. Ecient sensitivity analysis for structural optimization. Comp Meth Appl Mech Engng 1994;117: 143±56. [11] Kirsch U, Liu S. Exact structural reanalysis by a ®rst-order reduced basis approach. Struct Optimizat 1995;10:153±8. [12] Kirsch U, Liu S. Structural reanalysis for general layout modi®cation. Proc 6th AIAA/NASA/ISSMO Symposium on Multidisciplinary Analysis and Optimization. AIAA 1996. p. 1142±9. [13] Henderson HV, Searle R. On deriving the inverse of a sum of matrices. SIAM Rev 1981;23:53±60. [14] Kavlie D, Graham H, Powell GH. Ecient reanalysis of modi®ed structures. Proc ASCE 1971;97(ST1):377±92. [15] Sherman J, Morrison WJ. Adjustment of an inverse matrix corresponding to changes in one element of a given matrix. Ann Math Statist 1950;21:124±7. [16] Woodbury MA. Inverting modi®ed matrices. Memorandum Report 42, Princeton: NJ, Statistical Research Group, 1950.
M. Ohsaki / Computers and Structures 79 (2001) 673±679 [17] Arora JS, Elwakeil OA, Chahande AI, Hsier CC. Global optimization methods for engineering applications. Struct Optimizat 1995;9:137±59. [18] Sved G, Ginos Z. Structural optimization under multiple loading. Int J Mech Sci 1968;10:803±5.
679
[19] Cheng G, Guo X. e-relaxed approach in structural topology optimization. Struct Optimizat 1997;13:258± 66. [20] Ohsaki M. Genetic algorithm for topology optimization of trusses. Comput Struct 1995;57(2):219±25.