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.

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 [1720]. 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.

The grid points are defined in geodetic latitude and longitude , where and 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 , eccentricity , inclination , longitude of ascending node , argument of perigee , and mean anomaly .

The maximum Earth-centered half-angle and the maximum swath width of the satellite sensor occur when the satellite is at apogee . 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]:

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 is given by [31]where and the mean motion is given byDrag perturbation is described using the following equations [32]:

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]:where 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) becomesFor a noncircular orbit, the maximum swath width achieved on the perigee point is given byIt is recommended to use the time derivative of for a noncircular orbit with drag and perturbation included.

The latitude access time repeats for every period of mean argument of latitude or every . is defined as a set of latitude access times such that

3.2. Longitude Access

The correlation between satellite orbital elements and the area longitude is defined aswhere is the mean argument of latitude and 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 instead of the usual , since the mean argument of latitude lies between and . The analytical solution approach is described briefly by first defining several variables. Let , , and be defined as follows:Substituting (13)–(15) into (12) results in the following equation:

Equation (16) can be solved by using a solver in MATLAB®. is a quadrant-sensitive inverse of , which returns a value in the range . Using Taylor series expansion, can be defined asLet us put a value of into the equivalent Kepler equation in terms of the mean element and derive : is defined as a set of latitude access time for which

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 consists of access times, whereThe flowcharts of the algorithm are depicted in Figures 3 and 4.

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 parameter for the -th objective, as in [33]. where

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.

In the first generation, parent population of individuals creates an offspring population also in individuals through the common tournament Selection, Crossover, and Mutation operators. Next, population with a size of individuals goes through a nondominated sorting. After the first generation, crowded tournament selection is used to select the top individuals from to form the next generation . 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. denotes the number of fast nondominated target points which is the number of grid points located inside the polygon, for which 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.

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™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,  km, , 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.

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 , inclination , longitude of ascending node , argument of perigee , and mean anomaly , are used as optimization design variables, while eccentricity is kept as a constant. As a note, and 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.

The effect of drag is presented in Figure 7.

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(2)Maximum average area time coverage(3)Minimization of maximum coverage gap time(4)Average coverage gap time

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.

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 , 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 . However, Case suffers from longest computational time due to the increased number of cost function counting.

Case 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 , 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 , , and all employ and only as optimized variables. In comparison with Case , the number of populations and generations is increased in Case , and only two orbital planes are used in Case . It can be seen from Figures 1416 that the largest percent coverage from all four cases (a comparison with Case included) is achieved by Case on 88% coverage. Both the lowest ART and maximum revisit time are achieved by Case at 1.614 and 0.896 hours, respectively. Case 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 is compared with Case , 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 exhibits slightly better coverage with 5 minutes difference, while in terms of ART, Case achieves 16.8 minutes shorter time compared to Case . Figures 1116 show the results for various cases.

Case 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 , , and . Seven variables are used in total for three satellites and two orbits, two slots for , 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 , , , , and in Case . In Case , a constellation of three satellites distributed in two orbital planes is used. In comparison with Case , although the area coverage performance is similar, the time coverage achieved in Case is about 15 minutes longer. Also, it can be seen that Case 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 and , compared to variations of , , , , and ) within a considerably lower number of populations and generations. Figures 1719 show the results for Cases and .

Cases and were conducted to see the results if the number of populations and generations is increased separately. In comparison with Cases and , Case yields the best performance in all of the four optimized objectives, while the second-best performance is by Case . The results of the simulation are given in Table 3. From Table 4, similar or better performance with Case is achieved by Cases , , and 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 are similar to Case despite only half of the function evaluations in Case . Figures 2022 show the results for Cases , , , and .

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

In comparison with the best area coverage result in Case , it can be seen that, by introducing as an optimized variable in Case , 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 .

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.

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

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 , , , 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

:Reference cross-sectional area of the satellite, m2
:Orbit semimajor axis, km
:Drag coefficient of the satellite
:Eccentricity
:-th cost function
:Orbit altitude, km
:Maximum number of satellites in the constellation
:Orbit inclination, deg
:Earth oblateness gravity harmonic coefficient,
:Mean argument of latitude, deg
:Mean anomaly, deg
:Mass of the satellite, kg
:Total number of target points
:Orbit mean motion, deg/s
:Number of grids
:Number of individuals
naked:Number of uncovered grids by satellite denoted with in percent ratio of all total available grids
:Spherical radius of the Earth, 6378.1363 km
:Satellite’s apogee height, km
:Time, s
:Sensor half-angle of conical center, deg
:Elevation angle, deg
:Target location longitude, deg, where
:Gravitational parameter of the Earth,
:Longitude of ascending node, deg
:Atmospheric density, kg/m3
:Target location latitude, deg, where
:Earth-centered half-angle, deg
:Argument of perigee, deg
:The Earth’s rotation rate,  rad/s
:True anomaly, deg.

Conflicts of Interest

The authors declare that there are no conflicts of interest.