# Optimization of Interplanetary Trajectories Using the Colliding Bodies Optimization Algorithm

Authors: Marco Del Monte; Raffaele Meles; Christian Circi
Journal: International Journal of Aerospace Engineering
 (2020)

---

## Abstract

In this paper, a recent physics-based metaheuristic algorithm, the Colliding Bodies Optimization (CBO), already employed to solve problems in civil and mechanical engineering, is proposed for the optimization of interplanetary trajectories by using both indirect and direct approaches. The CBO has an extremely simple formulation and does not depend on any initial conditions. To test the performances of the algorithm, missions with remarkably different orbital transfer energies are considered: from the simple planar case, as the Earth-Mars orbital transfer, to more energetic ones, like a rendezvous with the asteroid Pallas.


---

## Body

## 1. Introduction

An important aspect of a space mission design is the obtaining of the nominal optimal trajectory, traditionally the one with minimum transfer time or maximum payload mass. Trajectory optimization is an old problem, its origins date back to the ancient Greeks, but its rigorous mathematical formulation, as an optimal control problem, arrives only with Pontryagin in the mid-1900s [1]. Optimal control is an issue concerning the determination of the inputs into a dynamical system that optimize (i.e., minimize or maximize) a specified performance index while satisfying several constraints [2]. These constraints can be differential, as the equations of motion, or algebraic, as departure, mid-course, and arrival constraints. Because of the complexity of most applications, optimal control problems are chiefly solved numerically. Numerical methods adopted are divided into two major classes: indirect methods and direct methods. In the former, the original problem is transcribed into a multiple-point boundary-value problem that is solved to determine candidate optimal trajectories, and the optimization phase consists in finding the optimal set of costate variables. In a direct method, the state and/or control of the optimal control problem is discretized, and the problem is transcribed into a nonlinear optimization problem [2]. Both approaches lead to a parametric optimization where a set of optimal parameters must be found. This optimization has often been conducted by means of gradient-based research (e.g., Newton-Raphson-based algorithms), but in the last decades, a new kind of optimization procedures has been proposed and developed: metaheuristics [3, 4]. The goal of metaheuristics is to efficiently explore the search space looking for near-optimal solutions. The fundamental characteristics are that they are problem independent and they possess a nondeterministic nature, useful to escape from local optima. This is achieved by either allowing worsening moves or generating new starting solutions for the local search in a more “intelligent” way than just providing random initial solutions. This stochastic nature is not employed blindly, but in an intelligent, biased manner, and is what truly differentiates them from gradient-based techniques, which are deterministic and strongly dependent on an initial guess of the solution [5, 6]. Gradient-based algorithms are largely employed in all fields of engineering, including space trajectories optimization. However, in recent years, metaheuristic algorithms have been increasingly adopted, especially in preliminary analysis of trajectories. Two of the most used algorithms, in this domain, are Particle Swarm Optimization [7, 8] and Differential Evolution [9, 10].This paper is focused on the optimization of space trajectories using the Colliding Bodies Optimization algorithm. This is a novel population-based metaheuristic inspired by the one-dimensional collision theory between bodies, where each candidate solution being considered as a body with mass. CBO utilizes a simple formulation to find extremals of functions and does not depend on any internal parameter [11]. In Section 2, CBO and its enhanced version will be described briefly. Test cases using indirect methods are studied in Section 3, while those using direct methods are studied in Section 4. Conclusions will be given in Section 5.


## 2. Colliding Bodies Optimization Algorithm

Colliding Bodies Optimization is a metaheuristic algorithm developed by Kaveh and Mahdavi [11–13], inspired by the one-dimensional collision theory. There are two versions of this algorithm, a basic one [11] and an enhanced one [12], that improves the basic version by means of a sort of Elitism and Crossover. In the next two sections, some basic statements are reported while details can be found in the cited references.


### 2.1. Basic CBO Formulation

Each search agent is modelled as a body with mass and velocity. The initial position of theith body is randomly provided in a j-dimensional search space set by the user:
(1)xij=xj,min+rand⋅xj,max−xj,min,where rand is a random number between 0 and 1. A collision occurs between two bodies, and their positions, after the impact, are updated based on the one-dimensional collision laws [11, 13]. Given the body Xk (also called particle or object), its mass is defined as follows:
(2)mk=1/Jk1/∑i=1n1/Ji,k=1,⋯,n,where Jk is the cost function value of the kth particle and n, which must be an even number, is the total number of bodies used in the optimization process (the population size). The n colliding bodies (CBs) are sorted into ascending order, according to their objective function values, and then divided into two equal groups: Stationary Objects (the lower half) and Moving Objects (the upper half). Objects of the MO group collide against members of the SO group to improve their position and push stationary objects towards better positions. In particular, the colliding pairs are established according to the ascending order with respect to the objective function. Hence, for instance, the best moving particle collides with the best stationary one. Bodies’ velocities before the collision are assigned as follows:
(3)Stationarybodies:vi=0,i=1,⋯,n2,(4)Movingbodies:vi=xi−n/2−xi,i=n2+1,⋯,n.As many other metaheuristic algorithms, velocities are not defined as the derivative of the position with respect to time, but they are expressed as displacements in the search space. According to the colliding bodies’ theory, velocities after the collision are calculated as follows:(5)Stationarybodies:vi′=mi+n/2+εmi+n/2vi+n/2mi+mi+n/2,i=1,⋯,n2,(6)Movingbodies:vi′=mi−εmi−n/2vimi+mi−n/2,i=n2+1,⋯,n,where ε is the Coefficient of Restitution, defined as the ratio of the relative velocity between two bodies after and before the collision:
(7)ε=vi+1′−vi′vi+1−vi.This coefficient is assumed varying linearly between 1 and 0 during the optimization process, in order to ensure the balance between exploration and exploitation. After the calculation of the displacement, it is possible to determine new positions of the stationary and moving bodies as follows:(8)Stationarybodies:xinew=xi+rand⋅vi′,i=1,⋯,n2,(9)Movingbodies:xinew=xi−n/2+rand⋅vi′,i=n2+1,⋯,n,where rand is a uniformly distributed random vector in the range [-1,1]. This iterative scheme, performed on all the particles at each iteration, is repeated until a given stopping criterion is fulfilled. A Pseudocode 1 of CBO is reported.Pseudocode 1
1. Initialize the CBO population in the search space (Equation (1))
2. Evaluate the objective functions and define the masses as in Equation (2)
3. Sort the population in order to identify stationary and moving groups and calculate the velocities as in Equations (3) and (4)
4. Calculate the velocity after the collisions by means of Equations (5) and (6)
5. The new positions can be determined by Equations (8) and (9)
6. If the terminating criterion is fulfilled, proceed to step 7; otherwise, go to step 2
7. Report the best solution found by the algorithm
8. END


