# Asteroid Rendezvous Mission Design Using Multiobjective Particle Swarm Optimization

Authors: Ya-zhong Luo; Li-ni Zhou
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

A new preliminary trajectory design method for asteroid rendezvous mission using multiobjective optimization techniques is proposed. This method can overcome the disadvantages of the widely employed Pork-Chop method. The multiobjective integrated launch window and multi-impulse transfer trajectory design model is formulated, which employes minimum-fuel cost and minimum-time transfer as two objective functions. The multiobjective particle swarm optimization (MOPSO) is employed to locate the Pareto solution. The optimization results of two different asteroid mission designs show that the proposed approach can effectively and efficiently demonstrate the relations among the mission characteristic parameters such as launch time, transfer time, propellant cost, and number of maneuvers, which will provide very useful reference for practical asteroid mission design. Compared with the PCP method, the proposed approach is demonstrated to be able to provide much more easily used results, obtain better propellant-optimal solutions, and have much better efficiency. The MOPSO shows a very competitive performance with respect to the NSGA-II and the SPEA-II; besides a proposed boundary constraint optimization strategy is testified to be able to improve its performance.


---

## Body

## 1. Introduction

The optimization of interplanetary trajectories to an asteroid continues to arouse a great deal of interest [1–4]. Although the low-thrust propulsion is employed in asteroid rendezvous missions, an impulsive trajectory is always assumed for preliminary mission design and optimization, in which the access feasibility is evaluated, the launch window is determined, the gravity-assist maneuvers scheme is planned, and so forth. In a preliminary asteroid rendezvous trajectory design, the Pork-Chop method is widely used. This method employs the two-impulse algorithms including the optimal two-impulse noncoplanar transfer [1, 2], the classical Hohmann transfer [3], and the two-impulse Lambert algorithm [4] to design the transfer trajectory. The total Δv corresponding to different departure time and arrival time is then obtained and the contours of minimum total Δv are plotted to assist the designer to find the best launch window and transfer trajectory.This type of method is very intuitionistic and easily executed. However, only the two-impulse trajectory is investigated in this method, and as demonstrated by Lawden’s theory [5, 6], the two-impulse trajectory is not the propellant-optimal solution under most conditions. Besides, this method is in essence of an exhaustive searching method, its computation cost increases exponentially as the search space increases, and much human intervention is required as no orderliness exists in most of the contours.In this paper, a new preliminary trajectory design method for asteroid rendezvous mission using multiobjective optimization techniques is proposed. This method can overcome the disadvantages of the Pork-Chop method. The multiobjective integrated launch window and multi-impulse transfer trajectory design model is formulated, which employs minimum-fuel cost and minimum-time transfer as two objective functions. In this model, the Earth departure date, hyperbolic velocity, and the interplanetary transfer impulses are all chosen optimization design variables. The multiobjective particle swarm algorithm is employed to locate the Pareto solution, by which the relationships characteristics among the overall mission parameters can be effectively revealed.The particle swarm optimization (PSO) algorithm firstly introduced by Kennedy and Eberhart [7] in 1995 has been successfully applied in many fields of research. However, it is only recently applied to aerospace trajectories optimization [8–11]. The single-objective PSO algorithm has been testified as one successful spacecraft trajectory optimizer in designing impulsive interplanetary trajectories [8], low-thrust trajectories for asteroid exploration [9], Lyapunov and Halo (periodic) orbits [10], and multiple-burn rendezvous trajectories [11].The implementation of the PSO algorithm adopts a population of particles, whose behavior is affected by either the best local (i.e., within a certain neighborhood) or the best global individual. The relative simplicity of PSO and the fact that it is a population-based technique have made it a natural candidate to be extended for multiobjective optimization. In a survey paper in 2006 on multiobjective particle swarm optimization (MOPSO) [12], it was reported that there were currently over twenty-five different proposals of MOPSO reported in the specialized literature. The studies on MOPSO remain a very active area of research, and the MOPSO has been successfully applied to many practical multiobjective optimization problems, recently applied to robotics [13], industrial management [14], and chemical engineering [15]. It also has been applied in the domain of aerospace including airfoil shape optimization, complex physics/shape optimization, and multidisplinary design optimization [16]. However, it was seldom applied in solving multiobjective spacecraft trajectory optimization problems. The recent studies show that the single-objective PSO algorithm is one effective, reliable, and accurate spacecraft trajectory optimizer [8–11]. In this study, we will show that the MOPSO could be one successful optimizer for multiobjective spacecraft trajectory optimization problems.In summary, the main contribution of this paper is twofold. (1) A novel asteroid rendezvous mission design method using the multiobjective techniques is proposed. Compared with the current widely employed Pork-Chop method, the proposed approach is demonstrated to be able to provide much more easily used results, obtain better propellant-optimal solutions, and have much better efficiency. (2) As far as we know, it is the first time to apply the MOPSO to spacecraft trajectory optimization. The MOPSO proves to be quite effective in finding the Pareto-optimal solutions to asteroid rendezvous multiobjective optimization problems. The MOPSO shows a very competitive performance with respect to two highly competitive multiobjective evolutionary algorithms: the nondominated sorting genetic algorithm-II (NSGA-II) [17] and the strength Pareto evolutionary algorithm-II (SPEA-II) [18].


## 2. Optimization Problem

### 2.1. Asteroid Rendezvous Design Problem

The interplanetary transfer trajectory is always divided into three different segments, that is, planet departure segment, heliocentric transfer segment, and capture segment by using the concept of influencing sphere. As the influencing sphere of the planet is much smaller than that of the sun, the flight path and flight time of the departure and capture segments are much small compared with that of the heliocentric transfer. Therefore, for a preliminary trajectory design of an asteroid exploration mission, the design emphasis is firstly focused on the heliocentric transfer, while the planet departure is assumed as instantaneous process with an impulsive maneuver and the capture segment is omitted.This paper studies the asteroid rendezvous problem departure from the Earth.t0 represents the spacecraft launch date which is to be designed, and V∞, u, and v define the departure hyperbolic velocity according to the formulas
(1)θ=2πu,φ=arccos(2v-1)-π2,v∞V∞=cos⁡(φ)cos⁡(θ)i+cos⁡(φ)sin(θ)j+sin(φ)k.LetrE(t0), vE(t0) be the heliocentric position and velocity of the Earth at t0; then the initial position and velocity vector of the spacecraft, r(t0) and v(t0), entering the heliocentric transfer trajectory are defined as
(2)r(t0)=rE(t0),v(t0)=vE(t0)+v∞.The heliocentric transfer trajectory is modeled by a two-body dynamic model with the following governing equations:(3)x˙=vx,y˙=vy,z˙=vz,v˙x=-μx(x2+y2+z2)3/2+ax,v˙y=-μy(x2+y2+z2)3/2+ay,v˙z=-μz(x2+y2+z2)3/2+az,
where x, y, and z are position components of spacecraft, vx, vy, and vz are velocity components of spacecraft, ax, ay, and az are acceleration component of spacecraft, and μ is the gravity parameter of the sun.The thrust accelerationΓ(t)=(ax,ay,az)T can be approximated as n impulses
(4)Γ(t)=∑i=1nΔviδ(t-ti),
where ti is the time where an impulse is applied, Δvi=(Δvxi,Δvyi,Δvzi)T(i=1,2.…,n), and they are all the variables to be designed. The final state conditions are defined
(5)r(tf)=rA(tf),v(tf)=vA(tf),
where tf is the rendezvous time of spacecraft with the target asteroid and rA(tf) and vA(tf) are the position and velocity vector of the target asteroid at the rendezvous time.


### 2.2. Feasible-Solution Iteration Model

In order to avoid dealing with the equality constraints described as (5), the Lambert algorithm is employed to establish the infeasible iteration optimization model. The chosen independent variables, that is, the optimization variables, are impulse times and the first n-2 impulses
(6)ti(i=1,2,…,n),Δvi(i=1,2,…,n-2).The last two impulses are determined by solving the Lambert problem constrained by (5). In this feasible iteration model, each evaluation of the objective function produces a feasible solution that satisfies implicitly the rendezvous conditions. Detail on this multi-impulse rendezvous optimization model using the Lambert algorithm can be found in [19, 20].


### 2.3. Multiobjective Optimization Model

The total velocity characteristic is chosen as the first objective function(7)min⁡f1(x)=Δv=V∞+∑i=1n|Δvi|.The heliocentric transfer time is chosen as the second objective function(8)min⁡f2(x)=tf-t0.The constraint on the time of impulse is considered. The general constraint onti(i=1,2,…,n) is
(9)t0≤t1<t2⋯<tn≤tf.As the transfer time is one of the objective functions,tf is chosen as an optimization variable. The total optimization variables include four parts: t0, V∞, u, v, ti(i=1,2,…,n), Δvi(i=1,2,…,n-2), and tf. In order to improve optimization performance, the variable-scaling method is imposed on ti. Let αi=(ti-t0)/(tf-t0)αi≤1. Therefore, the final optimization variables x are
(10)x=(t0,V∞,u,v,α1,…,αn,Δv1,…,Δvn-2,tf)T.


