Hierarchical Optimization Algorithm and Applications of Spacecraft Trajectory Optimization

: The pursuit of excellent performance in meta-heuristic algorithms has led to a myriad of extensive and profound research and achievements. Notably, many space mission planning problems are solved with the help of meta-heuristic algorithms, and relevant studies continue to appear. This paper introduces a hierarchical optimization frame in which two types of particles—B-particles and S-particles—synergistically search for the optima. Global exploration relies on B-particles, whose motional direction and step length are designed independently. S-particles are for ﬁne local exploitation near the current best B-particle. Two speciﬁc algorithms are designed according to this frame. New variants of classical benchmark functions are used to better test the proposed algorithms. Furthermore, two spacecraft trajectory optimization problems, spacecraft multi-impulse orbit transfer and the pursuit-evasion game of two spacecraft, are employed to examine the applicability of the proposed algorithms. The simulation results indicate that the hierarchical optimization algorithms perform well on given trials and have great potential for space mission planning.


Introduction
In the last twenty years, the rapid development of meta-heuristic optimization algorithms has promoted the extraordinary progress of widespread engineering applications, including variable selection in chemical modeling [1], pattern recognition [2], path planning of UAVs [3], feature selection [4] and data clustering [5]. Doubtlessly, the powerful capabilities of parameter optimization are beneficial to space missions planning, such as spacecraft rendezvous trajectory optimization [6], interplanetary trajectories by multiple gravity-assist [7], agile satellite constellation design [8] and spacecraft attitude maneuver path planning [9]. In order to effectively apply the meta-heuristic algorithms to the space mission planning, on the one hand, constructing an appropriate optimization problem model is of vital importance. On the other hand, it is valuable to improve the performance of algorithms, which is the focus of this paper. Compared with traditional mathematical programming methods, meta-heuristic algorithms attract the attention of a large number of scholars due to four characteristics: simplicity, flexibility, derivation-free mechanism and local optima avoidance [10]. Relevant research can be divided into three types: proposing new optimization algorithms, improving existing algorithms and applying existing algorithms to practical problems. Apparently, new algorithms are most noteworthy for their new mechanisms to provide inspiration for other methods.
To design a highly applicable meta-heuristic algorithm, researchers have abstracted phenomena and laws in the real world into mathematical descriptions, from which a great number of algorithms were born. Genetic Algorithm (GA) [11] imitates Darwinian evolution theory and is a pioneer of multitudinous evolutionary optimization algorithms, including Evolutionary Programming (EP) [12] and Differential Evolution (DE) [13]. Typical vided by the local search of the second hierarchy, in which a certain number of S-particles iteratively exploit the neighborhood, and the best result is returned to update the position of the best B-particle. The number of function evaluations is a constant predefined and identical to the total number of generated S-particles. Two specific algorithms, HOA-1 and HOA-2, are designed based on HOF, and the difference lies in the local search methods of S-particles. In both algorithms, B-particles, except the best one, iteratively search the space by referring to the best two B-particles. In HOA-1, S-particles perform a pattern search by iteratively searching the space dimension by dimension, and the search step length in each dimension is a uniformly distributed random number with stepped decreasing bounds with the increase in the iteration number of first hierarchy. In each iteration of the second hierarchy, S-particles in HOA-2 are generated simultaneously based on a set of Gaussian distributions, which are iteratively updated according to the positions of the current best two function values.
To test the effectiveness of algorithms, the benchmark functions test is a popular method. However, some existing algorithms perform well in solving benchmark functions whose solutions are the origin of the search space, and an obvious performance degradation appears when the solutions deviate a little from the origin. Considering this unreasonable phenomenon, twenty-three variant benchmark functions are designed. In contrast to regular benchmark functions repeated solving, we set a random position deviation to the solution of the benchmark function, while keeping the minimum value invariant. Therefore, these variant benchmark functions were solved by the two proposed algorithms and seven existing algorithms repeatedly. Besides the benchmark functions test, two spacecraft trajectory optimization problems were solved by proposed algorithms and compared algorithms. The first problem was the spacecraft multi-impulse orbit transfer between two coplanar orbits, and the target was to find the optimal velocity impulse vectors and corresponding time. The second problem is the two-spacecraft pursuit-evasion game solved by the multiple shooting method, and the target is to find a set of appropriate costate variable initial values.
The rest of this paper is arranged as follows: Section 2 gives a statement about the optimization problem to be solved. Section 3 describes HOF and two algorithms systematically. The performance tests of proposed algorithms on benchmark functions and trajectory optimization problems are respectively presented in Sections 4 and 5. The conclusion is given in Section 6 ultimately.

