Satellite Constellation Orbit Design Optimization with Combined Genetic Algorithm and Semianalytical Approach

Abstract

This paper focuses on maximizing the percent coverage and minimizing the revisit time for a small satellite constellation with limited coverage. A target area represented by a polygon defined by grid points is chosen instead of using a target point only. The constellation consists of nonsymmetric and circular Low Earth Orbit (LEO) satellites. A global optimization method, Genetic Algorithm (GA), is chosen due to its ability to locate a global optimum solution for nonlinear multiobjective problems. From six orbital elements, five elements (semimajor axis, inclination, argument of perigee, longitude of ascending node, and mean anomaly) are varied as optimization design variables. A multiobjective optimization study is conducted in this study with percent coverage and revisit time as the two main parameters to analyze the performance of the constellation. Some efforts are made to improve the objective function and to minimize the computational load. A semianalytical approach is implemented to speed up the guessing of initial orbital elements. To determine the best parametric operator combinations, the fitness value and the computational time from each study cases are compared.

[ProQuest: [...] denotes non US-ASCII text; see PDF]

Tania Savitri 1 and Youngjoo Kim 1 and Sujang Jo 2 and Hyochoong Bang 1

Academic Editor:Christian Circi

1, Department of Aerospace Engineering, Korea Advanced Institute of Science and Technology (KAIST), Daejeon 305-701, Republic of Korea

2, Korea Aerospace Research Institute (KARI), Daejeon 305-701, Republic of Korea

Received 27 February 2017; Revised 26 April 2017; Accepted 27 April 2017; 31 May 2017

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Compared to a single satellite, a constellation may provide better coverage and higher reliability under failure of some satellites, ensuring a higher rate of survival and mission success. A constellation can offer unique capabilities that are often difficult to achieve through different means, for instance, enhanced temporal coverage [1]. A satellite orbit constellation design usually relies on the Walker approach [2] or streets of coverage method [3, 4]. Besides, a ground track-based approach has been developed [5], and recently the sliding ground track concept applied to constellations composed of one or more orbital planes has been introduced [6].

Related to the coverage mission, the idea of satellite constellation designs for complex coverage was presented by Ulybyshev [7]. Indrikis and Cleave, in their work about SPACE EGGS or the Satellite Coverage Model for LEO (Low Earth Orbit) Constellations [8], also tried to access the effectiveness of global, regional, and area coverage for proliferated small satellite constellations in low altitude orbits, stressing the capability of conventional analytical techniques. The idea of staged deployments of satellite constellations in LEO was also proposed by including the uncertainty feature of the expected number of users and activity level [9]. In the process, a satellite's life cycle cost and its capacity is traded off to satisfy the minimum per-channel performance requirement. Another approach to optimize a satellite constellation design by using tiers of satellites with variations in the orbit's altitude and inclination parameters was proposed Razoumny et al. in [10]. This approach tries to reduce the redundancy in the Walker approach by dividing the regional coverage into several latitude regions that can be addressed by different pairs of satellites' altitude and inclination, allowing for a better revisit time or a reduced number of satellites.

Besides the analytic methods, the genetic algorithm (GA) is known for its robustness in obtaining a global optimum solution for nonlinear multivariable problems through its stochastic and heuristic search algorithms. In their thesis, Pegher and Parish [11] tried to compare the coverage optimization and the revisit time of sparse military satellite constellations using traditional approaches and GA. Another method of hybrid satellite constellation design called Genetic Satellite Constellation (GSC) was proposed in [12] by using single GA optimization. This approach demonstrates the performance of hybrid LEO/MEO, LEO/GEO, and MEO/GEO constellations. A study on Simulated Annealing and GA approaches to satellite constellation design for coverage of a limited latitude region was conducted by Crossley and Williams [13], where both methods outperformed the conventional Walker approach at low Earth central angles.

Multiobjective genetic algorithm (MOGA) is a variation of GA, which is known to be a robust technique to solve multiobjective optimization, resulting in a Pareto optimal set solution [14]. A class of fast and elite MOGA was introduced by Deb et al. in 2002 [15]. This new MOGA called NSGA-II stands for Nondominated Sorting GA, which uses a fast nondominated sorting and a crowding distance assignment in its algorithm. Several multiobjective optimization studies have been performed. For example, Ely used two-branch tournament GA for a constellation design with eccentric orbits, aiming to minimize the maximum constellation's orbit altitude and the total number of satellites [16]. Several works from Ferringer are also focused on multiobjective optimization of satellite constellation designs including the trade-off analysis [17-20]. Related to MOGA, Mason et al. introduced a type of MOGA approach, called the Modified Illinois Nondominated Sorting Genetic Algorithm (MINSGA), combined with Satellite Tool Kit (STK) software to produce several constellation designs that provides a continuous global coverage [21]. Another work was also done by Confessore et al. to optimize a satellite constellation's orbit design with GA aiming to minimize the number of satellites and maximize region coverage [22].

However, to the best of our knowledge, there is no specific unified design approach for a local continuous coverage or surveillance mission over a region. This paper addresses coverage and revisit time problems using a sparse satellite constellation design with limited coverage capability, as is often encountered in low-cost satellite missions. Some low-cost small satellites tend to use small cameras for image acquisition [23, 24]. Several recent developments made in this area are the BRITE mission, which is an explorer targeting bright stars, conducted by the Canadian Space Agency, partnering with the University of Vienna and Graz University of Technology from Austria and the Copernicus Astronomical Center from Poland [16], and the PRISM Earth-imaging validation mission developed by the University of Tokyo using a CMOS imager on a narrow-angle camera [25].

