Spatial Geometry Design of a Low Earth Orbit Constellation for Iranian Regional Navigation Satellite System

The regional navigation satellite system (RNSS) is recently used in some countries to cover or enhance their local navigation. The most important satellite navigation systems to date are in the Earth medium orbit (MEO) and Earth synchronous Earth orbit (GEO), which are characterized by big satellites, and high launch, construction, and operation costs. In contrast, low Earth orbit (LEO) small-satellite constellations have recently become attractive because of their advantages, such as a significant reduction in system cost, an increase in communication volume, and a reduction in latency. Therefore, in this study, the spatial geometry of a LEO constellation is designed for Iran to increase the required regional navigation performance. For this purpose, the optimal constellation configuration is obtained through a multi-objective genetic algorithm (MOGA) utilizing a cost function with a combination of the geometry dilution of precision (GDOP), the number of satellites, and orbital height in the form of a design procedure. Moreover, reducing the feasible region of longitude of the ascending nodes of orbit planes is applied in the design process to reduce the search space. The simulation results indicate the constellation designed performance.


INTRODUCTION
The development of navigation using satellite constellations has been a major area of research and development in recent decades. Especially the Global Navigation Satellite System (GNSS), which is widely designed and implemented by some countries.
The first operational GNSS was the GPS. It was followed by the Global Navigation Satellite System (GLONASS), launched in 1982, and is currently fully operational. Subsequently, other new GNSS, such as Galileo and BeiDou-2, are also deployed. Galileo is the first European GNSS that went live in 2016 (European GNSS Service Centre, www.gpsworld.com, www.gsc.europa.eu).
The second generation of the BeiDou Navigation Satellite System (BDS), also known as COMPASS or BeiDou-2, became operational in China in December 2011 with a partial constellation of 10 satellites in orbit. In 2015, China launched the third generation, BeiDou-3, for global coverage constellation. BeiDou-3 will eventually consist of 35 satellites and is expected to provide global services upon completion in 2020 (Wang 2016). In 2013, India launched the first satellite of the Indian Regional Navigation Satellite System (so called NavIC), it is currently operated with seven satellites, which are designed to provide accurate position information service to users in India and its neighborhoods (Government of India, https://www.isro.gov.in/irnss-programme).
The Japanese Quasi-Zenith satellite system (QZSS) is a four-satellite Regional Navigation Satellite System (RNSS), which is still in development. The first satellite was launched in 2010.
However, these constellations are characterized by high launching, building and operating costs, while the recent approach has been to reduce the total cost of deploying satellite systems. As a result, in recent years, the use of Low Earth Orbit (LEO) systems has become more popular because their overall costs have dropped significantly. The utilization of LEO constellations allows the use of less expensive transmitters, launchers, and mechanical design than their advanced GNSS and RNSS counterparts.
In the field of constellation design, many researchers have taken steps over the years to provide basic designs and models, as well as algorithms for optimizing systems, and to publish the results of their studies. Walker (1970) developed an algorithm for designing satellite constellations, which is commonly used as a starting point for any constellation design. Furthermore, in 1984, he designed the well-known Walker constellation, including the Delta, Sigma, and Omega constellations. The street of coverage (SOC) method, introduced by Lueders and Ginsberg (1974), solves the problem of continuous regional and global coverage by overlapping satellites on adjacent orbit planes. Mortari et al. (2004) designed a better performance constellation by setting up the coverage performance on the basis of the Flower constellation model compared to previous models, and its variants in Lattice Davis et al. 2013), Necklace (Casanova et al. 2014) and repeating ground-track constellations under the effect of the J2-perturbation (Arnas and Casanova 2020) formulations are developed. Ma et al. (2013) presented a LEO constellation for regional communication and considered symmetric and asymmetric structures of the constellation. Zhang et al. (2018) designed a LEO constellation to obtain a Chinese regional navigation augmentation with the aim of minimizing overall system cost. Shtark and Gurfil (2018) introduced a LEO constellation optimization method for RNSS, generated and propagated various constellations, and demonstrated their performance. Li et al. (2019) presented a LEO constellation to improve the performance of a current multi-GNSS real-time positioning service.
It should be noted that, in the LEO constellations, there is a potentially large Doppler shift due to the rapid motion of satellites.
There have been many studies that indicate efforts to address this challenge. Zantou et al. (2005) proposed a simple methodology for Doppler correction, which is adapted for low-cost transceivers. Xu et al. (2018)  Their experiments showed that the LoRa modulation has very high immunity to the Doppler effect. This immunity allows for the use of LoRa modulation in satellite radio communications in orbits above 550 km without any restrictions associated with the Doppler effect. Leyva-Mayorga et al. (2020) provide a comprehensive overview of the physical and logical links, along with the essential architectural and technological components that enable the full integration of LEO constellations into 5G and B5G systems. It seems that some solutions can be adopted to address this challenge.
On the other hand, the optimization of a constellation is a hard problem because of its nonlinearity, nonuniqueness and discontinuity nature. So, it is difficult to find an optimal solution using analytic approaches. For this purpose, modern algorithms, such as evolution-based or intelligence swarm-based approaches, are used to find the near-optimal solution. Crossley et al. (2000) proposed a simulated annealing and genetic algorithm to minimize the revisit time. Zhang et al. (2018) designed the LEO satellite constellation using differential evolution (DE) algorithm. Shtark and Gurfil (2018)  (MOGA) for optimizing LEO constellations for regional coverage. Quan and Chao (2007) applied the ant system algorithm (ASA) to optimize the constellation coverage performance.
In this paper, a design procedure was developed for designing the geometric part of a satellite constellation of possible future Iranian regional navigation satellite system (IRNSS), based on MOGA optimization.