Optimization Problem Formulation
Optimization problems cover a wide range of subproblem types, which can be singleobjective or multi-objective, unconstrained or constrained, static or dynamic. The focus of this work is to solve the static unconstrained single-objective optimization problem (SUSOP), since other types of problems can be studied by introducing more techniques and mechanisms into the achievement of SUSOP. Stochastic/heuristic optimization techniques have been extensively employed to handle optimization problems, and explicit reviews in [10,27] may help interested readers gain a comprehensive understanding of the development in this field.
This section gives a brief description and definition of the problem to be studied in the following work. Usually, a minimization problem is formulated to represent the optimization problem, as follows: Minimize : f (X)          X = (x 1 , x 2 , . . . , x N−1 , x N ) T lb i ≤ x i ≤ ub i , i = 1, 2, . . . , N Ub = (ub 1 , ub 2 , . . . , ub N ) T Lb = (lb 1 , lb 2 , . . . , lb N ) T (1) where lb i and ub i are boundary values of the ith element of state variable X due to the restrictions on an N-dimension search space in numerical calculation. The continuity of function f means that a best solution must exist in the search space, even though the global optima may not be located in this space, especially for those problems in which it is difficult to stipulate the search space, for instance, when it is necessary to search for the initial values of costate variables in two-point boundary problem solving.
To find the best solution of the SUSOP, abundant algorithms have been proposed. Meta-heuristic algorithms are one kind of method quite popular over recent decades. A typical feature of meta-heuristics is the random factor, which leads meta-heuristics to receive different results so as to increase the probability of finding the optimal solution. A meta-heuristic algorithm (MHA) can be formulated as follows: where f (X) is the optimization function with a given search space, X 0 is the initial values, {P} is the parameter set including constant and varied parameters, and X * is the solution of f (X) obtained by MHA. The performance of MHA is significantly influenced by the values of {P}. However, the focus of this work is to design the search strategy, while the values of parameters are intuitively set and manually adjusted by the simulation results in the following paper.

Hierarchical Optimization Algorithm
The structure of HOF is outlined in this section. Subsequently, two specific algorithms, HOA-1 and HOA-2, are provided and their differences are discussed.

Hierarchical Optimization Frame
In this paper, we used particles to represent the individuals in SI algorithms. For SI algorithms, a popular idea is to use the position of the current best particle (CBP) to guide other particles' motion, which makes the exploration and exploitation in the search space more effective and efficient. However, CBP itself seldom obtains benefit from this strategy. To alleviate this limitation, HOF was proposed to enhance the search ability of CBP by performing a local search and giving it additional function evaluations. HOF includes two hierarchies: H1 is the first hierarchy and H2 is the second hierarchy. In H1, a certain number of particles called B-particles are used to iteratively search for the global optimum. Therefore, the CBP in HOF is the current best B-particle. In each H1 iteration, B-particles cooperate to iteratively update their positions based on the current information of search results. Generally, only one function evaluation is conducted to complete the position update of each B-particle except CBP, whose position update depends on the local search of H2 iterations and requires more function evaluations. In H2, particles of another type called S-particles, initially distributed in the neighborhood of CBP, are assigned to iteratively exploit the space. After a predefined constant number of S-particles are randomly generated, the position of the best S-particle becomes the updated position of CBP based on the greedy strategy. The iterations of two hierarchies are executed in turn until the maximum number of H1 iterations is reached. Note that B-particles converge as the H1 iterative search progresses, while S-particles dispersedly exploit the local space from the initial position of CBP. I b and I s are the maximum number of iterations in H1 and H2, which means that I s times H2 iterations proceed after each H1 iteration, and the stopping criterion is I b times H1 iterations have been performed. The total number of function evaluations is determined according to the iteration numbers I b and I s , the B-particles number and the S-particles number. The whole frame is modularized and illustrated by Figure 1, in which the orange blocks represent adjustable sub-algorithms. A 2-D sphere function was used to show the CBP's update process in one H1 iteration with its subsidiary H2 iterations. In Figure 2, Figure 2a shows the positions of eight B-particles, in which the CBP is represented by a red marker, while other B-particles are blue markers. Four S-particles, marked by pentagrams, were used to perform the local search in the neighborhood of the CBP, and their initial positions and final positions are respectively given in Figure 2b,c. The best S-particle was selected to update the CBP's position after H2 iterations, and is denoted by the green marker in Figure 2d.  There is no specific formula given in the flowchart in Figure 1. To design specific algorithms, we designed an iterative method for B-particles used in this paper, while other methods are also applicable. In H1, we culled the best two B-particles in each iteration and set them as two leaders, L1 and L2, to guide the motion of other B-particles in the next iteration, except for L1's position, which was determined by the result of H2 iterations. Each SI algorithm has its own way of deriving particle propagation, including two aspects, the direction and the step length. The velocity computation of PSO provides a particle's motional direction and step length simultaneously. In comparison, the direction and step length of bacteria are computed respectively in Bacterial Foraging Algorithm [30], in which the direction is a unit vector randomly produced in the searching space, and the step length is a small predetermined value to meet search accuracy. Both algorithms have respective advantages and limitations, which inspired us to separate the determination of searching direction and step length to strengthen the algorithm's flexibility and adaptability for different optimization problems.
Since the search direction and step length are calculated respectively, we used positions of L1 and L2 to derive a particle's direction update formula by where B n i indicates ith B-particle's position in nth iteration in H1, L1 n and L2 n represent the positions of L1 and L2 in nth iteration respectively and the parameters c 1 and c 2 control the weights of L1 and L2 to B n i . The symbol · means the Euclidean norm of a vector. With the inconsistency of L1 and L2, the denominator of D n i is always larger than 0, which averts the singularity. Assuming L1 n exerts more influence on B n i than L2 n , we set c 1 > c 2 > 0. The step length of B n i is given by where α is the step ecoefficiency between 0 and 1 to directly tune the searching range of B n i . In total, the updating equation of B n i in H1 iteration is constructed as follows: note that the update mechanism of B-particles is deterministic, due to the constant values of c 1 , c 2 and α, which is not best, but effective. Further study on the autonomous adjustment of these parameters can be conducted by introducing randomness or adaptive laws. L1 n is actually the position of CBP in nth H1 iteration. Different ways to generate and update S-particles can construct many specific algorithms, one of which may better solve the focused problem. Two hierarchical optimization algorithms are introduced in the next subsection, and their difference lies in the local search method of S-particles.