The contribution of this paper is twofold: the first is a method to address a complex coverage area by using sparse satellite constellation design with a limited coverage capability and the second is an efficient approach to lower the computational burden when performing the GA optimization. In this paper, a case study is conducted on a sparse satellite constellation using multiple circular LEO satellites equipped with imaging sensors. A GA-based optimization is chosen due to its ability to locate a global optimum solution for a wide class of nonlinear problems that may be applicable to the case of this study. Several constellations' figures of merits [26], such as the area percent coverage and revisit time, are selected to evaluate the constellation performance using a multiobjective optimization algorithm. As much as GA seems to be a compatible optimization method to solve such a problem, it is often avoided due to its high computational load [27]. To minimize the computational load, a semianalytical approach is also proposed, which brings about a significant savings in computational time. The proposed approach, a combination of GA and an efficient semianalytical approach, will be a viable option for optimum satellite constellation designs.

2. Problem Statement

2.1. Mathematical Formulation of a Target Area

The area of interest is defined by inputting latitude-longitude coordinates of the main cities located in the country of interest. From these points, a convex hull polygon is generated using a computational geometry method, as depicted in Figure 1. The convex hull of a set of points is the smallest convex set that contains the points [28]. The area of interest is divided into several effective grid points, which are then evaluated to determine the target access time.

Figure 1: Example region (a) and the polygon made by convex hull mapping (b).

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

The grid points are defined in geodetic latitude _Φt and longitude _λt , where t=1,...,N and N is the total number of grid points. Furthermore, zero ellipsoidal height of the target points is assumed.

2.2. Orbit Design

Initial orbit design consists of six standard orbital elements: semimajor axis a, eccentricity e, inclination i, longitude of ascending node ohm, argument of perigee ω, and mean anomaly M.

The maximum Earth-centered half-angle θ and the maximum swath width of the satellite sensor occur when the satellite is at apogee _ra =a(1+e). As shown in Figure 2, a satellite footprint projection on the Earth's surface with elevation angle [straight epsilon] defines a coverage circle with a radius θ through the following equation [29]: [figure omitted; refer to PDF]

Figure 2: Satellite footprint projection.

[figure omitted; refer to PDF]

3. Semianalytical Initial Guess

Using a semianalytical initial guess, the time step for propagation could be adjusted in order to reduce the overall computational load. Comparison between this method and the numerical propagation method in terms of the area coverage difference is presented as well in this paper. This method is based upon the correlation between nonsingular orbital elements and the target Φ and Ψ. This approach is especially useful for observation over a repeating ground track orbit.

In comparison with numerical propagation, this approach allows for information about the satellite position at a specific time without having to propagate the satellite's orbital elements or the satellite's position and velocity vector over a time step. The accuracy of numerical propagation results depends on the size of the time step, which, on the other hand, conflicts with the size of available computing memory and computational load.

By the mean element theory [30], a satellite orbit is described using its mean orbital elements. The time rate of secular changes of mean Ω, ω, and M is given by [31] [figure omitted; refer to PDF] where K=(3/2)_J2^Re2 μ and the mean motion n is given by [figure omitted; refer to PDF] Drag perturbation is described using the following equations [32]: [figure omitted; refer to PDF]

3.1. Latitude Access

The correlation between the satellite mean orbital elements and the access to a target region's latitude is defined in the following form [32]: [figure omitted; refer to PDF] where _[varphi]s represents the angular distance along the ground track measured from the ascending node at time = 0. The equation is only applied to an orbit with small eccentricity. On the assumption that the satellite in this study is equipped with a conical sensor, (6) becomes [figure omitted; refer to PDF] For a noncircular orbit, the maximum swath width _θmax achieved on the perigee point is given by [figure omitted; refer to PDF] It is recommended to use the time derivative of θ for a noncircular orbit with drag and _J2 perturbation included. [figure omitted; refer to PDF]

The latitude access time repeats for every period of mean argument of latitude or every 2π/(ω +M ). _{Tlat i} is defined as a set of latitude access times such that [figure omitted; refer to PDF]

3.2. Longitude Access

The correlation between satellite orbital elements and the area longitude is defined as [figure omitted; refer to PDF] where l is the mean argument of latitude and - or +_θmax represents the entrance and exit of the access. A similar approach using a straightforward expansion to solve this equation is explained in [29] with a difference in the usage of atan[...]2 instead of the usual atan, since the mean argument of latitude lies between -π and π. The analytical solution approach is described briefly by first defining several variables. Let _A1 , _A2 , and _A3 be defined as follows: [figure omitted; refer to PDF] Substituting (13)-(15) into (12) results in the following equation: [figure omitted; refer to PDF]

Equation (16) can be solved by using a solver in MATLAB®. α=atan2(y,x) is a quadrant-sensitive inverse of tan[...]α=y/x, which returns a value in the range -π<=α<=π. Using Taylor series expansion, atan[...](_A3 sin[...]l,cos[...]l) can be defined as [figure omitted; refer to PDF] Let us put a value of l into the equivalent Kepler equation in terms of the mean element and derive Δt: [figure omitted; refer to PDF] _{Tlon j} is defined as a set of latitude access time for which [figure omitted; refer to PDF]

The time for intersection of a target point can be found from the time window intersection of each latitude and longitude access. The intersection time _Taccess consists of s access times, where [figure omitted; refer to PDF] The flowcharts of the algorithm are depicted in Figures 3 and 4.

Figure 3: Latitude and longitude target access determination flowchart.

[figure omitted; refer to PDF]

Figure 4: Target access time intersection flowchart.

[figure omitted; refer to PDF]

4. Multiobjective Optimization

4.1. Pareto Approach

