Open AccessArticle

Self-Organizing Control of Mega Constellations for Continuous Earth Observation

Yun Xu

^1,2,*

Yulin Zhang

^1,2

Zhaokui Wang

³,

Yunhan He

^1,2 and

Li Fan

^1,2

College of Control Science and Engineering, Zhejiang University, Hangzhou 310027, China

Huzhou Institute of Zhejiang University, Huzhou 313000, China

School of Aerospace Engineering, Tsinghua University, Beijing 100084, China

Author to whom correspondence should be addressed.

Remote Sens. 2022, 14(22), 5896; https://doi.org/10.3390/rs14225896

Submission received: 22 October 2022 / Revised: 16 November 2022 / Accepted: 16 November 2022 / Published: 21 November 2022

(This article belongs to the Special Issue CubeSats Applications and Technology)

Download

Browse Figures

Graphical abstract
"> Figure 1
Satellite coverage geometry. "> Figure 2
Street-of-coverage geometry of one orbital plane. "> Figure 3
Streets-of-coverage geometries of adjacent orbital planes. "> Figure 4
Street-of-coverage geometry of the perturbed coplanar satellites. D and E are two sub-satellite points. A and C are two intersections of two coverage circle. Point A, B, C form a right triangle. Point D, E, F form a right triangle. "> Figure 5
Illustration of the scale-independent relative orbital elements (<math display="inline"><semantics> <mrow> <mi>x</mi> <mo>−</mo> <mi>y</mi> </mrow> </semantics></math> plane). "> Figure 6
Illustration of the scale-independent relative orbital elements (<math display="inline"><semantics> <mrow> <mi>x</mi> <mo>−</mo> <mi>z</mi> </mrow> </semantics></math> plane). "> Figure 7
Relative motion-bound changing with time. "> Figure 8
Semi-major axis time histories without semi-major axis control. "> Figure 9
Semi-major axis time histories under semi-major axis control. "> Figure 10
Control effects changing with time. "> Figure 11
The initial distribution of the coplanar satellites in the geocentric inertial system. Red dots represent satellites that do not satisfy the continuous coverage constraint while green dots represent satellites that satisfy the continuous coverage constraint. "> Figure 12
The final distribution of the coplanar satellites in the geocentric inertial system. "> Figure 13
Number of unqualified satellites in the coplanar satellites. "> Figure 14
The differences in maximum bounds away from the corresponding constraints. Red line is zero-value standard. Blue lines are maximum bound differences. "> Figure 15
The differences in minimum bounds away from the corresponding constraints. Red line is zero-value standard. Blue lines are minimum bound differences. "> Figure 16
Semi-major axis differences away from the target value. "> Figure 17
RAAN differences away from the target RAAN under <math display="inline"><semantics> <msub> <mi>J</mi> <mn>2</mn> </msub> </semantics></math> perturbation. Red lines represent <math display="inline"><semantics> <mrow> <mi>δ</mi> <msub> <mi mathvariant="sans-serif">Ω</mi> <mi>min</mi> </msub> </mrow> </semantics></math> and <math display="inline"><semantics> <mrow> <mi>δ</mi> <msub> <mi mathvariant="sans-serif">Ω</mi> <mi>max</mi> </msub> </mrow> </semantics></math> respectively. Blue lines represent RAAN differences. "> Figure 18
Control effects of coplanar satellites in the geocentric inertial system. "> Figure 19
Initial distribution of the mega constellation. Red dots represent satellites that do not satisfy the continuous coverage constraint while green dots represent satellites that satisfy the continuous coverage constraint. "> Figure 20
Final distribution of the mega constellation. "> Figure 21
Number of unqualified satellites in the mega constellation. "> Figure 22
Control effects of constellation satellites. "> Figure 23
Maximum bound differences of constellation satellites. Red line is zero-value standard. Blue lines are maximum bound differences. "> Figure 24
Minimum bound differences of constellation satellites. Red line is zero-value standard. Blue lines are minimum bound differences. "> Figure 25
Semi-major axis differences of constellation satellites. "> Figure 26
RAAN differences of constellation satellites. Red lines represent <math display="inline"><semantics> <mrow> <mi>δ</mi> <msub> <mi mathvariant="sans-serif">Ω</mi> <mi>min</mi> </msub> </mrow> </semantics></math> and <math display="inline"><semantics> <mrow> <mi>δ</mi> <msub> <mi mathvariant="sans-serif">Ω</mi> <mi>max</mi> </msub> </mrow> </semantics></math> respectively. Blue lines represent RAAN differences. "> Figure 27
Coverage performance of the GW-2 sub-constellation. "> Figure 28
Multiplicity of coverage of the GW-2 sub-constellation at the initial and final moments. ">

Review Reports Versions Notes

Abstract

This work presents a novel self-organizing control method for mega constellations to meet the continuous Earth observation requirements. In order to decrease the TT&C pressure caused by numerous satellites, constellation satellites are not controlled according to the designed configurations but are controlled with respect to intersatellite constraints. By analyzing the street-of-coverage (SOC) of coplanar constellation satellites, the continuous coverage constraint of the mega constellation is transformed into constraints of the right ascension of ascending node (RAAN) and relative motion bound between every two adjacent coplanar satellites. The proposed continuous coverage constraint can be satisfied by most ongoing or planned mega constellations. Artificial potential functions (APFs) are used to realize self-organizing control. The scale-independent relative orbital elements (SIROEs) are innovatively presented as the self-organizing control variables. Using the Gaussian equations and Lyapunov’s theory, the stability of the APF control in quadratic form is proven, from which it can be concluded that the APF control variables of the controlled satellite should have the same time derivative as the target satellite states under two-body Keplerian motion condition, and SIROEs are ideal choices. The proposed controllers and self-organizing rules are verified in the sub-constellation of the GW-2 mega constellation by simulation. The results demonstrate the goodness in control effect and ground coverage performance.

Keywords:

mega constellation; continuous Earth observation; self-organizing control; artificial potential function

Graphical Abstract

1. Introduction

In recent years, mega constellations consisting of thousands of satellites have received widespread attention due to their potential to provide high-speed and full-coverage telecommunication services. Commercial companies, such as SpaceX and OneWeb, have made their plans public to deploy mega constellations in low Earth orbits [1]. SpaceX have launched nearly 2000 Starlink satellites and started to provide global communication and Internet services. China plans to construct two Guowang low-Earth-orbit constellations, totaling 12,992 satellites for telecommunication services [2].

Such giant constellations are wonderful choices for performing Earth observation missions. Since they evenly distribute at low Earth orbits all over the globe, they can be used to collect huge amounts of data on gravity [3] and atmosphere [4] for analysis, and they can obtain high–resolution Earth images in real time [5]. Among these remote sensing missions, the priority is to provide continuous global coverage. Limited by the computation load of analytic and grid methods, traditional coverage analysis methods are suitable for constellations with less than 100 satellites [6,7,8]. For the requirement of quick analysis of mega constellations, Dai et al. studied the geometric configuration of constellation projection points on Earth’s surface [9]. The multi-satellite coverage problem is decomposed into a one-satellite coverage problem to solve both continuous and discontinuous coverage problems. Gong et al. proposed a continuous coverage analysis method based on 2D map theory [10]. The proposed method avoids redundant calculation of intersatellite positions and positions between the satellites and ground targets so that the time complexity does not increase significantly with the increase in satellites. Huang et al. presented a systematic method for the design of continuous global coverage Walker and street-of-coverage constellations, by taking seven critical constellation properties (coverage, robustness, self-induced collision avoidance, launch, build-up, station-keeping, and end-of-life disposal) as design criteria [11]. However, the above studies rely on the overall motion information of the constellation, which will result in central ground control.

To that end, this paper innovatively introduces the self-organizing concept of mega constellation control, and the emphasis is on maintaining continuous coverage. Two traditional constellation maintenance methods are the absolute station-keeping method and the relative station-keeping method [12,13,14]. For mega constellations, the former will result in frequent on-orbit maneuvers, while the latter relies on the overall configuration information, which is also inefficient because of the huge number of satellites. In this paper, self-organizing control indicates a decentralized and distributed control mechanism that requires the relative motion information of neighbor satellites only. By that means, constellation control will not be affected by the failure or development of constellation satellites. Artificial potential function (APF) is a commonly used self-organizing control method for its clear mathematical descriptions. Trajectories based on APF control are not planned explicitly but rather are shaped by the potential field in reaction to the changing environment. APF control has shown its potential in many close proximity scenarios, including satellite cluster reconfiguration [15,16], collision avoidance [17], satellite rendezvous and docking [18,19], etc. Spencer improved the APF from Cartesian coordinate form to relative orbital elements (ROEs) form so that the relative orbital dynamics can be reflected in APF control [20]. Our previous work realized reconfiguration, bounded flight, and collision avoidance for satellite clusters by means of ROE-based APF control [21,22]. However, the relative motion distance between constellation satellites is much larger than that of close clusters. The above control methods cannot be applied to mega constellation scenarios.

For the requirement of continuous Earth observation of mega constellations, this paper analyzes the relationship between orbital planes, coplanar adjacent satellites, and coverage performance and proposes a continuous novel coverage constraint that relates to intersatellite motion only. The new constraint can be applied to the self-organizing control of mega constellations. The APF method is adopted to solve the self-organizing control problem in this work. Considering the large range of relative motion between constellation satellites, the scale-independent relative orbital elements (SIROEs) are proposed as the APF control variables. SIROEs have been proposed in our preliminary work [23,24], which have clear geometric interpretations, concise expressions, and high precision for relative motion with near-circular reference orbit and arbitrary intersatellite separation, which are particularly suitable for mega constellation satellites. The stability of the proposed self-organizing control is proven based on the Gaussian equations and Lyapunov’s theory. An ideal APF control variable of the controlled satellite should have the same time derivative as the target satellite states under a two-body Keplerian motion circumstance. Finally, the GW-2 constellation is simulated to verify the proposed self-organizing control method. The results indicate its goodness of control effect and ground-coverage performance.

The continuous coverage constraint is analyzed in Section 2.1. The configuration requirements of mega constellations are given here. The artificial potential functions for a single satellite and the self-organizing control rules for coplanar satellites are studied in Section 2.2. The detailed stability proofs of the proposed artificial potential functions are given in Section 2.3. The performance of self-organizing control for continuous coverage is evaluated based on a GW-2 constellation in Section 3. Conclusions are provided in Section 4.

2. Materials and Methods

2.1. Continuous Coverage Constraint of Mega Constellation

Considering the numerous orbital planes and satellites in one mega constellation, the continuous coverage constraint is redefined in this section. The configuration requirements for mega constellations that satisfy the coverage constraint are proposed. The intersatellite motion constraints during the lifetime are also analyzed.

2.1.1. Definition of Continuous Coverage

Walker constellation has been applied in many mega constellation missions because of its simple design and regular configuration. The spatial configuration of the Walker constellation can be represented by the total number of satellites N, the number of orbital planes P, and the phase parameters F, where

2 π F / N

represents the phase difference between two satellites in adjacent orbits. All satellites are placed in circular orbits at the same altitude and at the same inclination.

The latitude range that can be covered by the constellation is determined by the designed inclination and orbital altitude, and the distribution density of constellation satellites at higher latitudes is much larger than that at lower latitudes. Hence, it is convincible that continuous coverage for the designed latitude range can be achieved as long as low latitude regions can be covered continuously.

The coverage geometry of a single satellite is shown in Figure 1, where the angular radius of the coverage range is

α_{max}

, the orbital altitude is H, the Earth’s radius is

R_{e}

, and the minimum elevation angle is

E_{min}

