Mars Constellation Design and Low-Thrust Deployment Using Nonlinear Orbit Control

This research addresses the design of a Mars constellation composed of 12 satellites and devoted to telecommunications. While 3 satellites travel areostationary orbits, the remaining 9 satellites are placed in three distinct quasi-synchronous, inclined, circular orbits. The constellation at hand provides continuous global coverage, over the entire Martian surface. The use of 4 carrier vehicles, departing from a 4-sol orbit, is proposed as an affordable option for the purpose of deploying the entire constellation, even starting from dispersed initial conditions. Each carrier is driven toward the respective operational orbit using steerable and throttleable low-thrust propulsion, in conjunction with nonlinear orbit control. Lyapunov stability analysis leads to defining a feedback law that enjoys quasi-global stability properties. Orbit phasing concludes the constellation deployment, and is carried out by each satellite. The tradeoff between phasing time and propellant expenditure is characterized.


Introduction
Mars represents a major target for planetary exploration in the decades to come. Several scientific missions are already planned, and the prospect of carrying out the first human mission is becoming more and more concrete. A variety of 1 3 objects, ranging from rovers to scientific instrumentation for environmental analyses, were (and are being) released in different locations on the Martian surface. A satellite constellation aimed at providing continuous, global coverage over the great majority of the red planet could serve as a communications relay system and would certainly represent a valuable asset for future missions. In fact, such a constellation would significantly facilitate the communications among orbiters, probes and instruments located on the surface, and Earth ground stations. This would also lead to extending the existing capabilities of NASA's Deep Space Network, as remarked in Ref. [1]. Several goals were defined [1] for a satellite constellation about Mars: (a) global coverage over a prescribed time span, (b) communication support of the equatorial region, also as a support to human missions, (c) maximization of the communication and navigation performance across all latitudes and longitudes, (d) minimization of the variations of the navigation and communication performance, (e) redundant coverage in case of loss of any single satellite, and (f) minimization of the coverage variability due to long-term orbit perturbations, with special reference to the orbit plane precession.
Depending on the specific tasks, several types of orbits may be suitable for the purpose of designing and deploying a constellation around Mars [2,3]. Bell et al. [2] proposed 4 distinct constellation configurations of microsatellites that travel low-altitude, inclined orbits and ensure satisfactory global communication and navigation performance, with special regard to the equatorial region. Each configuration includes 6 satellites with different inclinations. Yet, for all configurations the maximum revisit time ranges from 1 to 8 h, and exceeds 2 h over the great majority of the Martian surface. Nann et al. [4] focused on precise characterization of the Martian atmosphere for future missions, and designed a constellation of 8 satellites suitable for the purpose of performing radio occultation measurements. The same constellation can be reconfigured, to provide (discontinuous) navigation services. Most recently, Kelly and Bevilacqua [5] proposed a constellation of 15 satellites that ensures global and continuous coverage of the entire Martian surface, using 5 equally spaced, inclined orbit planes.
Release of multiple satellites from a single carrier spacecraft would be an attractive and affordable option for constellation deployment. However, this objective is rather challenging, especially if the constellation configuration is designed to include different orbit planes with large angular displacements. Former contributions in the scientific literature were mainly focused on deployment strategies for Earth constellations [6][7][8][9]. In 2006, the FORMOSAT-3/COSMIC constellation [7] was released from a single launch vehicle. Each of the 6 satellites raised its orbit, and leveraged the differential precessional motion due to the Earth oblateness, with the final aim of achieving the desired orbit plane.