Multiobjective optimizations are usually expressed by introducing different objectives in a cost function, for which each objective importance is represented by a weight _wi parameter for the i-th objective, as in [33]. [figure omitted; refer to PDF] where [figure omitted; refer to PDF]

Using such a weighted sum approach, the objective functions should always be normalized or scaled so that their objective values are in similar magnitudes. However, according to [34], there is a certain drawback of minimizing the weighted sums of objectives for Pareto set generation in multicriteria optimization problems. The reason is that an evenly distributed set of weighting vectors cannot guarantee an evenly distributed representation of the Pareto optimal solutions.

4.2. Nondominated Sorting Genetic Algorithms Optimization (NSGA-II)

As described in Figure 5, NSGA-II, which stands for Nondominated Sorting GA, uses a fast nondominated sorting and a crowding distance assignment [8]. GA optimization itself is a heuristic, stochastic global search based on the Darwinian concepts of evolution and natural selection and was first introduced by John Holland in mid-1960s [33]. An elitist GA favors individuals with a better fitness measure (rank). A controlled elitist GA, which is used in this algorithm, facilitates the optimization process to converge into a set of diverse solutions, even if the solution might consist of several lower performing individuals. NSGA-II is selected for this research because it has been shown to exhibit a better spread and convergence, relatively close to the true Pareto optimal front compared to PAES and SPEA [34]. SPEA stands for Strength Pareto Evolutionary Algorithm, and PAES is Pareto Archived Evolution Strategy algorithm. PAES is a multiobjective optimizer approach, which uses an archive of previously found solutions in order to identify the dominance ranking of the current and candidate solution vectors. SPEA provides a form of elitism by using an archive of the nondominated set, which is maintained separately from the population of candidate solutions.

Figure 5: NSGA-II schema.

[figure omitted; refer to PDF]

In the first generation, parent population _Pt of I individuals creates an offspring population _Qt also in I individuals through the common tournament Selection, Crossover, and Mutation operators. Next, population _Rt with a size of 2I individuals goes through a nondominated sorting. After the first generation, crowded tournament selection is used to select the top I individuals from _Rt to form the next generation _Pt+1 . A solution that wins a tournament has either the highest rank or a better crowding distance parameter. The crowding distance parameter itself is defined as the largest cuboid surrounding a solution, for which no other solution is present [35].

Two parameters used in this algorithm to control the elitism are the Pareto fraction, which is the nondominated sorting using the Pareto set, and the distance function, which uses a crowding distance assignment. The Pareto fraction limits the number of individuals that can be put into the Pareto front (elite members). The distance function itself favors individuals who are located relatively far from the front since it can create a diverse, noncrowded set of solutions.

In this paper, the parametric operator is varied for the number of generations and population. The number of generations is usually set at the beginning of the optimization and is usually used as a stopping criterion for the optimization. Another useful parameter is the stall generation, where the optimization stops if the best fitness function value does not stop after hitting several generations until the stall generation is met. Population, which has a constant value during the optimization, is the number of individual sets for every generation.

5. Simulations

A target region for the case South Korea, with its 37 cities' latitude and longitude positions, is used as an input for simulation. The area of interest itself is represented as a polygon of boundary points that includes all of the 37 cities. This area of the polygon is then divided into a number of grid points. N denotes the number of fast nondominated target points which is the number of grid points located inside the polygon, for which 10×10 grids dividing the example area result in a total of 50 grid points. The basic design criteria and constraints are as follows. A maximum satellite constellation of three satellites with a 10-degree half-angle of a conical sensor is used. This conical sensor is set as well so that only one satellite cannot cover the whole area of interest. Although from (1) it was mentioned that the maximum swath width occurs at the apogee, due to the circular orbit used in this simulation, this parameter is a matter of the satellite altitude only, assuming that the pointing error is neglected.

The constellation is placed into a circular orbit at an epoch date on August 1, 2020. Up to six hours of the satellite's lifetime is studied. In comparison with the simulation done in [29], the satellite's orbiting time period is much smaller; therefore, the argument of perigee and mean anomaly are included as optimization variables to compensate for the shorter observation time.

The 6-hour slot of the satellite lifetime was studied as an example of applying the proposed optimization algorithm to satellite constellation orbit design. The part of the satellite lifetime was selected to compare performance of the proposed semianalytical algorithm to that of typical genetic algorithms in a perspective of computational time. To get optimized results for whole time period of repeating ground track seemed unnecessary since the satellites in the constellation are on the same configuration every about 30 days and therefore much more time for simulation is needed.

Imaging sensors is assumed for this study. The illuminating conditions will change at every pass because the orbit is not Sun-synchronous. However, we believe it is feasible to use imaging sensors in such constellation as existing systems equipped with CMOS or thermal infrared camera.

The constraints for initial orbit elements used as optimization variables are defined in Table 1.

Table 1: Parameters variation for constellation visibility GA simulation.

Variation parameter	Value
Inclination	0-90 degrees
Semimajor axis	6538 to 6978 km, orbit height 160 to 600 km (LEO)
Longitude of ascending node	0-360 degrees
Argument of latitude	0-360 degrees

A comparison between the 4th order Runge-Kutta (4RK) numerical integration and the semianalytical approach is performed as well. The results match within an error of 0.15% to 4.33% over a 1.5- to approximately a 13-day orbit propagation time. In a single application, the semianalytical approach can save up to four times computational time compared to the numerical integration approach, using a computer with an Intel ® Core(TM)2 Duo CPU E7500 2.93 GHz processor and only 2.0 GB RAM. This test was performed using a random target site at 37° latitude and 14°E longitude. A satellite located in orbit with 0.005 eccentricity, a=6812.2 km, i=^37° , and 30° Ω is propagated over five orbital periods. Figure 6 shows the error growth between the 4RK coverage results and the semianalytical approach, which can reach 4 minutes and 42 seconds during the longest 13-day propagation time. From the data, it is predicted that the difference will grow as the propagation time increases.

