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

Authors: Yaxiong Li; Xinglong Sun; Xinxue Liu; Jian Wu; Qingguo Liu
Journal: Modelling and Simulation in Engineering
 (2021)

---

## Abstract

On the basis that satellites given fixed count and orbit elements can be served in bounded time when an on-orbit serving mission order is set at any uncertain time in a given time interval, the deployment of on-orbit service vehicle (OSV) serving satellites becomes a complex multiple nested optimization problem, and the essence of deployment is to determine the count and orbit elements of OSVs. In consideration of the characteristics of this deployment problem, we propose a fuzzy adaptive particle swarm optimization (FAPSO) algorithm to solve this problem. First, on the basis of double pulse rendezvous hypothesis, a transfer optimization model of a single OSV serving multiple satellites is established based on genetic algorithm (GA), and this is used to compute the indexes of the subsequent two optimization models. Second, an assignment optimization model of OSVs is established based on the discrete particle swarm optimization (DPSO) algorithm, laying the foundation of the next optimization model. Finally, the FAPSO algorithm, which improves the performance of PSO algorithm by adjusting the inertia weight, is proposed to solve the deployment problem of multiple OSVs. The simulation results demonstrate that all optimization models in this study are feasible, and the FAPSO algorithm, which has a better convergence result than that obtained using the other optimization algorithms, can effectively solve the deployment problem of OSVs.


---

## Body

## 1. 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–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.Figure 1
Framework of the present study.Part A (Section2): 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 (Section3): 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 (Section4): 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 (Section5): 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.


## 2. 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 1st pulse thrust at some time which is called manoeuvre time, flies without thrust for some time, and manoeuvres because of the 2nd 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 1st 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 P1, and the position of the OSV at the rendezvous time is point P2. When the OSV manoeuvres from point P1 on the original orbit to point P2 on target orbit in Δt seconds to serve the satellite, r1 and v1 are its position vector and velocity vector, respectively, at point P1, and r2 and v2 are its position vector and velocity vector, respectively, at point P2. The determination of the 1st velocity increment Δv11 (Formula (1)) of the OSV serving the 1st satellite at point P1 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 2nd velocity increment Δv12 of the OSV serving the 1st satellite at point P2 on the target orbit is written as Formula (2). After serving the 1st satellite for a certain time, the OSV begins to serve the 2nd satellite in the same procedure as it served the 1st satellite.
(1)Δv11=v1−v0,(2)Δv12=v3−v2.Figure 2
The procedure of the double pulse rendezvous method.In Formula (1) and Formula (2), v0 is the velocity of the OSV at point P1 on the original orbit, v1 is the velocity of the OSV at point P1 on the transfer orbit after pushing by the 1st pulse thrust, v2 is the velocity of the OSV at point P2 on the transfer orbit, and v3 is the velocity of the OSV at point P2 on the target orbit after pushing by the 2nd 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.


### 2.1. Mathematical Model

From the viewpoint of energy, problems correlated to multiple OSVs serving multiple satellites are investigated in the present study. The weight of a single OSV includes the weight of the fuel it carries and the weight of its other parts.(3)m=mfuel+melse.In Formula (3), m is the weight of a single OSV, mfuel is the weight of the fuel carried by the OSV, and melse is the weight of the other parts.Given the sequence of serving satellites, the mathematical model can be written as(4)maxJ=mremains.t.t1<t2<t3<t4<⋯<t2n−1<t2n<Tmax−tservet3−t2≥tserve,⋯,t2n−1−t2n−2≥tserve∑i=1nΔmi≤mfuel−mkeep.In Formula (4), J is the optimization index A, mremain is the weight of the remaining fuel after the on-orbit mission is achieved, t1,t3,⋯,t2n−1 are the manoeuvre times, t2,t4,⋯,t2n are the rendezvous times, Tmax is the maximum time interval for finishing the on-orbit mission, tserve is the time interval for the OSV serving one satellite, n is the count of the transfer orbit segments (the ith segment corresponding to the ith satellite served by the OSV), Δmi is the weight of the consumed fuel in the ith segment, and mkeep is the weight of fuel used to keep the orbit. mremain can be written as Formula (5) and Δmi can be written as Formula (6).
(5)mremain=mfuel−∑i=1kΔmi,k≤K,(6)Δmi=mi−11−exp−Δvi1goIsp1−exp−Δvi2goIsp,i≤K.In Formulas (5) and (6), k is the count of satellites served, K is the count of satellites that single OSV can serve, mi−1 is the weight of the remaining fuel after a single OSV serves i−1 satellites, Isp is the specific impulse of the OSV, go is the acceleration of the earth, mi−1 is the weight of the remaining fuel after the OSV serves the i−1th satellite, Δvi1 is the 1st velocity increment of the OSV serving the ith satellite at the manoeuvre time, and Δvi2 is the 2nd velocity increment of the OSV serving the ith satellite at the rendezvous time.


### 2.2. Optimization Model

The GA is an intelligence algorithm which is widely used in many fields [17–19]. Taking the maximum value of the remaining fuel weight (max mremain) 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 t10,t21,⋯,t2n2n−1 are written as
(7)t2i2i−1=t2i−t2i−1≥3600,i=1,2,⋯,n,t2i−12i−2=t2i−1−t2i−2≥0,i=1,2,⋯,n.
In Formula (7), n is the count of the transfer orbit segments. Furthermore, t10,t21,⋯,t2n2n−1 are taken as genes encoded to a chromosome (Figure 3). The initial populations are randomly generated. The population size Mt, the crossover probability Pc, the mutation probability Pm, and the maximum number of iterations TA 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 (JRi−JRi−1≤10−6) or the maximum number of iterations is reached.Figure 3
Coding of the chromosome.Given the sequence of serving satellites, the maximumJ could be obtained through the GA.