Two Hierarchical Optimization Algorithms
Considering the performance differences caused by the different searching methods of S-particles, two hierarchical optimization algorithms, HOA-1 and HOA-2, are proposed to validate the effectiveness of HOF.

Formula of HOA-1
The searching method of S-particles can be divided into two steps: initialization and iterative update. The initialization is to determine the spatial domain and the position distribution law. The first step is to set a spatial domain of initial S-particle positions. A gradually shrinking domain makes a fine effect, because S-particles are used to search the neighborhood of CBP. In order to simplify the algorithm, we designed a function to control the boundary of this domain in Equation (6).
where F(·) is an operator to tune the range of deltaS n in nth H1 iterations, Ub and Lb are the boundary of the search space defined by (1), and T is a parameter to control the phased variation of search precision. In HOA-1, F(·) is defined in (7).
The symbol · indicates rounding down the variable. This function means the domain stays invariant for T times H1 iterations, and shrinks to one tenth in the next T iterations, and c is a constant used to determine the initial search precision. Generally speaking, we suppose that a better position exists within a hypercube centered on CBP, and take deltaS n to set the length of its edges.
In HOA-1, the number of S-particles is equal to the dimension N of the problem to be solved, which results from the idea that changing the value of CBP's position dimension by dimension serves to precisely search the space around CBP.
In H2 iterating process, the reference point P r is used to update the positions of Sparticles, and the initial P r is the CBP of the B-particles in the current H1 iteration. Once a better position is found, it becomes the new P r of the next H2 iteration. Let S m j be the position of jth S-particle in mth H2 iteration. The updating equation of S m j is constructed as follows: S m j = P r + r s · I N,j · deltaS n (j), m = 1, 2, . . . , I s (8) where I N,j is the jth column vector of a N-order identity matrix I N , deltaS n (j) is the jth element of deltaS n , and r s is a random number uniformly distributed in the interval of [-1,1]. After H2 iterations, the final P r is output to renew the position of CBP if its function value is smaller. The HOA-1 algorithm is summarized in Algorithm 1. In HOA-1, the way S-particles search for a better solution is by changing only one element of P r in a given domain to generate a new S-particle that is uniformly distributed in the hypercube. To design a different search strategy, a new algorithm HOA-2 is proposed in which the number of S-particles can be specified arbitrarily, every element of each S-particle is obtained from a Gaussian distribution, and all elements are changed simultaneously. To determine a suitable Gaussian distribution, the mean value µ and the standard deviation σ should be defined. Thus, the positions of L1 and L2 are input to the H2 optimization, and µ and σ can be given as follows: where µ 1 n and σ 1 n represent the first µ and σ in the H2 optimization of nth H1 iteration, respectively. The operator abs(·) is used to replace all elements of the vector with their absolute values. This specification stems from the assumption a better position exists closer to L1 than L2. Given that the range of this parameter setting is sometimes too large, another specification in Equation (10) is proposed, which is also subject to the above assumption.
The difference between Equations (9) and (10) is illustrated in Figure 3, in which the curves are the probability distribution density functions of two ways. It shows the interval constrained by Equation (9) is obviously larger. Due to the high probability of the generation of new S-particles between L1 and L2 by using Equation (10), this specification is preferred when the best solution is encircled by S-particles. To guarantee the validity of this estimation for better solutions, µ and σ are updated after every H2 iteration, which leads to the necessity of selecting and updating the best two S-particles, namely SL1 and SL2. Using SL1 and SL2 to severally replace L1 and L2 in Equations (9) or (10) can update µ and σ in each H2 iteration.
When the iterations in H1 proceed, if no other B-particle takes over these two particles, the difference between L1 and L2 becomes small, which makes the search space in H2 too narrow for S-particles to help CBP find a better position. With that in mind, a lower bound for σ is introduced to prevent S particles from being limited to a cramped space search. As with the restriction on deltaS n , Equation (11) gives the definition of σ and σ min .
where σ m n is σ in mth H2 iteration of nth H1 iteration, and SL1 m and SL2 m are positions of SL1 and SL2 in mth H2 iteration, respectively. Therefore, based on the setting of µ and σ in Equation (9), the updating equation of jth S-particle S m j in mth H2 iteration is constructed as follows: where the subscript k means kth element of a vector and σ m n (k) is kth element of σ m n . Likewise, the setting of µ and σ based on Equation (10) can be derived, but is omitted here. The HOA-2 algorithm is summarized in Algorithm 2. To maintain the relative independence between H1 and H2, only SL1's final position is output to update CBP's position. HOA-1 and HOA-2 share the same H1 and differ from the search mechanism in H2, which manifests more optimization algorithms that can be proposed based on HOF. By introducing more mechanisms into HOF, such as evolutionary operators, simulated annealing mechanism and other effective mechanisms, a new algorithm with better performance can be accomplished.

