MACS: AN AGENT-BASED MEMETIC MULTIOBJECTIVE

arXiv:1206.1305v1 [cs.CE] 6 Jun 2012

OPTIMIZATION ALGORITHM APPLIED TO SPACE TRAJECTORY

DESIGN

Massimiliano Vasile – Federico Zuiani
Space Advanced Research Team, University of Strathclyde
James Weir Building 75 Montrose Street G1 1XJ, Glasgow, United Kingdom
massimiliano.vasile@strath.ac.uk, f.zuiani@aero.gla.ac.uk

ABSTRACT
This paper presents an algorithm for multiobjective optimization that blends together a
number of heuristics. A population of agents combines heuristics that aim at exploring
the search space both globally and in a neighborhood of each agent. These heuristics
are complemented with a combination of a local and global archive. The novel agentbased algorithm is tested at first on a set of standard problems and then on three specific
problems in space trajectory design. Its performance is compared against a number of
state-of-the-art multiobjective optimisation algorithms that use the Pareto dominance
as selection criterion: NSGA-II, PAES, MOPSO, MTS. The results demonstrate that
the agent-based search can identify parts of the Pareto set that the other algorithms
were not able to capture. Furthermore, convergence is statistically better although the
variance of the results is in some cases higher.

1

KEYWORDS
Multiobjective optimisation; trajectory optimisation; memetic algorithms; multiagent
systems.

1. INTRODUCTION
The design of a space mission steps through different phases of increasing complexity.
In the first phase, a trade-off analysis of several options is required. The trade-off
analysis compares and contrasts design solutions according to different criteria and
aims at selecting one or two options that satisfy mission requirements. In mathematical
terms, the problem can be formulated as a multiobjective optimization problem.
As part of the trade-off analysis, multiple transfer trajectories to the destination need to
be designed. Each transfer should be optimal with respect to a number of criteria. The
solution of the associated multiobjective optimization problem, has been addressed, by
many authors, with evolutionary techniques. Coverstone et al. [1] proposed the use
of multiobjective genetic algorithms for the optimal design of low-thrust trajectories.
Dachwald et al. proposed the combination of a neurocontroller and of a multiobjective
evolutionary algorithm for the design of low-thrust trajectories [2]. In 2005 a study
by Lee et al. [3] proposed the use of a Lyapunov controller with a multiobjective evolutionary algorithm for the design of low-thrust spirals. More recently, Schütze et al.
proposed some innovative techniques to solve multiobjective optimization problems
for multi-gravity low-thrust trajectories. Two of the interesting aspects of the work
of Schütze et al. are the archiving of ǫ- and ∆-approximated solutions, to the known
best Pareto front [4], and the deterministic pre-pruning of the search space [5]. In
2009, Delnitz et al. [6] proposed the use of multiobjective subdivision techniques for
the design of low-thrust transfers to the halo orbits around the L2 libration point in the
Earth-Moon system. Minisci et al. presented an interesting comparison between an
EDA-based algorithm, called MOPED, and NSGA-II on some constrained and unconstrained multi-impulse orbital transfer problems [7].
In this paper, a hybrid population-based approach that blends a number of heuristics is

2

proposed. In particular, the search for Pareto optimal solutions is carried out globally
by a population of agents implementing classical social heuristics and more locally by
a subpopulation implementing a number of individualistic actions. The reconstruction
of the set of Pareto optimal solutions is handled through two archives: a local and a
global one.
The individualistic actions presented in this paper are devised to allow each agent to
independently converge to the Pareto optimal set. Thus creating its own partial representation of the Pareto front. Therefore, they can be regarded as memetic mechanisms
associated to a single individual. It will be shown that individualistic actions significantly improve the performance of the algorithm.
The algorithm proposed in this paper is an extension of the Multi-Agent Collaborative
Search (MACS), initially proposed in [8, 9], to the solution of multiobjective optimisation problems. Such an extension required the modification of the selection criterion,
for both global and local moves, to handle Pareto dominance and the inclusion of new
heuristics to allow the agents to move toward and along the Pareto front. As part of
these new heuristics, this papers introduces a dual archiving mechanism for the management of locally and globally Pareto optimal solutions and an attraction mechanisms
that improves the convergence of the population.
The new algorithm is here applied to a set of known standard test cases and to three
space mission design problems. The space mission design cases in this paper consider
spacecraft equipped with a chemical engine and performing a multi-impulse transfer.
Although these cases are different from some of the above-mentioned examples, that
consider a low-thrust propulsion system, nonetheless the size and complexity of the
search space is comparable. Furthermore, it provides a first test benchmark for multiimpulsive problems that have been extensively studied in the single objective case but
for which only few comparative studies exist in the multiobjective case [7].
The paper is organised as follows: section two contains the general formulation of the
problem, the third section starts with a general introduction to the multi-agent collaborative search algorithm and heuristics before going into some of the implementation
details. Section four contains a set of comparative tests that demonstrates the effec-

3

tiveness of the heuristics implemented in MACS. The section briefly introduces the
algorithms against which MACS is compared and the two test benchmarks that are
used in the numerical experiments. It then defines the performance metrics and ends
with the results of the comparison.

2. PROBLEM FORMULATION
A general problem in multiobjective optimization is to find the feasible set of solutions
that satisfies the following problem:

min f (x)

(1)

x∈D


where D is a hyperrectangle defined as D = xj | xj ∈ [blj buj ] ⊆ R, j = 1, ..., n
and f is the vector function:

f : D → Rm ,

f (x) = [f1 (x), f2 (x), ..., fm (x)]T

(2)

The optimality of a particular solution is defined through the concept of dominance:
with reference to problem (1), a vector y ∈ D is dominated by a vector x ∈ D if
fj (x) < fj (y) for all j = 1, ..., m. The relation x ≺ y states that x dominates y.
Starting from the concept of dominance, it is possible to associate, to each solution in
a set, the scalar dominance index:

Id (xj ) = |{i | i ∧ j ∈ Np ∧ xi ≺ xj }|

(3)

where the symbol |.| is used to denote the cardinality of a set and Np is the set of the
indices of all the solutions. All non-dominated and feasible solutions form the set:

X = {x ∈ D | Id (x) = 0}

(4)

Therefore, the solution of problem (1) translates into finding the elements of X. If X
is made of a collection of compact sets of finite measure in Rn , then once an element

4

of X is identified it makes sense to explore its neighborhood to look for other elements
of X. On the other hand, the set of non dominated solutions can be disconnected and
its elements can form islands in D. Hence, restarting the search process in unexplored
regions of D can increase the collection of elements of X.
The set X is the Pareto set and the corresponding image in criteria space is the Pareto
front. It is clear that in D there can be more than one Xl containing solutions that are
locally non-dominated, or locally Pareto optimal. The interest is, however, to find the
set Xg that contains globally non-dominated, or globally Pareto optimal, solutions.

3. MULTIAGENT COLLABORATIVE SEARCH
The key motivation behind the development of multi-agent collaborative search was to
combine local and global search in a coordinated way such that local convergence is
improved while retaining global exploration [9]. This combination of local and global
search is achieved by endowing a set of agents with a repertoire of actions producing
either the sampling of the whole search space or the exploration of a neighborhood
of each agent. More precisely, in the following, global exploration moves will be
called collaborative actions while local moves will be called individualistic actions.
Note that not all the individualistic actions, described in this paper, aim at exploring
a neighborhood of each agent, though. The algorithm presented in this paper is a
modification of MACS to tackle multiobjective optimization problems. In this section,
the key heuristics underneath MACS will be described together with their modification
to handle problem (1) and reconstruct X.

3.1. GENERAL ALGORITHM DESCRIPTION
A population P0 of npop virtual agents, one for each solution vector xi , with i =
1, ..., npop , is deployed in D. The population evolves through a number of generations.
At every generation k, the dominance index (3) of each agent xi,k in the population
Pk is evaluated. The agents with dominance index Id = 0 form a set Xk of nondominated solutions. Hence, problem (1) translates into finding a series of sets Xk
such that Xk → Xg for k → kmax with kmax possibly a finite number.
5

The position of each agent in D is updated through a number of heuristics. Some
are called collaborative actions because are derived from the interaction of at least two
agents and involve the entire population at one time. The general collaborative heuristic
can be expressed in the following form:

xk = xk + S(xk + uk )uk

(5)