## 2.1. Mathematical Model

From the viewpoint of energy, problems correlated to multiple OSVs serving multiple satellites are investigated in the present study. The weight of a single OSV includes the weight of the fuel it carries and the weight of its other parts.(3)m=mfuel+melse.In Formula (3), m is the weight of a single OSV, mfuel is the weight of the fuel carried by the OSV, and melse is the weight of the other parts.Given the sequence of serving satellites, the mathematical model can be written as(4)maxJ=mremains.t.t1<t2<t3<t4<⋯<t2n−1<t2n<Tmax−tservet3−t2≥tserve,⋯,t2n−1−t2n−2≥tserve∑i=1nΔmi≤mfuel−mkeep.In Formula (4), J is the optimization index A, mremain is the weight of the remaining fuel after the on-orbit mission is achieved, t1,t3,⋯,t2n−1 are the manoeuvre times, t2,t4,⋯,t2n are the rendezvous times, Tmax is the maximum time interval for finishing the on-orbit mission, tserve is the time interval for the OSV serving one satellite, n is the count of the transfer orbit segments (the ith segment corresponding to the ith satellite served by the OSV), Δmi is the weight of the consumed fuel in the ith segment, and mkeep is the weight of fuel used to keep the orbit. mremain can be written as Formula (5) and Δmi can be written as Formula (6).
(5)mremain=mfuel−∑i=1kΔmi,k≤K,(6)Δmi=mi−11−exp−Δvi1goIsp1−exp−Δvi2goIsp,i≤K.In Formulas (5) and (6), k is the count of satellites served, K is the count of satellites that single OSV can serve, mi−1 is the weight of the remaining fuel after a single OSV serves i−1 satellites, Isp is the specific impulse of the OSV, go is the acceleration of the earth, mi−1 is the weight of the remaining fuel after the OSV serves the i−1th satellite, Δvi1 is the 1st velocity increment of the OSV serving the ith satellite at the manoeuvre time, and Δvi2 is the 2nd velocity increment of the OSV serving the ith satellite at the rendezvous time.


## 2.2. Optimization Model

The GA is an intelligence algorithm which is widely used in many fields [17–19]. Taking the maximum value of the remaining fuel weight (max mremain) 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 t10,t21,⋯,t2n2n−1 are written as
(7)t2i2i−1=t2i−t2i−1≥3600,i=1,2,⋯,n,t2i−12i−2=t2i−1−t2i−2≥0,i=1,2,⋯,n.
In Formula (7), n is the count of the transfer orbit segments. Furthermore, t10,t21,⋯,t2n2n−1 are taken as genes encoded to a chromosome (Figure 3). The initial populations are randomly generated. The population size Mt, the crossover probability Pc, the mutation probability Pm, and the maximum number of iterations TA 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 (JRi−JRi−1≤10−6) or the maximum number of iterations is reached.Figure 3
Coding of the chromosome.Given the sequence of serving satellites, the maximumJ could be obtained through the GA.


## 3. Assignment Optimization Model

### 3.1. Mathematical Model

Assumption 1.
There areN OSVs that are remarked as O=O1,O2,⋯,ON, there are M satellites that are remarked as S=S1,S2,⋯,SM, a single OSV serves K satellites at most, and a single satellite can be served once at most. The assignment matrix X is written as
(8)xij=1,OiservingSj,0,OinotservingSj.In Formula (8), i=1,2,⋯,N, j=1,2,⋯,M. The mathematical model is written as
(9)maxfX=∑i=1N∑j=1Mxij+minmremainsms.t.∑j=1Mxij≤K,∀i=1,2,⋯,N∑i=1Nxij≤1,∀j=1,2,⋯,Mmremains≥mkeep,∀s=1,2,⋯,N.In Formula (9), maxfX 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, mremains is the weight of the remaining fuel after the sth OSV serves the satellites, and mkeep is the weight of the fuel used to keep the orbit.Optimization index B includes two parts.The first part (optimization index B1) is written as the formula ∑i=1N∑j=1Mxij, which refers to the count of satellites served. This part evaluates the use ratio of OSVs.The second part (optimization index B2) is written as the formula minmremains/m, which refers to the minimum weight (minmremains) 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 mremain 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 minmremains is larger, the overall ability of all OSVs to keep orbits would be better. The maximum value of the formula minmremains/m is mfuel/m, which the value is less than 1. In optimization index B, optimization B1 is the main optimization index. The value of increasing or decreasing one satellite (the definition of optimization index B1) is larger than the maximum value of optimization index B2.


### 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 isN×K. The coding of a particle is presented as Figure 4. In Figure 4, O1,O2,⋯,ON are N OSVs, and Sij refers to the jth satellite served by the ith OSV. The value of Sij is obtained from the set 0,1,2,⋯,M and 0 means that no satellite is served.Figure 4
Coding of particle.


#### 3.2.2. Calculation of the Fitness Function

The object function of Formula (9) is taken as the fitness function shown as
(10)maxfPsi=∑i=1N∑j=1Mxij+minmremainsm.In Formula (10), Psi is the ith particle.


