Dominance-Based Multi-Objective Simulated Annealing

20 Pages·2017·2.32 MB·English

IEEETRANSACTIONSONEVOLUTIONARYCOMPUTATION, DRAFT UNDERREVIEW 1
Dominance-Based
Multi-Objective Simulated Annealing
Kevin I. Smith, Richard M. Everson, Jonathan E. Fieldsend Member, IEEE,
Chris Murphy and Rashmi Misra
Abstract—Simulated annealing is a provably convergent opti- thereareseveralwelldevelopedgeneticalgorithmsandevolu-
miser for single-objective problems. Previously proposed multi- tionaryschemestoaddresssuchmulti-objectiveproblems(see,
objective extensions have mostly taken the form of a single-
[4], [5] for recent reviews), simulated annealing does not, in
objective simulated annealer optimising a composite function of
its usual formulation, provide a method for optimising more
the objectives. We propose a multi-objective simulated annealer
utilisingtherelativedominanceofasolutionasthesystemenergy than a single objective.Simulated annealing has been adapted
foroptimisation,eliminatingproblemsassociatedwithcomposite to multi-objective problems by combining the objectives into
objective functions. We also demonstrate a method for choosing a single objective function [6]–[10]; however, these methods
perturbationscalingspromotingsearch bothtowardsandacross
either damage the proof of convergence, or are limited (po-
the Pareto front.
tentially severely)in theirability to fullyexplorethe trade-off
We illustrate the simulated annealer’s performance on a
suite of standard test problems and provide comparisons with surface.
another multi-objective simulated annealer and the NSGA-II Weproposeamodifiedsimulatedannealingalgorithmwhich
genetic algorithm. The new simulated annealer is shown to maps the optimisation of multiple objectives to a single-
promote rapid convergence to the true Pareto front with a good objective optimisation using the true trade-off surface, main-
coverage of solutions across it comparing favourably with the
taining the convergence properties of the single-objective an-
other algorithms.
nealer while encouraging exploration of the full trade-offsur-
An application of the simulated annealer to an industrial
problem, the optimisation of a Code Division Multiple Access face. A method of practical implementation is also described,
(CDMA) mobile telecommunications network’s air interface, is usingtheavailablenon-dominateddatapointsfromthecurrent
presented and the simulated annealer is shown to generate optimisation to overcome the limitation that the true trade-off
non-dominated solutions with an even and dense coverage that
surface is unavailable for most real-world problems.
outperform single objective genetic algorithm optimisers.
In this paper, following some introductory material in sec-
Index Terms—Multiple objectives, simulated annealing, dom-
tion II, we start by briefly discussing methods that combine
inance, CDMA networks.
objectives into a single composite objective. In section III we
describe our dominance-based SA algorithm and, in section
I. INTRODUCTION IV, methods are described for improving the quality of the
Apopular androbustalgorithm for solving single-objective optimisation energy measure when the available data points
optimisation problems (those in which the user cares only are few. Choosing an efficient scale for perturbations is an
about a single dependant variable of the system) is simulated important component of scalar SA algorithms. The issue is
annealing (SA) [1], [2]. Geman & Geman [3] provided a further complicated in multi-objective algorithms because a
proofthatsimulatedannealing,ifannealedsufficientlyslowly, perturbation may not only move the current state closer to or
converges to the global optimum, and although the required furtherfromtheParetofront,butalsotransversally(i.e.,across
cooling rate is infeasibly slow for most purposes, simulated the front). In section V we describe a method for setting the
annealing often gives well convergedresults when run with a scale of perturbations and other run-time parameters. Results
fastercoolingschedule.Itisfrequentlythecaseinoptimisation showing that the algorithm converges on a range of standard
problems, however, that there are several objectives of the test problems are given in section VI, and we show that the
system which the user is interested in optimising simultane- algorithm compares favourably with both the popular NSGA-
ously.Clearly,simultaneousoptimisationofseveralobjectives IImulti-objectivegeneticalgorithm[31]andamulti-objective
is usually impossible and the curve (for two objectives) or simulated annealer suggested by Nam & Park [8]. In section
surface (for three or more objectives)that describes the trade- VII we present results demonstrating the simulated annealer’s
offbetweenobjectivesis knownas thePareto-front.Although performance on the optimisation of the air interface of a
CDMA network in the mobile telecommunications domain.
Kevin Smith is supported by the Engineering and Physical Sciences We drawconclusions in section VIII. A preliminary report on
Research Council (EPSRC) and Motorola. Jonathan Fieldsend is supported
this work appeared in [11].
bytheEPSRC,grantGR/R24357/01.
KIS, RME and JEF are with the Department of Computer
Science, University of Exeter, Exeter, EX4 4QF, UK. (e-mail: II. BACKGROUND
fK.I.Smith, R.M.Everson, [email protected]); CM and RM
A. Dominance and Pareto Optimality
are with Motorola, Thamesdown Drive, Swindon, SN25 4XY, UK. (e-mail:
fChris.Murphy,[email protected]) In multi-objective optimisation we attempt to simultane-
0000–0000/00$00.00ou(cid:13)csly200m5aIExEimEise or minimise D objectives, yi; which are
IEEETRANSACTIONSONEVOLUTIONARYCOMPUTATION, DRAFT UNDERREVIEW 2
functions of P variable parameters or decision variables, x= Algorithm 1 Simulated annealing
(x1;x2;:::;xP): Inputs:
fL gK Sequence of epoch durations
yi =fi(x); i=1;:::;D: (1) fTkkgKkk==11 Sequence of temperatures, Tk+1 <Tk
x Initial feasible solution
Withoutlossofgenerality,weassumethattheobjectivesareto
beminimised,sothatthemulti-objectiveoptimisationproblem
1: for k :=1;:::;K
may be expressed as:
2: for i:=1;:::;L
k
Minimise y=f(x)(cid:17)(f1(x);:::;fD(x)): (2) 3: x0 :=perturb(x)
4: (cid:14)E(x0;x):=E(x0)(cid:0)E(x)
The idea of dominance is generally used to compare two 5: u:=rand(0;1)
solutions f and g: If f is no worse for all objectives than
6: if u<min(1;exp((cid:0)(cid:14)E(x0;x)=T ))
k
g and wholly better for at least one objective it is said that f
7: x:=x0
dominates g; written f (cid:30)g. Thus f (cid:30)g iff:
8: end
f (cid:20)g 8i=1;:::;D and 9: end
i i (3)
f <g for at least one i: 10: end
i i
By a slight abuse of notation, dominance in objective space
is extended to parameter space; thus it is said that a(cid:30)b iff
sufficientlyslow annealing schedule is used, any sample from
f(a)(cid:30)f(b).
(cid:25) will almost surely lie at the minimum of E.
T
Clearly the dominates relation is not a total order and two
Sampling from the equilibrium distribution (cid:25) (x) at any
T
solutionsaremutuallynon-dominatingifneitherdominatesthe
particular temperature is usually achieved by Metropolis-
other. A set F of solutions is said to be a non-dominating set
Hastings sampling [2], which involves making proposals x0
if no element of the set dominates any other:
that are accepted with probability
a6(cid:30)b 8 a;b2F (4) A=min(1;expf(cid:0)(cid:14)E(x0;x)=Tg) (5)
A solution is said to be globally non-dominated, or Pareto- where
optimal, if no other feasible solution dominates it. The set
(cid:14)E(x0;x)(cid:17)E(x0)(cid:0)E(x): (6)
of all Pareto-optimal solutionsis knownas the Pareto-optimal
front,ortheParetoset,P;solutionsintheParetosetrepresent Intuitively, when T is high perturbations from x to x0 which
the possible optimal trade-offsbetween competing objectives. increase the energy are likely to be accepted (in addition
A human operator can select a solution with a knowledge to perturbations which decrease the energy, which are al-
of the trade-offs involved once this set has been revealed. ways accepted) and the samples can explore the state space.
Heuristic procedures, such as multiple objective evolutionary Subsequently, as T is reduced, only perturbations leading
algorithms and the multi-objective simulated annealing algo- to small increases in E are accepted, so that only limited
rithms discussed here, yield sets of mutually non-dominating explorationispossibleasthesystemsettleson(hopefully)the
solutions which will be only an approximation to the true global minimum. The algorithm is summarised in Algorithm
Paretofront.Somecarewithterminologyisthereforerequired, 1: during each of K epochs, the computational temperature is
and in this paper the set produced by such an algorithm is fixed at T and L samples are drawn from (cid:25) before the
k k Tk
referred to as the estimated Pareto front, which we denote by temperature is lowered in the next epoch. Each sample is a
F. perturbation (‘mutation’ in the nomenclature of evolutionary
algorithms) of the current state from a proposal density (line
3); the perturbed state x0 is accepted with probability given
B. Simulated Annealing
by (5), as shown in lines 4-8.
Simulated annealing, introduced by Kirkpatrick et al. [1]
As already alluded to, convergence is guaranteed if and
may be thought of as the computational analogue of slowly
only if the cooling schedule is sufficiently gradual [3], but
cooling a metal so that it adopts a low-energy, crystalline
experience has shown SA to be a very effective optimisation
state. At high temperatures particles are free to move around,
technique even with relatively rapid cooling schedules [12],
whereas as the temperature is lowered they are increasingly
[13].
confined due to the high energy cost of movement. It is
physically appealing to call the function to be minimised the
C. Multi-Objective SA with Composite Objective Functions
energy, E(x); of the state x; and to introduce a parameter T;
An attractive approach to multi-objectivesimulated anneal-
the computational temperature, which is lowered throughout
ing (MOSA), adopted by several investigators [7]–[10], [14]–
the simulation according to an annealing schedule. At each T
[17], is to combine the objectives as a weighted sum:
the SA algorithm aims to draw samples from the equilibrium
distribution (cid:25)T(x) / expf(cid:0)E(x)=Tg. As T ! 0 more and D
more of the probability mass of (cid:25)T, is concentrated in the E(x)=Xwifi(x): (7)
regionoftheglobalminimumofE;soeventually,assuminga i=1
IEEETRANSACTIONSONEVOLUTIONARYCOMPUTATION, DRAFT UNDERREVIEW 3
The composite objective is then used as the energy to be change probabilities. In this algorithm the relative dominance
minimised in a scalar SA optimiser. An equivalent alternative of the current and proposed solutions is tested and when the
[6] is to sum logf (x); and others (e.g., [8], [15]) have proposed solution dominates the current solution the proposal
i
investigated a number of non-linear and stochastic composite is accepted; this is analogous to the automatic acceptance
energies. of proposals with a lower state energy in single-objective
Itisclearthatsimulatedannealingwitha compositeenergy simulated annealing. In addition to the common practise
(7) will converge to points on the Pareto optimal front where of employing a state change probability which guarantees
the objectives have ratios given by w(cid:0)1, if such points exist. acceptance of strictly superior perturbations, Nam & Park
i
However, it is unclear how to choose the weights in advance, modify the acceptance rule so that proposals are accepted
indeed, one of the principal advantages of multi-objective with probability given by (5) and (7) if they are dominated
optimisation is that the relative importance of the objectives by x; but unconditionally accepted if x0 (cid:30) x or if x0 and
can be decided with the Pareto front on hand. Perhaps more x are mutually non-dominating. This promotes exploration
importantly, parts of the front are inaccessible with fixed of the search space and escape from local fronts but as the
weights [18]. Recognising this, investigators have proposed dimensionality increases so does the proportion of all moves
a variety of schemes for adapting the w during the annealing which are accepted unconditionally.This limits the behaviour
i
process to encourage exploration along the front. See for of the algorithm to that of a random walk through the search
example [19]. space when dealing with problems with high dimensionality.
It isnatural to keepan archive,F; of allthenon-dominated When the proposed solution is dominated by the current
solutions found so far, and this archive may be utilised to solution,Nam& Parkdefinetheenergydifferencecontrolling
furtherexplorationbyperiodicallyrestartingtheannealerfrom acceptance as the average differencein objectivevalues. Nam
a randomly chosen element of F [10]. & Park also employ 100 separate agents during optimisation,
A proposal x0 in scalar SA is either better or worse than whereeachagentisanindependentcopyofthealgorithm;this
thecurrentstatexdependingonthesignof(cid:14)E(x0;x);except serves a similar function to Suppapitnarm etal.’s return-to-
for pathological problems the probability that (cid:14)E = 0 is base approach to promoting diversity of the solutions located
vanishingly small. In multi-objective SA, however, x0 may by the algorithm.
dominate x or x0 may be dominated by x or they may Although it is clear that the assurance of a convergence
be mutually non-dominating: in fact, the probability that a proofcanbeprovidedforamulti-objectivesimulatedannealer
pair of randomly chosen points in D-dimensional space are using a scalar objective function and fixed weights (7), such
mutually non-dominating is 1(cid:0)2 1 D, so the mutually non- annealers are fundamentally limited in their coverage of the
(cid:0)2(cid:1)
dominating case becomes increasingly common with more Pareto front. On the other hand, it is difficult to see how
objectives.However,energiessuchas(7)mayleadtox0 being proofs of convergence might be obtained with the heuristic
accepted unconditionally ((cid:14)E < 0) even though x0 6(cid:30) x; modifications designed to promote exploration transversal to
because a large negative energy change from one objective the front. Given these difficulties, defining a multi-objective
may outweigh small positive changes on the other objectives. simulated annealer which utilises a composite objective func-
Each multi-objective simulated annealing algorithm which tion is undesirable. With this in mind, we investigate the
utilisesacompositeobjectivefunctionmustthereforedealwith efficacyof an energyfunctionbased on the defining notion of
this behaviour in some manner. dominance. The aim is the definition a single energyfunction
A good example of a composite objective function ap- appropriate to all cases of relative dominance between x and
proach to multi-objective simulated annealing is given by x0 withoutrequiringspecial cases for wherex0 (cid:30)x, orwhere
Suppapitnarm etal.[10]. Instead of weighting and summing x0 and x are mutually non-dominating, as has been the case
the objectives to produce a composite energy difference for in previous algorithms.
the acceptance criteria, this algorithm uses a multiplicative
function with individual temperatures for each objective each
III. ADOMINANCE BASEDENERGYFUNCTION
of which is adjusted independently by the algorithm. This
negates the need for a priori weighting of the objectives, In single objective optimisation problems the sign of the
and can be considered to function as a weighted compos- difference in energy (cid:14)E(x0;x) tells us whether the proposal
ite sum approach with algorithmically controlled weightings. x0 is a better, worse or (very rarely) equally good solution
This is vulnerable to the concentrated search properties of as the current solution x. Likewise the dominance relation
other composite objective techniques and Suppapitnarm etal. can be used to compare the relative merit of x0 and x in
employ a return-to-base scheme whereby the current solution multi-objectiveproblems,butnotethatitgivesessentiallyonly
is re-seeded with another solution from the non-dominated three values of quality—better, worse, equal—in contrast to
archive to promote a more even coverage. Suppapitnarm et the energydifference in uni-objectiveproblemswhich usually
al. promote exploration along the front by unconditionally gives a continuum.
accepting proposals that are not dominated by any member If the true Pareto front P were available, we could define a
of F; otherwise using (5). simple energyof x asthe measure ofthe front that dominates
Of the multi-objective simulated annealing techniques in f(x). Let Px be the portion of P that dominates f(x):
the current literature, perhaps the most promising is that of
Nam & Park [8] due to their use of dominance in state Px =fy2Pjy(cid:30)f(x)g: (8)
IEEETRANSACTIONSONEVOLUTIONARYCOMPUTATION, DRAFT UNDERREVIEW 4
in this case. Since all solutions lying on the front have
P
equal minimum energy, we would anticipate that a simulated
2 annealer using this energy would, having reached the front,
e
v perform a random walk exploration of the front.
cti
e We note that Fleischer [20] has proposed an alternative
bj
O measure of a non-dominated set, which may be loosely char-
acterised as being based on the volume dominated by the set
rather than the area of the dominating set.
Unfortunately,the true Pareto front P is unavailableduring
the course of an optimisation. We therefore propose to use an
energy defined in terms of the current estimate of the Pareto
front,F,whichisthesetofmutuallynon-dominatingsolutions
foundthus far in the annealing.Denote by F~ the union of the
F, the current solution x and the proposed solution x0, that is
Objective 1
F~ =F [fxg[fx0g: (10)
Fig.1. EnergyfromareaofthetrueParetofrontP dominatingasolution. Then, in a similar manner to (8), let F~x be the elements of F~
SolutionsaremarkedbycirclesandlinesindicatetheregionsofPdominating
that dominate x:
eachone.
F~x =fy2F~jy(cid:30)xg (11)
P so that we obtain an energy difference between the current
2 and proposed solutions of
e
v 1
ecti (cid:14)E(x0;x)= jF~j(cid:16)jF~x0j(cid:0)jF~xj(cid:17): (12)
bj
O
Division by jF~j ensures that (cid:14)E is always less than unity and
F
provides some robustness against fluctuations in the number
of solutions in F. If F~ is a non-dominating set the energy
difference between any two of its elements is zero. Note also
that (cid:14)E(x0;x) = (cid:0)(cid:14)E(x;x0): The inclusion of the current
solution and the proposal in F~ means that (cid:14)E(x0;x) < 0 if
x0 (cid:30)x; which ensures that proposals that move the estimated
fronttowardsthetruefrontarealwaysaccepted.Proposalsthat
Objective 1 aredominatedbyoneormoremembersofthecurrentarchive
are acceptedwith a probability dependingupon the difference
in the number of solutions in the archive that dominate x0
Fig.2. EnergyfromproportionoftheestimatedParetofrontF dominating
pointsdominatingasolution.ElementsofF areshownassmallgreycircles, and x. We emphasise that this probability does not depend
solutionsareshownaslargeropenorfilled circles. uponmetricinformationinobjectivespace;weputnoapriori
weighting on the objectives and the acceptance probability is
unaffected by rescalings of the objectives.
Then we define
Afurtherbenefitofthisenergymeasureisthatitencourages
explorationofsparselypopulatedregionsofthefront.Imagine
E(x)=(cid:22)(Px) (9)
two proposals, each dominated by some solutions in F; for
where (cid:22) is a measure defined on P. We shall be principally example, the solutions illustrated by the filled and unfilled
interested in finite sets approximating P and so shall take circles in Figure 2. The solution that is dominated by fewer
(cid:22)(Px)tobesimplythecardinalityofPx.IfP isacontinuous elements (the unfilled circle) has the lower energy and would
set, we can take (cid:22) to be the Lebesgue measure (informally, therefore be more likely to be accepted as a proposal.
thelength,areaorvolumefor2,3or4objectives);wefurther Defining the energy in this manner, unlike some pro-
discussmeasuresinducedonP insectionVI-E.Asillustrated posed multi-objective enhancements to simulated annealing
in Figure1, thisenergyhas thepropertieswe desire:if x2P discussedinsectionII-C,providesasingleenergyfunctionen-
then E(x)=0, and solutions more distant from the front are, couragingbothconvergencetoandcoverageoftheParetofront
in general, dominated by a greater proportion of P and so without requiring other modifications to the single-objective
have a higher energy; in Figure 1 the solution marked by an simulated annealing algorithm (beyond the obvious storage
open circle has a greater energy than the one marked by a of an archive of the estimated Pareto front). In particular no
filled circle. additionalrulesarerequiredforcasesinwhichthecurrentand
Clearly this formulation of an energy does not rely on proposed solutions are mutually non-dominating.
an a priori weighting of the objectives and the assurances Convergence to the true Pareto front is no longer an
of convergence [3] for uni-objective SA continue to hold immediate consequence of Geman & Geman’s work [3],
IEEETRANSACTIONSONEVOLUTIONARYCOMPUTATION, DRAFT UNDERREVIEW 5
because the energy based on F is only an approximation
to (9). However, Greening [21] offers proof of convergence,
2
albeit more slowly, even when the energy contains errors. e U
v
Current work is investigating the application of this work to cti
MOSA and in section VI we offer empirical evidence of the e
bj
convergence. O
An energy function based on (12) is straightforward to
calculate;countingthenumberofelementsofF~ thatdominate
x and x0 can be achieved in logarithmic time [22], [23].
S
Our proposed multi-objective algorithm closely follows the
standardSAalgorithm(Algorithm1),theonlyadditionthatis H
necessary is to maintain an archive, F of the current estimate
of the Pareto front and to calculate the energy difference
using (12). However, we postpone detailed description of the
Objective 1
algorithm until methods of increasing the empirical energy
resolution have been discussed.
Fig.3. AttainmentsurfaceSF istheboundaryoftheregion,U dominated
bythenon-dominatedsetF;whoseelementsaremarkedasdots.Dashedlines
IV. INCREASING ENERGYRESOLUTION denoteH theminimumrectanglecontainingF.
As mentioned earlier, the true Pareto-optimal front of so-
lutions is, in general, unavailable to us. While using the
proposebelowinsiststhatinterpolatingpointsareonlyweakly
archive of the estimated Pareto front F provides an estimate
dominated by the archive.
of solution energy, when F is small the resolution in the
energies can be very coarse. In fact, the difference in energy
betweentwosolutionsisanintegermultipleof1=jF~jbetween B. Linear Interpolation
0 and 1. Since the acceptance criterion (5) for new solutions Another apparently suitable method of augmenting F is
is determined by the difference in energy (cid:14)E(x;x0) between linear interpolation (in objective space) between the solutions
the current solution and the proposed solution, low resolu- in F. In this method, when the archive is smaller than some
tion of the energies leads to a low resolution in acceptance predefined size, new points in objective space are generated
probabilities. At low computational temperatures and with on the simplices defined by an element of F and its D (cid:0)
small archives it will become increasingly likely that this 1 nearest neighbours in F. This overcomes the limitations
granularity will make it almost impossible for even slightly of the previous method: Since new solutions are generated
detrimental moves (i.e., moves that increase E(x)) to be ‘on’ the current estimated Pareto front, the problems which
made. This is undesirable as, at its most severe, this effect could occur with using old, dominated elements of F in
reduces the algorithm to behaviour similar to a greedy search the energy calculations are avoided. The interpolated points
optimiser,andpreventstheexploratorybehaviourprovidedby generated can also be evenly distributed between the current
detrimental moves. estimated Pareto-optimal solutions, which is beneficial as it
This problem is alleviated by using a larger set for energy does not deter the algorithm from exploring any region of the
calculations. There are a couple of straightforward, but ulti- estimated front which is not already densely populated. The
mately inadequate, methods for artificially increasing the size principaldisadvantageofthismethodisthatproposalsmaybe
ofF whichwenowbrieflydiscussbeforedescribingamethod dominatedby an interpolated point, but not by anyof the real
using the attainment surface. elements of F; meaning that the proposal may erroneously be
disregarded.
A. Conditional Removal of Dominated Points
C. Attainment Surface Sampling
A straightforward method for increasing the size of the
archive is not to delete solutions known to be dominated Consideration of the previous two methods of augmenting
if deleting them would reduce jFj below some predefined theestimatedParetofrontsuggeststhattheaugmentingpoints
minimum. However, the existence of old solutions in F, may should have the following properties.
lead to desirable proposals (i.e., not non-dominated solutions) 1) The augmenting points must be sufficiently close to the
being rejected. In addition the old solutions may bias the currentestimationoftheParetofrontthattheycanaffect
search away from regions of the front that were previously the energy of solutions generated near to the current
well populated. estimated Pareto front.
A further disadvantage of this method is that the retained 2) Theymustbeevenlydistributedacrossthecurrentlyesti-
solutionsmaybepositionedsothat theyaredominatedbythe matedParetofrontsoastonotdiscouragethealgorithm
archive and indeed by the current point and the vast majority fromacceptingproposalsinpoorlypopulatedregionsof
of proposals. In this case they serve to increase the resolution the front.
of the energy at the expense of the range. By contrast the 3) Theymustnotdominateanyproposalwhichisnotdom-
interpolation method using the attainment surface that we inatedbyanymemberofF; so thatpotential entrantsto
IEEETRANSACTIONSONEVOLUTIONARYCOMPUTATION, DRAFT UNDERREVIEW 6
Algorithm 2 Sampling a point from the attainment surface
Inputs:
fL gD Elements of F; sorted by
d d=1
increasing coordinate d 1
0.9
1: for i:=1;:::;D Generate a random point, v 0.8
2: v :=rand(min(L );max(L )) 0.7
34:: edn:=dirandint(1;D) i Ci hoose a dimension, d Objective 300..56
0.4
5: for i:=1;:::;jFj Find smallest v s.t. v is
d 0.3 0
6: u:=L dominated by an element of F
d;i
0.2
7: vd :=ud 0.1 0.5
89:: ifrFet(cid:30)urvn v 0 0.2 0.4 0.6 0.8 1 1 Objective 1
Objective 2
10: end
11: end Fig. 4. 10000 samples from the attainment surface for an archive of 10
points,whicharemarkedwithheavydots.
the archive are not discarded. A consequence of this is
that they must all be dominated by at least one member points; 10000 samples are shown for visualisation. In the
of F. experiments reported in section VI F was augmented with
100 samples from S before calculating the energy of the
The attainment surface, which has previously been used F
proposal.
for estimated Pareto front visualisation [24] and is closely
related to the attainment function [25], is an interpolating
surface between the elements of F that has the requisite D. Multi-objective simulated annealing algorithm
properties. The attainment surface, S ; corresponding to F
F Having discussed sampling from the attainment surface to
is a conservative interpolation of the elements of F so that
increase the energy resolution, we are now in a position to
every point of S is dominated by an element of F. The
F summarise the main points of our proposed multi-objective
attainment surface for an F comprising four two-dimensional
simulated annealing algorithm. As shown in Algorithm 3,
elementsissketchedinFigure3.Moreformally,theattainment
the multi-objective algorithm differs from the uni-objective
surfaceistheboundaryoftheregioninobjectivespacewhich
algorithm in that an archive F of non-dominated solutions
is dominated by elements of F. If u;v 2 RD; we say that
foundso far is maintained,and the energydifferencebetween
u properly dominates v (denoted uCv) iff u < v 8i =
i i the proposed and current solution is calculated based on the
1;:::;D. Then if
current archive or its attainment surface.
F =fyju(cid:30)y for some u2Fg (13) The archive is initialised with the initial feasible point
(line 1 of Algorithm 3). At each stage the current solution
U =fyjuCy for some u2Fg (14)
x is perturbed to form the proposed solution x0. In the work
the attainment surface is S =F nU =@U. reported here, in which the parameters x are continuous and
F
We interpolate F with random samples uniformly dis- real valued, we perturb each element of x singly, drawing
tributed on S \H , the attainment surface restricted to the the perturbations from a Laplacian distribution centred on the
F F
minimumaxis-parallel hyper-rectanglecontainingF (see Fig- current value.
ure3).FromthedefinitionofS itisapparentthatinterpolated IftherearesufficientlymanysolutionsinF,theaugmented
F
points are dominated by an element of F; thus satisfying the archiveF~ isconstructedbyaddingxandx0 (line 9)toF and
third criterion. Uniform random sampling ensures that the the energy difference between x0 and x is calculated using
second criterion is met, as is the first criterion because S (12). If there are fewer than S solutions in then additional
F
interpolates F. samples are drawn from the attainment surface S using
F
Sampling from S may be performed using Algorithm 2, Algorithm 2 (line 6); the energy difference is then calculated
F
which works by sampling a point from a uniform distribution based on the sampled attainment surface, x and x0. In the
on the surface of H and then restricting one coordinate work reported here, we always augment F with 100 samples
F
so that the point is dominated by an element of F. This is from S . As even when there are a large number of solutions
F
facilitatedbytheuseoflistsL ,d=1;:::;Dwhichcomprise in the archive of the estimated Pareto front it is worthwhile
d
the elements of F sorted in increasing order of coordinate d. sampling from S as this samples evenly across the front,
F
Determining whether an element of F dominates v on line 8 providinggreater resolutionin sparselypopulated areas ofthe
may be efficiently implemented using binary searches of the front.
listsL ;in whichcasethealgorithmrequiresO(jFjlog(jFj)) If the proposal is accepted (line 14), the archive must be
d
time for the generation of each sample. Figure 4 illustrates updated.Ifxisnotdominatedbyanyofthearchivalsolutions,
the sampled attainment surface for a set of ten 3-dimensional allarchivalsolutionsthataredominatedbyxaredeletedfrom
IEEETRANSACTIONSONEVOLUTIONARYCOMPUTATION, DRAFT UNDERREVIEW 7
Algorithm 3 Multi-objective simulated annealing by using a short ‘burn-in’ period during which all solutions
Inputs: are accepted and setting the temperature equal to the average
fL gK Sequence of epoch durations positive change of energy divided by ln(2).
k k=1
fTkgKk=1 Sequence of temperatures, Tk+1 <Tk IntheworkreportedhereallepochsLk areofequallength,
x Initial feasible solution L = 100 and we adjust the temperature according to T =
k k
(cid:12)kT0, where (cid:12) is chosen so that Tk is 10(cid:0)5 after two thirds
1: F :=fxg Initialise archive of the evaluations are completed .
2: for k :=1;:::;K
3: for i:=1;:::;L
k B. Perturbation Scalings
4: x0 :=perturb(x)
5: if jFj<S If F is small Forsimplicitya proposalisgeneratedfromxbyperturbing
Construct attainment surface only one parameter or decision variable of x. The parameter
6: S :=interpolate(F) to be perturbed is chosen at random and perturbed with
F
7: F~ :=S [F [fxg[fx0g a random variable (cid:15) drawn from a Laplacian distribution,
F
8: else p((cid:15)) / e(cid:0)j(cid:27)(cid:15)j; where the scaling factor (cid:27) sets magnitude of
9: F~ :=F [fxg[fx0g theperturbation.TheLaplaciandistributionhastailsthatdecay
relativelyslowly,thus ensuringthat there isa high probability
10: end
Energy difference base on F~ of exploring regions distant from the current solutions.
We maintain two sets of scaling factors, since the per-
11: (cid:14)E(x0;x):=E(x0)(cid:0)E(x)
12: u:=rand(0;1) turbations generating moves to a non-dominated proposal
within a front (we call these traversals) may potentially be
13: if u<min(1;exp((cid:0)(cid:14)E(x0;x)=T ))
k
very different from those required to locate a front closer
14: x:=x0 Accept new current point
to P, which we call location moves. We maintain a scaling
If x is not dominated by any element of F
factor for each dimension of parameter space for each of
15: if z6(cid:30)x 8z2F
the location perturbations and the traversal perturbations, and
Remove dominated points from F
adjust these independently to increase the probability of such
16: F :=fz2Fjx(cid:6)zg
a move being generated. When perturbing a solution, it is
17: F :=F [x Add x to F
chosen randomly with equal probability whether the location
18: end
scaling set will be used, or the traversal scaling set. Statistics
19: end
are kept on perturbations generating traversal and location
20: end
moves; clearly these can be updated only after the proposal
21: end
has been generated so that the type of move is known. The
scalings are adjusted throughout the optimisation, whenevera
suitably large statistic set is available to reliably calculate an
the archive (line 16) and x is added to the archive (line 17).
appropriatescalingfactor.Thesescalingsareinitiallysetlarge
Clearly F is always a non-dominated set, although note that
enough to sample from the entire feasible space.
x0 may be dominated by members of F.
1) Traversal Scaling: The traversal rescaling for a particu-
lar decisionvariable x is performedwheneverapproximately
j
V. REALTIME ALGORITHMPARAMETEROPTIMISATION 50 traversal perturbations have been made to x since the last
j
The performance of this algorithm, in common with other rescaling.
simulatedannealingsystems,dependsuponparametersforthe In order to ensure wide coverage of the front we wish to
initial temperature, the annealing schedule and the size of maximise the distance (in objective space) covered by the
perturbations made to solutions when generating new propos- traversals to ensure the entire front is evenly covered. Gen-
als. Here we give details of methods which permit automatic erating traversals travelling a small distance will concentrate
setting of the initial temperature, and which adjust the scale theestimatedfrontaroundthe pointatwhichthe currentfront
of perturbations made to maximise the quality of proposed was discovered, an effect we aim to avoid.
solutions. We seek to generate proposals on approximately the scale
thathaspreviouslybeensuccessfulingeneratingwide-ranging
traversals. To achieve this, the perturbations are sorted by
A. Annealing Schedule
absolute size of perturbation in parameter space, and then
If the initial computational temperature is set too high, trisected in order, giving three groups, one of the smallest
all proposed solutions will be accepted, irrespective of their third of perturbations, the largest third of perturbations, and
relative energies, and if set too low proposals with a higher theremainingperturbations.Foreachgroupthemeantraversal
energythan the current solution will neverbe accepted, trans- size caused by the perturbations is calculated. The traversal
forming the algorithm into a greedy search. As a reasonable size is measured as the Euclidean distance travelled in objec-
starting point we set the initial temperature to achieve an tivespacewhenthecurrentsolutionandtheproposedsolution
initial acceptance rate of approximately 50% on derogatory aremutuallynon-dominating.Ifaperturbationandthecurrent
proposals.Thisinitialtemperature,T0;canbeeasilycalculated solution arenot mutuallynon-dominating,the traversalsize is
IEEETRANSACTIONSONEVOLUTIONARYCOMPUTATION, DRAFT UNDERREVIEW 8
TABLEI
made and accepted at all temperatures.
TESTPROBLEMDEFINITIONOFDTLZ1–DTLZ7OF[26]FOR3
The location perturbation scaling is recalculated for each
OBJECTIVES(USINGTHESUGGESTEDPARAMETERSIZES).(DEFINITION
parameter for which 20 proposals having (cid:14)E(x0;x) > 0
OFDTLZ5CORRECTED.)
have been generated, after which the count is reset. Location
perturbation rescaling is omitted in two cases: First, when the
Problem Definition
archive of the estimated Pareto front F has fewer than 10
DTLZ1 fffg123((((xxxx))))====122211100xx((111jxx((cid:0)12j(cid:0)x((cid:0)112x)+(21g)(+(1xg)+)(gx)()x)) bbmyyemtthhbeeerstsae.mmSppelecerosantfdurlroyem,dwothheeesnantthotaetincemoxmceenbetidnseu1dr.fsaTiczheeiswofihsFenbaemucgauumltsieepnltiweedde
0(cid:20)+xPi(cid:20)Pi=13;(fxoir(cid:0)i=0:51);22;(cid:0)::c:o;sP(;2P0(cid:25)=(x7i(cid:0)0:5))) daderjuosgtattohreysmcaolivnegssatpoparottxeimmpattetloykaeethpirtdh;ewachceenpttahniscevaraluteeoisf
f1(x)=cos(x1(cid:25)=2)cos(x2(cid:25)=2)(1+g(x))
f2(x)=cos(x1(cid:25)=2)sin(x2(cid:25)=2)(1+g(x)) too small, it becomes impossible to generate such a scaling,
DTLZ2 f3(x)=sin(x1(cid:25)=2)(1+g(x)) and so the scalings are kept at the most recent valid value.
fg01((cid:20)(xx)x)=i=(cid:20)Pc1oPi;s=(f3oxr(1xi(cid:25)i==(cid:0)2)10;c:o25s;):2(:x:2;(cid:25)P=;2P)(=1+12g(x)) parCtiocuunlatirngdimonenlysiomnoovfespagreanmeerateterdspfarcoem, thpeertaucrcbeapttiaonncsetroatea
f2(x)=cos(x1(cid:25)=2)sin(x2(cid:25)=2)(1+g(x)) of derogatory moves (cid:11) is the fraction of proposals to a
DTLZ3 f3(x)=sin(x1(cid:25)=2)(1+g(x)) greater energy which are accepted. If (cid:27) denotes the location
g(x)=100(jxj(cid:0)2
perturbation scaling for a particular dimension, the new (cid:27) is
0(cid:20)+xPi(cid:20)Pi=13;(fxoir(cid:0)i=0:51);22;(cid:0)::c:o;sP(;2P0(cid:25)=(x1i2(cid:0)0:5))) set as:
DTLZ4 f0gffff31122((cid:20)(((((xxxxxx)x)))))=i=====(cid:20)Pscccci1oooonPi;ssss=(cid:0)(((cid:0)(cid:0)f3ox(cid:18)(cid:18)xxr11((cid:11)1(cid:11)1(cid:11)1x(cid:25)(cid:25)i(cid:25)(cid:25)(cid:25)i======(cid:0)22222))1(cid:1)(cid:1)(cid:1)0cs;(csi:o21n5ios;n)s+(:(2(cid:0)(cid:18)(cid:0)(cid:18):x2xg2:)(cid:11)2;)(cid:11)2(P(((cid:25)x(cid:25)11;=)=+)+2P2(cid:1)(cid:1)gg=((1(1(xx1++)2)))gg((xx)))) Thi(cid:27)s u:=pd8><>:ate(cid:27)(cid:27)(cid:27)(=w1(o1+r+k2s(2b(cid:11)(e0(cid:0)c:3au0(cid:0)s:4e(cid:11),)=)i=0n0:6g:3)e)neraiiifffl,(cid:11)(cid:11)0s:m3<>a(cid:20)00ll::e(cid:11)34r(cid:20)pe0rt:u4rbati(o1n5s)
f3(x)=sin((cid:18)1(cid:25)=2)(1+g(x)) in parameter space are more likely to generate small changes
DTLZ5 g(cid:18)1(x=)x=1PPi=3(xi(cid:0)0:5)2 in objective space, resulting in smaller changes in energy.
(cid:18)2= 4(1+(cid:25)g(x))(1+2g(x)x2)
0(cid:20)xi(cid:20)1;fori=1;2;:::;P;P =12 VI. EXPERIMENTS
ff21((xx))==xx21 Weillustratetheperformanceofthisannealeronsomewell-
DTLZ6 f3(x)=(1+g(x))h(f1;f2;f3;g) known test functions from the literature, namely the DTLZ
g(x)= P(cid:0)92PPi=3xi test functions of Deb et al. [26], [30], and compare them to
h(f1;f2;f3;g)=3(cid:0)P3i=1(cid:16)1f+ig(1+sin(3(cid:25)fi))(cid:17) theperformanceofthewellestablishedNSGA-IIevolutionary
0(cid:20)xi(cid:20)1;fori=1;2;:::;P;P =22 algorithm[31](usingthePISAreferenceimplementation[32])
Ds.tT.LZ7 fgff3121((((xxxx))))====f1111113000(PPPx)1i2i3i+===000112411xfxx1iii(x)(cid:0)1(cid:21)0 wtaenhsdticfhNunawcmteiodn&issciuPssatsrhkia’ntstshemecuttilrotuin-eoIbP.jaTerchetteiovbefernosenimfit,tuPolaf;teiusdsiknangnonwthenea,lDesroTLw[8Ze]
s.t. g2(x)=f3(x)+4f2(x)(cid:0)1(cid:21)0 candiscoverhowcloseourestimatedarchiveF istoP;aswell
s.t. g3(x)=2f3(x)+f1(x)+f2(x)(cid:0)1(cid:21)0 as compare results from each algorithm. Note that we rectify
0(cid:20)xi(cid:20)1;fori=1;2;:::;P;P =30 a couple of minor typographical errors in the description
of DTLZ5 and DTLZ6 here, as the formulae published in
[26], [30] do not yield the Pareto fronts described.1 For
counted as being 0. The traversal perturbation scaling for this completeness we give the problem definitions in Table I; in
decision variableis then set to the averageperturbation of the all the experiments we use D=3 objectives.
group which generated the largest average traversal. Inthe work reportedhereall epochs L areof equallength
k
This heuristic is open to the criticism that it depends for the annealers, L = 100 and we adjust the temperature
k
uponmeasuringdistancesinobjectivespacewhiletherelative according to Tk = (cid:12)kT0, where (cid:12) is chosen so that Tk is
weighting ofthe D objectivefunctions is unknown.To allevi- 10(cid:0)5 after approximately two thirds of the evaluations are
atethisdifficulty,however,theobjectivesmayberenormalised completed; run lengths and the exact number of evaluations
during optimisation so the front has approximately the same before T is 10(cid:0)5 are given in table II. For MOSA, the
k
extent in each objective. We emphasise that, of course, the parameter perturbations are controlled using the scheme de-
use of metric information for setting the approximate scale of scribed in section V-B. The perturbations for Nam & Park’s
perturbations does not affect the dominance-based energy. annealer are performed using a scheme similar to that for
2) Location Scaling: Drawing from methods widely used MOSA but without the automatic rescaling feature novel to
in evolutionary algorithms (see [27]–[29] for recent work in
this area), we aim to adjust the scale of location perturbations 1In equation (25) of [26] only (cid:18)1 should be multiplied by (cid:25)=2 when
to keep the acceptance rate for x0 that have a higher energy calculating f1;:::;fM. In equation (27) the calculation of g(xM) is
inconsistent with the results provided, meaning all f3 values in the figure
thanxtoapproximately1=3,sothatexploratoryproposalsare in[30]arehalved.
IEEETRANSACTIONSONEVOLUTIONARYCOMPUTATION, DRAFT UNDERREVIEW 9
MOSA − DTLZ1
Nam & Park − DTLZ1 NSGA−II − DTLZ1
1
300 300
200 200
0.5
100 100
0 0 0
0 1 0 200 0 300
100
0.5 0.5 200 100 200 100 200
1 0 300 0 400 0
6 5
150
5
4
4
100 3
3
2
50 2
1
1
0 0 0
0 0.01 0.02 0.03 0.04 0 20 40 60 80 100 120 140 0 50 100 150 200
Distance from the true Pareto front Distance from the true Pareto front Distance from the true Pareto front
Fig.5. Top:ArchivesontestproblemDTLZ1after5000functionevaluationsforeachofthethreealgorithms.Bottom:Histogramsofthedistancesfrom
thetrueParetofrontofthearchive members(the5%mostdistanthavebeenomittedtoaidvisualisationinallsixfigures).
102 7000
101 6000
100 5000
ce10−1 size4000
Distan10−2 chive 3000
Ar
10−3 2000
10−4 1000
10−5 0
0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5 3
Iteration x 104 Iteration x 104
Fig. 6. Left: Distance of current point,x; andarchive F from the true Pareto front, P; versus iteration for DTLZ1. The dottedline shows median over
20runsofdistanceofxfromP;dashedlinesshowmaximumandminimum(overthe20runs)distancesateachiteration.Thethicklineshowsthemedian
(over20runs)ofthemediandistanceofarchivememberstoP.Right:Archivegrowthversusiteration.Thicklineshowsmedian(over20runs)archivesize
anddashedlinesshowmaximumandminimum.
MOSA;thescalingsare fixedat0:1(determinedfrom asmall each of the DTLZ test problems, after which we present
empirical study). The parameters for the NSGA-II algorithm statistical results summarising the performance over 20 runs.
used were those suggested as the default values in the PISA We use the non-parametric Mann-Whitney rank-sum test (at
[32]package2.Weuse100simultaneouschainsfortheNam& the 0.05 level) to test for significant differences between
Park implementation and a population of size 100 for NSGA- the algorithms in the hypervolume and true front distance
II. comparison measures. Files containing the final archives
We first discuss the performance of the algorithms on located by MOSA for each of these problems are avail-
able online at http://www.dcs.ex.ac.uk/people/
2The values for the PISA variator parameters are: individ- kismith/mosa/results/tec/.
ualmutationprobability=1, individual recombinationprobability=1,
variablemutationprobability=1, variableswapprobability=0.5,
variablerecombinationprobability=1,etamutation=15,etarecombination=5.
IEEETRANSACTIONSONEVOLUTIONARYCOMPUTATION, DRAFT UNDERREVIEW 10
TABLEII
apparent from the ‘steps’ in the median archive distance. The
ANNEALINGSCHEDULES
current solution clearly leads the archive, particularly at later
iterations when the computational temperature is low and the
Problem Run Timeto
Problem length Tk =10(cid:0)5 search is effectively a greedy search.
DTLZ1 5000 3000
DTLZ2 1000 500
B. DTLZ 2
DTLZ3 15000 10000
DTLZ4 5000 3000 Figure7presentsthearchiveresultingfromarepresentative
DTLZ5 1000 500
run of the algorithms on problem DTLZ2 for 1000 function
DTLZ6 5000 3000
DTLZ7 9000 6000 evaluations and a plot of the distances from the true front,
whichisthe eighthofa sphericalshellofradius1, centredon
the origin, lying in the positive octant. As the figure shows,
the archive lies close to the optimal front for each of the
A. DTLZ 1
algorithms, with MOSA significantly closer than the other
algorithms.
Figure 5 shows views in objective space of the archive
We remark that this problem, and several others of the
obtained from a single run of each of the algorithms on test
DTLZsuitewithoutaplethoraoflocalfronts,canbesuccess-
problems DTLZ1 after 5000 objective evaluations, together
fully treated with a rapid cooling schedule, as used here. Due
with plotsshowingthe distanceofthe membersof eachset to
to the ease of convergence to the true front on this problem,
thetrueParetofront.Foreachalgorithm,theplottedresultsare
we anticipate that any multi-objective optimiser will be able
those which have the median distance of solutions to the true
to produce a set of solutions close to the true front although
front out of a series of 20 runs; this ensures that the results
the density and coverage may vary significantly, as is the
presented are representative of the series. The true front for
casehere. Figure14 illustratesthat, while allthree algorithms
DTLZ1 is the segment of the plane passing through 0:5 on
have convergedclose to the true front, MOSA is significantly
each of the objective space coordinate axes, and it can be
closer than NSGA-II or Nam & Parks annealer. The volume
seen that the majority of solutions generated by MOSA lie
measure plot shows that the archive produced by MOSA also
very close to the front. This test problem has a large number
((cid:25) 115) of local fronts which lie as planes parallel to and has a greater coverage/density of solutions; even after only
1000 evaluations, the archive size plot clearly illustrates that
further from the origin than P; the existence of these fronts
MOSA has already converged very close to the true front
is evident from the histogram of the distances from P which
and is searching across the front improving the coverage and
shows solutions clustered at two distinct distances for MOSA
density.
and several for NSGA-II (this effect is less marked on the
While knowledge about the applicability of a short anneal-
Nam & Park front, where the solutions are distributed more
ing schedule would not be initially available for typical real-
evenlyacrossmanyfrontswhicharecloseinobjectivespace).
world problems, we anticipate that, for real-world problems,
It seems likely that it is these local fronts which preventNam
theannealerwouldberunwithaveryrapidannealingschedule
& Park’s annealer and NSGA-II from converging on the true
initially to discover if the problem were searchable in this
front,sinceinlaterproblemswithoutthisfeaturethedifference
manner.
in performance between the three algorithms is, while still
significant, much less extreme. Figure 14 provides, for each
test problem, box plots comparing the averagedistance of the C. DTLZ 3
archive to the true front, the volume measure of the archive
A striking example of the annealer’s performance is pro-
and the size of the archive(which is a fixed value for NSGA-
videdinFigure8,whereitsevaluationonDTLZ3isshownfor
II due to the constrained nature of the algorithm). For this
15000 function evaluations. The Pareto front here is again an
DTLZ1 problem, the figure clearly illustrates that MOSA has
eighth of a spherical shell, preceded by multiple local fronts,
notonlyconvergedtoasetveryclosetothetruefrontbutthat
of the same order as DTLZ1. The computational archive is
thefrontisalsowellcoveredasshownbythevolumemeasure
converged to within 0:01 of the true front, as supported by
results; size the MOSA archive is unconstrained in size it has
the histogram of solution distances from the front in Figure
beenabletogeneratealargearchiveclosetoandcoveringwell
8. Consistent with the findings by Deb et al. [30] NSGA-
the true front. We observe that the annealer on this problem
II had failed to converge (Deb et al. comment that in their
convergestoalocalfront,spreadsacrossituntilaperturbation
experiments that NSGA-II had still failed to converge after
‘breaksthrough’toafrontclosertoP afterwhichtheannealer
50000 function evaluations) and Nam and Park’s annealer
explores the nearer local front, adding solutions on this front
yieldsperformancesimilartoNSGA-II(asillustratedinFigure
to the archive and removing solutions on the previous local
14). Consistent with the previous problems, MOSA’s archive
front as they are dominated during the exploration. Figure 6
is shown to be large, dense and well covering in Figure 14.
shows the median, maximum and minimum (over 20 runs)
of the distance of the current point x to the true front P
D. DTLZ 4
versus iteration, together with the median (over 20 runs) of
the median distance of members of the archive F from P on The true Pareto front for this problem is again an eighth
a much longer set of runs. The presence of local fronts is of a spherical shell, but the solutions are unevenlydistributed