Figure 6: Coverage time comparison between 4th order Runge-Kutta propagation versus semianalytical approach.

[figure omitted; refer to PDF]

In overall application, the optimization's computational load can be reduced up to nine times using the semianalytical algorithm. If the satellite does not cover any target latitude, the algorithm automatically assigns a high value to the cost function and skips the access computation for longitude sections. Such bad performing individuals are automatically marked by the NSGA-II in the overall computation to exclude this individual in the next generation evaluation. This approach can speed up the entire computational time by skipping at least half of the calculation.

The fitness function of each individual is evaluated for every generation, and while the evaluation does not hit the stopping criteria, the selected parents continue to mate and reproduce new generations via Crossover and Mutation. The stopping criteria of the NSGA-II optimization algorithm used in this paper use MATLAB® default numbers, except for the TolFun or tolerance function being set at 10^-3 (driven by the total number of grids and the definition of the fitness function formula). Optimum constellation design then covers as many grid points as possible and has a minimum average revisit time (ART), also with the lowest maximum revisit time.

Five out of six orbital elements, semimajor axis a, inclination i, longitude of ascending node Ω, argument of perigee ω, and mean anomaly M, are used as optimization design variables, while eccentricity e is kept as a constant. As a note, ω and M are included as design variables, despite the circular orbit assumption, to maximize the range of individual diversity. Using (4) and (5), the drag perturbation effects can be modeled and show increasing θ over propagation time because drag tends to lower the semimajor axis of a satellite. The data for drag is obtained from [35] in Table 2.

Table 2: Drag parameters for a sample satellite (330 km LEO alt., e = 0, i = 79°).

Satellite parameter	Value
Mass (kg)	2.0
Area (m2 )	0.01
C D	2.0

The effect of drag is presented in Figure 7.

Figure 7: Effect of drag perturbation on the θ angle.

[figure omitted; refer to PDF]

There are four optimization objectives in this study, which is divided into two main case studies as follows: the coverage factor and the revisit time factor. The percent coverage for any point on the grid is described as the number of points covered, which is divided by the total number of points in the grid. The percent of coverage directly shows the effectiveness of satellite coverage. However, it cannot provide any information about the distribution of the gaps [32]. The coverage gap is defined here as the length of time where a point is not covered by any of the satellites in the constellation.

In the following, those two figures of merit of satellite constellations are represented in four objectives. Regarding the previous discussion about the convexity of the Pareto solution, the cost function in this study is represented in several objectives instead of using a weighted cost function. The area coverage is different from the time coverage since the time coverage is a summation of total satellite access time divided by the number of grids and the total simulation time. The following algorithm for GA-based optimization is depicted in Figure 8:

(1) Maximum area percent coverage [figure omitted; refer to PDF]

(2) Maximum average area time coverage [figure omitted; refer to PDF]

(3) Minimization of maximum coverage gap time [figure omitted; refer to PDF]

(4) Average coverage gap time [figure omitted; refer to PDF]

Figure 8: Schema for GA-based optimization algorithm.

[figure omitted; refer to PDF]

The reward is defined as a negative value that is added to the cost function if 100% area coverage can be achieved by less than the maximum initial three satellites. Several case studies that were conducted are presented in Figure 9.

Figure 9: Parameters variations for percent coverage and ART by NSGA-II simulation.

[figure omitted; refer to PDF]

Using the semianalytical technique, the time of access accuracy is independent of the time step since it does not use numerical integration such as the Runge-Kutta orbit propagation. In this simulation, test cases are analyzed by varying two types of parameters: the optimization algorithm parameters (or the NSGA-II parameters in this case) and the constellation orbit design parameters. The NSGA-II parameters that can be varied include the population, generation, and the stall criteria. Other adjustments that could be made for Selection, Crossover, and Mutation parameters for the population mating and generation are kept at their default values this time since the effect of changing those parameters has been discussed previously [33, 36].

In Case 2, a different number of populations and generations are tried. By increasing the number of populations and generations, there is a greater probability to achieve a global optimum for any nonlinear optimization problem as the number of search points increases. The analogy for this can be seen in Figure 10, where the increment in the number of populations looks like a mesh figure with an increasing number of mesh points. This results in a bigger search space throughout the overall optimization framework. In our case, it indicates that the Pareto front generated will have a better value, as can be seen in Case 2. However, Case 2 suffers from longest computational time due to the increased number of cost function counting.

Figure 10: Mesh figure for nonlinear problem.

[figure omitted; refer to PDF]

Case 3 iterates the constellation orbit design parameters by increasing the number of satellites in the constellation from 3 to 6. The increment in the number of satellites can increase the area percent coverage and coverage time as well while also reducing both ART and maximum revisit time. However, the increment achieved is still lower than the results achieved in Case 2, even though the computational load is about two times lower. The role of the inclination angle in this satellite constellation design is to maximize the launch vehicle capability by choosing the inclination that is closest to the target latitude.