#### 3.2.3. 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]
(11)vih+1=ω⋅vih+c1⋅r1⋅Pih−xih+c2⋅r2⋅Pgh−xih,xih+1=xih+vih+1.In Formula (11), Pih is the individual extreme in the hth iteration, Pgh is the global extreme in the hth iteration, xih is the fitness function value of the ith particle in the hth iteration, xih+1 is the fitness function value of the ith particle in the h+1th iteration, vih is the particle velocity value of the ith particle in the hth iteration, vih+1 is the particle velocity value of the ith particle in the h+1th iteration, ω is the inertia weight that reflects the extent that the particle maintains the present velocity, c1 is the cognitive coefficient that reflects the weight of the particle getting close to the individual extreme, c2 is the social coefficient which reflects the weight the particle getting close to the global extreme, r1 is a random number between 0 and 1, and r2 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:(12)Psih+1=c2⊗F3c1⊗F2ω⊗F1Psih,Pih,Pgh.In Formula (12), Psih is the fitness function value of the ith particle in the hth iteration, and Psih+1 is the fitness function value of the ith particle in the h+1th iteration. This formula includes three formulas which are written as
(13)Ψih=ω⊗F1Psih=F1Psih,Psih,rand<ω,rand≥ω,(14)Φih=c1⊗F2Ψih,Pih=F2Ψih,Pih,Ψih,rand<c1,rand≥c1,(15)Psih+1=c2⊗F3Φih,Pgh=F3Φih,Pgh,Φih,rand<c2,rand≥c2.Ψih is the value of the ith particle in the hth iteration after the calculation of ω⊗F1Psih. The formula rand means that a real number between 0 and 1 is randomly generated. When the condition rand<ω is met, the formula F1Psih is executed which means that two integers a and b are randomly generated in interval 0,N and the values on the ith particle at positions a and b are interchanged. When the condition rand≥ω is met, the formula Psih is executed which means that the coding values on the ith particle are unchanged.Φih is the value of the ith particle in the hth iteration after the calculation of c1⊗F2Ψih,Pih. Pih is the individual extreme in the tth iteration. When the condition rand<c1 is met, the formula F2Ψih,Pih is executed which means that two integers a and b are randomly generated in interval 0,N, and the values on Ψih at positions a and b as well as the values on Pih at positions a and b are interchanged. When the condition rand≥c1 is met, the formula Ψih is executed which means that the coding values on Ψih are unchanged.Pgh is the global extreme in the tth iteration. When the condition rand<c2 is met, the formula F3Φih,Pgh is executed which means that two integers a and b are randomly generated in interval 0,N, and the values on Φih at positions a and b as well as the values on Pgh at positions a and b are interchanged. When the condition rand≥c2 is met, the formula Φih is executed which means that the coding values on Φih are unchanged.The implementation steps of the DPSO algorithm are as follows.Step 1.
Set parameters: the population sizeMB of the particle swarm, the count N of OSVs, the maximum count K of satellites served by a single OSV, the maximum number of iterations TB, the inertia weight ω, the cognitive coefficient c1, the social coefficient c2, and h=1.Step 2.
The particle swarm population is initialized, the values on the particles are randomly generated, the fitness functionfPsi values of all particles are calculated, and the individual extreme Pih and global extreme Pgh are recorded.Step 3.
h=h+1, the position of the particles is updated through Formula (12), the fitness function fPsi values of all particles are calculated, and the individual extreme Pih and global extreme Pgh are updated.Step 4.
Step 3 is executed repeatedly until (Pgh−Pgh−1≤10−6) or the maximum number of iterations is reached.


## 3.1. Mathematical Model

Assumption 1.
There areN OSVs that are remarked as O=O1,O2,⋯,ON, there are M satellites that are remarked as S=S1,S2,⋯,SM, a single OSV serves K satellites at most, and a single satellite can be served once at most. The assignment matrix X is written as
(8)xij=1,OiservingSj,0,OinotservingSj.In Formula (8), i=1,2,⋯,N, j=1,2,⋯,M. The mathematical model is written as
(9)maxfX=∑i=1N∑j=1Mxij+minmremainsms.t.∑j=1Mxij≤K,∀i=1,2,⋯,N∑i=1Nxij≤1,∀j=1,2,⋯,Mmremains≥mkeep,∀s=1,2,⋯,N.In Formula (9), maxfX 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, mremains is the weight of the remaining fuel after the sth OSV serves the satellites, and mkeep is the weight of the fuel used to keep the orbit.Optimization index B includes two parts.The first part (optimization index B1) is written as the formula ∑i=1N∑j=1Mxij, which refers to the count of satellites served. This part evaluates the use ratio of OSVs.The second part (optimization index B2) is written as the formula minmremains/m, which refers to the minimum weight (minmremains) 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 mremain 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 minmremains is larger, the overall ability of all OSVs to keep orbits would be better. The maximum value of the formula minmremains/m is mfuel/m, which the value is less than 1. In optimization index B, optimization B1 is the main optimization index. The value of increasing or decreasing one satellite (the definition of optimization index B1) is larger than the maximum value of optimization index B2.


## 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 isN×K. The coding of a particle is presented as Figure 4. In Figure 4, O1,O2,⋯,ON are N OSVs, and Sij refers to the jth satellite served by the ith OSV. The value of Sij is obtained from the set 0,1,2,⋯,M and 0 means that no satellite is served.Figure 4
Coding of particle.


### 3.2.2. Calculation of the Fitness Function

The object function of Formula (9) is taken as the fitness function shown as
(10)maxfPsi=∑i=1N∑j=1Mxij+minmremainsm.In Formula (10), Psi is the ith particle.


