Shape optimization of 2D structures using simulated annealing

21 Pages·2007·2.44 MB·English

Comput.MethodsAppl.Mech.Engrg.196(2007)3279–3299
www.elsevier.com/locate/cma
Shape optimization of 2D structures using simulated annealing
Fazil O. Sonmez *
DepartmentofMechanicalEngineering,BogaziciUniversity,Istanbul,Bebek34342,Turkey
Received13September2005;receivedinrevisedform22May2006;accepted11January2007
Availableonline31March2007
Abstract
Thegoalofthisstudyistoobtaingloballyoptimumshapesfortwo-dimensionalstructuressubjecttoquasi-staticloadsandrestraints.
Forthispurposeatechniqueisproposed,usingwhichthevolume(orweight)ofastructurecanbeminimized.Theemphasisisonhow
onecan definethe shapeprecisely, andfind a shapethat accuratelyreflectsthe globally optimum shape.
As design constraints, stresses developed in the structure should not exceed the maximum allowable stress, and connectivity of the
structure should not be lost during shape changes. Optimization is achieved by a stochastic search algorithm called direct simulated
annealing(DSA),which seekstheglobalminimumthroughrandomlygeneratedconfigurations. Inordertoobtainrandomconfigura-
tions,aboundaryvariationtechniqueisused.Inthistechnique,asetofkeypointsischosenandconnectedbycubicsplinestodescribe
theboundaryofthestructure.Wheneverthepositionsofthekeypointsarechangedinrandomdirections,anewshapeisobtained.Thus,
coordinatesofthekeypointsserveasdesignvariables.Inordertoapplytheoptimizationprocedure,ageneralcomputercodewasdevel-
opedusingANSYSParametricDesignLanguage.Anumberofcaseswereexaminedtotestitseffectiveness.Theresultsshowthatthis
technique can be applied to two-dimensional shape optimization problems with high reliability even for cases where the entire free
boundaryis allowed tovary.
(cid:2)2007Elsevier B.V.All rightsreserved.
Keywords: Directsearchsimulatedannealing(DSA);Boundaryvariation;Precision;Accuracy;FEM;2Dstructures
1. Introduction sureandcentrifugalloads[69],uniaxiallyloadedplatewith
a hole [77], bridge under moving [76] and static [73] loads,
The increasing demand for lightweight, high-perfor- platescontainingcracks[37],plateswithanotch[35,59,85],
mance, low-cost structures, drives the considerable current 2Dstructuressubjecttofatigueloading [74],2Dstructures
research going on in the field of structural optimization. subjecttothermomechanicalboundaryconditions[57],2D
Studies on shape optimization [1–87] attest that shape structures containing a crack [55], various other 2D struc-
changes may lead to considerable savings in weight and tures [39,48,65,70], saw blades [23], cantilever beams
improvements in structural performance. Applications in [21,24,28,29,35,41,53,60,64,70,88], beams subject to tor-
this field include optimum shape designs of fillets [4,7,9, sional loading [49], beams with thin-walled sections [44],
12,20,25,27,30,43,50,51,53,66,71,80,85,87], cams [61], holes axially excitedbeams [72],thin-walledtubes subjecttotor-
in plates [3,7,9,10,12,19,22,25–28,32,39,41,42,60,66,73,76, sionalloading[77],shellsunderdynamicloading[67],com-
80,83,85,87] for single and multiple [68] load cases, rings ponents which are in contact [40,46,54,56,75,83], wheel
[30], torque arms [6,20,27,38,62,71,84], connecting rods profileofrailwayvehicles[81],thermalfins[17,58],axisym-
[31,53,74,79,83], brackets [6,38,62,79,84], chain links metricstructures[78],3Dstructureslikeenginebearingcap
[10,27,77], hooks [20,42], torque converter subject to pres- [11], cantilevered thick plate with a hole though the thick-
ness [16,47], aircraft turbine disk [16], control arm [63],
engine mount bracket [63], and slider bearings [82].
* Tel.:[phone];fax:[phone]. Thegeneralpurposeofshapeoptimizationistofindthe
E-mailaddress:[email protected] best profile for a structure under various constraints
0045-7825/$-seefrontmatter (cid:2)2007ElsevierB.V.Allrightsreserved.
doi:10.1016/j.cma.2007.01.019
3280 F.O.Sonmez/Comput.MethodsAppl.Mech.Engrg.196(2007)3279–3299
imposed by the requirements of the design. The optimum decreasingvalueofobjectivefunctionisiterativelycreated,
shape provides either the most efficient and effective use maygetstuckintoalocallyminimumpointotherthanthe
of material or the best performance. Accordingly, in some globally minimum solution. Its success depends on the
ofthestudiesthegoalwastominimizeweightofthestruc- choice of initial design. The usual approach is to employ
ture [10,11,16,17,20,24,25,30,31,35,38–41,48,49,51,53,62, thealgorithmmanytimesstartingfromdifferentconfigura-
69–71,77–79,84], in others to increase mechanical perfor- tions with the hope that one of the initial positions be suf-
mance, e.g. to minimize stress concentration [3,19,22,25, ficiently close to the globally minimum configuration, and
26,32,33,42,43,57,60,66,68,77,83,87], maximize fracture then to choose the lowest value as the globally minimum
strength [37,80,85], buckling strength [65], fatigue life solution. Another disadvantage is that if the starting point
[50,55,59,68], and heat flux [17,58], minimize peak contact is outside the feasible region, the algorithm may converge
stress [46,54,56,75,83], compliance [7,29,34,48,63,64,67, to a local minimum within the infeasible domain.
70,57,73,78], peak acceleration [67], and the probability If a downhill-proceeding local search algorithm is used,
of failure for brittle materials [23], and optimize dynamic an optimization procedure can only be successful if it is
behavior of structures [44,52,72,88]. usedtoimprovethecurrentdesign[4,69,83]oronlyasmall
In order to optimize a structure, some of its properties segment of the boundary is allowed to move [3,4,9,12,19,
affecting the cost function should be changed. As for a 26,28,35,39,43,53,69], or only a small number of dimen-
shape optimization problem, the boundary should be sional parameters are used to define its shape [6,26,29,30].
allowed to vary. In previous studies, boundaries of the The advantage of using these algorithms, on the other
structures were usually defined by spline curves [9,11,23– hand,istheircomputationalefficiency.Justasmallnumber
25,28,31,35,39–41,46,47,51,56–61,65,66,69–71,76,77,87,88] of iterations are needed to obtain an improved design.
passingthroughanumberofkeypoints.Positionsofthese Oneoftheoptimizationmethodsthathavebeenusedto
points thus become design variables. NURBS curves [27] obtain optimum structural designs is evolutionary struc-
defined by a polygon can also be used to define a shape. tural optimization. The approach is to gradually remove
Weights at each vertex of the polygon then become design inefficientmaterialsuntiltheshapeofthestructureevolves
variables. Some parameterized equations can be used for into an optimum. Although this method is not a global
boundary curve [26,30,34]; the parameters then control searchalgorithm,becauseofitssimplicityandeffectiveness,
the size, shape and aspect ratio of the boundary. The it has been applied to many structural optimization prob-
advantage of using these curves is that one may define lems [49,74–76,80]. This method is also suitable for topol-
the shape of the model using just a few design variables ogy optimization, where not only outer boundary but
andobtainasmoothboundary.Alternatively,inanumber also geometry of inner regions is allowed to change. There
ofstudies[3,32,43,64,68,83]allorsome ofthenodesatthe arealso similaralgorithmswhichcaneffectivelybeusedto
boundarieswereallowedtochangetheirpositions;ormag- find improved designs such as biological growth methods
nitudes of a set of fictitious loads on the boundary nodes [19,38,66] and metamorphic development methods [78].
werechanged[10];ormoresimplydimensionsofindividual In order to find the absolute minimum of an objective
parts of the structures were taken as design variables function without being sensitive to the starting position, a
[29,52]. Also, by removing or restoring materials within global optimization method has to be employed in struc-
the elements of the finite element model, the shape of the tural optimization problems. Stochastic optimization tech-
structure may be changed [20,26,49,74,75,80]. In that case, niques are quite suitable in this respect. Among their
presence of material becomes a design variable. advantages, they are not sensitive to starting point, they
During a shape optimization process, geometry of the cansearchalargesolutionspace,andtheycanescapelocal
structure may undergo substantial changes, such that it optimum points because they allow occasional uphill
can no longer be considered as a viable structure, e.g. geo- moves. The genetic algorithm (GA) and simulated anneal-
metric model may become unfeasible [20,21], stresses may ing (SA) algorithm are two of the most popular stochastic
exceed the admissible stress level [20,21,24,29–31,35,38– optimization techniques.
41,51,62,71,78,79,84], fatigue life may be shorter than the Genetic algorithms found many applications on shape
lowest allowable number of cycles [74], natural frequency optimization problems [15,21,31,41,47,48,61,64,70,57,77,
may be lower than a certain limit [24], displacements may 79]. They are based on the concepts of genetics and Dar-
betoolarge[11,16,24,35],finiteelementmeshmaybecome winian survival of the fittest [13]. The idea is to start the
too distorted [59], the area or volume may become too search process with a set of designs, called population.
large [7,29,46,65,67,72,73,88], manufacturing may become The search procedure is a simulation of the evolution pro-
too difficult [59]. In these cases, behavioral or side con- cess.Thegeneticalgorithmtransformsonepopulationinto
straints should be imposed. a succeeding population using operators of reproduction,
Locating globally optimum shape designs is a difficult recombination, and mutation. The convergence to the glo-
problem, requiring sophisticated optimization procedures. bal minimum depends on the proper choice of the design
In typical structural optimization problems, there may be parameters, rules of the reproduction and mutation.
many locally minimum configurations. For that reason, a Kirkpatricketal.[89]firstproposedsimulatedannealing
downhill-proceeding algorithm, in which a monotonically (SA) as a powerful stochastic search technique. The main
F.O.Sonmez/Comput.MethodsAppl.Mech.Engrg.196(2007)3279–3299 3281
advantageofSAisthatitisgenerallymorereliableinfinding in most of the previous studies, definition of shape was
the global optimum, i.e. the probability of locating the not precise enough to be capable of representing the best
global optimum is high even with large numbers of design possible shape. In this study, on the other hand, a proce-
variables.Applicationofsimulatedannealingtoshapeopti- dureisproposedtomonitorprecisionso that areasonable
mizationofthestructuresisquiterare.Anagnostouandhis compromisecanbemadebetweenabetterdefinitionanda
coworkers [17] used SA to design a thermal fin with mini- shorter computational time.
mummaterialuse.ShimandManoochehri[20]usedSAto Another concern in shape optimization is the accuracy,
obtainoptimaldesignsfortwo-dimensionalstructureswith i.e. how well the resulting optimum shape reflects the best
a minimum use of material. In their approach, a structure possible shape. In many of the previous studies, the goal
wasfirstdividedintoanumberofsmallfiniteelementblocks wastoimproveanintuitivelyfoundcurrentdesignthrough
andthentheseblockswererandomlyremovedorrestoredto alocalsearchalgorithmratherthantofindthebestpossible
obtainarandomlygeneratednewshape.Onedifficultywith shape.Insomeofthem,globalsearchalgorithmswereused,
thismethodisthatremovingandrestoringanelementmay but the accuracy remained questionable. In this study, a
violatemodelconnectivity.Wheneveranewconfigurationis numberofmethodsareproposedtoachievehighaccuracy.
generated its connectivity should be checked by a compli-
cated algorithm. Another difficulty is the roughness of the 2. Problem statement
resulting boundaries. If smooth boundaries are desired,
small elements should be used, resulting in long computa- Consideratwo-dimensionalstructureasinFig.1having
tionaltimes.However,theirstudydemonstratedtheprom- a geometry that can be defined by its boundary and thick-
iseofSAinsolvingshapeoptimizationproblems. ness.Thestructureissubjecttoin-planequasi-staticloads;
Theultimategoalofadesignerinusingashapeoptimiza- it is also restrained from moving at some points of the
tionalgorithmisjusttostatetheboundaryconditionsandlet boundary. Applied loads and restraints are considered as
the algorithm do some iterations without human interven- boundaryconditions.Thestructureshouldbeabletoresist
tion and automatically produce the best design. In this the loads without failure. This means the material should
respect, the previous studies had only a relative success. notyield.Besides,nopartofthestructureistoloseitscon-
Many of them relied too much on designer’s intuition, or nectiontotherestraints;i.e.thestructureshouldremainin
imposed tight restrictions on the movements of boundary, one piece. These are the two constraints imposed on the
or fixed some portions of boundary so that the number of structure. Our objective is to minimize the volume (or
possible designs among which the algorithm may search weight)ofthestructure,inotherwords,toobtainaconfig-
theoptimumonewasquitelimited.Thisstudyattemptsto uration(shape)withthemostefficientuseofmaterialwith-
come one step closer to our ultimate goal. The procedure out violating the behavioral and geometric constraints.
developed in this study may start from almost any initial The shape optimization problem for minimum weight
model thanks to the global search algorithm used. In fact can be formulated briefly as
the only interactive effort required from the designer is to
provideaninitialcoarseboundarymodelconsistingofspline Given: initial positions of key points, boundary
curvespassingthroughthegivenkeypoints,andtospecify conditions (applied forces and restraints),
theboundaryconditionsandwhichkeypointsofthebound- material properties
aryareallowedtomoveduringtheoptimizationprocess. Find: the optimum shape of the structure
Anotherrequirementforthesearchofgloballyoptimum Minimize: V=volume of the structure
configuration isthe precise definition of the optimized sys- Subject to: (i) the stress constraint (r 6r )
max allow
tem.Ifthesystemiscomplex,completedefinitionofitmay (ii) the geometric constraint (model
requirelargenumbersofdesignvariables.Ifallofthemare connectivity)
usedasoptimizationvariables,thatmeansiftheoptimized Design
system is precisely defined, long computational times are variables: the positions of the key points
needed,andalsolikelihoodoflocatingtheoptimumdesign
may become very low. If the definition of the system is
F
based on a limited number of design variables, that means 1
ifthedefinitionallowsforonlyalimitednumberofdistinct
configurations, the resulting optimum configuration may
F
not reflect the best possible design. Although using the 2
most precise definition of the system during optimization
y
isdesirable,thisisusuallynotpracticable;butifalesspre-
cisedefinitionisused,oneshouldbeawareofitsimpacton
F
how well the final design represents the best possible 3
x
design. To the author’s knowledge, there is no study on
P
structural optimization, or specifically on shape optimiza- 1
tion in which precision was a main concern. Apparently, Fig.1. Representationofa2Dshapeoptimizationproblem.
3282 F.O.Sonmez/Comput.MethodsAppl.Mech.Engrg.196(2007)3279–3299
This means that the maximum equivalent (Von Mises) DSA differs from SA basically in two aspects. Firstly,
stress of the structure, r , should be less than or equal DSA keeps a set of current configurations rather than just
max
tothemaximumallowablestress,r ,andthekeypoints one current configuration. Secondly, it always retains the
allow
which are used to define the boundaries of the structure best configuration. In a way, this property imparts a sort
should move without any loss in model connectivity. of memory to the optimization process. If a newly gener-
ated configuration is accepted, it just replaces the worst
3. Problem solution configuration.
In this study, direct search simulated annealing (DSA), 3.1. The cost function
proposedbyAlietal.[90],isadopted.Thisisanimproved
version of simulated annealing. SA gets its name from the Because the thickness of the structure is fixed and only
physicalprocesswherebythetemperatureofasolidisraised its lateral area is allowed to vary, the objective function
toameltingpoint,wheretheatomscanmovefreelyandthen tobeminimizedmaybetakenasitsareainsteadofitsvol-
slowlycooled.Themethodattemptstomodelthebehavior ume.Shapeoptimizationisaconstrainedoptimization;but
of atoms in forming arrangements in solid material during SA is only applicable to unconstrained optimization
annealing.Althoughtheatomsmoverandomly,astheirnat- problems. By integrating a penalty function for constraint
ural behavior they tend to form lower-energy configura- violationsintothecostfunction,theconstrainedoptimiza-
tions. However, this is a time-driven process. When a tion problem can be transformed into an unconstrained
material is crystallized from its liquid phase, it must be problem, for which the algorithm is suitable. A combined
cooled slowly if it is to assume its highly ordered, lowest- cost function may be constructed as
energy, perfect crystalline state. At each temperature level
during cooling, the material should reach equilibrium. As f ¼ A þc(cid:2)rmax (cid:2)1(cid:3)2: ð1Þ
the temperature decreases, the arrangement of the atoms A r
ini allow
getscloserandclosertothelowest-energystate.Thisprocess
Here, the bracketed term is defined as
continues until the temperature finally reaches freezing
point. (cid:2)r (cid:3) (cid:4)0 for r 6r ;
max (cid:2)1 ¼ max allow
Thereisananalogybetweenaphysicalannealingprocess r ðr =r Þ(cid:2)1 forr >r ;
allow max allow max allow
andanoptimizationprocess.Differentconfigurationsofthe
ð2Þ
problemcorrespondtodifferentarrangementsoftheatoms.
Thecostofaconfigurationcorrespondstotheenergyofthe where A is the area of the initial configuration; c is a
ini
system. Optimal solution corresponds to thelowest-energy weighingcoefficient,forwhichavalueabout0.9wasfound
state. Just in the same manner the atoms find their way to to be satisfactory. Any excursion into the infeasible region
build a perfect crystal structure through random move- (r >r ) results in an increase in the cost function.
max allow
ments, the global optimum is reached through a search On the other hand, if model connectivity is lost, the new
within randomlygenerated configurations. configuration is never accepted, and therefore there is no
InSA,anarbitrary initialdesign isselectedandsystem- need to calculate its cost function.
atically updated until a stopping criterion is satisfied.
Updatingisbased onan iterative procedure.In eachitera- 3.2. The boundary variation technique
tion, an arbitrary design is generated in the neighborhood
of the current configuration. If the cost function of the SA requires random generation of a new configuration
newdesignhasasmallervaluecomparedtothatofthecur- (in our case, a new shape) in each iteration. Definition of
rentdesign,itisaccepted.Thenewdesignthenreplacesthe a configuration should be made based on optimization
old one. On the other hand, if the new cost function has a variables. The algorithm tries to find the optimum values
larger value, the acceptability of the design is decided of these variables, which define the best configuration.
according to the probability of Boltzmann distribution. Shape of a 2D-structure is defined by its thickness and
The calculation of this probability depends on a tempera- boundary. Accordingly, its shape can easily be described
ture parameter, T, which is referred to as temperature, by spline curves passing through a number of key points.
because it plays a similar role in the optimization process Someofthesepointsmaybefixed,whileothersareallowed
as the temperature in the physical annealing process. The to move during optimization. As illustrated in Fig. 2,
temperature is kept constant for a number of trials and whenever the positions of the moving key points are chan-
then reduced. The rate of reduction should be slow so as gedinrandomdirectionsthroughrandomdistances,anew
not to get trapped at a locally minimum point. At initial boundary, thus a new configuration is obtained. Accord-
stages of the algorithm (at high temperatures), the proba- ingly, the x and y coordinates of the moving key points
bility of accepting worse designs is higher but later on at become optimization variables. Model connectivity is
low temperatures, this probability becomes lower and ensured by just checking whether or not the spline seg-
lower so that in the end the designs having a higher cost ments of the boundary cross over or touch each other,
are almost never accepted. which is automatically done by the FE software.
F.O.Sonmez/Comput.MethodsAppl.Mech.Engrg.196(2007)3279–3299 3283
f . The number of these configurations depends on the
S ‘
F dimension of the problem:
1
N ¼7ðnþ1Þ; ð3Þ
F wherenisthedimensionoftheproblem,i.e.thenumberof
2
design variables. Since x and y coordinates of the moving
y keypointsarethevariablesdefiningtheshapes,inourcase
n is given by
F
3 n¼2n ; ð4Þ
x m
P
1 where nm is the number of moving key points. This means
thatforeachadditionalmovingkeypoint,onemustgener-
Fig. 2. Generation of a random configuration through random move-
ate 14 more current configurations.
mentsofkeypoints.
3.4. Generation mechanism
Keypointsareallowedtomoveonlywithinaregion,S,
definedbythedesigner.Althoughsearchoftheglobalmin-
In each iteration, first, a configuration is randomly
imum without restricting the movements of the key points
selected from the set of current configurations, and then
is possible, computational time becomes unnecessarily
random displacements are given to the moving key points,
much longer. Search domain should not be too restrictive,
resulting in a new shape. If any of the moving key points
but should exclude the regions that are definitely expected
falls outside the search domain S (Fig. 2), this process is
to be far away from the boundary.
repeated until all of the moving key points fall within S.
For the direction of movement, a random real number
During movements of the key points if one segment of
between0and360,forthemovingdistancearandomnum-
the boundary crosses another one, or they touch each
ber between 0 and R is chosen. Here R is the maxi-
max max other,connectivityislost.Inthatcase,thegenerationpro-
mum moving distance. Our experience has shown that in
cess is repeated.
order to search the global optimum within a large region,
insteadofgivingsmallperturbations tothe current config-
3.5. Acceptability
urationtoobtainanewconfigurationinthenearneighbor-
hood of the current one, we should allow a large variance
Acceptability of a new configuration depends on the
in the current configurations. For this reason, the magni-
value of its cost function, f; its acceptability, A, is calcu-
tude of R is taken as 10% of the square root of the t t
max lated by
searchedarea,S;thuseachkeypointcanmoveuptoadis-
tanceofabout10%ofaveragedimension.Thismeansthat A ¼(cid:4)1 if ft 6fh; ð5Þ
the neighborhood of a current configuration where a new t expððf (cid:2)fÞ=T Þ if f >f
h t k t h
configuration is generated is initially quite large. This can
here f is the highest cost function value in A. This means
also be considered as a logical consequence of simulating h
every new design having a cost lower than the cost of the
the physical annealing process, where mobility of atoms
worst design is accepted.But,if the cost ishigher, the trial
is large at high temperatures. Also, as in the physical pro-
configuration may be accepted depending on the value of
cess,wheremobilityofatomsdecreasesasthetemperature
A. If it is greater than a number randomly generated be-
is lowered, moving distance is reduced as the temperature t
tween0and1,thetrialconfigurationisaccepted,otherwise
parameter is decreased; but the reduction scheme does
rejected.
notdirectlydependontemperature.Theconfigurationthat
Ifthetrialfunctionisaccepted,itreplacestheworstcon-
is worse than all current configurations except the worst
figuration. The best configuration, thus, always remains in
one is defined as the worse configuration, and if no
A.Ineachiteration,f andf areupdated.Athightemper-
improvementismadeontheworseshapeduringaMarkow h ‘
aturesitisunlikelythatAwillformadensecluster.Atlow
chain, R is reduced by 0.5%.
max temperatures, because the likelihood of accepting a worse
configuration is low, A will shrink to form a dense cluster.
3.3. The set of current configurations
3.6. Cooling schedule
DSA unlike ordinary SA starts with a set of N configu-
rations, A, rather than starting with only one configura- The simulated annealing process consists of first ‘‘melt-
tion. For the shape optimization problem, the initial set ing’’ the system being optimized at a high ‘‘temperature’’,
ofNdifferentshapesisselectedrandomlywithinthesearch T, then lowering the temperature slowly until the system
domain, S, without any regard to cost, but connectivity is ‘‘freezes’’andnofurtherchangesoccur.Here,temperature,
observed. Their costs are calculated and stored. The high- T, has no physical meaning; it is just a control para-
est and lowest cost function values are denoted as f and meter. Melting corresponds to the initial stage at which
h
3284 F.O.Sonmez/Comput.MethodsAppl.Mech.Engrg.196(2007)3279–3299
configurationsaregeneratedwithinthesolutiondomain,S, depends on howfast thevalue of the temperature parame-
without much regard to the cost. At high temperatures, ter is lowered. If the cooling rate is fast, the optimization
(high values of T) the probability of acceptance is high as process will probably end up with one of the higher-cost
Eq.(5)implies.Accordingly,configurationsthathaveeven local minima. If the cooling rate is slow, the reliability of
very high cost values may be accepted, just as in the phys- finding the global optimum solution will be high, but, the
ical annealing process the atoms in the melting state may optimizationprocesswilltakeexcessivelylongtime.There-
formconfigurationshavingveryhighenergy.Atlowvalues fore, the choice of the cooling schedule determines the
of temperature parameter, acceptability becomes low and effectiveness of the algorithm.
acceptance of worse configurations is then unlikely, just InDSA,atemperaturescalefactor,(a ),isspecifiedto
k+1
as the atoms become stable, and do not tend to change calculatethevalueofthetemperatureparameterinthenext
their arrangements at low temperatures. The cooling sche- Markow chain, T .
k+1
duleinSAcontrolstheconvergenceofthealgorithmtothe
T ¼a T ; ð8Þ
global minimum just as the cooling scheme in the physical kþ1 kþ1 k
annealing process controls the final microstructure. There- where T is the value of the current temperature. A value
k
fore, performance of SA depends on the cooling schedule. fora 2½a ;a (cid:3)iscalculatedusingthefollowingequa-
kþ1 min max
Inacoolingschedule,firstaninitialvalue,T ,shouldbe tion [90]:
0
specified for the temperature parameter. A scheme is then 8
a if L >L0;
required for reducing T and for deciding how many trials >< max k k
aretobeattemptedateachvalueofT.Lastly,thefreezing akþ1 ¼ ak(cid:2)ðak (cid:2)aminÞð1(cid:2)L0k(cid:2)1=LkÞ else if Lk >L0k(cid:2)1;
(or final value of the) temperature parameter is specified. >:a (cid:2)ða (cid:2)a ÞðL =L0 Þ else L 6L0 ;
max max k k k(cid:2)1 k k(cid:2)1
Initial value of the temperature parameter: Its initial
ð9Þ
value should be high enough to allow nearly all trials to
beaccepted.Thisprovidescomplete‘‘melting’’attheinitial where L0 is the actual number of trials executed in the kth
k
stages of the optimization process. Otherwise, the region Markowchain.IfnoconfigurationisfoundinthekthMar-
will not be thoroughly searched and thus the algorithm kow chain (inner loop) that is better than the best current
may become trapped in a local minimum. In the physical configurationL0 issettothevalueofL .Ifaconfiguration
k k
analogy mentioned earlier, choosing high T corresponds which is better than all is found, the inner loop is termi-
0
to heating up the solid until all particles are randomly nated,andL0 issettothenumberofiterationsactuallyexe-
k
arranged in the liquid phase such that atoms may freely cuted in this loop. In this study, a =0.999 and
max
arrange themselves. a =0.99 were chosen as appropriate values, which pro-
min
Length of the Markow chains: Iterations during which vided a compromise between the two effectiveness criteria,
thevalueofthetemperatureparameter,T,iskeptconstant i.e. short computational time and high reliability of locat-
are called Markov chains (or inner loops). Ali et al. [90] ingglobalminimum.Initialvalueofaistakenastheaver-
adopted the following equation to decide on the length of age of a and a .
max min
a Markow chain (the number of trials (or iterations)) for Because,anewconfigurationisgeneratedbygivingper-
the kth level of T: turbations to a configuration selected randomly among a
large number of current configurations, and once accepted
L ¼LþLð1(cid:2)ef‘(cid:2)fhÞ: ð6Þ
k the new configuration replaces only the worst one, cooling
Here rate hasless importance inDSA in comparisonto conven-
tional SA. The worst configuration seldom plays a role in
L¼10n; ð7Þ
the generation process. Again, it is important to start with
where n is the dimension of the problem given by Eq. (4). a temperature high enough to accept almost all generated
At high temperatures, the current configurations form a configurations.
sparse cluster; consequently f (cid:2)f is large and Markow Stopping criterion: The optimization process is termi-
h ‘
chain length is close to 2L. On the other hand, when they nated when no further appreciable improvement can be
form a dense cluster at low temperatures, it approaches achieved. This occurs when three conditions are satisfied.
to L. During the execution of kth Markow chain with Firstly,thecurrenttemperatureparametershouldbesmall,
lengthL ,ifa configuration isgenerated having a cost less in other words it should be close to the ‘‘freezing tempera-
k
than f , the best value in A, the current chain is stopped, ture’’(zero),i.e.nofurtherworsedesignsmaybeaccepted.
‘
andanewMarkowchainisstarted.Ifaconfigurationbet- Secondly, the current configurations, A, should form a
ter than the best configuration in A is not found, the com- densecluster;thusallalternativeswereeliminatedandopti-
plete chain of length L is executed. mization process converged to a single configuration.
k
The scheme for decreasing the temperature parameter: Thirdly, the maximum moving distance of key points,
After the initial value of the temperature parameter, T , R , should be small; this means the neighborhood of
0 max
isdetermined,adecrementrulemustbeestablishedtofind the current configurations within which optimum design
the subsequent values of the temperature parameter. The issoughtbecameverysmall.Accordingly,thestoppingcri-
probability of reaching the global optimum solution terion can be expressed as
F.O.Sonmez/Comput.MethodsAppl.Mech.Engrg.196(2007)3279–3299 3285
f (cid:2)f <e and of the global minimum. The value of e depends on the
h ‘ 1 3
dimensions of the structure.
T <e and ð10Þ
k 2
R <e ;
max 3
3.7. Solution procedure
wheree aresmallnumbers.Fortheoptimizationproblems
i
that were considered, it was more than sufficient to take As a first step in the optimization procedure, the
e =e =0.0001, and e =0.05 mm for a good estimate designer defines the force and displacement boundary
1 2 3
Fig.3. Thesolutionprocedure.
3286 F.O.Sonmez/Comput.MethodsAppl.Mech.Engrg.196(2007)3279–3299
conditions, chooses the material, and decides on the num- Anothersourceoflowaccuracyisduetoerrorsincalcu-
berofkeypoints,themovingboundariesandhowtheyare lating the cost of a configuration [35,86]. In many struc-
connected. Fig. 3 describes the procedure followed in this tural optimization problems, maximum stress value is
study to solve the optimization problem. First, N number used either in calculating the cost function or in checking
of initial configurations are randomly generated within constraint violations. Designers usually carry out a finite
the search domain, S; their costs are calculated and the element (FE) analysis to calculate the stress state in the
best, worse, and worst ones are determined. An initial structure,buttheytendtochoosearoughFEmeshtoalle-
value for the temperature parameter is chosen, which viatethecomputationalburden.However,thismayleadto
should be high enough to accept nearly all randomly gen- erroneous values of stress, and the resulting design, as
erated new shapes. In the figure, k and j show the inner shown before [86], will not even be similar to the globally
loop(Markowchain)numberandthenumberofgenerated optimum design.
configurations in the current inner loop, respectively. In Imprecise definition of shape also leads to optimal
each iteration of the Markow chain, a new configuration shapes not similar to the best possible shape, which will
isgeneratedwithinthesearchdomain,S,bygivingrandom be discussed next.
displacements to the key points of a randomly chosen cur-
rent configuration. If a key point falls outside S, random 5. Precision of an optimum design
repositioning of this key point is repeated. If connectivity
is lost, the generation process is repeated until a valid How well the optimized system is defined by the design
geometryisobtained.Iftheacceptabilityofthenewconfig- variablesisameasureofprecision.Someoftheparameters
uration calculated by Eq. (5) is greater than a randomly thatdefinethesystemareallowedtobechangedduringan
generated number, it is accepted, replacing the worst con- optimizationprocess.Thenumberoftheseparametersand
figuration. If its cost is lower than the best cost, the inner the range of values they may take determine the degree at
loop is terminated; otherwise the iterations are continued which the system can be tailored to the best performance.
until the iteration number, j, becomes equal to L calcu- By increasing the number of design variables and range
k
latedbyEq.(6).Whenaninnerloopends,thetemperature oftheirvalues,onemayobtainabetterdefinitionandalso
parameter is decreased. With a lower temperature, accept- a better optimum design. In the shape optimization prob-
ability decreases as Eq. (5) implies. This means that a con- lem, by using a larger number of moving key points, one
figuration having a cost higher than the cost of the worst may better describe the boundary, i.e. shape of the struc-
configuration is less frequently accepted. ture, and minimize its volume to a largerextent. However,
As in the physical annealing process where the mobility the chances of locating the globally optimum design
of the atoms diminishes during cool down, mobility of the become lowerandlowerasthe number ofdesign variables
key points is decreased as acceptability decreases. This is increases, at the same time computational times get longer
done when a configuration is not found during a Markow andlonger.Withalargenumberofoptimizationvariables,
loop which is better than the worse current configuration. locating the global minimum becomes difficult. As the
Inthisway,theneighborhoodofthecurrentconfigurations designer tries to obtain a more precise definition, he or
within which a new configuration is generated narrows she becomes less sure of the accuracy of the results. One
downasacceptabilitybecomeslowerandlower.Ifashape mayrepeattheoptimizationprocessmanytimesanddesig-
havingacostlowerthanthecostoftheworsecurrentcon- nate the design having the lowest cost as the global opti-
figurationisnotfoundduringtheinnerloopiterations,the mum design. However, because of long computational
maximum moving distance, R , is reduced. The optimi- times this approach is not feasible. Another way to over-
max
zationprocessiscontinueduntilthreestoppingcriteriasta- come this difficulty is to start optimization by choosing
ted in Eq. (10) are satisfied. only a few design variables and finding the optimum
design.Inthiscase,thereliabilityoflocatinggloballyopti-
4. Accuracy of an optimum design mumdesignwillbehigh,eventhoughprecisionwillbelow.
Then, the designer keeps increasing the number of vari-
The measure of accuracy for an optimum configuration ables and finding the optimum designs. With the assump-
is the degree of how well it represents the best possible tion that the discrepancy between the successively
design. obtained optimum designs will not be large as the number
As one of the sources detracting from the accuracy, of design variables is increased, one may observe conver-
search algorithm may get stuck into one of the local opti- gence to the most precise globally optimum design.
mums that fail to approximate the global optimum. One
may not ensure that the resulting configuration is globally 6. Results and discussions
optimal, but may use relatively reliable search algorithms
suchassimulatedannealing.Withagoodchoiceofoptimi- The optimization algorithm was implemented using the
zationparameters,onemayachievereliabilitygreaterthan parametric design language of ANSYS. By just executing
90%atthesametimeavoidexcessivelylongcomputational its commands, key points and spline curves are generated,
times. connectivity of the model is checked, the area bounded by
F.O.Sonmez/Comput.MethodsAppl.Mech.Engrg.196(2007)3279–3299 3287
these curves is calculated, a completely new finite element
mesh is automatically generated in each iteration, and the
maximum equivalent stress and cost are calculated.
Inalltheproblemsthatwereconsidered,thematerialof
thestructurewassteelwithanelasticmodulusof207MPa
and Poisson’s ratio of 0.28 and the structure was in plane-
stresscondition.Stresseswereassumedtobewithinthelin-
early elastic region of the stress–strain curve.
6.1. Optimal design of an eccentrically loaded plate
Oneoftheproblemsthatwereconsideredinthisstudyis
the optimal design of an eccentrically loaded plate
restrained at one end and loaded at the other (Fig. 4).
One segment of the border on the right is fixed in length
andsubjectto200MPapressure.Theleftborderisdefined
by a line between two key points; the lower one is fixed,
while the upper one is free but restrained to move in the
vertical direction. The search domain, S, within which
Fig.5. Theoptimalshapeoftheeccentricplatedefinedbyfivekeypoints.
thekeypointscanmove,isdefinedbyarectangleofdimen-
sion 10cm·12cm. The allowable stress of the material is
taken as 300MPa.
First,fivekeypointswereusedtodescribetheboundary
and using their coordinates as design variables the best
shape was found. Fig. 5 shows the optimal shape with its
finite element mesh. The cross-section becomes larger and
largerawayfromtheloadedregionsduetoincreasedbend-
ing moment. Because the number of design variables was
low, the optimum design could be obtained with high reli-
ability. Then, optimization process was repeated using
sevenandninemovingkeypointsresultingintheoptimum
designsshowninFigs.6and7.Themaximummovingdis-
tance, R , was initially 1.00cm, and at the end of the
max
optimization process, it decreased to a value of
0.036mm. In this way, the algorithm could search a large
region initially, and a close neighborhood towards the
end of the optimization. The lateral areas of the optimum
shapes defined by five, seven, and nine key points turned
out to be 37.100, 36.313, and 35.903cm2, respectively.
Consequently, with a higher number of key points, we
Fig. 6. The optimal shape of the eccentric plate defined by seven key
points.
could obtain a more precise definition of shape and also
Moving find an optimum configuration with a lower cost. Even
Fixed though reliability of the solution becomes low with a large
number of design variables, because it approximates the
more reliable optimum designs defined by smaller number
of key points, the reliability is ensured.
In order to ensure the accuracy at which the cost was
5 cm calculated, convergence of the FE solution was frequently
checked during the optimization process. If the change in
the magnitude of the maximum stress was greater than
0.5 cm
0.5% when FE analysis of currently best configuration
10 cm
wascarriedoutusingonefourthofthecurrentsizeoffinite
Fig.4. Optimumdesignproblemforaneccentricallyloadedplate. elements,thentheFEmeshwasrefined.Thesignificanceof
3288 F.O.Sonmez/Comput.MethodsAppl.Mech.Engrg.196(2007)3279–3299
Moving
10 cm
Fixed
Fig.7. Theoptimalshapeoftheeccentricplatedefinedbyninekeypoints.
accuracy in FE calculations is better seen in Fig. 8, which
shows the optimum shape obtained when a rough mesh is Fig.9. Optimumdesignproblemforahook.
used. Although the error in stress level is only 5%, the dis-
crepancy in shape was quite large.
move, but only in the horizontal direction and symmetri-
cally such that their midpoint retains its position. The rest
6.2. Optimal design of a hook of the key points on the border are free to move within a
search domain, S, of dimension 13cm·18cm. The allow-
Anotherproblemconsideredinthisstudyistheoptimal able stress of the material is taken as 800MPa.
designofaplaterestrainedatthetopedge,andloadedover At the start of the optimization, initial current configu-
an arc of surface as shown in Fig. 9, which corresponds to rationsarecreatedbydescribingsplinecurvesthroughkey
loading of a hook. The half circle has a radius of 1.5cm. points randomly generated within the search domain, S.
This curved border is fixed in length and subject to They are then ordered with respect to their cost; the ones
150MPa pressure at the center decreasing to zero magni- having the lowest and highest costs are named as the best
tudetowardthesides.The restrained topborder isdefined andworst configurations, respectively.Fig. 10showssome
by a line between two key points; both of them are free to these initial configurations created for the optimum hook
design problem including the best (Fig. 10a) and worst
(Fig. 10b) ones.
Fig.11showstheoptimalshapesobtainedusing6,8,10,
12,and13movingkeypoints.Thelateralareasofthemare
37.888,34.343,33.959,33.077,and33.412cm2,respectively.
When the number of key points is too large, locating the
global minimum becomes extremely difficult. More or less
the same curve can be generated by many different sets of
keypoints.Then,therecouldbenumerousnearglobalopti-
mumdesigns.Forthisreason,thealgorithmmaygetstuck
intooneofthemasinthiscase.Theoptimumshapedefined
by 13 key points is worse than the one defined by 12 key
points.Inthesecases,onemaystartatahighertemperature
and decrease the reduction ratio of the maximum moving
distancetoincreasetheprobabilityofconvergencetoabet-
terdesignattheexpenseoflongercomputationaltimes;but
there isa more effectivemethod to be discussedlater.
The maximum moving distance, R , was initially
max
1.60cm, and towards the end of the optimization process,
it decreased to a value of 0.092mm. The upper regions of
Fig.8. Inaccurateoptimalshapeobtainedwitharoughmesh. thehookaresubjecttopureaxialload;forthisreasonthey