Cases 5, 6, and 7 all employ i and Ω only as optimized variables. In comparison with Case 5, the number of populations and generations is increased in Case 6, and only two orbital planes are used in Case 7. It can be seen from Figures 14-16 that the largest percent coverage from all four cases (a comparison with Case 1 included) is achieved by Case 6 on 88% coverage. Both the lowest ART and maximum revisit time are achieved by Case 5 at 1.614 and 0.896 hours, respectively. Case 6 achieves the best overall performance in all four cases, although the CPU computational load significantly increases up to 4.5 times the average CPU computational load for all the three cases. To strengthen the results, Case 6 is compared with Case 2, both with 750 population and 100 generations. From these two cases, the results are almost identical in terms of the area percent coverage with only 2% difference. In terms of the time percent coverage, Case 6 exhibits slightly better coverage with 5 minutes difference, while in terms of ART, Case 2 achieves 16.8 minutes shorter time compared to Case 6. Figures 11-16 show the results for various cases.

Figure 11: NSGA-II result for f1 and f2 objectives (Cases 1, 2, 3, and 8).

[figure omitted; refer to PDF]

Figure 12: NSGA-II result for f1 and f3 objectives (Cases 1, 2, 3, and 8).

[figure omitted; refer to PDF]

Figure 13: NSGA-II result for f1 and f4 objectives (Cases 1, 2, 3, and 8).

[figure omitted; refer to PDF]

Figure 14: NSGA-II result for f1 and f2 objectives (Cases 1, 5, 6, and 7).

[figure omitted; refer to PDF]

Figure 15: NSGA-II result for f1 and f3 objectives (Cases 1, 5, 6, and 7).

[figure omitted; refer to PDF]

Figure 16: NSGA-II result for f1 and f4 objectives (Cases 1, 5, 6, and 7).

[figure omitted; refer to PDF]

Case 7 is conducted to seek an optimum design for a constellation of the same three satellites in two orbital planes in order to reduce the launch cost by using a combination of i, Ω, and ω. Seven variables are used in total for three satellites and two orbits, two slots for i, two for Ω, and three for ω variation. However, the separation between two satellites in the same plane is considered too small with only 13° difference in argument of latitude. Therefore, an additional case is attempted using five optimized variables in total, which are a, i, Ω, ω, and ν in Case 12. In Case 12, a constellation of three satellites distributed in two orbital planes is used. In comparison with Case 7, although the area coverage performance is similar, the time coverage achieved in Case 7 is about 15 minutes longer. Also, it can be seen that Case 7 has a greater number of solutions in the Pareto front with near optimum values. From the previous trends, regarding optimization results with respect to the number of optimized variables, one can see that the optimization using NSGA-II leads to a better set of solutions using a lower number of optimized variables (only i and Ω, compared to variations of a, i, Ω, ω, and ν) within a considerably lower number of populations and generations. Figures 17-19 show the results for Cases 7 and 12.

Figure 17: NSGA-II result for f1 and f2 objectives (Cases 7 and 12).

[figure omitted; refer to PDF]

Figure 18: NSGA-II result for f1 and f3 objectives (Cases 7 and 12).

[figure omitted; refer to PDF]

Figure 19: NSGA-II result for f1 and f4 objectives (Cases 7 and 12).

[figure omitted; refer to PDF]

Cases 9 and 10 were conducted to see the results if the number of populations and generations is increased separately. In comparison with Cases 9 and 10, Case 2 yields the best performance in all of the four optimized objectives, while the second-best performance is by Case 9. The results of the simulation are given in Table 3. From Table 4, similar or better performance with Case 9 is achieved by Cases 5, 7, and 12 with approximately 15% less computational load. By increasing the number of populations, the NSGA-II optimization approach allows a greater chance of finding the global optimum prior to converging at a cost of increased computational load. This argument explains why the results from Case 10 are similar to Case 9 despite only half of the function evaluations in Case 9. Figures 20-22 show the results for Cases 1, 2, 9, and 10.

Table 3: Results for the NSGA-II simulation.

Case number	Σ points in Pareto front	Computation time (hr)
1	105	4.104
2	199	16.134
3	125	8.308
4	175	4.279
5	16	4.116
6	263	18.864
7	175	4.098
8	175	3.716
9	175	12.584
10	263	5.386
11	175	3.913
12	175	3.942

Table 4: Best result for all objective functions (_f1 , _f2 , _f3 , _f4 ).

Case number	f 1 (%)	f 2 (min)	f 3 (hour)	f 4 (hour)
1	42	12.40	3.070	2.938
2	86	48.25	1.727	0.899
3	42	10.01	1.951	1.346
4	44	13.05	1.683	1.148
5	64	26.25	1.614	0.896
6	88	53.37	1.715	1.179
7	76	47.54	3.398	1.193
8	58	22.39	3.039	1.359
9	62	21.48	1.682	1.179
10	54	19.29	2.945	1.235
11	56	22.69	3.089	1.240
12	74	32.99	2.493	1.078

Figure 20: NSGA-II result for f1 and f2 objectives (Cases 1, 2, 9, and 10).

[figure omitted; refer to PDF]

Figure 21: NSGA-II result for f1 and f3 objectives (Cases 1, 2, 9, and 10).

[figure omitted; refer to PDF]

Figure 22: NSGA-II result for f1 and f4 objectives (Cases 1, 2, 9, and 10).

[figure omitted; refer to PDF]

Figures 23-25 present a comparison between Cases 1, 4, 5, and 11 results with a different number of orbital elements used for optimization variables. The comparison indicates that Case 5 yields the best results in terms of overall objectives (percent coverage and revisit time) followed by Case 11 with the best area coverage at 56% and then Case 4 at 44%. However, in terms of both the maximum revisit time and ART, the second rank is held by Case 4, while Case 11 offers a longer maximum revisit time at 3.089 hours compared to 1.683 hours in Case 11. The least performing constellation is Case 1 with a variation of a, i, Ω, and ω as optimized variables, while the best performing constellation is Case 5 with only i and Ω as optimized variables. Figures 26-29 show satellites' full ground tracks for various conditions and Tables 5-8 summarize orbit elements for the solutions.

