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

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 under J2 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.


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 from http://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][2][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][12][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 J 2 -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.

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. Figure 1 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 = 7396 to 7397, MJD2000 = 7402 to 7403, and MJD2000 = 7425 to 7426) with a fixed interval Δt = 120 s.
The primary index Obj 1 is aimed at maximizing the total navigation revenue of all key cities: where δ k is a 0-1 variable, which is used to measure whether the kth city can meet the navigation accuracy requirements, and w k is the weight of the kth city. GDOP k,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 2 International Journal of Aerospace Engineering sampling points are below 10, Obj 1 is equal to 229, which is the maximum value of Obj 1 . The secondary index Obj 2 is aimed at minimizing the entire cost of constellation construction. The construction phase is from MJD2000 = 7305 to 7396. During this phase, the navigation constellation can be constructed by new launches or piggyback launches.
where N launch , N carry , and N naviSat are the number of new launches, the number of piggyback launches, and the number of navigation satellites in the constellation, respectively. C launch = 1:2, C carry = 0:2, and C naviSat = 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 the J 2 perturbation are considered. The dynamic models are described as follows: where μ is the gravitational constant of the Earth, x, y, z are the position components of the satellite in the Earthcentered inertial (ECI) coordinate system, r = ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi x 2 + y 2 + z 2 p is the geocentric distance, R e is the equatorial radius, and J 2 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 under J 2 perturbation are described as follows: 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 A n = ffiffiffiffiffiffiffiffi ffi μ/a 3 p is the mean motion of the orbit. When Ω 0 , ω 0 , and M 0 are given at t 0 , the elements at time t can be calculated by During the construction phase, the satellites are propelled by a chemical engine, and the residual mass m f after an impulse is given by where m i is the satellite mass before the impulse, Δv is the magnitude of the velocity impulse, I sp is the specific impulse, and g 0 is the gravitational acceleration at sea level. All the parameters are set according to [21].

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: 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. Ω ij ð0Þ and M ij ð0Þ 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 M 0 , 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 where Sat ij ðt 0 Þ, Ω ij ðt 0 Þ, and M ij ðt 0 Þ 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 = t 0 , i = 1, 2, ⋯, P, j = 1, 2, ⋯, N. In this work, t 0 = 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.

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 N p , number of sidereal days to repeat the ground track N d , number of satellites N s , phasing parameters Ω 00 , ΔΩ, M 00 , 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 N p revolutions over N d days, the following relationship is required where T Ω and T ΩG are the nodal period of the orbit and nodal period of Greenwich, respectively. In this paper, N d is set to 1 according to the mission requirements. The orbits in a Flower constellation can be defined by the following equations [7]: International Journal of Aerospace Engineering where Sat i ðt 0 Þ, Ω i ðt 0 Þ, and M i ðt 0 Þ are the Kepler orbit elements, RAAN, and mean anomaly of the ith satellite at the initial time t = t 0 , i = 1, 2, ⋯, N p ; a, inc, Ω 00 , ΔΩ, and M 00 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 = ffiffiffiffiffiffiffiffi ffi μ/a 3 p is the mean motion; ω E is the Earth spin rate.

Problem Analysis
The deployment of a constellation with constrained piggyback launches brings additional complexity to the constella-tion 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 Figure 2. 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 3 shows how the primary index Obj 1 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, M 00 = 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 Obj 1 value increases with increasing orbit   5 International Journal of Aerospace Engineering radius, and the Obj 1 value decreases rapidly with an increase in the inclination when it is larger than the maximum inclination for the best Obj 1 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 Obj 1 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.

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 Table 2. 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.

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: where ϵ denotes the minimum elevation angle, which is set to 10°in this problem. Ulybyshev defined a twodimensional 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°. 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 with J 2 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 Þ.

