Deployment of On-Orbit Service Vehicles Using a Fuzzy Adaptive Particle Swarm Optimization Algorithm

,


Introduction
The on-orbit service of a spacecraft is a type of space operation that humans, robots, or both extend lives and improve their abilities in performing tasks of various spacecraft [1][2][3]. The main tasks of on-orbit service vehicles (OSVs) include on-orbit assembly, on-orbit maintenance, and logistic support [4]. As humans explore deeper in space, more and more diverse satellites are deployed in space. The main reasons why satellites are abandoned are that some components are broken and the fuel to maintain their orbits is not enough. The on-orbit service of a spacecraft is an effective method to reduce the waste of satellite resource and avoid producing space junk. Using OSVs to extend the lives of satellites deployed and improve the ability to perform tasks would be an important research direction.
The basic work of deployment of OSVs includes transfer orbit optimization and assignment optimization. In the existing literature, direct methods [5,6], indirect methods [7], and hybrid methods [8] are used to study the optimization of transfer trajectories of different types of spacecraft. OSV transfer orbit is one kind of transfer orbits. Although the direct methods suffer from the lack of high precision, its computation is efficient. Considering the deployment and assignment optimization need to calculate a lot of transfer orbits under different states, we use direct method to compute transfer orbits. The literature related to assignment optimization and the deployments of OSVs are as follows. A novel approach for the mission planning of on-orbit servicing such as visual inspection is proposed by Daneshjou et al. [9]. The on-orbit servicing task allocation is solved by a hybrid discrete particle swarm optimization algorithm which is proposed by Zhang and Zhang [10]. Pradeepmon et al. [11] solved quadratic assignment problems by using discrete particle swarm optimization (DPSO) algorithm.
Deployment of OSVs serving satellites, which is a complex multiple nested optimization problem, is almost not systematically studied in the existing literatures. Deployment of OSVs is the application of deployment problem. The optimization methods of other related deployment problems mainly include some intelligent algorithms. Tesarova & Vokalova [12] used genetic algorithm (GA) to deploy radio stations. Singh et al. [13] used particle swarm optimization (PSO) algorithm to deploy directional sensors. Ahmadian et al. [14] solved the large dimension electric distribution network problem by using a hybrid algorithm combining PSO and TS (Tabu search) algorithm (PSO+TS). Zhang et al. [15] obtained optimal deployment of flow monitors based on a hybrid algorithm combining GA and SA (simulated annealing) algorithm (GA+SA). The convergences of these intelligent algorithms are not good to compute complex multiple nested optimization problem. Thus, we propose a fuzzy adaptive particle swarm optimization (FAPSO) algorithm to improve the convergence of the results by balancing global and local search capability of PSO algorithm. The study background of the present study is that multiple satellites given fixed count and orbit elements need to be maintained by multiple OSVs in bounded time when an on-orbit serving mission order is set at any uncertain time t in a given time interval T (on-orbit mission). The main content of the present study consists of four parts (Sections 2-5). The framework of the present study is shown in Figure 1.
Part A (Section 2): on the basis of the double pulse rendezvous hypothesis and taking the maximum weight of remaining fuel as the optimization index (optimization index A), a transfer optimization model of single OSV serving multiple satellites is established based on the GA, which is used to calculate the indexes of the subsequent optimization models.
Part B (Section 3): at the time t of an on-orbit serving mission order set and taking the maximum value of both the count of satellites that can be served and the minimum value of a single OSV's optimization index A among all OSVs as the optimization index (optimization index B), an assignment optimization model of multiple OSVs serving multiple satellites is established based on DPSO algorithm.
Part C (Section 4): first of all, taking the time t in a given time interval T as the optimization variable and taking the minimum value of optimization index B as the optimization index (optimization index C), the optimization model to determine whether OSVs given initial count and orbit elements can achieve on-orbit mission or not is established based on the pattern search method and this is used to calculate the optimization index of the deployment model. Then, given the count of OSVs and taking the maximum value of optimization index C as the optimization index (optimization index D), the deployment optimization model of multiple OSVs serving multiple satellites is established based on the fuzzy adaptive particle swarm optimization (FAPSO) algorithm.
Part D (Section 5): the feasibility of all optimization models in the present study is verified through comparison with other intelligent algorithms, and the FAPSO algorithm is more efficient in solving the deployment problem of OSVs.