### 3.2.3. 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]
(11)vih+1=ω⋅vih+c1⋅r1⋅Pih−xih+c2⋅r2⋅Pgh−xih,xih+1=xih+vih+1.In Formula (11), Pih is the individual extreme in the hth iteration, Pgh is the global extreme in the hth iteration, xih is the fitness function value of the ith particle in the hth iteration, xih+1 is the fitness function value of the ith particle in the h+1th iteration, vih is the particle velocity value of the ith particle in the hth iteration, vih+1 is the particle velocity value of the ith particle in the h+1th iteration, ω is the inertia weight that reflects the extent that the particle maintains the present velocity, c1 is the cognitive coefficient that reflects the weight of the particle getting close to the individual extreme, c2 is the social coefficient which reflects the weight the particle getting close to the global extreme, r1 is a random number between 0 and 1, and r2 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:(12)Psih+1=c2⊗F3c1⊗F2ω⊗F1Psih,Pih,Pgh.In Formula (12), Psih is the fitness function value of the ith particle in the hth iteration, and Psih+1 is the fitness function value of the ith particle in the h+1th iteration. This formula includes three formulas which are written as
(13)Ψih=ω⊗F1Psih=F1Psih,Psih,rand<ω,rand≥ω,(14)Φih=c1⊗F2Ψih,Pih=F2Ψih,Pih,Ψih,rand<c1,rand≥c1,(15)Psih+1=c2⊗F3Φih,Pgh=F3Φih,Pgh,Φih,rand<c2,rand≥c2.Ψih is the value of the ith particle in the hth iteration after the calculation of ω⊗F1Psih. The formula rand means that a real number between 0 and 1 is randomly generated. When the condition rand<ω is met, the formula F1Psih is executed which means that two integers a and b are randomly generated in interval 0,N and the values on the ith particle at positions a and b are interchanged. When the condition rand≥ω is met, the formula Psih is executed which means that the coding values on the ith particle are unchanged.Φih is the value of the ith particle in the hth iteration after the calculation of c1⊗F2Ψih,Pih. Pih is the individual extreme in the tth iteration. When the condition rand<c1 is met, the formula F2Ψih,Pih is executed which means that two integers a and b are randomly generated in interval 0,N, and the values on Ψih at positions a and b as well as the values on Pih at positions a and b are interchanged. When the condition rand≥c1 is met, the formula Ψih is executed which means that the coding values on Ψih are unchanged.Pgh is the global extreme in the tth iteration. When the condition rand<c2 is met, the formula F3Φih,Pgh is executed which means that two integers a and b are randomly generated in interval 0,N, and the values on Φih at positions a and b as well as the values on Pgh at positions a and b are interchanged. When the condition rand≥c2 is met, the formula Φih is executed which means that the coding values on Φih are unchanged.The implementation steps of the DPSO algorithm are as follows.Step 1.
Set parameters: the population sizeMB of the particle swarm, the count N of OSVs, the maximum count K of satellites served by a single OSV, the maximum number of iterations TB, the inertia weight ω, the cognitive coefficient c1, the social coefficient c2, and h=1.Step 2.
The particle swarm population is initialized, the values on the particles are randomly generated, the fitness functionfPsi values of all particles are calculated, and the individual extreme Pih and global extreme Pgh are recorded.Step 3.
h=h+1, the position of the particles is updated through Formula (12), the fitness function fPsi values of all particles are calculated, and the individual extreme Pih and global extreme Pgh are updated.Step 4.
Step 3 is executed repeatedly until (Pgh−Pgh−1≤10−6) or the maximum number of iterations is reached.


## 3.2.1. Coding of a Particle

The particles are encoded in decimal system. The length of each particle isN×K. The coding of a particle is presented as Figure 4. In Figure 4, O1,O2,⋯,ON are N OSVs, and Sij refers to the jth satellite served by the ith OSV. The value of Sij is obtained from the set 0,1,2,⋯,M and 0 means that no satellite is served.Figure 4
Coding of particle.


## 3.2.2. Calculation of the Fitness Function

The object function of Formula (9) is taken as the fitness function shown as
(10)maxfPsi=∑i=1N∑j=1Mxij+minmremainsm.In Formula (10), Psi is the ith particle.