THE PROBLEM DEFINITION
Assuming that the region of Iran is covered by a quadrilateral, the vertices of the quadrilateral are considered as a set of four ground places PL = {pl j };j = 1...4. Each place pl j is defined by its latitude, longitude and altitude, λ plj , φ plj and H plj , respectively. Figure 1 depicts this quadrilateral region. Regarding three segments presented in every space mission, the space segment (in this case the satellite constellation), the ground segment (which are responsible for major control and management tasks of the space segment), and the user segment (which refers to the rest of the communication devices at ground level), the definition of the problem is the geometric design of the satellite constellation to achieve regional navigation in such a way that the positioning precision of each point of the given area at any time is accurate and cost-effective.

Mathematical formulation
This section formulates the related parameters for a satellite constellation model. For this purpose, only homogeneous satellite constellations, which include a number of orbit planes with equal satellites in each plane and the same orbits, are considered. In addition, circular orbits are assumed and the optimization of orbits (i.e., optimization of orbit altitude H P , inclination angle inc, the right ascension of the ascending node [RAAN] W and argument of perigee ω for each satellite) are combined into the number of orbital planes N P and the number of satellites per orbit N S .
It should be noted that the inclination inc is greater than the latitude of ground places because, otherwise, the satellite will not pass over the ground places, i.e., inc > φ plj , j = 1...4.
Following is a calculation of Ω and ω for any orbit at any time t, respectively (Eq. 1): (1) where t 0 is an epoch time, Ω • , ω • , and M • are the mean rates of the ascending node, the argument of perigee and mean anomaly due to effects of zonal harmonics perturbation, respectively, defined as Eq. 2: where J 2 is the Earth second zonal harmonic coefficient, R e is the Earth mean equatorial radius, a = R e + H p is the semi major, and μ is the standard gravitational parameter.
The positions of the available satellites in relation to the ground receivers play an important role in the accuracy of the position estimation, and this affect is measured by Geometry Dilution of Precision (GDOP). It is composite measure taking into account the geometrical impacts on both time and position. In general, GDOP can range from one to about 100. The mean value for the GPS constellation is about 2.7 (Wertz 2001). A larger GDOP value means poor performance for that measurement.
To compute the GDOP value, at least four satellites in view from any given places are required. Generally, a satellite is determined to be visible from a place if the elevation angle e (the angle measured at the ground place between the spacecraft and the local horizontal) is greater than a given minimum value. In this paper, ε min = 5° is assumed (Wertz 2001).
Calculations of GDOP values are performed. First, r ik is defined as the position vector of a satellite i with the argument of latitude θ ik = ω k + v ik (v ik is the satellite true anomaly) in a circular orbit k with an inclination angle inc, the altitude H p , the argument of perigee ω k , and the longitude of ascending node λ k . The following relation is used to calculate r ik in the ECEF frame (Wertz 2001 where λ k is obtained as in Eq. 4 (Vallado 2007): (4) In this relation, θ G is the Greenwich sidereal time that is obtained from Julian time of t. Next, ρ ij is obtained, which is the relative position of each satellite i with respect to the place pl j expressed the Earth centered Earth fixed (ECEF) frame, as follows (Eq. 5): It would be noted that, to calculate the elevation angle of the places, it is necessary to express δ ij in the SEZ frame. For any M visible satellites, the offset δr in the position and time bias of the receiver are related to the offset in the pseudorange measurements δρ through the linearized relation (Noureldin et al. 2013) (Eq. 6): where G is the geometry matrix with M × 4 dimensions, which characterizes the relative geometry of a satellite and the receiver at any time, written as Eq. 7: (7) Hence, the least squares solution for Eq. 6 is given by Eq. 8: (8) GDOP is characterized by the geometry matrix G, which relates the parameters of the user's position and time bias errors to those of the pseudorange errors. Thus, the GDOP related to place pl j is obtained at any time t using the following relation (Eq. 9): Obviously, the number of visible satellites from the relevant location must be at least four, so that the relation (9) is not singular.
In other words, for GDOP(t, pl j ) to be defined, G j must consist of at least four rows.

The Doppler shift
The Doppler shift is produced by the relative motion of the satellite with respect to the user. At the receiver antenna, the received frequency, f R , can be approximated by the classical Doppler equation as follows (Eq. 10): (10) where f T is the transmitted satellite signal frequency, V r is the satellite-to-user relative velocity vector in the ECEF frame, 1 ρ is the unit vector pointing along the line of sight from the user to the satellite, and c is the speed of propagation.
The Doppler offset (shift), due to the relative motion, is obtained from these relations as Eq. 11 (Kaplan and Hegarty 2006): The received frequency increases as the satellite approaches the receiver and decreases as it recedes from the user. At the GPS L1 frequency, the maximum Doppler frequency for a stationary user on the Earth is approximately 4 kHz, corresponding to a maximum line-of-sight velocity of approximately 800 m·s -1 (Kaplan and Hegarty 2006).

The nonlinear multi-objective cost functional
In this study, the GDOP values and the number of satellites, as well as orbits altitude, are considered as objects of the cost functional. Thus, the constellation design problem can be described as follows (Eq. 12): where m is the number of objective functions and p is the number of constraints. X ∈ n is a vector of design variables, where n is the number of independent variables. C(x) ∈ n is a vector of objective functions (criteria).
The predominant concept in defining an optimal point is that of Pareto optimality (Ehrgott 2005), which is defined as follows: Definition 1: Pareto optimal -A feasible solution x* ∈ X (X is the search space) is efficient (Pareto optimal) solution if there does not exist another vector x ∈ X, such that C(x) ≤ C(x*), and C i (x) < C i (x*), i = 1,..., m, for at least one function.
In this work, the components of C(x) are defined as (Eq. 13): where the criterion C 1 (x) represents the mean GDOP of the set of at least four visible satellites in the designed constellation from given places pl j ∈ PL at overall time. C 2 (x) = N T = N p × N S is the total number of satellites in the constellation, and C 3 (x) = H P is the orbit height.
The problem of constellation design is solved by searching the minimum of the defined objective function in the range of the design variables of the constellation.

Reduction of the search space
This section investigates the variations of the longitude of the ascending node (λ AN ) at any time for every orbit to reduce its range in the search space. Similar to the coverage region concept proposed by Ulybyshev (2015), the visibility conditions for the satellite in the space can be presented by an oval of boundary points defined in a two-dimensional space, wherein the x-axis is the λ k , and the y-axis is the time. The satellite trajectory map in the space is near a straight line along the y axis. If the line is intersected by the oval, then there is a visibility interval. Figure 2 shows the double connected fourfold coverage regions of the given four places for a day. Clearly, the grids area in the figure indicate, at any sampling time, at least four satellites in the constellation are simultaneous visible from any of the specified places around Iran. For the problem considered in this paper, Δl = Δl 1 + Δl 2 is about 221° as shown in Fig. 2. Consequently, the maximum value of the variable Δl is reduced to 221°, and the λ-range of the search space is reduced.

GENETIC ALGORITHM OPTIMIZATION SOLVER
Genetic algorithm (GA), as one of the evolutionary algorithms, has shown its effectiveness as an optimization method for complex multiobjective problems. The use of GA is a quickly evolving field of research, and there is much new to recommend. This method is developed based on biological evolution and Darwin' s theory. Similar to biological evolution, the GA starts with an initial population of random individuals (members). Then selects them with a good chance of reproduction (based on the best fitness performance) and reproduces a new generation of individuals using crossover and mutation operations. This process is repeated several times until at least one individual has a better fitness than the defined minimum fitness, or a maximum number of generations have been reached. Figure 3 shows the flowchart of a GA. An individual's parameters for optimization are displayed by a structure called a chromosome. These parameters are also called decision variables. To design a satellite constellation, each chromosome represents a constellation model, and the optimized parameters are the constellation design variables: number of orbit planes N P , number of satellites per orbit N S , the orbital height H p , the orbital inclination inc, the initial orbital argument of perigee ω 0 and initial orbital RAAN Ω 0 . For other orbital parameters, according to the predicted simplifications, the orbital eccentricity e = 0 and the mean anomaly M are evenly distributed for all constellations, that is used for obtaining the argument of perigee (Eq. 1). Thus, MOGA for the optimization problem consists of the vector of six decision variables x = [N p N S H P inc ω 0 Ω o ] and the multi-objective cost functional (fitness function), as presented in Eq. 13. GDOP values are also calculated using Eqs. 1 to 9 with the initial conditions presented in the next section.
It should be noted that, of these six parameters, the first two are of the integer type and the rest are of the real positive type. Therefore, the optimization problem is of the mixed-integer nonlinear programming type.

SIMULATION RESULTS
The optimization algorithm described in last section was implemented in MATLAB-R2019b and tested on a computer with 2.5 GHz Intel Core i5 and 4Gbytes RAM running on the 64bit-Windows 10 operating system. Table 1 contains the geographical coordinates of the four places (A, B, C, D) surrounding the region of Iran. To avoid the radiation emitted by the Van Allen belt, which may affect the satellite components, the upper bound of the altitude is considered 1000 km. Furthermore, to reduce the atmospheric drag and thus increase the mission lifetime, its lower bound is considered 500 km. The rest of upper bound and lower bound of decision variables of the optimization problem is as in Eq. 14: Parameters of the MOGA are chosen as follows: Number of population: N Population = 100; Number of generations: N Generation = 200; Probability of crossover operation: P Cross = 0.8; Probability of mutation operation: P Mutation = 0.05. MOGA simulation code is executed at the epoch of April 2 nd , 2020, 07:30 h. The sampling time step of 60 s is assumed. Figure 4 depicts Pareto front results of the MOGA. In this figure, three optimization objectives are illustrated in pairs. Table 2 presents results of the 13 favorite cases of the constellations based on optimal decision variables of the optimization problem.   Because this method searched the optimal constellations in terms of only the worst GDOP values of all the visible places at any sampling instant of time, it would be necessary to see the evolution in time of the minimum, mean and maximum GDOPs experienced by all the four places in the region of interest. Figure 5 shows these GDOP values of constellations presented in Table 2, which is sorted by number of satellites. It is shown the mean GDOP experienced by the regional places is below 5.0.   It should be noted that the number of satellites and orbital heights have a significant effect on the GDOP values of the constellations. Consequently, the designer is able to select economical to moderate constellations based on launch capacity and overall cost. For instance, the case No. 2 will be an admissible choice (as an economic case with GDOP much less than the first case and a smaller number of satellites than other cases), and case No. 6 will be another choice (as the first case from the top to the bottom of the table whose GDOP is less than 2). Figure 6 illustrates 3D schematic of optimal constellation No. 2. Figure 7 depicts resulting typical Doppler shift versus time for receivers overall 24 h. According to the figure, the maximum Doppler shift 33 kHz is achieved for all receivers based on the constellation No. 2 and L1 frequency (1575 MHz). Therefore, it is obvious that this value must be covered or compensated by the receivers. CONCLUSION LEO satellites are more accessible due to their lower altitudes and have advantages over a geostationary satellite, such as lower cost and low-latency with higher bandwidth communications. In this work, the spatial geometric part of a LEO satellite constellation for Iranian RNSS is designed to obtain the required navigation performance. For this purpose, a mathematical model of a satellite constellation with its design parameters is introduced and a MOGA optimization utilizing a cost function with a combination of the GDOP, the number of satellites and the orbit altitude is developed. The created model, along with the optimization module, provides a design procedure for the development of any system in the software environment. Moreover, the feasible region of the orbit planes is reduced in the design process to narrow the search space by choosing the visibility range of RAAN. Finally, a set of Pareto optimal solutions are obtained, as summarized in Table 2. These solutions are classified from economical to expensive constellations based on the orbital height, number of satellites and navigation accuracy (GDOP) so that one can be chosen by tradeoff. Given that the GDOP should also be the minimum value, it can be concluded that case No. 2 with 130 satellites in the height of 868 km will be an admissible choice (as an economic case with GDOP much less than the first case).
Regarding the results of the Doppler effect, depending on the type of modulation, signal bandwidth and power level, it is expected that the effect will be covered or compensated by the terminal receivers, which is beyond the scope of this research.
In addition, given that the region of the country is located within the assumed quadrilateral, any ground point of this region, such as the command and control centers, sensor stations, monitor stations, and mission center, are within the territory, and thus will be covered by the constellation.
Future work will continue with the development of a regional/global LEO constellation optimization design combination with GEO/IGSO/MEO.

DATA AVAILABILITY STATEMENT
All datasets were generated or analyzed in the current study.