Comparison with Other Algorithms
The concept of hierarchical optimization has been used in algorithm design and some results have been employed in single-objective optimization and multi-objective optimization. In [31], a hierarchical algorithm (HGA-PSO) was developed to construct two-layer optimization, in which the bottom layer, responsible for exploration, is based on GA, and the top layer, for exploitation, adopts PSO. The bottom layer includes many subgroups searching for the best solution independently. After each iteration in the bottom layer, some excellent individuals of every subgroup are brought together to form a new group to continue cooperative search in the top layer. Until the termination conditions are met, k individuals are selected randomly from the top layer to substitute random k individuals of every subgroup in the bottom layer before the next iteration starts. An external archive is established to store the results of the top layer search, and subgroups choose k individuals from the archive [32], which eliminates the negative influence of initialization and randomness of individuals in the top layer by outputting the final solution from the archive. Compared with the hierarchical optimization methods mentioned above, the HOF is based on a competition mechanism, which is analyzed in detail in Section 4.2.2 through simulation examples, instead of combining the best individuals from multiple subgroups to conduct a cooperative optimization. In [33], a multiagent collaborative search (MACS) method was proposed for local exploration of the subdomain generated by the branching method. Multiple agents were generated in the given subdomain and used to explore this subdomain by elaborate collaboration mechanisms. Afterwards, MACS was further developed based on a combination of Tchebycheff scalarization and Pareto dominance to solve multi-objective optimization problems [34]. Compared with MACS, the exploitation of S-particles around CBP is quite different. The most obvious difference is that S-particles locally search the space around CBP independently, while strong information interaction exists between agents in MACS. Also, a search subdomain is firstly determined, and then the population of agents is introduced and updated to search this invariant subdomain, while the local search space of S-particles varies with the position update of the current optimal solution.
Regarding the H2 optimization as the local search of H1 optimization, the frame of HOF is similar to memetic algorithms. Therefore, a further comparison of two proposed algorithms and existing algorithms is conducted from the point of view of the local search. An adaptive gradient descent-based local search was introduced in MA [35], in which the numbers of individuals considered in local search, the search steps and iterations decreased linearly during the iterations. However, the gradient calculation requires additional function evaluation to generate a new individual, which leads to the important task of determining the ratio of the function evaluations of gradient calculation to the total evaluation number. In HOA-1 and HOA-2, the local search of S-particles proceeds without gradient information, and accordingly, once function evaluation occurs when generating a new S-particle. The local search mechanism of HOA-1 can be classified as a pattern search method with random factors. The stepped decreasing upper bound of iterative step length of S-particles in each dimension is related to the iteration number of H1, and differentiates HOA-1 from other local pattern search methods, for example, the memetic algorithm combining PSO and pattern search refinement [36], and multi-coordinate search [37]. This stepped decreasing mechanism of search step guarantees the search ability of CBP, even if all B-particles converge in the early iterations of H1 and avoids the local minima.
Another point worth demonstration concerns the Gaussian distribution used in HOA-2. Kennedy proposed a Bare Bones Particle Swarms Optimization (BBPSO) in 2003 [38], in which every particle's position update is obtained by the Gaussian distribution defined by the best particle position and its historical optimal position. Then, many researchers improved the mechanism of Gaussian distribution parameter setting to increase its searching efficiency [39,40]. In these algorithms, the weight of best individual in defining the Gaussian distribution parameters is identical to other individuals, due to the balance of exploration and exploitation during the iterations. In HOA-2, the weight of L1 or SL1 is determinately bigger than that of L2 or SL2, because it adopts the hypothesis that the solution is closer to the best particle than other particles. Understandably, this hypothesis leads to closer distances between sampled particles and the best particle, and thus causes the premature convergence. Therefore, this hypothesis has a negative effect on the exploration and is rejected by existing algorithms. However, only the best one of sampled S-particles is selected to update the position of the CBP, and other S-particles have no influence on the iterative exploration of B-particles, which allows an increase in the weight of the best particle. In HOA-2, the S-particles are simultaneously sampled based on the same set of Gaussian distributions in each H2 iteration, and these Gaussian distributions are updated according to the positions of two best S-particles Sl1 and SL2. However, each particle in BBPSO has its own set of Gaussian distributions, which means that only one particle is sampled based on these distribution functions. The multiple sampling based on the Gaussian distributions can better exploit the neighborhood of the best particle, while the randomness of single sampling brings a negative effect to the local search. The local search methods based on covariance matrix adaption also utilize multiple dimension Gaussian distribution to describe the sampling subdomain around the center of selected individuals [41]. However, this type of algorithm emphasizes the adaptive laws of the covariance matrix and step size, without considering the weights of dominant individuals.

Experiments on Benchmark Functions
In this section, twenty-three benchmark functions are used to validate the performance of HOA-1 and HOA-2 compared with other optimization algorithms. Then, analysis and discussion about the results are elucidated.