3
The Journal of the Astronautical Sciences (2022) 69:  This well-known effect was also investigated by Crisp et al. [8] In addition, they designed a second strategy, based on the temporary release of the spacecraft in the proximity of the interior libration point of the Earth-Moon system. Most recently, Lee et al. [9] focused on a flexible deployment strategy in response to the possible evolution of the Earth regions of interest, whereas Pontani and Teofilatto [10] focused on low-thrust deployment of a low-Earth-orbit constellation tailored to polar ice monitoring.
This research considers a constellation composed of 12 satellites released in 4 distinct orbit planes about Mars. While 3 satellites travel an areostationary orbit, the remaining 9 satellites are released in three repeating, circular, quasi-synchronous, inclined orbits. This configuration is aimed at guaranteeing global and continuous coverage, while ensuring repetition and predictability of the visible passes as well as visibility of multiple satellites. This work also addresses the problem of deploying all the satellites using 4 carrier vehicles, released by the mother spacecraft. The latter orbits a highly elliptical path, i.e. the ESA 4-sol orbit [11], which is entered after planetary capture. Orbit dynamics about Mars is modeled with the inclusion of the most relevant perturbations, i.e. several harmonics of the areopotential, together with the gravitational pull due to the Sun as a third body. Each carrier embarks 3 satellites and employs steerable and throttleable low-thrust propulsion for the purpose of driving them toward the respective operational orbit. Nonlinear orbit control is proposed as a suitable strategy for autonomous real-time guidance, even in the presence of nonnominal flight conditions, e.g. if the initial highly elliptical orbit is significantly different from the expected one, due to injection errors at the planetary capture. Numerical simulations are intended to prove that the approach at hand is effective and accurate for releasing 4 carrier vehicles in the respective operational orbits. Phasing represents the last phase of deployment, and is performed by each satellite. The final goal is in gaining the correct angular position, without leaving the desired orbit plane. To this end, both internal phasing (using an intermediate orbit with shorter period) and external phasing (using an intermediate orbit with longer period) are being considered.
In short, the present research has the following objectives: (i) design a constellation of 12 satellites, capable of guaranteeing continuous global coverage of the entire Martian surface, with simultaneous visibility of multiple satellites, (ii) propose the use of 4 carrier vehicles for the release of the 12 satellites in 4 distinct orbit planes, (iii) use an effective low-thrust deployment strategy based on nonlinear orbit control, in the presence of nominal and nonnominal initial conditions, and (iv) describe the phasing maneuver, performed by each satellite, while providing a comprehensive analysis of the tradeoff between overall velocity change and phasing time.

Constellation Configuration
This research proposes a satellite constellation about Mars tailored to telecommunications. Global coverage of the Martian surface is sought, using 12 satellites placed in 4 distinct orbit planes, i.e. 3 satellites in areostationary orbit, and 1 3 9 satellites in 3 quasi-synchronous, circular, repeating-ground-track, inclined orbits. High-altitude orbits are considered, and no atmospheric drag is encountered. For the purpose of constellation design, the solar radiation pressure and the third body gravitational pull due to the Sun are negligible as well. In contrast, the Mars oblateness ( J 2 zonal harmonic) produces significant effects on the orbit elements and is considered in the model.
Desirable features of a satellite constellation are (a) repetitiveness, which implies predictability of the performance, and (b) invariance with respect to the J 2 perturbation. The use of circular orbits is a common and convenient requirement, specifically related to onboard sensors. In fact, if the altitude is constant, the sensor performance and operational modality does not change.

Coverage Analysis for 3 Satellites in Areostationary Orbit
The first three satellites that form the constellation travel areostationary orbits, at equally-spaced geographical longitudes. However, similarly to what occurs for geostationary satellites, areostationary spacecraft are subject to the perturbation related to harmonic J 2,2 of the areopotential. This term is responsible for the existence of two stable and two unstable geographical longitudes for an aerostationary satellite. The two unstable Martian longitudes equal − 105° and 75° [12]. Therefore, the following equally-spaced geographical longitudes are chosen for the three satellites placed in areostationary orbit: − 135°, − 15°, 105°. These values correspond to maximizing the angular separation of the subsatellite points from the unstable longitudes. The visible regions from each satellite are portrayed in Fig. 1, with different colors associated with distinct elevation angles. The black curves represent the bound of the visible regions with elevation ≥ min = 0 • , whereas the blue and red curves correspond to min = 20 • and 40°, respectively. Inspection of Fig. 1   areostationary orbits would certainly widen the latitudinal range associated with coverage with ≥ min = 20 • , but would not alter the latitudinal intervals (in the two hemispheres) where no satellite is visible with ≥ min = 20 • . In conclusion, 3 satellites placed in areostationary orbits, at equally-spaced longitudes, represent a viable option to continuously cover the equatorial region. Nevertheless, the use of inclined orbits for global continuous coverage is ineludible.

