Fuel-Efficient on-Orbit Service Vehicle Allocation Based on an Improved Discrete Particle Swarm Optimization Algorithm

Given the limited fuel capacity of an on-orbit service vehicle (OSV), proper OSV allocation to satellites during each service mission is critical for economic fuel consumption. +is allocation problem can be formulated as an optimization problem with many continuous and discrete design variables of wide domains. +is problem can be effectively handled through the proposed approach that combines the tabu search with the discrete particle swarm optimization algorithm (DPSO-TS). First of all, Pontryagin’s minimum principle and genetic algorithm (GA) are exploited to find the most fuel-efficient transfer trajectory. +is fuel efficiency maximization can then serve as the performance index of the OSV allocation optimization model problem. In particular, the maximization of the minimum residual fuel over individual OSVs is proposed as a performance index for OSV allocation optimization. +e optimization problem is numerically solved through the proposed DPSO-TS algorithm. Finally, the simulation results demonstrate that the DPSO-TS algorithm has a higher accuracy compared to the DPSO, the DPSO-PDM and the DPSO-CSA algorithms in the premise that these four algorithms have the basically same computational time. +e DPSO-TS algorithm can effectively solve the OSV allocation optimization problem.


Introduction
As humans explore space deeper, more and more diverse satellites are deployed into space. On-orbit service of satellites is a type of space operations performed by astronauts, robots, or both to extend the lives of satellites and improve their abilities to perform various tasks [1]. Onorbit service vehicles (OSV) are particularly useful as an effective tool to extend lives of satellites. ese vehicles perform several tasks including on-orbit assembly, maintenance, and logistic support [2]. Many important problems of the OSV are investigated in the existing literature. e study of transfer trajectory is the basic work of allocation of OSVs. Direct methods [3][4][5], indirect methods [6,7], and hybrid methods [8,9] are used to study the optimization of transfer trajectories of different types of space vehicles. e hybrid methods, one of which is adopted in this paper to optimize the transfer trajectories, combine the advantages of the direct and indirect methods.
In general, multiple satellites need to be maintained in one service mission. e allocation of OSVs to satellites is an important problem that has the characteristics of a general resource allocation problem. However, the OSV allocation problem is distinguished by two aspects: (a) there are many optimization variables with a wide search scope for each, and (b) the optimization variables are of both the continuous and discrete types. Allocation optimization problem is a complex combinatorial optimization problem and known to be an NP hard, and discrete particle swarm optimization (DPSO) is a commonly used algorithm to solve allocation optimization problems. e traditional DPSO algorithm [10][11][12] is easy to be trapped in the local optimum, and hence some improved DPSO algorithms are studied. Inertia weight of DPSO is used to balance the global and local search capacity, and then the way of designing the adjustment strategy of inertia weight to solve the problem of falling into local optimum is widely studied. Li et al. [13] propose DPSO-PDM algorithm by the hybridizing adjustment strategy based on new particle diversity and adaptive mutation strategy, which effectively solves the problem of complex network community detection. Xin et al. [14] propose an improved DPSO algorithm by designing the linearly decreasing inertia weight, which effectively solves the problem of scheduling MIMO radar tasks. Nagra et al. [15] propose an improved a hybrid self-inertia weight adaptive PSO algorithm. e methods mentioned above improve the performance of DPSO algorithm by adjusting the inertia weight. ere are some methods combining DPSO with other algorithms to improve the performance of DPSO algorithm. Wang et al. [16] propose SA-DPSO algorithm which combines the DPSO algorithm and simulated annealing algorithm. Vairam et al. [17] propose CHIDPSO algorithm by combining DPSO algorithm and cyber swarm algorithm, which improves the local search capacity of DPSO algorithm. Vasudevan and Sinha [18] proposed a hybrid optimization algorithm by combining GA and PSO. Shayeghi et al. [19] improve the mutation of DPSO based on similarity algorithm for optimization of transmission lines loading. Qi [20] proposes an improved DPSO algorithm by combing the iterated local search. Yang and Xu [21] propose DPSO-MFO algorithm, which is able to provide excellent performances as a natureinspired heuristic algorithm. e methods mentioned above improve the performance of DPSO by enhancing the exploration and exploitation abilities, but the balance of exploration and exploitation is not controlled. Too much exploitation makes the algorithm converge quickly often to a local optimum, and too much exploration slows down the convergence of the procedure although it increases the probability of finding the global or near-optimum solution of the optimization problem. e tabu search (TS) algorithm is a hill-climbing algorithm, which is simple and easy to implement [22][23][24]. In order to keep the balance of exploration and exploitation abilities, we introduce TS algorithm into DPSO algorithm to improve the performance of DPSO algorithm. e proposed method keeps the balance of exploration and exploitation abilities by controlling the inertia weight and the length of the tabu table. e background of this paper is that multiple satellites need to be maintained in one service mission. As fuel is extremely important for OSVs, the fuelefficient problem becomes the focus of allocation optimization problem. is paper is organized as follows. First, we establish in Section 2 an optimization model for the minimum duration transfer trajectory by using a hybrid method. en, we propose in Section 3 an algorithm that combines DPSO algorithm with TS algorithm. is DPSO-TS algorithm is used for solving the OSV allocation optimization problem. Finally, in Section 4, we carry on several simulations to demonstrate the superiority of the DPSO-TS algorithm over DPSO algorithm and some, and the results demonstrated have higher convergence accuracy than DPSO algorithm and two relevant representative algorithms including DPSO-PDM algorithm and CHIDPSO algorithm.
In order to make the paper readable, the nomenclature of variables is shown in Table 1 before the optimization models are established.