## 2.1. Asteroid Rendezvous Design Problem

The interplanetary transfer trajectory is always divided into three different segments, that is, planet departure segment, heliocentric transfer segment, and capture segment by using the concept of influencing sphere. As the influencing sphere of the planet is much smaller than that of the sun, the flight path and flight time of the departure and capture segments are much small compared with that of the heliocentric transfer. Therefore, for a preliminary trajectory design of an asteroid exploration mission, the design emphasis is firstly focused on the heliocentric transfer, while the planet departure is assumed as instantaneous process with an impulsive maneuver and the capture segment is omitted.This paper studies the asteroid rendezvous problem departure from the Earth.t0 represents the spacecraft launch date which is to be designed, and V∞, u, and v define the departure hyperbolic velocity according to the formulas
(1)θ=2πu,φ=arccos(2v-1)-π2,v∞V∞=cos⁡(φ)cos⁡(θ)i+cos⁡(φ)sin(θ)j+sin(φ)k.LetrE(t0), vE(t0) be the heliocentric position and velocity of the Earth at t0; then the initial position and velocity vector of the spacecraft, r(t0) and v(t0), entering the heliocentric transfer trajectory are defined as
(2)r(t0)=rE(t0),v(t0)=vE(t0)+v∞.The heliocentric transfer trajectory is modeled by a two-body dynamic model with the following governing equations:(3)x˙=vx,y˙=vy,z˙=vz,v˙x=-μx(x2+y2+z2)3/2+ax,v˙y=-μy(x2+y2+z2)3/2+ay,v˙z=-μz(x2+y2+z2)3/2+az,
where x, y, and z are position components of spacecraft, vx, vy, and vz are velocity components of spacecraft, ax, ay, and az are acceleration component of spacecraft, and μ is the gravity parameter of the sun.The thrust accelerationΓ(t)=(ax,ay,az)T can be approximated as n impulses
(4)Γ(t)=∑i=1nΔviδ(t-ti),
where ti is the time where an impulse is applied, Δvi=(Δvxi,Δvyi,Δvzi)T(i=1,2.…,n), and they are all the variables to be designed. The final state conditions are defined
(5)r(tf)=rA(tf),v(tf)=vA(tf),
where tf is the rendezvous time of spacecraft with the target asteroid and rA(tf) and vA(tf) are the position and velocity vector of the target asteroid at the rendezvous time.


## 2.2. Feasible-Solution Iteration Model

In order to avoid dealing with the equality constraints described as (5), the Lambert algorithm is employed to establish the infeasible iteration optimization model. The chosen independent variables, that is, the optimization variables, are impulse times and the first n-2 impulses
(6)ti(i=1,2,…,n),Δvi(i=1,2,…,n-2).The last two impulses are determined by solving the Lambert problem constrained by (5). In this feasible iteration model, each evaluation of the objective function produces a feasible solution that satisfies implicitly the rendezvous conditions. Detail on this multi-impulse rendezvous optimization model using the Lambert algorithm can be found in [19, 20].


## 2.3. Multiobjective Optimization Model

The total velocity characteristic is chosen as the first objective function(7)min⁡f1(x)=Δv=V∞+∑i=1n|Δvi|.The heliocentric transfer time is chosen as the second objective function(8)min⁡f2(x)=tf-t0.The constraint on the time of impulse is considered. The general constraint onti(i=1,2,…,n) is
(9)t0≤t1<t2⋯<tn≤tf.As the transfer time is one of the objective functions,tf is chosen as an optimization variable. The total optimization variables include four parts: t0, V∞, u, v, ti(i=1,2,…,n), Δvi(i=1,2,…,n-2), and tf. In order to improve optimization performance, the variable-scaling method is imposed on ti. Let αi=(ti-t0)/(tf-t0)αi≤1. Therefore, the final optimization variables x are
(10)x=(t0,V∞,u,v,α1,…,αn,Δv1,…,Δvn-2,tf)T.


## 3. Optimization Algorithms

### 3.1. Introduction of MOEA

A general multiobjective optimization problem is to find the design variables that optimize a vector objective function over the feasible design space. The objective functions are the quantities that the designer wishes to minimize, maximize, or attain a certain value. The problem formulation in standard form for a minimization is given here, which is similar for the other cases:(11)Minimizef(x)=(f1(x),f2(x),…,fm(x))T,(12)subjecttogi(x)≤0,i=1,2,…,p,hj(x)=0,j=1,2,…,q,
where
(13)x=(x1,x2,…,xv)T∈X⊂Rv.For the multiobjective asteroid rendezvous design problem, the two objective functions are described by (7) and (8), the constraints are described by (9), and the optimization variables are described as (10).The classical optimization method for a multiobjective optimization problem is the weighting method. In recent years, the multiobjective evolutionary algorithms (MOEA) have been greatly investigated in the domain of multiobjective optimization. There are many variants of MOEA reported in the literature; a recent survey on MOEA and their application in aeronautical and aerospace engineering has been made in [16].In the study, except for the MOPSO, we also test two other mostly popular MOEA. The first is the NSGA-II algorithm which is proposed by Deb et al. [17]. This algorithm uses the idea of transforming the m objectives to a single fitness measure by the creation of a number of fronts, sorted according to nondomination. During the fitness assignment, the first front is created as the set of solutions that is not dominated by any solutions in the population. These solutions are given the highest fitness and temporarily removed from the population, then a second nondominated front consisting of the solutions that are now nondominated is built and assigned the second-highest fitness, and so forth. This is repeated until each of the solutions has been assigned a fitness. After each front has been created, its members are assigned crowding distances (normalized distance to closest neighbors in the front in the objective space) later to be used for niching. The NSGA-II has been successfully applied in spacecraft trajectory optimization, for example, in designing a three-objective impulse rendezvous problem [20, 21] and a two-objective robust rendezvous problem with considering uncertainty [22].The second is the SPEA-II proposed by Zitzler et al. [18]. It uses an archive containing nondominated solutions previously found (the so-called external nondominated set). At each generation, nondominated individuals are copied to the external nondominated set, removing the dominated solutions. For each individual in this external set, a strength value is computed. Pareto dominance is adopted to ensure that the solutions are properly distributed along the Pareto front. It also uses a nearest neighbor density estimation technique and a fine-grained fitness assignment strategy which guide the search more efficiently.


### 3.2. Brief Description of the MOPSO

The MOPSO applied in this study is the algorithm proposed by Pulido and Coello [23], which was competitive against the most popular MOEA such as the NSGA-II, the PAES, and other MOPSO on typical benchmark problems, under common performance metrics [23]. The source code of the MOPSO is available from the EMOO repository located at http://delta.cs.cinvestav.mx/~ccoello/EMOO.The MOPSO is based on the use of Pareto ranking and a subdivision of decision variable space into several subswarms which is done using clustering techniques. The complete execution process of this algorithm can be divided into three stages: initialization, flight, and generation of results [23].At the first stage, every swarm is initialized. Each swarm creates and initializes its own particles and generates the leaders set among the particle swarm set by using Pareto ranking. In the second stage, it firstly performs the execution of the flight of every swarm; next, it applies a clustering algorithm to group the guide particles. This is performed until reaching a total ofGMax iterations. The execution of the flight of each swarm can be seen as an entire PSO process (with the difference that it will only optimize a specific region of the search space). First, each particle will select a leader to which it will follow. At the same time, each particle will try to outperform its leader and to update its position. If the updated particle is not dominated by any member of the leaders set, then it will become a new leader. The execution of the swarm will start again until a total ofsgmax iterations are reached. Once all the swarms have finished theirs flights, a clustering algorithm takes the control by grouping the closest particle guides into nswarms swarms. Each resulting group will be assigned to a different swarm. The third and final stage will present all the nondominated solutions found.Details of this algorithm can be found in [23] and its pseudocode code is shown in Algorithm 1.Algorithm 1
Begin
Foreach swarm
(1) Initialize its particles
(2) Initialize the set of global leaders:gleader
End For
DO
Foreach swarm
Do
For each particle
(3) Select a leader
(4) Perform the flight
(5) Update the value
If it is a leader then add to thegleader
End For
While (number of iterations ≤ sgmax)
(6) Store leaders ingleader in nswarms  groups
End For
(7) Assign each leader group to a swarm
While(number of iterations ≤ GMax)
EndThe MOPSO algorithm requires the following parameters: (1)GMax: the total number of generations that the algorithm will be executed; (2) sgmax: the number of internal generations that the particles of each swarm will run before sharing their leaders; (3)  nparticles: the total number of particles; and (4)  nswarms: the number of particle groups.


### 3.3. Constraints Optimization Method

