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

Authors: Tania Savitri; Youngjoo Kim; Sujang Jo; Hyochoong Bang
Journal: International Journal of Aerospace Engineering
 (2017)

---

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


---

## Body

## 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 Figure1. 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)
(b)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 axisa, eccentricity e, inclination i, longitude of ascending node Ω, 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=a1+e. As shown in Figure 2, a satellite footprint projection on the Earth’s surface with elevation angle ε defines a coverage circle with a radius θ through the following equation [29]:(1)θmax=-α+asin⁡a1+eResin⁡α.Figure 2
Satellite footprint projection.


## 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 Figure1. 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)
(b)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 axisa, eccentricity e, inclination i, longitude of ascending node Ω, 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=a1+e. As shown in Figure 2, a satellite footprint projection on the Earth’s surface with elevation angle ε defines a coverage circle with a radius θ through the following equation [29]:(1)θmax=-α+asin⁡a1+eResin⁡α.Figure 2
Satellite footprint projection.


## 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](2)Ω˙=-Ka7/21-e22cos⁡i,ω˙=Ka7/21-e222-52sin2⁡i,M˙=n+Ka7/21-e23/21-32sin2⁡i,where K=3/2J2Re2μ and the mean motion n is given by(3)n=μa3.Drag perturbation is described using the following equations [32]: (4)a˙d=-ρACdμam1-ωencos⁡i2.


### 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]:(5)ϕ-π2+cos-1⁡sin⁡icos⁡ϕs=0,(6)ϕs=3π2+ω0+ν˙Δt,where ϕ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(7)ϕ∓θmax-π2+cos-1⁡sin⁡icos⁡ϕs=0.For a noncircular orbit, the maximum swath width θmax achieved on the perigee point is given by(8)θmax=-αhalf+sin-1⁡a1+eREsin⁡αhalf.It is recommended to use the time derivative of θ for a noncircular orbit with drag and J2 perturbation included. (9)θ˙=sin⁡αhalf/Rer/aa˙+r2/pesin⁡ff˙1-rsin⁡αhalf/Re.The latitude access time repeats for every period of mean argument of latitude or every2π/ω˙+M˙. Tlati is defined as a set of latitude access times such that(10)Tlati=pini,pouti,i=1,…,m.


### 3.2. Longitude Access

The correlation between satellite orbital elements and the area longitude is defined as(11)Ψ-Ω0+ω0-GMST0∓θmax-Ω˙+ω˙-ωEll˙=atan⁡2cos⁡isin⁡l,cos⁡l,where l is the mean argument of latitude and ∓θ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:(12)A1=Ψ-Ω0+ω0-GMST0∓θmax,(13)A2=Ω˙+ω˙-ωEl˙,(14)A3=cos⁡i.Substituting (13)–(15) into (12) results in the following equation:(15)A1-A2l-atan⁡2A3sin⁡l,cos⁡l=0.Equation (16) can be solved by using a solver in MATLAB®. α=atan2y,x is a quadrant-sensitive inverse of tan⁡α=y/x, which returns a value in the range -π≤α≤π. Using Taylor series expansion, atan⁡A3sin⁡l,cos⁡l can be defined as(16)atan⁡x=x-x33+x55-x77+⋯,-1≤x≤1,(17)x=A3sin⁡l,cos⁡l.Let us put a value of l into the equivalent Kepler equation in terms of the mean element and derive Δt:(18)l=ω0+M0+ω˙+M˙Δt.Tlonj is defined as a set of latitude access time for which(19)Tlonj=qinj,qoutj,j=1,…,n.The time for intersection of a target point can be found from the time window intersection of each latitude and longitude access. The intersection timeTaccess consists of s access times, where(20)Taccess=⋃i=1mTlati,Tloni=Tink,Toutk,k=1,…,s.The flowcharts of the algorithm are depicted in Figures 3 and 4.Figure 3
Latitude and longitude target access determination flowchart.Figure 4
Target access time intersection flowchart.


## 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]:(5)ϕ-π2+cos-1⁡sin⁡icos⁡ϕs=0,(6)ϕs=3π2+ω0+ν˙Δt,where ϕ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(7)ϕ∓θmax-π2+cos-1⁡sin⁡icos⁡ϕs=0.For a noncircular orbit, the maximum swath width θmax achieved on the perigee point is given by(8)θmax=-αhalf+sin-1⁡a1+eREsin⁡αhalf.It is recommended to use the time derivative of θ for a noncircular orbit with drag and J2 perturbation included. (9)θ˙=sin⁡αhalf/Rer/aa˙+r2/pesin⁡ff˙1-rsin⁡αhalf/Re.The latitude access time repeats for every period of mean argument of latitude or every2π/ω˙+M˙. Tlati is defined as a set of latitude access times such that(10)Tlati=pini,pouti,i=1,…,m.


## 3.2. Longitude Access