Explanations of Performance Test
For comprehensive inspection of optimization algorithms, twenty-three benchmark functions, which consist of three types of functions-unimodal, multimodal and fixeddimension multimodal-were utilized in [10]. To better examine the performance of proposed HOA-1 and HOA-2, we designed variants of these twenty-three benchmark functions instead of directly using them to testify the overall effectiveness. The origin benchmark functions are listed in Tables A1-A3 in Appendix A. We introduced skills of variable substitution to shift the theoretical solution of each benchmark function in the searching space, while keeping the minima unchanged. The variable substitution can be expressed as ] which covers 1/20 of the range of ith dimension search space. Taking F1 sphere function as an example, the actual fitness function is , 100], and r i is a uniformly distributed constant in [-5, 5], which means the extreme point is no longer the coordinate origin but a random point near the origin. Solving the minima of F1 function repeatedly, the theoretical value is always 0, while the position is R = [r 1 , . . . , r 30 ] T , which is determined by the results of random sampling before the optimization starts. In this way, the efficacy of algorithms emerges more clearly.
In contrast to the proposed algorithms, seven algorithms were introduced as comparison algorithms, including three classical algorithms: GA, DE and PSO, and four recently proposed algorithms: GWO, Butterfly Optimization Algorithm (BOA) [42], HHO and AOA. The parameters of these algorithms are defined in Table 1. The numbers of iterations I max of all these contrastive algorithms were all set to 500, and the population size was 30. Since two hierarchies exist in HOA-1 and HOA-2, we set I b to 100 and I s to 4 for the sake of fairness, which guarantees a smaller number of function evaluations in HOA-1 and HOA-2 than in those of the compared algorithms. Because different benchmark functions require different search precision, the parameter T in HOA-1 and HOA-2 is defined correspondingly.

Algorithm Parameters
GA These algorithms were executed 30 times to solve each benchmark function on a laptop with an Intel(R) Core (TM) i7-10510U CPU at 1.8 GHz and 16.0 GB of RAM. All these algorithms are coded in MATLAB R2020a.

The Performance
The results of previous experiments are arranged meticulously, and the data are displayed in Tables 2-4, in which the best results obtained by nine algorithms for each benchmark function are presented in bold. In these three tables, B, A and S respectively stand for the best result, the average value and the standard deviation. The results of benchmark function test are given in Tables 2-4, in which the minimum values of these benchmark functions are bold. In the unimodal benchmark functions test, HOA-1 and HHO both found the best solution of three of the seven functions, and showed outstanding exploitation ability, according to Table 2. Furthermore, the results of the multimodal functions test in Table 3 illustrate that the exploration ability of HOA-1 is superior to other algorithms. For local minima avoidance, the data in Table 4 demonstrate that all these algorithms performed similarly, while HOA-1, DE and PSO had a weak advantage on the others. To statistically display the performance of tested algorithms, the Friedman test [43] was used to reflect the performance, and the result is listed in Table 5, which proves the advantage of HOA-1 over other algorithms. The p-value of the Friedman test was 8.8089 −15 , suggesting the existence of significant differences among the tested algorithms. Generally, HOA-1 outperformed other algorithms in the benchmark functions tests. Compared with HOA-1, HOA-2 performed better in F15 than the others. However, it should be noted that HOA-2 is also an effective algorithm, because the mean rank of HOA-2 was fifth among the nine algorithms according to the Friedman test result in Table 5. The reason HOA-1 performs better than HOA-2 in benchmark functions test can be explained by two factors. The first reason is the small number of I s which makes it difficult to fully embody the advantages of Gaussian distribution multiple dimensions sampling in HOA-2. The second reason is that the landscapes of benchmark functions are simpler than the actual engineering problems. Therefore, a one-dimensional position update used in HOA-1 can efficiently search the space.
The history of minimum-value searching of benchmark functions is diagrammed in Figure 4. The gentle decrease in the fitness value means the exploitation search was effectively conducted, while the sharp or even vertical decline reflects a better result found by the exploration search. From this point of view, the exploration abilities of GA and HHO are outstanding, and the exploitation ability of HOA-1 is conspicuous. The stepped attenuation of search step length of S-particles attributes more to the descent slope changes in HOA-1 and HOA-2. Note that the number of function evaluations, which is equal to the number of individuals generated in the search space, is used as the coordinate of x-axis. Accordingly, the coordinates of the x-axis of the final minimum values of HOA-1 and HOA-2 are not the same as other compared algorithms, due to the smaller function evaluations number, which is amplified and shown in the subplot of F1 function in Figure 4.
To dissect the search process and convergence of two proposed algorithms, we selected F1, F8 and F14 to represent three types of functions with which to analyze the characteristics of iterative search. The fitness history curves of solving three functions by HOA-1 and HOA-2 are delineated in Figure 5. The fitness history curves show that B-particles in HOA-1 converge faster than those in HOA-2, while B-particles in HOA-2 have better potential for local minima avoidance, which is reflected by the intermittent sharp decline during the iterative search in Figure 5b,c. The initial and final positions of all B-particles of the search in HOA-1 and HOA-2 are displayed in Figures 6 and 7, respectively. Note that the positions of B-particles are represented by the first two coordinates. The B-particles gathered in a subdomain of the search space after the iterative search, which demonstrated that both algorithms possess the ability to drive B-particles to converge and find a high-precision solution.