### 2.2. Enhanced CBO Formulation

The structure of the enhanced CBO (ECBO) algorithm is essentially the same as the basic CBO [12], with the difference that a Colliding Memory (CM) is introduced to save the best CBs’ positions obtained so far. In fact, the positions stored in the Colliding Memory substitute the worst positions occupied by the current bodies. In this way, the best positions are remembered and there is no global worsening of the objective function from one iteration to another. The number of the best CBs’ positions that are preserved, therefore the dimension of the Colliding Memory, is set by the user. Moreover, after the update of the CBs’ positions, the ECBO executes the following crossover instruction:
(10)ifrani<PRO,xij=xj,min+randxj,max−xj,min,otherwise,xij=xij,where xij is the jth variable of the ith CB, randomly selected; xj,min,xj,max are the lower and the upper bounds of the jth variable; PRO is the crossover probability that must be set by the user between [0,1]; rani is a random number uniformly distributed within [0,1], automatically generated for each particle, as well as rand. If rani is less than PRO, a crossover occurs. As PRO increases, the probability to perform a crossover increases. The Pseudocode 2 of the ECBO is reported.Pseudocode 2
1. Initialize the CBO population in the search space (Equation (1))
2. Evaluate the objective functions and define the masses as in Equation (2)
3. Update the Colliding Memory (CM) and population
4. Sort the population in order to identify stationary and moving groups and calculate the velocities as in Equations (3) and (4)
5. Calculate the velocity after the collisions by means of Equations (5) and (6)
6. The new positions can be determined by Equations (8), (9), and (10)
7. If the terminating criterion is fulfilled, proceed to step 8; otherwise, go to step 2
8. Report the best solution found by the algorithm
9. END


## 2.1. Basic CBO Formulation

Each search agent is modelled as a body with mass and velocity. The initial position of theith body is randomly provided in a j-dimensional search space set by the user:
(1)xij=xj,min+rand⋅xj,max−xj,min,where rand is a random number between 0 and 1. A collision occurs between two bodies, and their positions, after the impact, are updated based on the one-dimensional collision laws [11, 13]. Given the body Xk (also called particle or object), its mass is defined as follows:
(2)mk=1/Jk1/∑i=1n1/Ji,k=1,⋯,n,where Jk is the cost function value of the kth particle and n, which must be an even number, is the total number of bodies used in the optimization process (the population size). The n colliding bodies (CBs) are sorted into ascending order, according to their objective function values, and then divided into two equal groups: Stationary Objects (the lower half) and Moving Objects (the upper half). Objects of the MO group collide against members of the SO group to improve their position and push stationary objects towards better positions. In particular, the colliding pairs are established according to the ascending order with respect to the objective function. Hence, for instance, the best moving particle collides with the best stationary one. Bodies’ velocities before the collision are assigned as follows:
(3)Stationarybodies:vi=0,i=1,⋯,n2,(4)Movingbodies:vi=xi−n/2−xi,i=n2+1,⋯,n.As many other metaheuristic algorithms, velocities are not defined as the derivative of the position with respect to time, but they are expressed as displacements in the search space. According to the colliding bodies’ theory, velocities after the collision are calculated as follows:(5)Stationarybodies:vi′=mi+n/2+εmi+n/2vi+n/2mi+mi+n/2,i=1,⋯,n2,(6)Movingbodies:vi′=mi−εmi−n/2vimi+mi−n/2,i=n2+1,⋯,n,where ε is the Coefficient of Restitution, defined as the ratio of the relative velocity between two bodies after and before the collision:
(7)ε=vi+1′−vi′vi+1−vi.This coefficient is assumed varying linearly between 1 and 0 during the optimization process, in order to ensure the balance between exploration and exploitation. After the calculation of the displacement, it is possible to determine new positions of the stationary and moving bodies as follows:(8)Stationarybodies:xinew=xi+rand⋅vi′,i=1,⋯,n2,(9)Movingbodies:xinew=xi−n/2+rand⋅vi′,i=n2+1,⋯,n,where rand is a uniformly distributed random vector in the range [-1,1]. This iterative scheme, performed on all the particles at each iteration, is repeated until a given stopping criterion is fulfilled. A Pseudocode 1 of CBO is reported.Pseudocode 1
1. Initialize the CBO population in the search space (Equation (1))
2. Evaluate the objective functions and define the masses as in Equation (2)
3. Sort the population in order to identify stationary and moving groups and calculate the velocities as in Equations (3) and (4)
4. Calculate the velocity after the collisions by means of Equations (5) and (6)
5. The new positions can be determined by Equations (8) and (9)
6. If the terminating criterion is fulfilled, proceed to step 7; otherwise, go to step 2
7. Report the best solution found by the algorithm
8. END


## 2.2. Enhanced CBO Formulation