. The relationship between

α_{max}

, H, and

E_{\min}

is given by Equation (1).

α_{max} = arccos (\frac{R_{e}}{R_{e} + H} cos E_{min}) - E_{min}

(1)

The coverage range of adjacent coplanar satellites overlaps with each other, forming a street-of-coverage (SOC), shown in Figure 2. The width of the SOC

σ

depends on the coverage range

α_{max}

and the geocentric angle between adjacent coplanar satellites

θ

, written in Equation (2).

σ = arccos (\frac{cos 2 α_{max}}{cos θ})

(2)

Since mega Walker constellation satellites are located at the same inclination, the SOC of adjacent orbital planes are approximately parallel to each other at low latitudes, shown in Figure 3. If the SOC of adjacent orbital planes could overlap with each other, continuous coverage can be achieved at low latitudes. It is noted that adjacent orbital planes could be co-rotating planes or counter-rotating planes, where satellites in the co-rotating planes move in the same direction (

Δ Ω

is close to zero), while satellites in the counter-rotating planes move in the opposite direction (

Δ Ω

is close to

π

). Therefore, the relative motion relationships between co-rotating satellites are relatively fixed compared to counter-rotating satellites. From the perspective of self-organization, the adjacent orbital planes in the coverage constraint are specified as co-rotating planes.

To sum up, the continuous coverage constraint of the mega constellation is defined as follows: The continuous coverage can be guaranteed for a designed latitude range if the SOC of every two adjacent co-rotating orbital planes overlaps with each other. Theoretically, the existence of adjacent counter-rotating orbital planes can provide more coverage multiplicity.

2.1.2. Configuration Requirements for Mega Constellations

The mega constellation, containing a total of N satellites, is composed of S satellites evenly spaced on each of P orbital plane. The geocentric angle between adjacent coplanar satellites

θ

2 π / S

. According to Equation (2), the width of SOC of one orbital plane

σ

σ = arccos (\frac{cos 2 α_{max}}{cos (2 π / S)})

(3)

where

S \geq π / α_{max}

(4)

To ensure the SOC of every two adjacent co-rotating orbital planes overlaps with each other there exists

σ \geq 2 π / P

(5)

From Equations (3) and (5), it can be implied that

cos (\frac{2 π}{P}) cos (\frac{2 π}{S}) \geq cos (2 α_{max})

(6)

Therefore, the configuration requirements for mega constellations that satisfy the proposed continuous coverage constraint are given in Equations (4) and (6).

Several ongoing and planned mega constellation missions are listed in Table 1. Let the minimum elevation angle

E_{min}

takes 20

^{\circ}

. The coverage range

α_{\max}

can be calculated by Equation (1) together with the orbital altitude. Most of these mega constellations are proven to satisfy the configuration requirements for continuous coverage, i.e., Equations (4) and (6).

2.1.3. Intersatellite Motion Constraints

Due to perturbations and initial orbit errors, the coplanar satellites will drift from each other, and thus the RAAN differences between coplanar satellites will get larger over time without consideration of constellation control. The drift range of the RAAN difference between adjacent co-rotating orbital planes

δ Ω

, and the drift range of the geocentric angle between adjacent coplanar satellites

θ

under the continuous coverage constraint, are worth investigating.

The SOC geometry of the perturbed coplanar satellites is shown in Figure 4, where

∠ AC = arccos (\frac{cos 2 α_{max}}{cos θ})

(7)

∠ DE = θ

(8)

∠ EF = δ Ω

(9)

It can be calculated that the coverage range of two adjacent coplanar satellites along the track direction is

∠ DF = arccos (\frac{cos θ}{cos δ Ω})

(10)

The coverage range along the normal direction of the orbital plane is

∠ BC = \frac{sin θ cos 2 α_{max}}{\sqrt{{sin}^{2} θ {cos}^{2} θ - ({cos}^{2} θ - {cos}^{2} 2 α_{max}) {sin}^{2} δ Ω}}

(11)

To ensure the SOC of every two adjacent co-rotating orbital planes overlap with each other, there exists

∠ BC \geq \frac{2 π}{P}

(12)

To ensure the coverage range of every two adjacent coplanar satellites along the track direction overlap with each other, there exists

∠ D F_{max} \geq \frac{2 π}{S}

(13)

According to Equations (11) and (12), the drift range of

θ

can be obtained, which is related to the drift range of

δ Ω

θ_{min} \leq θ \leq θ_{max}

(14)

where

\begin{matrix} A_{1} = - {cos}^{2} (\frac{2 π}{P}) \\ B_{1} = - {cos}^{2} (\frac{2 π}{P}) (1 + {sin}^{2} δ Ω) + {cos}^{2} 2 α_{max} \\ C_{1} = {cos}^{2} (\frac{2 π}{P}) {sin}^{2} 2 α_{max} {sin}^{2} δ Ω \\ θ_{min} = arcsin (\sqrt{\frac{- B_{1} - \sqrt{B_{1}^{2} - 4 A_{1} C_{1}}}{2 A_{1}}}) \\ θ_{max} = arcsin (\sqrt{\frac{- B_{1} + \sqrt{B_{1}^{2} - 4 A_{1} C_{1}}}{2 A_{1}}}) \end{matrix}

(15)

According to Equations (10), (13) and (14), the drift range of

δ Ω

of every two adjacent satellites can be obtained as

0 \leq δ Ω \leq arcsin (\sqrt{\frac{B_{2}}{A_{2}}})

(16)

where

\begin{matrix} A_{2} = - 2 {sin}^{2} (\frac{2 π}{S}) {cos}^{2} (\frac{2 π}{S}) {cos}^{2} (\frac{2 π}{P}) - {cos}^{2} (\frac{2 π}{S}) {cos}^{2} 2 α_{max} + {cos}^{2} (\frac{2 π}{P}) {cos}^{2} 2 α_{max} \\ B_{2} = - {sin}^{2} (\frac{2 π}{S}) [{cos}^{2} (\frac{2 π}{S}) {cos}^{2} (\frac{2 π}{P}) - {cos}^{2} 2 α_{max}] \end{matrix}

(17)

The detailed derivation steps and the proofs of real roots are given in Appendix A.

We denote

Ω_{t}

as the target RAAN under

J_{2}

perturbation of the designed orbital plane. Then the drift range of

Ω_{i}

of all coplanar satellites away from

Ω_{t}

can be written as

\begin{matrix} δ Ω_{i} = Ω_{i} - Ω_{t} \geq - \frac{1}{2} arcsin (\sqrt{\frac{B_{2}}{A_{2}}}) = δ Ω_{min} \\ δ Ω_{i} = Ω_{i} - Ω_{t} \leq \frac{1}{2} arcsin (\sqrt{\frac{B_{2}}{A_{2}}}) = δ Ω_{max} \end{matrix}

(18)

We denote the position vectors of two adjacent coplanar satellites in the geocentric inertial frame as

r_{i}

and

r_{i + 1}

. The constraint of

θ

can be transformed to the relative distance constraint

d_{i, i + 1}

by means of cosine theorem.

\begin{matrix} d_{i, i + 1} \geq \sqrt{{∥r_{i}∥}^{2} + {∥r_{i + 1}∥}^{2} - 2 ∥r_{i}∥ ∥r_{i + 1}∥ cos θ_{min}} \\ d_{i, i + 1} \leq \sqrt{{∥r_{i}∥}^{2} + {∥r_{i + 1}∥}^{2} - 2 ∥r_{i}∥ ∥r_{i + 1}∥ cos θ_{max}} \end{matrix}

(19)

Considering that mega constellations are nearly circular orbits, the position vectors in the geocentric inertial frame are substituted by semi-major axis a, given in Equation (20).

\begin{matrix} d_{i, i + 1} \geq \sqrt{2 a^{2} (1 - cos θ_{min})} = d_{mint} \\ d_{i, i + 1} \leq \sqrt{2 a^{2} (1 - cos θ_{max})} = d_{maxt} \end{matrix}

(20)

To sum up, the intersatellite motion constraints that every two adjacent coplanar satellites should obey during their lifetime are presented in Equations (18) and (20).

2.2. Self-Organizing Control Using Artificial Potential Functions

In this section, the scale-independent relative orbital elements (SIROEs), as well as the corresponding bound formulations, are introduced first. The relative motion bound denotes the maximum or minimum value of the relative distance between two satellites. Based on SIROEs, the artificial potential functions for single satellite control are designed under the continuous coverage constraint. The self-organizing control rules for the coplanar satellites are presented.

2.2.1. Scale-Independent Relative Orbital Elements and Bound Formulations

Scale-independent relative orbital elements are high-precision configuration parameters for relative motion with near-circular reference orbit and arbitrary intersatellite separation. SIROEs have been proposed in our preliminary work [23]. Six independent SIROEs are

b_{1}, b_{2}, b_{3}, c, θ_{1}, θ_{2}

, shown in Figure 5 and Figure 6. The relationship between SIROEs and relative position is given in Equation (21). We denote

\tilde{x}, \tilde{y}, \tilde{z}

as the basic relative motion coordinates. The relative trajectory can be obtained by rotating the basic trajectory by

θ_{1}

, counterclockwise. The basic trajectory projected on the

x - y

plane is the combination of an ellipse (with a semi-major axis of

2 b_{1}

and a semi-minor axis of

b_{1}

), a circle (with a radius of

b_{2}

), and an additional offset (with x-offset of

b_{3} - a cos θ_{1}

and y-offset of

a sin θ_{1}

). The relative motion on the z-axis is harmonic.

\begin{matrix} [\begin{matrix} x \\ y \\ z \end{matrix}] & = [\begin{matrix} cos θ_{1} & - sin θ_{1} & 0 \\ sin θ_{1} & cos θ_{1} & 0 \\ 0 & 0 & 1 \end{matrix}] [\begin{matrix} - b_{1} cos (θ_{2} + ψ) + b_{2} cos 2 ψ + b_{3} - a cos θ_{1} \\ 2 b_{1} sin (θ_{2} + ψ) - b_{2} sin 2 ψ + a sin θ_{1} \\ c sin ψ \end{matrix}] \\ = [\begin{matrix} cos θ_{1} & - sin θ_{1} & 0 \\ sin θ_{1} & cos θ_{1} & 0 \\ 0 & 0 & 1 \end{matrix}] [\begin{matrix} \tilde{x} \\ \tilde{y} \\ \tilde{z} \end{matrix}] \end{matrix}

(21)

Without regard to the minor differences of semi-major axes between constellation satellites, the relative velocities represented by SIROEs are written as

[\begin{matrix} \dot{x} / n \\ \dot{y} / n \\ \dot{z} / n \end{matrix}] = [\begin{matrix} cos θ_{1} & - sin θ_{1} & 0 \\ sin θ_{1} & cos θ_{1} & 0 \\ 0 & 0 & 1 \end{matrix}] [\begin{matrix} b_{1} sin (θ_{2} + ψ) - 2 b_{2} sin 2 ψ \\ 2 b_{1} cos (θ_{2} + ψ) - 2 b_{2} cos 2 ψ \\ c cos ψ \end{matrix}]

(22)

where n is the mean angular velocity of the reference orbit.

According to Equations (21) and (22), SIROEs can be represented by translational states, given in Equation (23).