The Mechanism Explanation by Examples
The core idea of HOF is based on a competition mechanism that tilts computing resources to CBP to make it better, and compels all B-particles to compete for the rewards assigned to the CBP. There is one point to be emphasized: in both HOA-1 and HOA-2, the reward for CBP is granted before the corresponding H1 iteration starts. In other words, once a B-particle takes the title of CBP after one H1 iteration, the reward is conferred on it regardless of whether the alternation of CBP appears in the next iteration. This strategy guarantees the execution of one H2 optimization in every H1 iteration, and is named Normal-Mode, which implies that the number of H2 optimization implementation may not be fixed to 1. There are another two possible situations. One case, called Conservation-Mode, is when H2 optimization is implemented when the ith B-particle becomes CBP and keeps the place until its next personal H1 iteration comes. If other B-particles assume the place of CBP in this period, the ith B-particle will not gain extra assistance from Sparticles. The other case is defined as Redundance-Mode, in which H2 optimization is executed as soon as the ith B particle becomes CBP. In this case, there may be more than one H2 optimization in one H1 iteration. These three situations reveal the time delay between the reward of CBP winning and awarding. From the perspective of computational stability, we embed Normal-Mode in HOF to create HOA-1 and HOA-2. By the phenomenon above, this competition mechanism strongly stimulates the competition in B-particles, especially when the numbers of iterations and S-particles in H2 is increased, which exacerbates the imbalance between the CBP and other B-particles.
To make a specific impression on readers, simulated examples were carried out to provide more concrete specification. For better demonstration of the competition mechanism, we used HOA-1 to find the minimum of F8 function and set I b to 150 and I s to 3 in order to weaken the preponderance of CBP. To show the competition for the CBP place among B-particles, we adjusted the parameter T to 10 to preferably exhibit the stepped attenuation of search range.
The situation of CBP alternation is shown in Figure 8, in which Figure 8a gives the alternation of CBP in each H1 iteration and Figure 8b presents the actual execution of H2 optimization in the local region of the red square in Figure 8a. The differences between three situations can be easily observed: in Normal-Mode the H2 optimization is carried out in every H1 iteration, while the number of H2 optimization in Conservation-Mode and in Redundance-Mode decreases five times and increases five times, respectively (the number is 10 in the whole process). Note that in this example, at most, two B-particles become the CBP in one H1 iteration. If the competition is fiercer and more B-particles sequentially become the CBP in one H1 iteration, another example will reflect the result in Figure 9, in which the numbers of H2 optimization carried out in three situations are respectively 16, 14 and 20. In Conservation-Mode, H2 optimization is implemented only when one B-particle remains as CBP for at least two successive H1 iterations. For the performance comparison of different algorithms, Normal-Mode is preferred, due to its fixed amount of computation.  Another typical mechanism is the stepped attenuation of deltaS n , which is shown in Figure 10. Because the upper bound and lower bound of each dimension of F8 function is the same, the first element of deltaS n is qualified to embody its range variation. All elements of deltaS n are reduced to one tenth every 10 H1 iterations. However, it was found in the experiment that the algorithm may not converge to the global optima if the value of T is too small, while a T that is too large causes a loss of solution accuracy in a finite number of iterations. Therefore, further research is warranted to resolve this contradiction.

The Value of Parameter T
In the benchmark functions test, a significant feature of both HOA-1 and HOA-2 is that the parameter T is usually different and tuned according to the test function, which is due to the shrinking mechanism of search space. When the number of iterations in H1 is fixed, decreasing the value of T means to accelerate the process of space shrinking. If this rate is too fast, the minimum search by S-particles will be severely restricted and limited in its ability to help CBP update its position. Only when the optimization in H2 works properly can HOA-1 and HOA-2 seek out a good solution. Conversely, a slow shrinking rate reduces the search efficiency, due to the search space being too large. Meanwhile, a finite iterative number inevitably limits the amount of computation. Hence, in order to promote the smooth progress of optimization, we used different constant values to set the parameter T to control the search precision.
For better elucidation of the influence of different T on a given function, an additional example was provided. Without loss of generality, the F1 function was used to present the phenomenon of performance deterioration and the HOA-1 was assigned to execute the optimization mission. In HOA-1, we added I s from 4 to 10 and set the value of T as 4, 5, 6 and 12. These four group results were collected and compared to show the difference. Each group had 30 solutions, ranked in descending order of F1 function, and all results are displayed in Figure 11, in which Fbest represents the minimum after each run, and the common logarithm of Fbest is taken to conveniently reflect the difference. Three key features can be concluded:

1.
A smaller T has the potential to make the algorithm find a better solution.

2.
A destabilization of the search may occur when the value of T varies.

3.
A T that is too small may impede the algorithm search for the optimal solution. In this work, the value of T was set by trial and error to achieve a balance between stability and search efficacy, which leaves room for further research to design an adaptive law for parameter adjustment and to mitigate the trouble of manual parameter setting.

Applications of the Proposed Methods
In this section, two spacecraft trajectory optimization problems, the multi-impulse minimum fuel orbit transfer and the pursuit-evasion game of two spacecraft, are employed to reflect the capability of the proposed HOA-1 and HOA-2. The former problem is extensively solved by meta-heuristic algorithms [33,44]. The latter problem considers the requirement of tracking a non-cooperative target with maneuvering ability, and has been the subject of much recent attention [45]. The functions of both problems demonstrate strong nonlinearity and variable coupling, even without additional constraints.