Transfer Trajectory Optimization Based on a Hybrid Method
In this section, we provide the differential equations of motion of an OSV in an Earth-centred inertial (ECI) coordinate system. en, an optimization model of the OSV transfer trajectory is established.

OSV Equations of Motion.
e differential equations of motion [25] are where r is the position vector, v is the velocity vector, F is the thrust, m is the OSV mass, I sp is the specific impulse of the engine, and α � [α x , α y , α z ] T is the unit vector of in the thrust direction.

Transfer Trajectory Optimization.
Assuming that the maneuvering time of each OSV is known or has been already determined, an optimization model for the transfer trajectory between each pair of an OSV and a satellite is established based on a hybrid method, which exploits Pontryagin's minimum principle [26,27] and genetic algorithm (GA) [28]. e more fuel one OSV has retained, the longer it can stay in the orbit and the more satellites it can serve. erefore, the mass of the remaining fuel can be taken as the optimization index to be maximized, namely, where t f is the terminal time of each transfer trajectory and m is a function of the change in the OSV fuel mass. Based on the OSV equations of motion in formula (1), the Hamilton function can be derived as where λ r and λ v are the adjoint variables of the state variables and λ m is the adjoint variable of the OSV mass m. According to Pontryagin's minimum principle, the optimal thrust direction α * is obtained by sovling ‖α � 1‖ and zH/zα � 0. en, the optimal thrust direction is e partial derivatives of the Hamiltonian with respect to the adjoint variables are thus obtained. Consequently, the adjoint variables can be shown to satisfy the following differential equations: 2 Mathematical Problems in Engineering Initial boundary conditions are and the terminal boundary constraints are e optimal control under finite-thrust conditions is the switching mode of the engine strategy (bang-bang control). In this paper, the adopted switching mode is the on-off-on mode of the engine working sequence, and the switching moments are taken as the optimization variables. In the switching-on segment, the thrust is set to F max . In the switching-off segment, the thrust is set to zero.
For the GA-based hybrid method, the orbital optimization problem is converted into an optimization problem with parameters, including the initial values of the adjoint variables, the switching-off time t 1 off of the first thrust segment, the switching-on time t 2 on of the second thrust segment, and the switching-off time t 2 off of the second thrust segment. For maximizing the mass of the remaining OSV fuel under the given constraints, seventeen parameters β � [λ 1 (t 0 ), λ 2 (t 0 ), t 1 off , t 2 on , t 2 off ] are taken as optimization variables. e solution of this optimization problem using genetic algorithm is carried out as follows: Step 1: the design parameters β � [λ 1 (t 0 ), λ 2 (t 0 ), t 1 off , t 2 on , t 2 off ] are encoded. irty chromosomes are generated randomly to form the initial population. e population size P T , the crossover probability P c , the mutation probability P m , and the maximum number of iterations N t are set.
Step 2: for each chromosome, the equations are solved to get the optimal thrust direction under the conditions (6) and (7). en, the equations of motion (1) are solved to obtain the transfer trajectory. e fitness values of all chromosomes in the current generation are computed based on the fitness function in formula (2).
Step 3: for the next generation, chromosomes are selected using the roulette-wheel selection method. Some genes on two different chromosomes reciprocally cross according to the crossover probability and mutate according to the mutation probability. e execution of the selection, crossover, and mutation operations produces the next generation population.
Step 4: steps 2 and 3 are repeated until the GA fitness function values go below ε � 10 − 6 or the maximum number of iterations N is reached. e hybrid method results in the OSV transfer trajectory with the maximum mass of the remaining fuel at a given maneuvering moment.