\begin{matrix} c = \sqrt{z^{2} + {(\dot{z} / n)}^{2}} \\ sin ψ = z / c, cos ψ = \dot{z} / n c \\ b_{2} = \frac{1}{2} (a^{'} - \sqrt{a^{' 2} - c^{2}}) \\ b_{3} = \frac{1}{2} (a^{'} + \sqrt{a^{' 2} - c^{2}}) \\ cos θ_{1} = - \frac{2 b_{3} (2 \dot{y} + n x + n a) + 3 b_{2} sin 2 ψ (2 n y - \dot{x})}{(2 \dot{x} - n y) (2 n y - \dot{x}) - (2 n x + 2 n a + \dot{y}) (2 \dot{y} + n x + n a)} \\ sin θ_{1} = \frac{2 n b_{3} - (2 n x + 2 n a + \dot{y}) cos θ_{1}}{2 n y - \dot{x}} \\ b_{1} = \sqrt{{[- (x + a) cos θ_{1} - y sin θ_{1} + b_{2} cos 2 ψ + b_{3}]}^{2} + {(\frac{\dot{x}}{n} cos θ_{1} + \frac{\dot{y}}{n} sin θ_{1} + 2 b_{2} sin 2 ψ)}^{2}} \\ sin β_{1} = \frac{\dot{x} cos θ_{1} + \dot{y} sin θ_{1} + 2 n b_{2} sin 2 ψ}{n b_{1}} \\ cos β_{1} = \frac{- (x + a) cos θ_{1} - y sin θ_{1} + b_{2} cos 2 ψ + b_{3}}{b_{1}} \\ sin θ_{2} = sin β_{1} cos ψ - cos β_{1} sin ψ \\ cos θ_{2} = cos β_{1} cos ψ + sin β_{1} sin ψ \end{matrix}

(23)

According to Equation (21), the squared relative distance between satellites can be written as

w = x^{2} + y^{2} + z^{2} = \frac{A_{3} s^{6} + B_{3} s^{5} + C_{3} s^{4} + D_{3} s^{3} + E_{3} s^{2} + F_{3} s + G_{3}}{{(1 + s^{2})}^{3}}

(24)

where

s = tan \frac{ψ}{2}, \begin{matrix} ψ \in [0, 2 π) \end{matrix}

(25)

Parameters

A_{3}, B_{3}, C_{3}, D_{3}, E_{3}, F_{3}, G_{3}

are polynomials about SIROEs

b_{1}, b_{2}, b_{3}, c, θ_{1}, θ_{2}

. Therefore,

A_{3}, B_{3}, C_{3}, D_{3}, E_{3}, F_{3}, G_{3}

are time-invariant, while s is the only time variable. Differentiating w with respect to s results in

\frac{\partial w}{\partial s} = \frac{\begin{matrix} - B_{3} s^{6} + (6 A_{3} - 2 C_{3}) s^{5} + (5 B_{3} - 3 D_{3}) s^{4} + (4 C_{3} - 4 E_{3}) s^{3} + (3 D_{3} - 5 F_{3}) s^{2} \\ + (2 E_{3} - 6 G_{3}) s^{3} + F_{3} \end{matrix}}{{(1 + s^{2})}^{4}}

(26)

Let

\partial w / \partial s = 0

, there exists

\begin{matrix} - B_{3} s^{6} + (6 A_{3} - 2 C_{3}) s^{5} + (5 B_{3} - 3 D_{3}) s^{4} + (4 C_{3} - 4 E_{3}) s^{3} + (3 D_{3} - 5 F_{3}) s^{2} \\ + (2 E_{3} - 6 G_{3}) s^{3} + F_{3} = 0 \end{matrix}

(27)

s = \pm \infty

(28)

The real roots of Equation (27) together with

s = \pm \infty

constitute a set of s values that perhaps make the squared relative distance take the maximum or minimum value, which are relative motion bounds. Note that when

s = \pm \infty

ψ = \pm π

, the squared relative distance takes

w = A_{3}

Assuming that

s_{max}

is the s value that makes the squared relative distance w take the maximum value and

s_{min}

is the one that makes w take the minimum value, the maximum bound

{\overset{⌢}{d}}_{i j}

and the minimum bound

{\overset{⌣}{d}}_{i j}

for the current time can be written as