Table 5: Orbit elements for a solution in Case 7 (72% area coverage and 1 hour and 19 minutes of ART).

Orbit element	Sat #1	Sat #2	Sat #3
a (km)	6978	6978	6978
e	0	0	0
i (deg)	42	42	44
Ω (deg)	90	90	131
ω (deg)	324	326	284
ν (deg)	0	0	0

Table 6: Orbit elements for a solution in Case 7 with best maximum revisit time at 3 hours and 24 minutes and 32% area coverage.

Orbit element	Sat #1	Sat #2	Sat #3
a (km)	6978	6978	6978
e	0	0	0
i (deg)	42	42	53
Ω (deg)	154	154	136
ω (deg)	261	274	76
ν (deg)	0	0	0

Table 7: Orbit elements for a solution in Case 12 (best 74% area coverage and 9 hours of ART).

Orbit element	Sat #1	Sat #2	Sat #3
a (km)	6936	6936	6884
e	0	0	0
i (deg)	38	38	75
Ω (deg)	199	199	9
ω (deg)	251	159	158
ν (deg)	335	163	202

Table 8: Orbit elements for a solution in Case 12 with best maximum revisit time at 2 hours and 30 minutes and 44% area coverage.

Orbit element	Sat #1	Sat #2	Sat #3
a (km)	6938	6938	6871
e	0	0	0
i (deg)	39	39	78
Ω (deg)	198	198	118
ω (deg)	251	245	291
ν (deg)	333	342	23

Figure 23: NSGA-II result for f1 and f2 objectives (Cases 1, 4, 5, and 11).

[figure omitted; refer to PDF]

Figure 24: NSGA-II result for f1 and f3 objectives (Cases 1, 4, 5, and 11).

[figure omitted; refer to PDF]

Figure 25: NSGA-II result for f1 and f4 objectives (Cases 1, 4, 5, and 11).

[figure omitted; refer to PDF]

Figure 26: Satellites' full ground track (a) and close-up (b) for a solution in Case 7 with 72% area coverage and 1 hour and 19 minutes of ART. Red line: ground track of sat #1, pink line: ground track of sat #2, and blue line: ground track of sat #3. Green triangle: area covered by sat #1, pink triangle: area covered by sat #2, and blue star: area covered by sat #3.

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

Figure 27: Satellites' full ground track (a) and close-up (b) for a solution in Case 7 with best maximum revisit time at 3 hours and 24 minutes and 32% area coverage. Red line: ground track of sat #1, pink line: ground track of sat #2, and blue line: ground track of sat #3. Red star: area covered by sat #1 and blue triangle: area covered by sat #3. In this case, area covered by satellites #1 and #2 is similar, so the paths are stacked on top of each other.

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

Figure 28: Satellites' full ground track (a) and close-up (b) for a solution in Case 12 (best 74% area coverage and 9 hours of ART). Red line: ground track of sat #1, pink line: ground track of sat #2, and blue line: ground track of sat #3. Pink triangle: area covered by sat #1 and red star: area covered by sat #2.

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

Figure 29: Satellites' full ground track (a) and close-up (b) for a solution in Case 12 with best maximum revisit time at 2 hours and 30 minutes and 44% area coverage. Purple line: ground track of sat #1 and blue line: ground track of sat #2. Red star: area covered by sat #1.

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

In comparison with the best area coverage result in Case 7, it can be seen that, by introducing ν as an optimized variable in Case 12, the satellites' ascending node in the same orbit plane (satellite 1 and 2) is separated up to 264°. Nonetheless, this separation does not create significant changes in terms of percent coverage nor the revisit time performance compared to the results of Case 7.

In summary, the overall objective function can be improved using several methods that are described in Table 9. This table was made by taking the average of performances from the cases that have either the same orbital elements as optimized variables or cases represented by the same number of populations and/or generations.

Table 9: Summary of parameter variation optimization results.

Variation	Percent coverage	Average revisit time	Computational load	Pareto front
Increasing number of populations	Better up to 29%	Better up to 54%	Worse up to 31%	Increasing number of solutions in Pareto front
Increasing number of generations	Better up to 48%	Better up to 60%	Worse up to 206%	N/A
Increasing number of populations & generations	Better up to 71%	Better up to 69%	Worse up to 293%	Increasing number of solutions in Pareto front
Using a variation of a, i, Ω, ω, and ν	Better up to 33%	Better up to 82%	N/A	N/A
Using a variation of a, i, and Ω	Better up to 5%	Better up to 61%	N/A	N/A
Using a variation of i and Ω	Better up to 29%	Better up to 70%	N/A	N/A
Fuel penalty used	Better up to 38%	Better up to 54%	N/A	N/A
Semisparse constellation (Cases 7 & 12)	Better up to 25%	Better up to 61%	N/A	N/A
Increasing number of satellites in constellation	Same value	Better up to 54%	Worse up to 102%	N/A

The maximum coverage time is around 53 minutes on average per satellite orbit, which was achieved in Case 6. Also, 88% maximum area coverage was achieved in Case 6. In terms of the revisit time, the minimum ART is 54 minutes in Case 5, while the lowest maximum gap time is around 1 hour and 37 minutes for Case 5.

6. Conclusion and Future Work

