- Email: [email protected]
An algorithm for optimal design of steel frame structures
Applied Numerical Mathematics 47 (2003) 503–514 www.elsevier.com/locate/apnum
An algorithm for optimal design of steel frame structures ✩ V. Pereyra a,∗ , D. Lawver a , J. Isenberg b a Applied Sciences Division, Weidlinger Associates, 4410 El Camino Real #110, Los Altos, CA 94022, USA b Weidlinger Associates, 375 Hudson St., New York, NY 10014, USA
1. Introduction Structural optimization has been a topic of interest for over 100 years (see [14]). The 1960s showed simultaneously a flourishing due to the new finite element and mathematical programming algorithms, and a decay due to the lack of sufficient computer power to carry even the simplest optimization tasks. This led in the 1970s to approximate structural synthesis techniques, which coupled with advances in computers, allowed the solution of more significant problems, although mostly of academic interest. Vanderplaats states in his survey that for small and medium size problems the optimization costs are small relative to the total computational effort, but as the number of design variables and constraints increase, conventional methods will not work, and new algorithms need to be developed. He goes on to suggest the use of massive parallel or distributed computing systems as a possible practical solution. The work of Chan et al. [5,4] is the closest to the topic of this paper. The main differences are that full finite element analyses are used here, with no simplifying assumptions, and that it is assumed that derivatives of the objective function and constraints are not easily available because of the use of a commercial finite element code. This is a practical assumption, and it is chosen to recreate more closely the environment of an actual structural design shop, thus making the techniques more attractive to practicing designers. In [12] the problem of optimal design of steel structures has been considered by the authors using an optimization algorithm that was derived from the Optimality Criteria of Chan [4], which is also described in [5] and other publications. In that paper, reduction of the structure’s weight was chosen as the optimization goal, subject to load and other constraints and the algorithm was applied to some small structures. In this paper an improved, faster and more robust algorithm is presented and applied to larger structures of practical interest. Gravity, wind loads, and stress constraints according to more complete construction code requirements are considered. Other reasonable building constraints considered include grouping the ✩
This material is partly based upon work supported by the National Science Foundation under Grant No 0124874.
* Corresponding author.
E-mail address: [email protected] (V. Pereyra). 0168-9274/$30.00 2003 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/S0168-9274(03)00089-8
504
V. Pereyra et al. / Applied Numerical Mathematics 47 (2003) 503–514
members, so that an unique AISC W-section is assigned to all the members of a group. Also, the size or variety of W-sections that can be used for certain groups can be limited, thus reducing the number of parameters and the variety of sections. The authors understand that the total weight of the structure is not necessarily the main component of construction cost. One advantage of this algorithm is that many other different design goals and constraint types can be used with the same methodology. The resulting algorithm couples a commercial structural analysis package with an unconstrained optimization one. The nonlinear constraints are incorporated through an exact penalty method. Given an initial design that establishes the topology of the structure, the cross-sectional area (A) of groups of members are taken as the design parameters. Let NT be the total number of members, N the number of different groups of members, Nc the number of constraints, and Ns the number of currently active (i.e., violated) constraints. The weight factor for the group of members s is defined as ρLs j , (1) ws = j
where Ls j is the length of member j in group s, and ρ is the density of steel. The total weight is the goal functional: Φ(A) =
N
ws As ,
(2)
i=s
where As is the cross-sectional area of group member s. Thus the problem is: min Φ(A),
(3)
A
subject to gs (A) gsu ,
s = 1, . . . , Nc ,
(4)
where the gsu are the upper bounds for the constraints which, for example, can represent limit drifts of the different stories when subject to wind or earthquake loads, or they can represent limit values of the parameters, or they can represent limits on allowable stresses as defined in the AISC steel construction code [9]. The commercial finite element code SAP2000 [13] is used for analyzing the elastic response of the structure. The necessary information is extracted from its output in order to evaluate the goal functional and the constraint residuals.
2. Stress constraints The AISC ASD approach has been chosen for this work since this is still the most common method used and therefore it will be familiar to a larger group of readers. However, the method can be easily adapted to LRFD or metric approaches. The strength of members subject to combined stresses is determined according to the following provisions (see AISC 1989): Cmx fbx Cmy fby fa + + 1.0, )F Fa (1 − fa /Fex )Fbx (1 − fa /Fey by fbx fby fa + + 1.0. 0.6Fy Fbx Fby
(5) (6)
V. Pereyra et al. / Applied Numerical Mathematics 47 (2003) 503–514
505
Table 1 Definition of quantities used fa frc A Fa Cmx , Cmy fbx bm sx , sy Fbx Fby Fy
actual maximum axial stress of member = frc/A maximum force in member cross-sectional area allowable maximum axial stress 0.85 actual maximum bending stress = bm /sx maximum bending moment elastic section moduli allowable maximum bending stress = 0.66Fy = 23.76 ksi (163.71 MPa) (used for compact section) allowable maximum bending stress = 0.75Fy = 23.76 ksi (163.71 MPa) 36 ksi (248 MPa) for steel
Fex
12π 2 E = 149231(rx /K ∗ l)2 23(Klb /rx )2
Fey
12π 2 E = 149231(ry /K ∗ l)2 23(Klb /rx )2
l E K rx ry Ix Iy frcs bm2s , bm3s
unbraced length in the plane of bending modulus of elasticity of steel = 29000 ksi (199 810 MPa) effective length factor√ radius of gyration = Ix /A radius of gyration = Iy /A moment of inertia M from W-member properties moment of inertia M from W-member properties maximum force for member s bending moments for member s
The units used are inches and ksi. A description of all the quantities can be found in Table 1. The allowable maximum axial stress is calculated according to the following: CC = π 2 ∗ E/Fy = 126.1. slmax = max(Kx ∗ l/rx , Ky ∗ l/ry ),
(7)
If slmax > CC then, , Fa = Fex
(8)
else Fa =
[1 − (K l/rx )2 /(2CC2 )]Fy . 1.667 + 0.375K l/(CC rx ) − 0.125(K l/rx )3 /CC3
(9)
In other words, if slmax > CC then Fa =
12 ∗ π 2 ∗ E = 149331/ slmax2 , 23 ∗ slmax2
(10)
else (1 − 0.5 ∗ slmcc2 )Fy , 1.667 + 0.375 ∗ slmcc −0.125 ∗ slmcc3 where slmcc = slmax /CC . Putting these quantities into (5,6) one gets, for each member s, Fa =
(11)
506
V. Pereyra et al. / Applied Numerical Mathematics 47 (2003) 503–514
g1s (A) =
|frcs | 0.03577|bm2s | 0.03148|bm3s | + , + ) ) As ∗ Fa sx (1 − fa /Fex sy (1 − fa /Fey
0.0463|frcs | 0.0421|bm2s | 0.0421|bm3s | + + , g2s (A) = As sx sy
(12)
for a total of Nc = idrift +2 ∗ N constraints, where idrift is the number of drift constraints.
3. An exact penalty method In this section an exact penalty method is introduced for the constrained optimization problem above, which transforms it into an unconstrained one. The idea is to add to the weight functional a penalty term that is a function of the violated constraints. The larger the violation the larger this term sbecomes. There are many kinds of penalty and barrier methods. By defining, Ψ (A) = Φ(A) + κ N s=1 | ress (A)|, the exact penalty method from Gill et al. [6] can be stated as: minA Ψ (A),
(13)
where ress (A) = gsu − gs (A) 0,
(14)
correspond to the Ns active constraints. If A∗ is a Kuhn–Tucker point of the constrained problem, then for any κ sufficiently large, it can be proven that A∗ is also a minimum point of this unconstrained problem and viceversa. Other penalty methods require that the penalty parameter be driven to infinity in order to obtain convergence to the solution of the original problem. As the penalty parameter increases, the resulting problems become increasingly ill-conditioned and thus harder to solve. That is not the case with the exact penalty method chosen here. The resulting unconstrained minimization problem is solved using a derivative free iterative method due to Powell [12], and modified by Brent [2], who incorporated several important improvements to enhance efficiency and reliability. This algorithm remains as one of the best direct search methods yet written and it compares well with quasi-Newton methods that require derivatives. One of its main advantages is that the number of constraints is not an important factor in the optimization, contrary to what happens with other methods. The method searches along conjugate directions and it produces the unique minimum in N steps for a quadratic positive definite functional of N variables. For a general nonlinear functional the process is iterated until a convergence criterium is met. For the current problem, the appeal of this method comes from the fact that the derivatives of the goal functional with respect to the design parameters would be very hard to extract from a black box commercial code and they would be expensive to estimate by finite differences. In addition, the exact penalty functional itself is non-differentiable at the solution. If A0 is an initial design and ui are the columns of an orthogonal matrix, then the basic procedure can be described as follows: • For i = 1, . . . , N , calculate the step βˆi that minβi Ψ (Ai−1 + βi ui ) along the search direction ui . Then define Ai = Ai−1 + βˆi ui .
V. Pereyra et al. / Applied Numerical Mathematics 47 (2003) 503–514
507
• Replace un by uˆ n = An − A0 . • Calculate the step βˆ along the direction uˆ n by minβ Ψ (x 0 + β uˆ n ), and replace x 0 by x 0 + βˆ uˆ n . • Iterate until stopping criterion is satisfied. A periodic re-start solves the potential problem of linear dependence among the search directions that occurs if a β = 0. In ill-conditioned cases, or if the one-dimensional linear searches fail, random directions are also incorporated.
4. Design optimization Most intelligent search codes as the one selected here, or any of the stochastic search algorithms such as those used in [1,3,7,8], are not appropriate for problems with a large number of variables, since they require a very large number of function evaluations to converge. In order to extend the applicability of search methods to problems with many parameters and constraints a block decomposition approach is considered. In this approach the parameters are partitioned into blocks and the sub-problems are optimized for the parameters within the active block, while the remaining parameters are kept fixed. As soon as the optimization of a block is completed a master copy of the parameters is updated in a Gauss–Seidel type iteration. A full passage through all the blocks is called a sweep. Multiple sweeps are desirable and one can regulate the budget of the iterations in each block in order to make the full process efficient. In order to improve efficiency even further, the constraints associated with the active groups of members are controlled. If they are all satisfied but there are still external active constraints, then the iteration for this block is terminated. The rational being that constraints are given a large weight in the penalty functional and thus it does not pay to try to diminish the weight of a feasible block when there are violated constraints elsewhere. Observe that this block scheme is also suitable for parallel execution in a distributed network, as shown in [10], where sufficient conditions for convergence of the process are also stated. Thus the process is applicable to large scale engineering problems, an important feature. The optimization is an iterative process. Given an initial description of the structure, including a set of W-sections, the optimization code will trigger a structural analysis. With the results of this analysis, the goal functional and the constraints will be calculated, and a new set of cross-sectional areas Ai for the members in the active block will be generated. These new cross-sectional areas need not correspond to any existing W-section. A full description of the W-sections requires not only the cross-sectional area, but also the moments of inertia. Additionally, there is only a finite number of section types, and thus, in order to apply a continuous optimization process one needs to perform some additional steps. For each Ai values of the moments of inertia are calculated by linear interpolation within the W family of the previous value of that group. Thus, taking a generic cross-sectional area A, where the previous corresponding member was in family Wn , if A falls between the cross-sectional areas of sections Wnm and Wn,m+1 (where the sections have been ordered within each family Wn by increasing values of the cross-sectional area), then Ix = Ixnm +
An,m+1 − Anm nm A − A , Ixn,m+1 − Ixnm
(15)
508
V. Pereyra et al. / Applied Numerical Mathematics 47 (2003) 503–514
and similarly for Iy . For values of A outside the range of cross-sectional areas linear extrapolation is used instead. Once a trio (A, Ix , Iy ) has been calculated, the whole data base of W-sections (or the subset associated with that particular member) is searched, in order to find the closest W-section to it, always rounding above to avoid selecting too weak a section. As a measure of distance the square of the Euclidean norm of the relative deviations is used 2 2 2 (16) d = A − Aij /Aij + Ix − Ixij /Ixij + Iy − Iyij /Iyij . Ideally one would want the initial design to satisfy all the constraints, i.e., to be feasible (overdesigned). Then the optimization algorithm would reduce the size and weight of the structural members until a constraint is violated. Unfortunately this may be impractical, and as a matter of fact, the initial design for the 20 story multibay building below, produced by an experienced Structural Engineer, did violate several of the construction code constraints, since it was not optimized on purpose. In order to get closer to an ideal initial feasible design, a pre-processing analysis has been devised that takes advantage of the locality of the code constraints. In this pre-processing, the values of the member’s cross sectional areas are modified in proportion to the size of their code constraint functions. This brings the structure closer to compliance and produces a better starting design to optimize. This is done at the beginning of each block optimization. Also, the best calculated solution so far is used as a starting guess for the next block. In short, the algorithm works just as well with unfeasible as with feasible designs and no convergence difficulties have been observed. Also, the choice of the penalty parameter has been straightforward.
5. Eight story building with drift and design constraints A small example from Pereyra et al. [11], consisting of a one bay wide and one bay deep (10 × 20 feet2 ), eight story building is considered first. The columns of every two floors are grouped together and so are the beams. Wind drift (maximum allowable displacement of the top floor is 2 inches), stress constraints and two different load combinations are considered. Thus NT = 68,
N = 8,
Nc = 18.
(17)
The finite element package SAP2000 is used for structural analysis. This code calculates and reports maximum and minimum stresses and bending moments at a number of points for each member. SAP2000 is a fully interactive, graphic oriented package, but within the optimization algorithm is used as a batch code. Interaction with it occurs through the input and output files *.s2k and *.out, which are parsed and modified as necessary. Any other similar code could be used for this purpose, but it is emphasized that an important part of the approach is to show the possibility of coupling such a commercial black-box to an optimization code, even if the structural analysis package does not provides derivatives with respect to the design parameters. All the structural information that is needed by the code is extracted automatically from the *.s2k input file, thus eliminating any redundant and potentially conflicting input by the user. The calculation of the constraint values (12) is done for each member of a group, and for every load combination. Finally, the worst case (largest) is picked as the actual value of the constraint. Observe that members may be further subdivided into elements. In that case the check is at the element level, but the total length of the member has to be used. This information is also extracted from the *.s2k file.
V. Pereyra et al. / Applied Numerical Mathematics 47 (2003) 503–514
509
Table 2 Eight story building: Lightest feasible design Design
10 Weight 25627.57
Group
W-sect.
Cond
AXL
Ms
Mw
SUM
1
W14X68
1 2
0.81 0.79
0.17 0.20
0.00 0.00
0.99 0.98
2
W14X53
1 2
0.77 0.74
0.17 0.19
0.00 0.00
0.94 0.93
3
W14X53
1 2
0.92 0.70
0.00 0.20
0.00 0.00
0.92 0.90
4
W14X34
1 2
0.76 0.73
0.22 0.25
0.00 0.00
0.98 0.98
5
W14X61
1 2
0.00 0.00
0.35 0.41
0.00 0.00
0.35 0.41
6
W16X26
1 2
0.00 0.00
0.80 0.94
0.00 0.00
0.80 0.94
7
W6X9
1 2
0.01 0.01
0.62 0.73
0.00 0.00
0.63 0.74
8
W16X26
1 2
0.01 0.01
0.48 0.56
0.00 0.00
0.49 0.57
The engineer is required to design the structure in the usual way and produce a valid input file (*.s2k file). In general, it is advisable to be conservative and to produce a design that satisfies the constraints. In this example, the four first blocks correspond to the columns, while the last four correspond to the beams. A value of K = 1.2 was used in this calculation, but the value obtained by the structural analysis code could have been used as well. The lightest of the 17 different feasible designs found by the code using only 91 analyses is shown in Table 2. The total computer time on a 300 MHz Pentium 3 machine was just a few minutes for this small structure. For every group member the three individual terms in each of the two constraints (3) and their sum are also listed. This sum should be less than one for every group in order to have a feasible design. In this way the engineer can rapidly judge if the response of a particular group of members is dominated by axial forces or by bending moments. Observe that, since there are 268 different W-sections, the total number of all possible designs is 268!/260!, which is approximately 2.4 × 1019 . Obviously, this large number of possible combinations can be trimmed by limiting the types of sections that are reasonable for each member, but still, the combinatorics of the problem makes the exhaustive analysis of all reasonable structures impossible, even for such a small building. It also makes the eye-balling of an optimal design by trial and error a very hard, labor intensive task. In Fig. 2 the weights for the different feasible designs obtained by the algorithm are plotted. Observe that these are just some sampled values of a function of 8 variables, and should not be construed as the evolution of a one-dimensional minimization. The algorithm makes large excursions in an area surrounding the initial design, and thus offers the engineer an intelligent sampling of possible improved designs. The program does not eliminate the need for the designer. It just provides a computerized aid to generate alternative designs with smaller and smaller goal functionals, so that the engineer can be inspired
510
V. Pereyra et al. / Applied Numerical Mathematics 47 (2003) 503–514
Fig. 1. An eight story building: Front view.
Fig. 2. Eight story building: Feasible designs.
V. Pereyra et al. / Applied Numerical Mathematics 47 (2003) 503–514
511
by them to create some unexpectedly economic and good looking design. For the lightest feasible design the weight has been reduced by more than 30% with respect to the initial feasible design.
6. A multibay 20 story building A more realistic design of a multibay steel frame building that includes floors and walls is now considered. The initial design was done by a professional structural engineer. The building consists of 20 stories, with a story height of 13 –6 (4.11 m). Each floor has a layout of 6 bays wide by 4 bays deep, with each bay measuring 28 (8.53 m) square. Columns are numbered 1–7 in the E–W direction, and A–E in the N–S direction. Moment frames are placed along the perimeter, as well as along lines 3 and 5. All columns are oriented with their major axis in the N–S direction, except for columns A, E 2–6, which are rotated 90 deg, aligning their major axis in the E–W direction. The members are grouped into 11 different groups, and 9 of them are further divided into two story subdivisions, for a total of 92 different groups, as shown in Table 3. The floors are modeled as 3 × 3 shell element bays, with 6 (15 cm) standard modulus of elasticity concrete. The windows are modeled as minimum weight (1 lb/cu ft) 1/4 (0.635 cm) thick elements with the same modulus as the concrete. Loads are applied as indicated in Table 4. Some of the members Table 3 Twenty story building grouping of members Group numbers 1–10 11–20 21–30 31–40 41–50 51–60 61–70 71–80 81 82 83–92
Type
Column lines
Group name
Initial sizes
Corner columns E–W face columns N–S face columns, 2-way moment N–S face columns 1-way moment Inner moment columns Inner gravity columns N–S fascia beams E–W fascia beams Interior gravity beams Interior gravity girders Interior moment girders
A1,A7,E1,E7 B1,B7,C1,C7,D1,D7 A3,A5,E3,E5
corner col ew col nsm
W14 W14 W14
A2,A4,A6,E2,E4,E6
col nsg
W14
col inm col ing bm ns bm ew bm int gird grav gird m
W14 W14 W16,W24 W16,W24 W16,W24 W16,W24 W16,W24
B3,B5,C3,C5,D3,D5 B2,B4,B6,C2,C4,C6,D2,D4,D6 Lines A,E Lines 1,7 Lines B,C,D Lines 2,4,6 Lines 3,5
Table 4 Loads for twenty story building Dead load: Live load: Wind load:
Selfweight 720 plf 50 psf 26.0 psf 36.4 psf 42.9 psf 52.0 psf 61.1 psf
(75 psf from slab) perimeter line load (0–30 ft) (30–50 ft) (50–100 ft) (100–200 ft) (200–300 ft)
512
V. Pereyra et al. / Applied Numerical Mathematics 47 (2003) 503–514
Fig. 3. Plan view of 20 story building.
are subdivided into several elements. The wind loads are all applied to either the S or W face. All steel used was standard A36, but any other grade of steel could have been chosen. Observe that now the number of all possible combinations of W-sections for this structure is a staggering 268!/176!, which approximately equals 10215 . Considering that each analysis of this structure (on a Pentium III, 300 MHz computer) takes around 20 minutes, this strongly indicates the need for an intelligent search. In Fig. 4 the evolution of the weights for the feasible designs found are shown. The improvement, from the heaviest to the lightest design in the two runs is almost 10%, a significant reduction. The initial design, with a weight of 4 152 823 pounds (1 885 382 kg), is not in this figure since it was not feasible. In order to give some idea of the economic impact and where this design fits in current practice, the weight is corrected to include floor beams, which were not present in the initial design. From this corrected total weight, the number of pounds per square foot of construction is calculated. This correction is only approximated and pessimistic, since the average weight of the current structural E–W beams was used in order to calculate the weight of two additional floor beams per bay in that direction. Current weight of design #55 = 3 754 890 pounds (1 704 720 kg). Weight corrected to include 2 E–W floor beams per bay = 4 453 803 pounds (2 022 027 kg). Square footage: 376 320 (34 961 square m); Steel per sft = 11.84 p/sf (57.84 kg/square m).
V. Pereyra et al. / Applied Numerical Mathematics 47 (2003) 503–514
513
Fig. 4. Twenty story building feasible designs.
The weight of steel per square feet for this design is somewhat below of what is obtained by current best design practice, 12 p/sf (58.64 kg/square m).
7. Conclusions A procedure has been described that combines a finite element structural design code with an optimization one, in order to reduce the weight of a given structure subject to gravity and wind loads, under maximum drift and stress constraints. The algorithm uses a block approach and an exact penalty method for optimization within the blocks, that seems to be sufficiently robust for application to realistic, large scale structures. The unconstrained optimization problems are solved with a code that does not require derivatives, while the structural analyses are performed with SAP2000, a commercial code. The combined system has the ability of bringing an uncomplying design into a compliant one, and then work from there to reduce the weight while satisfying all the imposed constraints. A restart option permits to further optimize the design if that is desired. It presents to the structural engineer as output a list of the feasible designs found ordered by increasing weight. An additional auxiliary post-processing code can be used to generate SAP2000 *.s2k input files of selected feasible designs for visualization, checking or further processing. Thus, even if the optimal design is not reached, because of lack of time, useful results will be produced. Also, no a priori decision is made about what would be the most adequate design, to cover the realistic case when not all relevant building, economic or aesthetic goals or constraints have been included in the optimization, thus showing once again that the ultimate responsible is the engineer and that this is only a computer aid. Any other structural analysis program with SAP2000 capabilities could be used to produce a similar optimal design code. Because of the way in which the algorithm is crafted it would be straightforward
514
V. Pereyra et al. / Applied Numerical Mathematics 47 (2003) 503–514
to transform it into a parallel code, for which different blocks of variables would be processed asynchronously and simultaneously on different processors of a cluster of computers, thus improving the turn around time for large buildings.
Acknowledgement The 20 story building was designed by Eric Hoyt, Weidlinger Associates, New York. We thank Eric for several illuminating conversations and we also thank Tian-Fang Jing for directing the New York effort and for reviewing the manuscript.
References [1] J.S. Arora (Ed.), Guide to Structural Optimization, ASCE Manuals and Reports on Engineering Practice, Vol. 90, ASCE, New York, 1997. [2] R.P. Brent, Algorithms for Minimization Without Derivatives, Prentice-Hall, Englewood Cliffs, NJ, 1973. [3] C. Camp, S. Pezeshk, G. Cao, Optimized design of two-dimensional structures using a genetic algorithm, J. Struct. Engrg. (1998) 551–559. [4] C.-M. Chan, Optimal stiffness design to limit static and dynamic wind responses of tall steel buildings, Engrg. J. (1998) 94–105. [5] C.-M. Chan, D.E. Grierson, A.N. Sherbourne, Automatic optimal design of tall steel building frameworks, J. Struct. Engrg. ASCE 121 (5) (1995) 838–847. [6] P.E. Gill, W. Murray, M.H. Wright, Practical Optimization, Academic Press, New York, 1981. [7] M.-W. Huang, J.S. Arora, Optimal design of steel structures using standard sections, Struct. Optim. 14 (1997) 24–35. [8] N.S. Khot, V.B. Venkayya, L. Berke, Optimal structural design with stability constraints, Internat. J. Numer. Methods Engrg. (1976) 1097–1114. [9] Manual of Steel Construction. Allowable Stress Design, 9th Edition, AISC, Chicago, IL, 1989. [10] V. Pereyra, Asynchronous distributed solution of large scale nonlinear inversion problems, Appl. Numer. Math. 30 (1999) 31–40. [11] V. Pereyra, D. Lawver, J. Isenberg, Optimal design of steel frame structures, Appl. Numer. Math. 40 (2002) 59–77. [12] M.J.D. Powell, An efficient method for finding the minimum of a function of several variables without calculating derivatives, Comp. J. 7 (1964) 155–162. [13] SAP2000, Computers and Structures, Berkeley, CA, 1998. [14] G.N. Vanderplaats, Structural optimization: Where we’ve been and where we’re going, in: S. Hernandez, C.A. Brebbia (Eds.), Computer Aided Optimum Design of Structures V, Comp. Mechanics, W.I.T. Press, Southampton, 1997, pp. 45– 54.
An algorithm for optimal design of steel frame structures ✩ V. Pereyra a,∗ , D. Lawver a , J. Isenberg b a Applied Sciences Division, Weidlinger Associates, 4410 El Camino Real #110, Los Altos, CA 94022, USA b Weidlinger Associates, 375 Hudson St., New York, NY 10014, USA
1. Introduction Structural optimization has been a topic of interest for over 100 years (see [14]). The 1960s showed simultaneously a flourishing due to the new finite element and mathematical programming algorithms, and a decay due to the lack of sufficient computer power to carry even the simplest optimization tasks. This led in the 1970s to approximate structural synthesis techniques, which coupled with advances in computers, allowed the solution of more significant problems, although mostly of academic interest. Vanderplaats states in his survey that for small and medium size problems the optimization costs are small relative to the total computational effort, but as the number of design variables and constraints increase, conventional methods will not work, and new algorithms need to be developed. He goes on to suggest the use of massive parallel or distributed computing systems as a possible practical solution. The work of Chan et al. [5,4] is the closest to the topic of this paper. The main differences are that full finite element analyses are used here, with no simplifying assumptions, and that it is assumed that derivatives of the objective function and constraints are not easily available because of the use of a commercial finite element code. This is a practical assumption, and it is chosen to recreate more closely the environment of an actual structural design shop, thus making the techniques more attractive to practicing designers. In [12] the problem of optimal design of steel structures has been considered by the authors using an optimization algorithm that was derived from the Optimality Criteria of Chan [4], which is also described in [5] and other publications. In that paper, reduction of the structure’s weight was chosen as the optimization goal, subject to load and other constraints and the algorithm was applied to some small structures. In this paper an improved, faster and more robust algorithm is presented and applied to larger structures of practical interest. Gravity, wind loads, and stress constraints according to more complete construction code requirements are considered. Other reasonable building constraints considered include grouping the ✩
This material is partly based upon work supported by the National Science Foundation under Grant No 0124874.
* Corresponding author.
E-mail address: [email protected] (V. Pereyra). 0168-9274/$30.00 2003 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/S0168-9274(03)00089-8
504
V. Pereyra et al. / Applied Numerical Mathematics 47 (2003) 503–514
members, so that an unique AISC W-section is assigned to all the members of a group. Also, the size or variety of W-sections that can be used for certain groups can be limited, thus reducing the number of parameters and the variety of sections. The authors understand that the total weight of the structure is not necessarily the main component of construction cost. One advantage of this algorithm is that many other different design goals and constraint types can be used with the same methodology. The resulting algorithm couples a commercial structural analysis package with an unconstrained optimization one. The nonlinear constraints are incorporated through an exact penalty method. Given an initial design that establishes the topology of the structure, the cross-sectional area (A) of groups of members are taken as the design parameters. Let NT be the total number of members, N the number of different groups of members, Nc the number of constraints, and Ns the number of currently active (i.e., violated) constraints. The weight factor for the group of members s is defined as ρLs j , (1) ws = j
where Ls j is the length of member j in group s, and ρ is the density of steel. The total weight is the goal functional: Φ(A) =
N
ws As ,
(2)
i=s
where As is the cross-sectional area of group member s. Thus the problem is: min Φ(A),
(3)
A
subject to gs (A) gsu ,
s = 1, . . . , Nc ,
(4)
where the gsu are the upper bounds for the constraints which, for example, can represent limit drifts of the different stories when subject to wind or earthquake loads, or they can represent limit values of the parameters, or they can represent limits on allowable stresses as defined in the AISC steel construction code [9]. The commercial finite element code SAP2000 [13] is used for analyzing the elastic response of the structure. The necessary information is extracted from its output in order to evaluate the goal functional and the constraint residuals.
2. Stress constraints The AISC ASD approach has been chosen for this work since this is still the most common method used and therefore it will be familiar to a larger group of readers. However, the method can be easily adapted to LRFD or metric approaches. The strength of members subject to combined stresses is determined according to the following provisions (see AISC 1989): Cmx fbx Cmy fby fa + + 1.0, )F Fa (1 − fa /Fex )Fbx (1 − fa /Fey by fbx fby fa + + 1.0. 0.6Fy Fbx Fby
(5) (6)
V. Pereyra et al. / Applied Numerical Mathematics 47 (2003) 503–514
505
Table 1 Definition of quantities used fa frc A Fa Cmx , Cmy fbx bm sx , sy Fbx Fby Fy
actual maximum axial stress of member = frc/A maximum force in member cross-sectional area allowable maximum axial stress 0.85 actual maximum bending stress = bm /sx maximum bending moment elastic section moduli allowable maximum bending stress = 0.66Fy = 23.76 ksi (163.71 MPa) (used for compact section) allowable maximum bending stress = 0.75Fy = 23.76 ksi (163.71 MPa) 36 ksi (248 MPa) for steel
Fex
12π 2 E = 149231(rx /K ∗ l)2 23(Klb /rx )2
Fey
12π 2 E = 149231(ry /K ∗ l)2 23(Klb /rx )2
l E K rx ry Ix Iy frcs bm2s , bm3s
unbraced length in the plane of bending modulus of elasticity of steel = 29000 ksi (199 810 MPa) effective length factor√ radius of gyration = Ix /A radius of gyration = Iy /A moment of inertia M from W-member properties moment of inertia M from W-member properties maximum force for member s bending moments for member s
The units used are inches and ksi. A description of all the quantities can be found in Table 1. The allowable maximum axial stress is calculated according to the following: CC = π 2 ∗ E/Fy = 126.1. slmax = max(Kx ∗ l/rx , Ky ∗ l/ry ),
(7)
If slmax > CC then, , Fa = Fex
(8)
else Fa =
[1 − (K l/rx )2 /(2CC2 )]Fy . 1.667 + 0.375K l/(CC rx ) − 0.125(K l/rx )3 /CC3
(9)
In other words, if slmax > CC then Fa =
12 ∗ π 2 ∗ E = 149331/ slmax2 , 23 ∗ slmax2
(10)
else (1 − 0.5 ∗ slmcc2 )Fy , 1.667 + 0.375 ∗ slmcc −0.125 ∗ slmcc3 where slmcc = slmax /CC . Putting these quantities into (5,6) one gets, for each member s, Fa =
(11)
506
V. Pereyra et al. / Applied Numerical Mathematics 47 (2003) 503–514
g1s (A) =
|frcs | 0.03577|bm2s | 0.03148|bm3s | + , + ) ) As ∗ Fa sx (1 − fa /Fex sy (1 − fa /Fey
0.0463|frcs | 0.0421|bm2s | 0.0421|bm3s | + + , g2s (A) = As sx sy
(12)
for a total of Nc = idrift +2 ∗ N constraints, where idrift is the number of drift constraints.
3. An exact penalty method In this section an exact penalty method is introduced for the constrained optimization problem above, which transforms it into an unconstrained one. The idea is to add to the weight functional a penalty term that is a function of the violated constraints. The larger the violation the larger this term sbecomes. There are many kinds of penalty and barrier methods. By defining, Ψ (A) = Φ(A) + κ N s=1 | ress (A)|, the exact penalty method from Gill et al. [6] can be stated as: minA Ψ (A),
(13)
where ress (A) = gsu − gs (A) 0,
(14)
correspond to the Ns active constraints. If A∗ is a Kuhn–Tucker point of the constrained problem, then for any κ sufficiently large, it can be proven that A∗ is also a minimum point of this unconstrained problem and viceversa. Other penalty methods require that the penalty parameter be driven to infinity in order to obtain convergence to the solution of the original problem. As the penalty parameter increases, the resulting problems become increasingly ill-conditioned and thus harder to solve. That is not the case with the exact penalty method chosen here. The resulting unconstrained minimization problem is solved using a derivative free iterative method due to Powell [12], and modified by Brent [2], who incorporated several important improvements to enhance efficiency and reliability. This algorithm remains as one of the best direct search methods yet written and it compares well with quasi-Newton methods that require derivatives. One of its main advantages is that the number of constraints is not an important factor in the optimization, contrary to what happens with other methods. The method searches along conjugate directions and it produces the unique minimum in N steps for a quadratic positive definite functional of N variables. For a general nonlinear functional the process is iterated until a convergence criterium is met. For the current problem, the appeal of this method comes from the fact that the derivatives of the goal functional with respect to the design parameters would be very hard to extract from a black box commercial code and they would be expensive to estimate by finite differences. In addition, the exact penalty functional itself is non-differentiable at the solution. If A0 is an initial design and ui are the columns of an orthogonal matrix, then the basic procedure can be described as follows: • For i = 1, . . . , N , calculate the step βˆi that minβi Ψ (Ai−1 + βi ui ) along the search direction ui . Then define Ai = Ai−1 + βˆi ui .
V. Pereyra et al. / Applied Numerical Mathematics 47 (2003) 503–514
507
• Replace un by uˆ n = An − A0 . • Calculate the step βˆ along the direction uˆ n by minβ Ψ (x 0 + β uˆ n ), and replace x 0 by x 0 + βˆ uˆ n . • Iterate until stopping criterion is satisfied. A periodic re-start solves the potential problem of linear dependence among the search directions that occurs if a β = 0. In ill-conditioned cases, or if the one-dimensional linear searches fail, random directions are also incorporated.
4. Design optimization Most intelligent search codes as the one selected here, or any of the stochastic search algorithms such as those used in [1,3,7,8], are not appropriate for problems with a large number of variables, since they require a very large number of function evaluations to converge. In order to extend the applicability of search methods to problems with many parameters and constraints a block decomposition approach is considered. In this approach the parameters are partitioned into blocks and the sub-problems are optimized for the parameters within the active block, while the remaining parameters are kept fixed. As soon as the optimization of a block is completed a master copy of the parameters is updated in a Gauss–Seidel type iteration. A full passage through all the blocks is called a sweep. Multiple sweeps are desirable and one can regulate the budget of the iterations in each block in order to make the full process efficient. In order to improve efficiency even further, the constraints associated with the active groups of members are controlled. If they are all satisfied but there are still external active constraints, then the iteration for this block is terminated. The rational being that constraints are given a large weight in the penalty functional and thus it does not pay to try to diminish the weight of a feasible block when there are violated constraints elsewhere. Observe that this block scheme is also suitable for parallel execution in a distributed network, as shown in [10], where sufficient conditions for convergence of the process are also stated. Thus the process is applicable to large scale engineering problems, an important feature. The optimization is an iterative process. Given an initial description of the structure, including a set of W-sections, the optimization code will trigger a structural analysis. With the results of this analysis, the goal functional and the constraints will be calculated, and a new set of cross-sectional areas Ai for the members in the active block will be generated. These new cross-sectional areas need not correspond to any existing W-section. A full description of the W-sections requires not only the cross-sectional area, but also the moments of inertia. Additionally, there is only a finite number of section types, and thus, in order to apply a continuous optimization process one needs to perform some additional steps. For each Ai values of the moments of inertia are calculated by linear interpolation within the W family of the previous value of that group. Thus, taking a generic cross-sectional area A, where the previous corresponding member was in family Wn , if A falls between the cross-sectional areas of sections Wnm and Wn,m+1 (where the sections have been ordered within each family Wn by increasing values of the cross-sectional area), then Ix = Ixnm +
An,m+1 − Anm nm A − A , Ixn,m+1 − Ixnm
(15)
508
V. Pereyra et al. / Applied Numerical Mathematics 47 (2003) 503–514
and similarly for Iy . For values of A outside the range of cross-sectional areas linear extrapolation is used instead. Once a trio (A, Ix , Iy ) has been calculated, the whole data base of W-sections (or the subset associated with that particular member) is searched, in order to find the closest W-section to it, always rounding above to avoid selecting too weak a section. As a measure of distance the square of the Euclidean norm of the relative deviations is used 2 2 2 (16) d = A − Aij /Aij + Ix − Ixij /Ixij + Iy − Iyij /Iyij . Ideally one would want the initial design to satisfy all the constraints, i.e., to be feasible (overdesigned). Then the optimization algorithm would reduce the size and weight of the structural members until a constraint is violated. Unfortunately this may be impractical, and as a matter of fact, the initial design for the 20 story multibay building below, produced by an experienced Structural Engineer, did violate several of the construction code constraints, since it was not optimized on purpose. In order to get closer to an ideal initial feasible design, a pre-processing analysis has been devised that takes advantage of the locality of the code constraints. In this pre-processing, the values of the member’s cross sectional areas are modified in proportion to the size of their code constraint functions. This brings the structure closer to compliance and produces a better starting design to optimize. This is done at the beginning of each block optimization. Also, the best calculated solution so far is used as a starting guess for the next block. In short, the algorithm works just as well with unfeasible as with feasible designs and no convergence difficulties have been observed. Also, the choice of the penalty parameter has been straightforward.
5. Eight story building with drift and design constraints A small example from Pereyra et al. [11], consisting of a one bay wide and one bay deep (10 × 20 feet2 ), eight story building is considered first. The columns of every two floors are grouped together and so are the beams. Wind drift (maximum allowable displacement of the top floor is 2 inches), stress constraints and two different load combinations are considered. Thus NT = 68,
N = 8,
Nc = 18.
(17)
The finite element package SAP2000 is used for structural analysis. This code calculates and reports maximum and minimum stresses and bending moments at a number of points for each member. SAP2000 is a fully interactive, graphic oriented package, but within the optimization algorithm is used as a batch code. Interaction with it occurs through the input and output files *.s2k and *.out, which are parsed and modified as necessary. Any other similar code could be used for this purpose, but it is emphasized that an important part of the approach is to show the possibility of coupling such a commercial black-box to an optimization code, even if the structural analysis package does not provides derivatives with respect to the design parameters. All the structural information that is needed by the code is extracted automatically from the *.s2k input file, thus eliminating any redundant and potentially conflicting input by the user. The calculation of the constraint values (12) is done for each member of a group, and for every load combination. Finally, the worst case (largest) is picked as the actual value of the constraint. Observe that members may be further subdivided into elements. In that case the check is at the element level, but the total length of the member has to be used. This information is also extracted from the *.s2k file.
V. Pereyra et al. / Applied Numerical Mathematics 47 (2003) 503–514
509
Table 2 Eight story building: Lightest feasible design Design
10 Weight 25627.57
Group
W-sect.
Cond
AXL
Ms
Mw
SUM
1
W14X68
1 2
0.81 0.79
0.17 0.20
0.00 0.00
0.99 0.98
2
W14X53
1 2
0.77 0.74
0.17 0.19
0.00 0.00
0.94 0.93
3
W14X53
1 2
0.92 0.70
0.00 0.20
0.00 0.00
0.92 0.90
4
W14X34
1 2
0.76 0.73
0.22 0.25
0.00 0.00
0.98 0.98
5
W14X61
1 2
0.00 0.00
0.35 0.41
0.00 0.00
0.35 0.41
6
W16X26
1 2
0.00 0.00
0.80 0.94
0.00 0.00
0.80 0.94
7
W6X9
1 2
0.01 0.01
0.62 0.73
0.00 0.00
0.63 0.74
8
W16X26
1 2
0.01 0.01
0.48 0.56
0.00 0.00
0.49 0.57
The engineer is required to design the structure in the usual way and produce a valid input file (*.s2k file). In general, it is advisable to be conservative and to produce a design that satisfies the constraints. In this example, the four first blocks correspond to the columns, while the last four correspond to the beams. A value of K = 1.2 was used in this calculation, but the value obtained by the structural analysis code could have been used as well. The lightest of the 17 different feasible designs found by the code using only 91 analyses is shown in Table 2. The total computer time on a 300 MHz Pentium 3 machine was just a few minutes for this small structure. For every group member the three individual terms in each of the two constraints (3) and their sum are also listed. This sum should be less than one for every group in order to have a feasible design. In this way the engineer can rapidly judge if the response of a particular group of members is dominated by axial forces or by bending moments. Observe that, since there are 268 different W-sections, the total number of all possible designs is 268!/260!, which is approximately 2.4 × 1019 . Obviously, this large number of possible combinations can be trimmed by limiting the types of sections that are reasonable for each member, but still, the combinatorics of the problem makes the exhaustive analysis of all reasonable structures impossible, even for such a small building. It also makes the eye-balling of an optimal design by trial and error a very hard, labor intensive task. In Fig. 2 the weights for the different feasible designs obtained by the algorithm are plotted. Observe that these are just some sampled values of a function of 8 variables, and should not be construed as the evolution of a one-dimensional minimization. The algorithm makes large excursions in an area surrounding the initial design, and thus offers the engineer an intelligent sampling of possible improved designs. The program does not eliminate the need for the designer. It just provides a computerized aid to generate alternative designs with smaller and smaller goal functionals, so that the engineer can be inspired
510
V. Pereyra et al. / Applied Numerical Mathematics 47 (2003) 503–514
Fig. 1. An eight story building: Front view.
Fig. 2. Eight story building: Feasible designs.
V. Pereyra et al. / Applied Numerical Mathematics 47 (2003) 503–514
511
by them to create some unexpectedly economic and good looking design. For the lightest feasible design the weight has been reduced by more than 30% with respect to the initial feasible design.
6. A multibay 20 story building A more realistic design of a multibay steel frame building that includes floors and walls is now considered. The initial design was done by a professional structural engineer. The building consists of 20 stories, with a story height of 13 –6 (4.11 m). Each floor has a layout of 6 bays wide by 4 bays deep, with each bay measuring 28 (8.53 m) square. Columns are numbered 1–7 in the E–W direction, and A–E in the N–S direction. Moment frames are placed along the perimeter, as well as along lines 3 and 5. All columns are oriented with their major axis in the N–S direction, except for columns A, E 2–6, which are rotated 90 deg, aligning their major axis in the E–W direction. The members are grouped into 11 different groups, and 9 of them are further divided into two story subdivisions, for a total of 92 different groups, as shown in Table 3. The floors are modeled as 3 × 3 shell element bays, with 6 (15 cm) standard modulus of elasticity concrete. The windows are modeled as minimum weight (1 lb/cu ft) 1/4 (0.635 cm) thick elements with the same modulus as the concrete. Loads are applied as indicated in Table 4. Some of the members Table 3 Twenty story building grouping of members Group numbers 1–10 11–20 21–30 31–40 41–50 51–60 61–70 71–80 81 82 83–92
Type
Column lines
Group name
Initial sizes
Corner columns E–W face columns N–S face columns, 2-way moment N–S face columns 1-way moment Inner moment columns Inner gravity columns N–S fascia beams E–W fascia beams Interior gravity beams Interior gravity girders Interior moment girders
A1,A7,E1,E7 B1,B7,C1,C7,D1,D7 A3,A5,E3,E5
corner col ew col nsm
W14 W14 W14
A2,A4,A6,E2,E4,E6
col nsg
W14
col inm col ing bm ns bm ew bm int gird grav gird m
W14 W14 W16,W24 W16,W24 W16,W24 W16,W24 W16,W24
B3,B5,C3,C5,D3,D5 B2,B4,B6,C2,C4,C6,D2,D4,D6 Lines A,E Lines 1,7 Lines B,C,D Lines 2,4,6 Lines 3,5
Table 4 Loads for twenty story building Dead load: Live load: Wind load:
Selfweight 720 plf 50 psf 26.0 psf 36.4 psf 42.9 psf 52.0 psf 61.1 psf
(75 psf from slab) perimeter line load (0–30 ft) (30–50 ft) (50–100 ft) (100–200 ft) (200–300 ft)
512
V. Pereyra et al. / Applied Numerical Mathematics 47 (2003) 503–514
Fig. 3. Plan view of 20 story building.
are subdivided into several elements. The wind loads are all applied to either the S or W face. All steel used was standard A36, but any other grade of steel could have been chosen. Observe that now the number of all possible combinations of W-sections for this structure is a staggering 268!/176!, which approximately equals 10215 . Considering that each analysis of this structure (on a Pentium III, 300 MHz computer) takes around 20 minutes, this strongly indicates the need for an intelligent search. In Fig. 4 the evolution of the weights for the feasible designs found are shown. The improvement, from the heaviest to the lightest design in the two runs is almost 10%, a significant reduction. The initial design, with a weight of 4 152 823 pounds (1 885 382 kg), is not in this figure since it was not feasible. In order to give some idea of the economic impact and where this design fits in current practice, the weight is corrected to include floor beams, which were not present in the initial design. From this corrected total weight, the number of pounds per square foot of construction is calculated. This correction is only approximated and pessimistic, since the average weight of the current structural E–W beams was used in order to calculate the weight of two additional floor beams per bay in that direction. Current weight of design #55 = 3 754 890 pounds (1 704 720 kg). Weight corrected to include 2 E–W floor beams per bay = 4 453 803 pounds (2 022 027 kg). Square footage: 376 320 (34 961 square m); Steel per sft = 11.84 p/sf (57.84 kg/square m).
V. Pereyra et al. / Applied Numerical Mathematics 47 (2003) 503–514
513
Fig. 4. Twenty story building feasible designs.
The weight of steel per square feet for this design is somewhat below of what is obtained by current best design practice, 12 p/sf (58.64 kg/square m).
7. Conclusions A procedure has been described that combines a finite element structural design code with an optimization one, in order to reduce the weight of a given structure subject to gravity and wind loads, under maximum drift and stress constraints. The algorithm uses a block approach and an exact penalty method for optimization within the blocks, that seems to be sufficiently robust for application to realistic, large scale structures. The unconstrained optimization problems are solved with a code that does not require derivatives, while the structural analyses are performed with SAP2000, a commercial code. The combined system has the ability of bringing an uncomplying design into a compliant one, and then work from there to reduce the weight while satisfying all the imposed constraints. A restart option permits to further optimize the design if that is desired. It presents to the structural engineer as output a list of the feasible designs found ordered by increasing weight. An additional auxiliary post-processing code can be used to generate SAP2000 *.s2k input files of selected feasible designs for visualization, checking or further processing. Thus, even if the optimal design is not reached, because of lack of time, useful results will be produced. Also, no a priori decision is made about what would be the most adequate design, to cover the realistic case when not all relevant building, economic or aesthetic goals or constraints have been included in the optimization, thus showing once again that the ultimate responsible is the engineer and that this is only a computer aid. Any other structural analysis program with SAP2000 capabilities could be used to produce a similar optimal design code. Because of the way in which the algorithm is crafted it would be straightforward
514
V. Pereyra et al. / Applied Numerical Mathematics 47 (2003) 503–514
to transform it into a parallel code, for which different blocks of variables would be processed asynchronously and simultaneously on different processors of a cluster of computers, thus improving the turn around time for large buildings.
Acknowledgement The 20 story building was designed by Eric Hoyt, Weidlinger Associates, New York. We thank Eric for several illuminating conversations and we also thank Tian-Fang Jing for directing the New York effort and for reviewing the manuscript.
References [1] J.S. Arora (Ed.), Guide to Structural Optimization, ASCE Manuals and Reports on Engineering Practice, Vol. 90, ASCE, New York, 1997. [2] R.P. Brent, Algorithms for Minimization Without Derivatives, Prentice-Hall, Englewood Cliffs, NJ, 1973. [3] C. Camp, S. Pezeshk, G. Cao, Optimized design of two-dimensional structures using a genetic algorithm, J. Struct. Engrg. (1998) 551–559. [4] C.-M. Chan, Optimal stiffness design to limit static and dynamic wind responses of tall steel buildings, Engrg. J. (1998) 94–105. [5] C.-M. Chan, D.E. Grierson, A.N. Sherbourne, Automatic optimal design of tall steel building frameworks, J. Struct. Engrg. ASCE 121 (5) (1995) 838–847. [6] P.E. Gill, W. Murray, M.H. Wright, Practical Optimization, Academic Press, New York, 1981. [7] M.-W. Huang, J.S. Arora, Optimal design of steel structures using standard sections, Struct. Optim. 14 (1997) 24–35. [8] N.S. Khot, V.B. Venkayya, L. Berke, Optimal structural design with stability constraints, Internat. J. Numer. Methods Engrg. (1976) 1097–1114. [9] Manual of Steel Construction. Allowable Stress Design, 9th Edition, AISC, Chicago, IL, 1989. [10] V. Pereyra, Asynchronous distributed solution of large scale nonlinear inversion problems, Appl. Numer. Math. 30 (1999) 31–40. [11] V. Pereyra, D. Lawver, J. Isenberg, Optimal design of steel frame structures, Appl. Numer. Math. 40 (2002) 59–77. [12] M.J.D. Powell, An efficient method for finding the minimum of a function of several variables without calculating derivatives, Comp. J. 7 (1964) 155–162. [13] SAP2000, Computers and Structures, Berkeley, CA, 1998. [14] G.N. Vanderplaats, Structural optimization: Where we’ve been and where we’re going, in: S. Hernandez, C.A. Brebbia (Eds.), Computer Aided Optimum Design of Structures V, Comp. Mechanics, W.I.T. Press, Southampton, 1997, pp. 45– 54.