The correlation between satellite orbital elements and the area longitude is defined as(11)Ψ-Ω0+ω0-GMST0∓θmax-Ω˙+ω˙-ωEll˙=atan⁡2cos⁡isin⁡l,cos⁡l,where l is the mean argument of latitude and ∓θ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:(12)A1=Ψ-Ω0+ω0-GMST0∓θmax,(13)A2=Ω˙+ω˙-ωEl˙,(14)A3=cos⁡i.Substituting (13)–(15) into (12) results in the following equation:(15)A1-A2l-atan⁡2A3sin⁡l,cos⁡l=0.Equation (16) can be solved by using a solver in MATLAB®. α=atan2y,x is a quadrant-sensitive inverse of tan⁡α=y/x, which returns a value in the range -π≤α≤π. Using Taylor series expansion, atan⁡A3sin⁡l,cos⁡l can be defined as(16)atan⁡x=x-x33+x55-x77+⋯,-1≤x≤1,(17)x=A3sin⁡l,cos⁡l.Let us put a value of l into the equivalent Kepler equation in terms of the mean element and derive Δt:(18)l=ω0+M0+ω˙+M˙Δt.Tlonj is defined as a set of latitude access time for which(19)Tlonj=qinj,qoutj,j=1,…,n.The time for intersection of a target point can be found from the time window intersection of each latitude and longitude access. The intersection timeTaccess consists of s access times, where(20)Taccess=⋃i=1mTlati,Tloni=Tink,Toutk,k=1,…,s.The flowcharts of the algorithm are depicted in Figures 3 and 4.Figure 3
Latitude and longitude target access determination flowchart.Figure 4
Target access time intersection flowchart.


## 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 weightwi parameter for the i-th objective, as in [33]. (21)J=minx⁡∑i=1nwifix,x∈C,(22)J=minx⁡1-αf1x+αf2x,x∈C,where(23)α=sin⁡θcos⁡θ+sin⁡θ.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 Figure5, 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.In the first generation, parent populationPt 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.


## 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 weightwi parameter for the i-th objective, as in [33]. (21)J=minx⁡∑i=1nwifix,x∈C,(22)J=minx⁡1-αf1x+αf2x,x∈C,where(23)α=sin⁡θcos⁡θ+sin⁡θ.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 Figure5, 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.In the first generation, parent populationPt 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 Table1.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 degreesA 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™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.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 axisa, 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.0The effect of drag is presented in Figure7.Figure 7
Effect of drag perturbation on theθ angle.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 Figure8:(1)
Maximum area percent coverage(24)f1=1-∑i=1NcoverageidxN+reward.(2)
Maximum average area time coverage(25)f2=1-∑i=1Ncoverage  timeN·simulation  period.(3)
Minimization of maximum coverage gap time(26)f3=∑i=1Nmax⁡revisitN·60.(4)
Average coverage gap time(27)f4=∑i=1NrevisitN·60.Figure 8
Schema for GA-based optimization algorithm.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  Figure9.Figure 9
Parameters variations for percent coverage and ART by NSGA-II simulation.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 Case2, 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.Case3 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.Cases5, 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 forf1 and f2 objectives (Cases 1, 2, 3, and 8).Figure 12
NSGA-II result forf1 and f3 objectives (Cases 1, 2, 3, and 8).Figure 13
NSGA-II result forf1 and f4 objectives (Cases 1, 2, 3, and 8).Figure 14
NSGA-II result forf1 and f2 objectives (Cases 1, 5, 6, and 7).Figure 15
NSGA-II result forf1 and f3 objectives (Cases 1, 5, 6, and 7).Figure 16
NSGA-II result forf1 and f4 objectives (Cases 1, 5, 6, and 7).Case7 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 forf1 and f2 objectives (Cases 7 and 12).Figure 18
NSGA-II result forf1 and f3 objectives (Cases 7 and 12).Figure 19
NSGA-II result forf1 and f4 objectives (Cases 7 and 12).Cases9 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.942Table 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.078Figure 20
NSGA-II result forf1 and f2 objectives (Cases 1, 2, 9, and 10).Figure 21
NSGA-II result forf1 and f3 objectives (Cases 1, 2, 9, and 10).Figure 22
NSGA-II result forf1 and f4 objectives (Cases 1, 2, 9, and 10).Figures23–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 Case7 (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
0Table 6
Orbit elements for a solution in Case7 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
0Table 7
Orbit elements for a solution in Case12 (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
202Table 8
Orbit elements for a solution in Case12 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
23Figure 23
NSGA-II result forf1 and f2 objectives (Cases 1, 4, 5, and 11).Figure 24
NSGA-II result forf1 and f3 objectives (Cases 1, 4, 5, and 11).Figure 25
NSGA-II result forf1 and f4 objectives (Cases 1, 4, 5, and 11).Figure 26
Satellites’ full ground track (a) and close-up (b) for a solution in Case7 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)
(b)Figure 27
Satellites’ full ground track (a) and close-up (b) for a solution in Case7 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)
(b)Figure 28
Satellites’ full ground track (a) and close-up (b) for a solution in Case12 (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)
(b)Figure 29
Satellites’ full ground track (a) and close-up (b) for a solution in Case12 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)
(b)In comparison with the best area coverage result in Case7, 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 Table9. 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 ofa, i, Ω, ω, and ν
Better up to 33%
Better up to 82%
N/A
N/A
Using a variation ofa, i, and Ω
Better up to 5%
Better up to 61%
N/A
N/A
Using a variation ofi 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 (Cases7 & 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/AThe maximum coverage time is around 53 minutes on average per satellite orbit, which was achieved in Case6. 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 ofa, 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.



---

*Source: 1235692-2017-05-31.xml*