The multiobjective asteroid mission design problem is a highly constrained problem whose constraints are described in (9). The simulation experiments show that these constraints strongly influence the convergence. In the MOPSO, the constraints are always handled in checking Pareto dominance [23, 24]. When we compare two individuals, we first check their feasibility. If they are both feasible, then the comparison is done using Pareto dominance. If one is feasible and the other is infeasible, the feasible individual wins. If both are infeasible, then the individual with the lower amount of total constraint violation wins.The total constraint is calculated by making use of a nondifferentiable penalty function. For the general constrained problem in (10)–(13), the penalized, total constraint function C(x,M) is
(14)C(x,M)=M[∑i=1pmax⁡(0,gi(x))+∑j=1q|hj(x)|].In the present workM=1000000 is used as a penalty coefficient.For most optimization problems, each design variable has its own upper and low values. Thus, a strategy to maintain the particles within the search space in case they go beyond their boundaries is necessary for the MOPSO algorithm.In [23, 24], the following strategy was employed. When a decision variable goes beyond its boundaries, the decision variable takes the value of its corresponding boundary (either the lower or the upper boundary). The strategy will be likely effective when the Pareto solutions are located in the boundary of variables.However, our simulation experiments show that this strategy is not very effective in solving our problems. Therefore, another simple strategy is proposed. When a decision variable goes beyond its boundaries, the decision variable takes a random value from its feasible design space. The probability-based disposal to boundary constraint could enrich the diversity of swarm flight.The MOPSO with these two different boundary constraint optimization strategies is, respectively, called as MOPSO-I and MOPSO-II. Their performance is compared in Section4.3.


### 3.4. Algorithm Assessment Metrics

In order to allow a quantitative assessment of the performance of the MSOPO, we adopted the following two metrics.The first one is the epsilon indicator [25]. Given a reference Pareto front (ideally the true Pareto front, if available), the epsilon indicator measures the minimum amount ε necessary to translate all the points of the found Pareto front to weakly dominate the reference set. The epsilon indicator has two types, that is, the additive type and the multiplicative type, and the multiplicative type is used here.The second one is the hypervolume indicator [26]. The hypervolume of a set of solutions measures the size of the portion of objective space that is dominated by those solutions collectively. Generally, hypervolume is favored because it captures in a single scalar both the closeness of the solutions to the optimal set and, to some extent, the spread of the solutions across objective space.The statistical test chosen for result evaluation is the Mann-Whitney test [27]. This is a nonparametric rank-based test that can be used to compare two independent sets of sampled data. It outputs P values that estimate the probability of a failure to reject the null hypothesis of the study question. Here the (1-P) values can be interpreted as the probability that the performance of one algorithm is different (superior or junior) with statistical significance to that of the other.


## 3.1. Introduction of MOEA

A general multiobjective optimization problem is to find the design variables that optimize a vector objective function over the feasible design space. The objective functions are the quantities that the designer wishes to minimize, maximize, or attain a certain value. The problem formulation in standard form for a minimization is given here, which is similar for the other cases:(11)Minimizef(x)=(f1(x),f2(x),…,fm(x))T,(12)subjecttogi(x)≤0,i=1,2,…,p,hj(x)=0,j=1,2,…,q,
where
(13)x=(x1,x2,…,xv)T∈X⊂Rv.For the multiobjective asteroid rendezvous design problem, the two objective functions are described by (7) and (8), the constraints are described by (9), and the optimization variables are described as (10).The classical optimization method for a multiobjective optimization problem is the weighting method. In recent years, the multiobjective evolutionary algorithms (MOEA) have been greatly investigated in the domain of multiobjective optimization. There are many variants of MOEA reported in the literature; a recent survey on MOEA and their application in aeronautical and aerospace engineering has been made in [16].In the study, except for the MOPSO, we also test two other mostly popular MOEA. The first is the NSGA-II algorithm which is proposed by Deb et al. [17]. This algorithm uses the idea of transforming the m objectives to a single fitness measure by the creation of a number of fronts, sorted according to nondomination. During the fitness assignment, the first front is created as the set of solutions that is not dominated by any solutions in the population. These solutions are given the highest fitness and temporarily removed from the population, then a second nondominated front consisting of the solutions that are now nondominated is built and assigned the second-highest fitness, and so forth. This is repeated until each of the solutions has been assigned a fitness. After each front has been created, its members are assigned crowding distances (normalized distance to closest neighbors in the front in the objective space) later to be used for niching. The NSGA-II has been successfully applied in spacecraft trajectory optimization, for example, in designing a three-objective impulse rendezvous problem [20, 21] and a two-objective robust rendezvous problem with considering uncertainty [22].The second is the SPEA-II proposed by Zitzler et al. [18]. It uses an archive containing nondominated solutions previously found (the so-called external nondominated set). At each generation, nondominated individuals are copied to the external nondominated set, removing the dominated solutions. For each individual in this external set, a strength value is computed. Pareto dominance is adopted to ensure that the solutions are properly distributed along the Pareto front. It also uses a nearest neighbor density estimation technique and a fine-grained fitness assignment strategy which guide the search more efficiently.


## 3.2. Brief Description of the MOPSO

The MOPSO applied in this study is the algorithm proposed by Pulido and Coello [23], which was competitive against the most popular MOEA such as the NSGA-II, the PAES, and other MOPSO on typical benchmark problems, under common performance metrics [23]. The source code of the MOPSO is available from the EMOO repository located at http://delta.cs.cinvestav.mx/~ccoello/EMOO.The MOPSO is based on the use of Pareto ranking and a subdivision of decision variable space into several subswarms which is done using clustering techniques. The complete execution process of this algorithm can be divided into three stages: initialization, flight, and generation of results [23].At the first stage, every swarm is initialized. Each swarm creates and initializes its own particles and generates the leaders set among the particle swarm set by using Pareto ranking. In the second stage, it firstly performs the execution of the flight of every swarm; next, it applies a clustering algorithm to group the guide particles. This is performed until reaching a total ofGMax iterations. The execution of the flight of each swarm can be seen as an entire PSO process (with the difference that it will only optimize a specific region of the search space). First, each particle will select a leader to which it will follow. At the same time, each particle will try to outperform its leader and to update its position. If the updated particle is not dominated by any member of the leaders set, then it will become a new leader. The execution of the swarm will start again until a total ofsgmax iterations are reached. Once all the swarms have finished theirs flights, a clustering algorithm takes the control by grouping the closest particle guides into nswarms swarms. Each resulting group will be assigned to a different swarm. The third and final stage will present all the nondominated solutions found.Details of this algorithm can be found in [23] and its pseudocode code is shown in Algorithm 1.Algorithm 1
Begin
Foreach swarm
(1) Initialize its particles
(2) Initialize the set of global leaders:gleader
End For
DO
Foreach swarm
Do
For each particle
(3) Select a leader
(4) Perform the flight
(5) Update the value
If it is a leader then add to thegleader
End For
While (number of iterations ≤ sgmax)
(6) Store leaders ingleader in nswarms  groups
End For
(7) Assign each leader group to a swarm
While(number of iterations ≤ GMax)
EndThe MOPSO algorithm requires the following parameters: (1)GMax: the total number of generations that the algorithm will be executed; (2) sgmax: the number of internal generations that the particles of each swarm will run before sharing their leaders; (3)  nparticles: the total number of particles; and (4)  nswarms: the number of particle groups.


## 3.3. Constraints Optimization Method

The multiobjective asteroid mission design problem is a highly constrained problem whose constraints are described in (9). The simulation experiments show that these constraints strongly influence the convergence. In the MOPSO, the constraints are always handled in checking Pareto dominance [23, 24]. When we compare two individuals, we first check their feasibility. If they are both feasible, then the comparison is done using Pareto dominance. If one is feasible and the other is infeasible, the feasible individual wins. If both are infeasible, then the individual with the lower amount of total constraint violation wins.The total constraint is calculated by making use of a nondifferentiable penalty function. For the general constrained problem in (10)–(13), the penalized, total constraint function C(x,M) is
(14)C(x,M)=M[∑i=1pmax⁡(0,gi(x))+∑j=1q|hj(x)|].In the present workM=1000000 is used as a penalty coefficient.For most optimization problems, each design variable has its own upper and low values. Thus, a strategy to maintain the particles within the search space in case they go beyond their boundaries is necessary for the MOPSO algorithm.In [23, 24], the following strategy was employed. When a decision variable goes beyond its boundaries, the decision variable takes the value of its corresponding boundary (either the lower or the upper boundary). The strategy will be likely effective when the Pareto solutions are located in the boundary of variables.However, our simulation experiments show that this strategy is not very effective in solving our problems. Therefore, another simple strategy is proposed. When a decision variable goes beyond its boundaries, the decision variable takes a random value from its feasible design space. The probability-based disposal to boundary constraint could enrich the diversity of swarm flight.The MOPSO with these two different boundary constraint optimization strategies is, respectively, called as MOPSO-I and MOPSO-II. Their performance is compared in Section4.3.


## 3.4. Algorithm Assessment Metrics