Problem Formulation and Parameter Setting
Spacecraft orbit optimal transfer by impulsive thrusters is a classical and significant problem. However, an analytic solution of optimal impulsive control for orbit transfer in complicated situations is difficult to obtain when the number of impulses increases and additional complicated constraints are considered. Therefore, a direct approach was developed to construct a general nonlinear programming problem and solve it by advanced optimization algorithms [46]. Relevant research results can be referred to in [44].
Considering the goal of proving the effectiveness of the proposed methods, we chose minimum fuel transfer between two coplanar circular orbits as the problem to be optimized. This has an analytical solution known as the famous Hohmann transfer with two-impulse maneuver. Usually, the fuel consumed is positively correlated with the velocity increment, so the orbit transfer by minimum velocity increment is equivalent to the minimum fuel case. For additional details on the derivation, [47] is recommended. In accordance with the size limit of this paper, a brief introduction of Hohmann transfer is generalized below.
For two coplanar circular orbits around the celestial body, an ellipse connecting two orbits and tangent to them is shown in Figure 12. A spacecraft on the initial orbit can reach the final orbit after two-impulse maneuver ∆v1 and ∆v2. The transfer trajectory is half an ellipse from the perigee to the apogee of the transfer elliptic orbit. Two impulses ∆v1 and ∆v2 are given as follows: where µ is the gravitational constant of the celestial body, r1 and r2 are radii of the initial and final circular orbit. According to the characteristic of Hohmann transfer with the theoretical solution of two-impulse maneuver, we designed three cases to test the optimization performance of the proposed two algorithms. As with problem 1 in [44], they concerned two orbits around Mars (µ = 42, 830 km/s 2 ), with r1 and r2 being 8000 km and 15,000 km, respectively. The three cases are distinguished by search space as listed in Table 6. Table 6. Search space of three cases.

Case Search Space
Case 1  Without loss of generality, the start points of the three cases were all set to the (8000,0). In Case 1, we gave the first in-plane impulse to make the spacecraft leave the initial orbit and applied the second impulse to make it revolve on the final orbit when it reached one point of the final orbit. When the first impulse was provided, the second impulse, if the transfer trajectory intersected the final orbit, was able to be calculated. Therefore, the search space was a 2-dimensional vector whose components were the radial and tangential velocity increments. The landscape is depicted in Figure 13, in which ∆v x and ∆v y represent coordinate components of ∆v 1 , and the domain of the initial velocity impulse that cannot realize orbit transfer is marked with a red dot. Subsequently, we shrank the range of initial velocity increments so as to make the transfer impossible by two-impulse, by which a threeimpulse orbit transfer was constructed and the search space became a five-dimensional vector, added to the second 2-dimension impulse and the applied time. Different values of the upper bound of the applied time led to Case 2 and Case 3. In Case 2, the upper bound of applied time allowed the spacecraft to fly one revolution and return to the start point, by which the same total impulse increment with two-impulse Hohmann transfer could be achieved by three-impulse maneuver. However, the upper bound of applied time violated this condition in Case 3, making the result of three-impulse orbit transfer by minimum velocity increment different from Case 2.
DE, GA, PSO, GWO, HHO, AOA, HOA-1 and HOA-2 were used to solve the three cases according to the result of benchmark functions test. The numbers of iterations Imax of all these contrastive algorithms were all set to 400, and the population size was 50. Considering the inconsistency between the dimension of search space and the number of search individuals of these algorithms, we set both I b and I s to 30. The number of S-particles in HOA-1 is equal to N, while that can be changed to 4N in HOA-2 to ensure that the total number of function evaluations is the same as the comparison algorithms. The parameters of DE, GA, PSO, GWO, HHO and AOA were the same as those in Table 1. For more meticulous search, we set α to 0.1 and c to 0 in both HOA-1 and HOA-2.

Results and Discussion
For Case 1 and Case 2, adding in ∆v1 and ∆v2 according to Equations (13), the theoretical minimum velocity increment is 0.6091531 km/s 2 . Each algorithm was run 30 times and the best result have been selected and listed in Table 7, in which the best solutions found by these algorithms are bold. The results show that HOA-2 found the best solutions for the three cases, which evidently reflects the superiority of the proposed methods. Because the number of S-particles of HOA-1 should be equal to the problem dimension N, and that of HOA-2 can be properly set to 4N, the amount of computation of HOA-1 is noticeably less than other algorithms, which explains the performance degradation of HOA-1. The process of orbit transfer using the best solution of three cases is diagrammed in Figure 14, which is consistent with the previous analysis and proves the correctness of the results. Comparing the results of Case 1 and Case 2, it can be observed that the difficulty of finding the optimal solution increases when the dimensions of the search space increase. In a similar way, further reduction in velocity range can construct the orbit transfer by more impulses and increase the difficulty of spacecraft orbit maneuver optimization. Nevertheless, uncertain spacecraft trajectory is more likely with the increasing number of impulses [48].

