Autonomous Optimal Absolute Orbit Keeping through Formation Flying Techniques

Abstract

In this paper, the problem of autonomous optimal absolute orbit keeping for a satellite mission in Low Earth Orbit using electric propulsion is considered. The main peculiarity of the approach is to support small satellite missions in which the platform is equipped with a single thruster nozzle that provides acceleration on a single direction at a time. This constraint implies that an attitude maneuver is necessary before or during each thrusting arc to direct the nozzle into the desired direction. In this context, an attitude guidance algorithm specific for the orbit maneuver has been developed. A Model Predictive Control scheme is proposed, where the attitude kinematics are coupled with the orbital dynamics in order to obtain the optimal guidance profiles in terms of satellite state, reference attitude, and thrust magnitude. The proposed control scheme is developed exploiting formation flying techniques where the reference orbit is that of a virtual spacecraft that the main satellite is required to rendezvous with. In addition to the controller design, the closed-loop configuration is presented supported by numerical simulations. The efficacy of the proposed autonomous orbit-keeping approach is shown in several application scenarios.

1. Introduction

In recent years, the adoption of electric propulsion has become more popular among small satellite designers [1]. This trend is driven by the fact that electric propulsion systems are generally more sustainable and more efficient in terms of propellant usage in comparison to chemical propulsion ones. Moreover, electric thrusters could support large maneuvers while being smaller in size and lighter in weight than their chemical counterparts. This, in turn, can allow the allocation of more space/weight for the payload, and in many cases permits successful integration within some launchers that would be infeasible for the satellite if it were equipped with chemical thrusters.

One of the satellites designed to accommodate an electric propulsion system is Triton-X-Heavy [2] (referred to as Triton-X in the discussions that follow), which is a high-performance multi-mission microsatellite platform developed by LuxSpace for Low Earth Orbit (LEO) operations. To permit orbit correction and orbit-transfer maneuvers, Triton-X is equipped with one throttleable electric thruster aligned with its body frame’s z-axis.

Thanks to the simple design of the propulsion system, which requires minimal size and mass budgets, the satellite could be designed with a high level of flexibility, which permits the accommodation of a wide spectrum of payloads. However, having only one thruster onboard comes with its own challenges. For instance, the guidance and control schemes have to be optimized for the unidirectional thrust system, in the sense that attitude and orbit control systems have to be coupled to ensure that the thruster is always aligned with the desired firing direction. Moreover, this need becomes more relevant for orbit corrections that have to be performed on a limited time frame and/or when a very low thrust level is available.

To prepare Triton-X for its prospective missions, the AuFoSat project has been initiated to provide a ready-to-use Guidance, Navigation, and Control (GNC) toolbox specifically tailored for the specifications of Triton-X, which could be seamlessly used to also support formation flying (FF) applications. Previous AuFoSat research concentrated on orbit design [3], and on the relative navigation problem using the onboard sensors of Triton-X [4,5]. This paper, as part of AuFoSat, is focusing on optimal absolute orbit keeping for LEO satellites equipped with a single electric thruster. The problem is approached using formation flying techniques, meaning that Triton-X is assumed to be flying in formation with a virtual chief. Namely, when an orbit transfer maneuver is defined, the chief’s orbit is set as the target one, and the satellite is ultimately required to rendezvous with the virtual chief. On the one hand, opting for the exploitation of FF techniques enables the developed guidance and control schemes to be directly used in multi-satellite scenarios. On the other hand, the absolute orbit control problem can be efficiently tackled by linearizing the orbital dynamics around the reference desired trajectory. In this context, the relative orbit between the satellite and the virtual target can be formulated using the Relative Orbital Elements (ROEs), which render the linearized relative orbital dynamics with high levels of accuracy for neighboring orbits.

Many contributions on autonomous absolute orbit keeping for LEO satellites have been introduced over the last few decades [6,7,8,9]. Among the many, formation flying techniques were used for absolute orbit control in [10], where the spacecraft had to maintain its ground track close to a reference trajectory computed as the orbit of a virtual chief satellite. A specific set of variables was used to parameterize the relative motion and linear and quadratic optimal regulators were employed. In [11,12], formation reconfiguration with realistic mission constraints was considered. The proposed control schemes on these two articles relied on candidate Lyaponuv functions that comprised artificial potential around the prohibited areas (in the decision variables’ space), which guided the main satellite to follow feasible trajectories without having to solve an optimization problem. The main drawback of this approach is that the trajectories followed by the spacecraft are neither fuel- nor time-optimal. All the aforementioned contributions considered only satellites with 3D thrusting capabilities, and most of them adopted impulsive thrust. One mission that adopted a single impulsive thruster is the Autonomous Vision Approach Navigation and Target Identification (AVANTI) experiment [13], where a specific guidance scheme was adopted to include no-control windows that allowed the re-orientation of the available nozzle [14,15]. Furthermore, in [16], formation flying techniques were used to develop a control scheme that constantly corrected the orbit drifts of the satellite that comprised a single impulsive thruster. The satellite used a passive magnetic attitude control system to stabilize its longitudinal axis, along which the thruster was mounted, along the induction vector of the local geomagnetic field.

The only research, to the authors’ knowledge, that addresses the problem of absolute orbit control for unidirectional low-thrust satellites is [17]. Nonetheless, it does not consider the coupled attitude evolution, making such approach not suitable for prompt orbit corrections. Moreover, in [17], both the navigation and actuator errors were also not considered.

In this paper, a Model Predictive Control (MPC) scheme is proposed that couples the attitude and the relative translational equations, where the relative dynamics are formulated in the ROE space. In this setting, the MPC uses an attitude control surrogate model to emulate the functioning of Triton-X’s attitude control system as well as the ROE-based linearized relative dynamical model to predict the satellite’s future roto-translational state. The output of the MPC is a tuple that contains the optimal acceleration that should be provided by the thruster as well as the optimal desired attitude that the attitude control system needs to steer the satellite to. The overall cost function of the MPC at each receding horizon is a trade-off between delta-V optimality, time optimality, and attitude effort optimality, and is specifically designed for the purpose of small orbit changes (less than 1 km) as typically encountered in absolute orbit control problems around a reference orbit. The original contributions of the proposed approach are as follows:

• Considering the coupled attitude/relative orbit problem in the control loop;
• Proposing an MPC scheme for a satellite with a unidirectional propulsion system where the thrust can also be provided during slew maneuvers.

The following section of this paper introduces the formulation of the relative orbital dynamics in the ROE space as well as the adopted surrogate model for the attitude control system. In Section 3, the MPC problem is formally established in the sense that the employed cost function is introduced and the inherent constraints of the problem in hand are presented. The problem of MPC parameter tuning as well as the module execution logic of the closed control loop are discussed in the same section. The stability of the proposed MPC scheme is then tested by means of a Monte Carlo simulation in Section 4, and the proposed MPC is benchmarked against that introduced in [17].

2. Mathematical Models

In this section, the dynamical models used in developing the control schemes are presented. At first, the reference frames used in the research are identified, then the relative orbital dynamics of the deputy spacecraft with respect to the chief (virtual) satellite are shown, and finally the attitude kinematics are introduced.

The following reference frames are used in developing the dynamic models as well as the control scheme:

1.

Earth-centered inertial frame:

The employed Earth-Centered-Inertial (ECI) reference frame, ${\mathbb{F}}^i$, is defined by True-of-Date coordinate system with its origin at the center of the Earth, its x-axis along the true vernal equinox at the current epoch, its z-axis is aligned with the true rotation axis at the current epoch, while its y-axis completes the right-handed triad. A vector expressed in the ECI frame is signified by the ${\left(·\right)}^i$ superscript.

2.

Satellite-body-fixed frame:

The body-fixed reference frame, ${\mathbb{F}}^b$, is a reference frame attached to the spacecraft in question with the origin at the satellite’s center of mass. One common choice of the directions of the three axes is along the three principal axes of inertia of the satellite. A vector expressed in the body-fixed frame of the satellite will have the superscript ${\left(·\right)}^b$.

3.

Radial–transversal–normal frame:

The radial–transversal–normal (RTN) frame, ${\mathbb{F}}^r$, is centered on the center of mass of the chief satellite, where the x-axis of the RTN frame is defined to be along the position vector of the chief satellite, positive towards the Zenith direction, the z-axis is directed along the chief’s orbital angular momentum vector, and the y-axis completes the right-handed set. In the sequel, a vector expressed in the RTN frame is signified by the ${\left(·\right)}^r$ superscript.

Here, and in the discussions to follow, superscripts refer solely to frame labels.

2.1. Relative Orbital Dynamics

The orbital motion of a satellite under the gravitational influence of a major body (e.g., the Earth) can be parameterized in a planet-centered inertial frame using the following set of orbital elements:

$$\mathsf{\alpha }:=\left[\begin{array}{c}a\\ \theta \\ e_x\\ e_y\\ i\\ \Omega \end{array}\right],$$

(1)