Quasi-synchronous, Circular, Repeating-Ground-Track, Inclined Orbits
While 3 satellites are placed in an equatorial areostationary orbit, the remaining 9 satellites travel equally-inclined circular orbits with an identical altitude. On average, oblateness effects ( J 2 zonal harmonic perturbation) yield the following time derivatives [13] of the right ascension of the ascending node (RAAN) Ω , argument of periapse , and mean anomaly M: In Eqs. All the satellites are assumed to have identical inclination and SMA (or altitude H for circular orbits) in order to avoid differential actions due to the J 2 perturbation. This circumstance would cause a substantial alteration of the performance attainable by the constellation.
An orbit is termed repeating when phased with the planet rotation, i.e. when the ground track is periodically repeated. This occurs if the satellite completes N t orbits in m nodal days, In Eq. (5) T n = 2 ̇ is the nodal orbital period (i.e., the time interval between two consecutive ascending node crossings) and D n = 2 M −Ω is the nodal day (i.e. the time required for Mars to make a complete rotation with respect to the orbital plane; M denotes the Mars rotation rate). Both T n and D n depend on a and i (whereas e = 0 ). This fact implies that the ratio r t = N t m , which represents the number of orbits per nodal day, depends on a and i.
Quasi-synchronous satellites are stationary with respect to a specific location on the Martian surface. This location is selected as the one corresponding to the maximum latitude flown by the satellite, denoted with M . For direct orbits, the eastward component of the velocity at the maximum (and minimum) latitude coincides with the orbital velocity, i.e. √ M ∕R . Quasi-synchronism holds if the angular rate of an observer at latitude M equals the eastward angular velocity of the satellite, i.e.
In the last step, M ≡ i was used. Figure 2 portrays the radius of quasi-synchronous circular orbits as a function of inclination (in blue), together with the related number of orbits per nodal days q ∶= N t m (orange curve). Each point along the blue curve corresponds to a circular orbit whose motion is synchronous with the planet rotation rate, when the satellite is located at the maximum or minimum latitude. It is apparent that all these orbits have periods greater than one nodal day. The quasi-synchronous orbit associated with N t = 1 and m = 2 ( q = 0.5 ) is selected, and corresponds to R = 32427 km and i = 60.0 • . This choice corresponds to a supersynchronous orbit, and is motivated by the reasonable values of its inclination and Fig. 2 Quasi-synchronous circular orbits: radius as a function of inclination and related number of orbits per nodal day radius, together with the short repetition period, i.e. 2 nodal days. Figure 3 depicts the ground tracks of three satellites placed along 3 distinct orbits, with the preceding values of R and i and right ascension separation of 120°. It is worth noticing that each ground track has cusps at the extremal latitudes. This is a typical feature of orbits that are synchronous with the planet rotation at specific points, corresponding to the cusps of the ground track.

Coverage Analysis for 9 Satellites in Quasi-synchronous Orbits
Overall, 9 satellites are placed in the inclined quasi-synchronous circular orbits described in the preceding subsection. Each orbit is traveled by 3 satellites, with angular separation of 120°. As a result, any point along the ground track is flown over by satellite 2 and 3, with delays equal to T n 3 and 2T n 3 with respect to satellite 1.
While coverage of the 3 areostationary satellites is identified by the regions portrayed in Fig. 1, the analysis for the remaining 9 satellites is more complicated. In a repetition period, for a given location on the Martian surface, the minimum elevation of visibility is to be determined. To do this, the scheduling of the visible passes of all the 9 satellites must be obtained and inspected. However, due to the symmetry properties of the ground track, identifying the minimum elevation angle of a limited number of locations suffices to supply precise information on the global coverage properties. These sample locations are portrayed in Fig. 3, and lie along the symmetry meridians of the ground track. In fact, let min g , denote the minimum elevation angle over a repetition period, as a function of latitude, , and geographical longitude, g . If is set to a specific value , min g , becomes a function of a single variable, which is continuous and locally even with respect to six values of g ,  Table 1. It is apparent that the 9 satellites placed in inclined orbits can guarantee global and continuous coverage of the Martian surface with an elevation angle never lower than 20.1°. Figures 4 and 5 depict the elevation and azimuth time histories of all the satellites over two distinct sites on the planet surface. Inspection of these figures reveals that multiple satellites are visible during the entire repetition period.

