# Cost-Efficient LEO Navigation Augmentation Constellation Design under a Constrained Deployment Approach

Authors: Jing Ren; Dan Sun; Deng Pan; Mingtao Li; Jianhua Zheng
Journal: International Journal of Aerospace Engineering
 (2021)

---

## Abstract

The advantages of the Low Earth Orbit (LEO) satellite include low-latency communications, shorter positioning time, higher positioning accuracy, and lower launching, building, and maintenance costs. Thus, the introduction of LEO satellite constellation as a regional navigation augmentation system for the current navigation constellations is studied in this paper. To achieve the navigation performance requirement with the least system cost, a synthetic approach is presented to design and deploy a cost-efficient LEO navigation augmentation constellation over 108 key cities. To achieve lower construction costs, the constellation is designed to be deployed by constrained piggyback launches, which brings additional complexity to the constellation design. Two optimization models with discrete and continuous performance indices are established. They are solved by the genetic algorithm and differential evolution algorithm, and both Walker and Flower constellations are adopted. Results for 77 and 70 satellites are obtained. During the construction phase, a synthesis procedure containing five impulses is proposed by utilizing natural drift underJ2 perturbation. This work presents a method for designing the optimal LEO navigation constellation under a constraint deployment approach with the lowest construction cost and a strategy to deploy the constellation economically.


---

## Body

## 1. Introduction

Navigation constellation design has received interest for years. Most existing navigation constellations are deployed in Medium Earth orbits and geosynchronous Earth orbits, such as Global Positioning System (GPS), Global Navigation Satellite System (GLONASS), and GALILEO constellations (the orbital data is obtained fromhttp://delestrak.com). Resent research presents that a Low Earth Orbit (LEO) satellite constellation provides a significant advantage over a geostationary (GEO) satellite in terms of low-latency communications, shorter positioning time, higher positioning accuracy, and lower launching, building, and maintenance costs [1–3]. Thus, introducing Low Earth Orbit (LEO) constellation for navigation augmentation can reduce the cost of mission, shorten the position time, and improve the positioning accuracy. Therefore, LEO navigation augmentation has attracted a great deal of research interest.There has been much deep investigation focusing on LEO constellation design, and various constellation design techniques have been proposed, such as streets of coverage (SOC) [4], the Walker constellation theory [5, 6], and the Flower constellation theory [7]. In the Walker constellation, the orbital parameters are dependent on one another in a particular way, thereby reducing the complexity. The Walker technique provides the most symmetric geometry among all constellation design techniques. Thus, it is most suitable for navigation systems that are related to global coverage. To facilitate the use of on-orbit active spares, Patterson further developed the Walker theory and then proposed the Teledesic model, in which phasing between satellites in adjacent planes is not required [8]. Moreover, Palmade et al. combined two identical Walker constellations and proposed the SkyBridge constellation, in which the architecture of the constellation is made of 64 satellites distributed in two subconstellations shifted by -10° in the longitude of the ascending node and +14° in mean anomaly. This configuration offers better coverage performance with fewer satellites [9]. The Flower constellation theory [7, 10] was initially proposed as an alternative design strategy to the previous theory. By removing some design constraints, Flower constellations have been demonstrated to have better navigation performance than the existing GPS, GLONASS, and GALILEO using the same numbers of satellites or achieve the same performance using fewer satellites [11–13]. In recent years, this theory has also been applied to the constellation design problem domains including navigation, Earth observation space-based systems, regional coverage, and reconnaissance [14]. Remarkably, Lee and Ho presented an approach to design the Flower constellation patterns for regional coverage with predefined seed-satellite elements [15, 16]. In this paper, a design strategy for full constellations is presented, which consists the seed-satellite elements and configuration parameters under constrained deployment approaches.In addition to the constellation design, optimized deployment can significantly reduce the system cost. Budianto and Olds noted that the cost can be decreased by collaboratively optimizing the configuration, orbit design, spacecraft design, and launch manifest [17]. Furthermore, de Weck et al. presented an approach for deploying a constellation progressively and finding the best reconfigurable constellations within a given design space during deployment, thereby obtaining significant economic benefits [18]. Additionally, Lee et al. proposed a design method to optimize a multistage configuration for multiple possible scenarios with minimal lifecycle cost [19]. In addition, GMV analyzed the possible strategies for most of the fundamental phases of the constellation life cycle, including constellation launch, setup, replacement of failed satellites, and end-of-life policy [20]. The authors suggested that orbit transfer can be performed by taking advantage of the secular drift of the right ascension of the ascending node (RAAN) due to the J2-term of the terrestrial gravitational potential. This effect can be used to create a relative precession motion between two orbital planes to obtain the desired RAAN separation. This strategy allows the amount of fuel allocated to the orbit transfer to be reduced since it avoids out-of-plane maneuvers that are generally expensive.Based on Problem B of the 9th China Trajectory Optimization Competition (CTOC9 B) [21], several researchers made efforts to design a suitable constellation configuration and optimize the construction cost to satisfy the required navigation performance for selected key cities in China. This paper presents a synthetic approach to design and deploy a cost-efficient LEO regional navigation augmentation constellation. For the constellation design, in previous, the effectiveness of navigation constellation is generally considered a discrete index [15, 16, 21]. In this paper, an optimization model with a continuous performance index is developed to replace the discrete performance index and improve convergence. To validate this model, the Walker constellation and the Flower constellation are employed, and constellations that achieve the navigation performance requirement with fewer satellites are obtained. For the constellation deployment under predefined mission duration and limited fuel, a five-impulse transfer strategy is proposed to save fuel. A small tangential impulse at the start of the construction phase is applied to adjust the initial state of the satellite to achieve the proper phase after a period of natural drift. Moreover, an entire optimization procedure for the transfer process is presented. Overall, our investigation provides a method to design the optimal LEO navigation constellation under a constraint deployment approach with the fewest satellites and a strategy to deploy the constellation economically.


## 2. Background

### 2.1. Problem Description

Problem B of the 9th CTOC [21] focuses on the design and deployment of a regional augmentation system offering navigation services to 108 key cities in China. The goal is to design a suitable constellation configuration and optimize the navigation performance of key cities under the conditions of limited construction costs.Figure1 shows the positions of 108 key cities in China. Additionally, 100 available piggyback launches are provided; detailed information can be found in [21]. The constellation design process can be divided into two phases: construction and service. The construction phase is from Modified Julian Day 2000 (MJD2000) 7305 to 7396. During the construction phase, the navigation constellation can be constructed by new launches or piggyback launches. In the first approach, the distributer carries at most 16 navigation satellites and is located in a circular parking orbit with an altitude of 900 km. In the second approach, the distributer carries at most 8 navigation satellites and is located in any orbit of the piggyback launches. After that, for the two deployment methods stated above, the navigation satellites are separated from the distributer and are transferred to target orbits with chemical propulsion. If an elliptical orbit is employed, a critical inclination orbit at approximately 63.4° should be adopted. At MJD2000=7396, the constellation construction ends, and the service phase begins. The constellation service phase starts at MJD2000=7396. Then, the navigation performance of the constellation for all cities is evaluated on the 1st, 7th, and 30th days (i.e., MJD2000=7396to7397, MJD2000=7402to7403, and MJD2000=7425to7426) with a fixed interval Δt=120s.Figure 1
Locations of the 108 key cities.The primary indexObj1 is aimed at maximizing the total navigation revenue of all key cities:
(1)Obj1=max∑k=1Kδkwkδk=1,ifmaxGDOPk,l<100,ifmaxGDOPk,l>10,l=1,2,⋯,L;k=1,2,⋯,K,where δk is a 0-1 variable, which is used to measure whether the kth city can meet the navigation accuracy requirements, and wk is the weight of the kth city. GDOPk,l is the geometric dilution of precision (GDOP) of the kth city at the lth sampling point. Ideally, if the GDOP values of all the cities at all sampling points are below 10, Obj1 is equal to 229, which is the maximum value of Obj1. The secondary index Obj2 is aimed at minimizing the entire cost of constellation construction. The construction phase is from MJD2000=7305to7396. During this phase, the navigation constellation can be constructed by new launches or piggyback launches.
(2)Obj2=minNlaunchClaunch+NcarryCcarry+NnaviSatCnaviSat,where Nlaunch, Ncarry, and NnaviSat are the number of new launches, the number of piggyback launches, and the number of navigation satellites in the constellation, respectively. Claunch=1.2, Ccarry=0.2, and CnaviSat=0.05 are the cost of a new launch, the cost of a piggyback launch, and the cost of a navigation satellite, respectively, denoted in monetary units.During the construction and service phases, Earth’s central gravitational field and theJ2 perturbation are considered. The dynamic models are described as follows:
(3)x¨=−μxr31+32J2Rer21−5z2r2y¨=−μyr31+32J2Rer21−5z2r2z¨=−μzr31+32J2Rer23−5z2r2,where μ is the gravitational constant of the Earth, x,y,z are the position components of the satellite in the Earth-centered inertial (ECI) coordinate system, r=x2+y2+z2 is the geocentric distance, Re is the equatorial radius, and J2 is the Earth oblateness gravity coefficient.During the service phase, mean orbital elements are implemented to propagate the constellation and evaluate the GDOP. The dynamic models of mean elements underJ2 perturbation are described as follows:
(4)ȧ=0ė=0inc=0Ω̇=−3AnJ221−e22Rea2cosiω̇=3AnJ241−e22Rea25cos2i−1Ṁ=3AnJ241−e23/2Rea23cos2i−1,where a is the semimajor axis, e is the eccentricity, inc is the inclination, Ω denotes the RAAN, ω denotes the argument of the perigee, M denotes the mean anomaly, and An=μ/a3 is the mean motion of the orbit. When Ω0,ω0, and M0 are given at t0, the elements at time t can be calculated by
(5)Ω=Ω0+Ω̇t−t0ω=ω0+ω̇t−t0M=M0+An+Ṁt−t0.During the construction phase, the satellites are propelled by a chemical engine, and the residual massmf after an impulse is given by
(6)mf=miexp−ΔvIspg0,where mi is the satellite mass before the impulse, Δv is the magnitude of the velocity impulse, Isp is the specific impulse, and g0 is the gravitational acceleration at sea level. All the parameters are set according to [21].


### 2.2. Walker Constellation Theory

Based on the contributions of Walker [5, 6], a Walker constellation includes the star pattern, δ pattern, σ pattern, omega pattern, and rosette pattern. The σ pattern and rosette pattern are individual cases of Walker-δ patterns. Among all of these patterns, the Walker-δ pattern is the most widely used. According to the problem, if an elliptical orbit is employed, a critical inclination orbit should be adopted (63.4° for a prograde orbit; 116.6° for a retrograde orbit); however, there are few piggyback missions around this critical inclination, which can lead to difficulties in deployment. Thus, our work is based on a circular Walker constellation. The Walker constellation below refers to the Walker-δ pattern. The Walker-δ pattern consists of a group of satellites that share the same orbit radius and inclination, and these satellites are distributed evenly in several orbital planes. Walker developed a notation that is commonly used as a starting point for constellation design for labeling the orbits, and the Walker-δ pattern can be identified by T/P/F, where T is the total number of satellites, P is the total number of orbital planes, and F is the relative spacing between satellites in adjacent planes, which may assume any value between 0 and P−1. The number of satellites per plane, N, is given as N=T/P. The ascending nodes of all orbital planes are distributed evenly around the equator at intervals of 360°/P. Within each orbital plane, the N satellites are distributed at intervals of 360°/N. The only remaining issue is the relative phase between satellites in adjacent orbital planes. To identify this, Walker defined the phase difference (DF) in a constellation as the angle of the direction of motion from the ascending node of a satellite in one place to the nearest satellite in the next most westerly plane. The relative angular shift between satellites in adjacent orbital planes is equal to F·360°/T. The relevant equations are as follows:
(7)Ωij0=Ω0+i−12πPMij0=M0+i−1·F·2πT+j−1·P·2πT,where i=1,2,⋯,P;j=1,2,⋯,N, and Ω0,ω0 are the initial RAAN and the initial mean anomaly of the starting point, respectively. Ωij0 and Mij0 are the RAAN and mean anomaly, respectively, of the jth satellite in the ith orbital plane at the initial time t=0. Thus, a specific Walker-δ pattern can be uniquely determined by 7 parameters: orbital radius a, inclination inc, total number of satellites T, number of orbital planes P, initial RAAN value, and mean anomaly value of the starting point Ω0 and M0, and phase factor F.In this work, a quasi-Walker constellation model is established, which is a circular Walker-δ pattern that is distributed in a restricted region. The satellites are evenly deployed into several orbital planes with a certain separation angle in the RAAN between adjacent orbital planes instead of the orbital planes being evenly distributed around the equator. Therefore, a quasi-Walker constellation can be uniquely defined by
(8)Satijt0=a,0,inc,Ωijt0,0,Mijt0Ωijt0=Ω00+i−1ΔΩPMijt0=M00+i−1·F·ΔΩT+j−1·P·2πT,where Satijt0, Ωijt0, and Mijt0 denote the Kepler orbit elements, the RAAN, and the mean anomaly of the jth satellite in the ith orbital plane at the initial time t=t0,i=1,2,⋯,P,j=1,2,⋯,N. In this work, t0=MJD2000=7396. P,N,a,inc,Ω00,ΔΩ, and F denote the number of orbit planes, number of satellites in each plane, orbit radius, inclination, initial RAAN of the starting point, RAAN span of the constellation, and phase factor, respectively. In the typical Walker-δ pattern, the phase factor F is set to an integer value. In this work, F is defined as a continuous variable, which can enlarge the search space significantly and contribute to a better solution.


### 2.3. Flower Constellation Theory

The Flower constellation [7] is built by satellites that follow the same relative trajectory in the Earth-centered Earth-fixed (ECEF) coordinate system, so the orbits in a Flower constellation yield a perfectly repeated ground track. According to the above analysis, circular restricted Flower constellations [7] are adopted in this work. A circular Flower constellation can be identified by the following parameters: the number of petals Np, number of sidereal days to repeat the ground track Nd, number of satellites Ns, phasing parameters Ω00,ΔΩ,M00, orbital radius a, and orbit inclination inc (a and inc are orbit parameters and equal for all satellites). To obtain a repeating ground track, which will repeat after the satellite completes Np revolutions over Nd days, the following relationship is required
(9)NpTΩ=NdTΩG,where TΩ and TΩG are the nodal period of the orbit and nodal period of Greenwich, respectively. In this paper, Nd is set to 1 according to the mission requirements. The orbits in a Flower constellation can be defined by the following equations [7]:
(10)Satit0=a,0,inc,Ωit0,0,Mit0Ωit0=Ω00+i−1ΔΩNp−1Mit0=M00−i−1·n+Ṁ0ωE+Ω̇·ΔΩNp−1,where
(11)Ω̇=−2ξncosincṀ0=−ξn3sin2inc−2ξ=3Re2J24p2,where Satit0, Ωit0, and Mit0 are the Kepler orbit elements, RAAN, and mean anomaly of the ith satellite at the initial time t=t0,i=1,2,⋯,Np; a,inc,Ω00,ΔΩ, and M00 are the orbit radius, inclination, initial RAAN of the starting point, RAAN span of the constellation, and mean anomaly of the starting point, respectively; n=μ/a3 is the mean motion; ωE is the Earth spin rate.


## 2.1. Problem Description

Problem B of the 9th CTOC [21] focuses on the design and deployment of a regional augmentation system offering navigation services to 108 key cities in China. The goal is to design a suitable constellation configuration and optimize the navigation performance of key cities under the conditions of limited construction costs.Figure1 shows the positions of 108 key cities in China. Additionally, 100 available piggyback launches are provided; detailed information can be found in [21]. The constellation design process can be divided into two phases: construction and service. The construction phase is from Modified Julian Day 2000 (MJD2000) 7305 to 7396. During the construction phase, the navigation constellation can be constructed by new launches or piggyback launches. In the first approach, the distributer carries at most 16 navigation satellites and is located in a circular parking orbit with an altitude of 900 km. In the second approach, the distributer carries at most 8 navigation satellites and is located in any orbit of the piggyback launches. After that, for the two deployment methods stated above, the navigation satellites are separated from the distributer and are transferred to target orbits with chemical propulsion. If an elliptical orbit is employed, a critical inclination orbit at approximately 63.4° should be adopted. At MJD2000=7396, the constellation construction ends, and the service phase begins. The constellation service phase starts at MJD2000=7396. Then, the navigation performance of the constellation for all cities is evaluated on the 1st, 7th, and 30th days (i.e., MJD2000=7396to7397, MJD2000=7402to7403, and MJD2000=7425to7426) with a fixed interval Δt=120s.Figure 1
Locations of the 108 key cities.The primary indexObj1 is aimed at maximizing the total navigation revenue of all key cities:
(1)Obj1=max∑k=1Kδkwkδk=1,ifmaxGDOPk,l<100,ifmaxGDOPk,l>10,l=1,2,⋯,L;k=1,2,⋯,K,where δk is a 0-1 variable, which is used to measure whether the kth city can meet the navigation accuracy requirements, and wk is the weight of the kth city. GDOPk,l is the geometric dilution of precision (GDOP) of the kth city at the lth sampling point. Ideally, if the GDOP values of all the cities at all sampling points are below 10, Obj1 is equal to 229, which is the maximum value of Obj1. The secondary index Obj2 is aimed at minimizing the entire cost of constellation construction. The construction phase is from MJD2000=7305to7396. During this phase, the navigation constellation can be constructed by new launches or piggyback launches.
(2)Obj2=minNlaunchClaunch+NcarryCcarry+NnaviSatCnaviSat,where Nlaunch, Ncarry, and NnaviSat are the number of new launches, the number of piggyback launches, and the number of navigation satellites in the constellation, respectively. Claunch=1.2, Ccarry=0.2, and CnaviSat=0.05 are the cost of a new launch, the cost of a piggyback launch, and the cost of a navigation satellite, respectively, denoted in monetary units.During the construction and service phases, Earth’s central gravitational field and theJ2 perturbation are considered. The dynamic models are described as follows:
(3)x¨=−μxr31+32J2Rer21−5z2r2y¨=−μyr31+32J2Rer21−5z2r2z¨=−μzr31+32J2Rer23−5z2r2,where μ is the gravitational constant of the Earth, x,y,z are the position components of the satellite in the Earth-centered inertial (ECI) coordinate system, r=x2+y2+z2 is the geocentric distance, Re is the equatorial radius, and J2 is the Earth oblateness gravity coefficient.During the service phase, mean orbital elements are implemented to propagate the constellation and evaluate the GDOP. The dynamic models of mean elements underJ2 perturbation are described as follows:
(4)ȧ=0ė=0inc=0Ω̇=−3AnJ221−e22Rea2cosiω̇=3AnJ241−e22Rea25cos2i−1Ṁ=3AnJ241−e23/2Rea23cos2i−1,where a is the semimajor axis, e is the eccentricity, inc is the inclination, Ω denotes the RAAN, ω denotes the argument of the perigee, M denotes the mean anomaly, and An=μ/a3 is the mean motion of the orbit. When Ω0,ω0, and M0 are given at t0, the elements at time t can be calculated by
(5)Ω=Ω0+Ω̇t−t0ω=ω0+ω̇t−t0M=M0+An+Ṁt−t0.During the construction phase, the satellites are propelled by a chemical engine, and the residual massmf after an impulse is given by
(6)mf=miexp−ΔvIspg0,where mi is the satellite mass before the impulse, Δv is the magnitude of the velocity impulse, Isp is the specific impulse, and g0 is the gravitational acceleration at sea level. All the parameters are set according to [21].


## 2.2. Walker Constellation Theory

Based on the contributions of Walker [5, 6], a Walker constellation includes the star pattern, δ pattern, σ pattern, omega pattern, and rosette pattern. The σ pattern and rosette pattern are individual cases of Walker-δ patterns. Among all of these patterns, the Walker-δ pattern is the most widely used. According to the problem, if an elliptical orbit is employed, a critical inclination orbit should be adopted (63.4° for a prograde orbit; 116.6° for a retrograde orbit); however, there are few piggyback missions around this critical inclination, which can lead to difficulties in deployment. Thus, our work is based on a circular Walker constellation. The Walker constellation below refers to the Walker-δ pattern. The Walker-δ pattern consists of a group of satellites that share the same orbit radius and inclination, and these satellites are distributed evenly in several orbital planes. Walker developed a notation that is commonly used as a starting point for constellation design for labeling the orbits, and the Walker-δ pattern can be identified by T/P/F, where T is the total number of satellites, P is the total number of orbital planes, and F is the relative spacing between satellites in adjacent planes, which may assume any value between 0 and P−1. The number of satellites per plane, N, is given as N=T/P. The ascending nodes of all orbital planes are distributed evenly around the equator at intervals of 360°/P. Within each orbital plane, the N satellites are distributed at intervals of 360°/N. The only remaining issue is the relative phase between satellites in adjacent orbital planes. To identify this, Walker defined the phase difference (DF) in a constellation as the angle of the direction of motion from the ascending node of a satellite in one place to the nearest satellite in the next most westerly plane. The relative angular shift between satellites in adjacent orbital planes is equal to F·360°/T. The relevant equations are as follows:
(7)Ωij0=Ω0+i−12πPMij0=M0+i−1·F·2πT+j−1·P·2πT,where i=1,2,⋯,P;j=1,2,⋯,N, and Ω0,ω0 are the initial RAAN and the initial mean anomaly of the starting point, respectively. Ωij0 and Mij0 are the RAAN and mean anomaly, respectively, of the jth satellite in the ith orbital plane at the initial time t=0. Thus, a specific Walker-δ pattern can be uniquely determined by 7 parameters: orbital radius a, inclination inc, total number of satellites T, number of orbital planes P, initial RAAN value, and mean anomaly value of the starting point Ω0 and M0, and phase factor F.In this work, a quasi-Walker constellation model is established, which is a circular Walker-δ pattern that is distributed in a restricted region. The satellites are evenly deployed into several orbital planes with a certain separation angle in the RAAN between adjacent orbital planes instead of the orbital planes being evenly distributed around the equator. Therefore, a quasi-Walker constellation can be uniquely defined by
(8)Satijt0=a,0,inc,Ωijt0,0,Mijt0Ωijt0=Ω00+i−1ΔΩPMijt0=M00+i−1·F·ΔΩT+j−1·P·2πT,where Satijt0, Ωijt0, and Mijt0 denote the Kepler orbit elements, the RAAN, and the mean anomaly of the jth satellite in the ith orbital plane at the initial time t=t0,i=1,2,⋯,P,j=1,2,⋯,N. In this work, t0=MJD2000=7396. P,N,a,inc,Ω00,ΔΩ, and F denote the number of orbit planes, number of satellites in each plane, orbit radius, inclination, initial RAAN of the starting point, RAAN span of the constellation, and phase factor, respectively. In the typical Walker-δ pattern, the phase factor F is set to an integer value. In this work, F is defined as a continuous variable, which can enlarge the search space significantly and contribute to a better solution.


## 2.3. Flower Constellation Theory

The Flower constellation [7] is built by satellites that follow the same relative trajectory in the Earth-centered Earth-fixed (ECEF) coordinate system, so the orbits in a Flower constellation yield a perfectly repeated ground track. According to the above analysis, circular restricted Flower constellations [7] are adopted in this work. A circular Flower constellation can be identified by the following parameters: the number of petals Np, number of sidereal days to repeat the ground track Nd, number of satellites Ns, phasing parameters Ω00,ΔΩ,M00, orbital radius a, and orbit inclination inc (a and inc are orbit parameters and equal for all satellites). To obtain a repeating ground track, which will repeat after the satellite completes Np revolutions over Nd days, the following relationship is required
(9)NpTΩ=NdTΩG,where TΩ and TΩG are the nodal period of the orbit and nodal period of Greenwich, respectively. In this paper, Nd is set to 1 according to the mission requirements. The orbits in a Flower constellation can be defined by the following equations [7]:
(10)Satit0=a,0,inc,Ωit0,0,Mit0Ωit0=Ω00+i−1ΔΩNp−1Mit0=M00−i−1·n+Ṁ0ωE+Ω̇·ΔΩNp−1,where
(11)Ω̇=−2ξncosincṀ0=−ξn3sin2inc−2ξ=3Re2J24p2,where Satit0, Ωit0, and Mit0 are the Kepler orbit elements, RAAN, and mean anomaly of the ith satellite at the initial time t=t0,i=1,2,⋯,Np; a,inc,Ω00,ΔΩ, and M00 are the orbit radius, inclination, initial RAAN of the starting point, RAAN span of the constellation, and mean anomaly of the starting point, respectively; n=μ/a3 is the mean motion; ωE is the Earth spin rate.


## 3. Problem Analysis

The deployment of a constellation with constrained piggyback launches brings additional complexity to the constellation design. The search space should be reduced in the optimization of the constellation configuration to ensure that the constellation can be constructed by piggyback launches.


### 3.1. Constraint on Inclination

Specifically, to reduce the fuel consumption of the satellites, the inclination of the constellation is expected to be equal to that of some piggyback launches. Thus, the inclinations of the given 100 piggyback launches are rounded to the nearest whole number and counted in Figure2. As Figure 2 shows, among the 100 piggyback launches, the number of those inclined at angles of approximately 74°, 84°, and 98° with respect to the equator is the three largest.Figure 2
Histogram of the inclinations of the piggyback launches.Figure3 shows how the primary index Obj1 changes with increasing inclination. Four configurations are used for the experiments, which differ only in the orbit radius. Their basic configuration is 100/20/1,e=0,Ω00=0,M00=0, and the span of the RAAN ΔΩ is set to 360°. The results are shown in Table 1 and Figure 3. The maximum inclination for the best Obj1 value increases with increasing orbit radius, and the Obj1 value decreases rapidly with an increase in the inclination when it is larger than the maximum inclination for the best Obj1 value. Additionally, the range of the optimal inclinations grows broader with increasing orbit radius. Notably, the configuration with an orbital radius of 9250 km provides high-grade navigation performance at approximately 72°, while the other configurations that have lower orbital heights cannot achieve the maximum Obj1 at that angle. That is, if high-grade navigation performance is expected at a high inclination, the orbit radius should be sufficiently large. Therefore, 74° is chosen as the fundamental inclination in designing constellations, and the piggyback launches whose inclinations can be rounded to 74° are considered the candidate piggyback launches.Figure 3
Effect of the inclination onObj1.Table 1
Navigation performance of the four configurations.
Constellation indexOrbit radius (km)Minimum inclination for the bestObj1 value (deg)Maximum inclination for the bestObj1 value (deg)185005757287504168390004064492503872


### 3.2. Constraints on Orbital Height

The coplanar maneuverability of the satellite restricts the orbital height of the constellation. Based on the Hohmann transfer theory, the reachable maximum orbit radii of the candidate piggyback launches are estimated, assuming that all fuel is consumed to increase the height. The results are shown in Table2. Detailed information of the piggyback launches, including the mission indices and orbital elements, can be found in [21]. The values of the reachable maximum orbit radii of these piggyback launches are approximately 9400-9500 km except for mission 3 and mission 32. This range of values is used as the boundary for the search space of the orbit radius.Table 2
Maximum orbit radius estimations of the piggyback launches that are inclined at angles of approximately 74°.
Mission indexOrbit radius (km)Inclination (deg)Maximum orbit radius (km)137802.48948473.984000019462.087485927818.91111774.000400019484.031118447845.74114974.002200029519.899150987822.79526674.006200019489.222267107864.36017474.006700039544.800175537859.11106874.008699999537.779069337850.45956374.012200029526.208564617861.00530974.012600009540.313310757828.74211474.015099999497.172115487831.93553074.015900019501.440531307824.72411574.019000009491.80111687850.09894074.028499989525.725941477818.77659274.029099989483.851593827822.14526574.034599989488.35426637157.95139774.039099998606.443398327159.34382574.039099998608.27982667762.00892074.043999989408.024921


### 3.3. Constraints on RAANs

For the distribution of orbital planes, Ulybyshev developed a method to analyze the coverage requirements of satellite constellations [22, 23]. For a satellite located at orbital altitude h, the projection of its footprint onto Earth’s surface defines a circle of coverage with angle Θ as follows:
(12)Θ=cos−1ReRe+hcosϵ−ϵwhere ϵ denotes the minimum elevation angle, which is set to 10° in this problem. Ulybyshev defined a two-dimensional space associated with the satellite: the x axis is the RAAN, and the y axis is time. The visibility conditions for the satellite in this space can be represented by a polygon of boundary points. Figure 4 shows a map of coverage polygons for the 108 key cities, which are associated with the satellite at a=9450 km, inc=74°. Therefore, the satellites should be distributed in the span of the RAANs at 60°-350°.Figure 4
Coverage polygons for the 108 key cities,a=9450 km, inc=74°.Furthermore, the boundary of the RAANs of the orbital planes is also constrained within the RAAN drift range of the candidate piggyback launches, as the RAAN difference can be minimized withJ2 perturbation during construction. The RAAN drift range of each candidate piggyback launch during the construction phase is evaluated by Equation (5). The analysis results are shown in Figure 5. Each horizontal line segment represents the RAAN drift range of the corresponding piggyback launch during the construction phase, and the corresponding mission indices are labeled above the line segments. Detailed information of the piggyback launches, including the mission indices and orbital elements, can be found in [21]. The RAAN drift range of each orbit is approximately 60°, and an orbital plane at an RAAN of 100-380° can be covered by multiple piggyback launches. Thus, the initial RAAN value of the starting point Ω00 is set to 100-160°, and the orbital planes are expected to be distributed in Ω∈100°,375°. Therefore, the span of the RAAN ΔΩ is set to 120°,255°−Ω00.Figure 5
Accessible RAAN range of the piggyback launches that are inclined at angles of approximately 74° under theJ2 perturbation.


## 3.1. Constraint on Inclination

Specifically, to reduce the fuel consumption of the satellites, the inclination of the constellation is expected to be equal to that of some piggyback launches. Thus, the inclinations of the given 100 piggyback launches are rounded to the nearest whole number and counted in Figure2. As Figure 2 shows, among the 100 piggyback launches, the number of those inclined at angles of approximately 74°, 84°, and 98° with respect to the equator is the three largest.Figure 2
Histogram of the inclinations of the piggyback launches.Figure3 shows how the primary index Obj1 changes with increasing inclination. Four configurations are used for the experiments, which differ only in the orbit radius. Their basic configuration is 100/20/1,e=0,Ω00=0,M00=0, and the span of the RAAN ΔΩ is set to 360°. The results are shown in Table 1 and Figure 3. The maximum inclination for the best Obj1 value increases with increasing orbit radius, and the Obj1 value decreases rapidly with an increase in the inclination when it is larger than the maximum inclination for the best Obj1 value. Additionally, the range of the optimal inclinations grows broader with increasing orbit radius. Notably, the configuration with an orbital radius of 9250 km provides high-grade navigation performance at approximately 72°, while the other configurations that have lower orbital heights cannot achieve the maximum Obj1 at that angle. That is, if high-grade navigation performance is expected at a high inclination, the orbit radius should be sufficiently large. Therefore, 74° is chosen as the fundamental inclination in designing constellations, and the piggyback launches whose inclinations can be rounded to 74° are considered the candidate piggyback launches.Figure 3
Effect of the inclination onObj1.Table 1
Navigation performance of the four configurations.
Constellation indexOrbit radius (km)Minimum inclination for the bestObj1 value (deg)Maximum inclination for the bestObj1 value (deg)185005757287504168390004064492503872


## 3.2. Constraints on Orbital Height

The coplanar maneuverability of the satellite restricts the orbital height of the constellation. Based on the Hohmann transfer theory, the reachable maximum orbit radii of the candidate piggyback launches are estimated, assuming that all fuel is consumed to increase the height. The results are shown in Table2. Detailed information of the piggyback launches, including the mission indices and orbital elements, can be found in [21]. The values of the reachable maximum orbit radii of these piggyback launches are approximately 9400-9500 km except for mission 3 and mission 32. This range of values is used as the boundary for the search space of the orbit radius.Table 2
Maximum orbit radius estimations of the piggyback launches that are inclined at angles of approximately 74°.
Mission indexOrbit radius (km)Inclination (deg)Maximum orbit radius (km)137802.48948473.984000019462.087485927818.91111774.000400019484.031118447845.74114974.002200029519.899150987822.79526674.006200019489.222267107864.36017474.006700039544.800175537859.11106874.008699999537.779069337850.45956374.012200029526.208564617861.00530974.012600009540.313310757828.74211474.015099999497.172115487831.93553074.015900019501.440531307824.72411574.019000009491.80111687850.09894074.028499989525.725941477818.77659274.029099989483.851593827822.14526574.034599989488.35426637157.95139774.039099998606.443398327159.34382574.039099998608.27982667762.00892074.043999989408.024921


## 3.3. Constraints on RAANs

For the distribution of orbital planes, Ulybyshev developed a method to analyze the coverage requirements of satellite constellations [22, 23]. For a satellite located at orbital altitude h, the projection of its footprint onto Earth’s surface defines a circle of coverage with angle Θ as follows:
(12)Θ=cos−1ReRe+hcosϵ−ϵwhere ϵ denotes the minimum elevation angle, which is set to 10° in this problem. Ulybyshev defined a two-dimensional space associated with the satellite: the x axis is the RAAN, and the y axis is time. The visibility conditions for the satellite in this space can be represented by a polygon of boundary points. Figure 4 shows a map of coverage polygons for the 108 key cities, which are associated with the satellite at a=9450 km, inc=74°. Therefore, the satellites should be distributed in the span of the RAANs at 60°-350°.Figure 4
Coverage polygons for the 108 key cities,a=9450 km, inc=74°.Furthermore, the boundary of the RAANs of the orbital planes is also constrained within the RAAN drift range of the candidate piggyback launches, as the RAAN difference can be minimized withJ2 perturbation during construction. The RAAN drift range of each candidate piggyback launch during the construction phase is evaluated by Equation (5). The analysis results are shown in Figure 5. Each horizontal line segment represents the RAAN drift range of the corresponding piggyback launch during the construction phase, and the corresponding mission indices are labeled above the line segments. Detailed information of the piggyback launches, including the mission indices and orbital elements, can be found in [21]. The RAAN drift range of each orbit is approximately 60°, and an orbital plane at an RAAN of 100-380° can be covered by multiple piggyback launches. Thus, the initial RAAN value of the starting point Ω00 is set to 100-160°, and the orbital planes are expected to be distributed in Ω∈100°,375°. Therefore, the span of the RAAN ΔΩ is set to 120°,255°−Ω00.Figure 5
Accessible RAAN range of the piggyback launches that are inclined at angles of approximately 74° under theJ2 perturbation.


## 4. Configuration Optimization

In this section, two optimization models for the constellation configuration design are constructed, taking discrete and continuous performance indices. The genetic algorithm (GA) [17] is a stochastic optimization technique inspired by the evolution of living beings. The differential evolution algorithm (DE) [18] is a simple population-based, stochastic, and powerful function minimizer. Both algorithms are widely used global approaches to global numerical optimization. In this section, they are applied to solve the constellation design problem.


### 4.1. Optimization Model with Discrete Object

In this model, the inclination of the constellationinc is set to 74°, and the mean anomaly of the starting point M0 is set to 0. Thus, the quasi-Walker constellation configuration can be defined by Equation (8) in terms of x=P,N,a,Ω00,ΔΩ,F. Then, the performance index takes the form of the negative of Obj1 and can be calculated by Equation (1) with the predefined data of the 108 cities. The optimization model can be represented by
(13)findx=P,N,a,Ω0,ΔΩ,FminJDx=−Obj1s.t.xL<x<xU.GA [24] is one of the most widely used evolutionary algorithms. The mechanisms that govern the evolution are natural selection, reproduction, and mutation, which guarantee a robust convergence to a global numerical optimization. In this section, GA is employed to solve Equation (13) under the constraints on the design variables in Table 3. Among the design variables, P and N are integer variables, and N is set to be N≤8 because the distributer can carry at most 8 navigation satellites. The integer variables P and N are estimated by experience and then reduced gradually until there is no convergence. The parameters in GA are set as follows. The population size and maximum number of generations are set to 100 and 50, respectively. In detail, the crossover probability and mutation probability are 0.8 and 0.2, respectively.Table 3
Constraints on the design variables.
Design variablesLower boundaryUpper boundaryUnitP020—N08—a93009500KmΩ00100130DegΔΩ120320DegF0P−1—With the GA, a circular Walker constellation is obtained, which is a quasi-Walker constellation that consists of 77 satellites distributed evenly in 11 orbit planes;a=9337.001 km, inc=74°,ΔΩ=202.924°,Ω00=105.311°,F=5.914. This is a 77/11/5.914 quasi-Walker constellation. Seven orbital planes with a 20.29° separation angle in the RAAN are distributed evenly in the RAAN span of 105.31°-308.23°. With this constellation, the maximum primary index (Obj1=229) is obtained with the fewest satellites, which contributes to a system cost of Obj2=5.85. Figure 6 shows the configuration of the optimized constellation in ECI frames.Figure 6
77-satellite Walker constellation, which is the result for the optimization model with discrete object. (a) View fromazimuth=160° and elevation=45°; (b) view from azimuth=110° and elevation=90°.
(a)(b)


### 4.2. Optimization Model with Continuous Object

An improved optimization model with a continuous performance index is established here and adopted with two constellation configurations. According to Equation (1), the primary index (Obj1) is evaluated by the GDOP value at 2163 sampling points for each city. Therefore, there are 233,604 (2163×108) GDOP sampling values in total. To maximize the primary index (Obj1), the objective function in this model, which is to be minimized, can be defined as the rate of the unqualified sampling values:
(14)JCx=1−numSPGDOP≤10Nt·Ncity,where Nt=2163 is the number of sampling points of each city, Ncity=108 is the number of key cities, and numSPGDOP≤10 is the number of qualified sampling values for which the GDOP values are lower than 10.DE [25] originated from the Genetic Annealing algorithm and was initially presented to minimize continuous space functions. Later, as indicated by many studies, it exhibits better performance on many other problems than several other evolutionary algorithms [26]. Thus, for this optimization model with a continuous objective function, the DE/best/1/bin algorithm [25] is employed. In this section, the population size and maximum number of generations are set to 50 and 500, respectively. In detail, the scalar control parameter and crossover rate are 0.8 and 0.8, respectively.For a quasi-Walker constellation, the inclinationinc of the quasi-Walker constellation and mean anomaly of the starting point M0 are listed as the design variables. Thus, the constellation can be defined by Equation (8) in terms of x=P,N,a,inc,Ω00,ΔΩ,M00,F, which is the set of design variables. The constraints on the design variables are listed in Table 4.Table 4
Constraints on the design variables.
Design variablesLower boundaryUpper boundaryUnitP020—N08—a93009500kminc73.575.5degΩ00100160degΔΩ120255-Ω00degF0P−1—M000360degBy optimizing Equation (14) under the constraints in Table 4 with DE, an improved quasi-Walker constellation of 70 satellites, which are evenly distributed in 7 orbit planes, is obtained; a=9486.475km, inc=74.359°, ΔΩ=174.083°, Ω00=100.795°, F=0.004999, and M00=2.272°. Figure 7 shows the configuration of this constellation in ECI frames.Figure 7
70-satellite Walker constellation, which is the result for the optimization model with continuous object. (a) View fromazimuth=160° and elevation=45°; (b) view from azimuth=110° and elevation=90°.
(a)(b)For a restricted Flower constellation, the constellation can be defined by Equation (10) in terms of x=Np,a,inc,Ω00,ΔΩ,M00, which is the set of design variables. The constraints on the design variables are listed in Table 5.Table 5
Constraints on the design variables.
Design variablesLower boundaryUpper boundaryUnitNp6670—a93009500kminc73.575.5degΩ00100160degΔΩ120255-Ω00degM000360degBy optimizing Equation (14) under the constraints in Table 5 with DE, a restricted Flower constellation of 70 satellites is obtained; a=9487.106km, inc=74.377°, ΔΩ=219.399°, Ω00=100.554°, and M00=136.231°. Figure 8 shows the configuration of this constellation in ECI frames.Figure 8
70-satellite Flower constellation, which is the result for the optimization model with continuous object. (a) View fromazimuth=160° and elevation=45°; (b) view from azimuth=110° and elevation=90°.
(a)(b)


### 4.3. Simulation and Analysis

To analyze the navigation performance of the two constellations, the orbits are propagated on the 1st, 7th, and 30th days (i.e.,MJD2000=7396to7397, MJD2000=7402to7403, and MJD2000=7425to7426, respectively), and the time step is set to Δt=120s, which results in 2163 sampling points for each city. Then, the navigation performance for each city is evaluated. The GDOP values at all sampling points of each constellation are shown in Figures 9–11. The x scale denotes the sampling points in three days, so it is separately set to MJD2000=7396to7397, MJD2000=7402to7403, and MJD2000=7425to7426. As these figures show, all GDOP values of the three constellations are below 10, which indicates that both constellations can provide stable coverage for the target cities.Figure 9
Navigation performance of the 77-satellite Walker constellation.Figure 10
Navigation performance of the 70-satellite Walker constellation.Figure 11
Navigation performance of the 70-satellite Flower constellation.Among all 108 key cities, Qiqihar (123.95° E, 47.33° N) is located at the highest latitude, and Sanya (109.50° E, 18.25° N) is located at the lowest latitude. By taking Qiqihar, Beijing (116.40° E, 39.90° N), and Sanya as examples, Figures12–14 show the navigation performance features of the constellations for high-, middle-, and low-latitude areas, respectively. As these figures show, the GDOP values periodically exhibit variations. The 70-satellite constellations generate more GDOP values above 6 than the 77-satellite Walker constellation. Because of the high inclination, these three constellations show better navigation performance for higher-latitude areas. Although the three constellations showed different characteristics in terms of navigation performance, all of them satisfy the requirements of this regional navigation augmentation mission. Compared to the original discrete function, this new objective function, which is an approximate continuous function, effectively enlarges the convergence range. In addition, the improved optimization model generates constellations that can satisfy the navigation augmentation requirement with fewer satellites.Figure 12
Navigation performance for three typical cities of the 77-satellite Walker constellation.Figure 13
Navigation performance for three typical cities of the 70-satellite Walker constellation.Figure 14
Navigation performance for three typical cities of the 70-satellite Flower constellation.As indicated by [21, 27], the 77-satellite Walker constellation was ranked first in Problem B of the 9th China Trajectory Optimization Competition because this result can provide continuous high-accuracy navigation services for all key cities with the lowest system cost of 5.85. After the competition, the optimization model with continuous objective function was proposed, and better results have been obtained. Meanwhile, the quasi-Walker configuration was applied by Xi’an Satellite Control Center (XSCC) and compared with the Flower constellation [28]. The results obtained by XSCC [27] indicate that the optimal quasi-Walker constellation consists 75 satellites, which leads to a system cost of 5.75. In this paper, a promoted optimal quasi-Walker constellation with 70 satellites and the optimal Flower constellation are obtained by the optimization model with continuous objects. The system cost is reduced to 5.30, which verifies the effectiveness of the improved optimization model. Table 6 shows the optimization results in this paper and those obtained by XSCC.Table 6
Configuration model and optimization results in our studies and XSCC’s.
Configuration modelResults in this workResults obtained by XSCCQuasi-Walker constellationQuasi-Walker constellationRestricted Flower constellationQuasi-Walker constellationRestricted Flower constellationPerformance indexDiscreteContinuousContinuousPenalty functionSatellite7770707570P117702570N710131a (km)9337.009486.489487.119484.869487.39e00000inc (deg)74.0074.3674.3873.9074.39Ω00 (deg)105.31100.80100.5543.30152.94ΔΩ (deg)202.92174.08219.40218.69218.59M00 (deg)02.27136.23285.3612.63F5.910.005-9.35—-9.35System cost (monetary units)5.855.305.305.755.30


## 4.1. Optimization Model with Discrete Object

In this model, the inclination of the constellationinc is set to 74°, and the mean anomaly of the starting point M0 is set to 0. Thus, the quasi-Walker constellation configuration can be defined by Equation (8) in terms of x=P,N,a,Ω00,ΔΩ,F. Then, the performance index takes the form of the negative of Obj1 and can be calculated by Equation (1) with the predefined data of the 108 cities. The optimization model can be represented by
(13)findx=P,N,a,Ω0,ΔΩ,FminJDx=−Obj1s.t.xL<x<xU.GA [24] is one of the most widely used evolutionary algorithms. The mechanisms that govern the evolution are natural selection, reproduction, and mutation, which guarantee a robust convergence to a global numerical optimization. In this section, GA is employed to solve Equation (13) under the constraints on the design variables in Table 3. Among the design variables, P and N are integer variables, and N is set to be N≤8 because the distributer can carry at most 8 navigation satellites. The integer variables P and N are estimated by experience and then reduced gradually until there is no convergence. The parameters in GA are set as follows. The population size and maximum number of generations are set to 100 and 50, respectively. In detail, the crossover probability and mutation probability are 0.8 and 0.2, respectively.Table 3
Constraints on the design variables.
Design variablesLower boundaryUpper boundaryUnitP020—N08—a93009500KmΩ00100130DegΔΩ120320DegF0P−1—With the GA, a circular Walker constellation is obtained, which is a quasi-Walker constellation that consists of 77 satellites distributed evenly in 11 orbit planes;a=9337.001 km, inc=74°,ΔΩ=202.924°,Ω00=105.311°,F=5.914. This is a 77/11/5.914 quasi-Walker constellation. Seven orbital planes with a 20.29° separation angle in the RAAN are distributed evenly in the RAAN span of 105.31°-308.23°. With this constellation, the maximum primary index (Obj1=229) is obtained with the fewest satellites, which contributes to a system cost of Obj2=5.85. Figure 6 shows the configuration of the optimized constellation in ECI frames.Figure 6
77-satellite Walker constellation, which is the result for the optimization model with discrete object. (a) View fromazimuth=160° and elevation=45°; (b) view from azimuth=110° and elevation=90°.
(a)(b)


## 4.2. Optimization Model with Continuous Object

An improved optimization model with a continuous performance index is established here and adopted with two constellation configurations. According to Equation (1), the primary index (Obj1) is evaluated by the GDOP value at 2163 sampling points for each city. Therefore, there are 233,604 (2163×108) GDOP sampling values in total. To maximize the primary index (Obj1), the objective function in this model, which is to be minimized, can be defined as the rate of the unqualified sampling values:
(14)JCx=1−numSPGDOP≤10Nt·Ncity,where Nt=2163 is the number of sampling points of each city, Ncity=108 is the number of key cities, and numSPGDOP≤10 is the number of qualified sampling values for which the GDOP values are lower than 10.DE [25] originated from the Genetic Annealing algorithm and was initially presented to minimize continuous space functions. Later, as indicated by many studies, it exhibits better performance on many other problems than several other evolutionary algorithms [26]. Thus, for this optimization model with a continuous objective function, the DE/best/1/bin algorithm [25] is employed. In this section, the population size and maximum number of generations are set to 50 and 500, respectively. In detail, the scalar control parameter and crossover rate are 0.8 and 0.8, respectively.For a quasi-Walker constellation, the inclinationinc of the quasi-Walker constellation and mean anomaly of the starting point M0 are listed as the design variables. Thus, the constellation can be defined by Equation (8) in terms of x=P,N,a,inc,Ω00,ΔΩ,M00,F, which is the set of design variables. The constraints on the design variables are listed in Table 4.Table 4
Constraints on the design variables.
Design variablesLower boundaryUpper boundaryUnitP020—N08—a93009500kminc73.575.5degΩ00100160degΔΩ120255-Ω00degF0P−1—M000360degBy optimizing Equation (14) under the constraints in Table 4 with DE, an improved quasi-Walker constellation of 70 satellites, which are evenly distributed in 7 orbit planes, is obtained; a=9486.475km, inc=74.359°, ΔΩ=174.083°, Ω00=100.795°, F=0.004999, and M00=2.272°. Figure 7 shows the configuration of this constellation in ECI frames.Figure 7
70-satellite Walker constellation, which is the result for the optimization model with continuous object. (a) View fromazimuth=160° and elevation=45°; (b) view from azimuth=110° and elevation=90°.
(a)(b)For a restricted Flower constellation, the constellation can be defined by Equation (10) in terms of x=Np,a,inc,Ω00,ΔΩ,M00, which is the set of design variables. The constraints on the design variables are listed in Table 5.Table 5
Constraints on the design variables.
Design variablesLower boundaryUpper boundaryUnitNp6670—a93009500kminc73.575.5degΩ00100160degΔΩ120255-Ω00degM000360degBy optimizing Equation (14) under the constraints in Table 5 with DE, a restricted Flower constellation of 70 satellites is obtained; a=9487.106km, inc=74.377°, ΔΩ=219.399°, Ω00=100.554°, and M00=136.231°. Figure 8 shows the configuration of this constellation in ECI frames.Figure 8
70-satellite Flower constellation, which is the result for the optimization model with continuous object. (a) View fromazimuth=160° and elevation=45°; (b) view from azimuth=110° and elevation=90°.
(a)(b)


## 4.3. Simulation and Analysis

To analyze the navigation performance of the two constellations, the orbits are propagated on the 1st, 7th, and 30th days (i.e.,MJD2000=7396to7397, MJD2000=7402to7403, and MJD2000=7425to7426, respectively), and the time step is set to Δt=120s, which results in 2163 sampling points for each city. Then, the navigation performance for each city is evaluated. The GDOP values at all sampling points of each constellation are shown in Figures 9–11. The x scale denotes the sampling points in three days, so it is separately set to MJD2000=7396to7397, MJD2000=7402to7403, and MJD2000=7425to7426. As these figures show, all GDOP values of the three constellations are below 10, which indicates that both constellations can provide stable coverage for the target cities.Figure 9
Navigation performance of the 77-satellite Walker constellation.Figure 10
Navigation performance of the 70-satellite Walker constellation.Figure 11
Navigation performance of the 70-satellite Flower constellation.Among all 108 key cities, Qiqihar (123.95° E, 47.33° N) is located at the highest latitude, and Sanya (109.50° E, 18.25° N) is located at the lowest latitude. By taking Qiqihar, Beijing (116.40° E, 39.90° N), and Sanya as examples, Figures12–14 show the navigation performance features of the constellations for high-, middle-, and low-latitude areas, respectively. As these figures show, the GDOP values periodically exhibit variations. The 70-satellite constellations generate more GDOP values above 6 than the 77-satellite Walker constellation. Because of the high inclination, these three constellations show better navigation performance for higher-latitude areas. Although the three constellations showed different characteristics in terms of navigation performance, all of them satisfy the requirements of this regional navigation augmentation mission. Compared to the original discrete function, this new objective function, which is an approximate continuous function, effectively enlarges the convergence range. In addition, the improved optimization model generates constellations that can satisfy the navigation augmentation requirement with fewer satellites.Figure 12
Navigation performance for three typical cities of the 77-satellite Walker constellation.Figure 13
Navigation performance for three typical cities of the 70-satellite Walker constellation.Figure 14
Navigation performance for three typical cities of the 70-satellite Flower constellation.As indicated by [21, 27], the 77-satellite Walker constellation was ranked first in Problem B of the 9th China Trajectory Optimization Competition because this result can provide continuous high-accuracy navigation services for all key cities with the lowest system cost of 5.85. After the competition, the optimization model with continuous objective function was proposed, and better results have been obtained. Meanwhile, the quasi-Walker configuration was applied by Xi’an Satellite Control Center (XSCC) and compared with the Flower constellation [28]. The results obtained by XSCC [27] indicate that the optimal quasi-Walker constellation consists 75 satellites, which leads to a system cost of 5.75. In this paper, a promoted optimal quasi-Walker constellation with 70 satellites and the optimal Flower constellation are obtained by the optimization model with continuous objects. The system cost is reduced to 5.30, which verifies the effectiveness of the improved optimization model. Table 6 shows the optimization results in this paper and those obtained by XSCC.Table 6
Configuration model and optimization results in our studies and XSCC’s.
Configuration modelResults in this workResults obtained by XSCCQuasi-Walker constellationQuasi-Walker constellationRestricted Flower constellationQuasi-Walker constellationRestricted Flower constellationPerformance indexDiscreteContinuousContinuousPenalty functionSatellite7770707570P117702570N710131a (km)9337.009486.489487.119484.869487.39e00000inc (deg)74.0074.3674.3873.9074.39Ω00 (deg)105.31100.80100.5543.30152.94ΔΩ (deg)202.92174.08219.40218.69218.59M00 (deg)02.27136.23285.3612.63F5.910.005-9.35—-9.35System cost (monetary units)5.855.305.305.755.30


## 5. Constellation Deployment

In this section, a synthetic procedure to insert each satellite is proposed by applying the Hohmann transfer and the Lambert transfer and utilizing the natural drift underJ2 perturbation. A cross-plane deployment strategy is proposed, which is a creative method for economical constellation construction.


### 5.1. Cross-Plane Constellation Deployment Strategy

Figure15 shows the correspondence between the piggyback launches and target orbit planes. The vertical dashed lines represent the RAAN of the orbit planes at the initial time of the service phase in sequence, and each horizontal line segment represents the RAAN drift range of the corresponding piggyback launch during the construction phase. In cases where the vertical and horizontal lines cross, satellites in the piggyback launch can transfer to the target plane. If the RAAN drift range of a piggyback launch intersects with more than one orbit plane, the 8 satellites in that launch can be transferred to different planes. Ten piggyback launches are chosen to deploy a constellation consisting of 11 orbit planes. The piggyback launch indices are 6, 44, 33, 13, 8, 92, 61, 47, 10, and 82, and some groups of 8 satellites in one piggyback launch transfer to different planes. The deployment scheme is shown in Figure 16. As this figure illustrates, the 56 satellites of missions 33, 13, 8, 92, 61, 47, and 10 are deployed into 8 planes, and therefore, 1 piggyback launch is saved by this strategy. This cross-plane deployment strategy efficiently reduces the cost of constellation construction.Figure 15
Correspondence of piggyback launches and target orbital planes.Figure 16
The cross-plane deployment strategy.


### 5.2. Five-Impulse Orbit Transfer

J2 perturbation is expected to create a relative precession motion between two orbital planes to obtain the desired RAAN separation. Thus, a small tangential impulse Δv0 (<4 m/s) at the start of the construction phase (t=7305) is applied to adjust the initial state of the satellite in order to achieve the proper phase after a period of natural drift. Then, the satellites are maneuvered to the target height by two burns at the perigee and apogee of the transfer orbit. These two burns contribute to significant changes to the orbital height as well as moderate changes to the inclination and RAAN. Following these two burns, the satellites are propagated to the end of the construction phase. This drift period is also employed to allow the satellites to drift to the target phase. Finally, the satellites are precisely inserted through two burns by the Lambert transfer under J2 perturbation. As above, a five-impulse transfer procedure is applied to each satellite to achieve its target position. The entire procedure is performed as described below.In the first optimization model, grid search is used to findT0, which is the interval between the start of the deployment and the second impulse Δvp, for the minimal error of the RAAN at the insertion point. The initial values of Δvp and Δva are calculated by the Hohmann transfer [29, 30], and the initial value of Δv0 is set to 0. Then, taking the weighted sum of the error of the orbit radius, error of eccentricity, and error of the RAAN at the insertion point as the objective function, we optimize T0, the second and third impulses Δvp and Δva, and the small velocity increment at the start of the construction phase Δv0 by using the sequential quadratic programming algorithm. Notably, Δvp is expected to be performed at the perigee for a lower fuel consumption, but the optimization and grid search result for T0 may not be exactly the perigee. Thus, the orbit is propagated to the perigee after T0; then, Δvp is performed. The third impulse is performed at the apogee after Δvp. The optimization model is described as follows:
(15)findx=Δv0,T0,Δvp,ΔvaminJx=∆a∙w1+∆e∙w2+∆Ω∙w3s.t.xL<x<xU,where x=Δv0,T0,Δvp,Δvais the optimization variable; ∆a,∆e,and∆Ω are the error of the orbit radius, error of eccentricity, and error of RAAN at the insertion point, respectively; w1,w2,w3 are the corresponding weights; and xL=x−1.5×10−6km/s,3MJD2000,0.03km/s,0.03km/s and xU=x+1.5×10−6km/s,3MJD2000,0.03km/s,0.03km/s are the lower bound and upper bound of the optimization variable, respectively.The second optimization intends to minimize the error of the argument of latitude at the insertion pointΔu by optimizing the velocity increment at the start of the construction phase Δv0 in a small interval Δv0−0.0015km/s,Δv0−0.0015km/s using the sequential quadratic programming algorithm. Since the amplitude of latitude rapidly changes, the small change at the starting point causes a large difference in amplitude of latitude after a long propagation, while the RAAN and orbital height are less sensitive to small changes at the starting point. Thus, adding Δu in the objective function of the first optimization causes difficulties in convergence. By reoptimizing Δv0, Δu is significantly reduced, while ∆a,∆e, and ∆Ω do not greatly change. Therefore, optimizing Δv0 twice will facilitate the convergence and optimality of all parameters.Finally, the satellites inject att=7396 through the last two impulses ΔvL1,ΔvL2 calculated based on the Lambert transfer under the J2 perturbation. If the total fuel consumption is greater than 20 kg, which is the predefined fuel mass of the competition, the weights w1,w2,w3 should be modified until convergence. The sequence of the five-impulse transfer and the target errors of each impulse are shown in Table 7. Figure 17 shows the design process of the transfer of each satellite.Table 7
Velocity increments of the five-impulse transfer.
Impulse indexSymbolTimeTarget error1Δv0t=7305∆Ω,Δu2ΔvpPerigee afterT0∆a,∆e3ΔvaApogee afterT04,5ΔvL1,ΔvL2t=7396∆a,∆e,Δinc,∆Ω,Δu,ΔtI (ΔtI is the error of the insertion time)Figure 17
Flowchart of the transfer design process.Taking satellite #15 as an example, which is the first satellite in the third orbital plane, the five-impulse transfer procedure is presented in Table8. Satellite #15 departs from piggyback launch #33, and the total velocity increment is 0.632 km/s. The transfer process is shown in Figure 18, where the blue and red orbits are the piggyback launch and target orbit, respectively, the impulses are labeled by black dots, the arrows denote the directions of the velocity increments, and the grey parts are the natural drift periods of the satellite.Table 8
Transfer procedure of satellite #15.
Impulse indexTime (MJD2000)Velocity increment (km/s)Fuel consumption (kg)173050.0008160331332460.02773356074727377.3015283261500.3087904966734769.96114328537937377.3475928005710.2872459349474288.37296250756147395.9578833015830.0301758270006030.833076919784573960.0052407562472050.143814923541Figure 18
Transfer process of satellite #15.


## 5.1. Cross-Plane Constellation Deployment Strategy

Figure15 shows the correspondence between the piggyback launches and target orbit planes. The vertical dashed lines represent the RAAN of the orbit planes at the initial time of the service phase in sequence, and each horizontal line segment represents the RAAN drift range of the corresponding piggyback launch during the construction phase. In cases where the vertical and horizontal lines cross, satellites in the piggyback launch can transfer to the target plane. If the RAAN drift range of a piggyback launch intersects with more than one orbit plane, the 8 satellites in that launch can be transferred to different planes. Ten piggyback launches are chosen to deploy a constellation consisting of 11 orbit planes. The piggyback launch indices are 6, 44, 33, 13, 8, 92, 61, 47, 10, and 82, and some groups of 8 satellites in one piggyback launch transfer to different planes. The deployment scheme is shown in Figure 16. As this figure illustrates, the 56 satellites of missions 33, 13, 8, 92, 61, 47, and 10 are deployed into 8 planes, and therefore, 1 piggyback launch is saved by this strategy. This cross-plane deployment strategy efficiently reduces the cost of constellation construction.Figure 15
Correspondence of piggyback launches and target orbital planes.Figure 16
The cross-plane deployment strategy.


## 5.2. Five-Impulse Orbit Transfer

J2 perturbation is expected to create a relative precession motion between two orbital planes to obtain the desired RAAN separation. Thus, a small tangential impulse Δv0 (<4 m/s) at the start of the construction phase (t=7305) is applied to adjust the initial state of the satellite in order to achieve the proper phase after a period of natural drift. Then, the satellites are maneuvered to the target height by two burns at the perigee and apogee of the transfer orbit. These two burns contribute to significant changes to the orbital height as well as moderate changes to the inclination and RAAN. Following these two burns, the satellites are propagated to the end of the construction phase. This drift period is also employed to allow the satellites to drift to the target phase. Finally, the satellites are precisely inserted through two burns by the Lambert transfer under J2 perturbation. As above, a five-impulse transfer procedure is applied to each satellite to achieve its target position. The entire procedure is performed as described below.In the first optimization model, grid search is used to findT0, which is the interval between the start of the deployment and the second impulse Δvp, for the minimal error of the RAAN at the insertion point. The initial values of Δvp and Δva are calculated by the Hohmann transfer [29, 30], and the initial value of Δv0 is set to 0. Then, taking the weighted sum of the error of the orbit radius, error of eccentricity, and error of the RAAN at the insertion point as the objective function, we optimize T0, the second and third impulses Δvp and Δva, and the small velocity increment at the start of the construction phase Δv0 by using the sequential quadratic programming algorithm. Notably, Δvp is expected to be performed at the perigee for a lower fuel consumption, but the optimization and grid search result for T0 may not be exactly the perigee. Thus, the orbit is propagated to the perigee after T0; then, Δvp is performed. The third impulse is performed at the apogee after Δvp. The optimization model is described as follows:
(15)findx=Δv0,T0,Δvp,ΔvaminJx=∆a∙w1+∆e∙w2+∆Ω∙w3s.t.xL<x<xU,where x=Δv0,T0,Δvp,Δvais the optimization variable; ∆a,∆e,and∆Ω are the error of the orbit radius, error of eccentricity, and error of RAAN at the insertion point, respectively; w1,w2,w3 are the corresponding weights; and xL=x−1.5×10−6km/s,3MJD2000,0.03km/s,0.03km/s and xU=x+1.5×10−6km/s,3MJD2000,0.03km/s,0.03km/s are the lower bound and upper bound of the optimization variable, respectively.The second optimization intends to minimize the error of the argument of latitude at the insertion pointΔu by optimizing the velocity increment at the start of the construction phase Δv0 in a small interval Δv0−0.0015km/s,Δv0−0.0015km/s using the sequential quadratic programming algorithm. Since the amplitude of latitude rapidly changes, the small change at the starting point causes a large difference in amplitude of latitude after a long propagation, while the RAAN and orbital height are less sensitive to small changes at the starting point. Thus, adding Δu in the objective function of the first optimization causes difficulties in convergence. By reoptimizing Δv0, Δu is significantly reduced, while ∆a,∆e, and ∆Ω do not greatly change. Therefore, optimizing Δv0 twice will facilitate the convergence and optimality of all parameters.Finally, the satellites inject att=7396 through the last two impulses ΔvL1,ΔvL2 calculated based on the Lambert transfer under the J2 perturbation. If the total fuel consumption is greater than 20 kg, which is the predefined fuel mass of the competition, the weights w1,w2,w3 should be modified until convergence. The sequence of the five-impulse transfer and the target errors of each impulse are shown in Table 7. Figure 17 shows the design process of the transfer of each satellite.Table 7
Velocity increments of the five-impulse transfer.
Impulse indexSymbolTimeTarget error1Δv0t=7305∆Ω,Δu2ΔvpPerigee afterT0∆a,∆e3ΔvaApogee afterT04,5ΔvL1,ΔvL2t=7396∆a,∆e,Δinc,∆Ω,Δu,ΔtI (ΔtI is the error of the insertion time)Figure 17
Flowchart of the transfer design process.Taking satellite #15 as an example, which is the first satellite in the third orbital plane, the five-impulse transfer procedure is presented in Table8. Satellite #15 departs from piggyback launch #33, and the total velocity increment is 0.632 km/s. The transfer process is shown in Figure 18, where the blue and red orbits are the piggyback launch and target orbit, respectively, the impulses are labeled by black dots, the arrows denote the directions of the velocity increments, and the grey parts are the natural drift periods of the satellite.Table 8
Transfer procedure of satellite #15.
Impulse indexTime (MJD2000)Velocity increment (km/s)Fuel consumption (kg)173050.0008160331332460.02773356074727377.3015283261500.3087904966734769.96114328537937377.3475928005710.2872459349474288.37296250756147395.9578833015830.0301758270006030.833076919784573960.0052407562472050.143814923541Figure 18
Transfer process of satellite #15.


## 6. Conclusions

This work presents a synthetic method of designing a constellation for regional coverage under a constrained deployment approach. In summary, this work has three steps. First, the search space is reduced by estimating the coverage region for the target cities and the reachable position of the potential piggyback launches. Second, two optimization models are constructed for constellation configuration optimization. In particular, the continuous optimization model contributes to numerical convergence to the optimal solution with fewer satellites. Third, the five-impulse transfer and cross-plane deployment strategies make a significant contribution to economical constellation construction. The solution was ranked first in Problem B of the 9th China Trajectory Optimization Competition.



---

*Source: 5042650-2021-09-06.xml*