where a is the semi-major axis, $\theta =M+\omega$ is the argument of latitude with M being the mean anomaly and $\omega$ being the argument of periapsis, $\mathit{e}:={\left[\begin{array}{cc}{e}_x& e_y\end{array}\right]}^{⊺}={\left[\begin{array}{cc}ecos\omega & esin\omega \end{array}\right]}^{⊺}$ is the eccentricity vector with e being the orbital eccentricity, i is the orbital inclination, and $\Omega$ is the Right Ascension of the Ascending Node (RAAN). If a Cartesian state vector is used, then the motion of the satellite in the planet-centered inertial frame is parameterized by ${\mathit{x}}^i:={\left[\begin{array}{cc}{{\mathit{r}}^i}^{⊺}& {{\mathit{v}}^i}^{⊺}\end{array}\right]}^{⊺}$, where ${\mathit{r}}^i$ and ${\mathit{v}}^i$ are the absolute position and velocity vectors. These latter are related to the orbital elements through a set of nonlinear equations. Spacecraft position and velocity transform into osculating orbital elements, which in the remainder of this work will be denoted by $\tilde{\mathsf{\alpha }}$. Mean orbital elements, denoted by $\mathsf{\alpha }$, are to be intended as one-orbit averaged elements, where the short- and long-term oscillations generated by the $J_2$ harmonic of the Earth’s gravitational potential are removed. Mean/osculating elements’ conversions are performed through the transformations developed in [18]. To this end, the reader is warned not to confuse i, which refers to the orbital inclination, with vectors expressed in the ECI frame, e.g., ${\mathit{x}}^i$.

The relative motion between a deputy and a chief spacecraft can be described by the dimensionless quasi-nonsingular Relative Orbital Element (ROE) vector, which is a nonlinear transformation of the orbital elements vector introduced in (1),

$$\delta \mathsf{\alpha }:=\left[\begin{array}{c}\delta a\\ \delta \lambda \\ \delta e_x\\ \delta e_y\\ \delta i_x\\ \delta i_y\end{array}\right]=\left[\begin{array}{c}\Delta a/a_c\\ \Delta \theta +\Delta \Omega cos{i}_c\\ \Delta e_x\\ \Delta e_y\\ \Delta i\\ \Delta \Omega sin{i}_c\end{array}\right],$$

(2)

where $\delta a$ is the relative semi-major axis, $\delta \lambda$ is the relative mean longitude, $\delta \mathit{e}:={\left[\begin{array}{cc}\delta e_x& \delta e_y\end{array}\right]}^{⊺}$ is the relative eccentricity vector, and $\delta \mathit{i}:={\left[\begin{array}{cc}\delta i_x& \delta i_y\end{array}\right]}^{⊺}$ is the relative inclination vector. It is important to note that here, and in the sequel, the subscript ${\left(·\right)}_d$ denotes a quantity related to the deputy satellite, while the subscript ${\left(·\right)}_c$ is used for chief-related quantities. Moreover, $\delta \left(·\right)$ signifies a relative quantity between the deputy and the chief that is not necessarily the arithmetic difference between that of the deputy and that of the chief, while $\Delta \left(·\right)$ signifies the arithmetic difference between ${\left(·\right)}_d$ and ${\left(·\right)}_c$, i.e., $\Delta \left(·\right):={\left(·\right)}_d-{\left(·\right)}_c$. As for absolute elements, osculating ROEs are denoted by $\delta \tilde{\mathsf{\alpha }}$, whereas mean ROEs are expressed as $\delta \mathsf{\alpha }$. Both sets are computed from the corresponding osculating/mean orbital elements’ set by means of (2).

One of the advantages of the adopted ROEs is that they provide a direct insight on the shape of the relative orbit, unlike the relative position and velocity vectors, which need to be integrated over time in order to obtain a intuition of the relative orbit. Given the set of ROEs in (2), the shape of the relative orbit can be constructed as an ellipse, in both the transversal–radial and normal–radial planes (see Figure 1).

In this work, the adopted relative dynamics are the same as in [19], where the dynamics of the ROEs are linearized to the first order and the secular effect due to the $J_2$ zonal harmonic is included in the relative motion. The following assumptions apply:

• Neighboring orbits of the deputy and the chief,
• Near circular chief’s orbit.

The solution to the linear system under the influence of a constant input acceleration is given in the following form:

$$\delta \mathsf{\alpha }\left(t_{k+1}\right)=\mathsf{\Phi }\left(t_k,t_{k+1}\right)\delta \mathsf{\alpha }\left(t_k\right)+\mathsf{\Psi }\left(t_k,t_{k+1}\right){\mathit{u}}^r\left(t_k,t_{k+1}\right),$$

(3)

where $\mathsf{\Phi }\left(t_k,t_{k+1}\right)\in {\mathbb{R}}^{6\times 6}$ is the State Transition Matrix (STM) between the two time instants, $t_k$ and $t_{k+1}$, $\mathsf{\Psi }\left(t_k,t_{k+1}\right)\in {\mathbb{R}}^{6\times 3}$ is the convolution matrix between the same two time instants, and ${\mathit{u}}^r\left(t_k,t_{k+1}\right)\in {\mathbb{R}}^3$ is the acceleration vector in ${\mathbb{F}}^r$, which is constant over the period $\left[t_k,t_{k+1}\right)$. The reader is referred to [19] for detailed expressions of $\mathsf{\Phi }\left(t_k,t_{k+1}\right)$ and $\mathsf{\Psi }\left(t_k,t_{k+1}\right)$.

In the sequel, and in order to simplify the representation of equations, the following notations are used: ${\mathsf{\Phi }}_{k|k+1}\equiv \mathsf{\Phi }\left(t_k,t_{k+1}\right)$, ${\mathsf{\Psi }}_{k|k+1}\equiv \mathsf{\Psi }\left(t_k,t_{k+1}\right)$, $\delta {\mathsf{\alpha }}_k\equiv \delta \mathsf{\alpha }\left(t_k\right)$, and ${\mathit{u}}_{k|k+1}^r\equiv {\mathit{u}}^r\left(t_k,t_{k+1}\right)$

2.2. Attitude Kinematics