\{\begin{matrix} {\overset{⌢}{d}}_{i j} = \sqrt{\frac{A_{3} s_{max}^{6} + B_{3} s_{max}^{5} + C_{3} s_{max}^{4} + D_{3} s_{max}^{3} + E_{3} s_{max}^{2} + F_{3} s_{max} + G_{3}}{{(1 + s_{max}^{2})}^{3}}}, if s_{max} \neq \pm \infty \\ {\overset{⌢}{d}}_{i j} = \sqrt{A_{3}}, if s_{max} = \pm \infty \end{matrix}

(29)

\{\begin{matrix} {\overset{⌣}{d}}_{i j} = \sqrt{\frac{A_{3} s_{min}^{6} + B_{3} s_{min}^{5} + C_{3} s_{min}^{4} + D_{3} s_{min}^{3} + E_{3} s_{min}^{2} + F_{3} s_{min} + G_{3}}{{(1 + s_{min}^{2})}^{3}}}, if s_{min} \neq \pm \infty \\ {\overset{⌣}{d}}_{i j} = \sqrt{A_{3}}, if s_{min} = \pm \infty \end{matrix}

(30)

Relative motion bounds are complicated functions of SIROEs, and SIROEs are functions about translational states, given in Equation (23). Therefore,

{\overset{⌢}{d}}_{i j}

and

{\overset{⌣}{d}}_{i j}

are functions about translational states

x, y, z, \dot{x}, \dot{y}, \dot{z}

2.2.2. Design of Artificial Potential Functions for Single Satellite Control

Referring to the concept of relative motion bound, Equation (20) can be rewritten as

{\overset{⌢}{d}}_{i j} \leq d_{maxt}

{\overset{⌣}{d}}_{i j} \geq d_{mint}

. Together with Equation (18), the proposed continuous coverage constraint relates to RAAN differences and relative motion bounds. Therefore, the control effects will be performed in cases if

δ Ω_{i} > δ Ω_{max}

δ Ω_{i} < δ Ω_{min}

{\overset{⌢}{d}}_{i j} > d_{maxt}

{\overset{⌣}{d}}_{i j} < d_{mint}

The following artificial potential functions are all designed in quadratic form, given in Equation (31).

ϕ_{i} = {(X_{i} - X_{t})}^{T} K (X_{i} - X_{t})

(31)

where

X_{i}

is the control variable.

X_{t}

is the expected value of

X_{i}

K

is the control coefficient matrix, which is a positive definite symmetric matrix. In the geocentric inertial frame, the control effect generated from

ϕ_{i}

is given as

Δ V_{i} = - \frac{\partial X_{i}^{T}}{\partial V_{i}} \frac{\partial ϕ_{i}}{\partial X_{i}} = - [\begin{matrix} \frac{\partial X_{i}^{T}}{\partial \dot{X}} \frac{\partial ϕ_{i}}{\partial X_{i}} \\ \frac{\partial X_{i}^{T}}{\partial \dot{Y}} \frac{\partial ϕ_{i}}{\partial X_{i}} \\ \frac{\partial X_{i}^{T}}{\partial \dot{Z}} \frac{\partial ϕ_{i}}{\partial X_{i}} \end{matrix}]

(32)

where

V_{i} = {[\dot{X}, \dot{Y}, \dot{Z}]}^{T}

represents the geocentric velocity vector of satellite i. Hence, the velocity increments of

ϕ_{i}

are the negative gradient of

ϕ_{i}

with respect to satellite i’s geocentric velocity vector

V_{i}

. The stability proof will be given in Section 4.

(1): RAAN Control

To limit the RAAN difference, the artificial potential function is designed in Equation (33). The control variable takes

X_{i} = {[sin i_{i} sin Ω_{i}, sin i_{i} cos Ω_{i}]}^{T}

. The control coefficient matrix takes

K = diag (k_{i Ω}, k_{i Ω})

δ Ω_{i} = Ω_{i} - Ω_{t}

i_{t}

and

Ω_{t}

represent the target inclination and RAAN values under

J_{2}

perturbation of the designed orbital plane, respectively. Therefore,

ϕ_{i Ω}

limits the RAAN difference within the expected range and corrects the inclination meanwhile.

ϕ_{i, i Ω} = \{\begin{matrix} k_{i Ω} {(sin i_{i} sin Ω_{i} - sin i_{t} sin Ω_{t})}^{2} + k_{i Ω} {(sin i_{i} cos Ω_{i} - sin i_{t} cos Ω_{t})}^{2}, \\ if δ Ω_{i} > δ Ω_{max} \cup δ Ω_{i} < δ Ω_{min} \\ 0, if δ Ω_{min} \leq δ Ω_{i} \leq δ Ω_{max} \end{matrix}

(33)

According to Equation (32), the velocity increments generated by

ϕ_{i, i Ω}

in the geocentric inertial frame are

\begin{matrix} Δ V_{i, i Ω} = - \frac{\partial X_{i}^{T}}{\partial V_{i}} \frac{\partial ϕ_{i, i Ω}}{\partial X_{i}} \\ = - [\begin{matrix} \frac{\partial ϕ_{i, i Ω}}{\partial sin i_{i} sin Ω_{i}} \frac{\partial sin i_{i} sin Ω_{i}}{\partial \dot{X}} + \frac{\partial ϕ_{i, i Ω}}{\partial sin i_{i} cos Ω_{i}} \frac{\partial sin i_{i} cos Ω_{i}}{\partial \dot{X}} \\ \frac{\partial ϕ_{i, i Ω}}{\partial sin i_{i} sin Ω_{i}} \frac{\partial sin i_{i} sin Ω_{i}}{\partial \dot{Y}} + \frac{\partial ϕ_{i, i Ω}}{\partial sin i_{i} cos Ω_{i}} \frac{\partial sin i_{i} cos Ω_{i}}{\partial \dot{Y}} \\ \frac{\partial ϕ_{i, i Ω}}{\partial sin i_{i} sin Ω_{i}} \frac{\partial sin i_{i} sin Ω_{i}}{\partial \dot{Z}} + \frac{\partial ϕ_{i, i Ω}}{\partial sin i_{i} cos Ω_{i}} \frac{\partial sin i_{i} cos Ω_{i}}{\partial \dot{Z}} \end{matrix}] \end{matrix}

(34)

The relationship between

X_{i}

and

V_{i}

is given as

\begin{matrix} sin i_{i} sin Ω_{i} = \frac{Y \dot{Z} - Z \dot{Y}}{\sqrt{{(Y \dot{Z} - Z \dot{Y})}^{2} + {(Z \dot{X} - X \dot{Z})}^{2} + {(X \dot{Y} - Y \dot{X})}^{2}}} \\ sin i_{i} cos Ω_{i} = \frac{- (Z \dot{X} - X \dot{Z})}{\sqrt{{(Y \dot{Z} - Z \dot{Y})}^{2} + {(Z \dot{X} - X \dot{Z})}^{2} + {(X \dot{Y} - Y \dot{X})}^{2}}} \end{matrix}

(35)

In fact, the out-of-plane control can be avoided by setting an appropriate offset of the semi-major axis to ensure the time derivative of

Ω

is approximately zero under perturbations. Take the GW-2 constellation, for instance. The drift range of

δ Ω

\pm {12.14}^{\circ}

when

E_{min} = 20^{\circ}

according to Equation (18). The drift range of

δ Ω

is relatively large, where the RAAN control will basically not be implemented in two years.

(2): Relative Motion-bound Control

To realize the restriction of relative motion bound, the artificial potential function is designed in Equation (36). The control variable takes

X_{i} = {\overset{⌣}{d}}_{i j}

X_{i} = {\overset{⌢}{d}}_{i j}

, and the control coefficient matrix takes the

K = k_{b o u n d}

ϕ_{i, b o u n d} = \{\begin{matrix} k_{b o u n d} {({\overset{⌣}{d}}_{i j} - d_{mint})}^{2}, if {\overset{⌣}{d}}_{i j} < d_{mint} \\ k_{b o u n d} {({\overset{⌢}{d}}_{i j} - d_{maxt})}^{2}, if {\overset{⌢}{d}}_{i j} > d_{maxt} \\ 0, if {\overset{⌣}{d}}_{i j} \geq d_{mint} \cap {\overset{⌢}{d}}_{i j} \leq d_{maxt} \end{matrix}

(36)

Satellite j is designated as the reference satellite.

d_{mint}

and

d_{maxt}

are determined by Equation (20).

In the geocentric inertial frame, the velocity increments generated by

ϕ_{i, b o u n d}

can be obtained by calculating the negative gradient of

ϕ_{i, b o u n d}

with respect to the relative velocity vector

v_{i} = {[\dot{x}, \dot{y}, \dot{z}]}^{T}

, and then transforming them from satellite j’s orbital frame to the geocentric inertial frame.

Δ V_{i, b o u n d} = - Q_{j}^{- 1} \frac{\partial X_{i}^{T}}{\partial v_{i}} \frac{\partial ϕ_{i, b o u n d}}{\partial X_{i}} = - Q_{j}^{- 1} [\begin{matrix} \frac{\partial ϕ_{i, b o u n d}}{\partial {\overset{⌢}{d}}_{i j}} \frac{\partial {\overset{⌢}{d}}_{i j}}{\partial \dot{x}} + \frac{\partial ϕ_{i, b o u n d}}{\partial {\overset{⌣}{d}}_{i j}} \frac{\partial {\overset{⌣}{d}}_{i j}}{\partial \dot{x}} \\ \frac{\partial ϕ_{i, b o u n d}}{\partial {\overset{⌢}{d}}_{i j}} \frac{\partial {\overset{⌢}{d}}_{i j}}{\partial \dot{y}} + \frac{\partial ϕ_{i, b o u n d}}{\partial {\overset{⌣}{d}}_{i j}} \frac{\partial {\overset{⌣}{d}}_{i j}}{\partial \dot{y}} \\ \frac{\partial ϕ_{i, b o u n d}}{\partial {\overset{⌢}{d}}_{i j}} \frac{\partial {\overset{⌢}{d}}_{i j}}{\partial \dot{z}} + \frac{\partial ϕ_{i, b o u n d}}{\partial {\overset{⌣}{d}}_{i j}} \frac{\partial {\overset{⌣}{d}}_{i j}}{\partial \dot{z}} \end{matrix}]

(37)

\partial ({\overset{⌢}{d}}_{i j}, {\overset{⌣}{d}}_{i j}) / \partial \dot{x}

\partial ({\overset{⌢}{d}}_{i j}, {\overset{⌣}{d}}_{i j}) / \partial \dot{y}

\partial ({\overset{⌢}{d}}_{i j}, {\overset{⌣}{d}}_{i j}) / \partial \dot{z}

can be formulated based on Equations (23), (24), (29) and (30).

Q_{j}

represents the transformation matrix from the geocentric inertial frame to satellite j’s orbital frame.

Q_{j} = [\begin{matrix} - sin Ω_{j} cos i_{j} {sinu}_{j} + cos Ω_{j} {cosu}_{j} & cos Ω_{j} cos i_{j} {sinu}_{j} + sin Ω_{j} {cosu}_{j} & sin i_{j} {sinu}_{j} \\ - sin Ω_{j} cos i_{j} {cosu}_{j} - cos Ω_{j} {sinu}_{j} & cos Ω_{j} cos i_{j} {cosu}_{j} - sin Ω_{j} {sinu}_{j} & sin i_{j} {cosu}_{j} \\ sin Ω_{j} sin i_{j} & - cos Ω_{j} sin i_{j} & cos i_{j} \end{matrix}]

(38)

(3): Semi-Major Axis Control

The differences in semi-major axes between constellation satellites will lead to secular drift, and hence the relative motion-bound constraint will never be satisfied. Moreover, the semi-major axis will keep declining because of the influence of atmospheric drag. From this perspective, semi-major axis control is essential to mega constellation maintenance.

The artificial potential function is designed to be a quadratic function about

a_{i}

, given in Equation (39).

X_{i} = a_{i}

K = k_{a}

a_{t}

is the designed value of the orbital semi-major axis. Semi-major axis control will be implemented when the continuous coverage constraint has been satisfied.

ϕ_{i, a} = \{\begin{matrix} k_{a} {(a_{i} - a_{t})}^{2}, if {\overset{⌣}{d}}_{i j} \geq d_{mint} \cap {\overset{⌢}{d}}_{i j} \leq d_{maxt} \cap δ Ω \geq δ Ω_{min} \cap δ Ω \leq δ Ω_{max} \\ 0, if {\overset{⌣}{d}}_{i j} < d_{mint} \cup {\overset{⌢}{d}}_{i j} > d_{maxt} \cup δ Ω < δ Ω_{min} \cup δ Ω > δ Ω_{max} \end{matrix}

(39)

The velocity increment of

ϕ_{i, a}

in the geocentric inertial frame is

Δ V_{i, a} = - \frac{\partial X_{i}^{T}}{\partial V_{i}} \frac{\partial ϕ_{i, a}}{\partial X_{i}} = - [\begin{matrix} \frac{\partial ϕ_{i, a}}{\partial a_{i}} \frac{\partial a_{i}}{\partial \dot{X}} \\ \frac{\partial ϕ_{i, a}}{\partial a_{i}} \frac{\partial a_{i}}{\partial \dot{Y}} \\ \frac{\partial ϕ_{i, a}}{\partial a_{i}} \frac{\partial a_{i}}{\partial \dot{Z}} \end{matrix}]

(40)

The relationship between

a_{i}

and

V_{i}

is given as

a_{i} = - μ {[({\dot{X}}^{2} + {\dot{Y}}^{2} + {\dot{Z}}^{2}) - \frac{2 μ}{\sqrt{X^{2} + Y^{2} + Z^{2}}}]}^{- 1}

(41)

To sum up, the control potential functions for satellite i are given as follows.

ϕ_{i} = \{\begin{matrix} k_{i Ω} {(sin i_{i} sin Ω_{i} - sin i_{t} sin Ω_{t})}^{2} + k_{i Ω} {(sin i_{i} cos Ω_{i} - sin i_{t} cos Ω_{t})}^{2}, \\ if δ Ω > δ Ω_{max} \cup δ Ω < δ Ω_{min} \\ k_{b o u n d} {({\overset{⌣}{d}}_{i j} - d_{mint})}^{2}, if {\overset{⌣}{d}}_{i j} < d_{mint} \\ k_{b o u n d} {({\overset{⌢}{d}}_{i j} - d_{maxt})}^{2}, if {\overset{⌢}{d}}_{i j} > d_{maxt} \\ k_{a} {(a_{i} - a_{t})}^{2}, if δ Ω \geq δ Ω_{min} \cap δ Ω \leq δ Ω_{max} \cap {\overset{⌣}{d}}_{i j} \geq d_{mint} \cap {\overset{⌢}{d}}_{i j} \leq d_{maxt} \end{matrix}

(42)

2.2.3. Self-Organizing Control Rules for Coplanar Satellites

For three adjacent coplanar satellites,

i - 1

, i and

i + 1

, the self-organizing control rules are defined as follows.

ϕ_{i} = \{\begin{matrix} k_{i Ω} {(sin i_{i} sin Ω_{i} - sin i_{t} sin Ω_{t})}^{2} + k_{i Ω} {(sin i_{i} cos Ω_{i} - sin i_{t} cos Ω_{t})}^{2}, \\ if δ Ω > δ Ω_{max} \cup δ Ω < δ Ω_{min} \\ k_{b o u n d} {({\overset{⌣}{d}}_{i, i - 1} - d_{mint})}^{2}, if {\overset{⌣}{d}}_{i, i - 1} < d_{mint} \cap {\overset{⌣}{d}}_{i, i + 1} \geq d_{mint} \\ k_{b o u n d} {({\overset{⌣}{d}}_{i, i + 1} - d_{mint})}^{2}, if {\overset{⌣}{d}}_{i, i - 1} \geq d_{mint} \cap {\overset{⌣}{d}}_{i, i + 1} < d_{mint} \\ k_{b o u n d} {({\overset{⌢}{d}}_{i, i - 1} - d_{maxt})}^{2}, if {\overset{⌢}{d}}_{i, i - 1} > d_{maxt} \cap {\overset{⌢}{d}}_{i, i + 1} \leq d_{maxt} \\ k_{b o u n d} {({\overset{⌢}{d}}_{i, i + 1} - d_{maxt})}^{2}, if {\overset{⌢}{d}}_{i, i - 1} \leq d_{maxt} \cap {\overset{⌢}{d}}_{i, i + 1} > d_{maxt} \\ k_{a} {(a_{i} - a_{t})}^{2}, else \end{matrix}

(43)

To put it simply, the RAAN control will be performed once the RAAN difference between satellite i and the designed value does not fit the

δ Ω

range. The relative motion-bound control will work only if one of two adjacent coplanar satellites of satellite i (denoted as

i - 1

and

i + 1

) does not meet the boundary constraint while the other does. For other cases, the semi-major axis control will be implemented, even if

{\overset{⌢}{d}}_{i, i - 1}

{\overset{⌣}{d}}_{i, i - 1}

{\overset{⌢}{d}}_{i, i + 1}

, and

{\overset{⌣}{d}}_{i, i + 1}

do not all meet the boundary constraint.

2.3. Stability Proof of Satellite Control

In this section, the stability of the proposed self-organizing control is proved based on the Gaussian equations and Lyapunov’s theory. The selection criteria of potential function variables for satellite control are discussed.

2.3.1. Stability Proof of RAAN Control and Semi-Major Axis Control

If there exists one Lyapunov function

ϕ

that satisfies Equation (44), the control system can be proven to be stable.

\{\begin{matrix} ϕ (0) = 0, \\ ϕ (X) > 0, \forall X \neq 0 \\ \dot{ϕ} (0) = 0, \\ \dot{ϕ} (X) \leq 0, \dot{ϕ} (X) ≢ 0, \forall X \neq 0 \end{matrix}

(44)

In this paper,

ϕ_{i}

is designated as the Lyapunov function. Referring to Equation (31), it is obvious that

ϕ_{i}

is positive definite. If the time derivative of

ϕ_{i}

is negative definite, one can prove that the proposed APF control is stable.

For RAAN control and semi-major axis control,

X = X_{i} - X_{t} = [sin i_{i} sin Ω_{i} - sin i_{t} sin Ω_{t},

sin i_{i} cos Ω_{i} - sin i_{t} cos Ω_{t}, a_{i} - a_{t}]^{T}

. Hence,

X

relates to controlled satellite i only.

According to Equation (31), the time derivative of

ϕ_{i}

can be written as

\frac{d ϕ_{i}}{d t} = \frac{\partial ϕ_{i}}{\partial X_{i}^{T}} \frac{d X_{i}^{T}}{d t} + \frac{\partial ϕ_{i}}{\partial X_{t}^{T}} \frac{d X_{t}^{T}}{d t}

(45)

The time derivative of

X_{i}

\begin{matrix} \frac{d X_{i}}{d t} = \frac{\partial X_{i}}{\partial R_{i}^{T}} \frac{d R_{i}}{d t} + \frac{\partial X_{i}}{\partial V_{i}^{T}} \frac{d V_{i}}{d t} \\ = \frac{\partial X_{i}}{\partial R_{i}^{T}} \frac{d R_{i}}{d t} + \frac{\partial X_{i}}{\partial V_{i}^{T}} (- μ \frac{R_{i}}{R_{i}^{3}} + u) \\ = {(\frac{d X_{i}}{d t})}_{t w o - b o d y} + \frac{\partial X_{i}}{\partial V_{i}^{T}} u \end{matrix}

(46)

where

R_{i}

and

V_{i}

are the position vector and velocity vector in the geocentric inertial frame, respectively.

u

represents the control acceleration vector in the geocentric inertial frame.

{(d X_{i} / d t)}_{t w o - b o d y}

represents the time derivative of

X_{i}

under the two-body Keplerian motion condition.

Since

X_{t}

does not change under control, the time derivative of

X_{t}

\frac{d X_{t}}{d t} = {(\frac{d X_{t}}{d t})}_{t w o - b o d y}

(47)

According to Equation (32), the control effect generated by

ϕ_{i}

in geocentric inertial frame is the negative gradient of

ϕ_{i}

with respect to

V_{i} = {[\dot{X}, \dot{Y}, \dot{Z}]}^{T}

u = - \frac{\partial X_{i}^{T}}{\partial V_{i}} \frac{\partial ϕ_{i}}{\partial X_{i}} = - [\begin{matrix} \frac{\partial X_{1}}{\partial \dot{X}} & \frac{\partial X_{2}}{\partial \dot{X}} & \dots & \frac{\partial X_{n}}{\partial \dot{X}} \\ \frac{\partial X_{1}}{\partial \dot{Y}} & \frac{\partial X_{2}}{\partial \dot{Y}} & \dots & \frac{\partial X_{n}}{\partial \dot{Y}} \\ \frac{\partial X_{1}}{\partial \dot{Z}} & \frac{\partial X_{2}}{\partial \dot{Z}} & \dots & \frac{\partial X_{n}}{\partial \dot{Z}} \end{matrix}] [\begin{matrix} \frac{\partial ϕ_{i}}{\partial X_{1}} \\ \frac{\partial ϕ_{i}}{\partial X_{2}} \\ ⋮ \\ \frac{\partial ϕ_{i}}{\partial X_{n}} \end{matrix}]

(48)

Substituting Equations (46)–(48) into Equation (45), the time derivative of

ϕ_{i}

can be rewritten as

\begin{matrix} \frac{d ϕ_{i}}{d t} = \frac{\partial ϕ_{i}}{\partial X_{i}^{T}} ({(\frac{d X_{i}}{d t})}_{t w o - b o d y} + \frac{\partial X_{i}}{\partial V_{i}^{T}} u) + \frac{\partial ϕ_{i}}{\partial X_{t}^{T}} {(\frac{d X_{t}}{d t})}_{t w o - b o d y} \\ = \frac{\partial ϕ_{i}}{\partial X_{i}^{T}} ({(\frac{d X_{i}}{d t})}_{t w o - b o d y} - \frac{\partial X_{i}}{\partial V_{i}^{T}} \frac{\partial X_{i}^{T}}{\partial V_{i}} \frac{\partial ϕ_{i}}{\partial X_{i}}) + \frac{\partial ϕ_{i}}{\partial X_{t}^{T}} {(\frac{d X_{t}}{d t})}_{t w o - b o d y} \\ = - {(\frac{\partial X_{i}^{T}}{\partial V_{i}} \frac{\partial ϕ_{i}}{\partial X_{i}})}^{T} (\frac{\partial X_{i}^{T}}{\partial V_{i}} \frac{\partial ϕ_{i}}{\partial X_{i}}) + \frac{\partial ϕ_{i}}{\partial X_{i}^{T}} {(\frac{d X_{i}}{d t})}_{t w o - b o d y} + \frac{\partial ϕ_{i}}{\partial X_{t}^{T}} {(\frac{d X_{t}}{d t})}_{t w o - b o d y} \\ = - 4 {[\frac{\partial X_{i}^{T}}{\partial V_{i}} K (X_{i} - X_{t})]}^{T} [\frac{\partial X_{i}^{T}}{\partial V_{i}} K (X_{i} - X_{t})] + 2 {(X_{i} - X_{t})}^{T} K {(\frac{d X_{i}}{d t} - \frac{d X_{t}}{d t})}_{t w o - b o d y} \end{matrix}

(49)

For RAAN control and semi-major axis control,

X_{i} = {[sin i_{i} sin Ω_{i}, sin i_{i} cos Ω_{i}, a_{i}]}^{T}

stays constant under the two-body Keplerian motion condition without consideration of perturbations. Hence,

{(d X_{i} / d t)}_{t w o - b o d y}

is zero. Similarly,

{(d X_{t} / d t)}_{t w o - b o d y}

is also zero. The time derivative of

ϕ_{i}

\frac{d ϕ_{i}}{d t} = - 4 {[\frac{\partial X_{i}^{T}}{\partial V_{i}} K (X_{i} - X_{t})]}^{T} [\frac{\partial X_{i}^{T}}{\partial V_{i}} K (X_{i} - X_{t})]

(50)

The control coefficient matrix

K

is positive definite. Therefore, Equation (50) is a negative definite quadratic function when

X_{i} - X_{t} \neq 0

ϕ_{i}

is positive definite, and the time derivative of

ϕ_{i}

is negative definite. According to Lyapunov’s second theory, the RAAN and semi-major axis APFs are stable.

2.3.2. Stability Proof of Relative Motion-Bound Control

For relative motion-bound control,

X = X_{i} - X_{t} = {[{\overset{⌢}{d}}_{i j} - d_{maxt}, {\overset{⌣}{d}}_{i j} - d_{mint}]}^{T}

. Different from RAAN and semi-major axis control, the relative motion-bound control relates to the controlled satellite i and its adjacent coplanar satellite j. Therefore, the time derivative of

X_{i}

will be different. Here we give the stability proof of relative motion-bound control under self-organizing control rules.

Assuming that

ϕ_{i} = k_{bound} {({\overset{⌢}{d}}_{i, i + 1} - d_{maxt})}^{2}

. It can be inferred from Equation (43) that

{\overset{⌢}{d}}_{i, i - 1} \leq d_{maxt}

, and

{\overset{⌢}{d}}_{i, i + 1} > d_{maxt}

. In this case,

X_{i} = {\overset{⌢}{d}}_{i, i + 1}

ϕ_{i}

is positive definite. The time derivative of

X_{i}

can be written as

\begin{matrix} \frac{d X_{i}}{d t} = \frac{\partial X_{i}}{\partial R_{i}^{T}} \frac{d R_{i}}{d t} + \frac{\partial X_{i}}{\partial R_{i + 1}^{T}} \frac{d R_{i + 1}}{d t} + \frac{\partial X_{i}}{\partial V_{i}^{T}} \frac{d V_{i}}{d t} + \frac{\partial X_{i}}{\partial V_{i + 1}^{T}} \frac{d V_{i + 1}}{d t} \\ = {(\frac{d X_{i}}{d t})}_{t w o - b o d y} + \frac{\partial X_{i}}{\partial V_{i}^{T}} u_{i} + \frac{\partial X_{i}}{\partial V_{i + 1}^{T}} u_{i + 1} \\ = \frac{\partial X_{i}}{\partial V_{i}^{T}} u_{i} + \frac{\partial X_{i}}{\partial V_{i + 1}^{T}} u_{i + 1} \end{matrix}

(51)

The time derivative of

X_{t}

\frac{d X_{t}}{d t} = {(\frac{d X_{t}}{d t})}_{t w o - b o d y} = 0

(52)

The control effect working on satellite i is

u_{i} = - \frac{\partial X_{i}^{T}}{\partial V_{i}} \frac{\partial ϕ_{i}}{\partial X_{i}}

(53)

Since

{\overset{⌢}{d}}_{i, i - 1} \leq d_{maxt}

and

{\overset{⌢}{d}}_{i, i + 1} > d_{maxt}

, the control APF working on satellite

i + 1

ϕ_{i + 1} = \{\begin{matrix} k_{bound} {({\overset{⌢}{d}}_{i, i + 1} - d_{maxt})}^{2} = ϕ_{i}, i f {\overset{⌢}{d}}_{i, i + 1} > d_{maxt} a n d {\overset{⌢}{d}}_{i + 1, i + 2} \leq d_{maxt} \\ 0, i f {\overset{⌢}{d}}_{i, i + 1} > d_{maxt} a n d {\overset{⌢}{d}}_{i + 1, i + 2} > d_{maxt} \end{matrix}

(54)

The control effect working on satellite

i + 1

u_{i + 1} = \{\begin{matrix} - \frac{\partial X_{i}^{T}}{\partial V_{i + 1}} \frac{\partial ϕ_{i}}{\partial X_{i}}, i f {\overset{⌢}{d}}_{i, i + 1} > d_{maxt} a n d {\overset{⌢}{d}}_{i + 1, i + 2} \leq d_{maxt} \\ 0, i f {\overset{⌢}{d}}_{i, i + 1} > d_{maxt} a n d {\overset{⌢}{d}}_{i + 1, i + 2} > d_{maxt} \end{matrix}

(55)

Hence, the time derivative of

ϕ_{i}

\begin{matrix} \frac{d ϕ_{i}}{d t} = \frac{\partial ϕ_{i}}{\partial X_{i}^{T}} \frac{d X_{i}}{d t} + \frac{\partial ϕ_{i}}{\partial X_{t}^{T}} \frac{d X_{t}}{d t} = \frac{\partial ϕ_{i}}{\partial X_{i}^{T}} (\frac{\partial X_{i}}{\partial V_{i}^{T}} u_{i} + \frac{\partial X_{i}}{\partial V_{i + 1}^{T}} u_{i + 1}) \\ = \{\begin{matrix} - \frac{\partial ϕ_{i}}{\partial X_{i}^{T}} \frac{\partial X_{i}}{\partial V_{i}^{T}} \frac{\partial X_{i}^{T}}{\partial V_{i}} \frac{\partial ϕ_{i}}{\partial X_{i}} - \frac{\partial ϕ_{i}}{\partial X_{i}^{T}} \frac{\partial X_{i}}{\partial V_{i + 1}^{T}} \frac{\partial X_{i}^{T}}{\partial V_{i + 1}} \frac{\partial ϕ_{i}}{\partial X_{i}}, i f {\overset{⌢}{d}}_{i, i + 1} > d_{maxt} a n d {\overset{⌢}{d}}_{i + 1, i + 2} \leq d_{maxt} \\ - \frac{\partial ϕ_{i}}{\partial X_{i}^{T}} \frac{\partial X_{i}}{\partial V_{i}^{T}} \frac{\partial X_{i}^{T}}{\partial V_{i}} \frac{\partial ϕ_{i}}{\partial X_{i}}, i f {\overset{⌢}{d}}_{i, i + 1} > d_{maxt} a n d {\overset{⌢}{d}}_{i + 1, i + 2} > d_{maxt} \end{matrix} \\ = \{\begin{matrix} - 4 k_{bound}^{2} {({\overset{⌢}{d}}_{i, i + 1} - d_{maxt})}^{2} \frac{\partial {\overset{⌢}{d}}_{i, i + 1}}{\partial V_{i}^{T}} \frac{\partial {\overset{⌢}{d}}_{i, i + 1}}{\partial V_{i}} - 4 k_{bound}^{2} {({\overset{⌢}{d}}_{i, i + 1} - d_{maxt})}^{2} \frac{\partial {\overset{⌢}{d}}_{i, i + 1}}{\partial V_{i + 1}^{T}} \frac{\partial {\overset{⌢}{d}}_{i, i + 1}}{\partial V_{i + 1}}, \\ i f {\overset{⌢}{d}}_{i, i + 1} > d_{maxt} a n d {\overset{⌢}{d}}_{i + 1, i + 2} \leq d_{maxt} \\ - 4 k_{bound}^{2} {({\overset{⌢}{d}}_{i, i + 1} - d_{maxt})}^{2} \frac{\partial {\overset{⌢}{d}}_{i, i + 1}}{\partial V_{i}^{T}} \frac{\partial {\overset{⌢}{d}}_{i, i + 1}}{\partial V_{i}}, i f {\overset{⌢}{d}}_{i, i + 1} > d_{maxt} a n d {\overset{⌢}{d}}_{i + 1, i + 2} > d_{maxt} \end{matrix} \end{matrix}

(56)

Since

k_{bound} > 0

, the time derivative of

ϕ_{i}

is negative definite. To that end, the stability of the relative motion-bound control under self-organizing rules is proven.

2.3.3. Selection Criteria of Quadratic Potential Function Variables

In past studies, the control variables of quadratic APFs can be position states

X, Y, Z

[16,18,25], velocity states

\dot{X}, \dot{Y}, \dot{Z}

[19], relative orbital elements [20,21,22], classical orbital elements [26], etc. The stabilities of quadratic APFs in terms of different control variables are discussed as follows.

(1): Using $R$ and $V$ as Control Variables

In this case,

X_{i} = [R_{i}; V_{i}]

X_{t} = [R_{t}; V_{t}]

. The control effects aim to make the geocentric position and velocity of the controlled satellite converge to that of the target one. The difference between the time derivative of

X_{i}

and that of

X_{t}

under the two-body condition is

{(\frac{d X_{i}}{d t} - \frac{d X_{t}}{d t})}_{t w o - b o d y} = [\begin{matrix} \frac{d R_{i}}{d t} - \frac{d R_{t}}{d t} \\ \frac{d V_{i}}{d t} - \frac{d V_{t}}{d t} \end{matrix}] = [\begin{matrix} V_{i} - V_{t} \\ - \frac{μ R_{i}}{R_{i}^{3}} + \frac{μ R_{t}}{R_{t}^{3}} \end{matrix}]

(57)

When the relative distance between the controlled satellite and the target one is much closer than the orbital altitude, the difference in velocity and two-body gravity in Equation (57) is minor and is ignored. Referring to Equation (49),

{(d X_{i} / d t - d X_{t} / d t)}_{t w o - b o d y} = 0

, and the time derivative of

ϕ_{i}

is negative definite. The control effect is stable. However, the orbital dynamics differ a lot with the increase in the intersatellite distance.

{(d X_{i} / d t - d X_{t} / d t)}_{t w o - b o d y}

can be positive or negative, and the stability cannot be guaranteed anymore.

Therefore, quadratic APF control in terms of

R

and

V

is stable only if the relative distance is small, or if the gravity is weak, for example, on the ground or in deep space.

(2): Using Orbital Elements as Control Variables

In this case,

X_{i} = {[a_{i}, e_{i}, i_{i}, Ω_{i}, ω_{i}, M_{i}]}^{T}

X_{t} = {[a_{t}, e_{t}, i_{t}, Ω_{t}, ω_{t}, M_{t}]}^{T}

. Under two-body Keplerian motion condition, it can be inferred that

{(\frac{d X_{i}}{d t} - \frac{d X_{t}}{d t})}_{t w o - b o d y} = [\begin{matrix} d a_{i} / d t - d a_{t} / d t \\ d e_{i} / d t - d e_{t} / d t \\ d i_{i} / d t - d i_{t} / d t \\ d Ω_{i} / d t - d Ω_{t} / d t \\ d ω_{i} / d t - d ω_{t} / d t \\ d M_{i} / d t - d M_{t} / d t \end{matrix}] = [\begin{matrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ n_{i} - n_{t} \end{matrix}]

(58)

According to Equation (49), if the slow variables of orbital elements, i.e.,

a, e, i, Ω, ω

, are chosen as APF control variables, the control effect is stable.

(3): Using Relative Orbital Elements as Control Variables

In this case,

X_{i} = {[b_{1}, b_{2}, b_{3}, c, θ_{1}, β_{1}, ψ]}^{T}

X_{t} = {[b_{1 t}, b_{2 t}, b_{3 t}, c_{t}, θ_{1 t}, β_{1 t}, ψ_{t}]}^{T}

. There exists

{(\frac{d X_{i}}{d t} - \frac{d X_{t}}{d t})}_{t w o - b o d y} = [\begin{matrix} d b_{1} / d t - d b_{1 t} / d t \\ d b_{2} / d t - d b_{2 t} / d t \\ d b_{3} / d t - d b_{3 t} / d t \\ d c / d t - d c_{t} / d t \\ d θ_{1} / d t - d θ_{1 t} / d t \\ d β_{1} / d t - d β_{1 t} / d t \\ d ψ / d t - d ψ_{t} / d t \end{matrix}] = [\begin{matrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ n_{i} - n_{t} \\ n_{i} - n_{t} \end{matrix}]

(59)

Therefore, relative orbital elements

b_{1}, b_{2}, b_{3}, c, θ_{1}

are ideal control variables that ensure control stability. When

a_{i} = a_{t}

β_{1}, ψ

can be chosen as control variables.

To sum up, the selection criteria of quadratic APF control variables are given as follows. If the time derivative of

X_{i}

is the same as that of

X_{t}

under the two-body condition, the corresponding APF control is stable.

3. Results and Discussion

In this section, the self-organizing control under continuous coverage constraints is simulated. The GW-2 constellation is designated as the mega constellation example. The control consumption is assessed, and the coverage performance is analyzed.

3.1. Self-Organizing Control of Two Satellites

The two-satellite case is simulated, where one satellite, named the chief, performs semi-major axis control only, and the other, named the deputy, performs self-organizing control to keep the relative motion bound within the expected range.

The simulation parameters are listed in Table 2.

The relative motion-bound changing with time is shown in Figure 7, where the blue, black, red, and green curves represent the maximum bound under control, the minimum bound under control, the expected maximum bound, and the expected minimum bound, respectively. It can be found that both the maximum bound and the minimum bound are constrained within the expected range.

Due to the existence of atmospheric drag, the mean semi-major axis drops by 10 m every day without semi-major control, shown in Figure 8. It can be seen that the relative motion-bound control makes the semi-major axis of the deputy converge to that of the chief as well. That is because the expected bound values are fixed, and the fixed bounds exist only if the semi-major axis difference takes zero. Figure 9 shows the semi-major axes of the chief satellite and the deputy under semi-major axis control. The results indicate that the semi-major axes of the two satellites have been maintained around the expected value.

The control effect of the deputy is shown in Figure 10. The accumulated velocity increment is 1.98 m/s, and the single-step velocity increment is less than 0.02 m/s, which is minor enough for mega constellation control.

3.2. Self-Organizing Control of Coplanar Satellites

The coplanar satellites are selected as satellites of one orbit in GW-2 sub-constellation with an orbital altitude of 1145 km and inclination of

60^{\circ}

. The designed configuration of the GW-2 constellation has been presented in Table 1. There are 48 orbital planes and 36 coplanar satellites in the GW-2 sub-constellation total. The simulation parameters are listed in Table 3. The initial errors of

δ Ω, δ ω, δ f

are relatively large, which simulate the drift of

Ω, ω, f

under perturbations without control.

The angular radius of the coverage range is calculated as

α_{max} = {17.19}^{\circ}

. According to Equation (18), the RAAN drift should be within

δ Ω_{min} = - {12.14}^{\circ}

and

δ Ω_{max} = {12.14}^{\circ}

. The drift range of the relative motion bound can be calculated using Equations (15) and (20) based on the extreme value of

δ Ω

The initial distribution of the coplanar satellites in the geocentric inertial system is shown in Figure 11, while the final distribution is shown in Figure 12. The red dots represent satellites that do not satisfy the continuous coverage constraint, named unqualified satellites, whereas the green dots represent satellites that satisfy the continuous coverage constraint. The number of unqualified satellites changing with time is shown in Figure 13. The unqualified satellites are reduced from 11 to 0 due to the self-organizing control, which verifies the feasibility of the proposed control method.

The differences in the maximum relative motion bounds away from their corresponding constraints are shown in Figure 14. The differences should be negative theoretically and they actually are. Hence the maximum bound control has not been performed during the simulation time. The differences in the minimum relative motion bounds away from their corresponding constraints are shown in Figure 15, from which it can be found that all minimum bounds are controlled to be larger than the constraints. The occurrence of the unqualified satellite at 40 days in Figure 13 can be understood by referring to Figure 15. The small drift in the minimum bound away from the constraint due to perturbations is identified as unqualified.

The semi-major axis differences away from target value

a_{t}

are shown in Figure 16. The maximum drift of the semi-major axis at the final moment is 54.27 m. The RAAN differences away from the target RAAN under

J_{2}

perturbation of the designed orbital plane are shown in Figure 17. The maximum drift of RAAN at the final moment is

{9.73}^{\circ}

. The control effects of coplanar satellites in the geocentric inertial system are shown in Figure 18. The exact velocity increments have no value of reference because the initial orbit conditions are dramatically bad. However, Figure 16 reveals the convergence of the proposed self-organizing control.

3.3. Self-Organizing Control of Mega Constellation

This paper focuses on self-organizing control to keep continuous coverage within the designed latitude zones, whereas collision avoidance control has not been considered so far. Since the proposed continuous coverage constraint relates to coplanar satellites only, the self-organizing control of different orbital planes is independent of each other.

The initial simulation parameters are listed in Table 4. In order to demonstrate the effect of relative motion-bound control, the initial errors of

δ Ω, δ ω, δ f

are set as

δ Ω \in [- 5, 5] deg

δ ω \in [- 1, 1] deg

δ f \in [- 1, 1] deg

. The initial and final distributions of the mega constellation are shown in Figure 19 and Figure 20, respectively. Satellites that satisfy the continuous coverage constraint are colored green, while others are colored red. It can be observed that unqualified satellites turn into qualified satellites under self-organizing control. The number of unqualified satellites is reduced from 248 to 0, shown in Figure 21. The control effects in the geocentric inertial system changing with time are shown in Figure 22, which reveals the convergence of the proposed self-organizing control.

The time histories of the maximum bound difference, minimum bound difference, semi-major axis difference, and RAAN difference are given in Figure 23, Figure 24, Figure 25 and Figure 26. All the minimum relative motion bounds converge to expected constraints automatically under control. The semi-major axis differences are restricted within 98.02 m, and the maximum drift of RAAN at the final moment is

{5.46}^{\circ}

, which is far away from the

δ Ω

drift constraint.

Figure 27 shows the coverage performance of the GW-2 sub-constellation at the final moment, where yellow dots represent virtual ground stations that can be single-covered, and red dots represent ground stations that cannot be single-covered. The results indicate that the GW-2 sub-constellation can provide continuous coverage within

\pm 78^{\circ}

along latitudes. Figure 28 shows the multiplicity of coverage of the GW-2 sub-constellation, where the blue curve denotes the initial condition and the green curve denotes the final condition. It can be found that the GW-2 constellation always provides continuous coverage within

\pm 78^{\circ}

latitude zones during the simulation time.

4. Conclusions

In view of the numerous satellites in mega constellations, this paper innovatively introduces the self-organizing concept of mega constellation control. The emphasis is on realizing continuous Earth observation within the designed latitude zones so that the combination of multiple sub-constellations can provide global remote sensing services. A novel continuous coverage constraint has been presented that relates to intersatellite motion only, including the RAAN drift range and relative motion-bound drift range between two adjacent coplanar satellites. The configuration requirements for mega constellations have been studied, and most ongoing or planned mega constellations satisfy the proposed coverage constraint. Considering RAAN, the relative motion-bound and semi-major axis constraints, the corresponding artificial potential functions as well as the self-organizing rules have been designed. As a highlight, the scale-independent relative orbital elements are used as control variables in mega constellation control because these elements have high precision for near-circular orbits and arbitrary intersatellite ranges. The stability of satellite control using quadratic artificial potential functions has been proven. One key finding is that if the time derivative of the control variable is the same as that of the expected value under two-body Keplerian motion condition, the corresponding quadratic artificial potential function is stable. The self-organizing control simulations of coplanar satellites and the GW-2 sub-constellation have been carried out. The simulation results indicate the goodness of control effect for the small single-step velocity increment and low accumulated control consumption, and the mega constellation satisfies continuous coverage under control.

Author Contributions

Conceptualization: Y.X., Y.Z. and Z.W.; funding acquisition: Z.W. and L.F.; Methodology: Y.X.; validation verification: Y.X. and Y.H.; writing—original draft: Y.X.; writing—review and editing: Y.X. and Y.H. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Natural Science Foundation of China of funder grant number No. U20B2056 and the Huzhou Distinguished Scholar Program of funder grant number ZJIHI-KY0016. The APC was funded by Zhejiang University.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

According to Equations (11) and (12),

cos ∠ BC \leq cos (\frac{2 π}{P})

(A1)

{cos}^{2} 2 α_{max} {sin}^{2} θ \leq {cos}^{2} (\frac{2 π}{P}) ({sin}^{2} θ {cos}^{2} θ - {sin}^{2} δ Ω {cos}^{2} θ + {cos}^{2} 2 α_{max} {sin}^{2} δ Ω)

(A2)

Let

τ = {sin}^{2} θ

1 - τ = {cos}^{2} θ

\begin{matrix} - cos (\frac{2 π}{P}) τ^{2} + [{cos}^{2} (\frac{2 π}{P}) (1 + {sin}^{2} δ Ω) - {cos}^{2} 2 α_{max}] τ \\ - {cos}^{2} (\frac{2 π}{P}) {sin}^{2} 2 α_{max} {sin}^{2} δ Ω \geq 0 \end{matrix}

(A3)

Let

\begin{matrix} A_{1} = - {cos}^{2} (\frac{2 π}{P}) \\ B_{1} = - {cos}^{2} (\frac{2 π}{P}) (1 + {sin}^{2} δ Ω) + {cos}^{2} 2 α_{max} \\ C_{1} = {cos}^{2} (\frac{2 π}{P}) {sin}^{2} 2 α_{max} {sin}^{2} δ Ω \end{matrix}

(A4)

According to Equation (6), it can be proven that

cos (\frac{2 π}{P}) \geq cos (2 α_{max})

(A5)

Therefore,

\begin{matrix} A_{1} \geq 0 \\ B_{1} \leq - {cos}^{2} 2 α_{max} {sin}^{2} δ Ω \leq 0 \\ C_{1} \geq 0 \end{matrix}

(A6)

The range of

τ

can be obtained by Equation (A3)

\frac{- B_{1} - \sqrt{B_{1}^{2} - 4 A_{1} C_{1}}}{2 A_{1}} \leq τ \leq \frac{- B_{1} + \sqrt{B_{1}^{2} - 4 A_{1} C_{1}}}{2 A_{1}}

(A7)

According to Equation (A6), it is obvious that

\frac{- B_{1} + \sqrt{B_{1}^{2} - 4 A_{1} C_{1}}}{2 A_{1}} \geq \frac{- B_{1} - \sqrt{B_{1}^{2} - 4 A_{1} C_{1}}}{2 A_{1}} \geq 0

(A8)

Therefore,

τ = {sin}^{2} θ \geq 0

can be proven, yet

B_{1}^{2} - 4 A_{1} C_{1} \geq 0

depends on the

δ Ω

value.

The range of

sin θ

\sqrt{\frac{- B_{1} - \sqrt{B_{1}^{2} - 4 A_{1} C_{1}}}{2 A_{1}}} \leq sin θ \leq \sqrt{\frac{- B_{1} + \sqrt{B_{1}^{2} - 4 A_{1} C_{1}}}{2 A_{1}}}

(A9)

The range of

cos θ

\sqrt{\frac{2 A_{1} + B_{1} + \sqrt{B_{1}^{2} - 4 A_{1} C_{1}}}{2 A_{1}}} \leq cos θ \leq \sqrt{\frac{2 A_{1} + B_{1} - \sqrt{B_{1}^{2} - 4 A_{1} C_{1}}}{2 A_{1}}}

(A10)

According to Equations (10) and (13),

{(cos ∠ D F)}_{min} \leq cos (\frac{2 π}{S})

(A11)

{(cos θ)}_{min} \leq cos (\frac{2 π}{S}) cos δ Ω

(A12)

\sqrt{B_{1}^{2} - 4 A_{1} C_{1}} \geq - 2 A_{1} {sin}^{2} (\frac{2 π}{S}) - 2 A_{1} {cos}^{2} (\frac{2 π}{S}) {sin}^{2} δ Ω - B_{1}

(A13)

B_{1}^{2} - 4 A_{1} C_{1} \geq {[2 A_{1} {sin}^{2} (\frac{2 π}{S}) + 2 A_{1} {cos}^{2} (\frac{2 π}{S}) {sin}^{2} δ Ω + B_{1}]}^{2}

(A14)

Equation (A14) ensures that

θ

is a real value, given in Equations (A9) and (A10). Equation (A14) can be simplified as

A_{2} {sin}^{2} δ Ω - B_{2} \geq 0

(A15)

where

\begin{matrix} A_{2} = - 2 {sin}^{2} (\frac{2 π}{S}) {cos}^{2} (\frac{2 π}{S}) {cos}^{2} (\frac{2 π}{P}) - {cos}^{2} (\frac{2 π}{S}) {cos}^{2} 2 α_{max} + {cos}^{2} (\frac{2 π}{P}) {cos}^{2} 2 α_{max} \\ B_{2} = - {sin}^{2} (\frac{2 π}{S}) [{cos}^{2} (\frac{2 π}{S}) {cos}^{2} (\frac{2 π}{P}) - {cos}^{2} 2 α_{max}] \end{matrix}

(A16)

According to Equation (6),

\begin{matrix} A_{2} \leq - 2 {sin}^{2} (\frac{2 π}{S}) {cos}^{2} 2 α_{max} - {cos}^{2} (\frac{2 π}{S}) {cos}^{2} 2 α_{max} + {cos}^{2} (\frac{2 π}{P}) {cos}^{2} 2 α_{max} \\ = [{cos}^{2} (\frac{2 π}{P}) + {cos}^{2} (\frac{2 π}{S}) - 2] {cos}^{2} 2 α_{max} \\ = - [{sin}^{2} (\frac{2 π}{P}) + {sin}^{2} (\frac{2 π}{S})] {cos}^{2} 2 α_{max} \\ \leq 0 \\ B_{2} = - {sin}^{2} (\frac{2 π}{S}) [{cos}^{2} (\frac{2 π}{S}) {cos}^{2} (\frac{2 π}{P}) - {cos}^{2} 2 α_{max}] \leq 0 \end{matrix}

(A17)

δ Ω

is a positive variable. Therefore, the range of

δ Ω

0 \leq δ Ω \leq arcsin (\sqrt{\frac{B_{2}}{A_{2}}})

(A18)

References

Del Portillo, I.; Cameron, B.G.; Crawley, E.F. A technical comparison of three low earth orbit satellite constellation systems to provide global broadband. Acta Astronaut. 2019, 159, 123–135. [Google Scholar] [CrossRef]
Williams, A.; Hainaut, O.; Otarola, A.; Tan, G.H.; Rotola, G. Analysing the impact of satellite constellations and ESO’s role in supporting the astronomy community. arXiv 2021, arXiv:2108.04005. [Google Scholar]
Boggio, M.; Colangelo, L.; Virdis, M.; Pagone, M.; Novara, C. Earth Gravity In-Orbit Sensing: MPC Formation Control Based on a Novel Constellation Model. Remote Sens. 2022, 14, 2815. [Google Scholar] [CrossRef]
Wang, Z.; Zhang, Y.; Wen, G.; Bai, S.; Cai, Y.; Huang, P.; Han, D.; He, Y. Atmospheric Density Model Optimization and Spacecraft Orbit Prediction Improvements Based on Q-Sat Orbit Data. arXiv 2021, arXiv:2112.03113. [Google Scholar]
Mastro, P.; Masiello, G.; Serio, C.; Pepe, A. Change Detection Techniques with Synthetic Aperture Radar Images: Experiments with Random Forests and Sentinel-1 Observations. Remote Sens. 2022, 14, 3323. [Google Scholar] [CrossRef]
Lang, T. Low Earth orbit satellite constellations for continuous coverage of the mid-latitudes. In Proceedings of the Astrodynamics Conference, San Diego, CA, USA, 29 July–31 July 1996; p. 3638. [Google Scholar]
Santos, M.; Shapiro, B. Relating satellite coverage to orbital geometry. In Proceedings of the AIAA/AAS Astrodynamics Specialist Conference and Exhibit, Honolulu, HI, USA, 18–21 August 2008; p. 6609. [Google Scholar]
Hackett, T.M.; Bilén, S.G.; Bell, D.J.; Lo, M.W. Geometric approach for analytical approximations of satellite coverage statistics. J. Spacecr. Rocket. 2019, 56, 1286–1299. [Google Scholar] [CrossRef]
Dai, G.; Chen, X.; Wang, M.; Fernández, E.; Nguyen, T.N.; Reinelt, G. Analysis of satellite constellations for the continuous coverage of ground regions. J. Spacecr. Rocket. 2017, 54, 1294–1303. [Google Scholar] [CrossRef] [Green Version]
Gong, Y.; Zhang, S.; Peng, X. Quick coverage analysis of mega Walker Constellation based on 2D map. Acta Astronaut. 2021, 188, 99–109. [Google Scholar] [CrossRef]
Huang, S.; Colombo, C.; Bernelli-Zazzera, F. Multi-criteria design of continuous global coverage Walker and Street-of-Coverage constellations through property assessment. Acta Astronaut. 2021, 188, 151–170. [Google Scholar] [CrossRef]
Arnas, D.; Casanova, D.; Tresaco, E. Relative and absolute station-keeping for two-dimensional–lattice flower constellations. J. Guid. Control. Dyn. 2016, 39, 2602–2604. [Google Scholar] [CrossRef]
Chen, Y.; Zhao, L.; Liu, H.; Li, L.; Liu, J. Analysis of configuration and maintenance strategy of LEO walker constellation. J. Astronaut. 2019, 40, 1296–1303. [Google Scholar]
Li, J.; Hu, M.; Wang, X.; Li, F.; Xu, J. Analysis of configuration offsetting maintenance method for LEO Walker constellation. Chin. Space Sci. Technol. 2021, 41, 38. [Google Scholar]
Izzo, D.; Pettazzi, L. Autonomous and distributed motion planning for satellite swarm. J. Guid. Control. Dyn.s 2007, 30, 449–459. [Google Scholar] [CrossRef]
Nag, S.; Summerer, L. Behaviour based, autonomous and distributed scatter manoeuvres for satellite swarms. Acta Astronaut. 2013, 82, 95–109. [Google Scholar] [CrossRef]
Hu, Q.; Dong, H.; Zhang, Y.; Ma, G. Tracking control of spacecraft formation flying with collision avoidance. Aerosp. Sci. Technol. 2015, 42, 353–364. [Google Scholar] [CrossRef]
McInnes, C.R. Autonomous proximity manoeuvring using artificial potential functions. ESA J. 1993, 17, 159–169. [Google Scholar]
Cao, L.; Qiao, D.; Xu, J. Suboptimal artificial potential function sliding mode control for spacecraft rendezvous with obstacle avoidance. Acta Astronaut. 2018, 143, 133–146. [Google Scholar] [CrossRef]
Spencer, D.A. Automated trajectory control using artificial potential functions to target relative orbits. J. Guid. Control. Dyn. 2016, 39, 2142–2148. [Google Scholar] [CrossRef]
Wang, Z.; Xu, Y.; Jiang, C.; Zhang, Y. Self-organizing control for satellite clusters using artificial potential function in terms of relative orbital elements. Aerosp. Sci. Technol. 2019, 84, 799–811. [Google Scholar] [CrossRef]
Xu, Y.; Wang, Z.; Zhang, Y. Bounded flight and collision avoidance control for satellite clusters using intersatellite flight bounds. Aerosp. Sci. Technol. 2019, 94, 105425. [Google Scholar] [CrossRef]
Xu, Y.; Wang, Z.; Zhang, Y. Modeling and control of scale-independent relative orbital elements for near-circular orbits. Acta Astronaut. 2022, 198, 642–658. [Google Scholar]
Jiang, C.; Wang, Z.; Zhang, Y. Decomposition analysis of spacecraft relative motion with different inter-satellite ranges. Acta Astronaut. 2019, 163, 56–68. [Google Scholar]
An, M.; Wang, Z.; Zhang, Y. Self-organizing control strategy for asteroid intelligent detection swarm based on attraction and repulsion. Acta Astronaut. 2017, 130, 84–96. [Google Scholar] [CrossRef]
Sun, Y.; Shen, H. The control of mega-constellation at low earth orbit based on TLE. Acta Astronaut. 2020, 42, 156. [Google Scholar]

Figure 1. Satellite coverage geometry.

Figure 2. Street-of-coverage geometry of one orbital plane.

Figure 3. Streets-of-coverage geometries of adjacent orbital planes.

Figure 4. Street-of-coverage geometry of the perturbed coplanar satellites. D and E are two sub-satellite points. A and C are two intersections of two coverage circle. Point A, B, C form a right triangle. Point D, E, F form a right triangle.

Figure 5. Illustration of the scale-independent relative orbital elements (

x - y

plane).

Figure 5. Illustration of the scale-independent relative orbital elements (

x - y

plane).

Figure 6. Illustration of the scale-independent relative orbital elements (

x - z

plane).

Figure 6. Illustration of the scale-independent relative orbital elements (

x - z

plane).

Figure 7. Relative motion-bound changing with time.

Figure 8. Semi-major axis time histories without semi-major axis control.

Figure 9. Semi-major axis time histories under semi-major axis control.

Figure 10. Control effects changing with time.

Figure 11. The initial distribution of the coplanar satellites in the geocentric inertial system. Red dots represent satellites that do not satisfy the continuous coverage constraint while green dots represent satellites that satisfy the continuous coverage constraint.

Figure 12. The final distribution of the coplanar satellites in the geocentric inertial system.

Figure 13. Number of unqualified satellites in the coplanar satellites.

Figure 14. The differences in maximum bounds away from the corresponding constraints. Red line is zero-value standard. Blue lines are maximum bound differences.

Figure 15. The differences in minimum bounds away from the corresponding constraints. Red line is zero-value standard. Blue lines are minimum bound differences.

Figure 16. Semi-major axis differences away from the target value.

Figure 17. RAAN differences away from the target RAAN under

J_{2}

perturbation. Red lines represent

δ Ω_{\min}

and

δ Ω_{\max}

respectively. Blue lines represent RAAN differences.

Figure 17. RAAN differences away from the target RAAN under

J_{2}

perturbation. Red lines represent

δ Ω_{\min}

and

δ Ω_{\max}

respectively. Blue lines represent RAAN differences.

Figure 18. Control effects of coplanar satellites in the geocentric inertial system.

Figure 19. Initial distribution of the mega constellation. Red dots represent satellites that do not satisfy the continuous coverage constraint while green dots represent satellites that satisfy the continuous coverage constraint.

Figure 20. Final distribution of the mega constellation.

Figure 21. Number of unqualified satellites in the mega constellation.

Figure 22. Control effects of constellation satellites.

Figure 23. Maximum bound differences of constellation satellites. Red line is zero-value standard. Blue lines are maximum bound differences.

Figure 24. Minimum bound differences of constellation satellites. Red line is zero-value standard. Blue lines are minimum bound differences.

Figure 25. Semi-major axis differences of constellation satellites.

Figure 26. RAAN differences of constellation satellites. Red lines represent

δ Ω_{\min}

and

δ Ω_{\max}

respectively. Blue lines represent RAAN differences.

Figure 26. RAAN differences of constellation satellites. Red lines represent

δ Ω_{\min}

and

δ Ω_{\max}

respectively. Blue lines represent RAAN differences.

Figure 27. Coverage performance of the GW-2 sub-constellation.

Figure 28. Multiplicity of coverage of the GW-2 sub-constellation at the initial and final moments.

Table 1. Summary of several mega constellation missions.

Constellation	Altitude	Inclination	Planes	Satellites per Plane	Whether Continuous Coverage Constraint Is Satisfied
Starlink Genaration 1	550 km	53 $^{\circ}$	72	22	Yes
	540 km	53.2 $^{\circ}$	72	22	Yes
	570 km	70 $^{\circ}$	36	20	Yes
	560 km	97.6 $^{\circ}$	6	58	No
	560 km	97.6 $^{\circ}$	4	43	No
OneWeb	1120 km	87.9 $^{\circ}$	18	40	Yes
	1150 km	53 $^{\circ}$	32	50	Yes
	1110 km	53.8 $^{\circ}$	32	50	Yes
GW-A59	590 km	85 $^{\circ}$	16	30	No
	600 km	50 $^{\circ}$	40	50	Yes
	508 km	55 $^{\circ}$	60	60	Yes
GW-2	1145 km	30 $^{\circ}$	48	36	Yes
	1145 km	40 $^{\circ}$	48	36	Yes
	1145 km	50 $^{\circ}$	48	36	Yes
	1145 km	60 $^{\circ}$	48	36	Yes

Table 2. Simulation parameters for the two-satellites case.

Parameters	Values
Simulation Time	60 days
Step Length	600 s
Initial Conditions	$[a_{c}, e_{c}, i_{c}, Ω_{c}, ω_{c}, M_{c}] = [R_{e} + 550 km, 0, 53^{\circ}, 30^{\circ}, 40^{\circ}, 50^{\circ}]$
	$[a_{d}, e_{d}, i_{d}, Ω_{d}, ω_{d}, M_{d}] = [a_{c}, e_{c} + 0.01, i_{c}, Ω_{c} + 1^{\circ}, ω_{c} + 1^{\circ}, M_{c} + 10^{\circ}]$
Perturbations	$J_{2}$ , atmospheric drag, lunar gravity, solar gravity
Control Coefficients	$k_{bound} = 2 \times 10^{- 8}, k_{a} = 2 \times 10^{- 8}$
Ballistic Coefficient	$C_{D} A / m = 0.01$
Expected Range	$d_{maxt} = d_{max} - 1 km, d_{mint} = d_{min} + 1 km, a_{t} = R_{e} + 500 km$

Table 3. Simulation parameters for the coplanar satellites case.

Parameters	Values
Simulation Time	60 days
Step Length	600 s
Designed Configuration	$a_{t} = R_{e} + 1145 km$
	$e_{t} = 0$
	$i_{t} = 60^{\circ}$
	$Ω_{t} = 0^{\circ}$
	$ω_{t} = 0^{\circ}$
Initial Orbit Errors	$δ a = [- 100, 100] m, uniform$
	$δ e = [0, 0.1], uniform$
	$δ i = [- 0.1, 0.1] deg, uniform$
	$δ Ω = [- 10, 10] deg, uniform$
	$δ ω = [- 1, 1] deg, uniform$
	$δ f = [- 1, 1] deg, uniform$
Perturbations	$J_{2}$ , atmospheric drag, lunar gravity, solar gravity
Control Coefficients	$k_{bound} = 2 \times 10^{- 8}, k_{a} = 2 \times 10^{- 8}$
Ballistic Coefficient	$C_{D} A / m = 0.01$
Minimum Elevation Angle	$E_{min} = 20^{\circ}$

Table 4. Simulation parameters for the mega-constellation case.

Parameters	Values
Simulation Time	60 days
Step Length	600 s
Designed Configuration	$a_{t} = R_{e} + 1145 km$
	$e_{t} = 0$
	$i_{t} = 60^{\circ}$
	$P = 48$
	$S = 36$
Initial Orbit Errors	$δ a = [- 100, 100] m, uniform$
	$δ e = [0, 0.1], uniform$
	$δ i = [- 0.1, 0.1] deg, uniform$
	$δ Ω = [- 5, 5] deg, uniform$
	$δ ω = [- 1, 1] deg, uniform$
	$δ f = [- 1, 1] deg, uniform$
Perturbations	$J_{2}$ , atmospheric drag, lunar gravity, solar gravity
Control Coefficients	$k_{bound} = 2 \times 10^{- 8}, k_{a} = 2 \times 10^{- 8}$
Ballistic Coefficient	$C_{D} A / m = 0.01$
Minimum Elevation Angle	$E_{min} = 20^{\circ}$

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Xu, Y.; Zhang, Y.; Wang, Z.; He, Y.; Fan, L. Self-Organizing Control of Mega Constellations for Continuous Earth Observation. Remote Sens. 2022, 14, 5896. https://doi.org/10.3390/rs14225896

AMA Style

Xu Y, Zhang Y, Wang Z, He Y, Fan L. Self-Organizing Control of Mega Constellations for Continuous Earth Observation. Remote Sensing. 2022; 14(22):5896. https://doi.org/10.3390/rs14225896

Chicago/Turabian Style

Xu, Yun, Yulin Zhang, Zhaokui Wang, Yunhan He, and Li Fan. 2022. "Self-Organizing Control of Mega Constellations for Continuous Earth Observation" Remote Sensing 14, no. 22: 5896. https://doi.org/10.3390/rs14225896

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu

Self-Organizing Control of Mega Constellations for Continuous Earth Observation

Abstract

1. Introduction

2. Materials and Methods

2.1. Continuous Coverage Constraint of Mega Constellation

2.1.1. Definition of Continuous Coverage

2.1.2. Configuration Requirements for Mega Constellations

2.1.3. Intersatellite Motion Constraints

2.2. Self-Organizing Control Using Artificial Potential Functions

2.2.1. Scale-Independent Relative Orbital Elements and Bound Formulations

2.2.2. Design of Artificial Potential Functions for Single Satellite Control

2.2.3. Self-Organizing Control Rules for Coplanar Satellites

2.3. Stability Proof of Satellite Control

2.3.1. Stability Proof of RAAN Control and Semi-Major Axis Control

2.3.2. Stability Proof of Relative Motion-Bound Control

2.3.3. Selection Criteria of Quadratic Potential Function Variables

3. Results and Discussion

3.1. Self-Organizing Control of Two Satellites

3.2. Self-Organizing Control of Coplanar Satellites

3.3. Self-Organizing Control of Mega Constellation

4. Conclusions

Author Contributions

Funding

Data Availability Statement

Conflicts of Interest

Appendix A

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI