Adaptive Orbital Rendezvous Control of Multiple Satellites Based on Pulsar Positioning Strategy

Abstract

This paper addresses the orbital rendezvous control for multiple uncertain satellites. Against the background of a pulsar-based positioning approach, a geometric trick is applied to determine the position of satellites. A discontinuous estimation algorithm using neighboring communications is proposed to estimate the target’s position and velocity in the Earth’s Centered Inertial Frame for achieving distributed rendezvous control. The variables generated by the dynamic estimation are viewed as virtual reference trajectories for each satellite in the group, followed by a novel saturation-like adaptive control law with the assumption that the masses of satellites are unknown and time-varying. The rendezvous errors are proven to be convergent to zero asymptotically. Numerical simulations considering the measurement fluctuations validate the effectiveness of the proposed control law.

1. Introduction

In recent years, the cooperative rendezvous of orbital satellites has become a popular research topic among the academic community [1]. Advanced positioning and rendezvous control techniques are building a firm base for the potential applications, such as orbital maintenance, orbital refueling, and orbital assembly [2], via steering multiple orbital satellites to achieve rendezvous at a certain target. Particularly, research on the pulsar-based cooperative rendezvous control is of practical significance as the pulsar sources feature high-precision yet stable timing properties for determining the position coordinates of orbital vehicles [3,4].

The pulsar-based positioning technology allows the orbital vehicles to locate themselves by comparing the received signals from pulsar sources with a database of known pulsars and locations [5]. It is the next-generation navigation technology for orbiting or interplanetary spacecrafts [6] and an alternative calibrated source for GNSS (Global Navigation Satellite System) [7]. Compared with the positioning technique, the control technology plays a more important role as the mission success can only be achieved with a robust control design that provides orbital vehicles with robust properties towards the external disturbance, measurement noise, and modeling uncertainty [8,9,10,11,12,13]. In Ref. [8], the extended Kalman filter is applied to denoising the virtual noisy pulsar signal within Poisson distribution and the state observing criteria of a linearized pulsar model via using pulsar data is also proposed. With the help of orbit information and the long-term observation of a single pulsar, navigation algorithms developed based on the adaptive divided different Kalman filter under scenarios of 1–3 orbital satellites are reported in [9]. For improving the reliability, robust control design features the capacity of overcoming the problem of either parameter variation or external disturbances [10,11,14]. Considering the limitation of fuel storage of orbital spacecraft described by Clohessy–Wiltshire equations of motion, the robust $L_1$ control strategy shown in [10] has achieved better fuel efficiency during orbit transferring even with parameter variations. One other fuel optimization strategy based on the saw-like control updating algorithm can be found in [12]. As the command from the ground control station to the orbital satellites features an obvious delay phenomenon, advanced control schemes such as [15] can be applicable to fill this gap. Sometimes the trade-off problem between impulsive efficiency and task complexity should even be predetermined before launching the spacecraft into the space [13]. Though the results reported in [8,9,10,11,12,13,14,15] have solved various control problems, only solutions for a single satellite are developed without considering cooperative control of multiple orbital satellites.

As the complexity increases in space missions [16], the cooperative missions of multiple orbital vehicles have received lots of attention [17,18,19,20,21,22,23]. The formation control of numerous spacecrafts is studied in [17] under the J2 perturbation caused by the oblateness of the Earth. Robotic arms and wheeled mobile robots are generally used for ground algorithm investigations [18,19]. Two ground 6-DOF robotic arms are used for performing lidar-based rendezvous and docking control design for orbital satellites [18], in which the real-time target pose estimation and tracking algorithms are applied. A wheeled satellite simulator testbed with experimental docking discussions is reported in [19], illustrating the basic idea of maneuvering two spacecraft to achieve rendezvous in the final phase. In addition, applying atmosphere lift and drag force for satellites in low Earth orbits is proved to be an alternative to complete rendezvous and docking missions [20]. It is worth noting that the control schemes proposed in [17,18,19,20] convert the orbital formation into a local tracking problem and only achieve the formation control of two satellites. This, however, is not applicable for more general tasks such as orbital assembly that requires multiple orbital vehicles to rendezvous at the target orbit simultaneously [21,22,23]. In Ref. [21], the authors propose an adaptive fuzzy control law that ensures the spacecrafts’ rendezvous errors to be ultimately bounded by small numbers. In Ref. [22], the Lyapunov barrier functions are applied to solve the inner-agent collision avoidance problem of multiple orbital spacecrafts for safe maneuver. Considering the constrained field of view of vision-based devices, the feasible path generating algorithm reported in [23] ensures that the target spacecraft maintains the view of the camera during the rendezvous process. However, the cooperative algorithms developed for multiple satellites in [21,22,23] only achieve small attraction regions of the closed-loop system, i.e., the rendezvous errors must be initialized as tens of meters for maintaining the effectiveness of the control law or the model simplifications. In practice, the orbital vehicles launched for cooperative missions are orbiting in different orbits, which means that the initial distance between them might be hundreds of kilometers, which cannot be dealt with by the control laws reported in [21,22,23]. Additionally, the cooperative control schemes reported in [17,18,19,20,21,22,23] contain the assumption that all orbital vehicles can know the information of the reference target, which is not a necessity for achieving cooperative missions.

Motivated by the discussions above, this paper makes a further endeavor to solve the adaptive rendezvous control problem of multiple uncertain orbital satellites with the background of the pulsar-based positioning method. First, we introduce a direct geometric trick to determine orbital coordinates via assuming simultaneous observations of three pulsars for each satellite. Second, we consider the scenario that only part of the satellites in the group have access to the states of the orbital reference vehicle and propose a discontinuously distributed estimation algorithm to estimate the reference orbital vehicle. Specifically, the estimation errors are globally exponentially convergent to zero. Third, we design a novel saturation-like adaptive control strategy via viewing the estimation algorithm as a virtual reference trajectory for each networked satellite. It is proven in the Lyapunov sense that the rendezvous errors of all networked satellites converge to zero asymptotically.

The main contribution of this paper is presenting the design and analysis of a novel adaptive rendezvous strategy for multiple uncertain orbital satellites with an initial pulsar-based positioning method. In comparison with control schemes developed for the single orbital satellites in [8,9,10,11,12,13] or the rendezvous of two spacecrafts in [17,18,19,20], the cooperative control algorithm in this paper allows numerous satellites to achieve rendezvous at the reference target. The adaptive capacity towards unknown and time-varying masses provided by the presented control scheme also features more practical significance than that only considers constant mass in [23]. Additionally, the proposed rendezvous algorithm still works even if the initial rendezvous errors are initialized in hundreds of kilometers, largely expanding the efficient range of a region of attraction in comparison with the control laws reported in [21,22,23].

The remainder is organized as follows. Section 2 presents some useful preliminaries for control formulation. In particular, the pulsar-based position scheme is introduced. Section 3 elaborates on the designs of the estimation algorithm and adaptive control law. Numerical simulations are included in Section 4. Section 5 concludes the work briefly.

2. Preliminaries

2.1. Notations and Definitions

In this paper, $\mathbb{R}$ denotes the real number set, $\parallel ·\parallel$ is the Euclidean norm, $|·|$ represents the absolute value of a scalar, diag$\{·\}$ forms a diagonal matrix, ${\mathit{I}}_n$ is an n-dimensional identity matrix within a vector, $1_n$ is an n-dimensional identity row vector, $0_n$ denotes an n-dimensional zero row vector. For a square matrix, $\lambda (·)$, ${\lambda }_m(·)$ and ${\lambda }_M(·)$ represent the eigenvalue, the smallest and largest eigenvalue, respectively. The time t is sometimes omitted without making confusion.

2.2. Pulsar-Based Positioning

Position measurement for orbital satellites is essential for the onboard control systems to maintain stability. It highly relies on the Telemetry, Tracking, and Command (TT&C) stations built on the ground or other orbital GNSS satellites, which requires continuous yet enormous investment. In comparison with that, pulsars are natural high-precision timing sources and are suitable for determining the position coordinates of satellites in both orbital and interplanetary missions. Henceforth, using pulsars as positioning sources is of practical significance for rendezvous control.

In what follows, an initial pulsar-based positioning method is presented.

Assumption 1.

The solar system barycenter coincides with the center of mass of the Sun.

Assumption 2.

The axial precession of the Earth is neglected and the Earth orbit is a pure Keplerian orbit without being affected by other celestial bodies except the Sun.