In order to allow a quantitative assessment of the performance of the MSOPO, we adopted the following two metrics.The first one is the epsilon indicator [25]. Given a reference Pareto front (ideally the true Pareto front, if available), the epsilon indicator measures the minimum amount ε necessary to translate all the points of the found Pareto front to weakly dominate the reference set. The epsilon indicator has two types, that is, the additive type and the multiplicative type, and the multiplicative type is used here.The second one is the hypervolume indicator [26]. The hypervolume of a set of solutions measures the size of the portion of objective space that is dominated by those solutions collectively. Generally, hypervolume is favored because it captures in a single scalar both the closeness of the solutions to the optimal set and, to some extent, the spread of the solutions across objective space.The statistical test chosen for result evaluation is the Mann-Whitney test [27]. This is a nonparametric rank-based test that can be used to compare two independent sets of sampled data. It outputs P values that estimate the probability of a failure to reject the null hypothesis of the study question. Here the (1-P) values can be interpreted as the probability that the performance of one algorithm is different (superior or junior) with statistical significance to that of the other.


## 4. Simulation Results

### 4.1. Problem Configuration

In order to testify the effectiveness of the proposed method, two different asteroid mission designs are illustrated. Table1 lists the orbit elements of the two asteroids. The upper and lower space of the optimization variables are provided in Table 2. Three different cases with the number of impulses of 2, 3, and 4 are, respectively, tested for each asteroid mission.Table 1
Orbit elements of asteroid.
Index
Asteroid 1
Asteroid 2
Name
1999YR14
2340
a (AU)
1.65365126892224
0.84421076388332195
e
0.40069261757759106
0.44975834146342486
i (deg)
3.7221930161441943
5.8547882390182853
Ω (deg)
3.1338963493744654
211.50460158030430
ω (deg)
9.4143875285008676
39.994195753797953
M (deg)
114.73402134869427
240.44827444641544
Epoch (MJD2000)
3255.0000
3255.000000Table 2
Design space of variables.
Variables
Space
Units
t
0
[
4000,10000
]
MJD2000
V
∞
[
0,5
]
km/s
u
[
0,1
]
n/a
v
[
0,1
]
n/a
α
i
(
i
=
1,2
,
…
,
n
)
n
=
2: [0,0.6],[0.3,1]
n/a
n
=
3: [0,0.4],[0.1,0.8],[0.5,1]
n
=
4: [0,0.4],[0.1,0.8],[0.5,0.9],[0.7,1]
Δ
v
i
x
,
Δ
v
i
y
,
Δ
v
i
z  (i=1,2,…,n-2)
[
-
4,4
]
km/s
t
f
-
t
0
[
100,1500
]
days


### 4.2. Pareto Fronts Analysis

From our experiments, this multiobjective asteroid rendezvous design problem is very difficult to be solved, and the true Pareto fronts of this multiobjective problem are not known. In order to obtain the Pareto fronts as possible close to the true ones, the MOPSO is executed with a much larger number of function evaluations. The parameters of the MOPSO arenparticles=400, GMax=400, sgmax=5, and nswarms=20.Considering the stochastic characteristic of the MOPSO, 10 independent runs for each test case are completed. All the Pareto solutions of the 10 independent runs are compared, and the repeated and non-Pareto solutions are deleted, and the revised Pareto solutions are selected as the final solutions. In the following examples, the figured Pareto fronts are all obtained in the same method. The MOPSO with different boundary constraint optimization methods, that is, the MOPSO-I and MOPSO-II, is both tested. The performance comparisons between the MOPSO-I and MOPSO-II will be analyzed in the next section.Figure1 compares the Pareto solution front for two-impulse, three-impulse, and four-impulse of the Asteroid 1999YR14 mission. Figure 2 compares this for the Asteroid 2340 mission. The tradeoffs between the total Δv and the transfer time are clearly demonstrated by Figures 1 and 2, which will be useful for a mission designer. Several inherent peculiarities regarding the optimal multiobjective asteroid rendezvous trajectories have been observed from the optimization results. For the second asteroid, the corresponding total Δv reduces evidently as the transfer time increases when the transfer time is in the range of [100,140] day, but small change in other ranges; therefore the Pareto front concentrates on a narrow range of 40 days. For the first asteroid, the corresponding total Δv reduces evidently as the transfer time increases in a large range of [100,900] day; therefore the Pareto solutions distribute in a much larger space, but the Pareto front is discontinuous in small range as demonstrated in Figure 1.Pareto fronts of Asteroid 1.
(a)
Produced by MOPSO-I
(b)
Produced by MOPSO-IIPareto fronts of Asteroid 2.
(a)
Produced by MOPSO-I
(b)
Produced by MOPSO-IIIn order to explain this phenomenon, the propellant-optimal solutions are obtained by using the approach employed in [28]. In this optimization, the transfer time is fixed and Δv is calculated every 5 days in the transfer time of range [100,1100] days. The relations between the transfer time and Δv of the propellant-optimal solutions are illustrated in Figure 3. Some transfer time ranges with much-higher propellant cost are located; for example, the Δv is about 9.5 km/s for the transfer time of 300 days, while being only 6.2 km/s for the transfer time of 280 days. The former solution is larger than the latter solution by about 50% in Δv; besides its transfer time is also larger than the latter. This issue can explain why the Pareto fronts of two-impulse solutions discontinue near the transfer time of 300 days.Figure 3
Relations between transfer time and total characteristic velocity (propellant-optimal solutions, Asteroid 1).The influence of different number of impulses on the totalΔv has been also demonstrated by Figures 1 and 2. It is obvious that a three-impulse trajectory and a four-impulse trajectory cost less propellant in compared with a two-impulse trajectory. It is necessary to increase the number of impulse to find a much better trajectory in the preliminary asteroid rendezvous mission design.


### 4.3. Performance Analysis of the MOPSO

#### 4.3.1. Boundary Constraint Optimization Comparisons

As seen from Figures1 and 2, the boundary constraint optimization method of the MOPSO has a great influence on the Pareto fronts. Graphically, the MOPSO-II has much better performance in diversity, and the spread of its solutions found is much larger than that of the MOPSO-I. This can easily explain that probability-based disposal to boundary constraint can enrich the diversity of swarm flight.It is not easily to determine which one is much closer to the true Pareto fronts from Figure1, and the two-impulse comparison case is shown in Figure 4(a), while it is easy to see from Figure 2 that the Pareto fronts obtained by MOPSO-I are much closer to the true Pareto fronts, and the three-impulse comparison case is redrawn in Figure 4(b) to demonstrate this clearly.Pareto fronts produced by MOPSO with different boundary constraint optimization.
(a)
Asteroid 1, two-impulse
(b)
Asteroid 2, three-impulseTherefore, the quantitative metrics are calculated and the statistical results are provided in Table3. As seen from Table 3, the MOPSO-II has better performance with respect to the hybervolume indicator (larger value means better) for the two cases. While, for the epsilon indicator (smaller value means better), the MOPSO-II is much better in case of Asteroid 1, two-impulse problem and the MOPSO-I are much better in case of Asteroid 2, three-impulse.Table 3
MOPSO performance with different boundary constraint optimization.
Test problem
Algorithm
Epsilon
Hybervolume
Mean
std
P value
Mean
std
P value
Asteroid 1,two-impulse
MOPSO-I
0.011987
0.008378
0.1620
0.972618
0.006829
0.0257
MOPSO-II
0.006833
0.004066
0.977958
0.001301
Asteroid 2,three-impulse
MOPSO-I
0.025630
0.0278051
0.1620
0.87509
0.308357
0.6232
MOPSO-II
0.0399274
0.026254
0.968537
0.016548Through the comparisons provided in Table2 and Figure 4, our proposed boundary constraint demonstrates much better performance in total. Besides, these comparisons show that the same algorithm will have different performance in solving the same type of problem with different configurations.


#### 4.3.2. Algorithm Parameters Analysis