Configuration Optimization
In this section, two optimization models for the constellation configuration design are constructed, taking discrete and   [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.

Optimization Model with Discrete
Object. In this model, the inclination of the constellation inc is set to 74°, and the mean anomaly of the starting point M 0 is set to 0. Thus, the quasi-Walker constellation configuration can be defined by Equation (8) 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.
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 (Obj 1 = 229) is obtained with the fewest satellites, which contributes to a system cost of Obj 2 = 5:85. Figure 6 shows the configuration of the optimized constellation in ECI frames.

Optimization Model with Continuous
Object. An improved optimization model with a continuous performance index is established here and adopted with two where N t = 2163 is the number of sampling points of each city, N city = 108 is the number of key cities, and numðS

10
International Journal of Aerospace Engineering P GDOP ≤ 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 inclination inc of the quasi-Walker constellation and mean anomaly of the starting point M 0 are listed as the design variables. Thus, the constellation can be defined by Equation (8) in terms of x = fP, N, a, inc, Ω 00 , ΔΩ, M 00 , Fg, which is the set of design variables. The constraints on the design variables are listed in Table 4.
For a restricted Flower constellation, the constellation can be defined by Equation (10) in terms of x = fN p , a, inc, Ω 00 , ΔΩ, M 00 g, which is the set of design variables. The constraints on the design variables are listed in Table 5.

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 = 7396 to 7397,

12
International Journal of Aerospace Engineering MJD2000 = 7402 to 7403, and MJD2000 = 7425 to 7426, respectively), and the time step is set to Δt = 120 s, 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 = 7396 to 7397, MJD2000 = 7402 to 7403, and MJD2000 = 7425 to 7426. 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. 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, Figures 12-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.
As indicated by [21,27], the 77-satellite Walker constellation was ranked first in Problem B of the 9th China 13 International Journal of Aerospace Engineering 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.

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 under J 2 perturbation. A cross-plane deployment strategy is proposed, which is a creative method for economical constellation construction.    Figure 15 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.

5.2.
Five-Impulse Orbit Transfer. J 2 perturbation is expected to create a relative precession motion between two orbital planes to obtain the desired RAAN separation. Thus, a small tangential impulse Δv 0 (<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 J 2 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 find T 0 , which is the interval between the start of the deployment and the second impulse Δv p , for the minimal error of the RAAN at the insertion point. The initial values of Δv p and Δv a are calculated by the Hohmann transfer [29,30], and the initial value of Δv 0 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 T 0 , the second and third impulses Δv p and Δv a , and the small velocity increment at the start of the construction phase Δv 0 by using the sequential quadratic programming algorithm. Notably, Δv p is expected to be performed at the perigee for a lower fuel consumption, but the optimization and grid search result for T 0 may not be exactly the perigee. Thus, the orbit is propagated to the perigee after T 0 ; then, Δv p is performed. The third impulse is performed at the apogee after Δv p . The optimization model is described as follows: where x = fΔv 0 , T 0 , Δv p , Δv a g is 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; w 1 , w 2 , w 3 are the corresponding weights; and x L = x − ½1:5 × ð10 −6 km/sÞ, 3 MJD2000, 0:03 km/s, 0:03 km/s and x U = x + ½1:5 × ð10 −6 km/s, 3 MJD2000, 0:03 km/s, 0:03 km/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 Δv 0 in a small interval ðΔv 0 − 0:0015 km/s, Δv 0 − 0:0015 km/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 Δv 0 , Δu is significantly reduced, while Δa, Δe, and ΔΩ do not greatly change. Therefore, optimizing Δv 0 twice will facilitate the convergence and optimality of all parameters.
Finally, the satellites inject at t = 7396 through the last two impulses Δv L1 , Δv L2 calculated based on the Lambert transfer under the J 2 perturbation. If the total fuel consumption is greater than 20 kg, which is the predefined fuel mass of the competition, the weights w 1 , w 2 , w 3 should be modified until convergence. The sequence of the five-impulse transfer and the target errors of each impulse are shown in Piggy back mission index Target orbital plane Figure 16: The cross-plane deployment strategy. 16 International Journal of Aerospace Engineering Table 7. Figure 17 shows the design process of the transfer of each satellite. Taking satellite #15 as an example, which is the first satellite in the third orbital plane, the five-impulse transfer procedure is presented in Table 8. 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.

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.

Data Availability
The background material data can be provided if necessary.  Figure 17: Flowchart of the transfer design process.