The structure of the enhanced CBO (ECBO) algorithm is essentially the same as the basic CBO [12], with the difference that a Colliding Memory (CM) is introduced to save the best CBs’ positions obtained so far. In fact, the positions stored in the Colliding Memory substitute the worst positions occupied by the current bodies. In this way, the best positions are remembered and there is no global worsening of the objective function from one iteration to another. The number of the best CBs’ positions that are preserved, therefore the dimension of the Colliding Memory, is set by the user. Moreover, after the update of the CBs’ positions, the ECBO executes the following crossover instruction:
(10)ifrani<PRO,xij=xj,min+randxj,max−xj,min,otherwise,xij=xij,where xij is the jth variable of the ith CB, randomly selected; xj,min,xj,max are the lower and the upper bounds of the jth variable; PRO is the crossover probability that must be set by the user between [0,1]; rani is a random number uniformly distributed within [0,1], automatically generated for each particle, as well as rand. If rani is less than PRO, a crossover occurs. As PRO increases, the probability to perform a crossover increases. The Pseudocode 2 of the ECBO is reported.Pseudocode 2
1. Initialize the CBO population in the search space (Equation (1))
2. Evaluate the objective functions and define the masses as in Equation (2)
3. Update the Colliding Memory (CM) and population
4. Sort the population in order to identify stationary and moving groups and calculate the velocities as in Equations (3) and (4)
5. Calculate the velocity after the collisions by means of Equations (5) and (6)
6. The new positions can be determined by Equations (8), (9), and (10)
7. If the terminating criterion is fulfilled, proceed to step 8; otherwise, go to step 2
8. Report the best solution found by the algorithm
9. END


## 3. Numerical Simulations: Indirect Methods

In order to analyze the performances of this optimization algorithm, a total of five study cases will be presented. Two of them are solved by using an indirect strategy, while the remaining three cases adopt a direct approach.


### 3.1. Optimal Earth to Mars Orbital Transfer

The problem is to reach the orbit of Mars departing from the Earth’s orbit with the minimum transfer time by using a low thrust engine. This case has already been studied in literature with different techniques [7, 14, 15]. In this paper, the best trajectory is obtained by using the ECBO. As in the cited papers, the following hypotheses are established:
(1)
The orbits of the planets are coplanar and circular
(2)
The only attracting body is the Sun
(3)
The spacecraft’s initial position and velocity are the same as the Earth’s
(4)
The thrust magnitude of the spacecraft is constantThe spacecraft’s equations of motion are written in polar coordinates:(11)ṙ=vr,v̇r=−μs−rvθ2r2+Tmsinα,v̇θ=−vrvθr+Tmcosα,where r is the position vector; vr and vθ are, respectively, the radial and the horizontal velocity of the spacecraft, μs is the Sun’s gravitational parameter, T is the thrust magnitude; m is the spacecraft mass, and α is the thrust pointing angle relative to the local horizontal (Figure 1).Figure 1
Thrust pointing angleα.The control vector isut=α. The thrust-to-mass ratio is as follows [7]:
(12)Tm=Tm0−T/ct−t0=cn0c−n0t−t0,where m0 is the initial spacecraft mass, c=55.894km/s is the exhaust velocity, and n0=8.342×10−4m/s2 is the thrust-to-mass ratio at the initial time t0. A normalized set of units is as follows:
(i)
Distance unitDU=149.5⋅106km is the mean radius of the Earth’s orbit
(ii)
Time unitTU=5.018⋅106s, such that μs=DU3/TU2The objective function to minimize isJ=tf. The desired terminal conditions are expressed by the vector ω as follows:
(13)ω=ΔrΔvrΔvθ=rtf−RMarsvrtfvθtf−μsRMars=000.To write the necessary conditions for optimality, the Hamiltonian function has to be defined:(14)H=prvr+pvr−μs−rvθ2r2+Tmsinu+pvθ−vrvθr+Tmcosu,where the time-dependent set of costate variables pr,pvr,pvθ has been introduced. Furthermore, the costate differential equations are expressed as follows:
(15)ṗr=vθ2pvr−vrvθpvθr2−2μspvrr3,ṗvr=−pr+vθpvθr,ṗvθ=−2vθpvr+vrpvθr.From the Pontryagin principle, the optimal control law is as follows:(16)cosuopt=−pvθpvθ2+pvr2,sinuopt=−pvrpvθ2+pvr2.The set of necessary conditions for optimality can be completed with the transversality condition referred to the final time:(17)Htf=−∂J∂tf−λ⋅∂ω∂tf,where λ is another set of adjoint variables concerning the terminal conditions; hence, it has the same dimensions as ω. Rearranging Equation (17), the following condition is obtained:
(18)cn0c−n0tfpvrtf2+pvθtf2−1=0.In order to find the optimal trajectory that satisfies the necessary conditions, the problem reduces to the determination of four parameters:pr0,pvr0,pvθ0,tf. The CBO algorithm must find the solution exploring the following search space:
(19)1TU≤tf≤10TU,−1≤pr0,pvr0,pvθ0≤1.It is necessary to introduce a cost function that links the optimization algorithm to the problem. This function always includes the quantity that must be minimized, and if there is any constraint to respect, the most popular approach is to include them in the cost function. Then, the objective function is defined here as follows:(20)J=tf+100⋅Δr+100⋅Δvr+100⋅Δvθ.In the cost function, all the quantities (transfer duration and errors at final time) are dimensionless and the weights are chosen to scale them properly, in order to keep the balance between all terms throughout the simulation. The ECBO population is composed of 50 particles. The dimension of the Colliding Memory is set to 1 and crossover is not performed. The algorithm has to find the optimal set of costate initial values. The transversality condition can be neglected by the CBO, due to the homogeneity of the costate equations (Equation (15)). For this reason, if the CBO finds a set of initial costate variables that is proportional to the optimal one (p→0=bp→0opt), the same proportionality holds at any time. By writing the control as a function of the proportional set, it is possible to demonstrate that it coincides with the optimal control (Equation (16)). This means that the minimum time trajectory can be obtained also with an initial costate proportional to the optimal set. Nevertheless, the transversality condition will be violated, and by substituting the proportional set in Equation (18), it can be easily proved that Htfopt=−b. In order not to increase the number of equality constraints in the cost function, the transversality condition is verified at the end of the optimization process, to guarantee the optimality of the trajectory.The optimal solution found by the CBO consists in a 192.6 days transfer (Figure2).Figure 2
Earth-Mars optimal transfer trajectory.Figure3 shows the optimal thrust pointing angle αt and the costate evolution during the transfer trajectory. These results are in accordance with those in literature by using PSO [7] and are obtained quickly and easily thanks to the simplicity of the algorithm.Figure 3
Optimal costate and control evolution in the Earth-Mars transfer.


### 3.2. Optimal Earth to Mercury Orbital Transfer

A minimum time transfer between the orbits of the Earth and Mercury is studied by means of a solar sail. The problem consists in determining the optimal steering lawαt that minimizes the time of flight to reach Mercury’s orbit [16].As in the cited papers, a series of simplifying assumptions are made:(1)
The relative orbital inclination of Mercury and Earth is neglected, and the orbits are considered circular
(2)
The spacecraft’s initial position and velocity are the same as the Earth’s
(3)
The only attracting body is the SunA polar inertial reference frame is used, and the equations of motion for the solar sail are as follows:(21)ṙ=vr,θ̇=vθr,v̇r=vθ2r−μsr2+acb1+b2+b3REarthr2⋅cosαb1+b2cos2α+b3cosα,v̇θ=−vrvθr+acb1+b2+b3REarthr2⋅sinαcosαb2cosα+b3,where r is the position vector and θ is the polar angle measured anticlockwise from the axis that connects the Sun to the Earth at the initial instant. vr and vθ are, respectively, the radial and the tangential velocity. α is the angle between the direction that connects the Sun to the solar sail and the thrust direction (Figure 4); ac is the characteristic acceleration of the sail, and the terms b1, b2, and b3 are the coefficients that represent the optical properties of the sail. A sail with an aluminium-coated front side and a chromium-coated back side with b1=0.1728, b2=1.6544, and b3=−0.0109 is considered [16].Figure 4
Thrust pointing angleα.The minimum time trajectory is obtained with an indirect approach. The Hamiltonian function is as follows:(22)H=prvr+vθpθ−pvθvrr+pvrvθ2r−μsr2+accosαREarth/r2b1+b2+b3pvrb1+b2cos2α+b3cosα+pvθsinαb2cosα+b3,where pr,pθ,pvr,pvθ are the time-dependent costate variables. The Euler-Lagrange equations are as follows:
(23)ṗr=pθvθr2+pvrvθ2r2−2μsr3+2acREarth2r3b1+b2+b3cosαb1+b2cos2α+b3cosα−pvθvrvθr2−2acREarth2r3b1+b2+b3sinαcosαb2cosα+b3,ṗθ=0,ṗvr=−pr+vθpvθr,ṗvθ=−2vθpvr+vrpvθ−pθr.The optimal control will be obtained in the formα=αpvr,pvθ, but it is not possible to find an explicit solution in this form. The optimal steering law can be approximated by [17]:
(24)α=signpvθα~ifαp<αp∗,signpvθπ2ifαp≥αp∗,where
(25)cosαp=pvrpvθ2+pvr2withαp∈0,π,αp∗≅2.5392rad,α~≅0.008109αp6−0.05474αp5+0.1356αp4−0.1266αp3+0.08266αp2+0.3038αp+0.0008666withαp∈0,αp∗.The equations of motion have the following boundary conditions:(26)rt0=REarth,θt0=vrt0=0,vθt0=μsREarth.There are four unknown parameters:pr0,pvr0,pvθ0,tf, and the optimal values are found in the following search space:
(27)1TU≤tf≤20TU,−1≤pr0pvr0pvθ0≤1.To investigate the domain, a population of 50 particles is considered. As in the previous case, an ECBO that preserves the best position in the population at each iteration is employed. Therefore, the selected dimension of the Colliding Memory is 1 and the crossover is not considered. The cost function is written, as in the previous case:(28)J=0.01⋅tf+100⋅Δr+40⋅Δvr+40⋅Δvθ.Although the cost function has the same form as the previous one, coefficients must be chosen carefully and they are different for each problem. The characteristic acceleration is set to0.25mm/s2 and the resulting transfer lasts 2.8 years, in agreement with the results obtained in [16]. The optimal costate and α time history are shown in Figure 5. The CBO, also in this case, demonstrates a great ability in finding the optimal initial costate values within 6500 function evaluations with no parameter tuning at all.Figure 5
Optimal costate and control angle (alpha) for the Earth-Mercury transfer.


## 3.1. Optimal Earth to Mars Orbital Transfer

The problem is to reach the orbit of Mars departing from the Earth’s orbit with the minimum transfer time by using a low thrust engine. This case has already been studied in literature with different techniques [7, 14, 15]. In this paper, the best trajectory is obtained by using the ECBO. As in the cited papers, the following hypotheses are established:
(1)
The orbits of the planets are coplanar and circular
(2)
The only attracting body is the Sun
(3)
The spacecraft’s initial position and velocity are the same as the Earth’s
(4)
The thrust magnitude of the spacecraft is constantThe spacecraft’s equations of motion are written in polar coordinates:(11)ṙ=vr,v̇r=−μs−rvθ2r2+Tmsinα,v̇θ=−vrvθr+Tmcosα,where r is the position vector; vr and vθ are, respectively, the radial and the horizontal velocity of the spacecraft, μs is the Sun’s gravitational parameter, T is the thrust magnitude; m is the spacecraft mass, and α is the thrust pointing angle relative to the local horizontal (Figure 1).Figure 1
Thrust pointing angleα.The control vector isut=α. The thrust-to-mass ratio is as follows [7]:
(12)Tm=Tm0−T/ct−t0=cn0c−n0t−t0,where m0 is the initial spacecraft mass, c=55.894km/s is the exhaust velocity, and n0=8.342×10−4m/s2 is the thrust-to-mass ratio at the initial time t0. A normalized set of units is as follows:
(i)
Distance unitDU=149.5⋅106km is the mean radius of the Earth’s orbit
(ii)
Time unitTU=5.018⋅106s, such that μs=DU3/TU2The objective function to minimize isJ=tf. The desired terminal conditions are expressed by the vector ω as follows:
(13)ω=ΔrΔvrΔvθ=rtf−RMarsvrtfvθtf−μsRMars=000.To write the necessary conditions for optimality, the Hamiltonian function has to be defined:(14)H=prvr+pvr−μs−rvθ2r2+Tmsinu+pvθ−vrvθr+Tmcosu,where the time-dependent set of costate variables pr,pvr,pvθ has been introduced. Furthermore, the costate differential equations are expressed as follows:
(15)ṗr=vθ2pvr−vrvθpvθr2−2μspvrr3,ṗvr=−pr+vθpvθr,ṗvθ=−2vθpvr+vrpvθr.From the Pontryagin principle, the optimal control law is as follows:(16)cosuopt=−pvθpvθ2+pvr2,sinuopt=−pvrpvθ2+pvr2.The set of necessary conditions for optimality can be completed with the transversality condition referred to the final time:(17)Htf=−∂J∂tf−λ⋅∂ω∂tf,where λ is another set of adjoint variables concerning the terminal conditions; hence, it has the same dimensions as ω. Rearranging Equation (17), the following condition is obtained:
(18)cn0c−n0tfpvrtf2+pvθtf2−1=0.In order to find the optimal trajectory that satisfies the necessary conditions, the problem reduces to the determination of four parameters:pr0,pvr0,pvθ0,tf. The CBO algorithm must find the solution exploring the following search space:
(19)1TU≤tf≤10TU,−1≤pr0,pvr0,pvθ0≤1.It is necessary to introduce a cost function that links the optimization algorithm to the problem. This function always includes the quantity that must be minimized, and if there is any constraint to respect, the most popular approach is to include them in the cost function. Then, the objective function is defined here as follows:(20)J=tf+100⋅Δr+100⋅Δvr+100⋅Δvθ.In the cost function, all the quantities (transfer duration and errors at final time) are dimensionless and the weights are chosen to scale them properly, in order to keep the balance between all terms throughout the simulation. The ECBO population is composed of 50 particles. The dimension of the Colliding Memory is set to 1 and crossover is not performed. The algorithm has to find the optimal set of costate initial values. The transversality condition can be neglected by the CBO, due to the homogeneity of the costate equations (Equation (15)). For this reason, if the CBO finds a set of initial costate variables that is proportional to the optimal one (p→0=bp→0opt), the same proportionality holds at any time. By writing the control as a function of the proportional set, it is possible to demonstrate that it coincides with the optimal control (Equation (16)). This means that the minimum time trajectory can be obtained also with an initial costate proportional to the optimal set. Nevertheless, the transversality condition will be violated, and by substituting the proportional set in Equation (18), it can be easily proved that Htfopt=−b. In order not to increase the number of equality constraints in the cost function, the transversality condition is verified at the end of the optimization process, to guarantee the optimality of the trajectory.The optimal solution found by the CBO consists in a 192.6 days transfer (Figure2).Figure 2
Earth-Mars optimal transfer trajectory.Figure3 shows the optimal thrust pointing angle αt and the costate evolution during the transfer trajectory. These results are in accordance with those in literature by using PSO [7] and are obtained quickly and easily thanks to the simplicity of the algorithm.Figure 3
Optimal costate and control evolution in the Earth-Mars transfer.


## 3.2. Optimal Earth to Mercury Orbital Transfer

A minimum time transfer between the orbits of the Earth and Mercury is studied by means of a solar sail. The problem consists in determining the optimal steering lawαt that minimizes the time of flight to reach Mercury’s orbit [16].As in the cited papers, a series of simplifying assumptions are made:(1)
The relative orbital inclination of Mercury and Earth is neglected, and the orbits are considered circular
(2)
The spacecraft’s initial position and velocity are the same as the Earth’s
(3)
The only attracting body is the SunA polar inertial reference frame is used, and the equations of motion for the solar sail are as follows:(21)ṙ=vr,θ̇=vθr,v̇r=vθ2r−μsr2+acb1+b2+b3REarthr2⋅cosαb1+b2cos2α+b3cosα,v̇θ=−vrvθr+acb1+b2+b3REarthr2⋅sinαcosαb2cosα+b3,where r is the position vector and θ is the polar angle measured anticlockwise from the axis that connects the Sun to the Earth at the initial instant. vr and vθ are, respectively, the radial and the tangential velocity. α is the angle between the direction that connects the Sun to the solar sail and the thrust direction (Figure 4); ac is the characteristic acceleration of the sail, and the terms b1, b2, and b3 are the coefficients that represent the optical properties of the sail. A sail with an aluminium-coated front side and a chromium-coated back side with b1=0.1728, b2=1.6544, and b3=−0.0109 is considered [16].Figure 4
Thrust pointing angleα.The minimum time trajectory is obtained with an indirect approach. The Hamiltonian function is as follows:(22)H=prvr+vθpθ−pvθvrr+pvrvθ2r−μsr2+accosαREarth/r2b1+b2+b3pvrb1+b2cos2α+b3cosα+pvθsinαb2cosα+b3,where pr,pθ,pvr,pvθ are the time-dependent costate variables. The Euler-Lagrange equations are as follows:
(23)ṗr=pθvθr2+pvrvθ2r2−2μsr3+2acREarth2r3b1+b2+b3cosαb1+b2cos2α+b3cosα−pvθvrvθr2−2acREarth2r3b1+b2+b3sinαcosαb2cosα+b3,ṗθ=0,ṗvr=−pr+vθpvθr,ṗvθ=−2vθpvr+vrpvθ−pθr.The optimal control will be obtained in the formα=αpvr,pvθ, but it is not possible to find an explicit solution in this form. The optimal steering law can be approximated by [17]:
(24)α=signpvθα~ifαp<αp∗,signpvθπ2ifαp≥αp∗,where
(25)cosαp=pvrpvθ2+pvr2withαp∈0,π,αp∗≅2.5392rad,α~≅0.008109αp6−0.05474αp5+0.1356αp4−0.1266αp3+0.08266αp2+0.3038αp+0.0008666withαp∈0,αp∗.The equations of motion have the following boundary conditions:(26)rt0=REarth,θt0=vrt0=0,vθt0=μsREarth.There are four unknown parameters:pr0,pvr0,pvθ0,tf, and the optimal values are found in the following search space:
(27)1TU≤tf≤20TU,−1≤pr0pvr0pvθ0≤1.To investigate the domain, a population of 50 particles is considered. As in the previous case, an ECBO that preserves the best position in the population at each iteration is employed. Therefore, the selected dimension of the Colliding Memory is 1 and the crossover is not considered. The cost function is written, as in the previous case:(28)J=0.01⋅tf+100⋅Δr+40⋅Δvr+40⋅Δvθ.Although the cost function has the same form as the previous one, coefficients must be chosen carefully and they are different for each problem. The characteristic acceleration is set to0.25mm/s2 and the resulting transfer lasts 2.8 years, in agreement with the results obtained in [16]. The optimal costate and α time history are shown in Figure 5. The CBO, also in this case, demonstrates a great ability in finding the optimal initial costate values within 6500 function evaluations with no parameter tuning at all.Figure 5
Optimal costate and control angle (alpha) for the Earth-Mercury transfer.


## 4. Numerical Simulation: Direct Methods

In this case, the test cases proposed are minimum time rendezvous with three different celestial bodies: an outer planet, an inner planet, and a very inclined asteroid. The transfer is considered heliocentric and the mathematical model consists in a three-dimensional restricted two body problem, with the Sun as the attractor and the spacecraft with negligible mass. The equations of motion are as follows:(29)ẋ=vxẏ=vyż=vzv̇x=−μsr3rx+Txmv̇y=−μsr3ry+Tymv̇z=−μsr3rz+Tzmwhere r=x2+y2+z2 is the position vector of the spacecraft and T is the thrust vector. The mass of the spacecraft is indicated as m and its time varying law is m=m0−ṁ⋅t, where m0 is the initial spacecraft mass and ṁ is the mass flow ratio. This set of equations is conveniently normalized using this set of units:
(i)
Distance unitDU=149.5⋅106km is the mean radius of the Earth’s orbit
(ii)
Time unitTU=5.018⋅106s, such that μs=DU3/TU2
(iii)
Mass unitMU=m0A direct single shooting technique, with constant thrust magnitude, is proposed due to its simplicity and effectiveness [18]. The trajectory is divided into a finite number of arcs NA with variable duration ∆tk, with k=1,⋯,NA. During each arc, the angles αk and ψk represent the thrust vector direction in the local reference frame r^,θ^,h^ with respect to the inertial reference frame J2000 (Figure 6).Figure 6
Thrust vector control direction in the local reference frame and its orientation with respect to J2000.In a rendezvous problem, the departure time is another variable to compute during the optimization process. The number of arcsNA is selected after a preliminary study and chosen to reduce as much as possible the computational effort [19]. Globally, there will be three variables for each arc plus the departure time. The ephemerides of the planets are DE430 (provided by the NASA SPICE tool [20]). With the transfer time tf, we need to minimize the following arrival errors as well:
(30)Δx=xstf−xttf,Δy=ystf−yttf,Δz=zstf−zttf,Δvx=vxstf−vxttf,Δvy=vystf−vyttf,Δvz=vzstf−vzttf,where subscript s indicates the spacecraft state and t represents the target body. The desired final condition, hence, is that the arrival errors (Equation (30)) are zero.The objective function is defined as follows:(31)J=1000⋅Cttf+Cd⋅Δx+Δy+Δz+CvΔvx+Δvy+Δvz.The coefficientsCt,Cd and Cv (that weight of the different dimensionless contributes in the cost function) are chosen after a preliminary study to properly balance the different terms, and their values will be presented for each test case. In the next section, three rendezvous missions will be discussed and the selected celestial bodies are Mars, Venus, and Pallas. The stopping criterion is composed of two conditions. The first one concerns the accuracy of the mission: each transfer will be considered successful, and hence, the relative simulation will stop, if the errors on position and velocity, described in Equation (30), go, respectively, below the mean radius of the target body and the threshold velocity of 20 m/s. The second term of the stopping criterion is a computational condition that imposes a maximum number of cost function evaluations Feval. Each simulation will stop if Feval exceeds the threshold value set for the specific mission. If this threshold is reached, the algorithm automatically reinitializes itself with a different initial distribution of particles. Both versions of the CBO will be used. In order to perform a more trustworthy analysis, one thousand runs will be executed for each case. Each run departs from different initial distributions of particles in the search space because of the randomness of the initialization of the positions (Equation (1)). This is done to evaluate the average performances of the algorithm. The first performance index considered is the percentage of success %succ. A success occurs when the accuracy condition of the stopping criterion is met. Other performance indices are minimum and average values of transfer time tf and number of function evaluations Feval. In addition, the dispersion around the mean value of transfer time σtf will be shown.


### 4.1. Rendezvous with Mars

Here, the outer planet Mars has to be encountered by a spacecraft with the following characteristics:(i)
ThrustT=300mN
(ii)
Specific impulseIsp=3000s
(iii)
Initial massm0=1000kgThe search space, in terms of angles, arc duration, and departure time, must be chosen to let the optimization process start. The optimizer will look for the best solution out of the bounds here empirically defined:(1)
αk∈0°,180°
(2)
ψk∈−90°,90°
(3)
∆tk∈10,100days
(4)
tdep∈01/07/2019−01/07/2020The number of arcs chosen isNA=4; therefore, there are 13 variables. The coefficients of the cost function are so established: Ct=0.01,Cd=10,Cv=10. The maximum number of function evaluations selected for this problem is 20000. It is used a basic version of the CBO by employing 50 particles to investigate the search space.As can be seen from Table1, the minimum time transfer, on average, lasts 302 days with a minimum of 288 days. The standard deviation is small; therefore, although the search space is wide, the CBO finds similar trajectories and this demonstrates the capability of the algorithm to explore effectively the search space. In Figure 7, the behavior of the cost function is plotted for all the successful runs (77.8%). It can be seen that the cost function decreases similarly for each run, which can be interpreted as a sign of robustness of the algorithm at changes in the initial conditions.Table 1
Mars rendezvous results.
%succ
tfmindays
tfmean
σtf
77.8
288
302.08
8.418Figure 7
Cost function behavior for the Earth-Mars transfer.The best trajectory found so far carries the spacecraft to Mars in 288 days with a final mass of 746.35 kg (Figure8). Finally, the control law is plotted in Figure 9.Figure 8
Optimal trajectory for the Earth-Mars transfer.Figure 9
Optimal control law for the Earth-Mars transfer.


### 4.2. Rendezvous with Venus

Here, the inner planet Venus must be encountered by a spacecraft with the same characteristics as the previous one.The search space, in terms of angles, arc duration, and departure time, is defined as follows:(1)
αk∈−180°,0°
(2)
ψk∈−90°,90°
(3)
∆tk∈1,90days
(4)
tdep∈09/01/2019−09/01/2020The number of arcs chosen isNA=5; therefore, there are 16 variables. After some preliminary tests, the coefficients of the cost function are so established: Ct=0.01,Cd=100,Cv=50. In this problem, the number of maximum function evaluations is set to 50000. It is employed a basic CBO with 80 particles. The results are shown in Table 2. The ability of the CBO to both explore and exploit is confirmed here by the small dispersion of the duration transfer around the mean value.Table 2
Venus rendezvous results.
%succ
tfmindays
tfmean
σtf
38.8
232.8
241
5.82The balance between exploration and exploitation phase can also be deduced by the almost linear slope of the cost function trend, plotted in Figure10. The optimal control law and the optimal trajectory are shown in Figures 11 and 12.Figure 10
Cost function behavior for the Earth-Venus transfer.Figure 11
Optimal control law for the Earth-Venus transfer.Figure 12
Optimal trajectory for the Earth-Venus transfer.


### 4.3. Rendezvous with Pallas

Pallas is an asteroid belonging to the asteroid belt that describes a very inclined orbit around the Sun (34.8°) with an eccentricity of 0.2305, making it a difficult body to reach. The probe chosen for this mission has the following characteristics:(i)
ThrustT=80mN
(ii)
Specific impulseIsp=3000s
(iii)
Initial massm0=600kgThe search space is as follows:(1)
αk∈0°,180°
(2)
ψk∈−90°,90°
(3)
∆tk∈5,250days
(4)
tdep∈01/01/2019−01/01/2021The number of arcs chosen isNA=11 resulting in 34 variables. The coefficients are Ct=0.01,Cd=50,Cv=100. The maximum number of function evaluations allowed for this problem is 100000. It is employed an ECBO with one particle in the Colliding Memory without considering the crossover. The number of particles is 120. Despite the biggest energetic difference between the Earth and target orbits, the strictest arrival tolerances, the least powerful thruster, and the highest number of variables, the CBO can find optimal solutions for this rendezvous. The success percentage of 23% (Table 3) can be addressed to the concepts just mentioned, yet the value of standard deviations of tmin is very small (3% of the mean transfer time), revealing that the discovered trajectories are very close to one another. This can be interpreted as a sign of the reliability of the algorithm. Figure 13, Figure 14, and Figure 15 represent, respectively, the optimal control law, the cost function behavior, and the optimal trajectory.Table 3
Pallas rendezvous results.
%succ
tfmindays
tfmean
σtf
23
1429.3
1510.1
47.02Figure 13
Optimal control law for the Earth-Pallas transfer.Figure 14
Cost function behavior for the Earth-Pallas transfer.Figure 15
Optimal trajectory for the Earth-Pallas transfer.


### 4.4. Discussion of the Results

It is worth mentioning that the direct approach, utilized in the last three cases, does not contemplate optimality conditions; therefore, there is not a definite optimal trajectory. Since there is no reference solution, a large number of simulations is needed to characterize the behavior of the algorithm. Due to the many degrees of freedom offered, the developed strategies allow many solutions to satisfy the arrival constraints. Although the range of possible transfer times is wide, the solutions found by the CBO are not equally distributed on all the possible values oftf. In Figures 16, 17, and 18, the transfer durations of the successful trajectories are grouped into histograms. For all the cases tested above, they are centered around a mean value, biased towards the minimum time transfer, with very small standard deviation values (equal or less than 3% of the mean value), reported in Tables 1, 2, and 3.Figure 16
Transfer time histogram for the Earth-Mars transfer.Figure 17
Transfer time histogram for the Earth-Venus transfer.Figure 18
Transfer time histogram for the Earth-Pallas transfer.


## 4.1. Rendezvous with Mars

Here, the outer planet Mars has to be encountered by a spacecraft with the following characteristics:(i)
ThrustT=300mN
(ii)
Specific impulseIsp=3000s
(iii)
Initial massm0=1000kgThe search space, in terms of angles, arc duration, and departure time, must be chosen to let the optimization process start. The optimizer will look for the best solution out of the bounds here empirically defined:(1)
αk∈0°,180°
(2)
ψk∈−90°,90°
(3)
∆tk∈10,100days
(4)
tdep∈01/07/2019−01/07/2020The number of arcs chosen isNA=4; therefore, there are 13 variables. The coefficients of the cost function are so established: Ct=0.01,Cd=10,Cv=10. The maximum number of function evaluations selected for this problem is 20000. It is used a basic version of the CBO by employing 50 particles to investigate the search space.As can be seen from Table1, the minimum time transfer, on average, lasts 302 days with a minimum of 288 days. The standard deviation is small; therefore, although the search space is wide, the CBO finds similar trajectories and this demonstrates the capability of the algorithm to explore effectively the search space. In Figure 7, the behavior of the cost function is plotted for all the successful runs (77.8%). It can be seen that the cost function decreases similarly for each run, which can be interpreted as a sign of robustness of the algorithm at changes in the initial conditions.Table 1
Mars rendezvous results.
%succ
tfmindays
tfmean
σtf
77.8
288
302.08
8.418Figure 7
Cost function behavior for the Earth-Mars transfer.The best trajectory found so far carries the spacecraft to Mars in 288 days with a final mass of 746.35 kg (Figure8). Finally, the control law is plotted in Figure 9.Figure 8
Optimal trajectory for the Earth-Mars transfer.Figure 9
Optimal control law for the Earth-Mars transfer.


## 4.2. Rendezvous with Venus

Here, the inner planet Venus must be encountered by a spacecraft with the same characteristics as the previous one.The search space, in terms of angles, arc duration, and departure time, is defined as follows:(1)
αk∈−180°,0°
(2)
ψk∈−90°,90°
(3)
∆tk∈1,90days
(4)
tdep∈09/01/2019−09/01/2020The number of arcs chosen isNA=5; therefore, there are 16 variables. After some preliminary tests, the coefficients of the cost function are so established: Ct=0.01,Cd=100,Cv=50. In this problem, the number of maximum function evaluations is set to 50000. It is employed a basic CBO with 80 particles. The results are shown in Table 2. The ability of the CBO to both explore and exploit is confirmed here by the small dispersion of the duration transfer around the mean value.Table 2
Venus rendezvous results.
%succ
tfmindays
tfmean
σtf
38.8
232.8
241
5.82The balance between exploration and exploitation phase can also be deduced by the almost linear slope of the cost function trend, plotted in Figure10. The optimal control law and the optimal trajectory are shown in Figures 11 and 12.Figure 10
Cost function behavior for the Earth-Venus transfer.Figure 11
Optimal control law for the Earth-Venus transfer.Figure 12
Optimal trajectory for the Earth-Venus transfer.


## 4.3. Rendezvous with Pallas

Pallas is an asteroid belonging to the asteroid belt that describes a very inclined orbit around the Sun (34.8°) with an eccentricity of 0.2305, making it a difficult body to reach. The probe chosen for this mission has the following characteristics:(i)
ThrustT=80mN
(ii)
Specific impulseIsp=3000s
(iii)
Initial massm0=600kgThe search space is as follows:(1)
αk∈0°,180°
(2)
ψk∈−90°,90°
(3)
∆tk∈5,250days
(4)
tdep∈01/01/2019−01/01/2021The number of arcs chosen isNA=11 resulting in 34 variables. The coefficients are Ct=0.01,Cd=50,Cv=100. The maximum number of function evaluations allowed for this problem is 100000. It is employed an ECBO with one particle in the Colliding Memory without considering the crossover. The number of particles is 120. Despite the biggest energetic difference between the Earth and target orbits, the strictest arrival tolerances, the least powerful thruster, and the highest number of variables, the CBO can find optimal solutions for this rendezvous. The success percentage of 23% (Table 3) can be addressed to the concepts just mentioned, yet the value of standard deviations of tmin is very small (3% of the mean transfer time), revealing that the discovered trajectories are very close to one another. This can be interpreted as a sign of the reliability of the algorithm. Figure 13, Figure 14, and Figure 15 represent, respectively, the optimal control law, the cost function behavior, and the optimal trajectory.Table 3
Pallas rendezvous results.
%succ
tfmindays
tfmean
σtf
23
1429.3
1510.1
47.02Figure 13
Optimal control law for the Earth-Pallas transfer.Figure 14
Cost function behavior for the Earth-Pallas transfer.Figure 15
Optimal trajectory for the Earth-Pallas transfer.


## 4.4. Discussion of the Results

It is worth mentioning that the direct approach, utilized in the last three cases, does not contemplate optimality conditions; therefore, there is not a definite optimal trajectory. Since there is no reference solution, a large number of simulations is needed to characterize the behavior of the algorithm. Due to the many degrees of freedom offered, the developed strategies allow many solutions to satisfy the arrival constraints. Although the range of possible transfer times is wide, the solutions found by the CBO are not equally distributed on all the possible values oftf. In Figures 16, 17, and 18, the transfer durations of the successful trajectories are grouped into histograms. For all the cases tested above, they are centered around a mean value, biased towards the minimum time transfer, with very small standard deviation values (equal or less than 3% of the mean value), reported in Tables 1, 2, and 3.Figure 16
Transfer time histogram for the Earth-Mars transfer.Figure 17
Transfer time histogram for the Earth-Venus transfer.Figure 18
Transfer time histogram for the Earth-Pallas transfer.


## 5. Conclusions

The study conducted in this paper reveals that the CBO is a valid algorithm, capable of solving a broad range of trajectory optimization problems. In the indirect cases, the CBO finds the optimal set of costate variables, discovering the optimal trajectory in a straightforward way. The direct cases were made challenging on purpose, to test the effectiveness of the CBO thoroughly. Despite the increasing number of variables, the large domains, and the strict arrival tolerances, the algorithm succeeded in achieving satisfying trajectories for all the problems presented, exhibiting also a small dispersion around the suboptimal solution. Concerning the CBO performances, it is possible to state that, although it is straightforward in terms of computational effort and parameter settings, it shows great reliability and effectiveness in solving these kinds of problems. The remarkable advantage of CBO is that it is ready to use and it does not need an initial guess of the solution or insights of the problem. In addition, CBO has just one parameter to adjust (Elitism scheme); therefore, it does not need the fine-tuning preliminary operations required by metaheuristics in general.



---

*Source: 9437378-2020-01-03.xml*