Particle swarm approach for structural design optimization

Particle swarm approach for structural design optimizationR.E.Pereza,*,K.Behdinanba Institute for Aerospace Studies,University of Toronto,4925Dufferin Street,Toronto,Ont.,Canada M3H 5T6bDepartment of Aerospace Engineering,Ryerson University,350Victoria Street,Toronto,Ont.,Canada M5B 2K3Received 17February 2006;accepted 31October 2006Available online 6March 2007AbstractThis paper presents in detail the background and implementation of a particle swarm optimization algorithm suitable for constraint structural optimization tasks.Improvements,effect of the different setting parameters,and functionality of the algorithm are shown in the scope of classical structural optimization problems.The effectiveness of the approach is illustrated by three benchmark structural optimization tasks.Results show the ability of the proposed methodology to find better optimal solutions for structural optimization tasks than other optimization algorithms.Ó2007Elsevier Ltd.All rights reserved.Keywords:Structural optimization;Truss structures;Particle swarms;Constrained optimization;Soft Computing;Stochastic optimization1.IntroductionOver the past decade a number of optimization algo-rithms have been used extensively in structural optimiza-tion task,from gradient-based algorithms using convex and continuous design spaces,to non-gradient probabilis-tic-based search algorithms widely applied for global and non-convex design exploration.In the latter category many of these algorithms have been developed by mimicking nat-ural phenomena like simulated annealing [1],genetic algo-rithms [2],and bacterial foraging [3]among others.Recently,a family of optimization algorithms has been developed based on the simulation of social interactions among members of a specific species looking for food sources.From this family,the two most promising algo-rithms are ant colony optimization,which is based on the pheromone pathways used by ants to guide other ants in a colonies towards food sources [4],and particle swarm optimization or PSO,which is based on the social behav-iour reflected in flock of birds,bees,and fish that adjusttheir physical movements to avoid predators,and to seek the best food sources [5].The PSO algorithm was first proposed by Kennedy and Eberhart [6].It is based on the premise that social sharing of information among members of a species offers and evo-lutionary advantage.A number of advantages with respect to other algorithms make PSO an ideal candidate to be used in optimization tasks.The algorithm is robust,well suited to handle non-linear,non-convex design spaces with discontinuities.It can handle continuous,discrete and inte-ger variable types with ease.As compared to other robust design optimization methods PSO is more efficient,requir-ing fewer number of function evaluations,while leading to better or the same quality of results [7,8].In addition,it easiness of implementation makes it more attractive as it does not required specific domain knowledge information,internal transformation of variables or other manipula-tions to handle constraints.Furthermore,it is a popula-tion-based algorithm,so it can be efficiently parallelized to reduce the total computational effort.Recently,the PSO has been proved useful on diverse engineering design applications such as logic circuit design [9],control design [10–12],and power systems design [13–15]among others.Applications in structures had been done in the area of structural shape optimization [16,17],0045-7949/$-see front matter Ó2007Elsevier Ltd.All rights reserved.doi:10.1016/pstruc.2006.10.013*Corresponding author.E-mail addresses:rperez@utias.utoronto.ca (R.E.Perez),kbehdina@acs.ryerson.ca (K.Behdinan)./locate/compstrucComputers and Structures 85(2007)1579–1588and in topology optimization[18],with promising results in such structural design applications.This paper focuses on the implementation and application of PSO for structural optimization.The general PSO methodology as well as the enhancement of the basic algorithm is introduced. Application of the algorithm to classical constraint struc-tural optimization problems is shown.The development of the paper is as follows:Section2presents the general formulation of the particle swarm optimization approach. In Section3analytical properties of the algorithm are discussed and improvement of the basic algorithms are pre-sented.Section4presents different structural optimization case studies.They are used to analyze the behaviour and sensitivity of the PSO parameters,and demonstrate the effectiveness of the approach infinding optimal structural optimization solutions.Finally,the paper closes with some concluding remarks.2.Particle swarm algorithm2.1.Mathematical formulationThe particle swarm process is stochastic in nature;it makes use of a velocity vector to update the current posi-tion of each particle in the swarm.The velocity vector is updated based on the‘‘memory’’gained by each particle, conceptually resembling an autobiographical memory,as well as the knowledge gained by the swarm as a whole. Thus,the position of each particle in the swarm is updated based on the social behaviour of the swarm which adapts to its environment by returning to promising regions of the space previously discovered and searching for better positions over time.Numerically,the position x of a particle i at iteration k+1is updated as shown in Eq.(1)and illustrated in Fig.1.x i kþ1¼x ikþv ikþ1D tð1Þwhere v ikþ1is the corresponding updated velocity vector,and D t is the time step value[19].Throughout the present work a unit time step is used.The velocity vector of each particle is calculated as shown in Eq.(2),v i kþ1¼wv ikþc1r1ðp ikÀx ikÞD tþc2r2ðp gkÀx ikÞD tð2Þwhere v ikis the velocity vector at iteration k,r1and r2rep-resents random numbers between0and1;p ikrepresents thebest ever particle position of particle i,and p gkcorrespondsto the global best position in the swarm up to iteration k.The remaining terms are problem dependent parameters;for example,c1and c2represent‘‘trust’’parameters indicat-ing how much confidence the current particle has in itself(c1or cognitive parameter)and how much confidence ithas in the swarm(c2or social parameter),and w is the iner-tia weight.This later term plays an important role in thePSO convergence behaviour since it is employed to controlthe exploration abilities of the swarm.It directly impactsthe current velocity,which in turn is based on the previoushistory of rge inertia weights allow for widevelocity updates allowing to globally explore the designspace,while small inertia values concentrate the velocityupdates to nearby regions of the design space.2.2.Algorithm descriptionBased on the particle and velocity updates as explainedabove,the algorithm is constructed as follows:(1)Initialize a set of particles positions x iand velocitiesv irandomly distributed throughout the design spacebounded by specified limits.(2)Evaluate the objective function values fðx ikÞusing thedesign space positions x ik.(3)Update the optimum particle position p ikat currentiteration(k)and global optimum particle position p gk.(4)Update the position of each particle using its previousposition and updated velocity vector as specified inEqs.(1)and(2).(5)Repeat steps2–4until the stopping criteria ismet.For the current implementation the stopping cri-teria is defined based on the number of iterationsreached.Eqs.(3)and(4)are used to obtain the initial position andvelocity vectors randomly distributed throughout thedesign space for each particle of the swarm.The term rin both equations represents a random number between0and1,while x min and x max represent the design variableslower and upper bounds respectively.x i¼x minþrðx maxÀx minÞð3Þv i¼x minþrðx maxÀx minÞD tð4Þ2.3.Algorithm analysis2.3.1.PSO search direction and stochastic step sizeLet us replace the PSO velocity update equation(2)intoEq.(1)to get the following expression:x ikþ1¼x ikþwV ikþc1r1ðp ikÀx ikÞD tþc2r2ðp gkÀx ikÞD tD tð5Þ1580R.E.Perez,K.Behdinan/Computers and Structures85(2007)1579–1588Factorizing the cognitive and social terms from the above equation we obtain the following general equation:x i k þ1¼x i k þwV i k D t þðc 1r 1þc 2r 2Þc 1r 1p i kþc 2r 2p g k c 1r 1þc 2r 2Àx i k ð6Þwhich follows the general gradient line-search form x i k þ1¼^x i k þa p k such that:^x i k ¼x i k þwV i k D t a ¼c 1r 1þc 2r 2 p k ¼c 1r 1p i kþc 2r 2p g k1122Àx i kð7Þwhere a can be regarded as a stochastic step size,and p k as a stochastic search direction.The stochastic step size will be limited only by the selection of social and cognitive param-eters;knowing that r 1,r 22[0,1]it will belong to the inter-val [0,c 1+c 2]with a mean value of ½c 1þc 2.Similarly,the stochastic search direction will belong to the interval½Àx i k ;c 1p i kþc 2p g k12Àx i k ,which is dependent not only on the so-cial and cognitive parameters,but also in the individual and global best positions.2.3.2.PSO convergenceBy replacing the particle velocity equation (2)into the position equation (1)and re-arranging the position term,the general form for the i th particle position equation at iteration k +1can be obtained asx i k þ1¼x i k ð1Àc 1r 1Àc 2r 2ÞþwV i k D t þc 1r 1p i k þc 2r 2p g kð8ÞSimilarly,re-arranging the position term in the particle velocity equation (2)leads to Vi k þ1¼Àx i kðc 1r 1þc 2r 2ÞD t þwV ik þc 1r 1p i k D t þc 2r 2p g kD tð9ÞEqs.(8)and (9)can be combined and written in matrixform asx i k þ1V i k þ1"#¼1Àc 1r 1Àc 2r 2w D t Àðc 1r 1þc 2r 2ÞD tw 2435x i k V i k þc 1r 1c 2r 2c 1r 1D t c 2r 2D t "#p i kp g k ð10Þwhich can be considered as a discrete-dynamic system rep-resentation for the PSO algorithm where [x i ,V i ]T is the state subject to an external input [p i ,p g ]T ,and the first and second matrices correspond to the dynamic and input matrices,respectively.Assuming there is no external excitation in the dynamic system so [p i ,p g ]T is constant (other particles do not find better positions),then a convergent behaviour could be maintained.If such is the case,as the iterations k !1then ½x i k þ1;V i k þ1 T ¼½x i k ;V i k T which reduces the dynamic metric form to0 ¼Àðc 1r 1þc 2r 2Þw D t Àðc 1r 1þc 2r 2ÞD tw À12435x i k V i k þc 1r 1c 2r 2c 1r 1D t c 2r 2D t "#p i k p g kð11Þwhich is true only when V i k ¼0and both x i k and p i k coincide with p g k .Note that such position is not necessarily a local or global minimizer,it is instead an equilibrium point.Such point,however,will improve towards the optimum if there is external excitation in the dynamic system,so better indi-vidual-best and global-best positions are found from the optimization process.2.3.3.PSO stabilityThe dynamic matrix eigenvalues derived from Eq.(10)can provide additional insight in the system stability and dynamic behaviour.The dynamic matrix characteristic equation is derived ask 2Àðw Àc 1r 1Àc 2r 2þ1Þk þw ¼0ð12Þwhere the eigenvalues are given ask 1;2¼1þw Àc 1r 1Àc 2r 2Æffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffið1þw Àc 1r 1Àc 2r 2Þ2À4w q 2ð13ÞThe necessary and sufficient condition for stability of a dis-crete-dynamic system is that all eigenvalues (k )derived from the dynamic matrix lies inside a unit circle around the origin on the complex plane,so j k i =1,...,n j <1.An analysis on the derived eigenvalues (12)leads to the follow-ing set of stability conditions for the PSO:c 1r 1þc 2r 2>0ðc 1r 1þc 2r 2Þ2Àw <1w <1ð14ÞFrom the above conditions the inertia term lies in the inter-val ðc 1r 1þc 2r 2ÞÀ1<w <1from which the upper limit for c 1r 1+c 2r 2can be derived as (c 1r 1+c 2r 2)<4.Knowing that r 1,r 22[0,1]the following parameter selection heuristic can be derived:0<ðc 1þc 2Þ<4ðc 1þc 2Þ2À1<w <1ð15ÞTherefore,if w ,c 1,and c 2are selected with the heuristics specified in (15),and meet the conditions specified in (14)the system has guaranteed convergence to an equilibrium point as discussed in Section 2.3.2.3.Algorithm improvements 3.1.Inertia weight updateDue to the importance of the inertia weight in control-ling the global/local search behaviour of the PSOR.E.Perez,K.Behdinan /Computers and Structures 85(2007)1579–15881581algorithm,a dynamic improvement has proven useful by forcing an initial global search with a high inertia weight (w%1)and subsequently narrowing down the algorithm exploration to feasible areas of the design space by decreas-ing its value towards local search values(w<0.5).Two main approaches had been used to deal with such dynamic update.In Shi and Eberhart[19],a dynamic variation of inertia weight is proposed by linearly decreasing w with each algorithmic iteration as shown in Eq.(16)w kþ1¼w maxÀw maxÀw mink maxkð16Þwhile in Fourie and Groenwold[16]a dynamic decrease of w value has been suggested based on a fraction multiplier (k w)as shown in Eq.(17),when no improvements had been made for a predefined number of consecutive design iterationsw kþ1¼k w w kð17Þ3.2.Violated design points redirectionAn additional improvement introduced by Venter and Sobieszczanski-Sobieski[20]when dealing with violated constraint particles is included in the present optimization effort.The method restricts the velocity vector of a violated particle to a usable feasible direction that would reduce the objective function while pointing back towards feasible regions of the design space.The modification is made to the velocity vector of a particle i specified in Eq.(2)when there is one or more violated constraints by re-setting the velocity vector of the particle according to Eq.(18)which include only the obtained particle self-information of the best point and the current best point in the swarm.This strategy will introduce the new velocity vector pointing towards the feasible regions of the design space as shown in Fig.2.A new position for the violated constraint parti-cles is defined by using Eq.(1)with the velocity vector modified asv i kþ1¼c1r1ðp ikÀx ikÞD tþc2r2ðp gkÀx ikÞD tð18ÞThe velocity of particle i at iteration k+1is thus only influenced by the best point found so far for that particle,and the current best point in the swarm.If both of these best points are feasible,the new velocity vector will point back to a feasible region of the design space.Otherwise, the new velocity vector will point to a region of the design space that resulted in smaller constraint violations.The re-sult is to have the violated particle move back towards the feasible region of the design space,or at least closer to its boundary,in the next design iteration.3.3.Constraints handlingSimilar to other stochastic optimization methods,the PSO algorithm is defined for unconstrained problems.To accommodate the inclusion of constraints,a parameter-less adaptive penalty scheme is used as seen in Eq.(19).It uses the swarm information,such as the average of the objective function and the level of violation of each constraint during each iteration,in order to define different penalties for dif-ferent constraints:f0ðxÞ¼fðxÞif x is feasible;fðxÞþP mi¼1k i g iðxÞotherwise8<:ð19ÞThe penalty parameter is defined at each iteration ask i¼j fðxÞjg iðxÞPl¼1½ g lðxÞð20Þwhere f(x)is the objective function,m is the number of con-straints,g i(x)is a specific constraint value(with violated constraints having values larger than zero), fðxÞis the aver-age of the objective function values in the current swarm, and g iðxÞis the violation of the l th constraint averaged over the current population.The above formulation distributes the penalty coeffi-cients in a way that those constraints which are more diffi-cult to be satisfied will have a relatively higher penalty coefficient.Such distribution is achieved by making the i th coefficient proportional to the average of the violation of the i th constraint by the elements of the current popula-tion.An individual in the swarm whose i th violation equals the average of the i th violation in the current population for all i,has a penalty equal to the absolute value of the1582R.E.Perez,K.Behdinan/Computers and Structures85(2007)1579–1588average objective function of the population.Similarly,the average of the objective function equals fðxÞþj fðxÞj.4.Numerical examplesThe effectiveness of the implemented PSO algorithm on structural optimization is shown through the use of four classical truss optimization examples.4.1.Six-node trussThefirst example considers a well-known test problem corresponding to a10-bar truss non-convex optimization as describe by Sunar and Belegundu[21]and shown in Fig.3.In this problem the cross-sectional area for each of the10members in the structure are being optimized towards the minimization of total weight.The cross-sec-tional area varies between0.1in.2and35.0in.2.Con-straints are specified in terms of stress and displacement of the truss members.The allowable stress for each mem-ber is25,000psi for both tension and compression,and the allowable displacement on the nodes is±2in.,in the x and y directions.The density of the material is0.1lb/in.3, Youngs modulus is E=104ksi and vertical downward loads of100kips are applied at nodes2and4.In total the problem has a variable dimensionality of10and con-straint dimensionality of34(10tension constraints,10 compression constraints,and12displacement constraints). The effect of the different PSO parameters,specifically the effect of social and cognitive parameters and different dynamic inertia weight updates are presented,followed by general result of the sizing problem.4.1.1.Social and cognitive parameters variation effectA more detailed analysis of the effect of c1and c2param-eters is necessary to determine the sensitivity of the param-eters in the overall optimization procedure.The effect of the cognitive(c1)and social(c2)parameters on the general optimization history can be seen in Fig.4where objective function values vs.iteration number is shown.The results are representative of more than20trials for each tested case,with afixed inertia weight value of0.875,and the same initial particles and velocity values.The parameters were varied according to the formula c2=4Àc1,where c1was varied in the interval[4,0].It can be seen clearly that when only cognitive(c1=4,c2=0)or only social values (c1=0,c2=4)are used the resulting history converges very fast,within thefirst10iterations,but do not improve after the initial convergence indicating that the algorithm reaches a local suboptimal and is not able to escape from it,due to the lack of information exchange either from social or cognitive sources respectively.When higher emphasis is placed on the social exchange of information (c1=1,c2=3)the algorithm is able to gradually converge to better regions of the design space,but again if a local optimum is found it cannot escape from it.This fact indi-cates an overshooting of the algorithm where particle posi-tions are updated based mostly on social information not promoting the exchange of individual information which could be useful when an individual particle has found a bet-ter region of the design space.When the cognitive and social parameters are in balance,(c1=c2=2,so the mean of Eq.(2)is equal to1)or with slightly higher cognitive value(c1=3,c2=1),solutions converge near the global optimum solution;these two situations are shown by the gradual decrease of objective function values and subse-quent improvements over time.Notice that a slightly higher value for the cognitive parameter results in better solutions.This result is due to the fact that individuals con-centrate more in their own search regions avoiding over-shooting the best design space region,but promoting some global exchange of information to make the swarm to point toward the best regions.4.1.2.Inertia weight variation effectTo analyze the effect of the inertia weight variation, three different test cases were analyzed.Thefirst one con-siders afixed mean inertia weight of the typical w interval [1.00,0.35],the second one uses a linearly varying update over the same interval,while the third one uses a dynamic update as described by Eq.(4)with a fraction multiplier k w=0.975.Results are given in Fig.5.Notice thatfixedR.E.Perez,K.Behdinan/Computers and Structures85(2007)1579–15881583inertia weight allow for a fast convergence(within50iter-ations)to sub-optimal regions of the design space while lin-ear and dynamic weight updates converge slowly,due to the initial global search nature force by higher weight update values,and gradually converges to optimal regions of the design space approximately after iteration120;no clear distinction in the solution is apparent between both varying inertia weight updates,but the dynamic update tend to converge faster to the optimal solution which could be useful if a specific converge value is prescribed as the stop criteria of the algorithm.It is important to note that the optimal selection of the PSO parameters are in general problem-dependent.How-ever,the obtained results confirms the expectations as derived from the theoretical analysis(see i.e.[22–25])and experiments(see i.e.[26–28])regarding the sensitivity and behaviour of such tuning parameters.As long the stability conditions presented in(14)are met,it is observed that maintain an approximately equal weighting between the cognitive and social parameters in the interval of1.4–2lead to the optimal convergence for the PSO.Similarly,afixed inertia weight makes the PSO to stagnated with a high rate of convergence in local sub-optimal values,while dynami-cally decreasing such parameter allows for better explora-tion of the design space with similar rates of convergence.4.1.3.Structural optimization resultsTable1shows the PSO best and worst results of20inde-pendent runs for the10-bar truss problem.Based on the above inertia study,results were obtained using a dynamic inertia weight update to provide the best global/local opti-mization trade.Other published results,for the same case, are summarized as well in Tables1and2.They are foundTable1Optimization results for the10-bar trussTruss area PSO best PSO worst Gellatly andBerke[30]Schimit andMiura[29]Ghasemi et al.[37]Schimit andFarshi[34]Dobbs andNelson[31]0133.50033.50031.35030.57025.73033.43230.500 020.1000.1000.1000.3690.1090.1000.100 0322.76633.50020.03023.97024.85024.26023.290 0414.41713.30415.60014.73016.35014.26015.428 050.1000.1000.1400.1000.1060.1000.100 060.1000.1000.2400.3640.1090.1000.210 077.534 6.82638.3508.5478.7008.3887.649 0820.46718.93522.21021.11021.41020.74020.980 0920.39218.81422.06020.77022.30019.69021.818 100.1000.1000.1000.3200.1220.1000.100 Weight5024.215176.275112.005107.305095.655089.005080.00 Table2Optimization results for the10-bar truss(continuation)Truss area Rizzi[32]Haug andArora[39]Haftka andGrdal[40]Adeli andKamal[35]El-Sayed andJang[36]Galante[38]Memari andFuladgar[33]0130.73130.03130.52031.2832.96630.44030.561 020.1000.1000.1000.100.1000.1000.100 0323.93423.27423.20024.6522.79921.79027.946 0414.73315.28615.22015.3914.14614.26013.619 050.1000.1000.1000.100.1000.1000.100 060.1000.5570.5510.100.7390.4510.100 078.5427.4687.4577.90 6.3817.6287.907 0820.95421.19821.04021.5320.91221.63019.345 0921.83621.61821.53019.0720.97821.36019.273 100.1000.1000.1000.100.1000.1000.100 Weight5061.605060.9205060.805052.005013.244987.004981.1 1584R.E.Perez,K.Behdinan/Computers and Structures85(2007)1579–1588using different optimization approaches including gradient based algorithms both unconstrained [29],and constrain-ted [30–33],structural approximation algorithms [34],con-vex programming [35],non-linear goal programming [36],and genetic algorithms (global search optimization)[37,38].We can see that the PSO provide good results as compared with other methods for this problem.The opti-mal solutions found by the PSO meet all constraints and have only two active constraints including the displace-ments at nodes 3and 6.4.2.25-bar trussThe second example considers the weight minimization of a 25-bar transmission tower as described by Schmit and Fleury [41]and shown in Fig.6.The design variables are the cross-sectional area for the truss members,which are linked in eight member groups as shown in Table 3.Loading of the structure is presented in Table 4.Con-straints are imposed on the minimum cross-sectional area of each truss (0.01in.2),allowable displacement at each node (±0.35in.),and allowable stresses for the members in the interval [À40,40]ksi.In total this problem has a var-iable dimensionality of 8and a constraint dimensionality of 84.Table 5present the best and worst result of 20indepen-dent runs of the PSO after 200iterations with a swarm of 40individuals for the 25-bar truss optimization task.This table presents as well results for the same optimization task from different research efforts.Clearly,the PSO yield excel-lent solutions for both its best and worst results.The best PSO solution provides the lowest structural weight of any reported solution,while the worst PSO solution is in agree-ment with those reported in the literature.The optimal solutions obtained by the PSO have the same active con-straints as reported in other references as follows:the dis-placements at nodes 3and 6in the Y -direction for both load cases and the compressive stresses in members 19and 20for the second load case.Also,there are no con-straint violations.4.3.72-bar trussThe third example comprises the optimization of a four-story 72-bar space truss shown in Fig.7.The dimen-sions are indicated as well in the figure.The optimizationTable 3Truss members area grouping for the 25-bar truss Group Truss members A11A22–5A36–9A410,11A512,13A614–17A718–21A822–25Table 4Nodal loads for the 25-bar truss Node F x F y F z 11000À10,000À10,00020À10,000À10,0003500006600Table 5Optimization results for the 25-bar truss Area group PSO best PSO worst Zhou and Rozvany [42]Haftka and Grdal [40]Erbatur et al.[43]Zhu [44]Wu and Chow [45]A10.10000.10000.0100.0100.10.10.1A2 1.02270.9895 1.987 1.987 1.2 1.90.5A3 3.4000 3.4000 2.994 2.991 3.2 2.6 3.4A40.10000.10000.0100.0100.10.10.1A50.1000 3.40000.0100.012 1.10.1 1.5A60.63990.69990.6840.6830.90.80.9A7 2.0424 1.9136 1.677 1.6790.4 2.10.6A8 3.4000 3.4000 2.662 2.664 3.4 2.6 3.4Weight485.33534.84545.16545.22493.80562.93486.29R.E.Perez,K.Behdinan /Computers and Structures 85(2007)1579–15881585。