Given that the adopted satellite is equipped with only one thruster, continuous attitude maneuvers are necessary in order to direct the thruster to the desired direction during time-extended orbit maneuvers, which are likely to occur with the available low-thrust level. In this work, it is assumed that the attitude control system of Triton-X is capable of tracking a desired attitude profile. The specifications of the attitude control system are obtained from the publicly available Triton-X brochure (The brochure can be found through: https://luxspace.lu/resources/ (accessed on 1 August 2023) ). A surrogate model for the Attitude Determination and Control System (ADCS) is presented here, but first the attitude kinematics are introduced.

Let ${\mathfrak{q}}^{xy}\in \mathcal{Q}$, with $\mathcal{Q}:=\{\mathfrak{q}\in {\mathbb{R}}^4:∥\mathfrak{q}∥=1\}$ as the unit quaternion that rotates any vector, $\mathit{v}\in {\mathbb{R}}^3$, from one reference frame ${\mathbb{F}}^y$ to another frame ${\mathbb{F}}^x$, through

$${\mathit{v}}^x={\mathfrak{q}}^{xy}\circ {\mathit{v}}^y\circ {\tilde{\mathfrak{q}}}^{xy},$$

(4)

where ∘ is the quaternion multiplication operator, and ${\tilde{\mathfrak{q}}}^{xy}$ is the quaternion conjugate of ${\mathfrak{q}}^{xy}$. Equivalently, the attitude rotation can be performed using the Direction Cosine Matrix (DCM) $\mathbf{A}\left({\mathfrak{q}}^{xy}\right)\in \mathcal{A}$, with $\mathcal{A}:=\{\mathbf{A}\in {\mathbb{R}}^{3\times 3}:{\mathbf{A}}^{⊺}\mathbf{A}=\mathbf{I}\}$ (see Appendix A in [20]) as follows:

$${\mathit{v}}^x=\mathbf{A}\left({\mathfrak{q}}^{xy}\right){\mathit{v}}^y.$$

(5)

The attitude kinematics of a rigid body can be described as follows [21]:

$${\dot{\mathfrak{q}}}^{xy}=\frac{1}{2}{\mathfrak{q}}^{xy}\circ {\omega }^y,$$

(6)

where ${\omega }^y\in {\mathbb{R}}^3$ is the angular velocity of frame ${\mathbb{F}}^y$ with respect to frame ${\mathbb{F}}^x$, expressed in ${\mathbb{F}}^y$.

Letting $\mathfrak{q}\equiv \mathfrak{q}\left(t\right)={\mathfrak{q}}^{rb}\left(t\right)$ and $\mathsf{\omega }\equiv \mathsf{\omega }\left(t\right)={\mathsf{\omega }}^b\left(t\right)$, and letting ${\mathfrak{q}}_r\equiv {\mathfrak{q}}_r\left(t\right)$ be the desired/reference quaternion at time t, (6) can be written as

$$\dot{\mathfrak{q}}=\frac{1}{2}\mathfrak{q}\circ \mathsf{\omega },$$

(7)

and the error quaternion signal, ${\mathfrak{q}}_e$, is defined, in terms of the current attitude and the reference one, as

$${\mathfrak{q}}_e:=\left[\begin{array}{c}{q}_{e,0}\\ {\mathbf{Q}}_e\end{array}\right]:={\tilde{\mathfrak{q}}}_r\circ \mathfrak{q}.$$

(8)

The ultimate behaviour of an ADCS would be to force the error quaternion to asymptotically approach either ${\left[\begin{array}{cccc}1& 0& 0& 0\end{array}\right]}^{⊺}$ or ${\left[\begin{array}{cccc}-1& 0& 0& 0\end{array}\right]}^{⊺}$, which could be achieved through forcing the angular velocity, $\mathsf{\omega }$, profile to follow

$$\mathsf{\omega }=-K{q}_{e,0}{\mathbf{Q}}_e,$$

(9)

where K is a positive control gain. The rationale behind why the $\mathsf{\omega }$ profile is plausible is given in Appendix A.

An additional operational constraint is imposed on the ADCS that dictates that the angular rate around any axis does not exceed the maximum allowable angular rate, i.e., $∥\mathsf{\omega }∥\le {\omega }_{max}$. To be compliant with the attitude control specifications of Triton-X, the value of K is set to $0.2$, while the value of ${\omega }_{max}$ is set to 2 deg/s. Indeed, ${\omega }_{max}$ could be found directly in the data sheets. The numerical value of K, instead, has been assessed through numerical simulations in order to properly emulate the behavior of Triton-X ADCS, based on the time duration to perform a 60-degree slew.

Substituting from (9) into (7), the adopted surrogate model for the ADCS is obtained:

$$\dot{\mathfrak{q}}=-\frac{1}{2}K{q}_{e,0}\mathfrak{q}\circ {\mathit{q}}_e.$$

(10)

One thing that has to be acknowledged is that the attitude and angular velocity profiles of the actual satellite could differ slightly from the profiles generated by the surrogate model assumed here, depending on the exact quaternion feedback regulator used by the ADCS as well as the exact inertia tensor of the satellite. Nonetheless, changing the value of K could bring the attitude and angular velocity profiles of the surrogate model very close to those of the actual ADCS. This claim could be verified through benchmarking the surrogate model against a full attitude control system, which uses one of the commonly used feedback regulator proposed in [22], and the results are depicted in Figure 2. It is to be noted that the initial and reference quaternions are randomly chosen for the simulation in Figure 2. Moreover, the subscript ${\left(·\right)}_{\mathrm{full}}$ refers to a quantity relating to the full attitude control model, while the subscript ${\left(·\right)}_{\mathrm{surr}}$ refers to surrogate model quantities.

3. Model Predictive Control

3.1. MPC Problem Formulation

The Model Predictive Control problem is, in its essence, a recurring optimization problem that searches for the optimal control input over the receding control horizon, and, in most cases, also solves for the state variables over the prediction horizon. Only the first control input in the control horizon is realized by the actuators. In the context of our problem, attitude dynamics are coupled with the ROE-based relative translational motion, and the state vector is set to

$$\mathit{s}=\left[\begin{array}{c}\delta \mathsf{\alpha }\\ \mathfrak{q}\end{array}\right],$$

(11)

where the dynamics of the state vector are described by (3) and (10). The control input, on the other hand, is set to

$$\mathit{u}=\left[\begin{array}{c}u\\ {\mathfrak{q}}_r\end{array}\right],$$

(12)

where u is the magnitude of the piece-wise constant input acceleration provided by the single thruster, and ${\mathfrak{q}}_r$ is the desired quaternion. The state vector as well as the control vector are collated over the prediction and the control horizons in the two matrices $\mathbf{S}$ and $\mathbf{U}$ respectively, such that

$$\begin{array}{cc}\hfill \mathbf{S}& \hfill =\left[\begin{array}{cccc}{\mathit{s}}_0& {\mathit{s}}_1& \cdots & {\mathit{s}}_{n_p}\end{array}\right],\hfill \\ \hfill \mathbf{U}& \hfill =\left[\begin{array}{cccc}{\mathit{u}}_{0|1}& {\mathit{u}}_{1|2}& \cdots & {\mathit{u}}_{n_u-1|n_u}\end{array}\right],\hfill \end{array}$$

(13)

with $n_p$ being the length of the prediction horizon, and $n_u$ being the length of the control horizon, where $n_p\ge n_u$. It is to be noted that the notations ${\mathit{s}}_k:=\mathit{s}\left(t_k\right)$, where $k\in \left\{0,1,\dots ,n_p\right\}$, and ${\mathit{u}}_j:=\mathit{u}\left(t_j\right)$, where $j\in \left\{0,1,\dots ,n_u-1\right\}$, are used, and the time difference between any two consecutive instances is constant and is equal to the user-defined sampling time, i.e., $t_{k+1}-t_k=T_s$. Moreover, after the control horizon is over, the input acceleration is set to zero (i.e., $u_k=0,\forall k\in \left\{n_u,n_u+1,\dots ,n_p\right\}$) and the desired attitude is left to be decided by the mission-specific mode management logic of the ADCS. For the sake of simulating the attitude evolution, the reference attitude at an arbitrary time step beyond the control horizon is set to the quaternion from the previous step. Formally, ${\mathfrak{q}}_r\left(t_k\right)={\mathfrak{q}}_r\left(t_{k-1}\right)\forall k\in \left\{n_u,n_u+1,\dots ,n_p\right\}$.

It is also important to note that the adopted satellite is assumed, without loss of generality of the approach, to have its thruster directed along the z-direction of its body frame, and hence the input acceleration vectors in ${\mathbb{F}}^b$ and in ${\mathbb{F}}^r$ can be related as

$$\begin{array}{cc}\hfill {\mathit{u}}^r& \hfill =\mathbf{A}\left({\mathfrak{q}}^{rb}\right){\mathit{u}}^b,\hfill \\ \hfill {\mathit{u}}^b& \hfill ={\left[\begin{array}{ccc}0& 0& u\end{array}\right]}^{⊺}.\hfill \end{array}$$

(14)

After having introduced the necessary notations, the optimization problem could then be formulated as follows:

Problem 1. 

$$\begin{array}{cc}\hfill & \underset{\mathbf{S},\mathbf{U}}{min}{J}_{\delta \mathsf{\alpha }}+J_u+J_{\delta \mathfrak{q}}\hfill \\ \hfill & \mathit{Subject}\mathit{to}\hfill \end{array}$$

$$\begin{array}{cc}\hfill & \hfill {\mathit{s}}_0=\left[\begin{array}{c}\delta {\mathsf{\alpha }}_0\\ {\mathfrak{q}}_0\end{array}\right],\hfill \end{array}$$

(15)

$$\begin{array}{cc}& \delta {\mathsf{\alpha }}_{k+1}=\left\{\begin{array}{l}{\mathsf{\Phi }}_{k|k+1}\delta {\mathsf{\alpha }}_k+{\mathsf{\Psi }}_{k|k+1}{\overline{\mathit{u}}}_{k|k+1}^r,k\in \left\{0,1,\dots ,n_u-1\right\}\\ {\mathsf{\Phi }}_{k|k+1}\delta {\mathsf{\alpha }}_k,k\in \left\{n_u,n_u+1,\dots ,n_p-1\right\}\end{array}\right.\end{array}$$

(16)

$$\begin{array}{cc}\hfill & \hfill {\mathfrak{q}}_{k+1}=\mathfrak{f}\left(t_k,t_{k+1},{\mathfrak{q}}_k,{\mathfrak{q}}_r\left(t_k\right)\right),\forall k\in \left\{0,1,\dots ,n_p-1\right\}\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill & \hfill {∥{\mathfrak{q}}_k∥}^2=1,\forall k\in \left\{1,\dots ,n_p-1\right\},\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill & \hfill {∥{\mathsf{\omega }}_k∥}^2\le {\omega }_{max}^2,\forall k\in \left\{1,\dots ,n_p-1\right\},\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill & \hfill 0\le u_k\le u_{max},\forall k\in \left\{0,1,\dots ,n_u-1\right\},\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill & \hfill {∥{\mathfrak{q}}_r\left(t_k\right)∥}^2=1,\forall k\in \left\{0,1,\dots ,n_u-1\right\},\hfill \end{array}$$

(21)

where $t_0$, $\delta {\mathsf{\alpha }}_0$, and ${\mathfrak{q}}_0$, are the time, dimensionless ROE vector, and quaternion vector, respectively, at the beginning of each horizon. Note that these are receding quantities that are specific to each prediction horizon. Accordingly, they coincide with the quantities at the beginning of the simulation only once, at the time of the first prediction. Moreover, $u_{max}$ is the maximum applicable acceleration that can be provided by the thruster, $\mathit{\omega }$ is the angular velocity vector of the satellite, and ${\omega }_{max}$ is the maximum allowable angular rate around any axis.

Indeed, a constant value for the input acceleration vector, ${\mathit{u}}_{k|k+1}^r$, has to be fed to the ROE dynamics constraint (16). However, the attitude dynamics are much faster than those of the ROE, and since the vector ${\mathit{u}}_{k|k+1}^r$ is attitude-dependent, its value can change significantly from $t_k$ to $t_{k+1}$, even when ${\mathit{u}}_{k|k+1}^b$ is kept constant. Therefore, a prediction of ${\mathit{u}}_{k|k+1}^r$, ${\overline{\mathit{u}}}_{k|k+1}^r$ had to be fixed and fed to the dynamics constraint (16). The value of ${\overline{\mathit{u}}}_{k|k+1}^r$ is set to that of ${\mathit{u}}_{k|k+1}^r$ at time $\left(t_k+t_{k+1}\right)/2$, with its formal definition being as follows:

$$\begin{array}{cc}\hfill {\overline{\mathit{u}}}_{k|k+1}^r:=& \hfill \mathbf{A}\left(\mathfrak{q}\left(\frac{t_k+t_{k+1}}{2}\right)\right){\mathit{u}}_{k|k+1}^b,\hfill \\ \hfill {\mathit{u}}_{k|k+1}^b:=& \hfill {\left[\begin{array}{ccc}0& 0& u_k\end{array}\right]}^{⊺}.\hfill \end{array}$$

(22)

The attitude kinematics (10) are discretized through a symbolic Runge–Kutta fourth-order scheme and the discrete kinematics are imposed as constraint (17) in Problem 1. Moreover, the discrete angular velocity vector in constraint (19) in Problem 1, ${\omega }_k$, could be obtained by applying Formula (9) using the discrete unit quaternion signal ${\mathfrak{q}}_k$. It is also worth noting that the squared norms (${∥{\mathfrak{q}}_k∥}^2$, ${∥{\omega }_k∥}^2$, and ${∥{\mathfrak{q}}_r\left(t_k\right)∥}^2$) are constrained instead of the norms themselves. This formulation is adopted to facilitate the job of the optimizer since it allows the Jacobian of the constraints to be defined at all values of $\mathbf{S}$ and $\mathbf{U}$.

The cost function components of Problem 1, $J_{\delta \mathsf{\alpha }}$, $J_u$, and $J_{\delta \mathfrak{q}}$, are defined as

$$\begin{array}{cc}\hfill J_{\delta \mathsf{\alpha }}& ={\left(\delta {\mathsf{\alpha }}_{n_p}-\delta {\mathsf{\alpha }}_r\right)}^{⊺}\mathbf{P}\left(\delta {\mathsf{\alpha }}_{n_p}-\delta {\mathsf{\alpha }}_r\right)\hfill \\ \hfill & \phantom{=}+\sum _{k=0}^{n_p-1}{\left(\delta {\mathsf{\alpha }}_k-\delta {\mathsf{\alpha }}_r\right)}^{⊺}\mathbf{Q}\left(\delta {\mathsf{\alpha }}_k-\delta {\mathsf{\alpha }}_r\right),\hfill \\ \hfill J_u& =R_u\sum _{k=0}^{n_u-1}{u}_k^2,\hfill \\ \hfill J_{\delta \mathfrak{q}}& =R_{\delta \mathfrak{q}}\sum _{k=0}^{n_u-1}{\delta {\mathfrak{q}}_k}^{⊺}\delta {\mathfrak{q}}_k,\delta {\mathfrak{q}}_k:={\mathfrak{q}}_r\left(t_k\right)-\mathfrak{q}\left(t_k\right),\hfill \end{array}$$

(23)

where $\delta {\mathsf{\alpha }}_r$ is the reference ROE vector, while $\mathbf{P}\in {\mathbb{R}}^{6\times 6}$, $\mathbf{Q}\in {\mathbb{R}}^{6\times 6}$, $R_u\in {\mathbb{R}}^1$, and $R_{\delta \mathfrak{q}}\in {\mathbb{R}}^1$ are positive-definite MPC gains.

It is important to note that the ultimate goal of autonomous orbit keeping is to rendezvous with a virtual target then track its absolute orbit, i.e., $\delta {\mathsf{\alpha }}_r=\mathbf{0}$. Moreover, the matrices $\mathbf{P}$ and $\mathbf{Q}$ are chosen to be diagonal matrices, which renders $J_{\delta \mathsf{\alpha }}$ a summation over the squared weighted 2-norms of the error vector, $\delta {\mathsf{\alpha }}_k-\delta {\mathsf{\alpha }}_r$. This ultimately means that $J_{\delta \mathsf{\alpha }}$ is a measure of distance between $\delta {\mathsf{\alpha }}_k$ and $\delta {\mathsf{\alpha }}_r$, which is zero in our case, in a scaled-ROE space. Indeed, reasoning $J_{\delta \mathsf{\alpha }}$ in terms of the magnitude of the error vector in a scaled-ROE space does not only allow us to eventually drive $\delta \mathsf{\alpha }$ to $\delta {\mathsf{\alpha }}_r$, and hence for the satellite to rendezvous and then track the orbit of the virtual target, but also allows us to indirectly minimize the total delta-V cost, since the distance in the ROE space, or in a scaled-ROE space for that matter (e.g., the dimensional ROE space, which scales the ROE vector by the mean semi-major axis of the chief), can be used as a measure for the total delta-V cost [14,15]. A closer look at $J_{\delta \mathsf{\alpha }}$ reveals that while minimizing it indeed drives the error signal to zero, it is, in fact, an implicit trade-off between fuel and time optimality, depending on the ratio between the traces of $\mathbf{P}$ and $\mathbf{Q}$, respectively. The greater the value of this ratio, the more inclined towards fuel optimality the cost function becomes, while the lower the ratio, the more the scheme leans towards time optimality.

Adding $J_u$ to the compound cost function directly minimizes the total required thrust from the onboard single thruster. The adopted convex form of $J_u$ is the fuel-optimal form of a cost function, which also coincides with the delta-V-optimal form for single-thruster satellites [23].

Finally, the purpose of adding the last component of the cost function, $J_{\delta \mathfrak{q}}$, is to softly constraint the desired attitude to be as close as possible to the actual attitude in order to avoid unnecessary large slews, which in turn minimizes the attitude control effort. This soft constraint allows large slew maneuvers, depending on the choice of the MPC gains, only when it is meaningful from the ROE error (i.e., $J_{\delta \mathsf{\alpha }}$) or fuel cost (i.e., $J_u$) points of view.

3.2. Parameter Tuning

In the proposed formulation of MPC, there are many parameters and gains to be chosen, which can affect the overall performance and robustness of the MPC. In many industrial applications, such parameters are chosen engineering experience. In this application, the initial guess for the prediction horizon has been formulated based on the literature review. Since the $\Delta V$-optimal locations for ROE changes through impulsive delta-V increments are known to occur, at most, every half-orbit [15], and since the input acceleration of the low-thrust propulsion system should be distributed around these $\Delta V$-optimal locations (if $\Delta V$ optimality is sought), the prediction horizon, $n_p{T}_s$, is chosen to be at least half the orbital period. In this setting, the MPC is able to foresee all the $\Delta V$-optimal locations throughout a certain orbit, which hence leads to optimal allocation of the available thrust. Overly increasing the prediction horizon is expected to only overwhelm the onboard processor, while not having much effect on the optimality of the solution.

The choice of the control horizon, $n_u{T}_s$, on the other hand, is decided through trial and error, which is an iterative process as all the other parameters have to be fixed before the tuning process starts. Nonetheless, regarding the choice of the $n_u$ value, it has to be a small positive integer less than or equal to $n_p$. Common guidelines for choosing the prediction and control horizons can be found in [24].

Tuning of the sampling time, $T_s$, has also been performed through trial and error. In the context of Problem 1, although the dynamics of the ROE are slow and the sampling time could be chosen as a large value in order to not solve the optimization problem so frequently, the following important consideration has to be taken into account to carry out the tuning process. Since the attitude dynamics are much faster than that of the ROE, one cannot ignore that a fixed value for the input acceleration vector in the RTN frame, ${\overline{\mathit{u}}}_{k|k+1}^r$, is fed to the linearized ROE dynamics constraint (16), see (22), which renders the model used in the MPC not accurate unless a small value is chosen for the sampling time. Indeed, including $J_{\delta \mathfrak{q}}$ in the compound cost function helps in slowing down the attitude dynamics. Still, the need to choose a scant value for the sampling time is meaningful. A compromise has to be found between choosing a small $T_s$ and having a more accurate model for the MPC on one hand, and choosing a large $T_s$ that allows the optimization problem to be solved less frequently on the other.

The cost function gains (i.e., $\mathbf{P}$, $\mathbf{Q}$, $R_u$, and $R_{\delta \mathfrak{q}}$ in this paper) are more challenging to tune, since the cost function of the MPC is problem-dependent and no general guidelines exist for its weighting. Thus, more focus is put on analyzing the effect of choosing different cost function weights and on actually selecting an optimal set of these weights.

In order to reduce the number of tunable parameters, the components of the cost function, i.e., $J_{\delta \mathsf{\alpha }}$, $J_u$, and $J_{\delta \mathfrak{q}}$, are related to the main weighting matrix, $\mathbf{Q}$, such that

$$\begin{array}{cc}\hfill \mathbf{P}& \hfill =f_{\mathbf{P}}\mathbf{Q},\hfill \\ \hfill {\overline{J}}_u& \hfill =f_u{\overline{J}}_{\delta \mathsf{\alpha }},\hfill \\ \hfill {\overline{J}}_{\delta \mathfrak{q}}& \hfill =f_{\delta \mathfrak{q}}{\overline{J}}_{\delta \mathsf{\alpha }},\hfill \end{array}$$

(24)

where

$$\begin{array}{ll}{\overline{J}}_{\delta \mathsf{\alpha }}:=& \left(n_p+f_{\mathbf{P}}\right){\overline{\delta \mathsf{\alpha }}}^{⊺}\mathbf{Q}\overline{\delta \mathsf{\alpha }},\\ {\overline{J}}_u:=& n_u{R}_u{\overline{u}}^2,\\ {\overline{J}}_{\delta \mathfrak{q}}:=& N_N{R}_{\delta \mathfrak{q}}{\overline{\delta \mathfrak{q}}}^{⊺}\overline{\delta \mathfrak{q}},\end{array}$$

(25)

with $\overline{\delta \mathsf{\alpha }}$ being the expected value of $\delta \mathsf{\alpha }$ over the prediction horizon, which is approximated by the ROE vector at the beginning of each horizon, and $\overline{u}\approx 0.5u_{max}$ and $\overline{\delta \mathfrak{q}}\approx {\left[\begin{array}{cccc}0.1& 0.1& 0.1& 0.1\end{array}\right]}^{⊺}$ being predictions for the values of u and $\delta \mathfrak{q}$ over the control horizon. The approximated expected values, $\overline{\delta \mathsf{\alpha }}$, $\overline{u}$, and $\overline{\delta \mathfrak{q}}$, are illustrated graphically over one prediction horizon in Figure 3. For imaging purposes, each of them is taken to be one-dimensional. Indeed, the values assumed for $\overline{u}$ and $\overline{\delta \mathfrak{q}}$ are merely the authors’ predictions for the expected values of these quantities over the whole simulation; however, choosing other values for these variables would not affect the final MPC gains, as will be elaborated upon later in the text.

Substituting from (25) into (24), the MPC gains can be related to $\mathbf{Q}$, $\overline{\delta \mathsf{\alpha }}$, $\overline{u}$, and $\overline{\delta \mathfrak{q}}$ as

$$\begin{array}{cc}\hfill R_u& \hfill =\frac{f_u\left(n_p+f_{\mathbf{P}}\right)}{n_u{\overline{u}}^2}{\overline{\delta \mathsf{\alpha }}}^{⊺}\mathbf{Q}\overline{\delta \mathsf{\alpha }},\hfill \\ \hfill R_{\delta \mathfrak{q}}& \hfill =\frac{f_{\delta \mathfrak{q}}\left(n_p+f_{\mathbf{P}}\right)}{n_u{\overline{\delta \mathfrak{q}}}^{⊺}\overline{\delta \mathfrak{q}}}{\overline{\delta \mathsf{\alpha }}}^{⊺}\mathbf{Q}\overline{\delta \mathsf{\alpha }},\hfill \end{array}$$

(26)

where $R_u$ and $R_{\delta \mathfrak{q}}$ change for every prediction horizon.

It is conceivable that when $\delta \mathsf{\alpha }:=∥\delta \mathsf{\alpha }∥$ is very small (i.e., the rendezvous of the satellite with its virtual target has already taken place), the priority should be given to $J_{\delta \mathfrak{q}}$ in order for the ADCS not to hyper-react to the small errors in $\delta \mathsf{\alpha }$. Indeed, the phenomenon of hyper-reaction of the ADCS when $\delta \mathsf{\alpha }$ becomes too small has been verified via numerical simulation for a system using the definition of $R_{\delta \mathfrak{q}}$ in (26). It is for this reason that $R_{\delta \mathfrak{q}}$ has to be redefined in order to steer the priority to $J_{\delta \mathfrak{q}}$ when $\delta \mathsf{\alpha }$ is approaching zero, such that

$$R_{\delta \mathfrak{q}}=\mathrm{Sat}\left(\frac{f_{\delta \mathfrak{q}}\left(n_p+f_{\mathbf{P}}\right)}{n_u{\overline{\delta \mathfrak{q}}}^{⊺}\overline{\delta \mathfrak{q}}}{\overline{\delta \mathsf{\alpha }}}^{⊺}\mathbf{Q}\overline{\delta \mathsf{\alpha }},{R_{\delta \mathfrak{q}}}_{min},\infty \right),$$

(27)

where $\mathrm{Sat}\left(x,x_{min},x_{max}\right)$ is the saturation function which is defined as follows:

$$\mathrm{Sat}\left(x,x_{min},x_{max}\right):=\left\{\begin{array}{ccc}{x}_{min},& & x\le x_{min}\\ x,& & x_{min}\le x\le x_{max}\\ x_{max},& & x\ge x_{max}\end{array}\right..$$

(28)

Looking at Equations (24)–(27), it is clear that in order to define all the MPC gains, one has to choose $\mathbf{Q}$, $f_{\mathbf{P}}$, $f_u$, $f_{\delta \mathfrak{q}}$, and ${R_{\delta \mathfrak{q}}}_{min}$.

The choice of $\mathbf{Q}$ is rather simple since all the ROEs can be weighted equally. However, in order to enhance the stability of the final relative orbit (i.e., to keep the obtained orbit for as long as possible without it being affected much by the perturbations), an emphasis is put on minimizing $\delta a$ when the ROEs are close to the target ones, and $\mathbf{Q}$ is defined as

$$\mathbf{Q}=\left\{\begin{array}{ccc}Q{\mathbf{I}}_6,& & a\delta \alpha \le 1m\\ \\ Q\left[\begin{array}{cc}100& 0\\ 0& {\mathbf{I}}_5\end{array}\right],& & a_c\delta \alpha >1m\end{array}\right.,$$

(29)

where Q is a tuning parameter that controls the order of magnitude of the cost function, and ${\mathbf{I}}_n$ is the identity matrix with dimensions $n\times n$.

Fixing the value of Q to ${10}^5$, of $f_{\mathbf{P}}$ to 10, and of ${R_{\delta \mathfrak{q}}}_{min}$ to ${10}^{-5}$, a brute force sensitivity analysis over $f_u$ and $f_{\delta \mathfrak{q}}$ is carried out to choose adroit values for $f_u$ and $f_{\delta \mathfrak{q}}$. Thanks to the High-Performance Computer of the University of Luxembourg [25], 726 simulations could be run simultaneously where, in addition to interchanging the values of $f_u$ and $f_{\delta \mathfrak{q}}$, the seed of the Random Number Generator (RNG), which controls the initial state at the beginning of the simulation, could also be taken into account. In these simulations, and in the sequel, warm-starting is employed at the beginning of each prediction horizon, meaning that the optimized state and control profiles from the previous prediction horizon are utilized to construct the initial guess for the current one.

Before starting the simulations, it was noted that the admissible values for $f_u$ and $f_{\delta \mathfrak{q}}$ belong to the period $\left[0,1\right]$, since the most important component of the cost function is $J_{\delta \mathsf{\alpha }}$ and since the two other components are normalized with respect to $J_{\delta \mathsf{\alpha }}$ (see (24)).

One simulation would have its RNG seed, $f_u$, and $f_{\delta \mathfrak{q}}$ as one possible combination from the following sets:

$$\begin{array}{cc}\hfill & \mathrm{seed}\in \left\{0,1,2,3,4,5\right\},\hfill \\ \hfill & f_u,f_{\delta \mathfrak{q}}\in \left\{0,0.02,0.04,0.06,0.08,0.1,0.3,0.6,0.7,0.9,1\right\},\hfill \end{array}$$

(30)

To this end, a fitness function needs to be introduced in order to asses how good the output of the simulation is. Such fitness function in our context needed to address four criteria (performance metrics):

1.

Driving the ROE vector to zero, which is the main goal of the MPC scheme to achieve orbit keeping. This is assessed through observing the norm of the finale of the ROE time series, denoted as $\delta {\alpha }_{\mathrm{fin}}$.

2.

Enhancing the stability of the relative orbit finale which, in other words, is minimizing the Root Mean Square (RMS) of the last portion of the relative semi-major axis profile over time, denoted as $\delta a_{\mathrm{fin}}$.

3.

Minimizing the total delta-V, $\Delta V_{\mathrm{tot}}$.

4.

Reducing the attitude effort, which is quantified through the mean angular rate of the satellite, ${∥\mathsf{\omega }∥}_{\mathrm{mean}}$.

It is necessary to state that the term “finale” in this paper refers to the last $10\%$ of the simulation time span. Furthermore, although this research is concerned with low-thrust propellers, the total delta-V is considered instead of the total thrust in order to enable comparisons between the proposed scheme and those from the literature. The total delta-V in our context can be calculated using the following formula:

$$\Delta V_{\mathrm{tot}}={\int }_{t_{\mathrm{init}}}^{t_{fin}}udt,$$

(31)

where $t_{\mathrm{init}}$ and $t_{fin}$ are the initial and final times of the simulation, respectively. Given that the input thrust/acceleration from the single thruster is piece-wise-constant, the total delta-V can be rewritten as follows:

$$\Delta V_{\mathrm{tot}}=\sum _{i=1}^{n_h}{\left(\sum _{k=0}^{n_u-1}{u}_k{T}_s\right)}_j,$$

(32)

with $n_h$ being the number of receding horizons being processed within a simulation, j the horizon index, and $T_s$ being the fixed sampling time.

The adopted overall fitness function could then be written as

$$\begin{array}{c}\varphi =K_{\delta \mathsf{\alpha }}·\underset{{\varphi }_{\delta \mathsf{\alpha }}}{\underbrace{\mathrm{fitness}\left(\delta {\mathsf{\alpha }}_{\mathrm{fin}}\right)}}+K_{\delta a}·\underset{{\varphi }_{\delta a}}{\underbrace{\mathrm{fitness}\left(\delta a_{\mathrm{fin}}\right)}}\hfill \\ \hfill +K_{\Delta V}·\underset{{\varphi }_{\Delta V}}{\underbrace{\mathrm{fitness}\left(\Delta V_{\mathrm{tot}}\right)}}+K_{∥\mathsf{\omega }∥}·\underset{{\varphi }_{∥\mathsf{\omega }∥}}{\underbrace{\mathrm{fitness}\left({∥\mathsf{\omega }∥}_{\mathrm{mean}}\right)}},\end{array}$$

(33)

where $K_{\delta a}$, $K_{\delta \mathsf{\alpha }}$, $K_{\Delta V}$, and $K_{∥\mathsf{\omega }∥}$ are weights that determine the importance of each of the four criterion, and the function $\mathrm{fitness}\left(·\right)$ is defined as

$$\mathrm{fitness}\left(x\right)=\frac{1/x-min\left(1/x\right)}{max\left(1/x\right)-min\left(1/x\right)},$$

(34)

with $min\left(·\right)$ and $max\left(·\right)$ being the minimum and maximum functions over all the simulations. It is clear that this form of the fitness function only shoots out values between 0 and 1.

Running the simulations for the 726 combinations in (30) for a fixed initial chief orbit, defined in

Table 1. Chief’s initial orbit for all the 726 simulations used in the sensitivity analysis.

${\tilde{\mathit{a}}}_{\mathit{c}}$ [km]	${\tilde{\mathit{\theta }}}_{\mathit{c}}$ [${}^{\circ }$]	${\left.{\tilde{\mathit{e}}}_{\mathit{x}}\right|}_{\mathit{c}}$	${\left.{\tilde{\mathit{e}}}_{\mathit{y}}\right|}_{\mathit{c}}$	${\tilde{\mathit{i}}}_{\mathit{c}}$ [${}^{\circ }$]	${\tilde{\mathbf{\Omega }}}_{\mathit{c}}$ [${}^{\circ }$]
7121	0	${10}^{-5}$	0	45	0

, applying the definition of the fitness function in (34) to all the 726 simulations and to the four performance metrics, and fixing the values of the weights to $K_{\delta \alpha }=5$, $K_{\delta a}=1$, $K_{\Delta V}=1$, and $K_{∥\mathsf{\omega }∥}=1$, which reflect the high importance of minimizing $\delta {\alpha }_{\mathrm{fin}}$ in comparison to the other metrics, the heat maps that summarize all the simulations could be obtained and are presented in Figure 4.

The heat maps in Figure 4 could only be obtained after averaging the fitness over each RNG seed in (30), and after adopting a scattered data interpolant in order to smooth the heat maps.

It is clear from the figures that the smaller $f_u$ is, the better ${\varphi }_{\delta \alpha }$ and ${\varphi }_{\delta a}$ become, while this relation is conceivably inverted for ${\varphi }_{\Delta V}$. Furthermore, ${\varphi }_{\Delta V}$ and ${\varphi }_{∥\mathsf{\omega }∥}$ could be noticed to be positively correlated, and are both negatively correlated with ${\varphi }_{\delta \alpha }$ and ${\varphi }_{\delta a}$. Further investigations on the simulations where ${\varphi }_{\Delta V}$ and ${\varphi }_{∥\mathsf{\omega }∥}$ are both large revealed that minimal input acceleration is provided at the current attitude of the satellite, which is almost stagnant, while the effect of this acceleration on the orbit correction is, as well, minimal. Nonetheless, the ROE vector is approaching its set point, although very slowly.

The overall fitness function, $\varphi$, is depicted in Figure 5 and the fittest point, i.e., the point with maximum overall fitness, ($f_u=0$ and $f_{\delta \mathfrak{q}}=0.02$) is marked with the gray circle on the heat map. It is clear from the fittest $f_u-f_{\delta \mathfrak{q}}$ combination that, for the initial chief orbit in Table 1, it was not necessary to include the direct delta-V penalty, $J_u$, in the cost function, since $J_{\delta \mathsf{\alpha }}$ was already indirectly accounting for delta-V cost as previously mentioned. One has to acknowledge that the fittest combination of $f_u$ and $f_{\delta \mathfrak{q}}$ in Figure 5 is not necessarily universal for any arbitrary initial chief’s orbit, and the selection of the optimal $f_u-f_{\delta \mathfrak{q}}$ combination needs to be performed again for each initial orbit of the virtual chief. This does not represent a limitation of the approach, since the reference orbit of the chief is fixed once as soon as the scientific mission is designed.

The employed approach served the need to identify an adroit set of combinations of $f_u$ and $f_{\delta \mathfrak{q}}$. Future investigations may focus on the setup of a heuristic approach, e.g., genetic algorithms, to search for the combination that globally maximizes the overall fitness function. In this context, note that the definition of the fitness functions and of the related parameters impact the global optimum.

Finally, the values of the optimized $f_u$ and $f_{\delta \mathfrak{q}}$ are heavily dependant on the authors’ predictions of $\overline{u}$ and $\overline{\delta \mathfrak{q}}$, and another choice of these quantities would have definitely led to different optimal $f_u$ and $f_{\delta \mathfrak{q}}$; however, the values of $R_u$ and $R_{\delta \mathfrak{q}}$ will not change since $f_u$ is directly divided by ${\overline{u}}^2$ to obtain $R_u$ and the same dynamics apply for $f_{\delta \mathfrak{q}}$ with the squared 2-norm of $\overline{\delta \mathfrak{q}}$ to obtain $R_{\delta \mathfrak{q}}$ (see (26)).

3.3. Closing the Control Loop

The MPC described so far is one key module of the whole control loop, which, at practical implementation, requires feedback to operate properly, which is provided by the navigation system, which, in turn, relies on sensors as well as a state estimation filter. The full closed loop, which illustrates the module execution logic onboard the deputy spacecraft, is depicted in Figure 6. There, Osc2Mean is the function that transforms osculating orbital elements to mean ones [18], Body2Inert is the method that rotates any vector from the body frame to the inertial frame, and OE2ROE is the function that transforms the absolute orbital elements of both the chief and the deputy to a ROE vector (see (2)). The breve accent, $\breve{\left(·\right)}$, signifies a quantity which is affected by either one or a combination of the following sources of errors:

• Estimation errors (e.g., ${\breve{\mathsf{\alpha }}}_d$);
• Actuator errors (e.g., $\breve{\mathfrak{q}}$);
• Physical constraints (e.g., ${\breve{\mathit{u}}}^b$).

The proposed MPC is validated by means of numerical simulations, emulating the performance of the navigation module and the effects of the physical limitations of the adopted actuators.

The physical limitations are only present in the saturation block (see Figure 6), which could be implemented in the numerical simulations. In fact, this saturation is taken into account in the MPC implementation, but the saturation block after the MPC is implemented anyway as a safeguard in case the MPC computes an infeasible solution.

The “Sensors” and “Filter” blocks in Figure 6 are replaced by the following surrogate models, which treat the propagated state of the deputy as ground truth:

$$\begin{array}{cc}\hfill & \hfill {\breve{\tilde{\mathit{x}}}}_d^i={\tilde{\mathit{x}}}_d^i+\mathcal{N}\left(\mathbf{0},{\Sigma }_{\tilde{\mathit{x}}}\right),\hfill \\ \hfill & \hfill {\breve{\tilde{\mathsf{\alpha }}}}_d=\mathrm{Cart}2\mathrm{OE}\left({\breve{\tilde{\mathit{x}}}}_d^i\right),\hfill \end{array}$$

(35)

where “Cart2OE” is the method that transforms the Cartesian state to orbital elements, and $\mathcal{N}\left(\mu ,{\sigma }^2\right)$ is a normally distributed random variable with $\mu$ as its expected value and ${\sigma }^2$ as its variance. Hence, ${\Sigma }_{\tilde{\mathit{x}}}$ is the covariance matrix of the random noise affecting the estimation of the Cartesian state of the deputy satellite, which is defined as

$${\Sigma }_{\tilde{\mathit{x}}}=\left[\begin{array}{cccccc}{\sigma }_{\mathit{r}}^2& 0& 0& 0& 0& 0\\ 0& {\sigma }_{\mathit{r}}^2& 0& 0& 0& 0\\ 0& 0& {\sigma }_{\mathit{r}}^2& 0& 0& 0\\ 0& 0& 0& {\sigma }_{\mathit{v}}^2& 0& 0\\ 0& 0& 0& 0& {\sigma }_{\mathit{v}}^2& 0\\ 0& 0& 0& 0& 0& {\sigma }_{\mathit{v}}^2\end{array}\right],$$

(36)

where ${\sigma }_{\mathit{r}}^2$ and ${\sigma }_{\mathit{v}}^2$ are the variances of the one-dimensional position and velocity estimation errors, respectively, which, for the numerical simulations presented in this work, were extracted from the publicly available Triton-X brochure.

It is to be noted that the relative navigation is typically more accurate than the absolute one [4,5]. It is for this reason that a similar model is used to emulate the behaviour of the relative navigation system within the numerical simulations as follows:

$$\begin{array}{cc}\hfill & \hfill \delta \mathsf{\alpha }=\mathrm{OE}2\mathrm{ROE}\left({\mathsf{\alpha }}_c,{\mathsf{\alpha }}_d\right)\hfill \\ \hfill & \hfill \breve{\delta \mathsf{\alpha }}=\delta \mathsf{\alpha }+\mathcal{N}\left(\mathbf{0},{\Sigma }_{\delta \mathsf{\alpha }}\right),\hfill \end{array}$$

(37)

where ${\Sigma }_{\delta \mathsf{\alpha }}$ is the covariance matrix of the random noise affecting the estimation of the ROE vector, which is expanded as

$${\Sigma }_{\delta \mathsf{\alpha }}=\left[\begin{array}{cccccc}{\sigma }_{\delta a}^2& 0& 0& 0& 0& 0\\ 0& {\sigma }_{\delta \lambda }^2& 0& 0& 0& 0\\ 0& 0& {\sigma }_{\delta e_x}^2& 0& 0& 0\\ 0& 0& 0& {\sigma }_{\delta e_y}^2& 0& 0\\ 0& 0& 0& 0& {\sigma }_{\delta i_x}^2& 0\\ 0& 0& 0& 0& 0& {\sigma }_{\delta i_y}^2\end{array}\right],$$

(38)

with the diagonal elements being the variances of the random noise that affect each element of the ROE vector. These variances are taken to all be equal in the sequel, i.e., ${\sigma }_{\delta a}={\sigma }_{\delta \lambda }={\sigma }_{\delta e_x}={\sigma }_{\delta e_y}={\sigma }_{\delta i_x}={\sigma }_{\delta i_y}={\sigma }_{\delta \mathsf{\alpha }}$.

Lastly, the “Pointing error” block is replaced by the following surrogate model:

$$\begin{array}{cc}\hfill & \breve{\mathfrak{q}}=\mathfrak{q}\circ {\mathfrak{q}}_{\mathrm{pe}},\hfill \\ \hfill & {\mathfrak{q}}_{\mathrm{pe}}={\left[\begin{array}{cc}cos\left(\frac{{\zeta }_{{\mathfrak{q}}_{\mathrm{pe}}}}{2}\right)& sin\left(\frac{{\zeta }_{{\mathfrak{q}}_{\mathrm{pe}}}}{2}\right){\widehat{\mathit{q}}}_{\mathrm{pe}}^{⊺}\end{array}\right]}^{⊺},\hfill \end{array}$$

(39)

where ${\zeta }_{{\mathfrak{q}}_{\mathrm{pe}}}$ is the pointing error angle, which was obtained from the Triton-X brochure in the validation simulations, and ${\widehat{\mathit{q}}}_{\mathrm{pe}}$ is a random three-element unit vector.

It is important to note that although the chief, in this context, is a virtual point, it is propagated using the perturbed dynamics for two main reasons:

• In many applications, the tracking of the reference orbit and the absolute control of the reference orbit are handled separately;
• The proposed scheme can be directly applied to the rendezvous with an actual spacecraft.

As a result, the autonomous orbit keeping compensates errors and perturbations acting on the relative dynamics.

4. MPC Validation

In order to test the stability of the proposed MPC scheme, and thanks to the parallelization capabilities of the HPC facilities of the University of Luxembourg, the algorithm was run over 500 different simulations, each with a different initial ROE vector, such that

$$a_c·\delta \mathsf{\alpha }\left(t_{\mathrm{init}}\right)\left[\mathrm{m}\right]\in \left[\begin{array}{c}\left[-100,100\right]\\ \left[-1000,1000\right]\\ \left[-100,100\right]\\ \left[-100,100\right]\\ \left[-100,100\right]\\ \left[-100,100\right]\end{array}\right],$$

(40)

where $t_{\mathrm{init}}$ is the initial time of each simulation.

The simulation and MPC parameters used in the 500-run Monte Carlo simulation are summarized in

Table 2. Monte Carlo simulation and MPC parameters.

Chief’s initial orbit	${\tilde{a}}_c$ [km]	${\tilde{\theta }}_c$ [${}^{\circ }$]	${\left.{\tilde{e}}_x\right|}_c$	${\left.{\tilde{e}}_y\right|}_c$	${\tilde{i}}_c$ [${}^{\circ }$]	${\tilde{\Omega }}_c$ [${}^{\circ }$]
	7121	0	${10}^{-5}$	0	45	0
MPC parameters		$T_s$ [s]	$n_p$ [$T_s$]	$n_u$ [$T_s$]	$u_{max}$ [$\mathrm{m}/{\mathrm{s}}^2$]	${\omega }_{max}$ [${}^{\circ }/\mathrm{s}$]
		50 s	60	15	$3.5\times {10}^{-5}$	2
		Q	$f_{\mathbf{P}}$	$f_u$	$f_{\delta \mathfrak{q}}$	$R_{{\delta \mathfrak{q}}_{min}}$
		${10}^5$	10	$0.0$	$0.02$	${10}^{-5}$
Miscellaneous		K	${\sigma }_{\mathit{r}}$ [m]	${\sigma }_{\mathit{v}}$ [$\mathrm{m}/\mathrm{s}$]	$a_c{\sigma }_{\delta \mathsf{\alpha }}$ [m]	${\zeta }_{{\mathfrak{q}}_{\mathrm{pe}}}$ [arcs]
		$0.2$	10	0.5	1	25

.

It is to be noted that although the sampling time used in predicting the MPC states, $T_s$, is 50 s, the attitude dynamics are propagated each second in the simulations since the attitude changes much faster than the ROEs. Moreover, constraining the maximum allowable acceleration is a consequence of the physical constraint on the maximum thrust that can be provided by the onboard throttleable electric propulsion system. The value of the maximum allowable acceleration, $u_{max}$, was obtained for the following simulations from the information on Triton-X brochure on both, the maximum thrust, which is set equal to 7 mN, and the mass of the satellite, which is assumed to be constant, 200 kg, throughout the maneuver.

The main performance metrics over the 500 Monte Carlo runs are depicted in Figure 7, while the statistics of these metrics are presented in

Table 3. Performance metrics over 500 simulations.

	${\mathit{a}}_{\mathit{c}}\mathit{\delta }{\mathit{\alpha }}_{\mathbf{fin}}$ [m]	${\mathit{a}}_{\mathit{c}}\mathit{\delta }{\mathit{a}}_{\mathbf{fin}}$ [m]	$\mathbf{\Delta }{\mathit{V}}_{\mathbf{tot}}$ [$\mathit{m}/\mathit{s}$]	${∥\mathit{\omega }∥}_{\mathbf{mean}}$ [${}^{\circ }/\mathit{s}$]
Mean	$4.38$	$2.05$	$0.406$	$0.228$
Median	$4.29$	$2.02$	$0.409$	$0.226$
Max	$9.82$	$4.49$	$0.58$	$0.331$

.

In order to gain an insight into how the proposed scheme performs, the results of an arbitrary simulation out of the 500 simulations are presented. The randomly selected initial ROE (according to (40)) vector is given in

Table 4. Parameters of an arbitrary simulation out of the 500 Monte Carlo runs.

Initial ROE	$a_c\delta a$ [m]	$a_c\delta \lambda$ [m]	$a_c\delta e_x$ [m]	$a_c\delta e_y$ [m]	$a_c\delta i_x$ [m]	$a_c\delta i_y$ [m]
	$64.62$	$-947.77$	$-57.84$	$23.68$	$-80.35$	$24.03$
Reference ROE	$a_c\delta a$ [m]	$a_c\delta \lambda$ [m]	$a_c\delta e_x$ [m]	$a_c\delta e_y$ [m]	$a_c\delta i_x$ [m]	$a_c\delta i_y$ [m]
	0	0	0	0	0	0

alongside the reference ROE vector.

The ROE states are seen to approach their set points (zero for rendezvous with the virtual target) in Figure 8. It is of significance to mention that the chief’s mean orbital period, which is the unit of time in Figure 8, is $5974.46$ s.

The trajectory followed by the satellite is depicted in Figure 9 in both the transversal–radial and the normal–radial planes. This trajectory could be obtained by transforming the ROE to the relative Cartesian states in the RTN frame using the transformation in [26]. It can be seen from Figure 9a that the orbit maneuver is taking place gradually thanks to the continuous firing of the thruster, in contrast to the impulsive-thrust maneuvers where the velocity changes at a handful of points on the orbit.

The thrust exerted by the onboard thruster is depicted in Figure 10. This thrust is also projected on the RTN frame and the projection is presented in the same figure. It is interesting to see that the MPC does not follow a trajectory that minimizes the radial burns, which has been a usual constraint imposed on the MPC in earlier studies [17]. Although firing in the radial direction to correct $\delta \lambda$ errors is known to be more delta-V-expensive than firing in the transversal direction to build up $\delta a$ drift, which takes care of correcting the $\delta \lambda$ error by exploiting the natural dynamics, this is not the case for close-proximity maneuvers where the initial $\delta \lambda$ error is small. This can be verified by looking at the solution to the linearized forced Gauss variational equations [19], and it could also be verified via our preliminary experiments in which the projection of thrust on the radial direction was soft-constrained.

In Figure 11a, the error quaternion angle signal suggests that at each optimization step, the attitude starts from a value that does not coincide with the reference attitude. Nevertheless, by the end of each step, the attitude gradually converges towards the reference one, resulting in the error quaternion angle approaching a value close to zero. A very interesting feature of the proposed scheme is that the propulsion system is always turned on even when attitude redirection maneuvers are taking place. The optimizer calculates the reference attitude knowing that the satellite will not necessarily be aligned with it for the whole upcoming control step.

It is worth noting that the angle of any quaternion (e.g., the error quaternion, ${\mathfrak{q}}_e$, or the current quaternion, $\mathfrak{q}$) can be found as follows:

$${\zeta }_{\mathfrak{q}}=arctan\left(\frac{∥\mathit{q}∥}{q_0}\right),$$

(41)

where $∥\mathit{q}∥$ and $q_0$ are the vector and the scalar parts of the quaternion $\mathfrak{q}$, respectively. Lastly, the effect of adding the $J_{\delta \mathfrak{q}}$ term to the cost function of Problem 1 could be seen clearly in Figure 11b, where the quaternion angle rate is depicted. The quaternion angle rate, which can be thought of as a signed version of $∥\mathsf{\omega }∥$, is seen to be much lower than its maximum allowable values, $2^{\circ }/\mathrm{s}$. Our preliminary simulations, which did not contain the term $J_{\delta \mathfrak{q}}$ on their cost function, always had the attitude rate saturating at its maximum value, which is a hyper-reaction that is not desired in a real mission.

In order to test the efficiency and the optimality of the proposed MPC scheme, it is benchmarked against the approach proposed in [17] for a similar problem, where an out-of-plane (OOP) relative orbit correction maneuver had to be performed. The parameters used for this benchmark simulation are summarized in

Table 5. Benchmark simulation parameters.

Chief’s initial orbit	${\tilde{a}}_c$ [km]	${\tilde{\theta }}_c$ [${}^{\circ }$]	${\left.{\tilde{e}}_x\right|}_c$	${\left.{\tilde{e}}_y\right|}_c$	${\tilde{i}}_c{[}^{\circ }]$	${\tilde{\Omega }}_c$ [${}^{\circ }$]
	6828	0	${10}^{-5}$	0	78	0
Initial ROE	$a_c\delta a$ [m]	$a_c\delta \lambda$ [m]	$a_c\delta e_x$ [m]	$a_c\delta e_y$ [m]	$a_c\delta i_x$ [m]	$a_c\delta i_y$ [m]
	0	0	273	0	10	70
Reference ROE	$a_c\delta a$ [m]	$a_c\delta \lambda$ [m]	$a_c\delta e_x$ [m]	$a_c\delta e_y$ [m]	$a_c\delta i_x$ [m]	$a_c\delta i_y$ [m]
	0	0	273	0	400	120

. Any missing parameter in Table 5 is set to its value in Table 2, except for the maximum allowable acceleration, which is set to $u_{max}=3.2\times {10}^{-5}\mathrm{m}/{\mathrm{s}}^2$. Applying this maximum acceleration to the satellite employed in [17], which has a mass of 20 kg, renders the maximum thrust as $0.64$ mN.

The setting of this maneuver is interesting not only because it allows two different MPC schemes to be compared, but also because OOP maneuvers require only normal acceleration. In fact, the $\Delta V$-optimal locations for impulsive normal burns for such maneuvers could be analytically calculated [15]. For the case of electric propulsion, a controller is said to be $\Delta V$-optimal for out-of-plane maneuvers if the thrust in the normal direction is bang-bang-like and is distributed almost evenly around the $\Delta V$-optimal locations. The thrust profile of the proposed MPC is depicted in Figure 12a, which reveals that the proposed MPC behaves as a fuel-efficient scheme would be expected to behave.

It can also be seen from Figure 12a that the thrust components in the radial and the transversal directions are not exactly zero due to in-plane perturbation compensation and also owing to the fact that the thrust is provided during slew maneuvers. Therefore, the radial and transversal thrust components become more visible when the thrust direction is required to flip from the positive to the negative normal direction. In this situation, the satellite has to pass by many transient attitudes for which the radial and the transversal components are not necessarily zero.

Figure 12b depicts a comparison between the normal component of the thrust profile of the proposed MPC scheme and the MPC scheme proposed by [17].

The dimensional ROE profiles of the two MPCs are also shown in Figure 13, and a brief comparison between the performance of both schemes is presented in

Table 6. Comparison between the proposed and the reference MPCs.

	Convergence Time [orbits]	Terminal ${\mathit{a}}_{\mathit{c}}\left(\mathit{\delta }\mathit{\alpha }-\mathit{\delta }{\mathit{\alpha }}_{\mathit{r}}\right)\left[\mathbf{m}\right]$	$\mathbf{\Delta }{\mathit{V}}_{\mathbf{tot}}[\mathbf{m}/\mathbf{s}]$
Proposed MPC	$4.1$	${\left[\begin{array}{cccccc}0.9& -0.4& 1.9& -0.33.& 0.0& -0.5\end{array}\right]}^{⊺}$	$0.664$
Reference MPC	$6.6$	${\left[\begin{array}{cccccc}-3.6& 9.2& -1.4& 2.0& -2.9& 1.6\end{array}\right]}^{⊺}$	$0.501$

, where the convergence time is defined as the time it takes all elements in the error-dimensional ROE vector, $a_c\delta {\mathsf{\alpha }}_{\mathrm{err}}=a_c\left(\delta \mathsf{\alpha }-\delta {\mathsf{\alpha }}_r\right)$, to be less than $5\left[\mathrm{m}\right]$.

It is clear from Table 6 as well as Figure 13 that, for the simulation defined by Table 5, the proposed scheme is much faster, in fact $1.61$ times as fast as that of [17]. Furthermore, the proposed MPC appears to be much more precise than the reference one; however, it consumes $24.04\%$ more delta-V, as seen in Table 6. Indeed, consuming more delta-V than the reference MPC is conceivable since the proposed MPC is not designed to be delta-V-optimal in the first place. Again, the cost function of Problem 1 has two components (after eliminating $J_u$ as discussed in Section 3.2): the first one, $J_{\delta \mathsf{\alpha }}$, deals with gradually approaching the reference ROE vector, and implicitly implies a trade-off between fuel and time optimality (depending on $f_{\mathbf{P}}$), while the second component, $J_{\delta \mathfrak{q}}$, is added to minimize the attitude change throughout the maneuver. It is believed that a smaller total delta-V could be achieved using a different set of MPC gains; however, it would be at the cost of more abrupt attitude changes and/or more maneuvering time. The ROE profile in Figure 13 implies the evolution of the OOP variables, $\delta i_x$ and $\delta i_y$, from their initial values to their set points, while the in-plane variables fluctuate slightly around their initial/reference values as a response to the relative orbital perturbations as well as to the unavoidable radial and transversal accelerations during each attitude flip from the positive to the negative normal direction and vice-versa.

Lastly, Figure 14 relates to the attitude evolution of the system, where Figure 14a shows how the attitude is changing smoothly to recurrently flip the thrust direction from/to the positive to/from the negative normal direction, while the error quaternion angle is starting from a nonzero value only when the thrust direction is required to be flipped, and is seen to always approach zero. Moreover, although the objective is clearly to flip the thrust direction, the error quaternion angle is never ${180}^{\circ }$, which could only mean that the optimizer chooses the reference quaternion to evolve smoothly until it reaches the flipped attitude so that a hyper-reaction of angular velocity can be avoided. The quaternion angle rate is seen in Figure 14b to turn swiftly and almost reach its full potential, i.e., $\pm {\omega }_{max}$, when the thrust direction is set to change from/to the positive to/from the negative normal direction.

5. Conclusions

This paper proposes an MPC scheme to address the problem of optimal absolute orbit keeping for underactuated satellites that use electrical propulsion systems. Currently, several small-satellite platforms are equipped with unidirectional thrusters, requiring constant attitude slews to take place during time-extended orbit maneuvers in order to redirect the thruster to the desired propulsion direction. Although the proposed controller is concerning a single satellite, formation flying techniques have been utilized by assuming the satellite to be flying in formation with a virtual spacecraft (flying on the reference trajectory) with respect to which it is required to keep a predefined relative orbit (rendezvous in most cases). This approach allows astrodynamics’ insight into the classical linearized control problem to be leveraged. Moreover, this allows a direct implementation of the proposed algorithm in true multi-satellite relative-orbit-keeping contexts.

Designed to comprise a trade-off between fuel and time optimality through manipulating the controller gains, the proposed MPC couples the attitude and the relative orbital dynamics such that the control thrust is provided during the attitude redirection maneuvers. The stability of the control scheme is verified through Monte Carlo numerical simulations in which navigation errors, hardware errors, and physical constraints are emulated. In order to assess the controller’s performance, the proposed MPC has also been compared to a reference MPC from the literature. As for the compared scenario, the proposed MPC achieves more accurate orbit keeping in a shorter time frame. Such improved performance is achieved at the cost of a slightly larger delta-V cost, corresponding to approximately 24% of the total maneuver cost, and therefore is negligible at the mission level.