Our simulation experiments show thatnparticles and GMax are the two main parameters affecting the performance of the MOPSO in solving this complex multiobjective problem. In order to quantitatively evaluate their influence o, nparticles is, respectively, set as 400, 200, and 100, GMax is also, respectively, set as 400, 200, and 100, and a total of nine groups of parameters (in Table 4) are furthered tested. Other parameters are the same as used in Section 4.2.Table 4
MOPSO parameters configurations.
Index
Case  1
Case  2
Case  3
Case  4
Case  5
Case  6
Case  7
Case  8
Case  9
n
particles
400
400
400
200
200
200
100
100
100
GMax
400
200
100
400
200
100
400
200
100Ten independent runs for the MOPSO-II with each group of parameters in solving the three-impulse, Asteroid 1 problem, are executed.Figure5(a) compares the Pareto fronts of Cases 1, 5, and 9, Figure 5(b) compares that of Cases 2, 5, and 8, and Figure 5(c) compares that of Cases 7, 8, and 9. These comparisons are to show the influence of the total number of function evaluations. It is clearly seen from Figure 5 that the larger number of function evaluations evidently increases the performance of the MOPSO. Thus, a larger number of swarm size and iteration are necessary for obtaining the optimal Pareto fronts for this practical multiobjective optimization problem.Pareto fronts produced by MOPSO-II with different number of function evaluations.
(a)
Cases 1, 5, and 9
(b)
Cases 2, 5, and 8
(c)
Cases 7, 8, and 9Figure6(a) compares the Pareto fronts of Cases 2 and 4, Figure 6(b) compares that of Cases 6 and 8, and Figure 6(c) compares that of Cases 3, 5, and 7. These comparisons are to show the influence of nparticles and GMax with the same number of function evaluations. Evidently, we cannot determine which case is better by the graphical results provided in Figure 6. Therefore, the quantitative metrics are calculated and the statistical results are provided in Tables 5, 6, 7, and 8.Table 5
MOPSO-II performance with different parameters (Cases  2 and  4).
Test problem
Algorithm
Epsilon
Hybervolume
Mean
std
P value
Mean
std
P value
Asteroid 1,three-impulse
Case  2
0.010692
0.006797
0.6776
0.984063
0.002738
0.1405
Case  4
0.012918
0.012848
0.985643
0.002736Table 6
MOPSO-II performance with different parameters (Cases  6 and  8).
Test problem
Algorithm
Epsilon
Hybervolume
Mean
std
P value
Mean
std
P value
Asteroid 1,three-impulse
Case  6
0.006045
0.002645
0.4727
0.99544
0.000816
0.0757
Case  8
0.007806
0.004042
0.994456
0.001074Table 7
MOPSO-II performance with different parameters (Cases  3,  5, and  7).
Test problem
Algorithm
Epsilon
Hybervolume
Mean
std
Mean
std
Asteroid 1,three-impulse
Case  3
0.00583413
0.0036399
0.994062
0.00135846
Case  5
0.00438474
0.00129004
0.994755
0.00100777
Case  7
0.00564449
0.00274183
0.994917
0.00119178Table 8
P value of MOPSO (Cases  3,  5, and  7).
Case  3
Case  5
Case  7
Case  3
0.4727
0.9097
Case  5
0.3447
0.4727
Case  7
0.2123
0.7337
Bold indicates epsilon; the other is hybervolume.Pareto fronts produced by MOPSO-II with same number of function evaluations.
(a)
Cases 2, 4
(b)
Cases 6, 8
(c)
Cases 3, 5, and 7The comparisons between Case 2 and Case 4 show that larger size of swarm can obtain better average epsilon indicator with aP value of 67.76% but worse average hybervolume indicator with a P value of 14.05%. The comparisons between Case 6 and Case 8 show that larger size of swarm both increase average performance in epsilon and hybervolume, with a P value of 47.27% and 7.57%. By comparing Cases 3, 5, and 7, Case 5 demonstrates the best performance in total. However, the improvements are not evident as the P values are all a large value.Although the result obtained from this experiment seems to be inconclusive on which parameter is the best, we can argue that the MOPSO’s performance is not sensitive to the main PSO parameters, under the condition that the total number of function evaluation retains a large value. We have obtained competitive results in most cases without paying special attentions on PSO parameters.


### 4.4. Comparison of the MOPSO with Other MOEA

To demonstrate the performance of the employed MOPSO, two other popular algorithms are tested and compared: the NSGA-II and the SPEA-II.In the following examples, the total number of function evaluations was set to 20000 for all the algorithms compared. The NSGA-II and the SPEA-II were run using a population size of 200, a maximum number of generations of 100. The MOPSO-I and MOPSO-II used 200 particles, a maximum number of generations of 20, a maximum number of generations per swarm of 5, and a total of 5 swarms. The source code of the NSGA-II, and SPEA-II provided in the EMOO repository is also employed in this study.The two test problems are the same to those employed in Section4.3 to compare boundary constraint optimization. Ten independent runs for each test case are executed. The final Pareto fronts are compared in Figure 7 and the statistical results with respect to epsilon and hybervolume indicators are provided in Tables 9, 10, and 11.Table 9
Performance comparisons of MOPSO, NSGA-II, and SPEA-II.
Test problem
Algorithm
Epsilon
Hybervolume
Mean
std
Mean
std
Asteroid 1,two-impulse
MOPSO-I
0.035041
0.038848
0.876849
0.026781
NSGA-II
0.113438
0.094049
0.821205
0.080834
SPEA-II
0.100881
0.032094
0.811347
0.035008
Asteroid 2,three-impulse
MOPSO-II
0.159179
0.081796
0.875023
0.078886
NSGA-II
0.0594614
0.058546
0.838198
0.301149
SPEA-II
0.12172
0.055917
0.902724
0.072393Table 10
P value of MOPSO-I, NSGA-II, and SPEA-II.
MOPSO-I
NSGA-II
SPEA-II
MOPSO-I
0.0058
0.0113
NSGA-II
0.0073
0.5205
SPEA-II
0.0017
0.1620
Bold indicates epsilon; the other is hybervolume.Table 11
P value of MOPSO-II, NSGA-II, and SPEA-II.
MOPSO-II
NSGA-II
SPEA-II
MOPSO-II
0.0113
0.3847
NSGA-II
0.3447
0.0211
SPEA-II
0.4274
0.5205
Bold indicates epsilon; the other is hybervolume.Pareto fronts produced by MOPSO, NSGA-II, and SPEA-II.
(a)
Asteroid 1, two-impulse
(b)
Asteroid 2, three-impulseEvidently seen from Figure7(a), the MOPSO-I produces better Pareto fronts than the NSGA-II and SPEA-II for the first test problem. The statistical results of the two indicators also support this. The P value with respect to epsilon is calculated as 0.58% and 1.13%, respectively, and 0.73% and 0.17% for hybervolume. This shows that the MOPSO-I is absolutely super to the NSGA-II and SPEA-II.Also evidently seen from Figure7(b), the MOPSO-I produces much closer Pareto fronts than the NSGA-II and SPEA-II for the second test problem. However, its Pareto fronts are a little narrower, which may result in that its epsilon and hybervolume indicator are junior to the NSGA-II and SPEA-II. Furthermore, the fact that the Pareto fronts are highly discontinuous would make these metrics irrelevant. In the view of practical spacecraft mission design, obtaining much closer Pareto fronts from which the designer can choose one single satisfying solution for engineering design would be much desirable performance for a multiobjective optimizer. Thus, we can say that the MOPSO-I is a little better than the NSGA-II and SPEA-II; even the statistical results on quantitative metrics do not support this.In conclusion, the MOPSO showed a better performance with respect to the NSGA-II and the SPEA-II for the two test problems. For other test problems, the MOPSO is not always better, but its performance is very competitive.


### 4.5. Comparison with PCP Method

The widely used Pork-Chop (PCP) method is also tested here for comparison. TheΔv is calculated every 1 day in both Earth departure time and asteroid arrival time, and the two-impulse Lambert algorithm is employed. The contours of time of flight corresponding to the optimal transfers for Asteroid 1 are illustrated in Figure 8. For convenience, only the solutions with a Δv less than 18 km/s and Earth departure and asteroid arrive times in [5000,7000] (MJD2000) are presented.Relations between total characteristic velocity, transfer time, and earth departure time.
(a)
Asteroid 1, two-impulse
(b)
Asteroid 2, two-impulseWe analyze the comparisons between the proposed multiobjective optimization approach and the PCP method in the following three aspects.


#### 4.5.1. Uses of the Results

In the PCP method, the contours of theΔv are plotted to assist the designer to find the best launch window and transfer trajectory. However, the contours have many local peaks and it is not very convenient to determine which solution is the best one. Besides, the size of the region containing the global minimum is a very important factor in finding real trajectories in any given year. This is not easily observed from the Pork-Chop figures.In contrast to the PCP method, the proposed multiobjective optimization approach can provide friendly the designer of this information. The relations between total characteristic velocity, transfer time, and earth departure time of those obtained Pareto-optimal solutions for two test cases are provided in Figure8. From Figure 8, it can be clearly observed that the earth departure time for all the optimal rendezvous trajectories concentrates on one arrow domain for these two cases. The optimized Earth launch window and its size can be easily determined through the multiobjective optimization design, which will provide very useful information for engineering design.


#### 4.5.2. Global Optimality

The minimum-propellant solution searched by the PCP method for the first asteroid mission is calculated with aΔv of 6015.3 m/s and it is located in Figure 9 (* denotes). As seen in Figure 1, the MOPSO has located a set of Pareto solutions with a Δv about this value. This issue can also demonstrate the global convergence ability of the MOPSO. Besides, the PCP method obtains only the two-impulse trajectory. As is well known, the two-impulse trajectory is not the propellant-optimal solution under most conditions, and increasing the number of impulses will reduce the propellant cost. Our proposed approach can design the multi-impulse trajectory; thus it locates better solution than the PCP method.Figure 9
Contours of time of flight corresponding to the optimal transfers (Asteroid 1).


#### 4.5.3. Efficiency