Transfer Orbit Optimization Model
A single OSV can serve multiple satellites. The procedure of the double pulse rendezvous method is that a single OSV manoeuvres because of the 1 st pulse thrust at some time which is called manoeuvre time, flies without thrust for some time, and manoeuvres because of the 2 nd pulse thrust to get the same velocity as the first satellite at some time which is called rendezvous time when both the OSV and the first satellite reach the same position. The OSV begins to serve the first satellite after the double pulse rendezvous procedure. After serving the first satellite, the OSV serves other satellites in the same way as it served the first satellite.
The procedure of the OSV serving the 1 st satellite by double pulse rendezvous method is shown in Figure 2. The orbit of the OSV is original orbit, and the orbit of a satellite is target orbit. The position of the OSV at the manoeuvre time is point P 1 , and the position of the OSV at the rendezvous time is point P 2 . When the OSV manoeuvres from point P 1 on the original orbit to point P 2 on target orbit in Δt seconds to serve the satellite, r 1 and v 1 are its position vector and velocity vector, respectively, at point P 1 , and r 2 and v 2 are its position vector and velocity vector, respectively, at point P 2 . The determination of the 1 st velocity increment Δv 11 (Formula (1)) of the OSV serving the 1 st satellite at point P 1 on the original orbit is called the Lambert problem, and the method to solve the Lambert problem is studied by Schumacher et al. [16]. The 2 nd velocity increment Δv 12 of the OSV serving the 1 st satellite at point P 2 on the target orbit is written as Formula (2). After serving the 1 st satellite for a certain time, the OSV begins to serve the 2 nd satellite in the same procedure as it served the 1 st satellite.
In Formula (1) and Formula (2), v 0 is the velocity of the OSV at point P 1 on the original orbit, v 1 is the velocity of the OSV at point P 1 on the transfer orbit after pushing by the 1 st pulse thrust, v 2 is the velocity of the OSV at point P 2 on the transfer orbit, and v 3 is the velocity of the OSV at point P 2 on the target orbit after pushing by the 2 nd pulse thrust.
Given the sequence of serving satellites and taking the maximum weight of the remaining fuel as the optimization index, the transfer orbit optimization model of a single OSV serving multiple satellites is established by the GA. The manoeuvre and rendezvous times are searched globally using the GA. Optimization index A can be obtained by calculating the velocity increments of the procedures of an OSV serving multiple satellites.  Modelling and Simulation in Engineering Given the sequence of serving satellites, the mathematical model can be written as max J = m remain s:t: In Formula (4), J is the optimization index A, m remain is the weight of the remaining fuel after the on-orbit mission is achieved, t 1 , t 3 , ⋯, t 2n−1 are the manoeuvre times, t 2 , t 4 , ⋯, t 2n are the rendezvous times, T max is the maximum time interval for finishing the on-orbit mission, t serve is the time interval for the OSV serving one satellite, n is the count of the transfer orbit segments (the i th segment corresponding to the i th satellite served by the OSV), Δm i is the weight of the consumed fuel in the i th segment, and m keep is the weight of fuel used to keep the orbit. m remain can be written as Formula (5) and Δm i can be written as Formula (6).

Modelling and Simulation in Engineering
In Formulas (5) and (6), k is the count of satellites served, K is the count of satellites that single OSV can serve, m i−1 is the weight of the remaining fuel after a single OSV serves i − 1 satellites, I sp is the specific impulse of the OSV, g o is the acceleration of the earth, m i−1 is the weight of the remaining fuel after the OSV serves the ði − 1Þ th satellite, Δv i1 is the 1 st velocity increment of the OSV serving the i th satellite at the manoeuvre time, and Δv i2 is the 2 nd velocity increment of the OSV serving the i th satellite at the rendezvous time.