OSV Allocation Optimization Based on the DPSO-TS Algorithm
e OSV allocation optimization problem can be solved through the DPSO algorithm. e inertia weight ω is an important parameter of the DPSO algorithm. e larger the inertia weight ω is, the stronger the global search capability is. TS algorithm is a hill-climbing algorithm with good local where i � 1, 2, . . . , N, j � 1, 2, . . . , N and t m is maneuvering time. X includes the allocation plan of OSVs and their maneuevring times. e mathematical model of OSV allocation problem can be written as where m k remain is the remaining fuel mass of the k th OSV after maintaining, and m keep is the mass of the fuel used to keep the orbit. e inner min operator min[m 1 remain (X), m 1 remain (X), . . . , m s remain (X)] returns the minimum mass of the remaining fuel among individual OSVs used to maintain satellites according to X. e outer max operator max min[m 1 remain (X), m 1 remain (X), . . . , m s remain (X)] searches for the allocation plan X that maximizes the minimum mass value. e optimization index reflects the overall service capability of OSVs.

Optimization Problem Solution Based on the DPSO-TS
Algorithm.
e OSV allocation optimization problem can be solved using the proposed DPSO-TS algorithm. Details of the DPSO-TS algorithm are given for two parts: the implementation steps of DPSO-TS algorithm and DPSO algorithm and some steps and strategies of DPSO-TS algorithm.

3.2.1.
e Implementation Steps of DPSO-TS Algorithm. e implementation steps of DPSO-TS algorithm are Step 1: set the DPSO-TS parameters including the population size M B of the particle swarm, the OSV count N, the number of satellites K need to be maintained, the maximum number of iterations T D , the inertia weight ω, the cognitive coefficient c 1 , the social coefficient c 2 , and the iteration counter h � 1. Select a certain time step δ to discretize the mission duration t max .
Step 2: initialize the OSV particle swarm population. Set parameter values for each particle including an OSV serial number and a randomly generated maneuvering moment. For computational efficiency, maneuvering moments and OSV numbers of OSVs are separately used for the DPSO update by formula (12). e fitness function values f(Ps i ) of all particles are calculated, and the individual fitness extreme P i h and the global fitness extreme P g h are recorded.
Step 3: the particle positions are updated by formula (12). e fitness function values f(Ps i ) of all particles are calculated. e individual fitness extreme P h i and the global fitness extreme P h g are updated.
Step 4: the parameters of TS algorithm are set including the length l T of the tabu table, the size M T of the candidate set, and the maximum number of iterations T T .
Step 5: set the initial solution P T � P h g and the optimal solution P Tbest � P h g of the TS algorithm.
Step 6: according to the centrality and diversity search strategies, generate the neighborhood of the current solution, and select N T particles from the neighborhood as the candidate set Openlist(Y).
Step 7: according to the fitness function value f(Y i ), judge whether the amnesty criterion is met. If so, the corresponding candidate solution Y i will be placed in the tabu table Tabulist(Y) and the solution in the tabu table will be released. Otherwise, keep the tabu table intact.
Step 8: check whether the TS terminal conditions are satisfied. If so, output P h g � P Tbest . Otherwise, go to Step 6.
Step 9: increase the iteration number, h � h + 1. Check whether the DPSO termination conditions are satisfied. If so, output P h g . Otherwise, go to Step 3. e flow chart of DPSO-TS algorithm is shown in Figure 1.