The PCP is in essence of an exhaustive searching method. In our test, the search space for departure time and arrive time is[4000,10000] day. With the search step of 1 day, the PCP method calculates a total number of 6000*6000 trajectories, while only 400*400*5 trajectories are calculated in the proposed method for the most highly cost case. The proposed approach improves the calculation efficiency by about 40 times. Besides, the PCP method needs much human intention to determine the final solutions, while the proposed approach is an automated search method.


## 4.1. Problem Configuration

In order to testify the effectiveness of the proposed method, two different asteroid mission designs are illustrated. Table1 lists the orbit elements of the two asteroids. The upper and lower space of the optimization variables are provided in Table 2. Three different cases with the number of impulses of 2, 3, and 4 are, respectively, tested for each asteroid mission.Table 1
Orbit elements of asteroid.
Index
Asteroid 1
Asteroid 2
Name
1999YR14
2340
a (AU)
1.65365126892224
0.84421076388332195
e
0.40069261757759106
0.44975834146342486
i (deg)
3.7221930161441943
5.8547882390182853
Ω (deg)
3.1338963493744654
211.50460158030430
ω (deg)
9.4143875285008676
39.994195753797953
M (deg)
114.73402134869427
240.44827444641544
Epoch (MJD2000)
3255.0000
3255.000000Table 2
Design space of variables.
Variables
Space
Units
t
0
[
4000,10000
]
MJD2000
V
∞
[
0,5
]
km/s
u
[
0,1
]
n/a
v
[
0,1
]
n/a
α
i
(
i
=
1,2
,
…
,
n
)
n
=
2: [0,0.6],[0.3,1]
n/a
n
=
3: [0,0.4],[0.1,0.8],[0.5,1]
n
=
4: [0,0.4],[0.1,0.8],[0.5,0.9],[0.7,1]
Δ
v
i
x
,
Δ
v
i
y
,
Δ
v
i
z  (i=1,2,…,n-2)
[
-
4,4
]
km/s
t
f
-
t
0
[
100,1500
]
days


## 4.2. Pareto Fronts Analysis

From our experiments, this multiobjective asteroid rendezvous design problem is very difficult to be solved, and the true Pareto fronts of this multiobjective problem are not known. In order to obtain the Pareto fronts as possible close to the true ones, the MOPSO is executed with a much larger number of function evaluations. The parameters of the MOPSO arenparticles=400, GMax=400, sgmax=5, and nswarms=20.Considering the stochastic characteristic of the MOPSO, 10 independent runs for each test case are completed. All the Pareto solutions of the 10 independent runs are compared, and the repeated and non-Pareto solutions are deleted, and the revised Pareto solutions are selected as the final solutions. In the following examples, the figured Pareto fronts are all obtained in the same method. The MOPSO with different boundary constraint optimization methods, that is, the MOPSO-I and MOPSO-II, is both tested. The performance comparisons between the MOPSO-I and MOPSO-II will be analyzed in the next section.Figure1 compares the Pareto solution front for two-impulse, three-impulse, and four-impulse of the Asteroid 1999YR14 mission. Figure 2 compares this for the Asteroid 2340 mission. The tradeoffs between the total Δv and the transfer time are clearly demonstrated by Figures 1 and 2, which will be useful for a mission designer. Several inherent peculiarities regarding the optimal multiobjective asteroid rendezvous trajectories have been observed from the optimization results. For the second asteroid, the corresponding total Δv reduces evidently as the transfer time increases when the transfer time is in the range of [100,140] day, but small change in other ranges; therefore the Pareto front concentrates on a narrow range of 40 days. For the first asteroid, the corresponding total Δv reduces evidently as the transfer time increases in a large range of [100,900] day; therefore the Pareto solutions distribute in a much larger space, but the Pareto front is discontinuous in small range as demonstrated in Figure 1.Pareto fronts of Asteroid 1.
(a)
Produced by MOPSO-I
(b)
Produced by MOPSO-IIPareto fronts of Asteroid 2.
(a)
Produced by MOPSO-I
(b)
Produced by MOPSO-IIIn order to explain this phenomenon, the propellant-optimal solutions are obtained by using the approach employed in [28]. In this optimization, the transfer time is fixed and Δv is calculated every 5 days in the transfer time of range [100,1100] days. The relations between the transfer time and Δv of the propellant-optimal solutions are illustrated in Figure 3. Some transfer time ranges with much-higher propellant cost are located; for example, the Δv is about 9.5 km/s for the transfer time of 300 days, while being only 6.2 km/s for the transfer time of 280 days. The former solution is larger than the latter solution by about 50% in Δv; besides its transfer time is also larger than the latter. This issue can explain why the Pareto fronts of two-impulse solutions discontinue near the transfer time of 300 days.Figure 3
Relations between transfer time and total characteristic velocity (propellant-optimal solutions, Asteroid 1).The influence of different number of impulses on the totalΔv has been also demonstrated by Figures 1 and 2. It is obvious that a three-impulse trajectory and a four-impulse trajectory cost less propellant in compared with a two-impulse trajectory. It is necessary to increase the number of impulse to find a much better trajectory in the preliminary asteroid rendezvous mission design.


## 4.3. Performance Analysis of the MOPSO

### 4.3.1. Boundary Constraint Optimization Comparisons

As seen from Figures1 and 2, the boundary constraint optimization method of the MOPSO has a great influence on the Pareto fronts. Graphically, the MOPSO-II has much better performance in diversity, and the spread of its solutions found is much larger than that of the MOPSO-I. This can easily explain that probability-based disposal to boundary constraint can enrich the diversity of swarm flight.It is not easily to determine which one is much closer to the true Pareto fronts from Figure1, and the two-impulse comparison case is shown in Figure 4(a), while it is easy to see from Figure 2 that the Pareto fronts obtained by MOPSO-I are much closer to the true Pareto fronts, and the three-impulse comparison case is redrawn in Figure 4(b) to demonstrate this clearly.Pareto fronts produced by MOPSO with different boundary constraint optimization.
(a)
Asteroid 1, two-impulse
(b)
Asteroid 2, three-impulseTherefore, the quantitative metrics are calculated and the statistical results are provided in Table3. As seen from Table 3, the MOPSO-II has better performance with respect to the hybervolume indicator (larger value means better) for the two cases. While, for the epsilon indicator (smaller value means better), the MOPSO-II is much better in case of Asteroid 1, two-impulse problem and the MOPSO-I are much better in case of Asteroid 2, three-impulse.Table 3
MOPSO performance with different boundary constraint optimization.
Test problem
Algorithm
Epsilon
Hybervolume
Mean
std
P value
Mean
std
P value
Asteroid 1,two-impulse
MOPSO-I
0.011987
0.008378
0.1620
0.972618
0.006829
0.0257
MOPSO-II
0.006833
0.004066
0.977958
0.001301
Asteroid 2,three-impulse
MOPSO-I
0.025630
0.0278051
0.1620
0.87509
0.308357
0.6232
MOPSO-II
0.0399274
0.026254
0.968537
0.016548Through the comparisons provided in Table2 and Figure 4, our proposed boundary constraint demonstrates much better performance in total. Besides, these comparisons show that the same algorithm will have different performance in solving the same type of problem with different configurations.


### 4.3.2. Algorithm Parameters Analysis