Optimization Model.
The GA is an intelligence algorithm which is widely used in many fields [17][18][19]. Taking the maximum value of the remaining fuel weight (max m remain ) as the optimization index, the transfer optimization model of a single OSV serving multiple satellites is established based on the GA. The implementation steps are as follows.
Step 1. In order to meet the constraint conditions in Formula (4), the variables t 10 , t 21 , ⋯, t 2nð2n−1Þ are written as In Formula (7), n is the count of the transfer orbit segments. Furthermore, t 10 , t 21 , ⋯, t 2nð2n−1Þ are taken as genes encoded to a chromosome ( Figure 3). The initial populations are randomly generated. The population size M t , the crossover probability P c , the mutation probability P m , and the maximum number of iterations T A are set.
Step 2. The object function of Formula (4) is taken as the fitness function of the GA. Each chromosome has a fitness function value. The fitness function values of all chromosomes in the current generation are calculated.
Step 3. The chromosomes are selected by means of roulette wheel selection. Some genes on two different chromosomes reciprocally cross according to the crossover probability. Furthermore, some genes mutate according to the mutation probability.
Step 4. Steps 2 and 3 are repeatedly executed until the GA converged (jJðR i Þ − JðR i−1 Þj ≤ 10 −6 ) or the maximum number of iterations is reached.
Given the sequence of serving satellites, the maximum J could be obtained through the GA.

Assignment Optimization Model
In Formula (9), max f ðXÞ is the optimization index B. Furthermore, K is the count of satellites that a single OSV can serve, m is the weight of a single OSV, m s remain is the weight of the remaining fuel after the s th OSV serves the satellites, and m keep is the weight of the fuel used to keep the orbit.
Optimization index B includes two parts. The first part (optimization index B 1 ) is written as the formula ∑ N i=1 ∑ M j=1 x ij , which refers to the count of satellites served. This part evaluates the use ratio of OSVs.
The second part (optimization index B 2 ) is written as the formula min m s remain /m, which refers to the minimum weight (min m s remain ) of a single OSV's remaining fuel among all OSVs divided by the weight of the OSV (m). The purpose of the present study is to use the minimum count OSVs to serve the satellites given the count and their orbits based on the on-orbit mission. The formula m remain refers to the weight of the remaining fuel after a single OSV serves satellites, and this remaining fuel can be used to keep the orbit before the OSV execute its on-orbit mission, which evaluates the ability of a single OSV keeping its orbit. If the value of min m s remain is larger, the overall ability of all OSVs to keep orbits would be better. The maximum value of the formula min m s remain /m is m fuel /m, which the value is less than 1. In optimization index B, optimization B 1 is the main optimization index. The value of increasing or decreasing one satellite (the definition of optimization index B 1 ) is larger than the maximum value of optimization index B 2 .  Modelling and Simulation in Engineering 3.2. Optimization Model. The assignment optimization model of multiple OSVs serving multiple satellites is established based on the DPSO [20,21], and the main components of the DPSO algorithm are as follows.
3.2.1. Coding of a Particle. The particles are encoded in decimal system. The length of each particle is N × K. The coding of a particle is presented as Figure 4. In Figure 4, O1, O2, ⋯ , ON are N OSVs, and Si j refers to the j th satellite served by the i th OSV. The value of Si j is obtained from the set ð0, 1, 2 ,⋯,MÞ and 0 means that no satellite is served.

Calculation of the Fitness Function.
The object function of Formula (9) is taken as the fitness function shown as In Formula (10), Ps i is the i th particle.