Some Steps and Strategies of DPSO-TS Algorithm
(1) Coding of Particles. Particles of the OSV allocation problem are encoded in a decimal system where each particle has a length of 2 × N. Each particle is encoded as shown in Figure 2 Terminal conditions?
where Ps i is the i th particle.
(3) Particle Position and Velocity Update. e position of a particle is determined by the particle velocity, the individual fitness extreme, and the global fitness extreme [29]. e particle position and velocity in a PSO algorithm are updated as v h+1 where P h i is the individual fitness extreme of the i th particle in the h th iteration; P h g is the global fitness extreme in the h th iteration, x h i is the fitness function value of the i th particle in the h th iteration; x h+1 i is the fitness function value of the i th particle in the (h + 1) th iteration; v h i is the particle velocity value of the i th particle in the h th iteration; v h+1 i is the particle velocity value of the i th particle in the (h + 1) th iteration; ω is the inertia weight which reflects the extent to which a particle keeps the current velocity; c 1 is the cognitive coefficient, which measures how far a particle gets close to the individual fitness extreme; c 2 is the social coefficient, which measures how far a particle gets close to the global fitness extreme; and both r 1 and r 2 are random numbers in the interval [0, 1].
On the basis of the PSO, we propose the following DPSO update rules for the particle position and velocity: where Ps h i and Ps h+1 i are the fitness function values of the i th particle in the h th and (h + 1) th iterations, respectively. Formula (12) includes three formulas which are written as formulas (13)- (15).
In formula (13), Ψ h i is the value of the i th particle in the h th iteration and is obtained via the product ω ⊗ F 1 (Ps h i ). e function rand generates a random number from a uniform distribution on [0, 1]. If rand < ω, the function is computed leading to the random generation of two integers a and b, from the set (0, 1, 2, . . . , M), and hence the interchange of the values of the i th particle at positions a and b. If rand ≥ ω, the output is Ps h i which means that the coding values of the i th particle are kept unchanged.
In formula (14), Φ h i is the value of the i th particle in the h th iteration, and is obtained via the product In formula (15), P h g is the global fitness extreme position in the h th iteration. If rand < c 2 , function F 3 (Φ h i , P h g ) is computed leading to the random generation of two integers a and b, from the set (0, 1, 2, . . . , M), and hence the interchange of the values of Φ h i and P h g at positions a and b. If rand ≥ c 2 , the output is Φ h i which means that the coding values on Φ h i are kept unchanged.
(4) Some Strategies of TS Algorithm. In this paper, the initial solution of the TS algorithm is the global extremum. P g in each iteration of the DPSO algorithm. e fitness function is written as in formula (10). TS algorithm includes some strategies: if the solution is superior to the current optimal solution, the optimal solution should be amnesty; tabu tables which are first-in first-out (FIFO) queue are used to prevent repeated search for visited solutions and to avoid falling into local loops; the centralized search strategy refers to strengthening search in the neighborhood of good solutions. A diversity search strategy looks for solutions in unknown or highly unexplored regions.

