On the Design of an Adaptive Simulated Annealing Algorithm

15 Pages·2007·0.44 MB·English

by  Vincent A. Cicirello

On the Design of an Adaptive Simulated Annealing
Algorithm
VincentA.Cicirello
ComputerScienceandInformationSystems
TheRichardStocktonCollegeofNewJersey
Pomona,NJ08240
[email protected]
Abstract. Inthispaper,wedemonstratetheeaseinwhichanadaptivesimulated
annealingalgorithmcanbedesigned.Speci(cid:2)cally,weusetheadaptiveannealing
scheduleknownasthemodi(cid:2)edLamscheduletoapplysimulatedannealingto
theweightedtardinessschedulingproblemwithsequence-dependentsetups.The
modi(cid:2)edLamannealingscheduleadjuststhetemperaturetotrackthetheoretical
optimal rate of accepted moves. Employing the modi(cid:2)ed Lam schedule allows
ustoavoidtheoftentedioustuningoftheannealingschedule;asthealgorithm
tunes itselffor each instance during problem solving. Our resultsshow that an
adaptivesimulatedannealercanbecompetitivewhencomparedtohighlytuned,
hand crafted algorithms. Speci(cid:2)cally, we compare our results to a state-of-the-
artgeneticalgorithmforweightedtardinessschedulingwithsequence-dependent
setups.Ourstudyservesasanillustrationoftheeasewithwhichaparameter-free
simulatedannealercanbedesignedandimplemented.
1 Introduction
One of the often touted advantages of metaheuristic search algorithms such as simu-
latedannealing[1,2],tabusearch[3,4],geneticalgorithms(GA)[5,6],andthelike,is
thattheyareoftenarguedtoberelativelysimpletoimplementfornewproblems.While
insomesensethisisgenerallytrue,inanothersensethis(cid:147)easy toimplement(cid:148) advan-
tagelargelyignoresthenextstepoftuningparameterstoactuallygetthealgorithmto
performwellforyourproblem.
Mostmetaheuristicsarecontrolledbyseveralparameters.Forexample,thesimple
geneticalgorithm(SGA)hasapopulationsize,acrossoverrate,amutationrate.GAs
thatareslightlyfancierthantheSGAcanhaveadditionalparameterssuchasgeneration
gap,scalingwindow,elitismrate,etc.1 Likewise,algorithmssuchassimulatedanneal-
ing(SA)andtabusearch(TS)canalsohaveseveralcontrolparameters.Forexample,
arunofSAislargelycontrolledbythevalueofsomethingcalledthetemperature.The
temperatureinSAsearchgenerallybeginsatahighvalueanddecaysduringexecution
of the search. This temperature decay is controlled by something usually referred to
as the annealing schedule. All of the most commonly used annealing schedules have
parameters.Theseparameters,liketheparametersofaGA,mustbesetinsomeway.
1Forthesakeofthispaper,youdonotneedtoknowwhatanyoftheseGAparametersare.The
pointisthattherecanbeseveralofthem.
PerhapsthemostcommonapproachtotuningtheparametersofaGA,SA,TS,or
someothermetaheuristicis,unfortunately,atedioushand-tuning.Thatis,thealgorithm
implementersystematically(ornotsosystematically)triesoutvarioussetsofparame-
tervaluesonsomesetofprobleminstances,andselectsthesetofparametervaluesthat
appearstoperformwell.Fortunately,thereareotherswhohavetakenamorerigorous
approachtotheproblem.Forexample,theearliestformalexaminationoftheproblem
of parameter tuning in a GA was that of De Jong [7], who studied the performance
of the GA on a class of function optimization problems and empirically determined
an(cid:147)optimal(cid:148) setof parametervalues forthatclassofproblem.For decadessincethat
study,DeJong’sparametersetisoftenusedbyothersindesigningGAsfornewprob-
lems.However,DeJong’sparametersetisreallyonlyrelevantfortheclassofproblems
considered in his study. Others later followed De Jong’s lead and began looking for
waystoautomatetheparametertuningprocess.Grefenstette,forexample,introduced
theideaofusingaGAtooptimizetheparametersofaGA[8].Sincethen,ahostofre-
latedmeta-levelcontrolparameteroptimizationapproacheshavebeenintroduced,both
for tuning GA parameters [9(cid:150)12] as well as for parameters of other problem solving
systems [13(cid:150)16]. Others view the parameter tuning problem as one that should coin-
cidewithproblemsolving.Thatis,ratherthantuningcontrolparametersaheadoftime,
somearguethatproblemsolvingfeedbackcanbeusedtoadaptcontrolparametersto
theprobleminstanceathand(e.g.,[17(cid:150)19]).
In this paper, we too take the approach of adapting control parameters during the
search.Speci(cid:2)cally ,wefocusonsimulatedannealinganduseanexistingadaptivean-
nealing schedule to demonstrate the ease with which a parameter-free simulated an-
nealercanbeimplemented.OurempiricalresultsontheNP-Hardoptimizationproblem
knownasweightedtardinessschedulingwithsequence-dependentsetupsillustratesthat
suchanapproachcanbecompetitivewithevenhighly-tunedhand-craftedalgorithms.
Theremainderofthepaperproceedsasfollows.InSection2weoverviewsimulated
annealing,andspeci(cid:2)cally describetheadaptiveannealingscheduleknownasModi(cid:2)ed
Lam[19].NextinSection3weprovideanoverviewoftheweightedtardinessschedul-
ingproblemwithsequence-dependentsetupsthatweuseforourexperimentalstudy.In
Section4weprovidefurtherdetailsonthedesignofoursimulatedannealer.Wethen
presentempiricalresultsinSection5andconcludeinSection6.
2 TheModi(cid:2)edLamAnnealingSchedule
Figure1providespseudocodeforthetypicalSimulatedAnnealingalgorithm.Initially,
SAbeginswithsomearbitrarysolutiontotheproblem.Thisisoftenrandomlygener-
ated,althoughsometimesitcanbetheresultofsomeotherprocedure.SAthenconsid-
ersaseriesofrandommovesfromthecurrentstateS usinganeighborhoodfunction.
Forexample,iftheproblemSAissolvingisanorderingproblemsuchasthetraveling
salesperson problem, then the neighborhood function might be pairwise interchanges
of randomly selected cities from the current ordering. If the state that results is of a
lower cost than the current state, the move is accepted and this new state S0 becomes
thecurrentstateS.Otherwise,ifthenewstateisofhighercost(lesserqualitysolution),
SAacceptsitwithprobability,e(Cost(S)(cid:0)Cost(S0))=T.ThisisknownastheBoltzmann
SimulatedAnnealing
S GenerateInitialState
T SomeInitialHighTemperature,T0
forifrom1toEvalsmax
S0 PickRandomState(Neighborhood(S))
ifCost(S0)<Cost(S)
S S0fNote:acceptingamoveg
else
r Random(0;1)
ifr<e(Cost(S)(cid:0)Cost(S0))=T
S S0fNote:acceptingamoveg
end
end
T T0(cid:1)(cid:11)bRic
end
Fig.1.PseudocodeforastandardSimulatedAnnealingalgorithmusingthecommongeometric
annealingschedule.
distribution.Thatis,SAaccepts(cid:147)bad(cid:148) moveswithaprobabilitythatdependsonhow
(cid:147)bad(cid:148) thatmoveisandwhichalsodependsonthecurrentvalueofthetemperatureT.
The problem-solving effectiveness of SA largely depends on the design of what
is referred to as the annealing schedule. The annealing schedule refers to the way in
whichthetemperatureparameterT isupdatedduringthesearch.Themostcommonly
usedannealingscheduleisthegeometricschedule,whichisalsotheannealingschedule
suggestedbytheoriginatorsofSA[1].Thegeometricannealingscheduleisde(cid:2)ned as:
T T0 (cid:11)bRic; (1)
(cid:1)
where T is an initially (cid:147)high(cid:148) temperature, (cid:11) is referred to as the cooling rate, i is
0
simplytheindex ofthecurrentevaluationandcorrespondstotheiinthepseudocode
ofFigure1,andRistheroundlengthwhichcontrolshowmanyevaluationsaremade
foreachtemperatureT.ThevaluesofT ,(cid:11),andRareallparametersthatneedtobe
0
tunedduringthedesignoftheSA.
A well known theoretical result indicates that if the annealing schedule (cid:147)cools(cid:148)
thetemperatureatasuf(cid:2)ciently slowratetomaintainthesysteminastateofthermal
equilibrium,thenwithprobability1simulatedannealingwill(cid:2)nd thegloballyoptimal
solution. The problem with this result is that it requires such a slow cooling rate that
searchesrequireaprohibitivelylongannealingruntosettleuponasolution.
LamandDelosmeproposedanapproximatethermalequilibriumtheycallD-equi-
librium which balances the trade-off of required computation time and the quality of
thesolutionfoundbytherunofSA[20].Undercertainassumptionsabouttheformsof
thedistributionofthecostvaluesandthedistributionofcostvaluechanges,theyana-
lyzedtheirmodelanddeterminedtheannealingschedulethatmaintainsthesystemin
SimulatedAnnealingwithModi(cid:2)edLamAnnealingSchedule
S GenerateInitialState
T 0:5
AcceptRate 0:5
forifrom1toEvalsmax
S0 PickRandomState(Neighborhood(S))
ifCost(S0)<Cost(S)
S S0fNote:acceptingamoveg
AcceptRate 1 (499(cid:1)AcceptRate+1)
500
else
r Random(0;1)
ifr<e(Cost(S)(cid:0)Cost(S0))=T
S S0fNote:acceptingamoveg
AcceptRate 1 (499(cid:1)AcceptRate+1)
500
else
fNote:rejectingamoveg
AcceptRate 1 (499(cid:1)AcceptRate)
500
end
end
ifi=Evalsmax <0:15thenLamRate 0:44+0:56(cid:1)560(cid:0)i=Evalsmax=0:15
if0:15(cid:20)i=Evalsmax <0:65thenLamRate 0:44
if0:65(cid:20)i=EvalsmaxthenLamRate 0:44(cid:1)440(cid:0)(i=Evalsmax(cid:0)0:65)=0:35
ifAcceptRate>LamRate
T 0:999(cid:1)T
else
T T=0:999
end
end
Fig.2. Pseudocode for Simulated Annealing using Boyan’s version of Swartz’s Modi(cid:2)ed Lam
annealingschedule.
D-equilibrium (i.e., the annealing schedule that optimally balances the computational
cost / solution quality trade-off). This (cid:147)optimal(cid:148) annealing schedule adjusts the tem-
perature based on the parameter (cid:21) which controls the cost-quality trade-off and more
importantly based on the current rate of accepted moves. Analyzing their annealing
schedule, Lam and Delosme determined that the temperature is reduced the quickest
when the probability of accepting a move is equal to 0.44. Since a faster cooling rate
leads to a shorter annealing run, Lam and Delosme proposed allowing the size of the
neighborhoodconsideredformovestodynamically(cid:3)uctuate upordowntomatchthis
target move acceptance rate of 0.44 as closely as possible [20]. The general idea is
that you can increase the acceptance rate by decreasing the maximum distance from
yourcurrentstatethatthealgorithmconsidersfornextstates;andlikewisedecreasethe
acceptance rate by increasing this distance. The motivation for this is the assumption
that nearby search states are of similar quality. If you make shorter distance moves,
thedifferencein(cid:2)tness ofyourcurrentandpotentialfuturestatesshouldbesmalleron
averagethanifyouallowlargermoves.LamandDelosmevalidatedtheirmodelonthe
TravelingSalespersonProblemaswellasonacellplacementproblem,buttheanneal-
ing schedule itself is the result of a problem independent theoretical analysis of their
approximatethermalequilibriummodel.
Swartz later extensively studied Lam and Delosme’s annealing schedule and pre-
sentedamodi(cid:2)cation [21].First,heperformedadditionalempiricalanalysisoftheLam
scheduleanddeterminedthatatthebeginningofthesearch,themoveacceptancerate
isnear1.0(i.e.,allmovesatthebeginningareaccepted).Theacceptanceraterapidly
decreases (at an exponential rate) during the (cid:2)rst 15% of the run until it reaches the
targetof0.44.Itremainsnearlyconstantforthenext50%oftherun,andthenbeginsan
exponentialdeclinetoanacceptancerateof0%bytheendoftherun.Swartzthennoted
thatadjustingthesizeofthemoveneighborhoodisnottheonlywaytotrackthisthe-
oreticaltargetacceptancerate.Speci(cid:2)cally ,hepresentedanapproachfortrackingthis
targetacceptanceratebyadaptingthevalueofthetemperatureparameterupordown;
whereasLamandDelosmeusedamonotonicallydecreasingtemperaturescheduleand
adjusted the neighborhood function. Boyan presented a similar approach to adapting
thetemperatureparameter[19]totrackthetheoretical(cid:147)optimal(cid:148) acceptanceratewith-
out modifying the neighborhood function during the search. Swartz’s (and Boyan’s)
Modi(cid:2)ed Lam annealing schedule has the advantage over the original Lam schedule
in that it is problem independent; whereas adapting the neighborhood function would
requireaproblemdependentanalysis.Intheexperimentsofthispaper,weuseBoyan’s
recommendedimplementation[19]ofSwartz’sModi(cid:2)ed Lamannealingschedule[21].
PseudocodeforSAusingtheModi(cid:2)ed LamannealingscheduleisshowninFigure2.It
requiresnoparametertuningbytheimplementer.Thealgorithmusesproblem-solving
feedbacktoadjustthetemperatureT appropriately.
3 Weighted Tardiness Scheduling with Sequence-Dependent
Setups
The weighted tardiness scheduling problem with sequence-dependent setups can be
de(cid:2)ned asfollows.TheprobleminstanceconsistsofasetofjobsJ = j ;j ;:::;j .
0 1 N
f g
Eachofthejobsjhasaweightw ,duedated ,andprocesstimep .Furthermore,s is
j j j i;j
de(cid:2)ned astheamountofsetuptimerequiredimmediatelypriortothestartofprocessing
jobjifitfollowsjobionthemachine.Itisnotnecessarilythecasethats =s .The
i;j j;i
0-th(cid:147)job(cid:148) isthestartoftheproblem(p =0,d =0,s =0,w =0).Itspurposeis
0 0 i;0 0
forthespeci(cid:2)cation ofthesetuptimeofeachofthejobsifsequenced(cid:2)rst.
TheweightedtardinessobjectiveistosequencethesetofjobsJ tominimize:
T = w T = w max(c d ;0); (2)
j j j j j
(cid:0)
Xj2J Xj2J
whereT isthetardinessofjobj;andc ,d arethecompletiontimeandduedateofjob
j j j
j.Thecompletiontimeisequaltothesumoftheprocesstimesandsetuptimesofall
jobsthatcomebeforej inthesequenceplusthesetuptimeandprocesstimeofthejob
j.Speci(cid:2)cally ,let(cid:25)(j)bethepositioninthesequenceofjobj.De(cid:2)ne c as:
j
c = p +s : (3)
j i k;i
X
i;k2J;(cid:25)(i)<=(cid:25)(j);(cid:25)(i)=(cid:25)(k)+1
Single-machine scheduling to optimize the weighted tardiness objective function
is NP-Hard even if setups are independent of job ordering [22]. In the case that we
consider here, the challenge of the problem is greatly magni(cid:2)ed by the sequence-
dependentsetup constraints.Sen and Bagchi discuss the signi(cid:2)cance of the challenge
that sequence-dependent setups pose for exact solution procedures [23]. Speci(cid:2)cally ,
theydiscusshowsequence-dependentsetupsinduceanon-order-preservingpropertyof
theevaluationfunction.Atthetimeoftheirwriting,exactsolutionproceduressuchas
A*,Branch-and-Boundalgorithms,andtheirownGREC[23]forsequencingproblems
with sequence-dependent setups were limited to solving instances with no more than
approximately 25-30 jobs, even for easier objective functions. Problem instances of
largersizerequireturningtoinexactsolutionproceduressuchasmetaheuristicsearch.
4 DesigningAnAdaptiveSimulatedAnnealingAlgorithm
InthisSection,wedetailourdesignforasimulatedannealerfortheweightedtardiness
schedulingproblemwithsequence-dependentsetups.WeusetheModi(cid:2)ed Lamanneal-
ingscheduleasdiscussedearliertodesignaSAalgorithmfreeofanyparametersthat
requiretuning.Therearesimplytwodesigndecisionstomake:(1)whatneighborhood
functiontouse;and(2)howtogeneratetheinitialsolution.Weconsidereachofthese
inturnbelow.
Choosing the Neighborhood Function: The weighted tardiness scheduling problem
with sequence-dependent setups is a one machine sequencing problem. Any ordering
ofthejobsisavalidsolution.Theproblemisto(cid:2)nd agoodordering(i.e.,onewitha
lowvalueoftheweightedtardinessobjective).ThroughoutourSAsearch,thecurrent
search state is a complete ordering over the set of jobs. We consider three alternative
problem-independentneighborhoodoperators:
(cid:150) Insert: The insert operator randomly selects a job from the current sequence, re-
movesthatjob,andreinsertsitintoarandomlyselectednewposition.
(cid:150) Swap:Theswapoperatorrandomlyselectstwojobsfromthecurrentsequenceand
swapstheirpositions.
(cid:150) Scramble:Thescrambleoperatorrandomlyselectstwopositionsinthesequence
andrandomlyshuf(cid:3)es thejobslocatedbetweenthosetwopositions.
Thesethreealternativeoperatorsareapplicabletoanyorderingproblemandarecom-
mon choices of neighborhood operator. Later, we discuss empirical results using all
three. More sophisticated problem-dependent options are also possible here, but our
goalistoillustratethepowerofasimpleproblem-independentadaptiveSA.
Fig.3.Comparisonofneighborhoodoperators(scramble,swap,andinsert).Showsaverageper-
centdeviationfrombestknownsolutionstothebenchmarksetfordifferentlengthrunsofthe
adaptiveSAusingaheuristicallydeterminedinitialsolution.99%con(cid:2)denceintervalsareshown.
RandomStartvsHeuristicStart: Thenextdecisionwemustmakeishowtogenerate
theinitialcon(cid:2)guration ofthesearch.IfwewantedourSAapproachtobecompletely
problem-independent(asidefromtheoptimizationobjective),thentheSAsearchcould
simply begin with a random permutation of the jobs. We consider this option in our
experimentalresultsandrefertoitasRandom.
Wealso,however,considertheuseofaproblemspeci(cid:2)c heuristictoseedtheinitial
state of the search. Speci(cid:2)cally , we consider the use of the Apparent Tardiness Cost
with Setups (ATCS) heuristic [24]. ATCS was designed for use as a dispatch policy.
That is, its intended use is as a means of selecting the next job for the sequence in
a dynamic setting. It has also been used as a search heuristic for stochastic sampling
algorithms[25,26].ATCSisde(cid:2)ned asfollows:
w max(d p t;0) s
ATCS (t;l)= j exp( j (cid:0) j (cid:0) l;j); (4)
j p (cid:0) k p(cid:22) (cid:0) k s(cid:22)
j 1 2
wheretisthecurrenttime(orthesumoftheprocessandsetuptimesofthejobsalready
sequenced);listheindexofthejobmostrecentlyaddedtotheschedule;p(cid:22)istheaverage
processingtimeofalljobs;ands(cid:22)istheaveragesetuptime.
Fig.4.Comparisonofneighborhoodoperators(scramble,swap,andinsert).Showsaverageper-
centdeviationfrombestknownsolutionstothebenchmarksetfordifferentlengthrunsofthe
adaptiveSAusingarandominitialsolution.99%con(cid:2)denceintervalsareshown.
Thek andk areparametersthatareautomaticallytunedbasedoncharacteristics
1 2
oftheprobleminstanceasfollows.
4:5+R ifR 0:5
k1 =(cid:26)6:0 2Rother(cid:20)wise ; (5)
(cid:0)
where R is called the duedate range factor which is an indicator of how spread out
the duedates of the jobs are, and is de(cid:2)ned as: R = dmax(cid:0)dmin where d , d are
Cmax max min
the maximum and minimum duedates, and C is an estimate of the makespan (or
max
completiontimeofthelastjob).
(cid:28)
k = : (6)
2 2p(cid:17)
where(cid:28) iscalledtheduedatetightnessfactor,de(cid:2)ned by(cid:28) = 1 d(cid:22) ;and(cid:17) isthe
(cid:0) Cmax
setup time severity factor, de(cid:2)ned by (cid:17) = s(cid:22). The (cid:28) is an indicator of how urgent the
p(cid:22)
duedatesoftheproblemareand(cid:17)isanindicatorofhowlargeaneffectthesetuptimes
canhaveonproblemsolving.
In our experiments, we consider the ATCS heuristic as an alternative to random
initial search states. In this option, ATCS generates a single solution which becomes
theseedfortherunofSA.WerefertothisoptionintheexperimentsasHeuristic.
Fig.5.Comparisonofneighborhoodoperators(swapandinsert).Showsaveragepercentdevia-
tionfrombestknownsolutionstothebenchmarksetfordifferentlengthrunsoftheadaptivesim-
ulatedannealingalgorithmusingarandominitialsolution.99%con(cid:2)denceintervalsareshown.
5 ExperimentalResults
In this Section, we present experimental results. We begin with comparisons of the
threealternativeneighborhoodoperatorsandthetwoalternativesfortheinitialsearch
state. We use an existing set of benchmark instances [26]. This set consists of 120
probleminstances,eachwith60jobs,withvaryinglevelsofduedatetightness,duedate
range,andsetuptimeseverity.Nocurrentlyavailableapproachhasbeenabletoprovide
provably optimal solutions to the problem instances in this benchmark set. There are,
however,anumberofmetaheuristicsdesignedfortheproblem.
Figure3showsacomparisonofthethreeneighborhoodoperatorsfortheSAusing
theATCSheuristictoseedtheinitialstateforthesearch.Fivedifferentlengthsearches
areshownwithatotalnumberofevaluationsof1000,104,105,106,and107.Thetotal
numberofevaluationsisontheXaxisusinglogscale.TheYaxisshowstheaverage
percent deviation from the best known solutions to the benchmark instances. Speci(cid:2)-
cally,foreachlengthrunandforeachchoiceofoperator,wereportaveragesof1200
runs (10 runs on each of the 120 benchmark instances). 99% con(cid:2)dence intervals are
shown.Whenusingtheheuristictoseedthestartingcon(cid:2)guration, theclearlydominant
choiceofneighborhoodoperatoristheInsertoperator.
Figure4similarlyshowsacomparisonofthethreeneighborhoodoperators,butthis
timeusingrandompermutationstoseedtheSAsearch.Again,theresultsareaverages
Fig.6. Comparison of the adaptive simulated annealing algorithm using a heuristically deter-
minedinitialsolutionandusingarandominitialsolution.Bothversionsofthealgorithmusethe
insertneighborhoodoperator.Showsaveragepercentdeviationfrombestknownsolutionstothe
benchmarksetfordifferentlengthruns.99%con(cid:2)denceintervalsareshown.
of 1200 runs and 99% con(cid:2)dence intervals are shown. In this case, due to the partic-
ularly slow convergence when using the Scramble operator, it is dif(cid:2)cult to pick out
a clear winner. We can, however, see that Scramble is clearly dominated by either of
theothertwoneighborhoodoperators.Note,however,thatforespeciallylongsearches
(e.g.,107 totalnumberofevaluations),theperformanceoftheSAusingScrambleap-
pearstobecatchinguptotheothertwo.
Tohelpdistinguishtheperformanceoftheneighborhoodoperatorswhenusingran-
domstartingcon(cid:2)gurations, Figure5illustratestheresultsofjusttheSwapandInsert
operators,forfourdifferentlengthruns(totalnumberofevaluationsof104,105,106,
and107).Justlikeinthecaseofheuristicallydeterminedinitialstates,weseethatthe
Insertneighborhoodoperatoristheclearlydominatingperformerwhenusingrandom
startingcon(cid:2)gurations.
InFigure6,weconsidertheuseoftheATCSheuristictoseedthesearchvsusing
random initial con(cid:2)gurations. Both options in this (cid:2)gure use the Insert neighborhood
operator. Again, the results are averages of 1200 runs and 99% con(cid:2)dence intervals
are shown. For really short searches (e.g., 104 total evaluations), using the heuristic
to seed the search has a clear advantage. For slightly longer searches (e.g., 105 total
evaluations), the heuristically seeded SA still has an advantage, though not as strong.
Forlongersearchesoflength106evaluations,thereisnegligibledifference;andby107