## 3.2.3. 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]
(11)vih+1=ω⋅vih+c1⋅r1⋅Pih−xih+c2⋅r2⋅Pgh−xih,xih+1=xih+vih+1.In Formula (11), Pih is the individual extreme in the hth iteration, Pgh is the global extreme in the hth iteration, xih is the fitness function value of the ith particle in the hth iteration, xih+1 is the fitness function value of the ith particle in the h+1th iteration, vih is the particle velocity value of the ith particle in the hth iteration, vih+1 is the particle velocity value of the ith particle in the h+1th iteration, ω is the inertia weight that reflects the extent that the particle maintains the present velocity, c1 is the cognitive coefficient that reflects the weight of the particle getting close to the individual extreme, c2 is the social coefficient which reflects the weight the particle getting close to the global extreme, r1 is a random number between 0 and 1, and r2 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:(12)Psih+1=c2⊗F3c1⊗F2ω⊗F1Psih,Pih,Pgh.In Formula (12), Psih is the fitness function value of the ith particle in the hth iteration, and Psih+1 is the fitness function value of the ith particle in the h+1th iteration. This formula includes three formulas which are written as
(13)Ψih=ω⊗F1Psih=F1Psih,Psih,rand<ω,rand≥ω,(14)Φih=c1⊗F2Ψih,Pih=F2Ψih,Pih,Ψih,rand<c1,rand≥c1,(15)Psih+1=c2⊗F3Φih,Pgh=F3Φih,Pgh,Φih,rand<c2,rand≥c2.Ψih is the value of the ith particle in the hth iteration after the calculation of ω⊗F1Psih. The formula rand means that a real number between 0 and 1 is randomly generated. When the condition rand<ω is met, the formula F1Psih is executed which means that two integers a and b are randomly generated in interval 0,N and the values on the ith particle at positions a and b are interchanged. When the condition rand≥ω is met, the formula Psih is executed which means that the coding values on the ith particle are unchanged.Φih is the value of the ith particle in the hth iteration after the calculation of c1⊗F2Ψih,Pih. Pih is the individual extreme in the tth iteration. When the condition rand<c1 is met, the formula F2Ψih,Pih is executed which means that two integers a and b are randomly generated in interval 0,N, and the values on Ψih at positions a and b as well as the values on Pih at positions a and b are interchanged. When the condition rand≥c1 is met, the formula Ψih is executed which means that the coding values on Ψih are unchanged.Pgh is the global extreme in the tth iteration. When the condition rand<c2 is met, the formula F3Φih,Pgh is executed which means that two integers a and b are randomly generated in interval 0,N, and the values on Φih at positions a and b as well as the values on Pgh at positions a and b are interchanged. When the condition rand≥c2 is met, the formula Φih is executed which means that the coding values on Φih are unchanged.The implementation steps of the DPSO algorithm are as follows.Step 1.
Set parameters: the population sizeMB of the particle swarm, the count N of OSVs, the maximum count K of satellites served by a single OSV, the maximum number of iterations TB, the inertia weight ω, the cognitive coefficient c1, the social coefficient c2, and h=1.Step 2.
The particle swarm population is initialized, the values on the particles are randomly generated, the fitness functionfPsi values of all particles are calculated, and the individual extreme Pih and global extreme Pgh are recorded.Step 3.
h=h+1, the position of the particles is updated through Formula (12), the fitness function fPsi values of all particles are calculated, and the individual extreme Pih and global extreme Pgh are updated.Step 4.
Step 3 is executed repeatedly until (Pgh−Pgh−1≤10−6) or the maximum number of iterations is reached.


## 4. 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 servesK satellites at most. If M satellites need to be served, the initial count NN of satellites can be calculated:
(16)NN=ceilMK.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.


### 4.2. Optimization Index

Given the count and orbit elements of OSVs and taking the time intervalT 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 (minJ1). 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 (maxJ2).
(17)maxJ2=maxminJ1=maxminmaxfX.In Formula (17), maxfX is optimization index B.The implementation steps of the pattern search [24] are as follows.Step 1.
i=1. The function value J1Bi at base point Bi in the time interval T is calculated. The iteration step size δt and the convergence precision ε are set.Step 2.
The temporary pointTi is determined through
(18)Ti=Bi+δt,J1Bi+δt<J1Bi,Bi−δt,J1Bi−δt<J1Bi≤J1Bi+δt,Bi,J1Bi<minJ1Bi+δt,J1Bi−δt.
Formula (18) means that the temporary point Ti is determined in point Bi−δt, point Bi+δt, and point Bi which has a minimum value of J1. The temporary point Ti is defined as base point Bi+1. ith is the base point Bi pointing to base point Bi+1.Step 3.
If the conditionJ1Ti<J1Bi is met, the i+1th vector is twice longer than the ith vector. i=i+1. The temporary point Ti is written as
(19)Ti=Bi−1+2Bi−Bi−1.
The temporary pointTi is defined as the base point Bi+1 which is taken as the starting point of the ith vector. If the condition J1Ti≥J1Bi is met, the i+1th vector is half of the ith vector and its starting point becomes the base point Bi−1.Step 4.
Steps 2 and 3 are repeatedly executed. If the value ofTi remains unchanged, the iteration step size δt will be half until it reaches a certain degree of accuracy. The algorithm ends.


### 4.3. 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 Sections2.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.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.
(20)fit¯=fit−fitminfitmax−fitmin.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, fitmin is the minimum fitness function value in the present iteration, and fitmax 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].
(21)fs=1,x<aS,bS−xbS−aS,aS≤x≤bS,0,x>bS,(22)fM=0,x<aM,2x−aMbM−aM,aM≤x≤aM+bM2,2bM−xbM−aM,aM+bM2≤x≤bM,0,x>bM,(23)fL=0,x<aL,x−aLbL−aL,aL≤x≤bL,1,x>bL.Formulas (21), (22), and (23) are the membership functions that correspond to fuzzy language “small,” “medium,” and “large,” respectively. aS and bS are the constants when two input variables ω or fit¯ are within the range of fuzzy language “small.” aM and bM are the constants when two input variables ω or fit¯ are within the range of fuzzy language “medium.” aL and bL 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(24)μRX,Yx,y=μAx∧μBy.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, RX,Y is the implication relationship between X and Y, μAx is the membership function value of x, μBy is the membership function value of y, and “∧” is the operational symbol for choosing the smaller value between μAx and μBy.The fuzzy control rule base is described in Table1.
(25)vout=∫VvμVvdv∫VμVvdv.Step 3.
The fuzzy control variableΔω is translated into the precise value after the defuzzification procedure. The precise value vout of a particle is the centre point’s x-coordinate value of the area of the curve fΔω=μRX,Yx,y in Step 2 and the x-coordinate.Table 1
The fuzzy control rule base.
Numberfit¯ωΔω1SmallSmallMedium2SmallMediumSmall3SmallLargeSmall4MediumSmallLarge5MediumMediumMedium6MediumLargeSmall7LargeSmallLarge8LargeMediumMedium9LargeLargeSmallIn Formula (25), μVv is μRX,Yx,y, V is the area of the curve fΔω=μRX,Yx,y in Step 2 and the x-coordinate, and v is the x-coordinate value. The updated inertial weight ωupdated is written as
(26)ωupdated=ω+vout.The implementation steps of the FAPSO algorithm are as follows.Step 1.
Set parameters: the population sizeMC of the particle swarm, the count N of OSVs, the maximum number of iterations TC, the inertia weight ω, the cognitive coefficient c1, the social coefficient c2, and h=1. The coding of the ith particle is Psi=xi11,⋯,xi16,⋯,xij1,⋯,xij6,⋯,xiN6,⋯,xiN6, the expression xij6 refers to the kth orbit element of the jth OSV in the ith 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 Pih and the global extreme Pgh are recorded.Step 3.
The inertial weightω and the normalized fitness function values fit¯ in the hth 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 hth iteration. The fitness function values of all particles are calculated, and the individual extreme Pih and global extreme Pgh are recorded.Step 4.
Step 3 is repeatedly executed until (Pgh−Pgh−1≤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 on-orbit mission, the cycle calculation is ended.


## 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 servesK satellites at most. If M satellites need to be served, the initial count NN of satellites can be calculated:
(16)NN=ceilMK.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.


## 4.2. Optimization Index

Given the count and orbit elements of OSVs and taking the time intervalT 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 (minJ1). 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 (maxJ2).
(17)maxJ2=maxminJ1=maxminmaxfX.In Formula (17), maxfX is optimization index B.The implementation steps of the pattern search [24] are as follows.Step 1.
i=1. The function value J1Bi at base point Bi in the time interval T is calculated. The iteration step size δt and the convergence precision ε are set.Step 2.
The temporary pointTi is determined through
(18)Ti=Bi+δt,J1Bi+δt<J1Bi,Bi−δt,J1Bi−δt<J1Bi≤J1Bi+δt,Bi,J1Bi<minJ1Bi+δt,J1Bi−δt.
Formula (18) means that the temporary point Ti is determined in point Bi−δt, point Bi+δt, and point Bi which has a minimum value of J1. The temporary point Ti is defined as base point Bi+1. ith is the base point Bi pointing to base point Bi+1.Step 3.
If the conditionJ1Ti<J1Bi is met, the i+1th vector is twice longer than the ith vector. i=i+1. The temporary point Ti is written as
(19)Ti=Bi−1+2Bi−Bi−1.
The temporary pointTi is defined as the base point Bi+1 which is taken as the starting point of the ith vector. If the condition J1Ti≥J1Bi is met, the i+1th vector is half of the ith vector and its starting point becomes the base point Bi−1.Step 4.
Steps 2 and 3 are repeatedly executed. If the value ofTi remains unchanged, the iteration step size δt will be half until it reaches a certain degree of accuracy. The algorithm ends.


## 4.3. 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 Sections2.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.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.
(20)fit¯=fit−fitminfitmax−fitmin.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, fitmin is the minimum fitness function value in the present iteration, and fitmax 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].
(21)fs=1,x<aS,bS−xbS−aS,aS≤x≤bS,0,x>bS,(22)fM=0,x<aM,2x−aMbM−aM,aM≤x≤aM+bM2,2bM−xbM−aM,aM+bM2≤x≤bM,0,x>bM,(23)fL=0,x<aL,x−aLbL−aL,aL≤x≤bL,1,x>bL.Formulas (21), (22), and (23) are the membership functions that correspond to fuzzy language “small,” “medium,” and “large,” respectively. aS and bS are the constants when two input variables ω or fit¯ are within the range of fuzzy language “small.” aM and bM are the constants when two input variables ω or fit¯ are within the range of fuzzy language “medium.” aL and bL 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(24)μRX,Yx,y=μAx∧μBy.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, RX,Y is the implication relationship between X and Y, μAx is the membership function value of x, μBy is the membership function value of y, and “∧” is the operational symbol for choosing the smaller value between μAx and μBy.The fuzzy control rule base is described in Table1.
(25)vout=∫VvμVvdv∫VμVvdv.Step 3.
The fuzzy control variableΔω is translated into the precise value after the defuzzification procedure. The precise value vout of a particle is the centre point’s x-coordinate value of the area of the curve fΔω=μRX,Yx,y in Step 2 and the x-coordinate.Table 1
The fuzzy control rule base.
Numberfit¯ωΔω1SmallSmallMedium2SmallMediumSmall3SmallLargeSmall4MediumSmallLarge5MediumMediumMedium6MediumLargeSmall7LargeSmallLarge8LargeMediumMedium9LargeLargeSmallIn Formula (25), μVv is μRX,Yx,y, V is the area of the curve fΔω=μRX,Yx,y in Step 2 and the x-coordinate, and v is the x-coordinate value. The updated inertial weight ωupdated is written as
(26)ωupdated=ω+vout.The implementation steps of the FAPSO algorithm are as follows.Step 1.
Set parameters: the population sizeMC of the particle swarm, the count N of OSVs, the maximum number of iterations TC, the inertia weight ω, the cognitive coefficient c1, the social coefficient c2, and h=1. The coding of the ith particle is Psi=xi11,⋯,xi16,⋯,xij1,⋯,xij6,⋯,xiN6,⋯,xiN6, the expression xij6 refers to the kth orbit element of the jth OSV in the ith 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 Pih and the global extreme Pgh are recorded.Step 3.
The inertial weightω and the normalized fitness function values fit¯ in the hth 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 hth iteration. The fitness function values of all particles are calculated, and the individual extreme Pih and global extreme Pgh are recorded.Step 4.
Step 3 is repeatedly executed until (Pgh−Pgh−1≤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 on-orbit mission, the cycle calculation is ended.


## 5. Simulations and Analyses of Results

### 5.1. Simulations

The initial orbit elements of all satellites are presented in Table2. The search range orbit elements (A, e, I, Ω, α, and f) of the OSVs determined according to the orbit elements of all satellites are [6700 km,11000 km], [0, 1), [0°, 180°], [0°, 360°], [0°, 360°], and [0°, 360°], respectively. The initial orbit elements of OSVs are presented in 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 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.Table 2
Initial orbit elements of satellites.
SatelliteAkmeI°α°Ω°f°S193060.23855028.321224.138307.927135.862S270590.00022598.08476.815185.081283.332S392370.00016552.006171.386350.005188.701S478280.001219101.61887.562293.441282.681S577070.05408899.03156.562347.938308.882S675540.00288290.045245.84485.215113.968S784280.20164582.982262.55529.90144.341S895010.21320028.22634.607152.421325.392S994010.01655664.822302.192337.38756.304S1099390.32163956.917231.629218.498128.371S11106870.35591757.022347.224317.457122.776S1268740.00055597.420104.012316.284204.781S1371250.00392442.08451.815255.08130.200S1488380.00024998.17869.604155.872150.009S1591230.01379352.674109.39823.270200.020S1676850.00469337.47752.707158.01760.000S1792150.48930055.371128.871305.50180.000S1874530.00421823.61877.562293.44150.101S1996750.05408999.03156.562187.938160.032S2075340.00300291.145239.73584.325113.000S2169030.00054996.390110.033325.324203.070Table 3
Initial orbit elements of OSVs.
OSVAkmeI°α°Ω°f°O177510.00012032.22356.23276.223255.520O284000.00022578.002102.332253.010122.243O385260.00412076.924154.25773.020243.000O492010.37005655.223180.922223.00195.620O570020.00402187.235102.556330.76588.598O680020.529030145.667156.007200.321270.321O799230.052301120.023254.124122.23255.220O890010.00632278.336155.33695.66345.692In the simulation of the transfer orbit optimization model, the parameters of GA are taken from literature written by Tesarova & Vokalova [12]. The population size MA is 20, the weight of the fuel used to keep the orbit mkeep is 50 kg, the maximum number of iterations TA is 50, the crossover probability Pc is 0.85, and the mutation probability Pm 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 MB is 10, the maximum number of iterations TB is 50, the inertia weight ω is 0.75, the cognitive coefficient c1 is 2, and the social coefficient c2 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 MC is 10, the maximum number of iterations TC is 50, the inertia weight ω is 0.75, the cognitive coefficient c1 is 2, and the social coefficient c2 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 (2×365×24×60×60) seconds. The parameters of the fuzzy controller are shown in Table 4.Table 4
Parameters of the fuzzy controller.
ParameterFuzzy languageabfit¯Small00.06Medium0.050.4Large0.31ωSmall0.20.6Medium0.40.9Large0.61.1ΔωSmall-0.12-0.02Medium-0.040.04Large0.000.05We 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–15]. 30 independent runs are performed for each algorithm.


### 5.2. Analyses of Results

The results of the 10th 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 (maxJ2), the time to obtain maxJ2, 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.Table 5
Orbit elements of OSVs after simulations.
OSVAkmeI°α°Ω°f°O175090.00022128.10055.221125.122305.600O281440.00019882.998115.419320.38565.768O377920.000221100.325112.366102.981200.421O488790.25002255.348200.000250.887120.133O569980.00541196.951156.052349.870149.551O673310.45223938.002200.333250.55435.010O7100490.24000053.239300.256320.457225.000O895140.15100097.000265.107145.932102.000Table 6
Results of deployment.
Time (h)OSVSatelliteJ29336.450O1S8, S1, S1821.106O2S7, S20, S6O3S2, S4O4S10, S17, S15O5S12, S21, S5O6S13, S16O7S11, S9, S3O8S19, S14Table 7
OSVs’ manoeuvre time, rendezvous time, and the weight of remaining fuel.
OSVt1st2st3st4st5st6smremainskgO151692.55867975.303104345.873109193.762193120.182198267.663199.102O2175202.944187625.628219612.145225239.147250263.743255990.129111.257O317001.794148722.117218123.229239686.355——105.759O415887.2426912.82931971.05137423.46447426.58959549.761126.350O53605.539687.42951749.98163848.707——150.445O633806.93160783.39797927.029107616.136179462.553216913.012232.828O722523.51944028.819105397.328126423.851151691.717186117.106107.141O831522.04672061.83775662.262130570.746——164.904Figure 5
Parking orbits of all OSVs.Figure 6
Transfer orbit of O1 serving S8, S1, and S18.Figure 7
Transfer orbit of O2 serving S7, S20, and S6.Figure 8
Transfer orbit of O3 serving S10, S17, and S15.Figure 9
Transfer orbit of O4 serving S12, S21, and S5.Figure 10
Transfer orbit of O5 serving S11, S9, and S3.Figure 11
Transfer orbits of O6 serving S13 and S16.Figure 12
Transfer orbit of O7 serving S19 and S14.Figure 13
Transfer orbit of O8 serving S2 and S4.Figure 14
The convergence curves of the FAPSO algorithm and other algorithms.In Table6, the time to obtain optimization index D (maxJ2) is after 9336.450 hours. O1, O2, O4, O6, and O7 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, t1 and t2 are the manoeuvre time and rendezvous time to serve the 1st satellite, respectively, t3 and t4 are the manoeuvre time and rendezvous time to serve the 2nd satellite, respectively, t5 and t6 are the manoeuvre time and rendezvous time to serve the 3rd satellite, respectively, and mremains is the weight of remaining fuel after the sth 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.


## 5.1. Simulations

The initial orbit elements of all satellites are presented in Table2. The search range orbit elements (A, e, I, Ω, α, and f) of the OSVs determined according to the orbit elements of all satellites are [6700 km,11000 km], [0, 1), [0°, 180°], [0°, 360°], [0°, 360°], and [0°, 360°], respectively. The initial orbit elements of OSVs are presented in 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 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.Table 2
Initial orbit elements of satellites.
SatelliteAkmeI°α°Ω°f°S193060.23855028.321224.138307.927135.862S270590.00022598.08476.815185.081283.332S392370.00016552.006171.386350.005188.701S478280.001219101.61887.562293.441282.681S577070.05408899.03156.562347.938308.882S675540.00288290.045245.84485.215113.968S784280.20164582.982262.55529.90144.341S895010.21320028.22634.607152.421325.392S994010.01655664.822302.192337.38756.304S1099390.32163956.917231.629218.498128.371S11106870.35591757.022347.224317.457122.776S1268740.00055597.420104.012316.284204.781S1371250.00392442.08451.815255.08130.200S1488380.00024998.17869.604155.872150.009S1591230.01379352.674109.39823.270200.020S1676850.00469337.47752.707158.01760.000S1792150.48930055.371128.871305.50180.000S1874530.00421823.61877.562293.44150.101S1996750.05408999.03156.562187.938160.032S2075340.00300291.145239.73584.325113.000S2169030.00054996.390110.033325.324203.070Table 3
Initial orbit elements of OSVs.
OSVAkmeI°α°Ω°f°O177510.00012032.22356.23276.223255.520O284000.00022578.002102.332253.010122.243O385260.00412076.924154.25773.020243.000O492010.37005655.223180.922223.00195.620O570020.00402187.235102.556330.76588.598O680020.529030145.667156.007200.321270.321O799230.052301120.023254.124122.23255.220O890010.00632278.336155.33695.66345.692In the simulation of the transfer orbit optimization model, the parameters of GA are taken from literature written by Tesarova & Vokalova [12]. The population size MA is 20, the weight of the fuel used to keep the orbit mkeep is 50 kg, the maximum number of iterations TA is 50, the crossover probability Pc is 0.85, and the mutation probability Pm 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 MB is 10, the maximum number of iterations TB is 50, the inertia weight ω is 0.75, the cognitive coefficient c1 is 2, and the social coefficient c2 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 MC is 10, the maximum number of iterations TC is 50, the inertia weight ω is 0.75, the cognitive coefficient c1 is 2, and the social coefficient c2 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 (2×365×24×60×60) seconds. The parameters of the fuzzy controller are shown in Table 4.Table 4
Parameters of the fuzzy controller.
ParameterFuzzy languageabfit¯Small00.06Medium0.050.4Large0.31ωSmall0.20.6Medium0.40.9Large0.61.1ΔωSmall-0.12-0.02Medium-0.040.04Large0.000.05We 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–15]. 30 independent runs are performed for each algorithm.


## 5.2. Analyses of Results

The results of the 10th 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 (maxJ2), the time to obtain maxJ2, 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.Table 5
Orbit elements of OSVs after simulations.
OSVAkmeI°α°Ω°f°O175090.00022128.10055.221125.122305.600O281440.00019882.998115.419320.38565.768O377920.000221100.325112.366102.981200.421O488790.25002255.348200.000250.887120.133O569980.00541196.951156.052349.870149.551O673310.45223938.002200.333250.55435.010O7100490.24000053.239300.256320.457225.000O895140.15100097.000265.107145.932102.000Table 6
Results of deployment.
Time (h)OSVSatelliteJ29336.450O1S8, S1, S1821.106O2S7, S20, S6O3S2, S4O4S10, S17, S15O5S12, S21, S5O6S13, S16O7S11, S9, S3O8S19, S14Table 7
OSVs’ manoeuvre time, rendezvous time, and the weight of remaining fuel.
OSVt1st2st3st4st5st6smremainskgO151692.55867975.303104345.873109193.762193120.182198267.663199.102O2175202.944187625.628219612.145225239.147250263.743255990.129111.257O317001.794148722.117218123.229239686.355——105.759O415887.2426912.82931971.05137423.46447426.58959549.761126.350O53605.539687.42951749.98163848.707——150.445O633806.93160783.39797927.029107616.136179462.553216913.012232.828O722523.51944028.819105397.328126423.851151691.717186117.106107.141O831522.04672061.83775662.262130570.746——164.904Figure 5
Parking orbits of all OSVs.Figure 6
Transfer orbit of O1 serving S8, S1, and S18.Figure 7
Transfer orbit of O2 serving S7, S20, and S6.Figure 8
Transfer orbit of O3 serving S10, S17, and S15.Figure 9
Transfer orbit of O4 serving S12, S21, and S5.Figure 10
Transfer orbit of O5 serving S11, S9, and S3.Figure 11
Transfer orbits of O6 serving S13 and S16.Figure 12
Transfer orbit of O7 serving S19 and S14.Figure 13
Transfer orbit of O8 serving S2 and S4.Figure 14
The convergence curves of the FAPSO algorithm and other algorithms.In Table6, the time to obtain optimization index D (maxJ2) is after 9336.450 hours. O1, O2, O4, O6, and O7 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, t1 and t2 are the manoeuvre time and rendezvous time to serve the 1st satellite, respectively, t3 and t4 are the manoeuvre time and rendezvous time to serve the 2nd satellite, respectively, t5 and t6 are the manoeuvre time and rendezvous time to serve the 3rd satellite, respectively, and mremains is the weight of remaining fuel after the sth 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.


## 6. 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.



---

*Source: 6644339-2021-06-21.xml*