Our simulation experiments show thatnparticles and GMax are the two main parameters affecting the performance of the MOPSO in solving this complex multiobjective problem. In order to quantitatively evaluate their influence o, nparticles is, respectively, set as 400, 200, and 100, GMax is also, respectively, set as 400, 200, and 100, and a total of nine groups of parameters (in Table 4) are furthered tested. Other parameters are the same as used in Section 4.2.Table 4
MOPSO parameters configurations.
Index
Case  1
Case  2
Case  3
Case  4
Case  5
Case  6
Case  7
Case  8
Case  9
n
particles
400
400
400
200
200
200
100
100
100
GMax
400
200
100
400
200
100
400
200
100Ten independent runs for the MOPSO-II with each group of parameters in solving the three-impulse, Asteroid 1 problem, are executed.Figure5(a) compares the Pareto fronts of Cases 1, 5, and 9, Figure 5(b) compares that of Cases 2, 5, and 8, and Figure 5(c) compares that of Cases 7, 8, and 9. These comparisons are to show the influence of the total number of function evaluations. It is clearly seen from Figure 5 that the larger number of function evaluations evidently increases the performance of the MOPSO. Thus, a larger number of swarm size and iteration are necessary for obtaining the optimal Pareto fronts for this practical multiobjective optimization problem.Pareto fronts produced by MOPSO-II with different number of function evaluations.
(a)
Cases 1, 5, and 9
(b)
Cases 2, 5, and 8
(c)
Cases 7, 8, and 9Figure6(a) compares the Pareto fronts of Cases 2 and 4, Figure 6(b) compares that of Cases 6 and 8, and Figure 6(c) compares that of Cases 3, 5, and 7. These comparisons are to show the influence of nparticles and GMax with the same number of function evaluations. Evidently, we cannot determine which case is better by the graphical results provided in Figure 6. Therefore, the quantitative metrics are calculated and the statistical results are provided in Tables 5, 6, 7, and 8.Table 5
MOPSO-II performance with different parameters (Cases  2 and  4).
Test problem
Algorithm
Epsilon
Hybervolume
Mean
std
P value
Mean
std
P value
Asteroid 1,three-impulse
Case  2
0.010692
0.006797
0.6776
0.984063
0.002738
0.1405
Case  4
0.012918
0.012848
0.985643
0.002736Table 6
MOPSO-II performance with different parameters (Cases  6 and  8).
Test problem
Algorithm
Epsilon
Hybervolume
Mean
std
P value
Mean
std
P value
Asteroid 1,three-impulse
Case  6
0.006045
0.002645
0.4727
0.99544
0.000816
0.0757
Case  8
0.007806
0.004042
0.994456
0.001074Table 7
MOPSO-II performance with different parameters (Cases  3,  5, and  7).
Test problem
Algorithm
Epsilon
Hybervolume
Mean
std
Mean
std
Asteroid 1,three-impulse
Case  3
0.00583413
0.0036399
0.994062
0.00135846
Case  5
0.00438474
0.00129004
0.994755
0.00100777
Case  7
0.00564449
0.00274183
0.994917
0.00119178Table 8
P value of MOPSO (Cases  3,  5, and  7).
Case  3
Case  5
Case  7
Case  3
0.4727
0.9097
Case  5
0.3447
0.4727
Case  7
0.2123
0.7337
Bold indicates epsilon; the other is hybervolume.Pareto fronts produced by MOPSO-II with same number of function evaluations.
(a)
Cases 2, 4
(b)
Cases 6, 8
(c)
Cases 3, 5, and 7The comparisons between Case 2 and Case 4 show that larger size of swarm can obtain better average epsilon indicator with aP value of 67.76% but worse average hybervolume indicator with a P value of 14.05%. The comparisons between Case 6 and Case 8 show that larger size of swarm both increase average performance in epsilon and hybervolume, with a P value of 47.27% and 7.57%. By comparing Cases 3, 5, and 7, Case 5 demonstrates the best performance in total. However, the improvements are not evident as the P values are all a large value.Although the result obtained from this experiment seems to be inconclusive on which parameter is the best, we can argue that the MOPSO’s performance is not sensitive to the main PSO parameters, under the condition that the total number of function evaluation retains a large value. We have obtained competitive results in most cases without paying special attentions on PSO parameters.


## 4.3.1. Boundary Constraint Optimization Comparisons

As seen from Figures1 and 2, the boundary constraint optimization method of the MOPSO has a great influence on the Pareto fronts. Graphically, the MOPSO-II has much better performance in diversity, and the spread of its solutions found is much larger than that of the MOPSO-I. This can easily explain that probability-based disposal to boundary constraint can enrich the diversity of swarm flight.It is not easily to determine which one is much closer to the true Pareto fronts from Figure1, and the two-impulse comparison case is shown in Figure 4(a), while it is easy to see from Figure 2 that the Pareto fronts obtained by MOPSO-I are much closer to the true Pareto fronts, and the three-impulse comparison case is redrawn in Figure 4(b) to demonstrate this clearly.Pareto fronts produced by MOPSO with different boundary constraint optimization.
(a)
Asteroid 1, two-impulse
(b)
Asteroid 2, three-impulseTherefore, the quantitative metrics are calculated and the statistical results are provided in Table3. As seen from Table 3, the MOPSO-II has better performance with respect to the hybervolume indicator (larger value means better) for the two cases. While, for the epsilon indicator (smaller value means better), the MOPSO-II is much better in case of Asteroid 1, two-impulse problem and the MOPSO-I are much better in case of Asteroid 2, three-impulse.Table 3
MOPSO performance with different boundary constraint optimization.
Test problem
Algorithm
Epsilon
Hybervolume
Mean
std
P value
Mean
std
P value
Asteroid 1,two-impulse
MOPSO-I
0.011987
0.008378
0.1620
0.972618
0.006829
0.0257
MOPSO-II
0.006833
0.004066
0.977958
0.001301
Asteroid 2,three-impulse
MOPSO-I
0.025630
0.0278051
0.1620
0.87509
0.308357
0.6232
MOPSO-II
0.0399274
0.026254
0.968537
0.016548Through the comparisons provided in Table2 and Figure 4, our proposed boundary constraint demonstrates much better performance in total. Besides, these comparisons show that the same algorithm will have different performance in solving the same type of problem with different configurations.


## 4.3.2. Algorithm Parameters Analysis

Our simulation experiments show thatnparticles and GMax are the two main parameters affecting the performance of the MOPSO in solving this complex multiobjective problem. In order to quantitatively evaluate their influence o, nparticles is, respectively, set as 400, 200, and 100, GMax is also, respectively, set as 400, 200, and 100, and a total of nine groups of parameters (in Table 4) are furthered tested. Other parameters are the same as used in Section 4.2.Table 4
MOPSO parameters configurations.
Index
Case  1
Case  2
Case  3
Case  4
Case  5
Case  6
Case  7
Case  8
Case  9
n
particles
400
400
400
200
200
200
100
100
100
GMax
400
200
100
400
200
100
400
200
100Ten independent runs for the MOPSO-II with each group of parameters in solving the three-impulse, Asteroid 1 problem, are executed.Figure5(a) compares the Pareto fronts of Cases 1, 5, and 9, Figure 5(b) compares that of Cases 2, 5, and 8, and Figure 5(c) compares that of Cases 7, 8, and 9. These comparisons are to show the influence of the total number of function evaluations. It is clearly seen from Figure 5 that the larger number of function evaluations evidently increases the performance of the MOPSO. Thus, a larger number of swarm size and iteration are necessary for obtaining the optimal Pareto fronts for this practical multiobjective optimization problem.Pareto fronts produced by MOPSO-II with different number of function evaluations.
(a)
Cases 1, 5, and 9
(b)
Cases 2, 5, and 8
(c)
Cases 7, 8, and 9Figure6(a) compares the Pareto fronts of Cases 2 and 4, Figure 6(b) compares that of Cases 6 and 8, and Figure 6(c) compares that of Cases 3, 5, and 7. These comparisons are to show the influence of nparticles and GMax with the same number of function evaluations. Evidently, we cannot determine which case is better by the graphical results provided in Figure 6. Therefore, the quantitative metrics are calculated and the statistical results are provided in Tables 5, 6, 7, and 8.Table 5
MOPSO-II performance with different parameters (Cases  2 and  4).
Test problem
Algorithm
Epsilon
Hybervolume
Mean
std
P value
Mean
std
P value
Asteroid 1,three-impulse
Case  2
0.010692
0.006797
0.6776
0.984063
0.002738
0.1405
Case  4
0.012918
0.012848
0.985643
0.002736Table 6
MOPSO-II performance with different parameters (Cases  6 and  8).
Test problem
Algorithm
Epsilon
Hybervolume
Mean
std
P value
Mean
std
P value
Asteroid 1,three-impulse
Case  6
0.006045
0.002645
0.4727
0.99544
0.000816
0.0757
Case  8
0.007806
0.004042
0.994456
0.001074Table 7
MOPSO-II performance with different parameters (Cases  3,  5, and  7).
Test problem
Algorithm
Epsilon
Hybervolume
Mean
std
Mean
std
Asteroid 1,three-impulse
Case  3
0.00583413
0.0036399
0.994062
0.00135846
Case  5
0.00438474
0.00129004
0.994755
0.00100777
Case  7
0.00564449
0.00274183
0.994917
0.00119178Table 8
P value of MOPSO (Cases  3,  5, and  7).
Case  3
Case  5
Case  7
Case  3
0.4727
0.9097
Case  5
0.3447
0.4727
Case  7
0.2123
0.7337
Bold indicates epsilon; the other is hybervolume.Pareto fronts produced by MOPSO-II with same number of function evaluations.
(a)
Cases 2, 4
(b)
Cases 6, 8
(c)
Cases 3, 5, and 7The comparisons between Case 2 and Case 4 show that larger size of swarm can obtain better average epsilon indicator with aP value of 67.76% but worse average hybervolume indicator with a P value of 14.05%. The comparisons between Case 6 and Case 8 show that larger size of swarm both increase average performance in epsilon and hybervolume, with a P value of 47.27% and 7.57%. By comparing Cases 3, 5, and 7, Case 5 demonstrates the best performance in total. However, the improvements are not evident as the P values are all a large value.Although the result obtained from this experiment seems to be inconclusive on which parameter is the best, we can argue that the MOPSO’s performance is not sensitive to the main PSO parameters, under the condition that the total number of function evaluation retains a large value. We have obtained competitive results in most cases without paying special attentions on PSO parameters.


## 4.4. Comparison of the MOPSO with Other MOEA