where uk depends on the other agents in the population and S is a selection function
which yields 0 if the candidate point xk + uk is not selected or 1 if it is selected (see
Section 3.2). In this implementation a candidate point is selected if its dominance index
is better or equal than the one of xk . A first restart mechanism is then implemented to
avoid crowding. This restart mechanism is taken from [9] and prevents the agents from
overlapping or getting too close. It is governed by the crowding factor wc that defines
the minimum acceptable normalized distance between agents. Note that, this heuristic
increases the uniform sampling rate of D, when activated, thus favoring exploration.
On the other hand, by setting wc small the agents are more directed towards local
convergence.
After all the collaborative and restart actions have been implemented, the resulting
updated population Pk is ranked according to Id and split in two subpopulations: Pku
and Pkl . The agents in each subpopulation implement sets of, so called, individualistic
actions to collect samples of the surrounding space and to modify their current location.
In particular, the last npop − fe npop agents belong to Pku and implement heuristics that
can be expressed in a form similar to Eq. (5) but with uk that depends only on xk .
The remaining fe npop agents belong to Pkl and implement a mix of actions that aim
at either improving their location or exploring the neighborhood Nρ (xi,k ), with i =
1, ..., fe npop . Nρ (xi,k ) is a hyperectangle centered in xi,k . The intersection Nρ (xi,k )∩
D represents the local region around agent xi,k that one wants to explore. The size of
Nρ (xi,k ) is specified by the value ρ(xi,k ): the ith edge of Nρ has length 2ρ(xi,k ) max{buj −
xi,k [j], xi,k [j] − blj }.
The agents in Pkl generate a number of perturbed solutions ys for each xi,k , with s =

6

1, ..., smax . These solutions are collected in a local archive Al and a dominance index
is computed for all the elements in Al . If at least one element ys ∈ Al has Id = 0
then xi,k ← ys . If multiple elements of Al have Id = 0, then the one with the
largest variation, with respect to xi,k , in criteria space is taken. Figure 1(a) shows
three agents (circles) with a set of locally generated samples (stars) in their respective
neighborhoods (dashed square). The arrows indicate the direction of motion of each
agent. The figures shows also the local archive for the first agent Al1 and the target
global archive Xg .
The ρ value associated to an agent is updated at each iteration according to the rule
devised in [9]. Furthermore, a value s(xi,k ) is associated to each agent xi,k to specify
the number of samples allocated to the exploration of Nρ (xi,k ). This value is updated
at each iteration according to the rule devised in [9].
The adaptation of ρ is introduced to allow the agents to self-adjust the neighborhood
removing the need to set a priori the appropriate size of Nρ (xi,k ). The consequence
of this adaptation is an intensification of the local search by some agents while the
others are still exploring. In this respect, MACS works opposite to Variable Neighborhood Search heuristics, where the neighborhood is adapted to improve global exploration, and differently than Basin Hopping heuristics in which the neighborhood is
fixed. Similarly, the adaptation of s(xi,k ) avoids setting a priori an arbitrary number of
individualistic moves and has the effect of avoiding an excessive sampling of Nρ (xi,k )
when ρ is small. The value of s(xi,k ) is initialized to the maximum number of allowable individualistic moves smax . The value of smax is here set equal to the number
of dimensions n. This choice is motivated by the fact that a gradient-based method
would evaluate the function a minimum of n times to compute an approximation of
the gradient with finite differences. Note that the set of individualistic actions allows
the agents to independently move towards and within the set X, although no specific
mechanism is defined as in [10]. On the other hand the mechanisms proposed in [10]
could further improve the local search and will be the subject of a future investigation.
All the elements of Xk found during one generation are stored in a global archive Ag .
The elements in Ag are a selection of the elements collected in all the local archives.

7

Figure 1: Illustration of the a) local moves and archive and b) global moves and
archive.

Figure 1(b) illustrates three agents performing two social actions that yield two samples (black dots). The two samples together with the non-dominated solutions coming
from the local archive form the global archive. The archive Ag is used to implement an
attraction mechanism that improves the convergence of the worst agents (see Sections
3.4.1). During the global archiving process a second restart mechanism that reinitializes a portion of the population (bubble restart) is implemented. Even this second
restart mechanism is taken from [9] and avoids that, if ρ collapses to zero, the agent
keeps on sampling a null neighborhood.
Note that, within the MACS framework, other strategies can be assigned to the agents
to evaluate their moves in the case of multiple objective functions, for example a decomposition technique [11]. However, in this paper we develop the algorithm based
only on the use of the dominance index.

3.2. COLLABORATIVE ACTIONS
Collaborative actions define operations through which information is exchanged between pairs of agents. Consider a pair of agents x1 and x2 , with x1 ≺ x2 . One of the
two agents is selected at random in the worst half of the current population (from the
point of view of the property Id ), while the other is selected at random from the whole
population. Then, three different actions are performed. Two of them are defined by
adding to x1 a step uk defined as follows:
uk = αrt (x2 − x1 ),

(6)

and corresponding to: extrapolation on the side of x1 (α = −1, t = 1), with the further
constraint that the result must belong to the domain D (i.e., if the step uk leads out of
D, its size is reduced until we get back to D); interpolation (α = 1), where a random
point between x1 and x2 is sampled. In the latter case, the shape parameter t is defined

8

Algorithm 1 Main MACS algorithm
1:

2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:

Initialize a population P0 of npop agents in D, k = 0, number of function evaluations neval = 0, maximum number of function evaluations Ne , crowding factor
wc
for all i = 1, ..., npop do
xi,k = xi,k + S(xi,k + ui,k )ui,k
end for
Rank solutions in Pk according to Id
Re-initialize crowded agents according to the single agent restart mechanism
for all i = fe npop , ..., npop do
Generate np mutated copies of xi,k .
Evaluate the dominance of each mutated copy yp against xi,k , with p =
1, ..., np
if ∃p|yp ≺ xi,k then
p̄ = arg maxp kyp − xi,k k
xi,k ← yp̄
end if
end for
for all i = 1, ..., fe npop do
Generate s < smax individual actions us such that ys = xi,k + us
if ∃s|ys ≺ xi,k then
s̄ = arg maxs kys − xi,k k
xi,k ← ys̄
end if
Store candidate elements ys in the local archive Al
Update ρ(xi,k ) and s(xi,k )
end for
S S
S
Form Pk = Pkl Pku and Âg = Ag Al Pk
Compute Id of all the elements in Âg
Ag = {x|x ∈ Âg ∧ Id (x) = 0 ∧ kx − xAg k > wc }
Re-initialize crowded agents in Pk according to the second restart mechanism
Compute attraction component to Ag for all xi,k ∈ Pk \ Xk
k =k+1
Termination Unless neval > Ne , GoTo Step 2

as follows:
t = 0.75

s(x1 ) − s(x2 )
+ 1.25
smax

(7)

The rationale behind this definition is that we are favoring moves which are closer to
the agent with a higher fitness value if the two agents have the same s value, while in
the case where the agent with highest fitness value has a s value much lower than that
of the other agent, we try to move away from it because a small s value indicates that
improvements close to the agent are difficult to detect.

9

The third operation is the recombination operator, a single-point crossover, where,
given the two agents: we randomly select a component j; split the two agents into
two parts, one from component 1 to component j and the other from component j + 1
to component n; and then we combine the two parts of each of the agents in order
to generate two new solutions. The three operations give rise to four new samples,
denoted by y1 , y2 , y3 , y4 . Then, Id is computed for the set x2 , y1 , y2 , y3 , y4 . The
element with Id = 0 becomes the new location of x2 in D.
3.3. INDIVIDUALISTIC ACTIONS
Once the collaborative actions have been implemented, each agent in Pku is mutated
a number of times: the lower the ranking the higher the number of mutations. The
mutation mechanisms is not different from the single objective case but the selection is
modified to use the dominance index rather than the objective values.
A mutation is simply a random vector uk such that xi,k + uk ∈ D. All the mutated
solution vectors are then compared to xi,k , i.e. Id is computed for the set made of
the mutated solutions and xi,k . If at least one element yp of the set has Id = 0 then
xi,k ← yp . If multiple elements of the set have Id = 0, then the one with the largest
variation, with respect to xi,k , in criteria space is taken.
Each agent in Pkl performs at most smax of the following individualistic actions: inertia, differential, random with line search. The overall procedure is summarized in
Algorithm 2 and each action is described in detail in the following subsections.