Particle's Position and Velocity Update
Formulas. The position of particle swarm optimization (PSO) is determined through the particle velocity, individual extreme, and global extreme, and the particle's position and velocity update formulas of the PSO algorithm are written as [22,23] In Formula (11), P i h is the individual extreme in the h th iteration, P g h is the global extreme in the h th iteration, x i h is the fitness function value of the i th particle in the h th iteration, x i h+1 is the fitness function value of the i th particle in the ðh + 1Þ th iteration, v i h is the particle velocity value of the i th particle in the h th iteration, v i h+1 is the particle velocity value of the i th particle in the ðh + 1Þ th iteration, ω is the inertia weight that reflects the extent that the particle maintains the present velocity, c 1 is the cognitive coefficient that reflects the weight of the particle getting close to the individual extreme, c 2 is the social coefficient which reflects the weight the particle getting close to the global extreme, r 1 is a random number between 0 and 1, and r 2 is a random number between 0 and 1.
On the basis of the PSO, the particle position and velocity update formulas for the DPSO algorithm is written as follows: In Formula (12), Ps h i is the fitness function value of the i th particle in the h th iteration, and Ps h+1 i is the fitness function value of the i th particle in the ðh + 1Þ th iteration. This formula includes three formulas which are written as rand ðÞ ≥ c 2 : Ψ h i is the value of the i th particle in the h th iteration after the calculation of ω ⊗ F 1 ðPs h i Þ. The formula rand ðÞ means that a real number between 0 and 1 is randomly generated. When the condition rand ðÞ < ω is met, the formula F 1 ðPs h i Þ is executed which means that two integers a and b are randomly generated in interval ½0, N and the values on the i th particle at positions a and b are interchanged. When the condition rand ðÞ ≥ ω is met, the formula Ps h i is executed which means that the coding values on the i th particle are unchanged.
Φ h i is the value of the i th particle in the h th iteration after the calculation of c 1 ⊗ F 2 ðΨ h i , P h i Þ. P h i is the individual extreme in the t th iteration. When the condition rand ðÞ < c 1 is met, the formula F 2 ðΨ h i , P h i Þ is executed which means that two integers a and b are randomly generated in interval ½0, N, and the values on Ψ h i at positions a and b as well as the values on P h i at positions a and b are interchanged. When the condition rand ðÞ ≥ c 1 is met, the formula Ψ h i is executed which means that the coding values on Ψ h i are unchanged. P h g is the global extreme in the t th iteration. When the condition rand ðÞ < c 2 is met, the formula F 3 ðΦ h i , P h g Þ is executed which means that two integers a and b are randomly generated in interval ½0, N, and the values on Φ h i at positions a and b as well as the values on P h g at positions a and b are interchanged. When the condition rand ðÞ ≥ c 2 is met, the formula Φ h i is executed which means that the coding values on Φ h i are unchanged. The implementation steps of the DPSO algorithm are as follows.
Step 1. Set parameters: the population size M B of the particle swarm, the count N of OSVs, the maximum count K of satellites served by a single OSV, the maximum number of Figure 4: Coding of particle.

Modelling and Simulation in Engineering
iterations T B , the inertia weight ω, the cognitive coefficient c 1 , the social coefficient c 2 , and h = 1.
Step 2. The particle swarm population is initialized, the values on the particles are randomly generated, the fitness function f ðPs i Þ values of all particles are calculated, and the individual extreme P i h and global extreme P g h are recorded.
Step 3. h = h + 1, the position of the particles is updated through Formula (12), the fitness function f ðPs i Þ values of all particles are calculated, and the individual extreme P i h and global extreme P g h are updated.
Step 3 is executed repeatedly until (jP g h − P g h−1 j ≤ 10 −6 ) or the maximum number of iterations is reached.

Deployment Model
4.1. Problem Description. The deployment problem of OSVs includes two parts: one is to determine the count of OSVs, while the other is to determine the orbit elements of all OSVs.
A single OSV serves K satellites at most. If M satellites need to be served, the initial count NN of satellites can be calculated: In Formula (16), the function ceilðÞ means rounding up to the nearest integer. After the count of OSVs is determined by Formula (16), the orbit element orbit radius (A), orbit eccentricity ratio (e), orbit inclination angle (I), right ascension of the ascending node (Ω), argument of perigee (α), and true anomaly (f ) of NN OSVs are generated according to the orbit elements of given satellites.
The optimization model to determine whether OSVs' given count and orbit elements can achieve the on-orbit mission or not is established and its optimization index is called optimization index C, which is the foundation of the deployment optimization model. If the given OSVs could not achieve the on-orbit mission, the orbit elements of all OSVs are adjusted to determine whether the on-orbit mission can be achieved and obtain the best optimization index. If the on-orbit mission could still not be achieved, the count of OSVs is increased and the orbit elements of all OSVs are adjusted again to obtain the best optimization index.