With Assumptions 1 and 2, the geometric relationship among pulsars, Sun, Earth, and satellites can be shown in Figure 1. The Sun-Centered-Inertial (SCI) frame $O_s-X_s{Y}_s{Z}_s$ that Earth is rotating around is viewed as an inertial frame, where $O_s$ is the center of mass of the Sun. The $X_s$-axis coincides with the direction from the vernal equinox to the $O_s$, the $Z_s$-axis is perpendicular to the Earth’s orbit plane, and the $Y_s$-axis completes the right-hand law according to $Z\times X$. For the Earth-Centered-Inertial (ECI) plane $O_E-X_E{Y}_E{Z}_E$ that is centered by the Earth, the $X_E$-axis features the direction to the vernal equinox, the $Z_E$-axis points the rotating axis of the Earth upwards, and the $Y_E$ satisfies the right-hand law of $Z_E\times X_E$. In $O_s-X_s{Y}_s{Z}_s$ frame, the Earth orbit around the Sun satisfies the following Keplerian dynamics [24],

$${\dot{\overrightarrow{r}}}_{s,e}={\overrightarrow{v}}_{s,e},{\dot{\overrightarrow{v}}}_{s,e}=-\frac{{\mu }_s}{r_{s,e}^3}{\overrightarrow{r}}_{s,e},$$

(1)

where ${\overrightarrow{r}}_{s,e}$ and ${\overrightarrow{v}}_{s,e}$ denote the position and velocity of the Earth, respectively.

Let ${\alpha }_j$ and ${\lambda }_j$ denote the right ascension and declination in frame $O_s-X_s{Y}_s{Z}_s$ of the pulsar j, respectively. The radiation signal of the pulsar j has the following direction [9]:

$${\overrightarrow{q}}_j=[cos{\alpha }_{j}cos{\lambda }_j,cos{\alpha }_{j}sin{\lambda }_j,sin{\alpha }_j].$$

(2)

Suppose that there are n satellites and the coordinate of i-th satellite is given by row vector ${\overrightarrow{r}}_{s,i}=[x_{s,i},y_{s,i},z_{s,i}]$. Generally, the projection of the ${\overrightarrow{r}}_{s,i}$ on the vector ${\overrightarrow{q}}_j$ is supposed to be measurable [5,24], i.e., ${\displaystyle {\overrightarrow{s}}_{i,j}=\left({\overrightarrow{r}}_{s,i}·{\overrightarrow{q}}_j\right){\overrightarrow{q}}_j}$ is known, where $r_{s,i}=\parallel {\overrightarrow{r}}_{s,i}\parallel =c(t_{st}-t_{ssb})$ with c being the light speed, $t_{st}$ the time instant that the satellite receives the pulsar signal, and $t_{ssb}$ the time of solar system barycenter. Henceforth, the plane perpendicular to the vector ${\overrightarrow{s}}_{i,j}$ and passing the i-th satellite can be given by:

$${\overrightarrow{s}}_{i,j}·({\overrightarrow{r}}_{s,i}-{\overrightarrow{s}}_{i,j})=0.$$

(3)

Actually, the position vector of the satellite ${\overrightarrow{r}}_{s,i}$ in the $O_s-X_s{Y}_s{Z}_s$ frame can be uniquely determined via solving the linear equation set in the form of Equation (3) if there are more than three pulsars in the detecting area of the satellite. This is because three non-parallel two-dimensional planes would uniquely determine one single point in the three-dimensional surface [25]. For this sake, we make the following assumption:

Assumption 3.

For each satellite, at least three pulsars’ signals can be observed.

Remark 1.

Note that the state-of-the-art instruments designed for observing pulsar sources, such as those launched by NICER and XPNAV missions, can only perform observation on a single pulsar during certain time windows. Assumption 3 is built upon future technology wherein the onboard equipment allows the orbital satellites to observe multiple pulsar signals simultaneously, which would be achieved by improving sensor sensitivity and minimization packaging.

2.3. Relative Satellite Dynamics in ECI Frame

Considering n orbital satellites under ECI coordinates, let row vectors ${\overrightarrow{r}}_i=\left[x_{e,i},y_{e,i},z_{e,i}\right]$ and ${\overrightarrow{v}}_{e,i}={\dot{\overrightarrow{r}}}_{e,i}=\left[{\dot{x}}_{e,i},{\dot{y}}_{e,i},{\dot{z}}_{e,i}\right]$ denote the ECI position and velocity of the i-th satellite, respectively. The satellites can generate arbitrary thrust force to control their positions and the position control is decoupled from attitude loop. Via ignoring the J2 perturbation caused by Earth obslateness, the unperturbed orbital satellites can be described by [24]:

$${\dot{\overrightarrow{r}}}_{e,i}={\overrightarrow{v}}_{e,i},{\dot{\overrightarrow{v}}}_{e,i}=-\frac{{\mu }_e}{r_{e,i}^3}{\overrightarrow{r}}_{e,i}+\frac{{\tau }_i}{m_i\left(t\right)},i\in \mathcal{N}\triangleq \{1,2,\cdots ,n\},$$

(4)

where ${\mu }_e$ is the gravitational constant of the Earth, $m_i\left(t\right)$ denotes the mass of the i-the satellite, and ${\tau }_i\in {\mathbb{R}}^3$ represents the control input. We suppose that ${\dot{m}}_i=\alpha \parallel {\tau }_i\parallel$, where $\alpha <0$ is a negative constant related to specific impulse coefficient. Via observing Equations (1), (2) and (4), the position of the i-th satellite satisfies

$${\overrightarrow{r}}_{e,i}={\overrightarrow{r}}_{s,i}-{\overrightarrow{r}}_{s,e},$$

(5)

which reveals the basic principle of orbital satellite positioning by measuring pulsar signals. We assume that the satellites can obtain Earth’s position vector in the SCI frame. In reality, the satellite’s total mass is unknown and decreasing due to generating certain thrust forces and torques via burning the fuel. For this sake, the mass $m_i\left(t\right)$ and its ratio of change ${\dot{m}}_i\left(t\right)$ are supposed to satisfy:

Assumption 4.

For all $i\in \mathcal{N}$, the mass $m_i\left(t\right)$ is an unknown variable and satisfies

$$m_i\left(t\right)\ge {\underline{m}}_i,$$

(6)

where ${\underline{m}}_i$ is a known positive constant, and the ratio of change of mass ${\dot{m}}_i\left(t\right)$ is a known non-positive bounded variable and ${\ddot{m}}_i\left(t\right)$ is bounded.

Besides n satellites, we suppose that there is one single free-flying reference orbital vehicle indexed by 0 and described by:

$${\dot{\overrightarrow{r}}}_{e,0}={\overrightarrow{v}}_{e,0},{\dot{\overrightarrow{v}}}_{e,0}=-\frac{{\mu }_e}{r_{e,0}^3}{\overrightarrow{r}}_{e,0}.$$

(7)

Actually, the state vector of the satellite described by Equation (7) can be uniquely determined by its orbital elements that include the eccentricity, semimajor axis, inclination, longitude of the ascending node, the argument of perigee, and the true anomaly, as well as other equivalent expressions [24]. For endowing the control formulation with the distributed feature, we assume:

Assumption 5.

There is at least one satellite in $\mathcal{N}$ has access to the position and velocity of the 0-th satellite described by Equation (7).

Construct a local vertical and local horizontal (LVLH) coordinate frame $O_L-x_L{y}_L{z}_L$ according to position vector ${\overrightarrow{r}}_{e,ir}=[x_{e,ir},y_{e,ir},z_{i,r}]$ and velocity vector ${\overrightarrow{v}}_{e,ir}=[{\dot{x}}_{e,ir},{\dot{y}}_{e,ir},{\dot{z}}_{i,r}]$, where the local $z_L$ axis points towards the Earth center, the local $y_L$ axis perpendicular downwards to the orbital plane and the local $x_L$ axis completes the right hand law from $y_L$ to $z_L$. Without losing generality, the unit axis vectors can be given as ${\overrightarrow{x}}_L={\overrightarrow{y}}_L\times {\overrightarrow{z}}_L,{\overrightarrow{y}}_L=-\frac{{\overrightarrow{r}}_{e,ir}\times {\overrightarrow{v}}_{e,ir}}{\parallel {\overrightarrow{r}}_{e,ir}\times {\overrightarrow{v}}_{e,ir}\parallel },{\overrightarrow{z}}_L=-\frac{{\overrightarrow{r}}_{e,ir}}{\parallel {\overrightarrow{r}}_{e,ir}\parallel }$. The angular velocity of the LVLH frame $O_L-x_L{y}_L{z}_L$ with respect to the ECI frame is:

$${\overrightarrow{\omega }}_{ir}=\frac{{\overrightarrow{r}}_{e,ir}\times {\overrightarrow{v}}_{e,ir}}{r_{e,ir}^2},$$

(8)

where $r_{e,ir}=\parallel {\overrightarrow{r}}_{e,ir}\parallel$. Define the relative position between the i-th satellite and ${\overrightarrow{r}}_{e,ir}$ in $O_L-x_L{y}_L{z}_L$ frame by ${\overrightarrow{r}}_{ir}=({\overrightarrow{r}}_{e,i}-{\overrightarrow{r}}_{e,ir})·{[{\overrightarrow{x}}_L;{\overrightarrow{y}}_L;{\overrightarrow{z}}_L]}^T$, whose derivative is [24]:

$${\overrightarrow{v}}_{ir}\triangleq {\dot{\overrightarrow{r}}}_{ir}={\dot{\overrightarrow{r}}}_{e,i}-{\dot{\overrightarrow{r}}}_{e,ir}-{\overrightarrow{\omega }}_{ir}\times {\overrightarrow{r}}_{ir}.$$

(9)

Then, the associated relative acceleration can be given by:

$$\begin{array}{cc}\hfill {\dot{\overrightarrow{v}}}_{ir}& =-\frac{{\mu }_e}{r_{e,i}^3}{\overrightarrow{r}}_{e,i}+\frac{{\tau }_i}{m_i}+\frac{{\mu }_e}{r_{e,ir}^3}{\overrightarrow{r}}_{e,ir}-{\dot{\overrightarrow{\omega }}}_{ir}\times {\overrightarrow{r}}_{ir}-{\overrightarrow{\omega }}_{ir}\times ({\overrightarrow{\omega }}_{ir}\times {\overrightarrow{r}}_{ir})-2{\overrightarrow{\omega }}_{ir}\times {\overrightarrow{v}}_{ir},\hfill \end{array}$$

(10)

where ${\displaystyle {\dot{\overrightarrow{\omega }}}_{ir}=-2\frac{{\overrightarrow{v}}_{e,ir}{\overrightarrow{r}}_{e,ir}^T}{r_{e,ir}^2}{\overrightarrow{\omega }}_{ir}}$ denotes the angular acceleration.

2.4. Graph Theory

A graph $\mathcal{G}=\{\mathcal{N},\mathcal{E},\mathcal{A}\}$ can be used to describe the interaction among satellites, where $\mathcal{E}\subseteq \mathcal{N}\times \mathcal{N}$ is the edge set and $\mathcal{A}$ denotes the adjacency matrix [26]. An edge $\left\{\left(i,j\right):i\ne j\right\}\in \mathcal{E}$ denotes that the j-th satellite can send information to the i-th satellite via wireless devices. The adjacency matrix is defined by $\mathcal{A}=\left\{a_{ij}\right\}\in {\mathbb{R}}^{n\times n}$, where $a_{ij}=1$ if $\left(i,j\right)\in \mathcal{E}$, otherwise $a_{ij}=0$. Self connection is not allowed, i.e., $a_{ii}=0,\forall i\in \mathcal{N}$. For an undirected graph, $a_{ij}=1\iff a_{ji}=1$ holds, denoting that satellite i and satellite j can exchange information with each other. The path from satellite i to satellite j denotes an edge sequence $\{(i,j_1),(j_2,j_3),\cdots ,(j_{*},j)\}$. The graph $\mathcal{G}$ is called connected if there is a path between any two distinct nodes. The in-degree matrix is given by: $\mathcal{D}=\mathrm{diag}\left\{[l_{11},l_{22},\cdots ,l_{nn}]\right\},l_{ii}=\sum _{j=1}^n{a}_{ij},\forall i,j\in \mathcal{N},$ and the Laplacian matrix is: $\mathcal{L}=\mathcal{D}-\mathcal{A}.$ As reported, the matrix $\mathcal{L}$ is semi-positive definite and has only one zero eigenvalue and $n-1$ positive eigenvalues provided that the graph $\mathcal{G}$ is undirected and connected [26]. Define $a_{i0}=1$ if the there is a valid information flow from the spacecraft to the i-th satellite, otherwise $a_{i,0}=0$. Define:

$$\begin{array}{cc}\hfill \mathcal{H}& =\mathcal{L}+\mathcal{B},\hfill \end{array}$$

(11)

where $\mathcal{B}=\mathrm{diag}\left\{[a_{10},\cdots ,a_{n0}]\right\}$.

Assumption 6.

The graph $\mathcal{G}$ is undirected and connected, and $\mathcal{B}\ne 0_{n\times n}$.

Lemma 1.

([26]). Under Assumption 6, the matrix $\mathcal{H}$ is positive definite.

2.5. Problem Formulation

Until now, the control objective can be formally stated as: given Assumptions 1–6, find a control law ${\tau }_i,\forall i\in \mathcal{N}$ so that the networked satellites described by Equation (4) achieve rendezvous at the orbit trajectory generated by Equation (7), i.e.,

$$\underset{t\to \infty }{lim}{\overrightarrow{r}}_{e,i}-{\overrightarrow{r}}_{e,0}=0_3.$$

(12)

3. Control Design and Stability Analysis

3.1. Estimation Algorithm

Under Assumption 5, the states of the reference orbital vehicle are not available to all satellites, which might obstruct achieving the control objective Equation (12). To overcome this problem, we plan to construct an algorithm to estimate the $[{\overrightarrow{r}}_{e,0},{\overrightarrow{v}}_{e,0}]$ via the neighboring communications among satellites in the ECI frame.

Define ${\overrightarrow{\eta }}_{e,ir}=[x_{e,ir},y_{e,ir},z_{e,ir}]$,${\overrightarrow{\phi }}_{e,ir}=[{\dot{x}}_{e,ir},{\dot{y}}_{e,ir},{\dot{z}}_{e,ir}]$ and ${\overrightarrow{\rho }}_{e,ir}=[{\ddot{x}}_{e,ir},{\ddot{y}}_{e,ir},{\ddot{z}}_{e,ir}]$ as the estimated position, velocity and acceleration of the reference orbital vehicle for the i-th satellite. The updating algorithm for ${\overrightarrow{\eta }}_{e,ir},{\overrightarrow{\phi }}_{e,ir},{\overrightarrow{\rho }}_{e,ir}$ is designed as:

$$\left\{\begin{array}{cc}\hfill {\dot{\overrightarrow{\eta }}}_{e,ir}& ={\overrightarrow{\phi }}_{e,ir}\hfill \\ \hfill {\dot{\overrightarrow{\phi }}}_{e,ir}& ={\overrightarrow{\rho }}_{e,ir}\hfill \\ \hfill {\dot{\overrightarrow{\rho }}}_{e,ir}& =-g_1{\overrightarrow{\phi }}_{e,ir}-g_2{\overrightarrow{\rho }}_{e,ir}-g_3{s}_{e,ir}-g_4\mathrm{sign}{s}_{e,ir}\hfill \end{array}\right.,$$

(13)

where $g_1>0,g_2>0,g_3>0,g_4>\underset{t\ge 0}{sup}∥{\overrightarrow{r}}_{e,0}^{\left(3\right)}+g_1{\ddot{\overrightarrow{r}}}_{e,0}+g_2{\dot{\overrightarrow{r}}}_{e,0}∥$, $\mathrm{sign}(·)$ denotes sign function and

$$\begin{array}{cc}\hfill s_{e,ir}& =\sum _{j=1}^n{a}_{ij}\left({\overrightarrow{\rho }}_{e,ir}-{\overrightarrow{\rho }}_{e,jr}\right)+a_{i0}\left({\overrightarrow{\rho }}_{e,ir}-{\ddot{\overrightarrow{r}}}_{e,0}\right)+g_1\sum _{j=1}^n{a}_{ij}\left({\overrightarrow{\phi }}_{e,ir}-{\overrightarrow{\phi }}_{e,jr}\right)+g_1{a}_{i0}\left({\overrightarrow{\phi }}_{e,ir}-{\dot{\overrightarrow{r}}}_{e,0}\right)\hfill \\ \hfill & +g_2\sum _{j=1}^n{a}_{ij}\left({\overrightarrow{\eta }}_{e,ir}-{\overrightarrow{\eta }}_{e,jr}\right)+g_2{a}_{i0}\left({\overrightarrow{\eta }}_{e,ir}-{\overrightarrow{r}}_{e,0}\right).\hfill \end{array}$$

(14)

Lemma 2.

Given Assumption 6, the gains satisfying $g_1>0,g_2>0,g_3>0,g_4$

$>\underset{t\ge 0}{sup}∥{\overrightarrow{r}}_{e,0}^{\left(3\right)}+g_1{\ddot{\overrightarrow{r}}}_{e,0}+g_2{\dot{\overrightarrow{r}}}_{e,0}∥$ and any initial states, the estimation algorithm given by Equations (13) and (14) ensures that:

$$\underset{t\to \infty }{lim}{\overrightarrow{\eta }}_{e,ir}={\overrightarrow{r}}_{e,0},\underset{t\to \infty }{lim}{\overrightarrow{\phi }}_{e,ir}={\dot{\overrightarrow{r}}}_{e,0},\underset{t\to \infty }{lim}{\overrightarrow{\rho }}_{e,ir}={\ddot{\overrightarrow{r}}}_{e,0}.$$

(15)

Proof. 

Define $3n$-dimensional stacked row vectors,

$$\begin{array}{cc}\hfill {\eta }_{e,r}& ={[{\overrightarrow{\eta }}_{e,1r},\cdots ,{\overrightarrow{\eta }}_{e,nr}]}^T,{\phi }_{e,r}={[{\overrightarrow{\phi }}_{e,1r},\cdots ,{\overrightarrow{\phi }}_{e,nr}]}^T,\hfill \\ \hfill {\rho }_{e,r}& ={[{\overrightarrow{\rho }}_{e,1r},\cdots ,{\overrightarrow{\rho }}_{e,nr}]}^T,s_{e,r}={[s_{e,1r},\cdots ,s_{e,nr}]}^T.\hfill \end{array}$$

(16)

It then, by Equations (13) and (14) follows that:

$$\left\{\begin{array}{cc}\hfill {\dot{\eta }}_{e,r}& ={\phi }_{e,r}\hfill \\ \hfill {\dot{\phi }}_{e,r}& ={\rho }_{e,r}\hfill \\ \hfill {\dot{\rho }}_{e,r}& =-g_1{\phi }_{e,r}-g_2{\rho }_{e,r}-g_3{s}_{e,r}-g_4{P}_e{s}_{e,r}\hfill \end{array}\right.,$$

(17)

where

$$P_e=\mathrm{diag}\{\mathrm{sign}{s}_{e,1r},\cdots ,\mathrm{sign}{s}_{e,nr}\}\in {\mathbb{R}}^{3n\times 3n}.$$

(18)

Define estimation error

$${\epsilon }_{e,r}={\overrightarrow{\eta }}_{e,r}-{\mathbf{1}}_n\otimes {\overrightarrow{r}}_{e,0},$$

(19)

where ⊗ denotes the Kronecker product. Within Equation (19), the $s_{e,r}$ in Equation (16) can be arranged as:

$$s_{e,r}=\mathcal{H}\otimes I_3({\ddot{\epsilon }}_{e,r}+g_2{\dot{\epsilon }}_{e,r}+g_1{\epsilon }_{e,r}),$$

(20)

with time-derivative along the trajectory of Equation (14) being given by:

$$\begin{array}{cc}\hfill {\dot{s}}_{e,r}& =H\otimes I_3\left({\stackrel{⃛}{\overrightarrow{r}}}_{e,r}+g_2{\ddot{\epsilon }}_{e,r}+g_1{\dot{\epsilon }}_{e,r}\right)\hfill \\ \hfill & =H\otimes I_3\left[-g_3{s}_{e,r}-g_4{P}_e{s}_{e,r}-1_n\otimes \left({\stackrel{⃛}{\overrightarrow{r}}}_{e,0}+g_2{\ddot{\overrightarrow{r}}}_{e,0}+g_1{\dot{\overrightarrow{r}}}_{e,0}\right)\right].\hfill \end{array}$$

(21)

Given Assumption 6 and Lemma 1, we know that the matrix $\mathcal{H}$ is positive definite. Choose a Lyapunov candidate as:

$$V_1=0.5s_{e,r}^T{(\mathcal{H}\otimes I_3)}^{-1}{s}_{e,r},$$

(22)

which, by the relationship $\lambda (\mathcal{H}\otimes I_3)=\lambda \left(\mathcal{H}\right)$, leads to:

$$\frac{1}{2{\lambda }_M\left(\mathcal{H}\right)}\parallel s_{e,r}{\parallel }^2\le V_1\le \frac{1}{2{\lambda }_m\left(\mathcal{H}\right)}{\parallel s_{e,r}\parallel }^2.$$

(23)

Then, calculate the time-derivative of $V_1$ as:

$$\begin{array}{cc}\hfill {\dot{V}}_1& =-g_3{s}_{e,r}^T{s}_{e,r}-g_4{s}_{e,r}^T{P}_e{s}_{e,r}-s_{e,r}^T{1}_n\otimes \left({\stackrel{⃛}{\overrightarrow{r}}}_{e,0}+g_2{\ddot{\overrightarrow{r}}}_{e,0}+g_1{\dot{\overrightarrow{r}}}_{e,0}\right)\hfill \\ \hfill & =-g_3{s}_{e,r}^T{s}_{e,r}-\sum _{i=1}^n\left[g_4{s}_{e,ir}\mathrm{sign}{s}_{e,ir}^T+s_{e,ir}{\left({\stackrel{⃛}{\overrightarrow{r}}}_{e,0}+g_2{\ddot{\overrightarrow{r}}}_{e,0}+g_1{\dot{\overrightarrow{r}}}_{e,0}\right)}^T\right]\hfill \end{array}$$

(24)

Let $\vartheta \stackrel{\Delta }{=}\underset{t\ge 0}{sup}\left\{∥{\stackrel{⃛}{\overrightarrow{r}}}_{e,0}+g_2{\ddot{\overrightarrow{r}}}_{e,0}+g_1{\dot{\overrightarrow{r}}}_{e,0}∥\right\}$. The Equation (24) satisfies:

$$\begin{array}{cc}\hfill {\dot{V}}_1& \le -g_3{∥s_{e,r}∥}^2-\sum _{i=1}^n\left|s_{e,ir}\right|\left(g_4-\vartheta \right)\hfill \\ \hfill & \le -g_3{∥s_{e,r}∥}^2.\hfill \end{array}$$

(25)

Therefore, $s_{e,r}$ would decay to zero with exponential decaying rate $-0.5g_3$ [27]. We then transform Equation (20) into

$${\ddot{\epsilon }}_{e,r}+g_2{\dot{\epsilon }}_{e,r}+g_1{\epsilon }_{e,r}={(\mathcal{H}\otimes I_3)}^{-1}{s}_{e,r}.$$

(26)

The vectors ${\epsilon }_{e,r}$,${\dot{\epsilon }}_{e,r}$ and ${\ddot{\epsilon }}_{e,r}$ will converge to zero exponentially due to $s_{e,r}$ decays to zero exponentially and the matrix $\mathcal{H}\otimes I_3$ is positive definite [28]. By definition Equation (19), the claims in the lemma follows. This completes the proof. □

Via the information exchange with connected satellites, all satellites achieve the estimate of $[{\overrightarrow{r}}_{e,0},{\overrightarrow{v}}_{e,0}]$. One should notice that the generated velocity ${\overrightarrow{\phi }}_{e,ir}$ and acceleration ${\overrightarrow{\rho }}_{e,ir}$ are bounded because the associated components of the orbital reference vehicle are bounded, and the exponential convergence to zero of the estimation errors. The generated signals ${\overrightarrow{\eta }}_{e,ir},{\overrightarrow{\phi }}_{e,ir}$ and ${\overrightarrow{\rho }}_{e,ir}$ stated in Lemma 2 will be viewed as a reference signal for the i-th satellite for $i\in \mathcal{N}$. As can be seen, the efficacy of the proposed estimation algorithm relies only on partial access to $[{\overrightarrow{r}}_{e,0},{\overrightarrow{v}}_{e,0}]$ and undirected connected communications, which gains a more robust property of the group as a whole when compared with the centralized ones [17,18,19,20,21,22,23] that require all satellites in the group to know the information of the reference target.

3.2. Feedback Control

In this subsection, we will take the modeling uncertainty into account and design a control law for each satellite $i\in \mathcal{N}$ so as to track the reference signal generated by Equations (13) and (14).

Define tracking errors by:

$${\varpi }_{e,i}=[{\overrightarrow{r}}_{e,i}-{\overrightarrow{\eta }}_{e,ir}]·{\left[\begin{array}{c}\frac{{\overrightarrow{\eta }}_{e,ir}\times {\overrightarrow{\phi }}_{e,ir}}{\parallel {\overrightarrow{\eta }}_{e,ir}\times {\overrightarrow{\phi }}_{e,ir}\parallel }\times \frac{{\overrightarrow{\eta }}_{e,ir}}{\parallel {\overrightarrow{\eta }}_{e,ir}\parallel }\hfill \\ -\frac{{\overrightarrow{\eta }}_{e,ir}\times {\overrightarrow{\phi }}_{e,ir}}{\parallel {\overrightarrow{\eta }}_{e,ir}\times {\overrightarrow{\phi }}_{e,ir}\parallel }\hfill \\ -\frac{{\overrightarrow{\eta }}_{e,ir}}{\parallel {\overrightarrow{\eta }}_{e,ir}\parallel }\hfill \end{array}\right]}^T.$$

(27)

Following the derivations Equations (8)–(10), the first- and second-order derivatives of Equation (27) can be calculated as:

$$\begin{array}{cc}\hfill {\dot{\varpi }}_{e,i}& ={\dot{\overrightarrow{r}}}_{e,i}-{\phi }_{e,ir}-{\zeta }_{e,ir}\times {\varpi }_{e,i},\hfill \\ \hfill {\ddot{\varpi }}_{e,i}& =-\frac{{\mu }_e}{r_{e,i}^3}{\overrightarrow{r}}_{e,i}+\frac{{\tau }_i}{m_i\left(t\right)}+\frac{{\mu }_e}{{∥{\eta }_{e,ir}∥}^3}{\eta }_{e,ir}-{\dot{\zeta }}_{e,ir}\times {\varpi }_{e,i}-{\zeta }_{e,ir}\times ({\zeta }_{e,ir}\times {\varpi }_{e,i})-2{\zeta }_{e,ir}\times {\dot{\varpi }}_{e,i},\hfill \end{array}$$

(28)

where ${\displaystyle {\zeta }_{e,ir}=\frac{{\eta }_{e,ir}\times {\phi }_{e,ir}}{\parallel {\eta }_{e,ir}\parallel }}$ and ${\displaystyle {\dot{\zeta }}_{e,ir}=-2\frac{{\phi }_{e,ir}{\eta }_{e,ir}^T}{r_{e,ir}^2}\zeta }$ represent the angular velocity and angular acceleration of the frame formed by the reference signal $[{\overrightarrow{\eta }}_{e,ir},{\overrightarrow{\phi }}_{e,ir}]$.

Define an intermediate variable:

$$\begin{array}{cc}\hfill {\delta }_{e,i}& =-\frac{{\mu }_e}{r_{e,i}^3}{\overrightarrow{r}}_{e,i}+\frac{{\mu }_e}{{∥{\eta }_{e,ir}∥}^3}{\eta }_{e,ir}-{\dot{\zeta }}_{e,ir}\times {\varpi }_{e,i}-{\zeta }_{e,ir}\times ({\zeta }_{e,ir}\times {\varpi }_{e,i})\hfill \\ \hfill & -2{\zeta }_{e,ir}\times {\dot{\varpi }}_{e,i}+k_{3}tanh\left({\dot{\varpi }}_{e,i}+k_2{\varpi }_{e,i}\right)+k_{1}tanh\left({\dot{\varpi }}_{e,i}\right),\hfill \end{array}$$

(29)

where $k_1,k_2,k_3>0$. We then propose the following control law with a parameter updating algorithm,

$$\left\{\begin{array}{cc}\hfill {\tau }_i& =-{\widehat{m}}_i\left(t\right){\delta }_{e,i}\hfill \\ \hfill {\dot{\widehat{m}}}_i& ={\dot{m}}_i\left(t\right)+k_4\left[k_{3}tanh\left({\dot{\varpi }}_{e,i}+k_2{\varpi }_{e,i}\right)+k_{1}tanh{\dot{\varpi }}_{e,i}+k_2{\dot{\varpi }}_{e,i}\right]{\delta }_{e,i}^T\hfill \end{array}\right.,$$

(30)

where $k_4>0$ and ${\widehat{m}}_i\left(t\right)$ denotes the estimated value for $m_i$. The estimated value ${\widehat{m}}_i\left(t\right)$ is updated via the integration ${\widehat{m}}_i\left(t\right)={\widehat{m}}_i\left(0\right)+{\int }_0^t{\dot{\widehat{m}}}_i\left(\chi \right)\mathrm{d}\chi$.

Theorem 1.

If the control gains are selected such that $k_1>0,k_2>0,k_3>0,k_4>0$, the adaptive control law shown in Equation (30) guarantees that the tracking error ${\varpi }_{e,i}$ converges to zero globally asymptotically. In addition, the control objective Equation (12) is achieved.

Proof. 

Substituting Equation (29) into ${\ddot{\varpi }}_{e,i}$ in Equation (28) results in:

$${\ddot{\varpi }}_{e,i}=-k_{3}tanh\left({\dot{\varpi }}_{e,i}+k_2{\varpi }_{e,i}\right)-k_{1}tanh\left({\dot{\varpi }}_{e,i}\right)+\frac{{\tau }_i}{m_i\left(t\right)}+{\delta }_{e,i}.$$

(31)

Define estimation error as ${\overline{m}}_i\left(t\right)=m_i\left(t\right)-{\widehat{m}}_i\left(t\right)$ and taking the control law Equation (30) into Equation (31), we obtain:

$${\ddot{\varpi }}_{e,i}=-k_{3}tanh\left({\dot{\varpi }}_{e,i}+k_2{\varpi }_{e,i}\right)-k_{1}tanh\left({\dot{\varpi }}_{e,i}\right)+\frac{{\overline{m}}_i}{m_i\left(t\right)}{\delta }_{e,i}.$$

(32)

Choose a positive definite function $V_2=V_2({\varpi }_i,{\dot{\varpi }}_i,{\overline{m}}_i,t)$ as follows:

$$\begin{array}{cc}\hfill V_2& =k_3{m}_i\ln \cosh \left({\dot{\varpi }}_{e,i}+k_2{\varpi }_{e,i}\right)+k_1{m}_i\ln \cosh {\dot{\varpi }}_{e,i}+0.5k_2{m}_i{\dot{\varpi }}_{e,i}{\dot{\varpi }}_{e,i}^T+\frac{{\overline{m}}_i^2}{2k_4}.\hfill \end{array}$$

(33)

The time-derivative of Equation (33) along the solution trajectory of Equation (32) can then be calculated as:

$$\begin{array}{cc}\hfill {\dot{V}}_2& ={\dot{m}}_i\left[k_3\ln \cosh \left({\dot{\varpi }}_{e,i}+k_2{\varpi }_{e,i}\right)+k_1\ln \cosh {\dot{\varpi }}_{e,i}+0.5k_2{\dot{\varpi }}_{e,i}{\dot{\varpi }}_{e,i}^T\right]\hfill \\ \hfill & +k_3{m}_{i}tanh\left({\dot{\varpi }}_{e,i}+k_2{\varpi }_{e,i}\right){\left({\ddot{\varpi }}_{e,i}+k_2{\dot{\varpi }}_{e,i}\right)}^T+k_1{m}_{i}tanh{\dot{\varpi }}_{e,i}{\ddot{\varpi }}_{e,i}^T+k_2{m}_i{\dot{\varpi }}_{e,i}{\ddot{\varpi }}_{e,i}^T\hfill \\ \hfill & +\frac{{\overline{m}}_i}{k_4}\left({\dot{m}}_i-{\dot{\widehat{m}}}_i\right).\hfill \end{array}$$

(34)

As ${\dot{m}}_i\le 0,\forall t\ge 0$, the first row of Equation (34) is non-positive and hence the ${\dot{V}}_2$ satisfies:

$$\begin{array}{cc}\hfill {\dot{V}}_2& \le k_3{m}_{i}tanh\left({\dot{\varpi }}_{e,i}+k_2{\varpi }_{e,i}\right){\left({\ddot{\varpi }}_{e,i}+k_2{\dot{\varpi }}_{e,i}\right)}^T+k_1{m}_{i}tanh{\dot{\varpi }}_{e,i}{\ddot{\varpi }}_{e,i}^T+k_2{m}_i{\dot{\varpi }}_{e,i}{\ddot{\varpi }}_{e,i}^T\hfill \\ \hfill & +\frac{{\overline{m}}_i}{k_4}\left({\dot{m}}_i-{\dot{\widehat{m}}}_i\right).\hfill \end{array}$$

(35)

Substitute ${\ddot{\varpi }}_{e,i}$ of Equation (32) into Equation (35),

$$\begin{array}{cc}\hfill {\dot{V}}_2& \le m_i\left(t\right)\left\{-k_3^{2}tanh\left({\dot{\varpi }}_{e,i}+k_2{\varpi }_{e,i}\right){\left[tanh\left({\dot{\varpi }}_{e,i}+k_2{\varpi }_{e,i}\right)\right]}^T-k_1^{2}tanh{\dot{\varpi }}_{e,i}tanh{\left({\dot{\varpi }}_{e,i}\right)}^T\right\}\hfill \\ \hfill & +m_i\left(t\right)\left\{-k_1{k}_2{\dot{\varpi }}_{e,i}tanh{\left({\dot{\varpi }}_{e,i}\right)}^T-2k_3{k}_{1}tanh\left({\dot{\varpi }}_{e,i}+k_2{\varpi }_{e,i}\right)tanh{\left({\dot{\varpi }}_{e,i}\right)}^T\right\}\hfill \\ \hfill & +\left[k_{3}tanh\left({\dot{\varpi }}_{e,i}+k_2{\varpi }_{e,i}\right)+k_{1}tanh{\dot{\varpi }}_{e,i}+k_2{\dot{\varpi }}_{e,i}\right]{\overline{m}}_i{\delta }_{e,i}^T+\frac{{\overline{m}}_i}{k_4}\left({\dot{m}}_i-{\dot{\widehat{m}}}_i\right).\hfill \end{array}$$

(36)

Taking the updating algorithm for ${\widehat{m}}_i$ in Equation (30) into Equation (36) and performing some direct calculations, we obtain:

$$\begin{array}{cc}\hfill {\dot{V}}_2& \le -m_i\left(t\right){∥k_{3}tanh\left({\dot{\varpi }}_{e,i}+k_2{\varpi }_{e,i}\right)+k_{1}tanh{\dot{\varpi }}_{e,i}∥}^2-k_1{k}_2{m}_i\left(t\right){\dot{\varpi }}_{e,i}tanh{\left({\dot{\varpi }}_{e,i}\right)}^T\hfill \\ \hfill & \le -{\underline{m}}_i{∥k_{3}tanh\left({\dot{\varpi }}_{e,i}+k_2{\varpi }_{e,i}\right)+k_{1}tanh{\dot{\varpi }}_{e,i}∥}^2-k_1{k}_2{\underline{m}}_i{\dot{\varpi }}_{e,i}tanh{\left({\dot{\varpi }}_{e,i}\right)}^T\hfill \\ \hfill & \le 0.\hfill \end{array}$$

(37)

It is obvious that ${\dot{V}}_2$ is semi-negative definite, which means that ${\varpi }_{e,i},{\dot{\varpi }}_{e,i}$ and ${\overline{m}}_i$ are bounded. With direct calculation, one would find out that ${\delta }_{e,i}$ is bounded. Henceforth, ${\ddot{\varpi }}_{e,i}$ is bounded. All these bounded variables demonstrate that:

$$\begin{array}{cc}\hfill {\ddot{V}}_2& ={\ddot{m}}_i\left(t\right)\left[k_3\ln \cosh \left({\dot{\varpi }}_{e,i}+k_2{\varpi }_{e,i}\right)+k_{1}lncosh{\dot{\varpi }}_{e,i}+0.5k_2{\dot{\varpi }}_{e,i}{\dot{\varpi }}_{e,i}^T\right]\hfill \\ \hfill & -{\dot{m}}_i\left(t\right){∥k_{3}tanh\left({\dot{\varpi }}_{e,i}+k_2{\varpi }_{e,i}\right)+k_{1}tanh{\dot{\varpi }}_{e,i}∥}^2-k_1{k}_2{\dot{m}}_i\left(t\right){\dot{\varpi }}_{e,i}tanh{\left({\dot{\varpi }}_{e,i}\right)}^T\hfill \\ \hfill & -2m_i\left(t\right)\left[k_{3}tanh\left({\dot{\varpi }}_{e,i}+k_2{\varpi }_{e,i}\right)+k_{1}tanh{\dot{\varpi }}_{e,i}\right]\times {\left[\frac{k_3{\ddot{\varpi }}_{e,i}+k_3{k}_2{\dot{\varpi }}_{e,i}}{1+{∥{\dot{\varpi }}_{e,i}+k_2{\varpi }_{e,i}∥}^2}+\frac{k_1{\ddot{\varpi }}_{e,i}}{1+{∥{\dot{\varpi }}_{e,i}∥}^2}\right]}^T\hfill \\ \hfill & -k_1{k}_2{m}_i\left(t\right){\ddot{\varpi }}_{e,i}tanh{\left({\dot{\varpi }}_{e,i}\right)}^T-k_1{k}_2{m}_i\left(t\right)\frac{{\dot{\varpi }}_{e,i}{\ddot{\varpi }}_{e,i}^T}{1+{∥{\dot{\varpi }}_{e,i}∥}^2}\hfill \end{array}$$

(38)

is bounded. The boundedness of ${\ddot{V}}_2$ implies that ${\dot{V}}_2$ is uniformly continuous, which, according to Barbalat’s lemma [27], proves that ${\varpi }_{e,i}$ and ${\dot{\varpi }}_{e,i}$ converge to zero globally asymptotically. Rewrite Equation (27) as follows:

$${\overrightarrow{r}}_{e,i}-{\overrightarrow{r}}_{e,0}={\varpi }_{e,i}+({\overrightarrow{\eta }}_{e,ir}-{\overrightarrow{r}}_{e,0}).$$

(39)

Because $\underset{t\to \infty }{lim}{\varpi }_{e,i}=0$ and $\underset{t\to \infty }{lim}\left({\overrightarrow{\eta }}_{e,ir}-{\overrightarrow{r}}_{e,0}\right)=0$ hold simultaneously, we use a simple theory of a cascade system [28] and obtain that $\underset{t\to \infty }{lim}\left({\overrightarrow{r}}_{e,i}-{\overrightarrow{r}}_{e,0}\right)=0$. This completes the proof. □

The above analysis shows that the satellites in $\mathcal{N}$ would achieve rendezvous at the reference orbital vehicle. One might note that the safe orbit height is not taken into account, i.e., the satellite must be orbiting above a certain height with respect to the Earth’s surface. We here propose an initial solution for this issue without presenting the detailed proof. Choose two positive numbers $h_1$ and $h_2$ satisfying $h_2>h_1>0$, define a continuous and derivable function $f\left(x\right):(h_1,+\infty )\mapsto {\mathbb{R}}_{\ge 0}$ as:

$$f\left(x\right)=\left\{\begin{array}{lll}-{\alpha }_{1}ln{\alpha }_2\left(x-h_1\right)+{\alpha }_3& ,& \forall h_1<x\le 0.5\left(h_1+h_2\right);\\ {\alpha }_4{\left(x-h_2\right)}^2\hfill & ,& \forall 0.5\left(h_1+h_2\right)<x\le h_2;\\ 0\hfill & ,& \forall h_2<x,\end{array}\right.$$

(40)

where the constants are chosen as ${\displaystyle {\alpha }_1=0.5{\alpha }_4{\left(h_2-h_1\right)}^2,{\alpha }_2=\frac{1}{2\left(h_2-h_1\right)},{\alpha }_3=0.5{\alpha }_1,{\alpha }_4>0}$. It is direct to verify that $f\left(x\right)$ is derivable and continuous for all $x\in (h_1,+\infty )$, so we omit the proof here. To ensure the orbital satellites orbiting above the $h_1$ and $\parallel {\overrightarrow{r}}_{e,i}\parallel >h_1,\forall t\ge 0$, we borrow the basic idea from an artificial potential field approach [29] and modify the control law Equation (30) into:

$$\left\{\begin{array}{cc}\hfill {\tau }_i& =-{\widehat{m}}_i{\delta }_{e,i}+f(\parallel {\overrightarrow{r}}_{e,i}\parallel )\frac{{\overrightarrow{r}}_{e,i}}{\parallel {\overrightarrow{r}}_{e,i}\parallel }{\tau }_{i,c}\hfill \\ \hfill {\dot{\widehat{m}}}_i& ={\dot{m}}_i\left(t\right)+k_4\left[k_{3}tanh\left({\dot{\varpi }}_{e,i}+k_2{\varpi }_{e,i}\right)+k_{1}tanh{\dot{\varpi }}_{e,i}+k_2{\dot{\varpi }}_{e,i}\right]{\delta }_{e,i}^T\hfill \end{array}\right.,$$

(41)

with ${\tau }_{i,c}$ being a constant. The introduced term $f(\parallel {\overrightarrow{r}}_{e,i}\parallel )\frac{{\overrightarrow{r}}_{e,i}}{\parallel {\overrightarrow{r}}_{e,i}\parallel }{\tau }_{i,c}$ would produce a reversing force in the direction ${\overrightarrow{r}}_{e,i}$, pointing outwards from the Earth’s center, so that the $\parallel {\overrightarrow{r}}_{e,i}\parallel >h_1,\forall t\ge 0$. In addition, as can be seen from Equations (40) and (41), the term becomes efficient only when $\parallel {\overrightarrow{r}}_{e,i}\parallel \in (h_1,h_2)$. This means that the modified control law Equation (41) reduces to Equation (30) if $\parallel {\overrightarrow{r}}_{e,i}\parallel \ge h_2$.

4. Numerical Simulation

This section presents the numerical simulation to validate the proposed control design. Two parts are involved. First, we discuss how to propagate the projection vector ${\overrightarrow{s}}_{i,j}$ of the satellite position and calculate the satellite position for feedback control. Second, we validate the leader estimation algorithm and feed the disturbed state signals into the control loop to verify the effectiveness of the adaptive control algorithm.

4.1. Pulsar-Based Positioning and Simulation Setup

For simulation use, we select three pulsars as positioning sources. The pulsar positions are listed in

Table 1. Pulsar Position in Galactic Frame [30].

j	Name	Longitude	Latitude
1	PSR B0531+21	184.56${}^{\circ }$	−5.78${}^{\circ }$
2	PSR B1821−24	7.8${}^{\circ }$	−5.58${}^{\circ }$
3	PSR B1937+21	57.51${}^{\circ }$	−0.29${}^{\circ }$

and

Table 2. Pulsars Information in SCI Frame.

j	Name	Right Ascension	Declination
1	PSR B0531+21	84.102438550${}^{\circ }$	−1.29446370${}^{\circ }$
2	PSR B1821−24	275.56845533${}^{\circ }$	−1.54719148${}^{\circ }$
3	PSR B1937+21	301.97479547${}^{\circ }$	42.29726047${}^{\circ }$

. The galactic coordinates shown in Table 1 are transformed into SCI coordinates in Table 2 via the online coordinate calculator provided by https://ned.ipac.caltech.edu/coordinate_calculator, accessed on 11 November 2021.

Within the pulsar position in Table 2, we use Equation (2) to calculate the fixed vectors formed by pulsar radiations in the SCI frame as follows:

$$\begin{array}{cc}\hfill {\overrightarrow{q}}_1& =[0.1027,-0.0023,-0.0226],\hfill \\ \hfill {\overrightarrow{q}}_2& =[0.0970,-0.0026,-0.0270],\hfill \\ \hfill {\overrightarrow{q}}_3& =[0.3917,0.3564,0.6730].\hfill \end{array}$$

(42)

As discussed before, the projection of satellite position vector ${\overrightarrow{r}}_{s,i},i\in \mathcal{N}$ in the SCI frame on ${\overrightarrow{q}}_j,j=\{1,2,3\}$ is measurable. Given three pulsars, i.e., $j=3$, we define the measured row vectors for the i-th satellite by ${\overrightarrow{s}}_{i,1},{\overrightarrow{s}}_{i,2},{\overrightarrow{s}}_{i,3}$, which can be propagated by the following equations:

$${\overrightarrow{s}}_{i,j}=[({\overrightarrow{r}}_{e,i}+{\overrightarrow{r}}_{s,e})·{\overrightarrow{q}}_j]{\overrightarrow{q}}_j,j=\{1,2,3\},$$

(43)

where ${\overrightarrow{r}}_{e,i}$ and ${\overrightarrow{r}}_{s,e}$ are solved by the associated Keplerian equations of motion. Then we use the two-dimensional plane given in Equation (3) and obtain the following linear equations:

$$\left[\begin{array}{c}{\overrightarrow{s}}_{i,1}\\ {\overrightarrow{s}}_{i,2}\\ {\overrightarrow{s}}_{i,3}\end{array}\right]{\overrightarrow{r}}_{s,i}^T=[\parallel {\overrightarrow{s}}_{i,1}{\parallel }^2,\parallel {\overrightarrow{s}}_{i,2}{\parallel }^2,\parallel {\overrightarrow{s}}_{i,3}{{\parallel }^2]}^T.$$

(44)

Equation (44) suggests that the satellite position in the SCI frame using pulsar measurements can be solved by:

$${\overrightarrow{r}}_{s,i}={\left[{\left[\begin{array}{c}{\overrightarrow{s}}_{i,1}\\ {\overrightarrow{s}}_{i,2}\\ {\overrightarrow{s}}_{i,3}\end{array}\right]}^{-1}[\parallel {\overrightarrow{s}}_{i,1}{\parallel }^2,\parallel {\overrightarrow{s}}_{i,2}{\parallel }^2,\parallel {\overrightarrow{s}}_{i,3}{{\parallel }^2]}^T\right]}^T.$$

(45)

To investigate the robust property of the proposed algorithm regarding the disturbance of position and velocity measurement, we inject the unknown time-varying disturbance into the state variables for feedback control. The overall process can be summarized as the figure shown below (Figure 2), in which the satellite ECI coordinate II is fed into the control loop.

4.2. Estimate of ${\overrightarrow{r}}_{e,0}$ and ${\overrightarrow{v}}_{e,0}$

For validating the proposed algorithm, we suppose that four orbital satellites transporting supplies to the international space station (ISS) are indexed by 0, which accords with the rendezvous scenario. The initial orbital elements of ISS and satellites are listed in

Table 3. Initial Orbital Elements.

i	A	ECC	INCL	RAAN	AANP	APP
0	6792 km	0.005426	51.6438${}^{\circ }$	38.8886${}^{\circ }$	23.0560${}^{\circ }$	63${}^{\circ }$
1	6881 km	0.006340	50.3210${}^{\circ }$	40.0100${}^{\circ }$	20.2022${}^{\circ }$	60${}^{\circ }$
2	6922 km	0.005924	60.5380${}^{\circ }$	39.4500${}^{\circ }$	25.1991${}^{\circ }$	64${}^{\circ }$
3	7238 km	0.007020	52.6225${}^{\circ }$	43.1526${}^{\circ }$	28.5234${}^{\circ }$	58${}^{\circ }$
4	7055 km	0.009070	40.4819${}^{\circ }$	42.5128${}^{\circ }$	23.4040${}^{\circ }$	55${}^{\circ }$

Where the definitions for letters in Table 3 are listed below, A: semi-major axis, ECC: orbit eccentricity, INCL: Inclination, RAAN: right ascension of ascending node, AANP: angle between ascending node and periapsis, APP: angle between periapsis and the satellite.

.

According to the initial orbital elements in Table 3, we calculate ${\overrightarrow{r}}_{e,i}\left(0\right)$ and ${\overrightarrow{v}}_{e,i}\left(0\right),\forall i\in \mathcal{N}$ as follows:

$$\begin{array}{cc}\hfill {\overrightarrow{r}}_{e,1}\left(0\right)& =[-1880.5,4055.8,5201.9]\mathrm{km},{\overrightarrow{v}}_{e,1}\left(0\right)=[-6.3076,-4.1772,1.0317]\mathrm{km}/\mathrm{s},\hfill \\ \hfill {\overrightarrow{r}}_{e,2}\left(0\right)& =[-2082.9,2683.1,6010.5]\mathrm{km},{\overrightarrow{v}}_{e,2}\left(0\right)=[-5.9198,-4.7778,0.1278]\mathrm{km}/\mathrm{s},\hfill \\ \hfill {\overrightarrow{r}}_{e,3}\left(0\right)& =[-2669.4,3486.7,5719.6]\mathrm{km},{\overrightarrow{v}}_{e,3}\left(0\right)=[-5.6280,-4.8637,0.3940]\mathrm{km}/\mathrm{s},\hfill \\ \hfill {\overrightarrow{r}}_{e,4}\left(0\right)& =[-2493.6,4807.7,4463.1]\mathrm{km},{\overrightarrow{v}}_{e,4}\left(0\right)=[-6.2565,-4.1120,1.0216]\mathrm{km}/\mathrm{s}.\hfill \end{array}$$

The initial orbits of the satellites are depicted in Figure 3. As can be seen, the satellite $i,\forall i\in \mathcal{N}$ and the ISS are orbiting the Earth with independent orbits.

For achieving distributed estimation about $[{\overrightarrow{r}}_{e,0},{\overrightarrow{v}}_{r,0}]$, the satellites in $\mathcal{N}$ are supposed to exchange information via the communication topology shown in Figure 4.

The coefficients are chosen as $g_1=1,g_2=0,1,g_3=1,g_{4x}=g_{4y}=g_{4z}=1.2$. Set initial values by ${\overrightarrow{\eta }}_{e,ir}\left(0\right)={\overrightarrow{r}}_{e,i}\left(0\right),{\overrightarrow{\phi }}_{e,ir}\left(0\right)={\overrightarrow{v}}_{e,i}\left(0\right)$ and ${\overrightarrow{\rho }}_{e,ir}\left(0\right)=0_3$. The simulation time length is the same as the next subsection. However, we only plot the error trajectory in the time interval $t\in [0,300]$ for showing the convergence clearly. The simulation results are depicted in Figure 5. It can be found out that all estimation errors converge to zero asymptotically. Within $t\ge 50 \mathrm{s}$, the signals ${\overrightarrow{\eta }}_{e,1r},\forall i\in \mathcal{N}=\{1,2,3,4\}$ feature stable shape without any chattering phenomenon occurring.

4.3. Control Validation

For achieving cooperative tasks, we plan to steer all satellites to the orbit ${\overrightarrow{r}}_{e,0}$ via the proposed control law and the pulsar-based positioning method. Within the initial condition shown in Table 3, we set the control coefficients by $k_1=0.2,k_2=0.4,k_3=0.2,k_4=0.01$. The real satellite masses are set as $m_1=450\mathrm{kg},m_2=423\mathrm{kg},m_3=467\mathrm{kg}$ and $m_4=433\mathrm{kg}$. These masses are unknown for each satellite in $\mathcal{N}$. The ratio of change of mass is chosen as ${\alpha }_i=-0.1,\forall i\in \{1,2,3,4\}$. The initial estimated mass values are set as ${\widehat{m}}_i\left(0\right)$ = (1 + 6%)$m_i$. Set the Sun gravitational constant as ${\mu }_s=1.32712439935\times {10}^{11}{\mathrm{kg}}^3/{\mathrm{s}}^2$ and the Earth gravitational constant as ${\mu }_e=3.986005\times {10}^5{\mathrm{kg}}^3/{\mathrm{s}}^2$. The semi-major axis of Earth rotating round the Sun is chosen as $a_{es}=1.496\times {10}^8\mathrm{km}$ within eccentricity given by $ec{c}_e=0.0167$. The initial anomaly of the Earth is set as zero. The Earth’s radius is set as $R_e=6371$ (km). Accordingly, we set $h_1=R_e+300\left(\mathrm{km}\right)=6671\left(\mathrm{km}\right)$, $h_2=R_e+350\left(\mathrm{km}\right)=6721\left(\mathrm{km}\right)$, ${\alpha }_1=125,{\alpha }_2=0.04,{\alpha }_3=62.5$, ${\alpha }_4=0.1$ and ${\tau }_{i,c}=100$. This setting means that we want all satellites are always orbiting 300 km above the Earth’s surface. In addition, the position and velocity measurements used for feedback control are supposed to suffer from disturbances ${\overrightarrow{\Delta }}_{e,r,i}$ and ${\overrightarrow{\Delta }}_{e,v,i}$ respectively, defined as follows:

$$\begin{array}{cc}\hfill {\overrightarrow{\Delta }}_{e,r,1}& =50[sin0.01t,cos0.05t,sin0.02t]\left(\mathrm{m}\right),{\overrightarrow{\Delta }}_{e,r,2}=50[sin0.03t,sin0.04t,cos0.02t]\left(\mathrm{m}\right),\hfill \\ \hfill {\overrightarrow{\Delta }}_{e,r,3}& =50[cos0.02t,sin0.01t,cos0.03t]\left(\mathrm{m}\right),{\overrightarrow{\Delta }}_{e,r,4}=50[sin0.015t,cos0.02t,sin0.012t]\left(\mathrm{m}\right),\hfill \\ \hfill {\overrightarrow{\Delta }}_{e,v,1}& =2[sin0.023t,cos0.03t,sin0.012t](\mathrm{m}/\mathrm{s}),{\overrightarrow{\Delta }}_{e,v,2}=2[sin0.016t,sin0.035t,cos0.022t](\mathrm{m}/\mathrm{s}),\hfill \\ \hfill {\overrightarrow{\Delta }}_{e,v,3}& =2[cos0.026t,sin0.013t,cos0.03t](\mathrm{m}/\mathrm{s}),{\overrightarrow{\Delta }}_{e,v,4}=2[sin0.018t,cos0.02t,sin0.032t](\mathrm{m}/\mathrm{s}).\hfill \end{array}$$

More precisely, the position measurement errors are fluctuated within $\pm 75.7\mathrm{m}$, while the fluctuation for velocity measurement is ±3.4 m/s. We set simulation time as $5000\mathrm{s}$ and start the control update after $t\ge 100\mathrm{s}$ according to the convergence of the estimation errors stated previously. The simulation results are depicted in Figure 6, Figure 7, Figure 8 and Figure 9.

As depicted in Figure 6a, all satellites in $\mathcal{N}$ move to the orbit of ${\overrightarrow{r}}_{e,0}$, achieving a rendezvous control scenario from initial orbits. In Figure 6a, the Earth is rotating around the Sun while the satellites are rotating with the Earth. The collective helix curves demonstrate that the control algorithm is effective even though the controller Equation (30) is developed in ECI coordinate frame. Figure 6b further supports this point of view and demonstrates that the position rendezvous errors converge into small regions enclosing the zero.

Due to position and velocity measurement disturbances during the transient phase, the rendezvous errors are not steered converging to zero. We plot the final relative position and velocity with respect to $[{\overrightarrow{r}}_{e,0},{\overrightarrow{v}}_{e,0}]$ in Figure 7. The simulation resulting in the steady-state of $t\ge 3000$ s shows that the final relative position is maintained to be less than 90 m and the relative velocity to be less than 3 m/s. More precisely, the final relative position and velocity mean values are calculated and summarized in

Table 4. The mean values of final relative states for $t\ge 3000$ s.

i	1	2	3	4
mean$\{\parallel {\overrightarrow{r}}_{e,i}-{\overrightarrow{r}}_{e,0}\parallel \}$	60.45 m	59.65 m	61.64 m	62.03 m
mean$\{\parallel {\overrightarrow{v}}_{e,i}-{\overrightarrow{v}}_{e,0}\parallel \}$	1.83 m/s	1.84 m/s	1.36 m/s	1.01 m/s

. As can be seen, the relative position and velocity in the steady-state phase are maintained at a reasonable level. The robust property of the proposed control algorithm towards measurement disturbances is validated. For such small ranges of final relative position and velocity difference, high-precision onboard measurement units can then be used to guide the satellites to accomplish tasks such as docking and circumventing.

By observing Figure 8, we find that the satellites are always over the safe height $h_1$, i.e., they would not collide with the Earth’s surface and guarantee $\parallel {\overrightarrow{r}}_i\parallel >h_1,\forall i\in \mathcal{N},t\ge 0$.

The change of masses due to fuel consumption is plotted in Figure 8. The mass losses are related to the initial rendezvous error. Via simple calculation, we obtain the initial rendezvous error for the first to the fourth satellites as $661.13$ km, $1134.40$ km, $569.04$ km, and $1529.90$ km, respectively. These initial errors explain the fuel consumption difference depicted in the figure. Additionally, the masses of satellites for $t\ge 3000$ s decrease slightly to overcome the measurement disturbances.

5. Conclusions

This paper presents the design and analysis of an adaptive rendezvous control strategy for multiple orbital satellites. Nonlinear tools, including a discontinuous approach, the reduced-order method, and saturation-like adaptive control design, are applied to derive a novel control law that steers multiple uncertain orbital satellites to achieve rendezvous at a certain reference orbital vehicle. It is proven that the rendezvous errors are convergent to zero asymptotically. In the future, the authors will consider more problems for the pulsar-based rendezvous control, such as intermitted observation of pulsar sources, statistical measurements, and fuel optimizations.