3.3.1. Inertia
If agent i has improved from generation k − 1 to generation k, then it follows the
direction of the improvement (possibly until it reaches the border of D), i.e., it performs
the following step:
ys = xi,k + λ̄∆I
where λ̄ = min{1, max{λ : ys ∈ D}} and ∆I = (xi,k − xi,k−1 ).

10

(8)

3.3.2. Differential
This step is inspired by Differential Evolution [12]. It is defined as follows: let xi1 ,k , xi2 ,k , xi3 ,k
be three randomly selected agents; then

ys = xi,k + e [xi1 ,k + F (xi3 ,k − xi2 ,k )]

(9)

with e a vector containing a random number of 0 and 1 (the product has to be intended
componentwise) with probability 0.8 and F = 0.8 in this implementation. For every
component ys [j] of ys that is outside the boundaries defining D then ys [j] = r(buj −
blj ) + blj , with r ∈ U (0, 1). Note that, although this action involves more than one
agent, its outcome is only compared to the other outcomes coming from the actions
performed by agent xi,k and therefore it is considered individualistic.
3.3.3. Random with Line Search
This move realizes a local exploration of the neighborhood N (xi,k ). It generates a first
random sample ys ∈ N (xi,k ). Then if ys is not an improvement, it generates a second
sample ys+1 by extrapolating on the side of the better one between ys and xi,k :


ys+1 = xi,k + λ̄ α2 rt (ys − xi,k ) + α1 (ys − xi,k )

(10)

with λ̄ = min{1, max{λ : ys+1 ∈ D}} and where α1 , α2 ∈ {−1, 0, 1}, r ∈
U (0, 1) and t is a shaping parameter which controls the magnitude of the displacement.
Here we use the parameter values α1 = 0, α2 = −1, t = 1, which corresponds to
extrapolation on the side of xi,k , and α1 = α2 = 1, t = 1, which corresponds to
extrapolation on the side of ys .
The outcome of the extrapolation is used to construct a second order one-dimensional
model of Id . The second order model is given by the quadratic function fl (σ) =
a1 σ 2 +a2 σ+a3 where σ is a coordinate along the ys+1 −ys direction. The coefficients
a1 , a2 and a3 are computed so that fl interpolates the values Id (ys ), Id (ys+1 ) and
Id (xi,k ). In particular, for σ = 0 fl = Id (ys ), for σ = 1 fl = Id (ys+1 ) and for

11

σ = (xi,k − ys )/kxi,k − ys k fl = Id (xi,k ). Then, a new sample ys+2 is taken at the
minimum of the second-order model along the ys+1 − ys direction.
Algorithm 2 Individual Actions in Pkl
1:
2:

3:
4:

5:
6:
7:

8:
9:
10:

s = 1, stop = 0
if xi,k ≺ xi,k−1 then
ys = xi,k + λ̄(xi,k − xi,k−1 )
with λ̄ = min{1, max{λ : ys ∈ D}}.
end if
if xi,k ≺ ys then
s=s+1
ys = e[xi,k − (xi1 ,k + (xi3 ,k − xi2 ,k ))]
∀j|ys (j) ∈
/ D, ys (j) = r(bu (j) − bl (j)) + bl (j),
with j = 1, ..., n and r ∈ U (0, 1)
else stop = 1
end if
if xi,k ≺ ys then
s=s+1
Generate ys ∈ Nρ (xi,k ).
Compute ys+1 = xi,k + λ̄rt (xi,k − ys )
with λ̄ = min{1, max{λ : ys ∈ D}}
and r ∈ U (0, 1).
Compute ys+2 = λ̄σmin (ys+1 − ys )/k(ys+1 − ys )k,
with σmin = arg minσ {a1 σ 2 + a2 σ + Id (ys )},
and λ̄ = min{1, max{λ : ys+2 ∈ D}}.
s=s+2
else stop = 1
end if
Termination Unless s > smax or stop = 1 , GoTo Step 4

The position of xi,k in D is then updated with the ys that has Id = 0 and the longest
vector difference in the criteria space with respect to xi,k . The displaced vectors ys
generated by the agents in Pkl are not discarded but contribute to a local archive Al ,
one for each agent, except for the one selected to update the location of xi,k . In order
to rank the ys , the following modified dominance index is used:

Iˆd (xi,k ) = j | fj (ys ) = fj (xi,k ) κ+

j | fj (ys ) > fj (xi,k )

(11)

where κ is equal to one if there is at least one component of f(xi,k ) = [f1 , f2 , ..., fm ]T
which is better than the corresponding component of f(ys ), and is equal to zero other-

12

wise.
Now, if for the sth outcome, the dominance index in Eq. (11) is not zero but is lower
than the number of components of the objective vector, then the agent xi,k is only
partially dominating the sth outcome. Among all the partially dominated outcomes
with the same dominance index we select the one that satisfies the condition:

min
s


f(xi,k )) − f(ys ) , e

where e is the unit vector of dimension m, e =

[1,1,1,...,1]T
√
m

(12)

. All the non-dominated and

selected partially dominated solutions form the local archive Al .
3.4. THE LOCAL AND GLOBAL ARCHIVES AL AND AG
Since the outcomes of one agent could dominate other agents or the outcomes of other
agents, at the end of each generation, every Al and the whole population Pk are added
to the current global archive Ag . The global Ag contains Xk , the current best estimate of Xg . The dominance index in Eq.(3) is then computed for all the elements
S
S
in Âg = Ag l Al Pk and only the non-dominated ones with crowding distance
kxi,k − xAg k > wc are preserved (where xAg is an element of Ag ).
3.4.1. Attraction
The archive Ag is used to direct the movements of those agents that are outside Xk .
All agents, for which Id 6= 0 at step k, are assigned the position of the elements in Ag
and their inertia component is recomputed as:

∆I = r(xAg − xi,k )

(13)

More precisely, the elements in the archive are ranked according to their reciprocal distance or crowding factor. Then, every agent for which Id 6= 0 picks the least crowded
element xAg not already picked by any other agent.

13

3.5. STOPPING RULE
The search is stopped when a prefixed number Ne of function evaluations is reached.
At termination of the algorithm the whole final population is inserted into the archive
Ag .

4. COMPARATIVE TESTS
The proposed optimization approach was implemented in a software code, in Matlab,
called MACS. In previous works [8,9], MACS was tested on single objective optimization problems related to space trajectory design, showing good performances. In this
work, MACS was tested at first on a number of standard problems, found in literature,
and then on three typical space trajectory optimization problems.
This paper extends the results presented in Vasile and Zuiani (2010) [13] by adding
a broader suite of algorithms for multiobjective optimization to the comparison and a
different formulation of the performance metrics.

4.1. TESTED ALGORITHMS
MACS was compared against a number of state-of-the-art algorithms for multiobjective
optimization. For this analysis it was decided to take the basic version of the algorithms
that is available online. Further developments of the basic algorithms have not been
considered in this comparison and will be included in future works. Note that, all the
algorithms selected for this comparative analysis use Pareto dominance as selection
criterion.
The tested algorithms are: NSGA-II [14], MOPSO [15], PAES [16] and MTS [17]. A
short description of each algorithm with their basic parameters follows.

4.1.1. NSGA-II
The Non-Dominated Sorting Genetic Algorithm (NSGA-II) is a genetic algorithm
which uses the concept of dominance class (or depth) to rank the population. A crowding factor is then used to rank the individuals within each dominance class. Optimiza-

14

tion starts from a randomly generated initial population. The individuals in the population are sorted according to their level of Pareto dominance with respect to other individuals. To be more precise, a fitness value equal to 1 is assigned to the non-dominated
individuals. Non-dominated individuals form the first layer (or class). Those individuals dominated only by members of the first layer form the second one and are assigned
a fitness value of 2, and so on. In general, for dominated individuals, the fitness is
given by the number of dominating layers plus 1. A crowding factor is then assigned
to each individual in a given class. The crowding factor is computed as the sum of
the Euclidean distances, in criteria space, with respect to other individuals in the same
class, divided by the interval spanned by the population along each dimension of the
objective space. Inside each class, the individuals with the higher value of the crowding
parameter obtain a better rank than those with a lower one.
At every generation, binary tournament selection, recombination, and mutation operators are used to create an offspring of the current population. The combination of
the two is then sorted according to dominance first and then to crowding. The nondominated individuals with lowest crowding factor are then used to update the population.
The parameters to be set are the size of the population, the number of generations, the
crossover and mutation probability, pc and pm , and distribution indexes for crossover
and mutation, ηc and ηm , respectively. Three different ratios between population size
and number of generations were considered: 0.08, 0.33 and 0.75. The values pc and
pm were set to 0.9 and 0.2 respectively and kept constant for all the tests. The values,
5, 10 and 20, were considered for ηc , while tests were run for values of ηm equal to 5,
25 and 50.