Optimization Index.
Given the count and orbit elements of OSVs and taking the time interval T as the searching region for optimization variable t, the optimization model to determine whether OSVs' given count and orbit elements can achieve on-orbit mission or not is established based on the pattern search method and its optimization index is called optimization index C (min J 1 ). In the entire searching region, optimization index B of the assignment optimization model at any time point can be calculated and minimum optimization index B is called optimization index C. The optimization index of the deployment optimization model is called optimization index D (max J 2 ).
In Formula (17), max f ðXÞ is optimization index B. The implementation steps of the pattern search [24] are as follows.
Step 1. i = 1. The function value J 1 ðB i Þ at base point B i in the time interval T is calculated. The iteration step size δt and the convergence precision ε are set.
Step 2. The temporary point T i is determined through Formula (18)  Step 3. If the condition J 1 ðT i Þ < J 1 ðB i Þ is met, the ði + 1Þ th vector is twice longer than the i th vector. i = i + 1. The temporary point T i is written as The temporary point T i is defined as the base point B i+1 which is taken as the starting point of the i th vector. If the condition J 1 ðT i Þ ≥ J 1 ðB i Þ is met, the ði + 1Þ th vector is half of the i th vector and its starting point becomes the base point B i−1 .
Step 4. Steps 2 and 3 are repeatedly executed. If the value of T i remains unchanged, the iteration step size δt will be half until it reaches a certain degree of accuracy. The algorithm ends.