Problem Formulation and Parameters Setting
Spacecraft pursuit-evasion is a typical game problem, in which the pursuing spacecraft attempts to capture the target while the evading spacecraft strives for escape. Accordingly, this pursuit-evasion problem is a zero-sum differential game. As a result of the great difficulty of an analytical solution derivation, the numerical methods have been extensively used in recent decades. The indirect method of solving this problem is to derive the necessary optimal conditions, transform it into a two-point boundary value problem (TPBVP) and find a saddle point by numerical methods. In this paper, the TPBVP solution is used to test the performance of the proposed methods. The dynamic model and necessary conditions derivation are adopted from [49], and a transcription is provided below.
The differential equations of the pursuit-evasion game are written as follows: .
where X = [x P , y P , z P , .
x P , . y P , z E ] T is the state variable of positions and velocities of two spacecraft in the LVLH coordinate system of a virtual spacecraft in circular reference orbit, and the subscripts P and E represent the pursuing spacecraft and the evading spacecraft, respectively. U P = [0 1×3 , u Px , u Py , u Pz , 0 1×6 ] T and U E = [0 1×9 , u Ex , u Ey , u Ez ] T represent the direction vector of control exerted on two spacecraft, and satisfy the restraint conditions u 2 Ix + u 2 Iy + u 2 Iz ≤ 1, in which I means P or E. T P and T E are the value of maximum thrust per unit mass of two spacecraft. Clohessy-Wiltshire equations [50] are used to describe the relative motion of two spacecraft with respect to the reference rotating coordinate system, and the matrix A(t) can be expressed by Equation (15).
where ω is the constant rotational angular velocity of reference coordinate system in the gravitational field of the Earth. Set the terminal time as t f , so the payoff function is as follows: Therefore, the spacecraft pursuit-evasion game can be described thus: the goal of pursuing spacecraft is to change U P to minimize J, while the evading spacecraft aims at maximum by J by U E .
The adjoint variable λ = [λ P , λ E ] T is introduced corresponding to X, in which z ] T and I = P, E. Construct the Hamiltonian function as H = λ T X, and the control of two spacecraft can be derived according to the minimum principle and Cauchy-Schwarz inequality.
The transversal condition is derived by (19).
So far, the TPBVP has been derived completely. Input the guess of initial adjoint variables to the optimization algorithms and minimize J to test the performance of the proposed methods. Relevant parameters are set as follows: A reasonable initial value of adjoint variable λ 0 is given in (21), and the corresponding J is 0.489126783247148.
According to λ 0 , two cases were designed, as displayed in Table 8, to test the performance of the same eight algorithms as in Section 5.1. Note that in Table 8, [1] 12 is a 12-dimension column vector with all elements being 1. The parameters of these algorithms are almost the same as in Section 5.1, except that we changed I b to 40 for HOA-1 to improve its performance, and the number of S-particles of HOA-1 and HOA-2 were respectively set to 12 and 20, to ensure a smaller amount of computation compared with other algorithms.

Results and Discussion
For both cases, each algorithm was run 10 times and the best result was selected and listed in Table 9, in which the best solutions are bold. The results show that HOA-2 performed best in this test, while the performance of HOA-1 is poorer than GWO and HOA-2, but better than other algorithms. Considering the difference between the two cases, it can be concluded that Case 1 examined the exploration of algorithms, while exploitation was more important in Case 2. Thus, the result is consistent with the previous statement that hierarchical optimization is designed to assign more computation resources to the CBP position update. It is universally accepted that the solution of TPBVP is difficult to acquire due to the strongly sensitive characteristic of the initial guess value of the adjoint variable that lacks physical meaning. However, a better initial guess may help us solve the TPBVP by the multiple shooting method, whose limitation lies in the high requirement of initial value. To illustrate the diversion caused by poorer initial adjoint variables, an example is provided in which two initial values of adjoint variable are used as input to the multiple shooting method to check whether the evading spacecraft could be captured by the designed control (17). The first initial value is λ 0 , and the second one λ 0 is obtained by decreasing the first element of λ 0 by 20%. The payoff function J with λ 0 is 36.448. The result is diagrammed in Figure 15, which shows that a tiny change in the initial value of TPBVP had a great impact on the result. Therefore, the proposed methods can effectively help to guess the initial value of TPBVP. In these two spacecraft trajectory optimization problems, HOA-2 performed better than HOA-1. This can be explained by two factors. The first is that the function evaluation number of HOA-1 is significantly less than other algorithms. The second factor is that the S-particle multidimensional sampling, based on the Gaussian distribution function in HOA-2, performs more efficiently when I s increases. Therefore, the proper allocation of the numbers of B-particles and S-particles can make better use of hierarchical optimization algorithms. We suggest that HOA-1 is suitable when the problem dimension is relatively large, and otherwise HOA-2 is preferred, especially when the S-particles can be set to a much larger number than the problem dimension.

Conclusions
For improving the search ability of the best individual compared with others, a hierarchical optimization frame is proposed, in which B-particles synergistically explore the search space, and a local search performed by S-particles is conducted to update the position of the best B-particle. Based on this framework, two algorithms (HOA-1 and HOA-2) were designed. The local search in HOA-1 is a type of pattern search with random factors, and that of HOA-2 is based on the sampling in the iteratively updated Gaussian distributions. Considering the limitations of regular benchmark functions, twentythree variant benchmark functions were designed and solved to accurately reflect the performance of the proposed algorithms and compared algorithms. The experiment results show that HOA-1 outperforms other algorithms in the benchmark function test. In order to further verify the feasibility and superiority of the proposed algorithms, two spacecraft trajectory optimization problems, multiple-impulse orbit transfer and spacecraft pursuitevasion game, were introduced, due to their complexity and challenges. The HOA-2 algorithm excels outstandingly in solving trajectory optimization problems. In HOA-1, the number of S-particles should be identical to the dimension of the problem, which partially limits its performance and adaptivity.
The parameter T greatly influences the performance of proposed algorithms, especially on the convergence speed. Generally, a larger T is suitable to solve the unimodal function, while a smaller T is preferred in solving the multimodal function. Also, the allocation of function evaluation numbers between the B-particles and S-particles is important to determine the search efficiency. For future research, searching strategy design and techniques of parameter setting are both promising research directions for promoting the development of hierarchical optimization.

Conflicts of Interest:
The authors declare no conflict of interest. Table A1. Unimodal benchmark functions.