4.1.2. MOPSO
MOPSO is an extension of Particle Swarm Optimization (PSO) to multiobjective problems. Pareto dominance is introduced in the selection criteria for the candidate solutions to update the population. MOPSO features an external archive which stores all the
non-dominated solutions and at the same time is used to guide the search process of the

15

swarm. This is done by introducing the possibility to direct the movement of a particle
towards one of the less crowded solutions in the archive. The solution space is subdivided into hypercubes through an adaptive grid. The solutions in the external archive
are thus reorganized in these hypercubes. The algorithm keeps track of the crowding
level of each hypercube and promotes movements towards less crowded areas. In a
similar manner, it also gives priority to the insertion of new non-dominated individuals
in less crowded areas if the external archive has already reached its predefined maximum size. For MOPSO three different ratios between population size and number of
generations were tested: 0.08 0.33 0.75. It was also tested with three different numbers
of subdivisions of the solution space: 10, 30, and 50. The inertia component in the
motion of the particles was set to 0.4 and the weights of the social and individualistic
components were se to 1

4.1.3. PAES
PAES is a (1 + 1) Evolution Strategy with the addition of an external archive to store
the current best approximation of the Pareto front. It adopts a population of only a
single chromosome which, at every iteration, generates a mutated copy. The algorithm
then preserves the non-dominated one between the parent and the candidate solution. If
none of the two dominates the other, the algorithm then checks their dominance index
with respect to the solutions in the archive. If also this comparison is inconclusive,
then the algorithm selects the one which resides in the less crowded region of the
objective space. To keep track of the crowding level of the objective space, the latter is
subdivided in an n-dimensional grid. Every time a new solution is added to (or removed
from) the archive, the crowding level of the corresponding grid cell is updated. PAES
has two main parameters that need to be set, the number of subdivisions in the space
grid and the mutation probability. Values of 1,2 and 4 were used for the former and 0.6
0.8 and 0.9 were used for the latter.

16

4.1.4. MTS
MTS is an algorithm based on pattern search. The algorithm first generates a population
of uniformly distributed individuals. At every iteration, a local search is performed by
a subset of individuals. Three different search patterns are included in the local search:
the first is a search along the direction of each decision variable with a fixed step length;
the second is analogous but the search is limited to one fourth of all the possible search
directions; the third one also searches along each direction but selects only solutions
which are evenly spaced on each dimension within a predetermined upper bound and
lower bound. At each iteration and for each individual, the three search patterns are
tested with few function evaluations to select the one which generates the best candidate
solutions for the current individual. The selected one is then used to perform the local
search which will update the individual itself. The step length along each direction,
which defines the size of the search neighbourhood for each individual, is increased if
the local search generated a non-dominated child and is decreased otherwise. When
the neighbourhood size reaches a predetermined minimum value, it is reset to 40% of
the size of the global search space. The non-dominated candidate solutions are then
used to update the best approximation of the global Pareto front. MTS was tested with
a population size of 20, 40 and 80 individuals.

4.2. PERFORMANCE METRICS
Two metrics were defined to evaluate the performance of the tested multiobjective optimizers:

Mspr =

Mconv =

Mp
1 X
fj − g i
min 100
Mp i=1 j∈Np
gi

Np
gj − fi
1 X
min 100
Np i=1 j∈Mp
gj

(14)

(15)

where Mp is the number of elements, with objective vector g, in the true global Pareto
front and Np is the number of elements, with objective vector f , in the Pareto front that
a given algorithm is producing. Although similar, the two metrics are measuring two
17

different things: Mspr is the sum, over all the elements in the global Pareto front, of the
minimum distance of all the elements in the Pareto front Np from the the ith element in
the global Pareto front. Mconv , instead, is the sum, over all the elements in the Pareto
front Np , of the minimum distance of the elements in the global Pareto front from the
ith element in the Pareto front Np .
Therefore, if Np is only a partial representation of the global Pareto front but is a
very accurate partial representation, then metric Mspr would give a high value and
metric Mconv a low value. If both metrics are high then the Pareto front Np is partial
and poorly accurate. The index Mconv is similar to the mean Euclidean distance [15],
although in Mconv the Euclidean distance is normalized with respect to the values of the
objective functions, while Mspr is similar to the generational distance [18], although
even for Mspr the distances are normalized.
Given n repeated runs of a given algorithm, we can define two performance indexes:
pconv = P (Mconv < tolconv ) or the probability that the index Mconv achieves a value
less than the threshold tolconv and pspr = P (Mspr < tolspr ) or the probability that
the index Mspr achieves a value less than the threshold tolconv .
According to the theory developed in [7, 19], 200 runs are sufficient to have a 95%
confidence that the true values of pconv and pspr are within a ±5% interval containing
their estimated value.
Performance index (14) and (15) tend to uniformly weigh every part of the front. This
is not a problem for index (14) but if only a relatively small portion of the front is
missed the value of performance index (15) might be only marginally affected. For this
reason, we slightly modified the computation of the indexes by taking only the MP∗ and
NP∗ solutions with a normalised distance in criteria space that was higher than 10−3 .
The global fronts used in the three space tests were built by taking a selection of about
2000 equispaced, in criteria space, nondominated solutions coming from all the 200
runs of all the algorithms.

18

4.3. PRELIMINARY TEST CASES
For the preliminary tests, two sets of functions, taken from the literature [14, 15], were
used and the performance of MACS was compared to the results in [14] and [15].
Therefore, for this set of tests, MTS was not included in the comparison. The function
used in this section can be found in Table 1.
Table 1: Multiobjective test functions

The first three functions were taken from [15]. Test cases Deb and Scha are two
examples of disconnected Pareto fronts, Deb2 is an example of problem presenting
multiple local Pareto fronts, 60 in this two dimensional case.
The last three functions are instead taken from [14]. Test case ZDT 2 has a concave
Pareto front with a moderately high-dimensional search space. Test case ZDT 6 has
a concave Pareto front but because of the irregular nature of f1 there is a strong bias
in the distribution of the solutions. Test case, ZDT4, with dimension 10, is commonly
recognized as one of the most challenging problems since it has 219 different local
Pareto fronts of which only one corresponds to the global Pareto-optimal front.
As a preliminary proof of the effectiveness of MACS, the average Euclidean distance
of 500 uniformly spaced points on the true optimal Pareto front from the solutions
stored in Ag by MACS was computed and compared to known results in the literature.
MACS was run 20 times to have a sample comparable to the one used for the other
algorithms. The global archive was limited to 200 elements to be consistent with [14].
The value of the crowding factor wc , the threshold ρtol and the convergence ρmin were
kept constant to 1e-5 in all the cases to provide good local convergence.
To be consistent with [15], on Deb, Scha and, Deb2, MACS was run respectively for
4000, 1200 and 3200 function evaluations. Only two agents were used for these lower
dimensional cases, with fe = 1/2. On test cases ZDT 2, ZDT 4 and ZDT 6, MACS
was run for a maximum of 25000 function evaluations to be consistent with [14], with
three agents and fe = 2/3 for ZDT 2 and four agents and fe = 3/4 on ZDT 4 and
ZDT 6.

19

The results on Deb, Scha and, Deb2 can be found in Table 2, while the results on
ZDT 2, ZDT 4 and ZDT 6 can be found in 3.
On all the smaller dimensional cases MACS performs comparably to MOPSO and better than PAES. It also performs than NSGA-II on Deb and Deb2. On Scha MACS
performs apparently worse than NSGA-II, although after inspection one can observe
that all the elements of the global archive Ag belong to the Pareto front but not uniformly distributed, hence the higher value of the Euclidean distance. On the higher
dimensional cases, MACS performs comparably to NSGA-II on ZDT2 but better than
all the others on ZDT4 and ZDT6. Note in particular the improved performance on
ZDT4.
Table 2: Comparison of the average Euclidean distances between 500 uniformly
space points on the optimal Pareto front for various optimization algorithms:
smaller dimension test problems.

