Surrogate-enhanced evolutionary annealing simplex algorithm for effective and efficient

21 Pages·2015·4.21 MB·English

by  Ioannis Tsoukalas

EnvironmentalModelling&Software77(2016)122e142
ContentslistsavailableatScienceDirect
Environmental Modelling & Software
journal homepage: www.elsevier.com/locate/envsoft
Surrogate-enhanced evolutionary annealing simplex algorithm for
effective and efficient optimization of water resources problems on a
budget
*
Ioannis Tsoukalas , Panagiotis Kossieris, Andreas Efstratiadis, Christos Makropoulos
DepartmentofWaterResourcesandEnvironmentalEngineering,NationalTechnicalUniversityofAthens,HeroonPolytechneiou5,GR-15780,Zographou,
Greece
a r t i c l e i n f o a b s t r a c t
Articlehistory: Inwaterresourcesoptimizationproblems,theobjectivefunctionusuallypresumestofirstrunasimu-
Received10April2015 lationmodelandthenevaluateitsoutputs.However,longsimulationtimesmayposesignificantbarriers
Receivedinrevisedform to the procedure. Often, to obtain a solution within a reasonable time, the user has to substantially
1December2015
restrict the allowable numberof function evaluations, thus terminating the search much earlier than
Accepted7December2015
required. A promising strategy to address these shortcomings is the use of surrogate modeling tech-
Availableonlinexxx
niques.HereweintroducetheSurrogate-EnhancedEvolutionaryAnnealing-Simplex(SEEAS)algorithm
thatcouplesthestrengthsofsurrogatemodelingwiththeeffectivenessandefficiencyoftheevolutionary
Keywords:
annealing-simplex method. SEEAS combines three different optimization approaches (evolutionary
Globaloptimization
Meta-models search, simulated annealing, downhill simplex). Its performance is benchmarked against other
Radialbasisfunctions surrogate-assisted algorithms in several test functions and two water resources applications (model
Hydrologicalcalibration calibration,reservoirmanagement).ResultsrevealthesignificantpotentialofusingSEEASinchallenging
Multi-reservoirmanagement optimizationproblemsonabudget.
Syntheticdata ©2015ElsevierLtd.Allrightsreserved.
Softwareavailability space to maximize the system performance; at each trial, new
valuesareassignedtothecontrolvariablesofthesimulationmodel,
EASandSEEAS(RandMATLABversions)areavailableonlineat:: which runs automatically to update the value of the objective
http://www.itia.ntua.gr/en/softinfo/29/ function.
Combinedsimulation-optimizationschemesforwaterresource
systems can be generally classified into two categories: (a)
1. Introduction
decision-making problems, in which the system properties and
associated processes are known a priori, but either some of its
Couplingofsimulationandoptimizationisapowerfultechnique
thathasgainedsignificantattentioninwaterresourcesscienceand designquantitiesoritsmanagementpolicyareunknown;and(b)
calibration problems, in which some internal properties of the
technology, sinceitensures greatadvantages overthetraditional
system,eitherphysicalorconceptual,areunknownandhavetobe
individual implementation of the two approaches (e.g.,
invertedbyminimizingthedeparturesofthesimulatedresponses
KoutsoyiannisandEconomou,2003).Inthiscontext,asimulation
against the observed ones. Despite their different rationale, both
model is used to faithfully represent the dynamics of the system
types of problems suffer from significant uncertainties and com-
under study in subsequent time steps and next to evaluate its
overall performance against one or more user-specified criteria. plexities, and theyare subject to multiple (and often conflicting)
criteriaaswellasnumerousconstraints.
Provided that these criteria are expressed in terms of objective
Forconvenience,weconsiderthatallcriteriaareaggregatedina
function,simulationcanbedrivenbyanoptimizationmodel,which
singleobjectivefunctionrepresentingaglobalperformancemea-
employs systematic search through the parameter (or decision)
sure of the system (an alternative approach would require the
formulation of a multiobjective function and the identification of
* Correspondingauthor.Tel.:þ302107722842. acceptable tradeoffs among conflicting criteria, which is not the
E-mailaddress:[email protected](I.Tsoukalas). case here). We also assume that all “internal” constraints (i.e.,
http://dx.doi.org/10.1016/j.envsoft.2015.12.008
1364-8152/©2015ElsevierLtd.Allrightsreserved.
I.Tsoukalasetal./EnvironmentalModelling&Software77(2016)122e142 123
constraints associated with the system dynamics) are handled 2006; Feyen et al., 2007; He et al., 2007; Regis and Shoemaker,
through the simulation model (Koutsoyiannis and Economou, 2009; Dias et al., 2013); (2) computationally efficient optimiza-
2003),whileanyadditional“external”constraints,whichareusu- tionalgorithms(e.g.,TolsonandShoemaker,2007;Kuzminetal.,
allyassociated with decision-making problems, are embedded in 2008;Tanetal.,2008;Tolsonetal.,2009);(3)strategiestoavoid
the objective function, typically as penalty terms. Under this opportunistically(expensive)modelevaluations(e.g.,Ostfeld and
premise, the combined simulation-optimization problem is Salomons, 2005; Razavi et al., 2010; Matott et al., 2012); and (4)
formalized as the determination of the global optimum (for con- surrogatemodelingtechniques,alsoreferredtoasmeta-modeling,
venience,minimum)ofanonlinearobjectivefunctionf(x),where function approximation, response surface modeling and model
f(,)representsthesimulationmodelandxisthevectorofcontrol emulation (Razavi et al., 2012a), where surrogate approaches are
variables. The search space is a hypervolume, since the unique used to approximate the responses of the original simulation
constraintsoftheproblemarethelowerandupperboundsofpa- model.Parallelcomputing,allowingtheexecutionofindependent
rameters.Asf(x)isablack-boxfunction,itsanalyticalexpressionas simulations by multiple processors, requires significant in-
wellasitsderivativesarenotavailable,whichprohibitstheuseof vestmentsinhardwareinfrastructure,whichmakesitimpractical
gradient-basedoptimization.Givenalsothat,duetouncertainties forcommonuse.Weremarkthatinordertoreducetheentiretime
and complexities of the system, f(x) is non-convex, and thus of computations three orders of magnitude e a reasonable
multimodal(i.e.,itcontainsmultiplelocaloptima),derivative-free requirement when dealing with complex simulation models e
methods combined with stochastic search approaches are essen- 1,000parallelprocessorsshouldbeused,whichisfarfromrealistic.
tialtosolvethisso-calledglobaloptimizationproblem. The other two options, i.e., the improvement of efficiency of
The need for advanced global optimization tools (e.g., evolu- existing algorithms, as well as the interruption of the function
tionaryalgorithms)hasbeenearlyrecognizedbythehydrological evaluationprocedure, when the model performance seems to be
community,whichhassignificantexperienceintheiruseandalso verypoorfromearlystepsofsimulation,maysavesometimebut
remarkablecontributionintheirdevelopment.Intheliteratureare notasmuchasrequired.Ontheotherhand,surrogatemodelsdo
found numerous reviews of optimization approaches in such nothaveanyspecificrequirementsincomputerresourcesandalso
problems.Forinstance,inthecontextofwaterresourcesplanning ensureveryfastcomputations,sincetheyreplace,tosomecontext,
and management, we distinguish the works by Labadie (2004), the(expensive)simulationmodel.Theirkeyobjectiveistogenerate
Fowler et al. (2008), Nicklow et al. (2010), Reed et al. (2013) modelsthatareaccurateinacertainregionofthesearchspace(i.e.,
(emphasis to multiobjective applications) and Ahmad et al. around a potential optimum) and thus intelligently guide the
(2014). The literature for hydrological calibration is even more optimization(Couckuytetal.,2013).
extended.Forconvenience,wehighlighttherecentworksbyDuan Although response surface approaches go back to 70's
(2013) and Efstratiadis and Koutsoyiannis (2010), who provide a (Blanning,1975),surrogate-basedoptimizationmethodshavebeen
comprehensive review of global and multiobjective calibration popularizedsincethepioneeringworkbyJonesetal.(1998),who
approaches,respectively.Itisalsoworthmentioningthearticleby developedtheEfficientGlobalOptimization(EGO)algorithm.EGO
Maier et al. (2014), who summarize the current status of evolu- uses Kriging as surrogate model and an acquisition function
tionary algorithms and other metaheuristics, and highlight new (namedExpectedImprovement),inordertolocatepotentialgood
directionsforfutureresearchacrosswaterresourcesapplications. samples that should be evaluated through expensive simulation
Apparently,inthewholecomputationalprocedure,simulation functions(Sacksetal.,1989;Jonesetal.,1998).Later,Sasenaetal.
isbyfarthemosttime-consumingcomponent.Asmodelsbecome (2002) implemented and investigated various acquisition func-
more complex and data-demanding, their requirements in tions for EGO. Literature also reports multi-objective versions of
computationaltimeand/orCPUincreasesubstantially(e.g.,Tolson EGO(e.g.,Knowles,2005;Ponweiseretal.,2008;Couckuytetal.,
and Shoemaker, 2007; Keating et al., 2010; Razavi et al., 2010; 2013).
Tsoukalas and Makropoulos, 2015a). Typical example is the case OthercommonlyusedsurrogatemodelsareRadialBasisFunc-
of physically-based hydrological models of fine spatial and tem- tions(RBFsePowell,1992;Buhmann,2003),polynomials(Myers
poral resolution, in contrast to lumped conceptual rainfall-runoff and Montgomery, 1995), artificial neural networks, and support
models.Inotherapplications,referredtoasstochasticsimulation vectormachines(CortesandVapnik,1995;Dibikeetal.,2001).The
problems,thecomputationaleffortincreasestwoorthreeordersof use RBFs within the context of evolutionary algorithms was
magnitudeduetotheuseofsynthetic(insteadofhistorical)time popularized after the publication Regis and Shoemaker (2004).
series of very large length (e.g., thousands of years), in order to Other typical examples of RBFs are the Multistart Local Metric
provide estimations for probabilistic quantities (e.g., reliability, StochasticRBF(MLMSRBF)andtheConstrLMSRBF,whichhandles
risk)withsatisfactoryaccuracy.Ontheotherhand,dependingon inequalityconstraints(RegisandShoemaker,2007b;Regis,2011).
the number of parameters and the irregularity of the response Additionally,Regis(2014)andTangetal.(2012)proposedhybrid-
surface,theoptimizationalgorithmmayneedtocallthesimulation izations of the particle swarm optimization algorithm (Kennedy
modelhundredsorthousandsof times,inordertoconvergetoa andEberhart,1995)thatuseRBFstoassistthesearch.Shoemaker
good solution. Therefore, the time effort of simulation imposes a et al. (2007) developed an evolutionary algorithm that uses an
practical barrier to optimization, which is necessary to run with RBF approximation and benchmarked its performance against
significantlyrestricted“budget”,bymeansofmaximumallowable severaltestproblems,withdimensionsrangingfrom8-Dto14-D.
numberoffunctionevaluations.Considerasimulationmodelthat Finally, Regis and Shoemaker (2013) developed the DYnamic CO-
requiresapproximately1.5minforasinglesimulationrunandan ordinate Search (DYCORS) that uses Response Surface models to
optimization algorithm that requires 10,000 function evaluations handle high-dimensional expensive optimization problems.
(iterations)toapproximatetheglobalminimum.Consequently,the DYCORSwasbenchmarkedagainstotherRBF-basedalgorithmsina
procedurewouldlastmorethantendays,whichmakesitpracti- varietyoftestproblems,rangingfrom14-Dto200-D.
callyinfeasible. Comprehensive reviews of surrogate-based optimization
AccordingtoRazavietal.(2010),theapproachestoalleviatethe methodscanbefoundinthebroaderoptimizationliterature(e.g.,
computational burden imposed by time-consuming simulation Jin,2005;ForresterandKeane,2009;Jin,2011).Therearealsore-
models are classified into four main categories: (1) parallel ported several successful applications in time-demanding hydro-
computing(e.g.,Schutteetal.,2004;Chengetal.,2005;Vrugtetal., logicalproblems(e.g.,Broadetal.,2005;Mugunthanetal.,2005;
124 I.Tsoukalasetal./EnvironmentalModelling&Software77(2016)122e142
Mugunthan and Shoemaker, 2006; Regis and Shoemaker, 2007a; 2. Optimizationmethodology
Zou et al., 2007; Kourakos and Mantoglou, 2009; Tsoukalas and
Makropoulos, 2015a). Razavi et al. (2012b) summarize the use of 2.1. EvolutionaryAnnealing-Simplex
surrogate modeling techniques in water resource systems, also
classifyingtheexistingmeta-modelingframeworks. EAS1isaheuristic,population-basedglobaloptimizationtech-
It is important to remark that in the context of combined nique, originally developed by Efstratiadis and Koutsoyiannis
simulation-optimization schemes, surrogate models play the role (2002),thatcouplesthestrengthofsimulatedannealinginrough
ofblack-boxapproachesthataimestablishingadata-drivenrela- search spaces along with the efficiency of the downhill simplex
tionshipbetweenthecontrolvariablesofthesimulationmodel(i.e., method(NelderandMead,1965)insmootherspaces.Itskeyideais
explanatoryvariables)andtheobjectivefunctionoftheoptimiza- theintroductionofanexternalvariableT,whichplaysarolesimilar
tionmodel(i.e.,responsevariable).Therefore,theyclearlydonot totemperatureinareal-worldannealingprocess,anddetermines
intendtoreproducethedynamicbehavioroftheoriginalsimula- the degree of randomness of the search procedure. This is
tionmodel(Razavietal.,2012b).Infact,thetaskofreproducingthe expressedthroughastochastictermthatisrelativetotemperature
dynamicbehaviorofthesimulationmodelisperformedbyaquite and is added to the initial objective function f(x), thus getting a
different surrogate modeling approach, generally referred to as modifiedfunctiong(x)¼f(x)±uT(whereuisavectorofuniformly
model reduction or reduced-order modeling. This yields a low- distributed random numbers). Search is based on an evolving
order, dynamic “equivalent” of the simulation model, by preser- populationoffeasiblepoints,wherecriticaldecisionsaredrivenby
ving,tosomeextent,thestate-spacerepresentationoftheoriginal the modified function. The genetic operators are either quasi-
model and allowing a physical interpretation of its structure stochastic geometric transformations, inspired by the downhill
(Castellettietal.,2012a,2012b). simplex method, or fully-probabilistic transitions (mutations). As
This paper introduces the Surrogate-Enhanced Evolutionary searchproceeds,thesystemtemperaturereducesaccordingtoan
Simplex-Annealing approach (SEEAS), which is a novel global adaptive annealing cooling schedule, and all transitions become
optimizationalgorithm,focusedontime-expensivefunctions.Our moredeterministic.
motivationarisesfromchallengingsimulation-optimizationprob- EAS has been successfully employed in several hydrological
lemsthatarecommonlyfoundinwaterresources,andtheyimpose, applications(e.g.,Rozosetal.,2004;Nalbantisetal.,2011;Kossieris
intheeverydaypractice,verylimitedcomputationalbudgets,e.g., etal.,2013;Efstratiadisetal.,2015).Ithasbeenalsoincorporated
offewhundredfunctionevaluations.SEEAShasbeendesignedfor within advanced modeling tools, i.e., Hydronomeas (Efstratiadis
bothtypesofsuchproblems,i.e.,decision-makingandcalibration, et al., 2004), Hydrogeios (Efstratiadis et al., 2008) and HyetosR
sufferingfromdifferentpeculiaritiesandcomplexities,whicharein (Kossierisetal.,2012)tosolvechallengingsimulation-optimization
turnreflectedinthedifferentgeometryoftheassociatedresponse problems.Theoriginalalgorithmhasbeenalsoadaptedtohandle
surfaces. multiobjectiveproblems(EfstratiadisandKoutsoyiannis,2008)and
SEEASisbuiltupontheEvolutionaryAnnealing-Simplex(EAS) stochastic (i.e., noisy) objective functions (Kossieris et al., 2013).
method (Efstratiadis and Koutsoyiannis, 2002), which is a hybrid HereweintroduceanimprovedversionofEAS,calledSurrogate-
schemecombiningglobalandlocalsearchstrategiesandassisted Enhanced Evolutionary Annealing-Simplex (SEEAS) algorithm,
by a RBF surrogate model. SEEAS uses an external archive to whichispresentedindetailherein.
maintain all visited solutions in order to formulate, update and
exploit the surrogate model during search. There are also some 2.2. Surrogate-EnhancedEvolutionaryAnnealing-Simplex
improvementsinthekeycoreofEAS,regardingthesimplextran-
sitionsandthemutationoperator.SEEASiscomparedandbench- 2.2.1. OverviewofSEEASalgorithm
markedagainsttheoriginalversionofEASandthreestate-of-the- Thealgorithmisasurrogate-enhancedextensionofEASinaway
art optimization algorithms that are mentioned before, i.e., DDS thatbuilds,maintainsandexploitssurrogatemodeling(SM)tech-
(Tolson and Shoemaker, 2007), MLMSRBF (Regis and Shoemaker, niquesthatgenerateapproximatedresponsesurfaces,whichallow
2007b), and DYCORS-LMSRBF (Regis and Shoemaker, 2013). Eval- effectively guiding search towards promising areas of the real
uationsaremadeonthebasisof12mathematicalproblems(i.e.,six responsesurface.ThemodelusedistheRBF,whichisawell-known
test functions for two alternative dimensions,15-D and 30-D), a interpolationtechnique(Fig.1,left).Duringtheiterativeprocedure,
hydrological calibration problem with 11 parameters, configured the algorithm maintains an external archive of all visited points,
withbothrealandsyntheticdata,andamulti-reservoirmanage- alreadyevaluatedthroughthe(expensive)objectivefunction.This
mentproblemwith20decisionvariables,usingsyntheticinflowsof archive is used to update the SM, in an attempt to progressively
500yearslength.Theuseofsyntheticdataisoneofthenoveltiesof provide more accurate approximations of the current region of
our testing framework. Moreover, most of the known surrogate- interest(i.e.theareaaroundthecurrentbestpoint).InSEEAS,the
based schemes have been onlyevaluated in calibrationproblems surrogatemodelhasadoublerole.Thefirstisprovidingnewpoints
andnotintime-demandingwatermanagementapplications,with thatareaddedtothecurrentpopulation,andthesecondisassisting
few exceptions (e.g., Razavi et al., 2012b; Tsoukalas and thegenetic operatorsof thedownhillsimplexschemetoidentify
Makropoulos, 2015a). The results of this extended analysis are suitable directions across the search space (e.g., favorable slopes
very encouraging, since the proposed method is effective and andnewareasofattraction).
efficient, in terms of locating a satisfactory solution as close as In order to balance exploration (i.e., detailed sampling) and
possibletotheglobaloptimum,withinreasonablecomputational exploitation(i.e.,blinduseofSM),SEEASusesaweightedmetric,
time,andclearlyoutperformstheotherexaminedapproaches,in termed acquisition function (AF), which accounts for the pre-
almostalltests. dictionsprovidedbytheSMaswellasthespreadofallpreviously
evaluatedpoints(bymeansofadistancequantity).Inoppositeto
commonpractices that use a standard expression of the AF with
constant weights, in our approach the weights are dynamically
adjusted, thus improving the efficiency of the algorithm. Details
abouttheacquisitionfunction(AF)aregiveninSection2.2.3.
1 EASandSEEASareavailableonlineat:http://www.itia.ntua.gr/en/softinfo/29/. SEEASfollowsaniterativesearchprocedure.Attheendofeach
I.Tsoukalasetal./EnvironmentalModelling&Software77(2016)122e142 125
Fig.1. Approximatedsurface(RBF)ina2-Dexample(Ackleyfunction)usingallavailablesamplepoints(leftpanel).Therightpaneldemonstratesarandomlyselectedsimplexand
themodifiedsurrogate-enhancedreflectionmovementusingcandidatepointsonthelineformedfromthesimplexcentroidandthemaximumreflectionpoint.Thesimplexis
reflectedatthecandidatepointwiththeminimumfunctionvalue.
iteration cycle (or generation, according to the terminology of
evolutionarytheory),weobtainatleastonenewpointthatenters XNs
sðxÞ¼ l 4ðjjx(cid:2)xjjÞþpðxÞ (1)
thepopulationandreplacesoneofitsexistingmembers.Atypical i i
iterationcycleofSEEASstartsbyfittingthesurrogatemodeltothe i¼1
cthurroreungthpLoaptiunlaHtiyopne(ricnuitbiealSlya,mthpilsinpgo,pLuHlaSt)i.oNneixstr,awnedormunlyagneinnetreartneadl whereli∊R,4isabasisfunctionoftheform4(r)¼r3,jj.jjisthe
Euclideandistance(norm)andp(x)isapolynomialtailoftheform
gElAoSb)aalcorpotsismtihzeatsiounrraolggaotreitrhemspo(pnasretiscuurlfaarclye,,tuhseinogriagsinoabljveectrisvioenthoef rpa(mx)e¼terbsΤxl,þb,aa,nwdhaearreebde¼te(rbm1,in…e,dbbny)Tsaonlvdinag∊thRe.lTinheeamrsoydsetelmpa:-
acquisitionfunction(AF),inordertolocateacandidatesolutionto
(cid:2) (cid:3)(cid:2) (cid:3) (cid:2) (cid:3)
enterthepopulation(providedthatthissolutionoutperformsthe
F P l y
currentworstpoint).Thereafter,wefollowasearchprocedurethat PT 0 c ¼ 0 (2)
is mostly based on the genetic operators of EAS, enhanced by
sinurorTorhdgeaegtree-tnaoesrseainslthiedadenascteiesptthsoeiuntcsiuliimrzreepntlehtxek-ibnnaofsowerdlmedtargtaienonsifnogratmihneaetdisoebnlyesc.tthioenSMof, Nwsh(cid:3)er(enFþis1)anmNatsr(cid:3)ixN,tshmeaitthrixrowwitohfewlehmicehnitss(41,ijx¼Ti)4,(ljj¼xie(l1x,j…jj),,PlNiss)Ta,
simplextransitions.Acharacteristicexampleinvolvingthereflec- c¼(b1,…,bn,a)T,andy¼(y1,…,yNs)T.Wementionthatthematrix
ofEq.(2)isinvertibleifandonlyifRank(P)¼nþ1(Powell,1992).
tionstepisillustratedinFig.1,right(forsimplicity,wedemonstrate
the predictions of the surrogate model and not the AF). In the
originalversionofEAS,afterspecifyingthedirectionofreflection
(definedbythedifferencebetweentheworstvertexofthesimplex 2.2.3. Acquisitionfunction
andthecentroidofallrestvertices),thealgorithmemploysablind Acquisition functions (AF) are well-established techniques,
trial-and-error procedure, i.e., it generates subsequent random aiming to balance exploration-exploitation in surrogate-based
pointsalongthisdirectionandevolvesaccordingtotheirvalues.In optimization algorithms (e.g., Sasena et al., 2002; Forrester and
this scheme, the original objective function is called whenever a Keane, 2009). SEEAS implements a novel scheme, in which the
newtrialpointisgenerated.Sincetheexpansioncontinuesaslong weights are automatically adjusted during the iterative process,
as the function value improves, this procedure may be quite according tothe current number of function evaluations and the
expensive,intermsoffunctionevaluations.Inopposite,inSEEAS maximumallowednumberofevaluations.
weemployacandidatescreeningprocedureusingtheSM,which ConsiderasetofNspoints,xjs,withknownresponsevalue,f(xjs),
allows making multiple trials with negligible computational cost andanothersetofNcpointsxic,withapproximatedresponsevalues
andguidingsearchusingallpriorinformation.Similarscreeningis s(xi).Thelatterareconventionallycalledcandidatepoints,inthe
c
employed within all simplex transformations (except shrinkage), sense that they are used within infilling or internal search pro-
thusprovidingsignificantaidtotheassociateddecisions. cedures,e.g.,selectionofthemostappropriatereflectionpointin
the graphical example of Fig. 1. The acquisition function is esti-
matedasfollows:
2.2.2. Surrogatemodel(RBF) Step A: Standardize the approximated response values of all
SEEASimplementstheRadialBasisFunction(RBF)interpolation candidate solutions by setting s*(xi) ¼ [s(xi) e smin]/
c c
method(Powell,1992;Buhmann,2003),andmorespecificallythe [smaxesmin],wheresminandsmaxarethecorrespondingmini-
RBFwithcubicbasisfunctionsandlinearpolynomialtail.Thisisa mumandmaximumvalues.
commonly used surrogate model of proven effectiveness, as re- Step B: Calculate the minimum Euclidean distance of each
portedinnumerousstudies(e.g.,Mugunthanetal.,2005;Regisand candidatepointxi fromallpreviouslyevaluatedpoints,xj,i.e.,
c s
Shoemaker, 2007a, b; Shoemaker et al., 2007; Regis and di ¼ d*(xic) ¼ min1 (cid:4) j (cid:4) NSjjxic e xjsjj, and standardize them by
Shoemaker,2013;MüllerandShoemaker,2014). settingdi*¼(diedmin)/(dmaxedmin),wheredminanddmaxare
ThecomputationalprocedureofRBFisthefollowing.GivenNs thecorrespondingminimumandmaximumdistances.
samplesx∊Rnwithresponsey,wegetthepairs(xi,yi).Thepre- StepC:CalculatetheweightedvalueofAFforeverycandidate
dictions(x)ofRBFmodelatsamplepointxisgivenby: pointusingtheformula:
126 I.Tsoukalasetal./EnvironmentalModelling&Software77(2016)122e142
(cid:4) (cid:5) (cid:4) (cid:5) Step3:AnewpointxpisgeneratedbyminimizingAF,usingthe
AF ¼ws* xi þð1(cid:2)wÞd* xi (3) originalversionofEASforinternaloptimization.Thenewpoint
i c c
is evaluated through f(x) and replaces the worst point of the
wherewisadimensionlessweightingcoefficient,ensuringbalance currentpopulation,ifthelatterisworse(higher)thanf(xp).
Step 4: A set of n þ 1 points is randomly selected from the
between exploitation and exploration. To finalize the infilling
currentpopulation,inordertoformulatetheverticesofasim-
routine,thecandidatewiththeminimumAFvaluewillbeselected
andassessedthroughtheobjectivefunction.Asmentionedbefore, plexinthen-dimensionalsearchspace,symbolizedS¼[x1,x2,
the minimization of the AF across the surrogate search space is …,xnþ1].TheelementsofSaresortedsuchasf(x1)corresponds
carriedoutthroughtheoriginalEASalgorithm. to the best (lowest) and f(xn þ 1) to the worst value of the
objectivefunction.
Step5:Fromthesubset[x2,…,xnþ1]weselectacandidatepoint
xw to be replaced in the population, based on the modified,
quasi-stochasticobjectivefunction:
2.2.4. DetaileddescriptionofSEEAS
Let f(x) be anonlinearobjectivefunction in the feasible space gðxÞ¼fðxÞþuT (6)
xL (cid:4) x (cid:4) xU, where x is an n-dimensional vector of continuous
control variables (in practice, f(x) represents the performance
where u is a uniform random number in the interval [0, 1]. By
measureofasimulationmodel).Forconvenience,wesearchforthe
addingthestochasticcomponentuTtotheobjectivefunctionf(x),
globalminimumoff(x),allowingabudgetofMFEfunctionevalu-
thealgorithmbehavesasinbetweenrandomanddownhillsearch.
ations.Thealgorithmusestwoarchives.Thefirstisthepopulation
Attheearlystagesofoptimization,whentemperatureisstillhigh,
P[t],whichisevolvedduringthesearchprocedure(wheretdenotes
any point except for the best one can be replaced. On the other
tehxeteirtnearaltaiorcnhcivyecleA[ot]r, gwehniecrhatcioonn)t,aainnsdatlhlevsiseictoedndpiosinthtsefsroo-mcaltlhede hand,inthelimitingcaseT/0,theactuallyworstpoint,i.e.,xnþ1,
isreplaced,asconsideredintheoriginaldownhillsimplexmethod.
beginningoftheoptimization(t¼0),includingthemembersofthe
cthuerroebnjtecptoivpeulfautnioctni.oWnfh(xe)n,eitveenrtaernsetwhepaoricnhtivxeiAs[et]v(atlhueaatercdhtihvreomugahy Step6:AsetofNrtrialpointsxkcraregeneratedbyreflectingthe
simplexaccordingtheformula:
beupdatedseveraltimeswithinageneration).Atthebeginningof
ceoanchsidneerwinggtehneercautriorenntt,eltehmeesnutrsroogfaAt[et].mThoedseilzeisorfet-heevaploupauteladtiobny xkcr ¼gþð0:5þdkÞðg(cid:2)xwÞ (7)
ism(cid:5)nþ1(i.e.,theminimumnumberofpointsrequiredtofita wheregisthecentroidofthesubset[x2,…,xnþ1]anddkisascale
RBFwithlinearpolynomialaswellastoformulateasimplexinthe coefficientequallyspreadintheinterval[0,1],thusdk¼(ke1)/
n-dimensionalspace),andremainsconstant,whilethesizeofthe (Nre1),fork¼1,…,Nr.Amongallcandidates,weselecttheone
external archive progressively increases, thus ensuring more ac- thatminimizesAF,whichwewillnextcallthereflectionpoint,xr.
curateapproximationsoftheresponsesurfaceand,consequently, The reflection point is evaluated on the basis of the objective
more reliable predictions. The initial population P[0] is generated
functionandenterstheexternalarchive.
viatheLatinHypercubeSampling(LHS)technique,whichensures
Asaptpisafraecntotlryy,tshpereianditiaaclraorscshtihveeAfe[0a]siisbliedesnptaicceal(tGoiuPn[0t]a. et al., 2003). Step7:Iff(xr)<f(xw),wereplacexwbyxrinthepopulationand
move to steps 8a or 8b, according tothe outcome of its com-
auxSiliimariylaprlayratomEeAteSr,,tTh[et],scuarlrleodgatteem-epnehraantucreed.Tahlgeocroitnhcmepatlsooriugsineastaens parisonwiththecurrentbestvertex,i.e.,f(xr)<f(x1).Otherwise,
wemovetostep9,todecidewhetherxrshouldbeacceptedor
from simulated annealing, where the key role of temperature is
withdrawn,thusseekinganothercandidate.
ensuringbalancebetweenrandomnessanddeterminism.InSEEAS, Step8a:Iff(xr)<f(x1),thevectorxrex1definesadirectionof
temperatureisdynamicallyadjusted(i.e.,reduced)usingempirical
rules,consideringtheextremevalues,f½t(cid:6) andf½t(cid:6) ,ofthecurrent minimization.Weremarkthatthedetectionofdownhillslopes
min max inhigh-dimensionalspacesofcomplexgeometryisnotanoften
populationP[t],andadimensionlessprogressindex,definedas
case.Thismakesessentialtotakeadvantageinordertoaccel-
PI¼logðFEÞ=logðMFEÞ (4) eratethesearchprocedure,byemployingasequenceofNetrial
expansionstepsthroughtherecursiveformula:
whereFEisthecurrentnumberoffunctionevaluationsandMFEis
the maximum allowable number of FE, which is a user-specified xkce¼gþdkðxr(cid:2)gÞ (8)
termAitnyaptiicoanlcitreitreartiioonn.cycle ofSEEAS,an outlineofwhichis illus- wheredkisascalecoefficientgivenbydk¼dke1þ(ke1)/(Nee1),
tratedinFig.2,comprisesthefollowingsteps(generationindextis fork¼1,…,Ne.TheexpansioncontinuesaslongastheAFvalueis
improved(oruntilreachingtheboundsofthefeasiblespace).The
omittedforsimplicity):
optimal(intermsofAF)trialpoint,xe,iskeptintheexternalarchive
Step1:Theinterpolationsurfaces(x)isupdatedusingthecur- andreplacesxrinthecurrentpopulation,providedthatf(xe)<f(xr).
Inthatcase,thealgorithmmovestostep12tofinalizethecycle.
rentinformationstoredintheexternalarchiveA(i.e.,allpoints
evaluatedsofarthroughtheoriginalobjectivefunction).
Step2:TheweightingcoefficientoftheAFisupdatedusingthe Step8b:Iff(xr)>f(x1),weattemptdetectingapromisingsolu-
empiricalformula: tion in theneighborhoodofx1,byemploying Nc trialcontrac-
tionsofthesimplexintheintervalbetweenthecentroidandthe
reflectionpoint,accordingtotheformula:
w¼max½0:75; minðPI; 0:95Þ(cid:6) (5)
The above formula ensures that at the early stages of optimi- xkcc¼gþð0:25þ0:5dkÞðxr(cid:2)gÞ (9)
zation,moreweightisgiventoexploration(upto0.25),butgrad-
uallyitscontributiondiminishesthusnotexceeding0.05. wheredk¼(ke1)/(Nce1),fork¼1,…,Nc.Theoptimal(intermsof
I.Tsoukalasetal./EnvironmentalModelling&Software77(2016)122e142 127
Fig.2. OutlineofSEEASalgorithmfollowingthestepsexplainedinSection2.2.4(*denotestheuseofthesurrogatemodelwithintheassociatedsimplextransformations).
AF)trialpoint,xc,iskeptintheexternalarchiveandreplacesxrin Thethresholdof0.50prohibitsafastreductionoftemperature
thecurrentpopulation,providedthatf(xc)<f(xr).Inthatcase,the and therefore maintains enough randomness within decisions,
algorithmmovestostep12tofinalizethegenerationcycle. which in turn prohibits early convergence to local optima. After
reducingT,theiterationcycleisfinalized(step12).
Step9:Iff(xr)>f(xw),weusethemodifiedobjectivefunction(7)
todecidewhetheremployinginsidecontractionofthesimplex, Step 10b: The reflection point xr is accepted although being
thus seeking for a potential local optimum, or expanding to- worsethanxw.Next,Nuuphill(i.e.,maximization)movements
wards a non-optimal (i.e., uphill) direction, in an attempt to areperformedusingthesameformulawithmultipleexpansion
escape from the current area of attraction. In this respect, if (Eq. (9)), in an attempt to pass the hill and discover adjacent
g(xr)>g(xw)wemovetostep10a,otherwisewemovetostep regions of attraction. This geometrical transformation was
10b. introduced by Pan and Wu (1998), to facilitate the simplex
Step10a:WerejectxrandimplementNctrialinsidecontractions escapingfromlocalminima.Similarlytoprevioussteps,weuse
of the simplex in the interval between the centroid and the theAFtodeterminetheoptimumuphillpoint,xu.Iff(xu)<f(xr),
worstpoint,accordingtotheformula: thispointiskeptintheexternalarchiveandreplacesxrinthe
current population, while the algorithm moves to step 12 to
xkcc¼g(cid:2)ð0:25þ0:5dkÞðg(cid:2)xrÞ (10) finalize the generation cycle. Otherwise, none of the simplex
transformations results to a better solution than the worst
wheredk¼(ke1)/(Nce1),fork¼1,…,Nc.Theoptimal(intermsof vertex xw, thus the last option is to attempt a pure stochastic
AF)trialpoint,xc,iskeptintheexternalarchiveandreplacesxwin generator,referredtoasmutation(step11).
thecurrentpopulation,providedthatf(xc)<f(xw).Otherwise,the Step11:Weseekarandompointoutofthetypicalrangeofthe
simplexshrinkstowardsthebestvertexx1,suchas: current population, defined on the basis of the mean, mP, and
standarddeviation,sP,ofallmembersofP.Inthisrespect,we
xs;i¼0:5ðx1þxiÞ for i¼2;…;nþ1 (11) [gmePneerastPe,maPnþorsmP]a,lwlyh-dicihstrisibauctceedptpeodiniftf(xxmm)o<utf(oxfr).thOethienrtwerivsael,
We remarkthat the above transformation is the sole evolving weaccountforauser-specifiedmutationprobabilitypminorder
mechanismofthealgorithmallowingthesimultaneousgeneration toacceptornottherandomlygeneratedpoint,xm,andreplacing
of multiple points; particularly, n new points are generated that xr in the current population. Anyway, since xm is evaluated
replaceallpreviousverticesinthecurrentpopulation.Thiscanbe throughtheobjectivefunction,itenterstheexternalarchive.
consideredasmilestoneofthesearchprocedure,inthesensethata Step 12: Considering the new member (or members, in the
local minimum, lying in the neighborhood of x1, has been sur- particularcaseofsimplexshrinkage)ofthepopulation,were-
rounded. This is the time to reduce the temperature of the opti- evaluatethecurrentminimum,xmin,andmaximum,xmax,and
mizationsystembyareductionfactorj.IncontrasttoEAS,wherej their function values, fmin and fmax. We also re-evaluate the
isaconstantparameteroftheannealingcoolingschedule,usually currentnumberoffunctionevaluations,FE,andcheckwhether
takingvaluesintotheinterval0.90e0.99,initssurrogate-enhanced thishasn'texceededtheterminationcriterion,MFE.Finally,we
versionjisautomaticallyadjustedtoalsoaccountfortheprogress re-evaluate the temperature so that T (cid:4) x(fmax e fmin), where
indexPI,usingthefollowingexpression: x(cid:5) 1 is a user-specified parameterof the annealingschedule,
usuallysetbetween2and5.ThisrestrictionpreventsTtaking
j¼maxð1(cid:2)PI; 0:50Þ (12)
128 I.Tsoukalasetal./EnvironmentalModelling&Software77(2016)122e142
Table1
Configurationofbenchmarkingsuite.
Problem Algorithms Numberofcontrol Max.Function Independentrunswithrandom Population Surrogatemodel(metamodel)
variables,n evaluations(MFE) initialpopulations size
Testfunctions All 15 500,1000 30 32 RBFwithcubicbasisfunctionsand
Testfunctions All 30 500,1000 30 62 linearpolynomialtail
Modelcalibrationwith All 11 500,1000 30 24
realdata
Toycalibrationwith All 11 500,1000 30 24
syntheticdata
Multireservoir SEEAS,DYCORS, 20 500 10 42
management MLMSRBF
problem
extremelyhighvalues,whichwoulddeterioratetheefficiencyof dependingonMFE(inparticular,SEEASusestheprogressindexPI,
SEEAS,asfarassearchwouldbecometoorandom. definedinEq.(4),withintheannealingcoolingschedule).Finally,
forthethreesurrogate-basedmethods(SEEAS,DYCORS,MLMSRBF)
Torunthealgorithm,itisessentialprovidingvaluesforallinput weemployedthesamemetamodel(RBFwithcubicbasisfunctions
arguments, which are the number of desirables steps within and linear polynomial tail), thus ensuring similar computational
different simplex transitions (Nr, Ne, Nc, Nu), the mutationproba- effortforbuilding,updatingandexploitingtheRBF(Razavietal.,
bility pm, and the adjusting factor x of the annealing cooling 2012a). We remark that in real-world problems the effort of the
schedule.Recommendedvalues,alsousedinallnextbenchmarking optimization routines (including metamodel fitting) is much less
tests,areNr¼Ne¼Nc¼Nu¼20,pm¼0.10andx¼2.Thesevalues than the effort of simulation, and therefore the runtime of the
weredeterminedonthebasisofextendedinvestigationswithinthe overallsearchprocedureispracticallyrelativetoMFE.
developmentofSEEAS,andtheyhavebeenalsovalidatedthrough AllcomputationswereimplementedinMATLABmathematical
thesensitivityanalysisofSection4.4. environmentusinga3.0GHzIntelCorei5processorwith4GBof
RAM, running on Windows 8 OS. For the SEEAS method we
3. Benchmarkingmethodology employedthetypicalinputargumentsgiveninSection2.2.4,while
fortheotheralgorithms,i.e.,EAS,DDS,MLMSRBFandDYCORS,we
3.1. Benchmarkingprotocol used the default values suggested in the associated articles
(Efstratiadis and Koutsoyiannis, 2002; Regis and Shoemaker,
To assess the performance of SEEAS we compared it with the 2007b;TolsonandShoemaker,2007;RegisandShoemaker,2013).
originalversionofEASaswellasthreestate-of-the-artoptimiza-
tion algorithms, which are synoptically presented in Section 3.3. 3.2. Performanceevaluationapproach
Twoofthebenchmarkalgorithms,i.e.,DYCORSandMLMSRBF,are
surrogate-assisted, while EAS and DDS do not employ surrogate Following the ideas of Razavi et al. (2012a) and Matott et al.
modelsthroughsearch. (2012),afterimplementingallrunsforeachspecificoptimization
Avarietyoftestproblemswereexamined,theoreticalaswellas problemsolvedwithaspecificalgorithm,weplottedthecumula-
real-world. Briefly, the hereafter called benchmarking “suite” in- tive distribution function (CDF) of the optimal values of f(x) ob-
cludessixmathematicaltestfunctions,formulatedwith15and30 tained within the specific budget. In order to quantify the
controlvariables,ahydrologicalcalibrationproblemwithrealand probability of attaining an equal or better solution, we used the
syntheticdata,andatime-expensivemulti-reservoirmanagement conceptofstochasticdominance(SD),introducedbyLevy(1992),to
problem(6(cid:3)2þ1(cid:3)2þ1¼15problems,intotal). comparetheCDFsofthealgorithms.LetFАandFВbetheCDFsof
To ensure fair comparison and safely infer about the perfor- algorithmsAandB,respectively.Assumingtheminimizationofa
manceofthealgorithmsweattemptedtoensureasmuchassimilar randomquantityq,weassumethatAdominatesBifFА(q)>FB(q)
configurations, as summarized in Table 1. In all problems we for all q, and vice versa. On the contrary, if the two CDFs are
employed multipleindependent runs, using the same population intersectedatsomepointqu,thenSDisnotapplicable.Inthiscase,
size and the same random generation technique, i.e., LHS. The weevaluatedthemedianpoint,i.e.,theonewith50%probabilityof
populationsizewassetequaltom¼2(nþ1),asrecommendedby exceedance,andconsideredasbetterthealgorithmwiththebest
Regis and Shoemaker (2007a, 2013), where n is the problem performanceatthispoint.Infact,toensurethatthedifferenceof
dimension(i.e.,thenumberofcontrolvariables).Weremarkthat thetwoalgorithmsatthepointofinterestisstatisticallysignificant,
otherresearchersrelatetheinitialpopulationsize(alsoreferredto we employed the non-parametric ManneWhitney U-test (MWU;
as design of experiment, DoE) to the available computational MannandWhitney,1947).ThenullhypothesisoftheMWUtestis
budget, quantified in terms of MFE, in order to design a more thatdatainFАandFВaresamplesfromcontinuousdistributions
detailedmetamodel;forinstance,Razavietal.(2012b)suggestthat withequalmedians.TheconfidenceleveloftheMWUtestwasset
m ¼ max[2(n þ 1), 0.1 MFE]. However, in our tests we avoided to95%.
associatingmwithMFE,inordertoinvestigatetheimpactsofthe
problem dimension to the performance of the examined algo- 3.3. Briefdescriptionofbenchmarkoptimizationalgorithms
rithms.Furthermore,wepreferredsavingresourcesfortheevolu-
tionaryprocedure,insteadofspendinganon-negligiblepartofour 3.3.1. DynamicallyDimensionedSearch(DDS)
budgettotheinitialDoE. Dynamically Dimension Search2 (DDS), is a stochastic, single-
Eachproblembutthelastwassolvedconsideringtwoalterna- solution based algorithm, developed by Tolson and Shoemaker
tivecomputationalbudgets,MFE(500and1000).Werunalltests
withtwodifferentbudgets(insteadofthemaximumofthem)since
all examined algorithms (except EAS) involve parameters 2 https://github.com/akamel001/Dynamic-Dimension-Search.
I.Tsoukalasetal./EnvironmentalModelling&Software77(2016)122e142 129
(2007)tolocatenear-optimalsolutionswithfewfunctionevalua- 4. Testfunctions
tions. DSS is designed to search globally at the early stages and
more locally when approaching a user-specified number of 4.1. Setupofoptimizationproblems
maximumfunctionevaluations(MFE).Itevolvesbyperturbingthe
current best solution in randomly selected dimensions, using an Thefirstsuiteofbenchmarkproblemsinvolvestheoptimization
evolutionary operator based on the normal distribution. The of six well-known mathematical problems (test functions),
probability of selecting a dimension toperturb is proportional to combining two alternative formulations in terms of number of
thecurrentnumberoffunctionevaluationsandMFE.Thetransition variables (n ¼ 15 and 30), and two algorithmic configurations in
from global to local search is employed by dynamically reducing termsofMFE(500and1000).Thissettingallowedforassessingthe
thenumberofperturbeddimensions.Intheliteraturearereported performanceofthealgorithmsagainstincreasinglevelsofdimen-
several successful applications of DDS (e.g., Tolson et al., 2009; sionality and increasing computational budget. Considering two
Razavi et al., 2010; Matottetal., 2012; Razavi etal., 2012a; Regis alternative dimensions and two computational budgets, we
andShoemaker,2013). configured four different problems for each test function, i.e., 24
optimization problems in total. According to the benchmarking
protocolexplainedinSection3.1,forallproblems,weemployed30
independentruns,thusrandomlychangingtheinitialpopulationof
3.3.2. MultistartLocalMetricStochasticRBFalgorithm(MLMSRBF) eachsearchexperiment.Thepopulationsizeofallalgorithmswe
Regis and Shoemaker (2007b) developed the Multistart Local set equal to 32 and 62, for the 15-D and 30-D formulations,
Metric Stochastic RBF3 (MLMSRBF), which is surrogate-assisted respectively.
optimization algorithm that can be considered as extension of Table 2 summarizes the main characteristics of the examined
DDS.ThefirststepistheimplementationoftheinitialDoEtofitthe test functions, which represent search spaces of different
surrogate model (particularly, RBF), which evolves by perturbing complexity. Two of them (Sphere and Zakharov) are unimodal,
thecurrentbestpoint(similartoDDS),usingnormaldistribution whiletherestaremultimodal(Ackley,Griewank,Rastrigin,Levy).
withzeromeanandaspecifiedcovariancematrix.Additionally,in Inallcasestheglobalminimumisknownandequaltozero.The
ordertolocatepromisingcandidates,thealgorithmusesametric analyticalexpressionofthetestfunctionsandtheboundsoftheir
thatbalancestheRBFpredictionandtheminimumdistancefrom variablesaregivenintheAppendix.
previously evaluated points (this is similar to the acquisition
functionintroducedin2.2.3,butwithconstantweights).Theglobal 4.2. Statisticalevaluationofoptimalsolutions
character of the algorithm is further enhanced by implementing
multipleDoEs.Thismultistartstrategyisenabledonlyifthealgo- Aninitialassessmentoftheperformanceofthefiveexamined
rithmappearstohavebeentrappedtoalocalminimum.Regisand algorithmswasmadeonthegroundsofmeanandstandarddevi-
Shoemaker (2007b) demonstrated the efficiency of MLMSRBF in ationofthebestfunctionvaluesobtainedfromeachoptimization
several benchmark problems, including 17 multimodal test func- set(i.e.,30independentrunsofthealgorithm).Theclosesttozero
tionsanda12-dimensionalgroundwaterbioremediationproblem. is themean and the lowest the standard deviation indicates that
Intheliteraturearealsoreportedothersuccessfulapplicationsof thealgorithmreachesthetheoreticaloptimumwithhighaccuracy
the method (e.g., Mugunthan et al., 2005; Mugunthan and andreliability.
Shoemaker,2006;RegisandShoemaker,2013). The statistical superiorityof SEEAS is exhibited in all problem
configurations,asshowninTables3and4,forproblemdimensions
n¼15and30,respectively.Specifically,forthe15-Dformulation
(Table 3), SEEAS achieves the best performance (i.e., the lowest
mean) in three out of six (OF1, OF3, OF6) and four out of six
3.3.3. DYnamicCOordinateSearch-MultistartLocalMetric
problems(OF1,OF2, OF3,OF6), forMFE¼ 500 and1000,respec-
StochasticRBF(DYCORS-LMSRBF)
tively. By doubling the dimensionality of the test functions to
The DYCORS framework was recently proposed by Regis and
n¼30,thussignificantlyincreasingthecomplexityoftheassoci-
Shoemaker (2013) for surrogate-based optimization of high-
ated optimization problems, SEEAS outperforms the other algo-
dimensionalexpensivefunctions.Theauthorspresentedtwover-
sions, DYCORS-LMSRBF and DYCORS-DDSRBF.4 The former is rithmsinfouroutofsix(OF1,OF2,OF3,OF6)andthreeoutofsix
problems (OF1, OF3, OF6), for MFE ¼ 500 and 1000, respectively
extension of LMSRBF and the latter is a surrogate-assisted DDS
(Table 4). Considering all alternative configurations, SEEAS is
(hereweuseDYCORS-LMSRBFthatperformedslightlybetterthan
optimalfor14outof24problems,DYCORSandEASareoptimalfor
DYCORS-DDSRBF). DYCORS employs a strategy similar to DDS by
4outof24,andDDSisoptimalfor3outof24.MLMSRBFdoesnot
dynamically and probabilistically reducing the number of per-
outperforminnoneofthe24testproblems.
turbeddimensionsuntilreachingtheMFE.Inordertogeneratetrial
Asexpected,theincreaseofcomputationalbudgetfrom500to
candidate points (on the selected/perturbed dimensions) the al-
1000 improves the performance of all algorithms. In general, the
gorithmusesanormaldistributionwithzeromeanandstandard
deviation sn, but this does not remain constant, since sn is mostsignificantimprovementisachievedbyEASandDDS,whichis
reasonablesincethesealgorithmsarenotsurrogate-assisted,thus
dynamicallyadjusted tocontrol the range of perturbation. More-
theyarebydefinitiondesignedtoproceed slowerthantheother
over,DYCORS-LMSRBFiscyclingthroughasetofweightsinorderto
schemes. The convergence behavior of the algorithms is further
balanceexploration and exploitation of the surrogatemodel. The
investigatedinnextsection.
authors assessed the performance of the two algorithms against
It is also worth mentioning that all algorithms exhibit poor
severaloptimizationschemesinavarietyoftestproblems,among
performanceagainstfunctionsOF4(Zakharov)andOF5(Rastrigin),
whicha14-Dhydrologicalcalibrationproblem.
sincetheyfaillocatingsatisfactorysolutionsforthegivenbudgets.
Inparticular,theplate-shaped valleyofZakharovfunctionmakes
extremely difficult fitting metamodels, which degenerates to hy-
3 https://courses.cit.cornell.edu/jmueller/or http://people.sju.edu/~rregis/pages/
software.html. perplanewithpracticallyzeroslopes.ItisnotsurprisingthatEAS
4 https://courses.cit.cornell.edu/jmueller/. ensures the best solutions, although these are still far from the
130 I.Tsoukalasetal./EnvironmentalModelling&Software77(2016)122e142
Table2
Summarycharacteristicsoftestfunctions(theirreferencesaregiveninAppendix).
Problem Testfunction Responsesurfaceproperties
OF1 Sphere Unimodalandconvex
OF2 Ackley Multimodalwithmanylocalminima
OF3 Griewank Multimodalwithmanyregularlydistributedlocalminima
OF4 Zakharov Unimodalwithaplate-shapedvalley
OF5 Rastrigin Multimodalwithmanylocalminima
OF6 Levy Multimodalwithmanylocalminimaandparabolicvalleys
Table3
Meanandstandarddeviationofbestsolutionsin15-Dtestproblems(optimalresultsarehighlighted).
MFE Testfunction EAS DDS SEEAS DYCORS MLMSRBF
Mean StDev Mean StDev Mean StDev Mean StDev Mean StDev
500 OF1 1.938 0.978 0.852 0.479 0.002 0.001 0.002 0.001 0.019 0.014
OF2 7.159 1.723 6.025 1.314 0.812 0.233 0.809 0.372 2.231 0.658
OF3 7.682 2.997 2.626 1.269 0.538 0.118 0.885 0.084 1.085 0.052
OF4 39.434 14.894 137.447 52.366 59.144 28.023 158.669 47.788 150.411 49.875
OF5 86.245 14.148 24.887 7.081 46.268 15.359 38.958 12.340 45.920 18.803
OF6 1.905 0.877 0.681 0.314 0.203 0.105 1.208 1.406 1.344 2.129
1000 OF1 0.378 0.177 0.150 0.079 0.001 0.001 0.001 0.000 0.011 0.007
OF2 3.523 0.936 3.847 0.528 0.437 0.208 0.607 0.092 1.862 0.556
OF3 2.444 1.061 1.505 0.299 0.368 0.140 0.809 0.082 1.040 0.037
OF4 26.828 17.895 97.541 38.226 41.290 26.639 121.266 36.925 121.359 37.730
OF5 59.735 17.012 11.233 3.136 29.733 12.838 33.585 13.490 35.784 11.031
OF6 0.767 0.292 0.234 0.104 0.124 0.060 0.536 0.860 0.524 0.863
Table4
Meanandstandarddeviationofbestsolutionsin30-Dtestproblems(optimalresultsarehighlighted).
MFE Testfunction EAS DDS SEEAS DYCORS MLMSRBF
Mean StDev Mean StDev Mean StDev Mean StDev Mean StDev
500 OF1 4.305 1.163 9.516 2.737 0.019 0.006 0.083 0.034 0.739 0.708
OF2 9.923 1.160 12.872 1.329 1.878 0.301 4.297 3.721 6.193 4.362
OF3 17.866 3.455 38.398 12.050 0.782 0.118 1.265 0.079 3.459 1.927
OF4 117.821 28.757 562.145 113.230 173.240 44.185 472.815 90.897 575.424 174.073
OF5 228.693 18.442 132.149 24.567 122.658 19.427 112.046 23.076 165.437 46.846
OF6 6.338 2.652 15.823 5.481 0.659 0.184 3.407 2.540 7.326 10.944
1000 OF1 2.529 0.933 2.112 0.791 0.006 0.004 0.011 0.004 0.358 0.177
OF2 6.516 0.845 7.670 0.924 1.206 0.297 1.085 0.168 3.643 1.103
OF3 8.836 2.617 8.273 2.679 0.549 0.093 1.020 0.026 2.420 0.713
OF4 94.598 20.317 412.238 118.573 151.472 54.097 403.812 93.081 491.425 146.097
OF5 198.335 16.587 71.598 15.028 98.371 19.505 85.267 22.956 134.864 39.193
OF6 2.683 0.736 3.921 2.215 0.443 0.126 4.213 5.440 2.865 4.583
theoretical optimum. EAS has been designed to also handle flat allow implementing steep downhill transitions. In general, the
response surfaces, which are often met in water management greatadvantageofthesimplex-basedtransitionsistheindirectuse
optimizationproblems,asfurtherexplainedinSection6.3.Onthe oftheconceptofgradient,whichfavorsquicklocationofregionsof
other hand, DDS is the algorithm that generallyensures the best attraction of local optima. This is of particular importance in
solutionoftheRastriginproblem.Again,thisisnotsurprising,since computational expensive problems, where the algorithm should
thesearchspaceofthisfunctionisextremelyrough,withmultiple quicklydetectpromisingdescentdirections.Infact,SEEASisclearly
localminima,thusthemoststochasticofallschemesisexpectedto superior tothe other two surrogate-assisted algorithms (DYCORS
bethemostefficient. and MLMSRBF) in all problems, except for Rastrigin. The most
impressive case is the Levy problem, whereSEEAS locates avery
good solution after the first one hundred of function evaluations
4.3. Evaluationofconvergencebehavior (Fig.9a),whilethemeanbestvaluefoundbyotheralgorithmssofar
iseventwoordersofmagnitudehigher.Similararetheresultsfor
Inordertofurtherinvestigatetheconvergencebehaviorofthe the Griewank function (Fig. 9b), which could be interpreted as a
algorithms, we plotted the average (out of 30 trials) value of the rough, multimodal version of sphere. A plausible explanation for
bestpointfoundsofaragainstthenumberoffunctionevaluations thisisthecombinedeffectoftheknowledgegainedbythemeta-
(Figs.3e8).Eachfigurereferstoaspecifictestfunctionandcom- model,whicheasilyrecognizesthesphericalstructureofGriewank,
prises four charts, for the alternative configurations (two and the simplex-based operators, using approximations of the
dimensions(cid:3)twoMFE). gradientofthefunction.
Inmostcases,SEEASexhibitsthefasterconvergence,evidently Aninterestingconclusionisthat,regardingSEEAS,theincrease
becausetheexpansionmechanismssupportedbythemetamodel ofthecomputationbudgethasmildeffectsintheimprovementof
(which provides enhanced overview of the surface geometry),
I.Tsoukalasetal./EnvironmentalModelling&Software77(2016)122e142 131
Fig.3. ConvergencecurvesfortestfunctionOF1(Sphere)with15(a,b)and30variables(c,d),withMFE¼500(a,c)andMFE¼1000(b,d).
themeanbestsolution.Thisisanotherevidenceofthesuitabilityof occasionally called in order to avoid making search too random.
SEEAS for extremely time-demanding optimization problems, in Finally,thesetupwithx¼2providessystematicallybetterresults
which the desirable number of function evaluations should be compared to x ¼ 1, while it exhibits either better or similar per-
minimal. formance when the annealing cooling parameter increases up to
x¼4(Table7).Nevertheless,acommonoutcomeformtheabove
4.4. SensitivityanalysisagainstinputparametersofSEEAS investigationsistherelativelylowsensitivityofSEEASagainstthe
examinedconfigurations,formostoftestproblems.
As mentioned in Section 2.2.4, SEEAS requires determining
severalinputarguments,intermsofstepparametersNr,Ne,Ncand 4.5. Suitabilityassessmentbasedonstochasticdominance
Nu,mutationprobabilitypm,andadjustingfactorxoftheannealing
cooling schedule. In order to investigate the sensitivity of SEEAS Foreachproblemandeachalgorithmweillustratetheempirical
against the default values adopted so far (i.e., CDFs,usingthesampleof30bestsolutionsfoundafterthetermi-
Nr ¼ Ne ¼ Nc ¼ Nu ¼ 20, pm ¼ 0.10 and x ¼ 2), weemployed 30 nationofthecorrespondingsearchprocedures.Basedonthem,we
independentrunsofthetestsfunctions(forn¼15variablesand calculatedthemediansof theCDFs(Table8),andnextemployed
MFE¼500),assigningdifferentvaluestoitsinputparameters.The theMWUtestbetweenthealgorithmsprovidingthebest(lower)
configurations and summary statistics, in terms of means and medians,toassesswhethertheobtaineddifferencesaresignificant.
standard deviation of the optimal solution of each set of optimi- TheresultsofalltestsaresummarizedinTable9.
zations,aregiveninTables5e7. The outcomes of the MWU test are in full accordance with
The analysis justifies our recommendations for the input pa- previousconclusions,andclearlyprovethestatisticalsuitabilityof
rameters of SEEAS. As shown in Table 5, the performance of the SEEAS.Consideringthefullsetofproblems,SEEASisevaluatedas
algorithm is significantly improved by increasing the common “preferred” or “equally good” in 18 out of 24 cases. Next best
value of the step parameters from 5 to 20, while it is slightly methodisDYCORS,whichispreferredorequallygoodin6outof24
improvedbyfurtherincreasingthisvalueto50.Actually,thesim- cases.Ifweisolate the lessbeneficialsubset,i.e., the formulation
plextransitionsareconsiderablyassistedbyusingtheoutcomesof with30decisionvariablesunderthelowercomputationalbudget
thesurrogatemodelwithinlocalsearch;however,itdoesnotmake (lowerleftpanelofTable9),thesuperiorityofSEEASisevenmore
sense calling the SM too many times, which introduces unnec- evident.
essary computations with marginally only benefit. Regarding the
mutation probability (Table 6), the algorithm provides almost 5. Hydrologicalcalibration
identical results for pmvalues as lowas 0.05 or 0.10, yet its per-
formanceisclearlydeterioratedbyincreasingthisprobabilityupto 5.1. Studyarea,simulationmodelandcalibrationsetup
0.30.Thisisalsoanon-surprisingconclusion,sinceitiswell-known
that in evolutionary algorithms the mutation operator should be Hydrological calibration is probably the most typical global