Optimization Model.
Given the count of OSVs, the optimized orbit elements are obtained by adjusting orbit elements of all OSVs. The deployment problem of OSVs is a high-dimension nonlinear problem which nests the optimization problems in Sections 2.2, 3.2, and 4.2.
If the inertia weight ω is bigger, the global search capability of the PSO algorithm will be stronger. If the inertia weight ω is smaller, the local search capability of the PSO algorithm will be stronger. Considering that the searching space of deployment problem of OSVs is large and the optimization variables are sensitive to the initial values, a larger inertial weight ω will be chosen to enhance the global search capability. Considering that the fuel taken by the OSV is expensive, enhancing the local search capability to obtain the better local results will also be very important. Hence, the fuzzy control system is introduced to adjust the inertia weight ω, which can balance the global search capability and local search capability. 6 Modelling and Simulation in Engineering The essence of the FAPSO algorithm is updating the inertia weight ω through the fuzzy controller. The updated inertial weight ω updated is determined through the fuzzy controller in three steps.
Step 1. The input variables of the fuzzy controller are translated into fuzzy variables. The two input variables are the inertial weight ω and the normalized fitness function values fit of particles in the present iteration.
In Formula (20), fit is the fitness function value of one particle in the present iteration, fit min is the minimum fitness function value in the present iteration, and fit max is the maximum fitness function value in the present iteration.
The fuzzy variables of these two input variables are calculated through the membership functions [25,26].
Formulas (21), (22), and (23) are the membership functions that correspond to fuzzy language "small," "medium," and "large," respectively. a S and b S are the constants when    7 Modelling and Simulation in Engineering two input variables ω or fit are within the range of fuzzy language "small." a M and b M are the constants when two input variables ω or fit are within the range of fuzzy language "medium." a L and b L are the constants when two input variables ω or fit are within the range of fuzzy language "large." Step 2. The values of two input variables ω and fit and their fuzzy variables are used to calculate the fuzzy control variable Δω according to the fuzzy reasoning algorithm and the fuzzy control rule base.
The fuzzy reasoning algorithm can be expressed through In Formula (24), x is the input variable ω, y is the input variable fit, X is the set where x is in, Y is the set where y is in, RðX, YÞ is the implication relationship between X and Y , μ A ðxÞ is the membership function value of x, μ B ðyÞ is the membership function value of y, and "∧" is the operational symbol for choosing the smaller value between μ A ðxÞ and μ B ðyÞ.
The fuzzy control rule base is described in Table 1.
Step 3. The fuzzy control variable Δω is translated into the precise value after the defuzzification procedure. The precise value v out of a particle is the centre point's x-coordinate value of the area of the curve f ðΔωÞ = μ RðX,YÞ ðx, yÞ in Step 2 and the x-coordinate.
In Formula (25), μ V ðvÞ is μ RðX,YÞ ðx, yÞ, V is the area of the curve f ðΔωÞ = μ RðX,YÞ ðx, yÞ in Step 2 and the x-coordinate, and v is the x-coordinate value. The updated inertial weight ω updated is written as The implementation steps of the FAPSO algorithm are as follows.
Step 1. Set parameters: the population size M C of the particle swarm, the count N of OSVs, the maximum number of iterations T C , the inertia weight ω, the cognitive coefficient c 1 , the social coefficient c 2 , and h = 1. The coding of the i th particle is the expression x i j6 refers to the k th orbit element of the j th OSV in the i th particle, and k = 1, 2, 3, 4, 5, 6 corresponds to the six orbit elements A, e, I, Ω, α, and f , respectively.
Step 2. The particle swarm population is initialized, the values on the particles are randomly generated, the fitness function values of all particles are calculated using Formula (17), and the individual extreme P i h and the global extreme P g h are recorded.
Step 3. The inertial weight ω and the normalized fitness function values fit in the h th iteration are calculated and the updated inertial weight ω updated is determined through the fuzzy controller in three steps. h = h + 1. The updated inertial weight ω updated is placed into Formula (11) to calculate the fitness function value of the particle in the h th iteration. The fitness function values of all particles are calculated, and the individual extreme P i h and global extreme P g h are recorded.
Step 3 is repeatedly executed until (jP g h − P g h−1 j ≤ 10 −6 ) or the maximum number of iterations is reached.
Step 5. If the optimization result cannot achieve the on-orbit mission, one OSV is added to the present OSVs and Steps 1-4 are executed. If the optimization result can achieve the onorbit mission, the cycle calculation is ended.  Table 3. The weight of a single OSV is 1000 kg. A single OSV serves three satellites at most and carries 600 kg of fuel. The specific  S2, S4 O4 S10, S17, S15 O5 S12, S21, S5 O6 S13, S16 O7 S11, S9, S3 O8 S19, S14 8
impulse is 3000 seconds. The duration of time to serve a satellite is more than an hour. The on-orbit mission should be achieved within 72 hours.
In the simulation of the transfer orbit optimization model, the parameters of GA are taken from literature written by Tesarova & Vokalova [12]. The population size M A is 20, the weight of the fuel used to keep the orbit m keep is 50 kg, the maximum number of iterations T A is 50, the crossover probability P c is 0.85, and the mutation probability P m is 0.1.
In the simulation of assignment optimization model, the parameters of DPSO are taken from literature written by Pradeepmon et al. [11]. The population size M B is 10, the maximum number of iterations T B is 50, the inertia weight ω is 0.75, the cognitive coefficient c 1 is 2, and the social coefficient c 2 is 2.
In the simulation of deployment model, the parameters of PSO are taken from literature written by Singh et al. [13]. The population size M C is 10, the maximum number of iterations T C is 50, the inertia weight ω is 0.75, the cognitive coefficient c 1 is 2, and the social coefficient c 2 is 2. The particle velocities corresponding to the six orbit elements A, e, I, Ω, α, and f are 50000, 0.00001, 10, 10, 10, and 10, respectively. In the optimization model to determine whether OSVs' given count and orbit elements can achieve on-orbit mission or not, the iteration step size δt of the pattern search is (24 × 60 × 60) seconds and the time interval T is    Figure 8: Transfer orbit of O3 serving S10, S17, and S15. 10 Modelling and Simulation in Engineering (2 × 365 × 24 × 60 × 60) seconds. The parameters of the fuzzy controller are shown in Table 4. We design the comparison experiments between the proposed FAPSO algorithm and GA, PSO, PSO+TS, and GA+SA algorithms. The parameters of these algorithms are successively taken from the existing literatures mentioned in the present paper [12][13][14][15]. 30 independent runs are performed for each algorithm.