Table 3: Comparison of the average Euclidean distances between 500 uniformly
space points on the optimal Pareto front for various optimization algorithms:
larger dimension test problems.

On the same six functions a different test was run to evaluate the performance of different variants of MACS. For all variants, the number of agents, fe , wc , ρtol and ρmin was
set as before, but instead of the mean Euclidean distance, the success rates pconv and
pspr were measured for each variant. The number of function evaluations for ZDT 2,
ZDT 4, ZDT 6 and Deb2 is the same as before, while for Scha and Deb it was reduced respectively to 600 and 1000 function evaluations given the good performance
of MACS already for this number of function evaluations. Each run was repeated 200
times to have good confidence in the values of pconv and pspr .
Four variants were tested and compared to the full version of MACS. Variant MACS no
local does not implement the individualistic moves and the local archive, variant MACS
ρ = 1 has no adaptivity on the neighborhood, its size is kept fixed to 1,variant MACS
ρ = 0.1 has the size of the neighborhood fixed to 0.1, variant MACS no attraction has
the attraction mechanisms not active.

20

The result can be found in Table 4. The values of tolconv and tolspr are respectively
0.001 and 0.0035 for ZDT 2, 0.003 and 0.005 for ZDT 4, 0.001 and 0.025 for ZDT 6,
0.0012 and 0.035 for Deb, 0.0013 and 0.04 for Scha, 0.0015 and 0.0045 for Deb2.
This thresholds were selected to highlight the differences among the various variants.
The table shows that the adaptation mechanism is beneficial in some cases although,
in others, fixing the value of ρ might be a better choice. This depends on the problem
and a general rule is difficult to derive at present. Other adaptation mechanisms could
further improve the performance.
The use of individualistic actions coupled with a local archive is instead fundamental,
so is the use of the attraction mechanism. Note, however, how the attraction mechanism penalizes the spreading on biassed problems like ZDT 6, this is expected as it
accelerates convergence.
Table 4: Comparison of different variants of MACS.

4.4. APPLICATION TO SPACE TRAJECTORY DESIGN
In this section we present the application of MACS to three space trajectory problems: a two-impulse orbit transfer from a Low Earth Orbit (LEO) to a high-eccentricity
Molniya-like orbit, a three-impulse transfer from a LEO to Geostationary Earth Orbit
(GEO) and a multi-gravity assist transfer to Saturn equivalent to the transfer trajectory
of the Cassini mission. The first two cases are taken from the work of Minisci et al. [7].
In the two-impulse case, the spacecraft departs at time t0 from a circular orbit around
the Earth (the gravity constant is µE = 3.9860105 km3 s−2 ) with radius r0 = 6721
km and at time tf is injected into an elliptical orbit with eccentricity eT = 0.667
and semimajor axis aT = 26610 km. The transfer arc is computed as the solution
of a Lambert’s problem [20] and the objective functions are the transfer time T =
tf − t0 and the sum of the two norms of the velocity variations at the beginning and
at the end of the transfer arc ∆vtot . The objectives are functions of the solution vector
x = [t0 tf ]T ∈ D ⊂ R2 The search space D is defined by the following intervals
t0 ∈ [0 10.8], and tf ∈ [0.03 10.8].
21

In the three-impulse case, the spacecraft departs at time t0 from a circular orbit around
the Earth with radius r0 = 7000 km and after a transfer time T = t1 + t2 is injected
into a circular orbit with radius rf = 42000. An intermediate manoeuvre is performed
at time t0 + t1 and at position defined in polar coordinates by the radius r1 and the
angle θ1 . The objective functions are the total transfer time T and the sum of the
three impulses ∆vtot . The solution vector in this case is x = [t0 , t1 , r1 , θ1 , tf ]T ∈
D ⊂ R5 . The search space D is defined by the following intervals t0 ∈ [0 1.62],
t1 ∈ [0.03 21.54], r1 ∈ [7010 105410], θ1 ∈ [0.01 2π − 0.01], and t2 ∈ [0.03 21.54].
The Cassini case consists of 5 transfer arcs connecting a departure planet, the Earth,
to the destination planet, Saturn, through a sequence of swing-by’s with the planets:
Venus, Venus, Earth, Jupiter. Each transfer arc is computed as the solution of a Lambert’s problem [21] given the departure time from planet Pi and the arrival time at
planet Pi+1 . The solution of the Lambert’s problems yields the required incoming
and outgoing velocities at each swing-by planet vin and vrout . The swing-by is modeled through a linked-conic approximation with powered maneuvers [22], i.e., the mismatch between the required outgoing velocity vrout and the achievable outgoing velocity vaout is compensated through a ∆v maneuver at the pericenter of the gravity
assist hyperbola. The whole trajectory is completely defined by the departure time t0
and the transfer time for each leg Ti , with i = 1, ..., 5. The normalized radius of the
pericenter rp,i of each swing-by hyperbola is derived a posteriori once each powered
swing-by manoeuvre is computed. Thus, a constraint on each pericenter radius has to
be introduced during the search for an optimal solution. In order to take into account
this constraint, one of the objective functions is augmented with the weighted violation
of the constraints:

f (x) = ∆v0 +

4
X

∆vi + ∆vf +

4
X

wi (rp,i − rpmin,i )2

(16)

i=1

i=1

for a solution vector x = [t0 , T1 , T2 , T3 , T4 , T5 ]T . The objective functions are, the toP
tal transfer time T = 5i Ti and f (x). The minimum normalized pericenter radii are
rpmin,1 = 1.0496, rpmin,2 = 1.0496, rpmin,3 = 1.0627 and rpmin,4 = 9.3925. The

22

search space D is defined by the following intervals: t0 ∈ [−1000, 0]MJD2000, T1 ∈
[30, 400]d, T2 ∈ [100, 470]d, T3 ∈ [30, 400]d, T4 ∈ [400, 2000]d, T5 ∈ [1000, 6000]d.
The best known solution for the single objective minimization of f (x) is fbest =
4.9307 km/s, with xbest = [−789.753, 158.2993, 449.3859, 54.7060, 1024.5896, 4552.7054]T .
4.4.1. Test Results
For this second set of tests, each algorithm was run for the same number of function
evaluations. In particular, consistent with the tests performed in the work of Minisci et
al., [7] we used 2000 function evaluations for the two-impulse case and 30000 for the
three-impulse case. For the Cassini case, instead, the algorithms were run for 180000,
300000 and 600000 function evaluations.
Note that the version of all the algorithms used in this second set of tests is the one that
is freely available online, written in c/c++. We tried in all cases to stick to the available
instructions and recommendations by the author to avoid any bias in the comparison.
The thresholds values for the two impulse cases was taken from [7] and is tolconv =
0.1, tolspr = 2.5. For the three-impulse case instead we considered tolconv = 5.0,
tolspr = 5.0. For the Cassini case we used tolconv = 0.75, tolspr = 5, instead. These
values were selected after looking at the dispersion of the results over 200 runs. Lower
values would result in a zero value of the performance indexes of all the algorithms,
which is not very significant for a comparison.
MACS was tuned on the three-impulse case. In particular, the crowding factor wc , the
threshold ρtol and the convergence ρmin were kept constant to 1e-5, which is below
the required local convergence accuracy, while fe and npop were changed. A value of
1e-5 is expected to provide good local convergence and good density of the samples
belonging to the Pareto front. Table 5 reports the value of performance indexes pconv
and pspr over 200 runs of MACS with different settings. The index pspr and the index
pconv have different, almost opposite, trends. However, it was decided to select the
setting that provides the best convergence,i.e. npop = 15 and fe = 1/3. This setting
will be used for all the tests in this paper.
On top of the complete algorithm, two variants of MACS were tested: one without

23

individualistic moves and local archive, denoted as no local in the tables, and one with
no attraction towards the global archive Ag , denoted as no att in the tables. Only these
two variants are tested on these cases as they displayed the most significant impact in
the previous standard test cases and more importantly were designed specifically to
improve performances.
NSGA-II, PAES, MOPSO and MTS were tuned as well on the three-impulse case. In
particular, for NSGA-II the best result was obtained for 150 individuals and can be
found in Table 6. A similar result could be obtained for MOPSO, see Table 7. For
MTS only the population was changed while the number of individuals performing
local moves was kept constant to 5. The results of the tuning of MTS can be found in
Table 9. For the tuning of PAES the results can be found in Table 8.
All the parameters tuned in the three impulse case were kept constant except for the
population size of NSGA-II and MOPSO. The size of the population of NSGA-II and
MOPSO was set to 100 and 40 respectively on the two impulse case and was increased
with the number of function evaluations in the Cassini case. In particular for NSGAII the following ratios between population size and number of function evaluations
was used: 272/180000, 353/300000, 500/600000. For MOPSO the following ratios
between population size and number of function evaluations was used: 224/180000,
447/300000, 665/600000. This might not be the best way to set the population size for
these two algorithms but it is the one that provided better performance in these tests.
Note that the size of the global archive for MACS was constrained to be lower than
the size of the population of NSGA-II, in order to avoid any bias in the computation of
Mspr .
The performance of all the algorithms on the Cassini case can be found in Table 12 for
a variable number of function evaluations.
Figure 2: Three-impulse test case: a) Complete Pareto front, b) close-up of the
Pareto fronts

Table 5: Indexes pconv and pspr for different settings of MACS

24

Figure 3: Cassini test case: a) Complete Pareto front, b) close-up of the Pareto
fronts

For the three-impulse case, MACS was able to identify an extended Pareto front (see
Fig.2(a) and Fig. 2(b) where all the non-dominated solutions from all the 200 runs are
compared to the global front), compared to the results in [7]. The gap in the Pareto
front is probably due to a limited spreading of the solutions in that region. Note the
cusp due the transition between the condition in which 2-impulse solutions are optimal
and the condition in which 3-impulse solutions are optimal.
Table 11 summarizes the results of all the tested algorithms on the two-impulse case.
The average value of the performance metrics is reported with, in brackets, the associated variance over 200 runs. The two-impulse case is quite easy and all the algorithms
have no problems identifying the front. However, MACS displays a better convergence
than the other algorithms while the spreading of MOPSO and MTS is superior to the
one of MACS.
The three-impulse case is instead more problematic (see Table 10). NSGA-II is not
able to converge to the upper-left part of the front and therefore the convergence is 0
and the spreading is comparable to the one of MACS. All the other algorithms perform
poorly with a value of almost 0 for the performance indexes. This is mainly due to the
fact that no one can identify the vertical part of the front.
Note that the long tail identified by NSGA-II is actually dominated by two points belonging to the global front (see 2(b)).
Table 12 reports the statistics for the Cassini problem. On top of the performance of the
three variants tested on the other two problems, the table reports also the result for 10
agents and fe = 5. For all numbers of function evaluations MACS has better spreading
than NSGA-II because NSGA-II converges to a local Pareto front. Nonetheless NSGAII displays a more regular behavior and a better convergence for low number of function
evaluations although it never reaches the best front. The Pareto fronts are represented
in Fig. 3(a) and Fig. 3(b) (even in this case the figures represent all the non-dominated
solutions coming from all the 200 runs). Note that, the minimum f returned by MACS
25

is the best known solution for the single objective version of this problem (see [9]). All
the other algorithms perform quite poorly on this case.
Of all the variants of MACS tested on these problems the full complete one is performing the best. As expected, in all three cases, removing the individualistic moves
severely penalizes both convergence and spreading. It is interesting to note that removing the attraction towards the global front is also comparatively bad. On the two
impulse case it does not impact the spreading but reduces the convergence, while on
the Cassini case, it reduces mean and variance but the success rates are zero. The observable reason is that MACS converges slower but more uniformly to a local Pareto
front.
Finally, it should be noted that mean and variance seem not to capture the actual performance of the algorithms. In particular they do not capture the ability to identify the
whole Pareto front as the success rates instead do.
Table 6: NSGAII tuning on the 3-impulse case

Table 7: PAES tuning on the 3-impulse case

5. CONCLUSIONS
In this paper we presented a hybrid evolutionary algorithm for multiobjective optimization problems. The effectiveness of the hybrid algorithm, implemented in a code
called MACS, was demonstrated at first on a set of standard problems and then its performance was compared against NSGA-II, PAES, MOPSO and MTS on three space
trajectory design problems. The results are encouraging as, for the same computational
effort (measured in number of function evaluations, MACS was converging more accuTable 8: MOPSO tuning on the 3-impulse case.
Table 9: MTS tuning on the 3-impulse

26

rately than NSGA-II on the two-impulse case and managed to find a previously undiscovered part of the Pareto front of the three-impulse case. As a consequence, on the
three-impulse case, MACS, has better performance metrics than the other algorithms.
On the Cassini case NSGA-II appears to converge better to some parts of the front although MACS yielded solutions with better f and identifies once more a part of the
front that NSGA-II cannot attain. PAES and MTS do not perform well on the Cassini
case, while MOPSO converges well locally but, with the settings employed in this
study, yielded a very poor spreading.
From the experimental tests in this paper we can argue that the following mechanisms
seem to be particularly effective: the use of individual local actions with a local archive
as they allow the individuals to move towards and within the Pareto set; the use of an
attraction mechanism as it accelerates convergence.
Finally it should be noted that all the algorithms tested in this study use the Pareto
dominance as selection criterion. Different criteria, like the decomposition in scalar
subproblems, can be equality implemented in MACS, without disrupting its working
principles, and lead to different performance results.

ACKNOWLEDGEMENTS
The authors would like to thank Dr. Edmondo Minisci and Dr. Guilio Avanzini for the
two and three impulse cases and the helpful advice.

REFERENCES
[1] V. Coverstone-Caroll, J.W. Hartmann, and W.M. Mason. Optimal multi-objective
low-thrust spacecraft trajectories. Computer Methods in Applied Mechanics and
Engineering, 186:387–402, 2000.
[2] B. Dachwald. Optimization of interplanetary solar sailcraft trajectories using evolutionary neurocontrol. Journal of Guidance, and Dynamics, February 2004.
[3] S. Lee, P. von Allmen, W. Fink, A.E. Petropoulos, and R.J. Terrile. Multiobjective evolutionary algorithms for low-thrust orbit transfer optimization. In
27

Proceedings of the Genetic and Evolutionary Computation Conference (GECCO
2005), Washington DC, USA, June 25–29 2005.
[4] O. Schütze, M. M. Vasile, and C.A. Coello Coello. Approximate solutions in
space mission design. In Parallel Problem Solving from Nature (PPSN 2008),
Dortmund, Germany, September 13–17 2008.
[5] O. Schütze, M. Vasile, O. Junge, M. Dellnitz, and D. Izzo. Designing optimal
low-thrust gravity-assist trajectories using space pruning and a multi-objective
approach. Engineering Optimization, 41(2), February 2009.
[6] M. Dellnitz, S. Ober-Blbaum, M. Post, O. Schtze, and Bianca Thiere. A multiobjective approach to the design of low thrust space trajectories using optimal
control. Celestial Mechanics and Dynamical Astronomy, 105(1), February 2009.
[7] E. Minisci and G. Avanzini. Optimisation of orbit transfer manoeuvres as a
test benchmark for evolutionary algorithms. In Proceedings of the 2009 IEEE
Congress on Evolutionary Computation (CEC2009), Trondheim, Norway, May
2009.
[8] M. Vasile. A behavioral-based meta-heuristic for robust global trajectory optimization. In Proceedings of the 2007 IEEE Congress on Evolutionary Computation (CEC2007), Singapore, September 2007.
[9] M. Vasile and M. Locatelli. A hybrid multiagent approach for global trajectory
optimization. Journal of Global Optimization, 44(4):461–479, August. 2009.
[10] Oliver Schütze, Adriana Lara, Gustavo Sanchez, and Carlos A. Coello Coello.
HCS: A new local search strategy for memetic multi-objective evolutionary algorithms. IEEE Transactions on Evolutionary Computation (to appear), 2010.
[11] Q. Zhang and H. Li. Moea/d: A multiobjective evolutionary algorithm based on
decomposition. IEEE Transactions on Evolutionary Computation, 11(6), November.

28

[12] K.V. Price, R.M. Storn, and J.A. Lampinen. Differential Evolution. A Practical
Approach to Global Optimization. Natural Computing Series. Springer, 2005.
[13] M. Vasile and F. Zuiani. A hybrid multiobjective optimization algorithm applied
to space trajectory optimization. In Proceedings of the IEEE International Conference on Evolutionary Computation, Barcelona, Spain, July 2010.
[14] K. A. Deb, A. Pratap, and T. Meyarivan. Fast elitist multi-objective genetic algorithm: Nga-ii. Kangal report no. 200001, KanGAL, 2000.
[15] C. Coello and M. Lechuga. Mopso: A proposal for multiple objective particle
swarm optimization. Technical report evocinv-01-2001., CINVESTAV. Instituto
Politecnico Nacional. Col. San Pedro Zacatenco. Mexico., 2001.
[16] J.D. Knowles and D.W. Corne. The pareto archived evolution strategy : A new
baseline algorithm for pareto multiobjective optimisation. In Proceedings of the
IEEE International Conference on Evolutionary Computation, Washington DC,
US, 1999.
[17] Lin-Yu Tseng and Chun Chen. Multiple trajectory search for multiobjective optimization. In Proceedings of the IEEE International Conference on Evolutionary
Computation, Singapore, 25-28 September 2007.
[18] D.A. Van Veldhuizen and G.B. Lamont. Evolutionary computation and convergence to a pareto front. In Late Breaking papers at the Genetic Programming,
1998.
[19] M. Vasile, E. Minisci, and M. Locatelli. On testing global optimization algorithms
for space trajectory design. In AIAA/AAS Astrodynamic Specialists Conference,
Honolulu, Hawaii, USA, Aug 2008.
[20] G. Avanzini. A simple lambert algorithm. Journal of Guidance, Control, and
Dynamics, 31(6):1587–1594, Nov.–Dec. 2008.
[21] R. Battin. An Introduction to the Mathematics and Methods of Astrodynamics.
AIAA, 1999.
29

[22] D. R. Myatt, V.M. Becerra, S.J. Nasuto, and J.M. Bishop.
mization tools for mission analysis and design.

Final rept. esa ariadna itt

ao4532/18138/04/nl/mv,call03/4101, ESA/ESTEC, 2004.

30

Global opti-

Table 10: Summary of metrics Mconv and Mspr , and associated performance
indexes pconv and pspr , on the three impulse test cases.

Table 11: Metrics Mconv and Mspr , and associated performance indexes pconv
and pspr , on the two impulse test cases.

Table 12: Metrics Mconv and Mspr , and associated performance indexes pconv
and pspr , on the Cassini case.

31

NOMENCLATURE

Ag

global archive

Al

local archive

aT

semimajor axis

D

search space

eT

eccentricity

f

cost function

fe

fraction of the population doing local moves

Id

dominance index

Mconv

convergence metrics

Mspr

spreading metrics

Nρ

neighborhood of solution x

Ne

maximum number of allowed function evaluations

neval

number of function evaluations

npop

population size

Pi

i-th planet

Pk

population at generation k

pconv

percentage of success on convergence

pspr

percentage of success on spreading

r

random number

rp

pericentre radius

rpmin

minimum pericentre radius

S

selection function

s

resource index

T

transfer time

t0

departure time

ti

manoeuvre time

tf

final time

U

uniform distribution

32

u

variation of the solution x

X

Pareto optimal set

x

solution vector

y

mutate individual

wc

tolerance on the maximum allowable crowding

Greek symbols
∆v

variation of velocity

ρ

size of the neighborhood Nρ

µE

gravity constant

θi

true anomaly of manoeuvre i

33

LIST OF TABLES
Table 1: Multiobjective test functions
Scha
x ∈ [−5, 10]
Deb
x1 , x2 ∈ [0, 1]
α = 2; q = 4
Deb2
x1 ∈ [0, 1]
x2 ∈ [−30, 30]
ZDT2
xi ∈ [0, 1];
i = 1, . . . , n
n = 30
ZDT4
x1 ∈ [0, 1];
xi ∈ [−5, 5];
i = 2, . . . , n
n = 10
ZDT6
xi ∈ [0, 1];
i = 1, . . . , n
n = 10

2
f2 = (x
 − 5)
−x



−2 + x
f1 =
4−x



−4 + x
f1 = x1
h

if
if
if
if

f2 = (1 + 10x2 ) 1 −
i
x1
sin(2πqx
)
1
1+10x2

x≤1
1<x<3
3<x≤4
x>4


x1
1+10x2

α

−

f1 = x1
f2 = g(x1 , x2 )h(x1 ,q
x2 ); g(x1 , x2 ) = 11 + x22 − 10 cos(2πx2 )
n
1 − fg1 iff1 ≤ g
h(x1 , x2 ) =
0 otherwise
9 Pn
g = 1 + n−1
i=2 xi
f1 2
h=1−( g )
f1 = x1 ; f2 = gh
g = 1 + 10(n
q − 1) +
f1
g

h=1−
f1 = x1 ; f2 = gh

Pn

2
i=2 [xi

− 10 cos(2πqxi )];

q Pn
i=2 xi
g = 1 + 9 4 n−1
h = 1 − ( fg1 )2
f1 = 1 − exp(−4x1 ) sin6 (6πx1 ); f2 = gh

Table 2: Comparison of the average Euclidean distances between 500 uniformly
space points on the optimal Pareto front for various optimization algorithms:
smaller dimension test problems.
Approach
MACS
NSGA-II
PAES
MOPSO

Deb2

Scha

Deb

1.542e-3
(5.19e-4)
0.094644
(0.117608)
0.259664
(0.573286)
0.0011611
(0.0007205)

3.257e-3
(5.61e-4)
0.001594
(0.000122)
0.070003
(0.158081)
0.00147396
(0.00020178)

7.379e-4
(6.36e-5)
0.002536
(0.000138)
0.002881
(0.00213)
0.002057
(0.000286)

34

Table 3: Comparison of the average Euclidean distances between 500 uniformly
space points on the optimal Pareto front for various optimization algorithms:
larger dimension test problems.
Approach
MACS
NSGA-II
PAES

ZDT2

ZDT4

ZDT6

9.0896e-4
(4.0862e-5)
0.000824
(<1e-5)
0.126276
(0.036877)

0.0061
(0.0133)
0.513053
(0.118460)
0.854816
(0.527238)

0.0026
(0.0053)
0.296564
(0.013135)
0.085469
(0.006644)

Table 4: Comparison of different versions of MACS.
Approach

Metric

ZDT2

ZDT4

ZDT6

Scha

Deb

Deb2

MACS

pconv
pspr
pconv
pspr
pconv
pspr
pconv
pspr
pconv
pspr

83.5%
22.5%
14%
1%
84%
22%
56%
5%
21%
0%

75%
28%
0%
0%
22%
7%
42%
15%
0.5%
0.5%

77%
58.5%
45%
34%
78%
63%
57%
61%
14.5%
78.5%

73%
38.5%
0.5%
0%
37%
0%
78%
88%
0.5%
0%

70.5%
83%
72.5%
4%
92%
54%
85%
94.5%
88.5%
0%

60%
67.5 %
11%
15%
21%
38%
42%
74%
0%
0%

MACS
no local
MACS
ρ=1
MACS
ρ = 0.1
MACS
no attraction

Table 5: Indexes pconv and pspr for different settings of MACS on the 3-impulse
case
pconv

npop = 5

npop = 10

npop = 15

fe = 1/3
fe = 1/2
fe = 2/3

45.5%
48.0%
45.0%

55.5%
51.0%
52.5%

61.0%
55.5%
43.0%

pspr
fe = 1/3
fe = 1/2
fe = 2/3

npop = 5
68.5%
65.0%
67.5%

npop = 10
62.0%
57.0%
51.0%

npop = 15
56.0%
46.0%
36.5%

35

Table 6: NSGAII tuning on the 3-impulse case
Mean Mconv
ηc / ηm
5
10
20
Mean Mspr
ηc / ηm
5
10
20
pconv
ηc / ηm
5
10
20

5
36.1
32.3
31.7

25
38.3
39.4
39.6

50
43.0
40.6
42.5

5
6.77
5.91
5.78

25
7.24
7.50
7.50

50
8.08
7.81
8.16

5
0.0%
0.0%
0.0%

25
0.0%
0.0%
0.0%

50
0.0%
0.0%
0.0%

Var Mconv
ηc / ηm
5
10
20
Var Mspr
ηc / ηm
5
10
20
pspr
ηc / ηm
5
10
20

5
201.0
182.0
175.0

25
202.0
172.0
183.0

50
185.0
182.0
169.0

5
9.97
9.74
9.75

25
9.47
8.68
8.53

50
7.25
8.34
8.04

5
44.8%
57.8%
61.0%

25
37.7%
33.1%
32.5%

50
23.4%
29.2%
24.7%

Table 7: MOPSO tuning on the 3-impulse case.
Mean Mconv
Particles/Subdivisions
50
100
150
Mean Mspr
Particles/Subdivisions
50
100
150
pconv
Particles/Subdivisions
50
100
150

10
59.1
50.0
47.8

30
50.2
43.3
41.4

50
41.5
41.3
39.4

10
16.1
14.7
14.8

30
13.5
12.2
11.9

50
12.3
11.6
11.4

10
0%
0%
0%

30
0%
0%
0%

50
0%
0%
0%

36

Var Mconv
Particles/Subdivisions
50
100
150
Var Mspr
Particles/Subdivisions
50
100
150
pspr
Particles/Subdivisions
50
100
150

10
1080.0
591.0
562.1

30
1010.2
721.1
550.2

50
713.3
778.1
608.2

10
41.3
32.8
30.0

30
37.3
25.8
22.2

50
24.0
24.5
22.5

10
0.5%
0.5%
0.5%

30
2.5%
4.5%
2%

50
3.5%
4%
4%

Table 8: PAES tuning on the 3-impulse case
Mean Mconv
Subdivisions/Mutation
1
2
4
Mean Mspr
Subdivisions/Mutation
1
2
4
pconv
Subdivisions/Mutation
1
2
4

0.6
53.7
52.8
53.0

0.8
70.6
70.2
70.2

0.9
70.2
70.0
70.1

0.6
14.2
13.6
13.8

0.8
27.7
27.6
27.7

0.9
36.7
36.6
36.6

0.6
0.0%
0.0%
0.0%

0.8
0.0%
0.0%
0.0%

0.9
0.0%
0.0%
0.0%

Var Mconv
Subdivisions/Mutation
1
2
4
Var Mspr
Subdivisions/Mutation
1
2
4
pspr
Subdivisions/Mutation
1
2
4

0.6
525.0
479.0
453.0

0.8
275.0
266.0
266.0

0.9
297.0
305.0
311.0

0.6
20.0
17.5
17.2

0.8
14.3
14.6
15.8

0.9
17.3
16.8
17.0

0.6
0.0%
0.0%
0.0%

0.8
0.0%
0.0%
0.0%

0.9
0.0%
0.0%
0.0%

Table 9: MTS tuning on the 3-impulse
3imp
Mconv
Mspr

Population
Mean
Var
Mean
Var
pconv
pspr

20
17.8
97.6
12.6
34.7
1.0%
0.5%

40
22.6
87.8
19.9
26.2
0.0%
0.0%

80
23.6
73.2
18.4
18.6
0.0%
0.0%

Table 10: Summary of metrics Mconv and Mspr , and associated performance
indexes pconv and pspr , on the three impulse test cases.
Metric

MACS

MACS
no local

MACS
no att

NSGA-II

PAES

MOPSO

MTS

Mconv

5.53
(15.1)
5.25
(3.73)
61.0%
56.0%

7.58
(26.3)
6.03
(3.95)
40.5%
36%

154.7
(235.0)
9.16
(2.07)
0.0%
0.0%

31.7
(175.0)
5.78
(9.75)
0.0%
61.0%

53.0
(453.0)
13.8
(17.2)
0.0%
0.0%

39.4
(608.1)
11.4
(22.5)
0%
4.0%

17.8
(97.6)
12.6
(34.7)
1.0%
0.5%

Mspr
pconv
pspr

37

Table 11: Metrics Mconv and Mspr , and associated performance indexes pconv
and pspr , on the two impulse test cases.
Metric

MACS

MACS
no local

MACS
no att

NSGA-II

PAES

MOPSO

MTS

Mconv

0.0077
(1e-3)
2.89
(0.49)
98.5%
29.5%

0.0039
(1.4e-2)
6.87
(7.36)
91.5%
0%

0.0534
(9.5e-3)
3.08
0.943
84%
24%

0.283
(4.35e-3)
2.47
(0.119)
0%
62.5%

0.198
(0.332)
332.0
(2.61e4)
75.5%
0.5%

0.378
(0.0636)
2.11
(2.65)
9%
94.5%

0.151
(0.083)
1.95
(1.41)
57.5%
85.5%

Mspr
pconv
pspr

38

Table 12: Metrics Mconv and Mspr , and associated performance indexes pconv
and pspr , on the Cassini case.
Approach

Metric

180k

300k

600k

MACS
5/15

Mconv
Mspr
pconv
pspr

6.50 (229.1)
12.7 (126.0)
6%
27.5%

4.48 (107.1)
11.1 (83.1)
14%
31%

3.91 (62.6)
8.64 (35.5)
21.5%
41%

MACS
5/10

Mconv
Mspr
pconv
pspr

6.74 (217.9)
12.1 (83.0)
8%
26.0%

5.56 (136.4)
10.3 (63.7)
10%
31.5%

3.14 (50.1)
8.11 (30.0)
25.5%
45.5%

MACS
no local

Mconv
Mspr
pconv
pspr

13.5 (436.2)
31.1 (274.3)
1.0%
1.0%

10.2 (350.1)
27.9 (278.9)
1.5%
2.5%

7.86 (190.2)
21.9 (226.7)
2.5%
6%

MACS
no att

Mconv
Mspr
pconv
pspr

2.62 (1.46)
27.8 (83.2)
0.0%
0.0%

2.2 (0.88)
22.9 (58.0)
0.0%
0.0%

1.82 (0.478)
18.2 (28.0)
1.0%
0.0%

NSGA-II

Mconv
Mspr
pconv
pspr

2.43 (18.0)
11.6 (71.4)
17.5%
15.5%

1.99 (16.8)
11.0 (47.5)
24.0%
12.5%

1.24 (1.62)
8.78 (28.2)
29.0%
25.0%

MOPSO

Mconv
Mspr
pconv
pspr

2.62 (7.33)
28.0 (308.3)
0.5%
0.0%

2.4 (2.57)
24.6 (260.4)
1.0%
0.0%

2.14 (0.94)
21.8 (231.3)
1.0%
0.5%

PAES

Mconv
Mspr
pconv
pspr

24.0 (54.5)
30.1 (47.3)
0.0%
0.0%

19.8 (32.9)
26.0 (33.9)
0.0%
0.0%

15.2 (16.6)
21.4 (19.5)
0.0%
0.0%

MTS

Mconv
Mspr
pconv
pspr

3.71 (1.53)
18.1 (18.2)
0.0%
0.0%

3.39 (1.67)
15.6 (13.4)
0.0%
0.0%

3.02 (1.69)
13.1 (8.46)
0.0%
0.0%

39

LIST OF FIGURES

40

(a)

(b)

Figure 1: Illustration of the a) local moves and archive and b) global moves and
archive.