Low-Thrust Orbit Dynamics
In recent years, low thrust propulsion has attracted a strong interest in the scientific community, and has already been used in several space missions, e.g. the NASA Deep Space 1 and the ESA Smart-1 missions [14,15], to name a few. Low-thrust propulsion leads to substantial propellant savings, at the price of increasing (even considerably) the time of flight. Recently, nonlinear control was proven to be an effective option for real-time feedback guidance in Earth orbit transfers, as well as for orbit maintenance [10,[16][17][18] This section is focused on low-thrust orbit dynamics, using the modified equinoctial elements to describe the dynamics of 4 carrier vehicles, modeled as point masses. Orbital motion of each carrier is mainly affected by the Martian gravitational field. Therefore, the spacecraft dynamics can be investigated appropriately by employing a perturbed two-body-problem model. The most relevant orbit perturbations are included in this dynamical framework, i.e. (i) some harmonics of the areopotential (i.e., J 2 , J 3 , and J 4 ) and (ii) third body gravitational pull due to the Sun. These 4 spacecraft depart from a specified parking orbit, and drive the 12 satellites that form the constellation toward their orbits. More specifically, each carrier vehicle contains 3 satellites, and delivers all of them in the respective orbit, either the aerostationary orbit or one of the 3 quasi-synchronous, inclined, circular orbits. Steerable, throttleable low-thrust propulsion actuates the orbit injection maneuvers of the 4 carriers. After delivery from the carrier vehicle, each satellite performs phasing maneuvers, using its own chemical propulsion system, with the final aim of reaching the correct argument of latitude along the operational orbit.
The spacecraft dynamics can be described using the osculating orbit elements, i.e. semimajor axis a, eccentricity e, inclination i, right ascension of the ascending node (RAAN) Ω , argument of periapse , and true anomaly f. However, the Gauss equations [19], which govern the time evolution of the orbit elements, become singular while approaching a circular or equatorial orbit and when the eccentricity tends to 1. For these reasons, the modified equinoctial elements (MEE) [19] l, m, n, s, and q are employed, together with the semilatus rectum (parameter) p, used in place of a. The five elements l, m, n, s, and q are defined as [19] These elements are nonsingular for all trajectories, with the only exception of equatorial retrograde orbits ( i = ). If ∶= 1 + l cos q + m sin q , the orbit (instantaneous) radius equals r = p∕ . The governing equations for the modified equinoctial elements are [17,19,20] where the Mars gravitational parameter and Vector a is a (3 × 1)-vector that collects the components of the non-Keplerian acceleration that affects the spacecraft motion. These are denoted with a r , a , a h and refer to the local vertical local horizontal (LVLH) rotating frame aligned with r,̂,ĥ , where unit vector r is directed toward the instantaneous position vector r (taken from the Mars center), whereas ĥ is aligned with the spacecraft angular momentum. Vector a includes both the thrust acceleration a T and the perturbing acceleration a P inherent to the space environment. Thus a = a T + a P . The perturbing acceleration (due to harmonics of the areopotential and third body gravitational attraction) are then projected in the LVLH-frame. It is worth remarking that the preceding equations of motion reduce to ̇x 1 = … =̇x 5 = 0 and ̇x 6 3 1 if a = a T + a P = 0 (Keplerian motion). This means that 5 out of 6 MEE remain constant along Keplerian paths, similarly to what occurs if classical orbit elements are used.
Let T max and m 0 denote the maximum available thrust magnitude and the initial mass. If x 7 represents the mass ratio and T the thrust magnitude, the following equation holds for x 7 : where c represents the effective exhaust velocity of the propulsion system. The magnitude of the instantaneous thrust acceleration is a T = u T m 0 m = u T ∕x 7 and is con- In conclusion, the spacecraft dynamics is described in terms of the state vector

Nonlinear Orbit Control for Low-Thrust Deployment
In this section, a nonlinear technique is described, aimed at identifying the thrust direction and magnitude for a carrier spacecraft tailored to constellation deployment and equipped with a low-thrust propulsion system. The initial parking orbit is a typical ESA 4-sol orbit [11], associated with the following orbit elements: At the initial epoch t 0 , set to 19 April 2023 at 0:00 UTC, the true anomaly is assumed to equal 180°.

Deployment Conditions
Each carrier must deploy 3 satellites in the respective desired orbit, identified by semimajor axis, eccentricity, inclination, and RAAN. On average, only the latter element changes in time due to J 2 (cf. Eq. (1)). This means that for satellite k the desired RAAN The remaining conditions regard the (time-varying) orbit plane orientation, defined by i d and Ω (k) d . Omitting the superscript (k) hence forward, these conditions are enforced by requiring that the instantaneous angular momentum, with direction aligned with ĥ , point toward the desired direction, denoted with ĥ d . Unit vector ĥ can be written in terms of Ω and i, i.e. ĥ =ĉ 1 sin Ω sin i +ĉ 2 cos Ω sin i +ĉ 3 cos i . [13] The alignment condition is After some steps, omitted for the sake of brevity, the preceding relation leads to the following equality: The final conditions for deployment (13) and (15) can be incorporated as the three components 1 , 2 , 3 of vector and finally expressed in compact form as It is worth remarking that is time-varying due to Ω (k) d . For the 3 inclined orbits, the initial values of Ω k t 0 are

Feedback Law and Stability Analysis
In a preceding section, the spacecraft motion was shown to be governed by Eqs. (8)- (10). In particular, Eq. (8) can be rewritten as where the perturbing acceleration a P includes several contributions, related to the space environment, i.e. harmonics of the areopotential and third body gravitational attraction due to the Sun. In this research, orbit perturbations are assumed to be perfectly known, for the purpose of applying the feedback laws described in this section. For systems governed by Eq. (18) with a P = 0 , the Jurdjevic-Quinn theorem [21] provides a feedback control law that drives the dynamical system to an arbitrary target state, making the controlled system Lyapunov-stable. This subsection outlines the fundamental steps needed to identify a saturated feedback control law for lowthrust orbit injection of the four carrier vehicles. Further details and proofs are contained in a recent publication [10], which adopts a similar approach for a completely different constellation about Earth.
The three (scalar) deployment conditions are written in the compact form (16), which defines the target set. Due to Eqs. (13) and (15), this is a connected and differentiable manifold. Nonlinear orbit control is aimed at defining a feedback control law capable of driving the dynamical system at hand (associated with Eqs. (9), (11), and (18)) toward the target conditions identified by Eq. (16). To do this, the following candidate Lyapunov function is introduced: The Journal of the Astronautical Sciences (2022) 69:1691-1725 where K denotes a diagonal matrix with constant, positive elements, which play the role of arbitrary weights. These are selected a priori. It is obvious that V > 0 unless = 0 . Yet, further conditions are required in order that V be an actual Lyapunov function [22]. As a preliminary steps, the following auxiliary vectors are introduced: Two propositions establish the conditions for V to be a Lyapunov function, and lead to identifying a saturated feedback law. Proof of Proposition 1 is reported in Ref. [10]. This proposition includes the assumption u (max) T . In this case, in place of (21), an alternative feedback law can be used. Proof of Proposition 2 is reported in Ref. [10]. Actually, Eqs. (21) and (22) provide a straightforward indication on the value of u T that minimizes V . In fact, selection of u T = u (max) T corresponds to the least value of V . In short, due to the last consideration, the two feedback laws (21) and (22) can be written in compact form as Equation (23) incorporates the saturation condition on u T , i.e. | | u T | | ≤ u (max) T , and provides a control law that can be actuated using steerable and throttleable propulsive thrust (with time-varying magnitude and direction).
It is worth remarking that Propositions 1 and 2 provide some sufficient conditions for stabilizing the dynamical system of interest. This circumstance implies that the assumptions of Propositions 1 and 2 can be violated (in some time intervals), without necessarily compromising asymptotic convergence to the desired final condition identified by Eq. (16).
Because the attracting set contains other subsets other than the target set, the asymptotic convergence toward the desired conditions is only local, based on Lyapunov's stability theorem. However, the equality b = 0 can be considered again, to rule out, if possible, subsets 1 and 2. The condition b = 0 implies ̇b = 0 , i.e. ̇b 1 =̇b 2 =̇b 3 = 0 , where Equation (24) yields three relations. Inspection of their closed-form expressions [10] leads to ruling out subset 1. Therefore, only subsets 2 and 3 correspond to equilibria. Actually, the possible convergence toward subset 2 is only theoretical. In fact, using x 2 2 + x 2 3 = e 2 in the expression of 2 , the Lyapunov function can be rewritten as where e is the instantaneous eccentricity. Because 1 and 3 are independent of x 2 and x 3 , the partial derivative of V with respect to e is V∕ e = 2k 2 e e 2 − e 2 d . It is apparent that V∕ e = 0 (i.e. V is stationary) at e = 0 , which is consistent with the fact that subset 2 corresponds to an equilibrium condition. However, if e = e > 0 (with e arbitrarily small), then V∕ e < 0 , and the reduction of V leads e to increasing up to the desired value e d , which is associated with subset 3, i.e. the target set. This means that the latter corresponds to a stable equilibrium, unlike subset 2. This circumstance has the very interesting practical consequence that-from the numerical point of view-the dynamical system of interest enjoys global convergence The Journal of the Astronautical Sciences (2022) 69:1691-1725 toward the desired operational conditions, provided that the control law (23) is adopted, while holding the assumptions of Propositions 1 or 2. As a final step, vector ( ∕ z) −1 ( ∕ t) , which appears in the definition of d, can be obtained for the deployment problem at hand. Its first two components are identically 0, whereas the third component is proven to be finite almost everywhere [10]. This circumstance, in conjunction with Eq. (23) (that includes the case |b + d| > u (max) T ), ensures that the feedback law (23) is feasible.

Low-Thrust Constellation Deployment
Electric propulsion requires onboard electrical power in order to operate. In some cases, this can be provided only when the space vehicle is illuminated. However, in this subsection, electric propulsion is assumed to be able to operate regardless of the spacecraft lighting conditions. On the other hand, it is worth remarking the carrier vehicles encounter eclipse conditions only in the last portion of the transfer path, therefore the preceding assumption is reasonable. The overall propulsive performance is identified by the following parameters: The numerical simulations are performed using canonical units. The distance unit (DU) equals the Martian radius, whereas the time unit (TU) is such that The reference epoch is set to April 19, 2023. Specific tolerances are allowed on the desired final conditions, associated with Eqs. (13) and (15), and this allows avoiding an excessive propellant expenditure. In particular, if the three constraints are met to a prescribed accuracy, then the respective weight is set to zero, i.e.
For each satellite, orbit injection is achieved when the three conditions of Eq. (26) are met. Four distinct final orbits are considered. For each of them, numerical search on the weighting coefficients leads to identifying the respective values that minimize the time of flight (cf. Appendix).
In the numerical simulations, the following perturbations are modeled: (i) some harmonics of the areopotential (i.e., J 2 , J 3 , and J 4 ) and (ii) third body gravitational pull due to the Sun. As a preliminary step, the inequality u (max) T > |d| is checked. Using the spacecraft data (i.e. propulsion and mass), the minimum available thrust acceleration u (max) T turns out to exceed the maximal magnitude of the perturbing acceleration. As a result, the feedback law guarantees the convergence toward the desired operational conditions, associated with the four final orbits.

Nominal Initial Conditions
This subsection considers the low-thrust orbit transfers using nonlinear orbit control, leading the carrier spacecraft to inject into the four operational orbits, i.e. the areostationary orbit and the 3 distinct quasi-synchronous, inclined orbits. Subsequent phasing is performed by each satellite that forms the constellation (cf. Section 6). The initial (nominal) orbit is a typical ESA 4-sol orbit [11], associated  6, 7, and 8 portray the transfer paths leading to delivering the carrier vehicles in quasi-synchronous, inclined orbits 1, 2, and 3, together with the corresponding time histories of SMA, inclination, eccentricity, and RAAN. It is worth remarking that the latter has a desired value that is time-varying. Figure 9 depicts the transfer trajectory toward the areostationary orbit, accompanied by

Nonnominal Initial Conditions
Nonlinear orbit control appears as particularly appropriate as a feedback strategy for autonomous guidance when nonnominal flight conditions occur. This section considers the problem of driving the carrier spacecraft toward the respective operational orbits in the presence of large dispersions on the initial orbit elements. These displacements are simulated stochastically in a Monte Carlo campaign composed of 100 simulations, by assuming Gaussian distribution for all the orbit elements, with mean value equal to the nominal value of the 4-sol orbit (cf. Eq. (12)) and standard deviations reported in Table 2. However, to avoid unrealistic values and hyperbolic trajectories, the limiting bounds reported in Table 2 are assumed for each orbit element.  Table 3 reports the average values (denoted with the upper bar) and the standard deviation (superscript ( ) ) for the time of flight and the final mass ratio, as well as the correlation coefficients between either the time of flight and the initial inclination ( i ) or the time of flight and the initial RAAN ( Ω ). Figures 10,11,12,13,15,16,17,18,20,21,22, and 23 portray the time histories of the semimajor axis, eccentricity, inclination, and RAAN for the 3 quasi-synchronous inclined orbits. Figures 14, 19, and 24 illustrate the time of flight as a function of the initial inclination and RAAN, for the 3 orbits of interest. It is apparent that a correlation Ω exists for deployment into 2 quasi-synchronous orbits out of 3. Instead, the transfer paths leading to quasi-synchronous orbit 1 exhibit weak correlation Ω . This is due to the time evolution of the RAAN, shown in Fig. 13. In fact, in this case two subsets of transfer trajectories exist: in subset 1, the RAAN increases, whereas in subset 2 the RAAN reduces, with convergence toward the final desired value in both cases. For this orbit, correlation i is greater than for the remaining 2 quasi-synchronous orbits. Figures 25, 26, 27, and 28 depict the time histories of the semimajor axis, eccentricity, and inclination referred to the final areostationary orbit, as well as the time of flight as a function of the initial inclination and RAAN. In this case, a correlation emerges between the time of flight and the initial inclination, testified by the value of i . This means that the transfer time needed to converge to the final areostationary (equatorial) orbit is substantially affected by the initial inclination. In all cases, in spite of largely dispersed initial conditions, the time of flight has an average value close to the nominal one and standard deviations between 5.5 and 13.5 days. The final mass ratio has average values ranging from 0.931 (quasi-synchronous orbit 2) to 0.950

Orbit Phasing
The first phase of constellation deployment ends with orbit injection of each carrier, which corresponds to achieving the desired orbit plane, with zero eccentricity and the prescribed orbit radius. Then, the satellites are released, and each of them performs corrective phasing maneuvers, with the use of chemical propulsion. Usually, this implies very short ignitions, and the impulsive thrust approximation can be adopted as a result [23,24]. The methodology for phasing described in this section resembles the one described in Ref. [10] for Earth orbits. Orbit phasing is performed by placing the space vehicle in a parking elliptic orbit. Two options are investigated: (1) internal phasing: a braking horizontal velocity change injects the spacecraft into an inner elliptic orbit; after a certain phasing time, a second velocity variation (at apoapse) inserts the satellite into the final orbit; (2) external phasing: an accelerating horizontal velocity change injects the spacecraft into an outer elliptic orbit; after a certain phasing time, a second velocity variation (at periapse) inserts the satellite into the final orbit.
In both cases the two velocity changes have the same magnitude. The initial and final orbit are identical, but the spacecraft is shifted along the orbit, to reach the desired position, as a result of the phasing maneuver. In the context of orbit phasing, harmonic J 2 of the areopotential (related to the Martian oblateness) is the only perturbing action included in the model, because it is largely dominant if compared to the remaining orbit perturbations. After the first velocity variation Δv̂ (either in case (a) or in case (b)), the parking elliptic orbit has semimajor axis a PH and eccentricity e PH . The satellite completes n PH parking orbits (from either apoapse to apoapse, case (1), or periapse to periapse, case (2)) in a time interval where ̇f represents the average time derivative of the true anomaly [25]. Along the parking orbit, the RAAN is subject to natural drift due to Mars oblateness (harmonic J 2 ) that differs from the desired one. The displacement ΔΩ PH between the desired change of RAAN and the actual variation is where Δt PH is given by Eq. (29). To compensate this undesired drift, an out-of-plane velocity change Δv hĥ is applied when the argument of latitude equals 90°, to avoid any alteration of the orbit inclination i. [19] To do this, the out-of-plane velocity change has component However, the out-of-plane velocity variation has the undesired effect of altering the argument of latitude along the circular orbit [19], To address this issue, a virtual vehicle is introduced, with argument of latitude d corresponding to the desired angular position; d,0 denotes the value of d when the phasing maneuver begins. Instead, the actual argument of latitude along the parking orbit is , with initial value 0 (when the first velocity change is applied). Let Δ ∶= d,0 − 0 denote the phasing angle, i.e. the difference between the desired and the actual argument of latitude when the phasing maneuver starts. After two inplane velocity changes and the final out-of-plane velocity variation (cf. Eq. (32)) the final argument of latitude equals Equation (33) includes also the variation due to the out-of-plane velocity change Δ (cf. Equation (32)). In the last relation the average time derivatives ̇ and ̇f are given by Eqs. (2) and (29), respectively. The condition that ensures correct phasing is with either k = 0 (internal phasing) or k = 1 (external phasing). The symbol ̇d represents the average time derivative of the argument of latitude of the virtual vehicle. Using the definition of Δ ∶= d,0 − 0 , the phasing Eq. (34) becomes In short, three velocity changes (two in-plane changes and a single out-of-plane velocity variation) allow performing the phasing maneuver. For a specified n PH and a specified (desired) angular displacement Δ , the following steps lead to identifying Δv and Δv h : In the end, the preceding steps allow placing the satellite at the correct angular position, while taking into account the coupled effects of in-plane and out-of-plane velocity changes. Figure 29 portrays the overall velocity change needed for phasing as a function of Δ , along the quasi-synchronous inclined orbit, whereas Fig. 30 depicts the corresponding in-plane and out-of-plane velocity variations. Different numbers of orbits (i.e. 10, 20, 30, 40, and 50) along the parking path are considered, and correspond to the durations reported in the right inset. Only the most convenient option between internal and external phasing is depicted. It is apparent that for values of Δ larger than 180°, external phasing is more convenient in terms of overall velocity change. Conversely, internal transfers are preferable if the satellite must be moved ahead of its initial position by an angle that does not exceed 180°. Not surprisingly, larger durations correspond to reduced propellant budgets. In particular, Fig. 29 points out that the overall velocity change reduces to modest values if the time spent in the parking orbit exceeds 80 days. As a further remark, the comparison of in-plane and out-of-plane velocity changes reveals that the latter have negligible values if compared to in-plane velocity variations. Moreover, the out-of-plane velocity change turns out to be relatively independent of the phasing duration. Figure 31 depicts the overall velocity change as a function of Δ , along the areostationary orbit. It is worth remarking that in this case no out-of-plane velocity change is required. Also in this case, longer durations correspond to reduced velocity variations needed to complete orbit phasing, with modest values if the time spent in the parking orbit exceeds 50 days.

Concluding Remarks
In this research, a satellite constellation that guarantees global continuous coverage over Mars is designed. This could serve as a communications relay system and would certainly represent a valuable asset for future missions. While 3 satellites travel an areostationary orbit, the remaining 9 satellites are released in three repeating-ground-track, circular, quasi-synchronous, inclined orbits. This configuration guarantees global and continuous coverage, while ensuring repetition and predictability of the visible passes and simultaneous visibility of multiple satellites. Moreover, this work addresses the problem of deploying the 12 satellites using 4 carrier vehicles, released by the mother spacecraft. The latter orbits a highly elliptical path, i.e. the ESA 4-sol orbit, which is entered after planetary capture. Orbit dynamics about Mars is modeled with the inclusion of the most relevant perturbations, i.e. several harmonics of the areopotential, together with the gravitational pull due to  the Sun as a third body. Each carrier embarks 3 satellites and employs steerable and throttleable low-thrust propulsion for the purpose of driving them toward the respective operational orbit. A saturated feedback law for the low-thrust magnitude and direction-with an upper bound on magnitude-is described and used, and enjoys global stability properties. In particular, some sufficient conditions for the asymptotic convergence toward the operational conditions are stated. These conditions include the perturbing acceleration term, and lead to identifying three types of transfer arcs: (a) maximum-thrust arcs, (b) time-varying, intermediate-thrust arcs, and (c) coast arcs. The numerical simulations demonstrate that the low-thrust deployment strategy based on nonlinear orbit control is effective, and needs a limited amount of propellant. Moreover, the feedback control law at hand does not require any offline reference trajectory. Therefore, it is effective as an autonomous real-time guidance strategy, even in the presence of nonnominal flight conditions. With this regard, a Monte Carlo campaign proves that orbit injection is successfully completed even starting from dispersed initial conditions, which can be regarded as representative of initial elliptical orbits significantly different from the expected one, due to injection errors at the planetary capture. Orbit phasing represents the last phase of deployment, and is performed by each satellite. Both internal and external phasing are considered. A simple and effective strategy is described that includes out-of-plane and in-plane velocity changes and avoids leaving the operational orbit plane. The tradeoff between overall propellant expenditure and phasing time is identified. and k 3 . Overall, 16 simulations were run for each final orbit, under the assumption of neglecting orbit perturbations, with the results collected in Tables 4, 5 Tables 4, 5 , 6, 7, 8, 9, 10, and 11 leads to choosing the most appropriate gains for each case, i.e., • areostationary orbit: k 1 = 1, k 2 = 10 4 , k 3 = 10 6 • quasi-synchronous orbit 1: k 1 = 1, k 2 = 10 4 , k 3 = 10 6 • quasi-synchronous orbit 2: k 1 = 1, k 2 = 10 4 , k 3 = 10 5 • quasi-synchronous orbit 3: k 1 = 1, k 2 = 10 4 , k 3 = 10 6 For the case of final areostationary orbit, the result attained after gain tuning was compared to the minimum-time orbit transfer, detected with the use of the indirect heuristic method [26]. It is worth remarking that the latter approach is computationally intensive and cannot be performed in real time. Moreover, optimal control does not provide a feedback law capable of compensating for nonnominal flight conditions and orbit perturbations. Nevertheless, the minimum-time transfer is useful for the purpose of comparing the numerical results found with nonlinear orbit control and optimal control. For the case at hand (final areostationary orbit), the following results were obtained: • optimal control: t f = 33.4 days, x 7,f = 0.953 • nonlinear feedback control: t f = 39.1 days, x 7,f = 0.946 These results provide a clear indication on the capability of approaching the performance obtained with optimal control, through proper selection of k 1 , k 2 , and k 3 .