Analyses of Results.
The results of the 10 th experiment presented in Tables 5-7 and Figures 5-13 are obtained using the FAPSO algorithm. The obtained orbit elements of OSVs are presented in Table 5, and the parking orbits of OSVs are presented in Figure 5. The deployment results of OSVs are presented in Table 6 including the value of optimization D (max J 2 ), the time to obtain max J 2 , and the assignment results of OSVs. On the basis of results presented in Table 6, OSVs' manoeuvre time, rendezvous time, and the weight of remaining fuel are presented in Table 7. The transfer orbits of OSVs are presented in Figures 6-13. The convergence curves (the average values of 30 independent experiments) of the FAPSO algorithm, GA, PSO algorithm, PSO+TS algorithm, and GA+SA algorithm are presented in Figure 14.
In Table 6, the time to obtain optimization index D (max J 2 ) is after 9336.450 hours. O1, O2, O4, O6, and O7 O4 Parking orbit S10 Parking orbit S17 Parking orbit S15 Parking orbit O4 -S10 Transfer orbit O4-S17 Transfer orbit O4-S15 Transfer orbit  Figure 9: Transfer orbit of O4 serving S12, S21, and S5.   Figure 10: Transfer orbit of O5 serving S11, S9, and S3. each serve 3 satellites. O3, O5, and O8 each serve 2 satellites. O1 serves S8, S1, and S18 in sequence. O2 serves S7, S20, and S6 in sequence. O3 serves S2 and S4 in sequence. O4 serves S10, S17, and S15 in sequence. O5 serves S12, S21, and S5 in sequence. O6 serves S13 and S16 in sequence. O7 serves S11, S9, and S3 in sequence. O8 serves S19 and S14 in sequence. In Table 7, t 1 and t 2 are the manoeuvre time and rendezvous time to serve the 1 st satellite, respectively, t 3 and t 4 are the manoeuvre time and rendezvous time to serve the 2 nd satellite, respectively, t 5 and t 6 are the manoeuvre time and rendezvous time to serve the 3 rd sat-ellite, respectively, and m s remain is the weight of remaining fuel after the s th OSV serving the satellites. In Figure 14, the convergence results of the FAPSO algorithm, PSO algorithm, PSO+TS algorithm, GA, and GA+SA algorithm are 21.106, 21.083, 21.084, 21.080, and 21.094, respectively. After the convergence results are entered into the objection function of Formula (9), the weight of the remaining fuel obtained through the FAPSO algorithm is 23 kg, 22 kg, 26 kg, and 12 kg more than that obtained by the PSO algorithm, PSO+TS algorithm, GA, and GA+SA algorithm, respectively. O6 Parking orbit S13 Parking orbit S16 Parking orbit O6-S13 Transfer orbit O6-S16 Transfer orbit Figure 11: Transfer orbits of O6 serving S13 and S16. Km O7-S9 Transfer orbit O7-S3 Transfer orbit O7-S11 Transfer orbit O7 Parking orbit S11 Parking orbit S9 Parking orbit S3 Parking orbit Figure 12: Transfer orbit of O7 serving S19 and S14.

Conclusion and Further Research
The deployment of multiple OSVs is a complex multiple nested optimization problem, which is systematically solved by a series of intelligent algorithms in the present study. The deployment model includes the transfer orbit optimization model and the assignment optimization model. The transfer orbit optimization is a continuous optimization problem, which is solved by the intelligence algorithm GA. The assignment optimization is a discrete optimization problem, which is solved by the intelligence algorithm the DPSO algorithm. The deployment is a continuous optimization problem, which is solved by the FAPSO algorithm. The FAPSO algorithm, which has a better convergence result than that obtained using the other similar optimization algorithms, can effectively solve the deployment problem of OSVs.
In further studies, some aspects should be considered such as the finite thrust engine of the OSV, the influence on the parking orbits because of the disturbing force, and optimization algorithms getting the simulation results faster. Km Km Km O8 Parking orbit S19 Parking orbit S14 Parking orbit O8-S19 Transfer orbit O8-S14 Transfer orbit