Simulations.
e parking orbit elements of satellites (at 0:00 on July 1, 2099) in need of service are shown in Table 2.
e parking orbit elements of the OSVs (at 0:00 on July 1, 2099) are shown in Table 3. e maintaining time for a satellite is an hour ,and the on-orbit service mission should be achieved in 24 hours (From 0:00 on July 1, 2099 to 0:00 on July 2, 2099), which means t m + t tra (t tra is the duration of OSV and satellite from maneuvering time to rendezvous time) is in [0, 23 × 3600 s]. e mass of each OSV is 1000 kg, each OSV carries 600 kg fuel the engine thrust is 490 N, and the gas jet velocity of 3000 m/s. m keep is 100 kg. e population size P T is 20, the crossover probability P c is 0.7, the mutation probability P m is 0.1, and the maximum number of iterations N t is 30. e population size M B is 20, maximum number of iterations T D is 30, the inertia weight ω is 0.8, the cognitive coefficient c 1 is 0.6, and the social coefficient c 2 is 0.6. e length l T of the tabu table is 2, the size M T of candidate is 5, and the maximum number of iterations T T is 8. Time step δ is 1 s, and I � t m /δ. e simulation performs ten independent experiments. In order to reduce computational time, each trajectory is computed by one independent blade of blade servers.
We design comparison experiments with DPSO algorithm, DPSO-PDM algorithm, and CHIDPSO algorithm. e parameters of DPSO algorithm are the same as those of DPSO-TS algorithm. e parameters of DPSO algorithm, part of DPSO-PDM algorithm, are the same as those of DPSO-TS algorithm, and the other parameters are the same as those in the article written by Li et al. [13]. e parameters of the DPSO algorithm part of CHIDPSO algorithm are the same as those of DPSO-TS algorithm, and the other parameters are the same as those in the article written by Vairam et al. [17].  Tables 4-7, and the best results can be obtained at these maneuvering times. It can be seen from Tables 4-7 that the average mass of the remaining fuel obtained by DPSO-TS algorithm, DPSO algorithm, DPSO-PDM algorithm, and CHIDPSO algorithm is 538.5 kg, 467.8 kg, 488.9 kg, and 508.8 kg; the average mass of the remaining fuel obtained by DPSO-TS algorithm is better than that of the other algorithms. Except the ninth experiment of DPSO-TS algorithm, the mass of the remaining fuel obtained by DPSO-TS algorithm in the other experiments is superior to that obtained by the other algorithms in any experiment; the average computational time obtained by DPSO-TS algorithm, DPSO algorithm, DPSO-PDM algorithm, and CHIDPSO algorithm is 1792.3 s, 1623.2 s, 1769.1 s, and 1797.3 s. e average computational time obtained by DPSO-TS algorithm differs from the shortest average computational time by 169.1 s, which has no effect on the service mission because there is a long preparation time before the mission is executed. e computational time of these four algorithms can be considered basically the same.    Figures 3-7 in turn. Figures 8-12 show the comparison of best fitness value curves obtained by DPSO-TS algorithm, DPSO algorithm, DPSO-PDM algorithm, and CHIDPSO algorithm in the experiments with top 5 optimization results in turn, from which we can see that the DPSO-TS algorithm's fitness value of each iteration is superior to the other algorithms after the eleventh iteration. Figures 8-12 show that the TS algorithm has strong local search capacity.

Analyses of
e simulation results demonstrate that the DPSO-TS algorithm has a higher accuracy compared to the DPSO, the DPSO-PDM, and the DPSO-CSA algorithms and can effectively solve the OSV allocation optimization problem.

Conclusions and Future Research
Fuel-efficient OSV allocation is a basic work to accomplish the on-orbit service mission in the background of deploying a large number of satellites in space in the future. We establish a transfer trajectory optimization model to obtain the mass of remaining fuel, which lays a foundation of allocation model. en, we propose DPSO-TS algorithm to solve the allocation model, and the simulation results demonstrate that the proposed algorithm is more effective than DPSO algorithm, DPSO-PDM algorithm, and CHIDPSO algorithm. e DPSO-TS algorithm can be applied to other allocation optimization problems. Although the proposed method has high accuracy, it still has some limitations. Further accelerating the computational time and investigating more algorithms with strong local search capacity in future research would be better.

Data Availability
e simulation data used to support the findings of this study are included within the article.

Conflicts of Interest
e authors declare that they have no conflicts of interest.