41

Pareto Front Comparison
30
Global
MACS
NSGA−II

25
20
T [h]

Part of the front not identified by NSGAII
15
10

Dominated tail

5
0
0

50

100

150
∆v

tot

200
[km/s]

250

300

350

(a)
Pareto Front Comparison
8
Global
MACS
NSGA−II

7
6

T [h]

5
4
Dominated tail

3
2
1
0
10

20

30
40
∆ vtot [km/s]

50

60

(b)

Figure 2: Three-impulse test case: a) Complete Pareto front, b) close-up of the
Pareto fronts

42

Pareto Front Comparison
7000
Global
MACS
NSGA−II

6000

T [day]

5000
4000
3000
2000
1000
0

50

100

150

200

250

f [km/s]

(a)
Pareto Front Comparison
Global
MACS
NSGA−II

6000
5500

T [day]

5000
4500
4000
3500
3000
2500
2000
4

5

6

7
f [km/s]

8

9

10

(b)

Figure 3: Cassini test case: a) Complete Pareto front, b) close-up of the Pareto
fronts

43

Pareto Front Comparison
7000
Global
NSGA−II
MACS

6000

T [day]

5000
4000
3000
2000
1000
0

50

100

150
f [km/s]

200

250

Close−up of the Cassini Pareto fronts
7000
Global
NSGA−II
MACS

T [day]

6000

5000

4000

3000

2000
5

5.5

6

6.5
f [km/s]

7

7.5

8