By implementing a semianalytical approach, this study showed that computational load was reduced by over nine times within 0.5% error without having to numerically integrate satellite position or orbital elements over time to get the area coverage. Several figures of merits were chosen to analyze the performance of the optimized constellation orbit design, including maximum coverage, average coverage time, ART, and maximum revisit time over a target area. Circular LEO satellite constellation was used as a case study with a combination of a, i, Ω, and ω as optimization design variables. The satellites in this constellation design case were assumed to have a small conical sensor angle, which means that the satellites were subjected to a limited coverage. Therefore, the study cases were defined as well to achieve maximum coverage efficiency from the available number of satellites. Several cases were analyzed to see how the NSGA-II parameter or the orbital parameter changes might impact the optimizations results. To determine the best parametric operator combination, NSGA-II employed the Pareto concept to give a set of points that are able to satisfy all objectives.

From the conducted study cases, it was shown that results could be improved by increasing the number of populations and generations in NSGA-II or by increasing the number of satellites in the constellation. The latter option was considered inefficient since the results of those two cases were not much different, and the option of adding more satellite in the constellation may cost more. Also, by minimizing the number of optimization variables, better results were obtained since the complexity of the problem was decreased with the decreased number of variables. The percent of coverage was also improved by introducing a fuel penalty. Future work, including analysis of sensor pinpointing accuracy, may yield more meaningful results and make this problem closer to a real case.

Nomenclature

[figure omitted; refer to PDF] :

Reference cross-sectional area of the satellite, m²

[figure omitted; refer to PDF] :

Orbit semimajor axis, km

[figure omitted; refer to PDF] :

Drag coefficient of the satellite

[figure omitted; refer to PDF] :

Eccentricity

[figure omitted; refer to PDF] :

[figure omitted; refer to PDF] -th cost function

[figure omitted; refer to PDF] :

Orbit altitude, km

[figure omitted; refer to PDF] :

Maximum number of satellites in the constellation

[figure omitted; refer to PDF] :

Orbit inclination, deg

[figure omitted; refer to PDF] :

Earth oblateness gravity harmonic coefficient, [figure omitted; refer to PDF]

[figure omitted; refer to PDF] :

Mean argument of latitude, deg

[figure omitted; refer to PDF] :

Mean anomaly, deg

[figure omitted; refer to PDF] :

Mass of the satellite, kg

[figure omitted; refer to PDF] :

Total number of target points

[figure omitted; refer to PDF] :

Orbit mean motion, deg/s

[figure omitted; refer to PDF] :

Number of grids

[figure omitted; refer to PDF] :

Number of individuals

naked [figure omitted; refer to PDF] :

Number of uncovered grids by satellite denoted with [figure omitted; refer to PDF] in percent ratio of all total available grids

[figure omitted; refer to PDF] :

Spherical radius of the Earth, 6378.1363 km

[figure omitted; refer to PDF] :

Satellite's apogee height, km

[figure omitted; refer to PDF] :

Time, s

[figure omitted; refer to PDF] :

Sensor half-angle of conical center, deg

[figure omitted; refer to PDF] :

Elevation angle, deg

[figure omitted; refer to PDF] :

Target location longitude, deg, where [figure omitted; refer to PDF]

[figure omitted; refer to PDF] :

Gravitational parameter of the Earth, [figure omitted; refer to PDF]

[figure omitted; refer to PDF] :

Longitude of ascending node, deg

[figure omitted; refer to PDF] :

Atmospheric density, kg/m³

[figure omitted; refer to PDF] :

Target location latitude, deg, where [figure omitted; refer to PDF]

[figure omitted; refer to PDF] :

Earth-centered half-angle, deg

[figure omitted; refer to PDF] :

Argument of perigee, deg

[figure omitted; refer to PDF] :

The Earth's rotation rate, [figure omitted; refer to PDF] rad/s

[figure omitted; refer to PDF] :

True anomaly, deg.

References

[1] A. da Silva Curiel, A. Cawthorne, M. Sweeting, "Progress in small satellite technology for Earth observation missions,", Small Satellites for Earth Observation: Selected Proceedings of the 5th International Symposium of the International Academy of Astronautics, Berlin, April 4-8 2005 , pp. 50, 2005.

[2] J. G. Walker, "Satellite constellations,", Journal of the British Interplanetary Society , vol. 37, pp. 559-572, 1984.

[3] R. D. Luders, "Satellite networks for continuous zonal coverage,", ARS Journal , vol. 31, no. 2, pp. 179-184, 1961.

[4] T. J. Lang, W. S. Adams, "A comparison of satellite constellations for continuous global coverage,", Mission Design & Implementation of Satellite Constellations , pp. 51-62, Springer, 1998.

[5] J. C. King, Quantization and Symmetry in Periodic Coverage Patterns with Applications to Earth Observation , 1975.

[6] C. Circi, E. Ortore, F. Bunkheila, "Satellite constellations in sliding ground track orbits,", Aerospace Science and Technology , vol. 39, pp. 395-402, 2014.

[7] Y. Ulybyshev, "Satellite constellation design for complex coverage,", Journal of Spacecraft and Rockets , vol. 45, no. 4, pp. 843-849, 2008.

[8] J. Indrikis, R. Cleave, 'SPACE EGGS'-Satellite Coverage Model for Low Earth Orbit Constellations , 1991.

[9] O. L. De Weck, R. D. Neufville, M. Chaize, "Staged deployment of communications satellite constellations in low earth orbit,", Journal of Aerospace Computing, Information and Communication , vol. 1, pp. 119-136, 2004.

[10] Y. N. Razoumny, P. G. Kozlov, V. Y. Razoumny, A. A. Moshnin, "On optimization of Earth coverage characteristics for compound satellite constellations based on orbits with synchronized nodal regression," in Proceedings of the 2nd IAA Conference on Dynamics and Control of Space Systems, Roma, Italy, March 2014.

[11] D. J. Pegher, J. A. Parish, Optimizing Coverage and Revisit Time in Sparse Military Satellite Constellations: A Comparison of Traditional Approaches and Genetic Algorithms , of DTIC Document, 2004.