To demonstrate the performance of the employed MOPSO, two other popular algorithms are tested and compared: the NSGA-II and the SPEA-II.In the following examples, the total number of function evaluations was set to 20000 for all the algorithms compared. The NSGA-II and the SPEA-II were run using a population size of 200, a maximum number of generations of 100. The MOPSO-I and MOPSO-II used 200 particles, a maximum number of generations of 20, a maximum number of generations per swarm of 5, and a total of 5 swarms. The source code of the NSGA-II, and SPEA-II provided in the EMOO repository is also employed in this study.The two test problems are the same to those employed in Section4.3 to compare boundary constraint optimization. Ten independent runs for each test case are executed. The final Pareto fronts are compared in Figure 7 and the statistical results with respect to epsilon and hybervolume indicators are provided in Tables 9, 10, and 11.Table 9
Performance comparisons of MOPSO, NSGA-II, and SPEA-II.
Test problem
Algorithm
Epsilon
Hybervolume
Mean
std
Mean
std
Asteroid 1,two-impulse
MOPSO-I
0.035041
0.038848
0.876849
0.026781
NSGA-II
0.113438
0.094049
0.821205
0.080834
SPEA-II
0.100881
0.032094
0.811347
0.035008
Asteroid 2,three-impulse
MOPSO-II
0.159179
0.081796
0.875023
0.078886
NSGA-II
0.0594614
0.058546
0.838198
0.301149
SPEA-II
0.12172
0.055917
0.902724
0.072393Table 10
P value of MOPSO-I, NSGA-II, and SPEA-II.
MOPSO-I
NSGA-II
SPEA-II
MOPSO-I
0.0058
0.0113
NSGA-II
0.0073
0.5205
SPEA-II
0.0017
0.1620
Bold indicates epsilon; the other is hybervolume.Table 11
P value of MOPSO-II, NSGA-II, and SPEA-II.
MOPSO-II
NSGA-II
SPEA-II
MOPSO-II
0.0113
0.3847
NSGA-II
0.3447
0.0211
SPEA-II
0.4274
0.5205
Bold indicates epsilon; the other is hybervolume.Pareto fronts produced by MOPSO, NSGA-II, and SPEA-II.
(a)
Asteroid 1, two-impulse
(b)
Asteroid 2, three-impulseEvidently seen from Figure7(a), the MOPSO-I produces better Pareto fronts than the NSGA-II and SPEA-II for the first test problem. The statistical results of the two indicators also support this. The P value with respect to epsilon is calculated as 0.58% and 1.13%, respectively, and 0.73% and 0.17% for hybervolume. This shows that the MOPSO-I is absolutely super to the NSGA-II and SPEA-II.Also evidently seen from Figure7(b), the MOPSO-I produces much closer Pareto fronts than the NSGA-II and SPEA-II for the second test problem. However, its Pareto fronts are a little narrower, which may result in that its epsilon and hybervolume indicator are junior to the NSGA-II and SPEA-II. Furthermore, the fact that the Pareto fronts are highly discontinuous would make these metrics irrelevant. In the view of practical spacecraft mission design, obtaining much closer Pareto fronts from which the designer can choose one single satisfying solution for engineering design would be much desirable performance for a multiobjective optimizer. Thus, we can say that the MOPSO-I is a little better than the NSGA-II and SPEA-II; even the statistical results on quantitative metrics do not support this.In conclusion, the MOPSO showed a better performance with respect to the NSGA-II and the SPEA-II for the two test problems. For other test problems, the MOPSO is not always better, but its performance is very competitive.


## 4.5. Comparison with PCP Method

The widely used Pork-Chop (PCP) method is also tested here for comparison. TheΔv is calculated every 1 day in both Earth departure time and asteroid arrival time, and the two-impulse Lambert algorithm is employed. The contours of time of flight corresponding to the optimal transfers for Asteroid 1 are illustrated in Figure 8. For convenience, only the solutions with a Δv less than 18 km/s and Earth departure and asteroid arrive times in [5000,7000] (MJD2000) are presented.Relations between total characteristic velocity, transfer time, and earth departure time.
(a)
Asteroid 1, two-impulse
(b)
Asteroid 2, two-impulseWe analyze the comparisons between the proposed multiobjective optimization approach and the PCP method in the following three aspects.


### 4.5.1. Uses of the Results

In the PCP method, the contours of theΔv are plotted to assist the designer to find the best launch window and transfer trajectory. However, the contours have many local peaks and it is not very convenient to determine which solution is the best one. Besides, the size of the region containing the global minimum is a very important factor in finding real trajectories in any given year. This is not easily observed from the Pork-Chop figures.In contrast to the PCP method, the proposed multiobjective optimization approach can provide friendly the designer of this information. The relations between total characteristic velocity, transfer time, and earth departure time of those obtained Pareto-optimal solutions for two test cases are provided in Figure8. From Figure 8, it can be clearly observed that the earth departure time for all the optimal rendezvous trajectories concentrates on one arrow domain for these two cases. The optimized Earth launch window and its size can be easily determined through the multiobjective optimization design, which will provide very useful information for engineering design.


### 4.5.2. Global Optimality

The minimum-propellant solution searched by the PCP method for the first asteroid mission is calculated with aΔv of 6015.3 m/s and it is located in Figure 9 (* denotes). As seen in Figure 1, the MOPSO has located a set of Pareto solutions with a Δv about this value. This issue can also demonstrate the global convergence ability of the MOPSO. Besides, the PCP method obtains only the two-impulse trajectory. As is well known, the two-impulse trajectory is not the propellant-optimal solution under most conditions, and increasing the number of impulses will reduce the propellant cost. Our proposed approach can design the multi-impulse trajectory; thus it locates better solution than the PCP method.Figure 9
Contours of time of flight corresponding to the optimal transfers (Asteroid 1).


### 4.5.3. Efficiency

The PCP is in essence of an exhaustive searching method. In our test, the search space for departure time and arrive time is[4000,10000] day. With the search step of 1 day, the PCP method calculates a total number of 6000*6000 trajectories, while only 400*400*5 trajectories are calculated in the proposed method for the most highly cost case. The proposed approach improves the calculation efficiency by about 40 times. Besides, the PCP method needs much human intention to determine the final solutions, while the proposed approach is an automated search method.


## 4.5.1. Uses of the Results

In the PCP method, the contours of theΔv are plotted to assist the designer to find the best launch window and transfer trajectory. However, the contours have many local peaks and it is not very convenient to determine which solution is the best one. Besides, the size of the region containing the global minimum is a very important factor in finding real trajectories in any given year. This is not easily observed from the Pork-Chop figures.In contrast to the PCP method, the proposed multiobjective optimization approach can provide friendly the designer of this information. The relations between total characteristic velocity, transfer time, and earth departure time of those obtained Pareto-optimal solutions for two test cases are provided in Figure8. From Figure 8, it can be clearly observed that the earth departure time for all the optimal rendezvous trajectories concentrates on one arrow domain for these two cases. The optimized Earth launch window and its size can be easily determined through the multiobjective optimization design, which will provide very useful information for engineering design.


## 4.5.2. Global Optimality

The minimum-propellant solution searched by the PCP method for the first asteroid mission is calculated with aΔv of 6015.3 m/s and it is located in Figure 9 (* denotes). As seen in Figure 1, the MOPSO has located a set of Pareto solutions with a Δv about this value. This issue can also demonstrate the global convergence ability of the MOPSO. Besides, the PCP method obtains only the two-impulse trajectory. As is well known, the two-impulse trajectory is not the propellant-optimal solution under most conditions, and increasing the number of impulses will reduce the propellant cost. Our proposed approach can design the multi-impulse trajectory; thus it locates better solution than the PCP method.Figure 9
Contours of time of flight corresponding to the optimal transfers (Asteroid 1).


## 4.5.3. Efficiency

The PCP is in essence of an exhaustive searching method. In our test, the search space for departure time and arrive time is[4000,10000] day. With the search step of 1 day, the PCP method calculates a total number of 6000*6000 trajectories, while only 400*400*5 trajectories are calculated in the proposed method for the most highly cost case. The proposed approach improves the calculation efficiency by about 40 times. Besides, the PCP method needs much human intention to determine the final solutions, while the proposed approach is an automated search method.


## 5. Conclusions

The paper formulates the asteroid rendezvous preliminary trajectory design as a multiobjective optimization problem and employs the multiobjective particle swarm optimization (MOPSO) algorithm to locate the Pareto-optimal solution set. Compared with the widely employed Pork-Chop method, the proposed approach is demonstrated to be able to provide much more easily used results, obtain better propellant-optimal solutions, and have much better efficiency. The results show that the proposed approach can effectively and efficiently demonstrate the relations among the mission characteristic parameters such as launch time, transfer time, propellant cost, and number of maneuvers, which will provide useful reference for practical asteroid mission design. The MOPSO proves to be quite effective in finding the Pareto-optimal solutions and its performance can be improved by a proposed boundary constraint optimization strategy. The MOPSO is found to be very competitive with respect to two highly competitive multiobjective evolutionary algorithms: the NSGA-II and the SPEA-II.



---

*Source: 823659-2014-03-31.xml*