[12] M. Asvial, R. Tafazolli, B. G. Evans, "Genetic hybrid satellite constellation design,", Parameters , vol. 1, pp. 100, 2003.

[13] W. A. Crossley, E. A. Williams, "Simulated annealing and genetic algorithm approaches for discontinuous coverage satellite constellation design,", Engineering Optimization , vol. 32, no. 3, pp. 353-371, 2000.

[14] C. M. Fonseca, Fleming P. J. "Genetic algorithms for multiobjective optimization: formulation discussion and generalization," in Proceedings of the 5th International Conference on Genetic Algorithms (ICGA '93), pp. 416-423, 1993.

[15] K. Deb, A. Pratap, S. Agarwal, T. Meyarivan, "A fast and elitist multiobjective genetic algorithm: NSGA-II,", IEEE Transactions on Evolutionary Computation , vol. 6, no. 2, pp. 182-197, 2002.

[16] T. Ely, "Satellite constellation design for zonal coverage using genetic algorithms," in Proceedings of the 8th AAS/AIAA Space Flight Mechanics Meeting, pp. 443-460

[17] M. P. Ferringer, D. B. Spencer, "Satellite constellation design tradeoffs using multiple-objective evolutionary computation,", Journal of Spacecraft and Rockets , vol. 43, no. 6, pp. 1404-1411, 2006.

[18] M. P. Ferringer, R. S. Clifton, T. G. Thompson, "Efficient and accurate evolutionary multi-objective optimization paradigms for satellite constellation design,", Journal of Spacecraft and Rockets , vol. 44, no. 3, pp. 682-691, 2007.

[19] M. P. Ferringer, R. S. Clifton, T. G. Thompson, "Constellation design with parallel multi-objective evolutionary computation," in Proceedings of AIAA/AAS Astrodynamics Specialist Conference and Exhibit, Keystone, Colo, USA, 2006.

[20] M. P. Ferringer, D. B. Spencer, P. M. Reed, R. S. Clifton, T. G. Thompson, "Pareto-hypervolumes for the reconfiguration of satellite constellations," in Proceedings of the AIAA/AAS Astrodynamics Specialist Conference and Exhibit, pp. 1-31, August 2008.

[21] W. J. Mason, V. Coverstone-Carroll, J. Hartmann, Optimal Earth Orbiting Satellite Constellations via A Pareto Genetic Algorithm , University of Illinois at Urbana-Champaign, 2001.

[22] G. Confessore, M. Di Gennaro, S. Ricciardelli, "A Genetic Algorithm to Design Satellite Constellations for Regional Coverage,", Operations Research Proceedings , vol. 2001, pp. 35-41

[23] N. Orr, J. Eyer, B. Larouche, R. Zee, "Precision formation flight: the CanX-4 and CanX-5 dual nanosatellite mission," in Proceedings of the 21st Annual AIAA/USU Conference on Small Satellites, 2007.

[24] H. Ashida, K. Fujihashi, S. Inagawa, Y. Miura, K. Omagari, N. Miyashita, S. Matunaga, T. Toizumi, J. Kataoka, N. Kawai, "Design of Tokyo Tech nano-satellite Cute-1.7+APD II and its operation,", Acta Astronautica , vol. 66, no. 9-10, pp. 1412-1424, 2010.

[25] M. Komatsu, S. Nakasuka, "University of Tokyo Nano Satellite Project 'PRISM',", Transactions of Space Technology Japan , vol. 7, 2009.

[26] J. R. Wertz, R. J. Wertz, Mission Geometry; Orbit and Constellation Design and Management: Spacecraft Orbit and Attitude Systems , vol. 13, of Space Technology Library, Microcosm, Boston, Mass, USA; Kluwer Academic Publishers, El Segundo, Calif, USA, 2001.

[27] E. S. H. Hou, N. Ansari, H. Ren, "Genetic algorithm for multiprocessor scheduling,", IEEE Transactions on Parallel and Distributed Systems , vol. 5, no. 2, pp. 113-120, 1994.

[28] D. Sunday, The Convex Hull of a 2D Point Set or Polygon , vol. 109, 2006.

[29] P. Sengupta, S. R. Vadali, K. T. Alfriend, "Satellite orbit design and maintenance for terrestrial coverage,", Journal of Spacecraft and Rockets , vol. 47, no. 1, pp. 177-187, 2010.

[30] Y. Kozai, "The motion of a close earth satellite,", The Astronomical Journal , vol. 64, pp. 367-377, 1959.

[31] H. Schaub, J. L. Junkins, Analytical Mechanics of Space Systems , AIAA, 2003.

[32] J. R. Wertz, Mission Geometry; Orbit and Constellation Design and Management , Space Technology Library, 2001.

[33] H.-D. Kim, O.-C. Jung, H. Bang, "A computational approach to reduce the revisit time using a genetic algorithm," in Proceedings of the International Conference on Control, Automation and Systems (ICCAS '07), pp. 184-189, IEEE, Seoul, Republic of Korea, October 2007.

[34] I. Das, J. E. Dennis, "A closer look at drawbacks of minimizing weighted sums of objectives for Pareto set generation in multicriteria optimization problems,", Structural Optimization , vol. 14, no. 1, pp. 63-69, 1997.

[35] D. Ilyas, Orbital Propagation and Formation Flying of CubeSats within QB50 Constellation , Science and Technology (SpaceMaster), 2011.

[36] T. Savitri, O. Jung, H. Bang, "Optimization of grid points coverage by satellite constellation using genetic algorithm," in Proceedings of the Asia-Pacific International Symposium, Takamatsu, Japan, 2013.