Piecewise Parallel Optimal Algorithm

Abstract

This chapter studies a new optimal algorithm that can be implemented in a piecewise parallel manner onboard spacecraft, where the capacity of onboard computers is limited. The proposed algorithm contains two phases. The predicting phase deals with the open-loop state trajectory optimization with simplified system model and evenly discretized time interval of the state trajectory. The tracking phase concerns the closed-loop optimal tracking control for the optimal reference trajectory with full system model subject to real space perturbations. The finite receding horizon control method is used in the tracking program. The optimal control problems in both programs are solved by a direct collocation method based on the discretized Hermite–Simpson method with coincident nodes. By considering the convergence of system error, the current closed-loop control tracking interval and next open-loop control predicting interval are processed simultaneously. Two cases are simulated with the proposed algorithm to validate the effectiveness of proposed algorithm. The numerical results show that the proposed parallel optimal algorithm is very effective in dealing with the optimal control problems for complex nonlinear dynamic systems in aerospace engineering area.

1. Introduction

Space tether system is a promising technology over decades. It has wide potential applications in the space debris mitigation & removal, space detection, power delivery, cargo transfer and other newly science & technic missions. Recently, there is continuous interest in the space tether systems, in leading space agencies such as, NASA’s US National Aeronautics and Space Administration, ESA’s European Space Agency, and JAXA’s Japan Aerospace Exploration Agency [1]. Their interest technologies include the electrodynamic tether (EDT) propulsion technology, retrieval of tethered satellite system, multibody tethered system and space elevator system. Compared with existing technologies adopted by large spacecraft such as the rocket or thruster, the space tether technology has the advantages of fuel-efficiency (little or no propellant required), compact size, low mass, and ease-of-use [2]. These advantages make it reasonable to apply the space tethered system for deorbiting the fast-growing low-cost micro/nano-satellites and no-fuel cargo transfer. The difficulty associated with space tether system is to control & suppress its attitudes during a mission process for the technology to be functional and practical. Many works have been devoted to solving this problem, and one effort is to use the optimal control due to its good performances in the complex and unstable nonlinear dynamic systems. In this chapter, a new piecewise onboard parallel optimal control algorithm is proposed to control and suppress the attitudes of the space tether system. To test its validity, two classical space tether systems, the electrodynamic tether system (EDT) and partial space elevator (PSE) system are considered and tested.

An EDT system with constant tether length is underactuated. The electric current is the only control input if there are no other active forces such as propulsion acting on the ends of an EDT. The commonly adopted control strategy in the literature is the current regulation using energy-based feedback in this underactuated control problem. Furthermore, many efforts have been done to solve this problem with optimal control. Stevens and Baker [3] studied the optimal control problem of the EDT libration control and orbital maneuverer efficiency by separating the fast and slow motions using an averaged libration state dynamics as constraints instead of instantaneous dynamic constraints in the optimal control algorithm. The instantaneous states are propagated from the initial conditions using the optimal control law in a piecewise fashion. Williams [4] treated the slow orbital and fast libration motions separately with two different discretization schemes in the optimal control of an EDT orbit transfer. The differential state equations of the libration motion are enforced at densely allocated nodes, while the orbital motion variables are discretized by a quadrature approach at sparsely allocated nodes. The two discretization schemes are unified by a specially designed node mapping method to reflect the coupling nature of orbital and libration motions. The control reference, however, is assumed known in advance.

A PSE system is consisted with one main satellite and two subsatellites (climber & end body) connected to each other by tether(s). The difficulty associated to such a system is to suppress the libration motion of the climber and the end body. This libration is produced by the moving climber due to the Coriolis force, which will lead the system unstable. While the climber is fast moving along the tether, the Coriolis force will lead to the tumbling of the PSE system. Thus, the stability control for suppressing such a system is critical for a successful climber transfer mission. To limit the fuel consuming, tension control is widely used to stable the libration motion of the space tethered system due to it can be realized by consuming electric energy only [5]. Many efforts have been devoted to suppressing the libration motion of space tethered system such as, Wen et al. [6] stabled the libration of the tethered system by an analytical feedback control law that accounts explicitly for the tension constraint. The study shows good computational effect, and the proposed method requires small data storage ability. Ma et al. [7] used adaptive saturated sliding mode control to suppress the attitude angle in the deployment period of the space tethered system. Optimal control [8, 9] is also proved as a way to overcome the libration issue. The above tension control schemes are helpful for both two-body and three-body tethered system. Up to data, limited devotions have been done on the libration suppression of a PSE system using tension control only. Williams used optimal control to design the climber’s speed function of a climber for a full space elevator [10]. Modeled by simplified dynamic equations, an optimal control problem is solved, and the solution results in zero in-plane libration motion of the ribbon in the ending phase of climber motion. The study shows that to eliminate the in-plane oscillations by reversing the direction of the elevator is possible. Kojima et al. [11] extended the mission function control method to eliminate the libration motion of a three body tethered system. The proposed method is effective when the total tether length is fixed and the maximum speed of the climber no more than 10 m/s. Although these efforts are useful to suppress the libration motion of the PSE system, it still difficult to control the attitudes of such a system in the transfer period.

To overcome the challenges in aforementioned works, we propose a parallel onboard optimal algorithm contains two phases. Phase 1 concerns the reference state trajectory optimization within a given time interval, where an optimal control model is formulated based on the timescale separation concept [3, 12] to simplify the dynamic calculations of the EDT & PSE system. An open-loop optimal state trajectory is then obtained by minimizing a cost function subject to given constraints. The state trajectory of paired state and control input variables is solved approximately by the direct collocation method [13] that is based on the Hermite-Simpson method [14]. In this phase, the simplified dynamic model is by used. Phase 2 concerns the tracking of the open-loop optimal state trajectory within the same interval. A closed-loop optimal control problem is formulated in a quadrature form to track the optimal state trajectory obtained in phase 1. Unlike phase 1, all the major perturbative forces are included, and more realistic geomagnetic and gravitational field models are considered. While the system is running the process in phase 2 with one CPU, the next phase 1 calculation is running in another CPU with data modification based on the errors obtained in the last calculation program. The simulation results demonstrate the effectiveness of the approach in fast satellite deorbit by EDTs in equatorial orbit. Furthermore, for fast transfer period of the partial space elevator, the propose method also shows good effect on suppression the libration angles of the climber and the end body with tension control only.

2. Optimal control algorithm

2.1. Control scheme

Assume two CUPs are used to process the calculation. CUP-1 is used to determine the open-loop optimal control trajectory of dynamic states employing the simple dynamic equations. The obtained optimal state trajectory will be tracked by CPU-2 using closed-loop RHC. While the system is tracking the i-th interval, the (i + 1)-th optimal trajectory is being calculated in CPU-1. Once the tracking for the i-th interval is finished (implemented by CPU-2), the real finial state S_i will be stored in the memory and the (i + 1)-th optimal trajectory can be tracked. By repeating the above process, the optimal suppression control problem is solved in a parallel piecewise manner until the transfer period is over.

The calculation of the optimal trajectory in CPU-1 is a prediction state trajectory whose initial state is estimated as ${\tilde{\mathit{S}}}_i^{i+1}={\tilde{\mathit{S}}}_i^i+\frac{1}{2}{\mathit{e}}_{i-1}$, where ${\tilde{\mathit{S}}}_i^i$ denotes the estimated finial state of the i-th interval obtained by CUP-1, the superscript and subscript denote the internal number and the states node number. ${\tilde{\mathit{S}}}_i^{i+1}$ is the estimated initial state of the (i + 1)-th interval, and ${\mathit{e}}_{i-1}={\mathit{S}}_{i-1}-{\tilde{\mathit{S}}}_{i-1}^{i-1}$ is the error between the real and the estimated final state of the (i-1)-th interval, the data used to calculate the error is picked up from the memory. For the first interval, i = 1, ${\mathit{S}}_0={\tilde{\mathit{S}}}_0^1$ and ${\mathit{e}}_0=0$. The computational diagram of the entire control strategy is given as shown in Figure 1.

Figure 1.

Control scheme.

2.2. Open-loop control trajectory

The libration angles of the climber and the end body are required to be kept between the desired upper/lower bounds in the climber transfer process. To make the calculation convenient and simple. The accessing process is divided into a series of intervals, such that, the transfer process $\left[t_i,t_{i+1}\right]$ is discretized into n intervals, where $t_i$ and $t_{i+1}$ are the initial and final time, respectively. $t_{i+1}$ can be obtained in terms of the transfer length of the climber and its speed. To make calculation convenient to be realized in practical condition, the transfer process is divided evenly. The optimal trajectory should be found to satisfy the desired cost functions for each time intervals as

$$J_i={\int }_{t_i}^{t_{i+1}}\Pi \left(\mathit{x},\mathit{u}\right)\mathit{dt}$$E1

subject to the simplified dynamic equations. All the errors between simple model and the entire model are regarded as perturbations. The above cost function minimization problem is solved by a direct solution method, which uses a discretization scheme to transform the continuous problem into a discrete parameter optimization problem of nonlinear programming within the interval, to avoid the difficulty usually encountered when standard approaches are used on derivation of the required conditions for optimality [15]. There are a number of efficient discretization schemes, such as, Hermite-Legendre-Gauss-Lobatto method [16] and Chebyshev pseudospectral method [8], in the literature for discretizing the continuous problem. In the current work, a direct collocation method, based on the Hermite-Simpson scheme [14, 17], is adopted because of its simplicity and accuracy.

Assume that the time interval $\left[t_i,t_{i+1}\right]$ is discretized into n subintervals with n + 1 nodes at the discretized time ${\tau }_k$ ($k=0,1,.\dots ,n$).

$${\gamma }_k={\tau }_{k+1}-{\tau }_k,\sum _{k=1}^n{\gamma }_k=t_{i+1}-t_i$$E2

The state vectors and control inputs are discretized at n + 1 nodes, ${\mathit{x}}_0$, ${\mathit{x}}_1$, ${\mathit{x}}_2$, …, ${\mathit{x}}_n$ and ${\upsilon }_0$, ${\upsilon }_1$, ${\upsilon }_2$, …, ${\upsilon }_n$. Further, denote the state vectors and the control inputs at mid-points between adjacent nodes by ${\mathit{x}}_{0.5}$, ${\mathit{x}}_{1.5}$, ${\mathit{x}}_{2.5}$, …, ${\mathit{x}}_{n-0.5}$ and ${\upsilon }_{0.5}$, ${\upsilon }_{1.5}$, ${\upsilon }_{2.5}$, …, ${\upsilon }_{n-0.5}$ and the mid-point state vectors, ${\mathit{x}}_{k+0.5}$, can be derived by the Hermite interpolation scheme,

$${\mathit{x}}_{k+0.5}=\frac{1}{2}\left({\mathit{x}}_k+{\mathit{x}}_{k+1}\right)+\frac{\gamma }{8}\left[\mathrm{\Gamma }\left({\mathit{x}}_k,{\upsilon }_k,{\tau }_k\right)-\mathrm{\Gamma }({\mathit{x}}_{k+1},{\upsilon }_{k+1},{\tau }_{k+1})\right]$$E3

Accordingly, the cost function in Eq. (1) can be discretized by the Simpson integration formula as

$$J\cong \frac{\gamma }{6}\sum _{k=0}^{n-1}\left[\mathrm{\Pi }\left({\mathit{x}}_k,{\upsilon }_k,{\tau }_k\right)+4\mathrm{\Pi }({\mathit{x}}_{k+0.5},{\upsilon }_{k+0.5},{\tau }_{k+0.5})+\mathrm{\Pi }({\mathit{x}}_{k+1},{\upsilon }_{k+1},{\tau }_{k+1})\right]$$E4

The nonlinear constraints based on the tether libration dynamics, the first-order states, can also be denoted by discretized equations using the Simpson integration formula, such that

$$\frac{\gamma }{6}\left[\mathrm{\Pi }\left({\mathit{x}}_k,{\upsilon }_k,t_k\right)+4\mathrm{\Pi }({\mathit{x}}_{k+0.5},{\upsilon }_{k+0.5},{\tau }_{k+0.5})+\mathrm{\Pi }({\mathit{x}}_{k+1},{\upsilon }_{k+1},{\tau }_{k+1})\right]+{\mathit{x}}_k-{\mathit{x}}_{k+1}=\mathit{0}$$E5

The left-hand side of Eq. (5) is also named as the Hermite-Simpson Defect vector in the literature. Finally, the discretization process is completed by replacing the constraints for the initial states and the continuous box constraints with the discretized constraints,

$${\mathit{x}}_0={\mathit{x}}_{\mathit{\text{start}}},{\mathit{x}}_{\min }\le {\mathit{x}}_k\le {\mathit{x}}_{\max },{\upsilon }_{\min }\le {\upsilon }_k\le {\upsilon }_{\max },{\upsilon }_{\min }\le {\upsilon }_{k+0.5}\le {\upsilon }_{\max }$$E6

The minimization problem of a continuous cost function is now transformed to a nonlinear programming problem. It searches optimal values for the programming variables that minimize the discretized form of cost function shown in Eq. (4) while satisfying the constraints of Eqs. (5) and (6). The subscript index “k” will be refreshed in the next time interval.

2.3. Closed-loop optimal control for tracking open-loop optimal state trajectory

The RHC is implemented by converting the continuous optimal control problem into a discrete parameter optimization problem that can be solved analytically. Like the open-loop trajectory optimization problem, the same direct collocation based on the Hermite-Simpson method is used to discretize the RHC problem.

By using the similar notations of discretization above, the cost function is discretized using the Simpson integration formula as

$$\begin{array}{c}G=\frac{1}{2}\delta {\mathit{x}}_{i+1}^T\mathit{S}\delta {\mathit{x}}_{i+1}\\ +\frac{\gamma }{12}\sum _{k=0}^{n-1}\left[\delta {\mathit{x}}_k^T\mathit{Q}\delta {\mathit{x}}_k+4\delta {\mathit{x}}_{k+0.5}^T\mathit{Q}\delta {\mathit{x}}_{k+0.5}+\delta {\mathit{x}}_{k+1}^T\mathit{Q}\delta {\mathit{x}}_{k+1}+\mathit{R}\left(\delta {\mathbf{\upsilon }}_k^2+4\delta {\mathbf{\upsilon }}_{k+0.5}^2+\delta {\mathbf{\upsilon }}_{k+1}^2\right)\right]\end{array}$$E7

and the constraints are discretized into

$$\delta {\mathit{x}}_k-\delta {\mathit{x}}_{k+1}+\frac{\gamma }{6}\left[{\mathit{A}}_k\delta {\mathit{x}}_k\mathit{+}{B}_k{\mathit{\delta \upsilon }}_k+4{\mathit{A}}_{k+0.5}\delta {\mathit{x}}_{k+0.5}+4{\mathit{B}}_{k+0.5}{\mathit{\delta \upsilon }}_{k+0.5}+{\mathit{A}}_{k+1}\delta {\mathit{x}}_{k+1}\mathit{+}{\mathit{B}}_{k+1}{\mathit{\delta \upsilon }}_{k+1}\right]=\mathit{0}$$E8

$$\delta {\mathit{x}}_k=\mathit{x}\left({\tau }_k\right)-{\mathit{x}}_{\mathit{opt}}\left({\tau }_k\right)$$E9

$$\delta {\mathit{x}}_{k+0.5}=\frac{1}{2}\left(\delta {\mathit{x}}_k+\delta {\mathit{x}}_{k+1}\right)+\frac{\gamma }{8}\left({\mathit{A}}_k\delta {\mathit{x}}_k+{\mathit{B}}_k{\mathit{\delta \upsilon }}_k-{\mathit{A}}_{k+1}\delta {\mathit{x}}_{k+1}-{\mathit{B}}_{k+1}{\mathit{\delta \upsilon }}_{k+1}\right)$$E10

where, ${\mathit{A}}_k=\mathit{A}\left({\tau }_k\right)$, ${\mathit{A}}_{k+0.5}=\mathit{A}\left({\tau }_{k+0.5}\right)$, ${\mathit{B}}_k=\mathit{B}\left({\tau }_k\right)$, ${\mathit{B}}_{k+0.5}=\mathit{B}\left({\tau }_{k+0.5}\right)$.

The derivation of Eq. (8a) finally leads to a quadratic programming problem to find a programming vector $\mathit{Z}={\left[\delta {\mathit{x}}_0^T\delta {\mathit{x}}_1^T\dots \delta {\mathit{x}}_n^T{\mathit{\delta \upsilon }}_0{\mathit{\delta \upsilon }}_1\dots {\mathit{\delta \upsilon }}_n{\mathit{\delta \upsilon }}_{0.5}{\mathit{\delta \upsilon }}_{1.5}\dots {\mathit{\delta \upsilon }}_{n-0.5}\right]}^T$, which minimizes the cost function:

$$G=\frac{1}{2}{\mathit{Z}}^T\mathit{MZ}$$E11

subject to.

$$\mathit{CZ}\mathit{=}\mathit{X} \mathit{X}={\left[\delta {\mathit{x}}_0^{T}00\dots 0\right]}^T$$E12

where the matrices C and M are given in the Appendix.

It is easy to find the solution analytically to this standard quadratic programming problem by [24].

$${\mathit{Z}}^{\ast }\mathit{=}{\mathit{M}}^{-1}{\mathit{C}}^T{\left({\mathit{CM}}^{-1}{\mathit{C}}^T\right)}^{-1}\mathit{X}$$E13

and the control correction at the current time can be obtained as

$$\mathit{\delta \upsilon }\left(t_i\right)\mathit{=}{\mathit{VZ}}^{\ast }\mathit{=}{\mathit{VM}}^{-1}{\mathit{C}}^T{\left({\mathit{CM}}^{-1}{\mathit{C}}^T\right)}^{-1}\mathit{X}\triangleq \mathit{K}\left(t_i,n,t_h\right)\left[\mathit{x}\left(t_i\right)-{\mathit{x}}_{\mathit{opt}}\left(t_i\right)\right]$$E14

where the row vector $\mathit{V}$ is defined to “choose” the target value from the optimal solution, and the position of “1” in the row vector $\mathit{V}$ is the same as the position of ${\mathit{\delta \upsilon }}_0$ in the column vector $\mathit{Z}$. Finally, the control input of the closed-loop control, $\upsilon \left(t_i\right)$, is

$$\upsilon \left(t_i\right)={\upsilon }_{\mathit{opt}}\left(t_i\right)+\mathit{\delta \upsilon }\left(t_i\right)={\upsilon }_{\mathit{opt}}\left(t_i\right)+\mathit{K}\left(t_i,n,t_h\right)\left[\mathit{x}\left(t_i\right)-{\mathit{x}}_{\mathit{opt}}\left(t_i\right)\right]$$E15

It is apparent that the closed-loop control law derived here is a linear proportional feedback control law, and the feedback gain matrix $\mathit{K}$ is a function of time. Without any explicit integration of differential equations, $\mathit{K}$ can either be determined offline or online depending on the computation and restoration ability onboard the satellite.

It is worth to point out some advantages of this approach. Firstly, the matrices M and C are both formulated by the influence matrices at certain discretization notes (A_k, B_k, A_k + 0.5, B_k + 0.5). If t_i and t_i are both set to be coincident with the discretization nodes used in the open-loop control problem, then most of the influence matrices calculated previously can be used directly in the tracking control process to reduce computational efforts. This is the advantage of using the same discretization method in the current two-phased optimal control approach. Secondly, the matrix $\mathit{M}$ is unchanged if treating the terminal horizon time $t_i+t_h$ as the time-to-go and keep the future horizon interval $\left[t_i,t_{i+1}\right]$ unchanged. This means the inverse of $\mathit{M}$ could be calculated only once in the same interval. It is attractive for the online implementation of HRC, where the computational effort is critical. As the entire interval can be discretized into small intervals, if these intervals are sufficiently small relative to the computing power of the satellite, then the calculation process can be carried by the onboard computer. Furthermore, for small intervals, the computation for the open-loop optimal trajectory of the next interval can be done by CPU-1 while the tracking is still in process, see Figure 1. This makes of the proposed optimal suppression control a parallel online implementation, which is another advantage of this control scheme.

3. Cases study

3.1. Parallel optimal algorithm in attitudes control of EDT system

In order to test the validity of the proposed optimal algorithm, a case study of the attitudes control of EDT system in aerospace engineering is used. The obtained results are compared with some existing control methods.

3.1.1. Problem formulation

The EDT system’s orbital motion is generally described in an Earth geocentric inertial frame (OXYZ) with the origin O at the Earth’s centre, see Figure 2(a). The X-axis directs to the point of vernal equinox, the Z-axis aligns with the Earth’s rotational axis, and the Y-axis completes a right-hand coordinate system, respectively. The equation of orbital dynamic motion can be written in the form of Gaussian perturbation [18], which is a set of ordinary differential equations of six independent orbital elements ($a$, $\Omega$, $i$, $e_x$, $e_y$, $\phi$)

$$\frac{\mathit{da}}{\mathit{dt}}=\frac{2}{n\sqrt{1-e^2}}\left({\sigma }_{z}esin\nu +{\sigma }_x\frac{p}{r}\right)$$E16

$$\frac{d\mathrm{\Omega }}{\mathit{dt}}=\frac{{\sigma }_{y}rsinu}{{\mathit{na}}^2\sqrt{1-e^2}sini}$$E17

$$\frac{\mathit{di}}{\mathit{dt}}=\frac{{\sigma }_{y}rcosu}{{\mathit{na}}^2\sqrt{1-e^2}}$$E18

$$\frac{{\mathit{de}}_x}{\mathit{dt}}=\frac{\sqrt{1-e^2}}{\mathit{na}}\left\{{\sigma }_{z}sinu+{\sigma }_x\left[\left(1+\frac{r}{p}\right)cosu+\frac{r}{p}{e}_x\right]\right\}+\frac{d\mathrm{\Omega }}{\mathit{dt}}{e}_{y}cosi$$E19

$$\frac{{\mathit{de}}_y}{\mathit{dt}}=\frac{\sqrt{1-e^2}}{\mathit{na}}\left\{-{\sigma }_{z}cosu+{\sigma }_x\left[\left(1+\frac{r}{p}\right)sinu+\frac{r}{p}{e}_y\right]\right\}-\frac{d\mathrm{\Omega }}{\mathit{dt}}{e}_{x}cosi$$E20

$$\frac{\mathit{d\phi }}{\mathit{dt}}=n-\frac{1}{\mathit{na}}\left[{\sigma }_z\left(\frac{2r}{a}+\frac{\sqrt{1-e^2}}{1+\sqrt{1-e^2}}ecos\nu \right)-{\sigma }_x\left(1+\frac{r}{p}\right)\frac{\sqrt{1-e^2}}{1+\sqrt{1-e^2}}esin\nu \right]-\frac{{\sigma }_{y}rcosisinu}{{\mathit{na}}^2\sqrt{1-e^2}sini}$$E21

Figure 2.

Illustration of coordinate system for the EDT’s orbital (a) and libration (b) motion.

The components of perturbative accelerations are defined in a local frame. The components ${\sigma }_x$ and ${\sigma }_z$ are in the orbital plane, and ${\sigma }_z$ is the radial component pointing outwards. The out-of-plane component ${\sigma }_y$ completes a right-hand coordinate system. The components of perturbative accelerations depend on the tether attitude and the EDT’ orbital dynamics is coupled with the tether libration motion.

The libration motion of a rigid EDT system is described in an orbital coordinate system shown in Figure 2(b). The z-axis of the orbital coordinate system points from the Earth’s center to the CM of the EDT system, the x-axis lies in the orbital plane and points to the direction of the EDT orbital motion, perpendicular to the z-axis. The y-axis completes a right-hand coordinate system. The unit vectors along each axis are expressed as ${\overrightarrow{\mathit{e}}}_{\mathit{ox}}$, ${\overrightarrow{\mathit{e}}}_{\mathit{oy}}$, and ${\overrightarrow{\mathit{e}}}_{\mathit{oz}}$, respectively. Then, the instantaneous attitudes of the EDT system are described by an in-plane angle $\alpha$ (pitch angle, rotating about the y-axis) followed by an out-of-plane angle $\beta$ (roll angle, rotating about the x’-axis, the x-axis after first rotating about the y-axis). Thus, the equations of libration motion of the EDT system can be derived as,

$$\ddot{\alpha }+\ddot{\nu }-2\left(\dot{\alpha }+\dot{\nu }\right)\dot{\beta }tan\beta +3\mu r^{-3}sin\alpha cos\alpha =\frac{Q_{\alpha }}{\tilde{m}{L}^2{cos}^2\beta }$$E22

$$\ddot{\beta }+{\left(\dot{\alpha }+\dot{\nu }\right)}^{2}sin\beta cos\beta +3\mu r^{-3}{cos}^2\alpha sin\beta cos\beta =\frac{Q_{\beta }}{\tilde{m}{L}^2}$$E23

where $\tilde{m}=\left[m_1{m}_2+\left(m_1+m_2\right)m_t/3+m_t^2/12\right]m_{\mathit{EDT}}^{-1}$ is the equivalent mass, ($Q_{\alpha }$, $Q_{\beta }$) are the corresponding perturbative torques by the perturbative forces to be discussed below.

The perturbative accelerations (${\sigma }_x$, ${\sigma }_y$, ${\sigma }_z$) and torques ($Q_{\alpha }$, $Q_{\beta }$) in Eq. (15) are induced by multiple orbital perturbative effects, namely, (i) the electrodynamic force exerting on a current-carrying EDT due to the electromagnetic interaction with the geomagnetic field, (ii) the Earth’s atmospheric drag, (iii) the Earth’s non-homogeneity and oblateness, (iv) the lunisolar gravitational perturbations, and (v) the solar radiation pressure, respectively. The EDT system is assumed thrust-less during the deorbit process, while the atmosphere, geomagnetic and ambient plasma fields are assumed to rotate with the Earth at the same rate. The geodetic altitude, instead of geocentric altitude, should be used in the evaluation of the environmental parameters, such as, atmospheric and plasma densities, to realistically account for the Earth’s ellipsoidal surface, such that,

$$h_g=r-r_{\mathit{po}}{\left(1-e_E^2{cos}^2\theta \right)}^{-1/2}$$E24

where the polar radius $r_{\mathit{po}}$ and the Earth’s eccentricity $e_E$ are provided by NASA [19].

Moreover, the local strength of geomagnetic field is described by the IGRF2000 model [20, 21, 22] in a body-fixed spherical coordinates of the Earth, such that

$$\begin{array}{l}{B}_{\varphi }=\frac{1}{sin\theta }\sum _{n=1}^{\infty }{\left(\frac{r_0}{r}\right)}^{n+2}\sum _{m=0}^{n}m\left[g_n^{m}sin\left(\mathit{m\varphi }\right)-h_n^{m}cos\left(\mathit{m\varphi }\right)\right]P_n^m\left({\theta }_c\right)\\ B_{\theta }=\sum _{n=1}^{\infty }{\left(\frac{r_0}{r}\right)}^{n+2}\sum _{m=0}^n\left[g_n^{m}cos\left(\mathit{m\varphi }\right)+h_n^{m}sin\left(\mathit{m\varphi }\right)\right]\frac{\partial P_n^m\left({\theta }_c\right)}{\partial {\theta }_c}\\ B_r=\sum _{n=1}^{\infty }{\left(\frac{r_0}{r}\right)}^{n+2}\left(n+1\right)\sum _{m=0}^n\left[g_n^{m}cos\left(\mathit{m\varphi }\right)+h_n^{m}sin\left(\mathit{m\varphi }\right)\right]P_n^m\left({\theta }_c\right)\end{array}$$E25

where r₀ = 6371.2 × 10³ km is the reference radius of the Earth, respectively.

The average current in the EDT is defined as

$$I_{\mathit{ave}}=\frac{1}{L}{\int }_0^{L}I\left(s\right)\mathit{ds}$$E26

The open-loop optimal control problem for EDT deorbit can be stated as finding a state-control pair $\left\{\mathit{x}\left(t\right),\upsilon \left(t\right)\right\}$ over each time interval $\left[t_i,t_{i+1}\right]$ to minimize a cost function of the negative work done by the electrodynamic force

$$J={\int }_{t_i}^{t_{i+1}}{\overrightarrow{\mathit{F}}}_e·\overrightarrow{\mathit{v}}\mathit{dt}\triangleq {\int }_{t_t}^{t_{i+1}}\mathrm{\Pi }\left(\mathit{x},\upsilon ,t\right)\mathit{dt}$$E27

subject to the nonlinear state equations of libration motion

$$\begin{array}{ll}{x}_1^{\prime }& =x_2\\ x_2^{\prime }& =2{\eta }^{3}esin\nu +2\left(x_2+{\eta }^2\right)x_{4}tan{x}_3-3{\eta }^{3}sin{x}_{1}cos{x}_1\\ & +\left[sinitan{x}_3\left(2sinucos{x}_1-cosusin{x}_1\right)-cosi\right]\left(\varsigma -\lambda \right)I_{\mathit{ave}}\frac{{\mu }_m}{\mu \tilde{m}}{\eta }^3\end{array}$$E28

$$\begin{array}{cc}{x}_3^{\prime }& =x_4\\ x_4^{\prime }& =-{\left(x_2+{\eta }^2\right)}^{2}sin{x}_{3}cos{x}_3-3{\eta }^3{cos}^2{x}_{1}sin{x}_{3}cos{x}_3\\ & -sini\left(2sinusin{x}_1+cosucos{x}_1\right)\left(\varsigma -\lambda \right)I_{\mathit{ave}}\frac{{\mu }_m}{\mu \tilde{m}}{\eta }^3\end{array}$$E29

where $\left(x_1,x_2,x_3,x_4\right)=\left(\alpha ,{\alpha }^{\prime },\beta ,{\beta }^{\prime }\right)$, $\eta =1+ecos\nu$, $\dot{\nu }={\left(\mu /p^3\right)}^{0.5}{\left(1+ecos\nu \right)}^2$, $r=p{\left(1+ecos\nu \right)}^{-1}$, $\lambda =\left(m_1+0.5m_t\right)/m_{\mathit{EDT}}$ is determined by the mass ratio between the end-bodies, and $\varsigma$ is determined by the distribution of current along the EDT, such that, $\varsigma =I_{\mathit{ave}}^{-1}{L}^{-2}{\int }_0^L\mathit{sI}\left(s\right)\mathit{ds}$. Accordingly, $\varsigma =0.5$ is used for the assumption of a constant current in the EDT. The initial conditions $\mathit{x}\left(t_i\right)={\mathit{x}}_{\mathit{\text{start}}}$ and the box constraint $\left|\alpha \right|\le {\alpha }_{\max },\left|\beta \right|\le {\beta }_{\max },I_{\min }\le I_{\mathit{ave}}\le I_{\max }$. The environmental perturbations are simplified by considering only the electrodynamic force with a simple non-tilted dipole model of geomagnetic field, such that,

$$\overrightarrow{\mathbf{B}}=\frac{{\mu }_m}{r^3}cosusini{\stackrel{\rightharpoonup }{\mathit{e}}}_{\mathit{ox}}+\frac{{\mu }_m}{r^3}cosi{\stackrel{\rightharpoonup }{\mathit{e}}}_{\mathit{oy}}-\frac{2{\mu }_m}{r^3}sinusini{\stackrel{\rightharpoonup }{\mathit{e}}}_{\mathit{oz}}$$E30

Accordingly, the electrodynamic force ${\overrightarrow{\mathit{F}}}_e$ exerting on the EDT can be obtained as,

$${\overrightarrow{\mathbf{F}}}_e=-{\int }_0^L\stackrel{\rightharpoonup }{\mathbf{B}}\times I\stackrel{\rightharpoonup }{\mathit{l}}\mathit{ds}=-I_{\mathit{ave}}L\stackrel{\rightharpoonup }{\mathbf{B}}\times \stackrel{\rightharpoonup }{\mathit{l}}$$E31

3.1.2. Results and discussion

The initial and boundary conditions of box constraints of the case are shown in Tables 1 and 2.

Parameters	Values
Mass of the main satellite	5 kg
Mass of subsatellite	1.75 kg
Mass of the tether	0.25 kg
Dimensions of main satellite	0.2 × 0.2 × 0.2 m
Dimensions of subsatellite	0.1 × 0.17 × 0.1 m
Tether length	500 m
Tether diameter	0.0005 m
Tether conductivity (aluminum)	3.4014 × 10^7 Ω^−1 m^−1
Tether current lower/upper limits	0 ~ 0.8 A for the equatorial orbit
Orbital altitudes	700 ~ 800 km

Table 1.

Parameters of an EDT system.

Parameters	Values
$I_{\max }$ (equatorial)	0.4 A
$I_{\max }$ (inclined)	$0.1\left[-sinisin{\overline{\beta }}_{\mathit{B}}sin\left({\mathrm{\Omega }}_G+{\overline{\alpha }}_{\mathit{B}}-\mathrm{\Omega }\right)+cosicos{\overline{\beta }}_{\mathit{B}}\right]{cos}^{-1}\left(i-{\overline{\beta }}_{\mathit{B}}\right)$ A
$I_{\min }$ (equatorial & inclined)	0 A
${\alpha }_{\max }$ (equatorial & inclined)	45 degrees
${\beta }_{\max }$ (equatorial & inclined)	45 degrees

Table 2.

Boundary Values of Box Constraints in the Open-Loop Trajectory Optimization.

Firstly, the validity of the proposed optimal control scheme in the equatorial orbit where the EDT system gets the highest efficiency is demonstrated. The solid line in Figure 4 shows the time history of the EDT’s average current control trajectory obtained from the open-loop optimal control problem. It is clearly shown that the average current in the open-loop case reaches the upper limit most of the time, which indicates the electrodynamic force being maximized for the fast deorbit. As expected, the current is not always at the upper limit in order to avoid the tumbling of the EDT system. This is evident that the timing of current reductions coincides with the peaks of pitch angles shown in Figure 5. The effectiveness of the proposed control scheme in terms of keeping libration stability is further demonstrated by the solid lines in Figure 5, where the trajectory of libration angles is no more than 45 degrees. It is also found that the amplitude of the pitch angle nearly reaches 45 degrees, the maximum allowed value, whereas the roll angle is very small in the whole deorbit process. As a comparison, tracking optimal control with the non-tilted dipole model and the IGRF 2000 model of the geomagnetic fields are conducted respectively.

The dashed lines in Figures 3 and 4 show the tracking control simulations where all perturbations mentioned before are included with the non-tilted dipole geomagnetic field model. As expected, the closed-loop tracking control works well in this case since the primary electrodynamic force perturbation is the same as the one used in the open-loop trajectory optimization. It is shown clearly in Figure 4 that the pitch angle under the proposed closed-loop control tracks the open-loop optimal trajectory very closely with this simple environment model. Figure 4 also shows the roll angle is almost zero even if it is not tracked. At the same time, Figure 3 shows that the current control modification to the optimal current trajectory is relative small, i.e., 12% above the maximum current, for the same reason. Now the same cases are analyzed again using a more accurate geomagnetic field model – the IGRF 2000 model with up to 7th order terms (Figures 5 and 6). The solid line in Figure 5 is the open-loop current control trajectory while the dashed line is the modified current control input obtained by the receding horizon control. Compared with Figure 3, it shows more current control modifications are needed to track the open-loop control trajectory because of larger differences in dynamic models between the open-loop and closed-loop optimal controls, primarily due to the different geomagnetic field models. Because of the same reason, it is noticeable in Figure 6 that the instantaneous states of the EDT system, controlled by the closed-loop optimal control law, are different from the open-loop reference state trajectory at the end of the interval. The instantaneous states are used for the next interval as initial conditions for the open-loop optimal control problem to derive the optimal control trajectory in that interval. This is reflected in Figure 6 that the solid lines are discontinuous at the beginning of each interval. The dashed line in Figure 7 shows that the pitch and roll angle under the closed-loop control has been controlled to the open-loop control trajectory, indicating the effectiveness of the proposed optimal control law. The roll angle is not controlled in this case as mentioned before. Compared Figure 4 with Figure 6, it shows that the roll angle increases significantly since there is an out-of-plane component of the electrodynamic force resulting from the IGRF 2000 geomagnetic model. However, the amplitude of the roll angle is acceptable within the limits and will not lead to a tumbling of the EDT system.

Figure 3.

Time history of average current in the equatorial orbit (non-tilted dipole geomagnetic model). Solid line: Open-loop state trajectory. Dashed line: Close-loop tracking trajectory.

Figure 4.

Time history of pitch and roll angles in the equatorial orbit (non-tilted dipole geomagnetic model). Solid line: Open-loop state trajectory. Dashed line: Close-loop tracking trajectory.

Figure 5.

Time history of average current in the equatorial orbit (IGRF 2000 model). Solid line: Open-loop state trajectory. Dashed line: Close-loop tracking trajectory.

Figure 6.

Time history of pitch and roll angles in the equatorial orbit (IGRF 2000 model). Solid line: Open-loop control trajectory. Dashed line: Close-loop tracking trajectory.

Figure 7.

Comparison of EDT deorbit rates using different control laws and geomagnetic field models.

Finally, we make a comparison to show the performance of the proposed onboard parallel optimal control law from the aspect of deorbit rate. A simple current on–off control law from a previous work of Zhong and Zhu [21] is used here as baselines for the comparison of EDT deorbit efficiency. The current on–off control becomes active only if the libration angles exceed the maximum allowed values. Furthermore, it will turn on the current only in the condition that the electrodynamic force does negative work in both pitch and roll directions. In this paper, the maximum allowed amplitude for pitch and roll angles was set to 20° and the turned-on current was assumed to be 0.4 A, roughly the average value of the current control input into the closed-loop optimal control. Besides, a minimum interval of 10 minutes for the switching was imposed to avoid equipment failure that might happen due to the frequent switching. Figure 8 shows the comparisons of the deorbit rates in different cases (the present optimal control and the current switching control with the non-tilted dipole or the IGRF 2000 model of geomagnetic field). It is shown that the EDT deorbit under the proposed optimal control scheme is faster than the current on–off control regardless which geomagnetic field model is used. The deorbit time of proposed optimal control based on the IGRF2000 model is about 25 hours, which equals approximately 15 orbits, whereas the deorbit time of simple current on–off control based on the same geomagnetic field model is about 55 hours, which equals approximately 33 orbits. The results also indicate that in the optimal control scheme, the effect is mostly shown in the current control input, instead of the deorbit rate, where Figure 6 shows much more current control effort is required due to the different magnetic field models were used in the open-loop control trajectory optimization and the closed-loop optimal tracking control.

Figure 8.

PSE system.

3.2. Parallel optimal algorithm in libration suppression of partial space elevator

For further test of the effect of the onboard parallel algorithm, the proposed control method is used to suppress the libration motions of the partial space elevator system. As studied in [24], this system is a non-equilibrium nonlinear dynamic system. It is difficult to suppress such a system in the mission period by using the common control design methods. In this case, we mainly concern obtaining the local time optimization.

3.2.1. Problem formulation

Consider an in-plane PSE system in a circular orbit is shown in Figure 8, where the main satellite, climber and the end body are connected by two inelastic tethers $L_1$ and $L_2$, respectively. The masses of the tethers are neglected. Assuming the system is subject to a central gravitational field and orbiting in the orbital plane. All other external perturbations are neglected. The main satellite, climber and end body are modeled as three point masses ($M$, $m_1$ and $m_2$) since the tether length is much greater than tethered bodies [5, 23]. Thus, the libration motions can be expressed in an Earth inertial coordinate system OXY with its origin at the centre of Earth. Denoting the position of the main satellite ($M$) by a vector $\mathbf{r}$ measuring from the centre of Earth. The climber $m_1$ is connected to the main satellite $M$ by a tether 1 with the length of $L_1$ and a libration angle ${\theta }_1$ measured from the vector $\mathbf{r}$. The distance between them is controlled by reeling in/out tether 1 at main satellite. The end body $m_2$ is connected to $m_1$ by a tether 2 with the length of $L_2$ and a libration angle ${\theta }_2$ measured from the vector $\mathbf{r}$. The length of tether 2 $L_2$ is controlled by reeling in or out tether 2 at end body. The mass of the main satellite is assumed much greater than the masses of the climber and the end body. Therefore, the CM of the PSE system can be assumed residing in the main satellite that moves in a circular orbit. Based on the aforementioned assumptions, the dynamic equations can be written as.

$${\ddot{\theta }}_1=-\frac{3{\omega }^{2}sin2{\theta }_1}{2}-\frac{2\left(\omega +{\dot{\theta }}_1\right){\dot{L}}_1}{L_1}-\frac{sin\left({\theta }_1-{\theta }_2\right)T_2}{L_1{m}_1}$$E32

$${\ddot{\theta }}_2=\frac{-3{\omega }^{2}sin2{\theta }_2}{2}-\frac{2\left(\omega +{\dot{\theta }}_2\right)\left({\dot{L}}_c-{\dot{L}}_1\right)}{L_0-L_1+L_c}+\frac{sin\left({\theta }_1-{\theta }_2\right)T_1}{\left(L_0-L_1+L_c\right)m_1}$$E33

$${\ddot{L}}_1=3{\omega }^2{L}_1{cos}^2{\theta }_1+2\omega L_1{\dot{\theta }}_1+L_1{{\dot{\theta }}_1}^2-\frac{T_1}{m_1}+\frac{cos\left({\theta }_1-{\theta }_2\right)T_2}{m_1}$$E34

$$\begin{array}{l}{\ddot{L}}_c=3{\omega }^2\left(L_0-L_1+L_c\right){cos}^2{\theta }_2+3{\omega }^2{L}_1{cos}^2{\theta }_1+\left(2\omega +{\dot{\theta }}_2\right)\left(L_0-L_1+L_c\right){\dot{\theta }}_2\\ +\left(2\omega +{\dot{\theta }}_1\right)L_1{\dot{\theta }}_1+\frac{\left[cos\left({\theta }_1-{\theta }_2\right)-1\right]T_1}{m_1}-\frac{\left[m_1-m_{2}cos\left({\theta }_1-{\theta }_2\right)+m_2\right]T_2}{m_1{m}_2}\end{array}$$E35

where $L_0$ is the initial total length of two pieces of the tethers and $L_c$ is the length increment relates to $L_0$.

The libration angles are required to be kept between the desired upper/lower bounds in the climber transfer process. The accessing process is divided into a series of intervals. In this case, we modified the aforementioned parallel optimal algorithm. The total transfer length $L_{10}-L_{1f}$ of the climber is discretized evenly. The optimal trajectory should be found to make the transfer time minimize in each equal tether transferring length. Then the cost function can be rewritten as

$$J_i=t_f$$E36

subject to the simplified dynamic equations

$${\ddot{\theta }}_1=-3{\omega }^2{\theta }_1-\frac{2\left(\omega +{\dot{\theta }}_1\right){\dot{L}}_1}{L_1}-\frac{\left({\theta }_1-{\theta }_2\right)T_2}{L_1{m}_1}$$E37

$${\ddot{\theta }}_2=-3{\omega }^2{\theta }_2-\frac{2\left(\omega +{\dot{\theta }}_2\right)\left({\dot{L}}_c-{\dot{L}}_1\right)}{L_0-L_1+L_c}+\frac{\left({\theta }_1-{\theta }_2\right)T_1}{\left(L_0-L_1+L_c\right)m_1}$$E38

$${\ddot{L}}_1=3{\omega }^2{L}_1+2\omega L_1{\dot{\theta }}_1+L_1{{\dot{\theta }}_1}^2+\frac{T_2-T_1}{m_1}$$E39

$${\ddot{L}}_c=3{\omega }^2\left(L_0+L_c\right)+\left(2\omega +{\dot{\theta }}_2\right)\left(L_0-L_1+L_c\right){\dot{\theta }}_2+\left(2\omega +{\dot{\theta }}_1\right)L_1{\dot{\theta }}_1-\frac{T_2}{m_2}$$E40

where $\mathit{u}=\left(T_1,T_2\right)$, $\mathit{x}=\left({\theta }_1,{\dot{\theta }}_1,{\theta }_2,{\dot{\theta }}_2,L_1,{\dot{L}}_1\right)$ and i denotes the interval number. In (26) the gravitational perturbations and the trigonometric functions are ignored, then they can be simplified following the assumptions: $sin{\theta }_j\sim {\theta }_j,cos{\theta }_j\sim 1\left(j=1,2\right)$. All the errors between simple model and the entire model are regarded as perturbations.

To ensure the availability and the suppression of the libration angles, following constrains are also required to be subjected $0\le T_1\le T_{1\text{max}},0\le T_2\le T_{2\text{max}},\left|{\theta }_1\right|\le {\theta }_{1\text{max}},\left|{\theta }_2\right|\le {\theta }_{2\text{max}},\left|L_c\right|\le L_{\mathit{\text{cLimit}}}$${\dot{L}}_{1m}\le {\dot{L}}_1\le {\dot{L}}_{1M},{\dot{L}}_{\mathit{cm}}\le {\dot{L}}_c\le {\dot{L}}_{\mathit{cM}}$, where T_1max and T_2max are the upper bounds of the tension control inputs T₁ and T₂, respectively. θ_1max, θ_2max and L_cLimit are the magnitudes of libration angles θ₁, θ₂ and the maximum available length scale of $L_c$, respectively. ${\dot{L}}_{1m}$ and ${\dot{L}}_{1M}$ are the lower and upper bounds of climber’s moving speed ${\dot{L}}_1$, respectively. ${\dot{L}}_{\mathit{cm}}$ and ${\dot{L}}_{\mathit{cM}}$ are the lower and upper bounds of end-bodies’ moving speed ${\dot{L}}_c$, respectively. It should be noting that, to avoid the tether slacking, the control tensions are not allowed smaller than zero. Dividing the time interval $\left[t_i,t_{i+1}\right]$ evenly into n subintervals. The cost function minimization problem for each time interval can be solved by Hermite–Simpson method, due its simplicity and accuracy [17]. Then nonlinear programming problem is to search optimal values for the programming variables that minimize the cost function for each interval shown in (25). The closed-loop optimal tracking control method, is same as that in case 1. Direct transcription methods are routinely implemented with standard nonlinear programming (NLP) software. The sparse sequential quadratic programming software SNOPT is used via a MATLAB-executable file interface.

3.2.2. Results and discussion

The proposed control scheme is used to suppress the libration angles of the PSE system in the ascending process with following system parameters and initial conditions: r = 7100 km, m₁ = 500 kg, m₂ = 1000 kg, θ₁(0) = θ₂(0) = 0, L₀ = 20 km, L₁(0) = 19,500 m, L_c (0) = 0, ${\dot{\theta }}_1\left(0\right)={\dot{\theta }}_2\left(0\right)=0$, ${\dot{L}}_1\left(0\right)=-20m/s$, and ${\dot{L}}_c\left(0\right)=0$ for the ascending process. The whole transfer trajectory is divided into 50 intervals. The climber’s ascending speed along tether 1 is allowed to be controlled to help suppress the libration angles and keep the states of the system in an acceptable area. The constrains are set as $T_{1\text{max}}=T_{2\text{max}}=200N$, ${\theta }_{1\text{max}}={\theta }_{2\text{max}}=0.3\mathit{rad}$, ${\dot{L}}_{1m}=-15m/s$, ${\dot{L}}_{1M}=-25m/s$, ${\dot{L}}_{\mathit{cm}}=-10m/s$ and ${\dot{L}}_{\mathit{cM}}=10m/s$.

The simulation results of this case are shown in Figures 9–11. The climber’s open-loop libration angle approaches its upper bound at 850 s. After 850 s${\theta }_1$ is kept at 0.3 rad by the end of the ascending period, see the dashed line in Figure 10. Using the closed-loop control, the tracking trajectory of ${\theta }_1$ matches the open-loop trajectory very well overall, see solid line in Figure 9. A short gap appears between 875 s – 880 s, this is caused by the errors of the model and computation. Figure 9 also shows the changes of the trajectories of ${\theta }_2$.The trajectory of ${\theta }_2$ obtained by closed-loop control tracks the open-loop trajectory well and reaches 0.1 rad by the end of the ascending period. The closed-loop trajectories of $L_1$ and $L_c$ are shown in Figure 10. They are the reflections of the control inputs. Both $L_1$ and $L_c$ show smooth fluctuations between 40s and 350 s. Figure 10 shows the time history of trajectories of ${\dot{L}}_1$ and ${\dot{L}}_c$, respectively. In the first 140 s the trajectory of ${\dot{L}}_c$ increases continuously until reaches its upper bound. Then, it keeps at 10 m/s by 270 s with some slight fluctuations. From 270 s to 355 s, it reduces continuously to −3 m/s. After that, ${\dot{L}}_c$ fluctuates around −3 m/s by the end of the transfer period. As a reflection of the control input, ${\dot{L}}_1$ also shows fluctuation during the transfer phase with obvious small-scale fluctuations appear in the period of 120 s – 260 s and 360 s – 750 s. This impacts during the whole transfer period, the changeable speed of the climber has the ability to help the suppression of the libration angles and states trajectory tracking. This time history of the control inputs is shown in Figure 11 with frequent changes between its lower bound and upper bound.

Figure 9.

Libration angle of ${\theta }_1$ and ${\theta }_2$ with variable climber speed.

Figure 10.

Length and its changing ratio of $L_1$, $L_c$ with variable climber speed.

Figure 11.

Control inputs for the closed-loop control with changing climber speed.

4. Conclusions

This chapter investigated a piecewise parallel onboard optimal control algorithm to solve the optimal control issues in complex nonlinear dynamic systems in aerospace engineering. To test the validity of the proposed two-phase optimal control scheme, the long-term tether libration stability and fast nano-satellite deorbit under complex environmental perturbations and the libration suppression for PSE system are considered. For EDT system, instead of optimizing the control of fast and stable nano-satellite deorbit over the complete process, the current approach divides the deorbit process into a set of intervals. For the PSE system, each time interval is set depends on the minimize transfer time for equal transfer length interval. Within each interval, the predicting phase simplifies significantly the optimal control problem. The dynamic equations of libration motion are further simplified to reduce computational loads using the simple dynamic models. The trajectory of the stable libration states and current control input is then optimized for the fast deorbit within the interval based on the simplified dynamic equations. The tracking optimizes the trajectory tracking control using the finite receding horizon control theory within the time interval corresponding to the open-loop control state trajectory with the same interval number. By applying the close-loop control modification, the system motions are integrated without any simplification of the dynamics or environmental perturbations and the instantaneous states of the orbital and libration motions. The i-th time interval’s closed-loop tracking is processed in tracking phase while the (i + 1)-th time interval’s optimal state trajectory is predicted in the predicting phase. This prediction is based on the error between the real state and the predicting state in the (i-1)-th time interval. By repeating the process, the optimal control problem can be achieved in a piecewise way with dramatically reduced computation effort. Compared with the current on–off control where the stable libration motion is the only control target, numerical results show that the proposed optimal control scheme works well in keeping the libration angles within an allowed limit.

Acknowledgments

This work is funded by the Discovery Grant of the Natural Sciences and Engineering Research Council of Canada, the National Natural Science Foundation of China, Grand No. 11472213 and the Chinese Scholarship Council Scholarship No. 201606290135.

A. Appendix

The detailed expressions for the matrixes C and M are shown as followed,

$$\mathit{C}=\left[{\mathit{C}}_1{\mathit{C}}_2{\mathit{C}}_3\right],\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathit{M}=\left[\begin{array}{lll}{\mathit{M}}_{11}& {\mathit{M}}_{12}& 0\\ {\mathit{M}}_{12}^T& {\mathit{M}}_{22}& 0\\ 0& 0& {\mathit{M}}_{33}\end{array}\right]$$

$$\begin{array}{l}{C}_1=\left[\begin{array}{lllll}E\hfill & \hfill & \hfill & \hfill & 0\hfill \\ {\chi }_0^1\hfill & {\chi }_0^2\hfill & \hfill & \hfill & \hfill \\ \hfill & {\chi }_1^1\hfill & {\chi }_1^2\hfill & \hfill & \hfill \\ \hfill & \hfill & \ddots \hfill & \ddots \hfill & \hfill \\ \hfill & \hfill & \hfill & {\chi }_{n-1}^1\hfill & {\chi }_{n-1}^2\hfill \\ 0\hfill & \hfill & \hfill & \hfill & E\hfill \end{array}\right],C_2=\left[\begin{array}{lllll}0\hfill & \hfill & \hfill & \hfill & 0\hfill \\ {\chi }_0^3\hfill & {\chi }_0^4\hfill & \hfill & \hfill & \hfill \\ \hfill & {\chi }_1^3\hfill & {\chi }_1^4\hfill & \hfill & \hfill \\ \hfill & \hfill & \ddots \hfill & \ddots \hfill & \hfill \\ \hfill & \hfill & \hfill & {\chi }_{n-1}^3\hfill & {\chi }_{n-1}^4\hfill \\ 0\hfill & \hfill & \hfill & \hfill & 0\hfill \end{array}\right],C_3=\frac{2}{3}\overline{\gamma }\left[\begin{array}{lll}0\hfill & \hfill & 0\hfill \\ B_{0.5}\hfill & \hfill & \hfill \\ \hfill & \ddots \hfill & \hfill \\ \hfill & \hfill & B_{n-0.5}\hfill \\ 0\hfill & \hfill & 0\hfill \end{array}\right]\hfill \\ {\mathit{\chi }}_j^1=\mathit{E}+\frac{\overline{\gamma }}{6}{\mathit{A}}_j+\frac{\overline{\gamma }}{3}{\mathit{A}}_{j+0.5}+\frac{{\overline{\gamma }}^2}{12}{\mathit{A}}_{j+0.5}{\mathit{A}}_j,{\mathit{\chi }}_j^2=-\mathit{E}+\frac{\overline{\gamma }}{6}{\mathit{A}}_{j+1}+\frac{\overline{\gamma }}{3}{\mathit{A}}_{j+0.5}-\frac{{\overline{\gamma }}^2}{12}{\mathit{A}}_{j+0.5}{\mathit{A}}_{j+1}\hfill \\ {\mathit{\chi }}_j^3=\frac{{\overline{\gamma }}^2}{12}{\mathit{A}}_{j+0.5}{\mathit{B}}_j+\frac{\overline{\gamma }}{6}{\mathit{B}}_j,{\mathit{\chi }}_j^4=-\frac{{\overline{\gamma }}^2}{12}{\mathit{A}}_{j+0.5}{\mathit{B}}_{j+1}+\frac{\overline{\gamma }}{6}{\mathit{B}}_{j+1}\hfill \end{array}$$

$$\begin{array}{l}{\mathit{M}}_{11}=\left[\begin{array}{lllll}{\varpi }_{00}^{11}& {\varpi }_{01}^{11}& & & 0\\ {\left({\varpi }_{01}^{11}\right)}^T& {\varpi }_{11}^{11}& & & \\ & & \ddots & & \\ & & & {\varpi }_{n-1n-1}^{11}& {\varpi }_{n-1n}^{11}\\ 0& & & {\left({\varpi }_{n-1n}^{11}\right)}^T& {\varpi }_{\mathit{nn}}^{11}\end{array}\right]{\mathit{M}}_{22}=\left[\begin{array}{lllll}{\varpi }_{00}^{22}& {\varpi }_{01}^{22}& & & 0\\ {\left({\varpi }_{01}^{22}\right)}^T& {\varpi }_{11}^{22}& & & \\ & & \ddots & & \\ & & & {\varpi }_{n-1n-1}^{22}& {\varpi }_{n-1n}^{22}\\ 0& & & {\left({\varpi }_{n-1n}^{22}\right)}^T& {\varpi }_{\mathit{nn}}^{22}\end{array}\right]\hfill \end{array}$$

$$\begin{array}{l}{\mathit{M}}_{12}=\left[\begin{array}{lllll}{\varpi }_{00}^{12}& {\varpi }_{01}^{12}& & & 0\\ {\varpi }_{10}^{12}& {\varpi }_{11}^{12}& & & \\ & & \ddots & & \\ & & & {\varpi }_{n-1n-1}^{12}& {\varpi }_{n-1n}^{12}\\ 0& & & {\varpi }_{\mathit{nn}-1}^{12}& {\varpi }_{\mathit{nn}}^{12}\end{array}\right]{\mathit{M}}_{33}=\frac{2}{3}\overline{\gamma }\left[\begin{array}{lll}R& & 0\\ & \ddots & \\ 0& & R\end{array}\right]\hfill \end{array}$$

$$\begin{array}{l}{\varpi }_{\mathit{jj}}^{11}=\frac{\overline{\gamma }}{6}\left(4\mathit{Q}+\frac{{\overline{\gamma }}^2}{8}{\mathit{A}}_j^T{\mathit{QA}}_j\right),{\varpi }_{\mathit{jj}+1}^{11}=\frac{\overline{\gamma }}{6}\left(\mathit{Q}-\frac{{\overline{\gamma }}^2}{16}{\mathit{A}}_j^T{\mathit{QA}}_{j+1}+\frac{\overline{\gamma }}{4}{\mathit{A}}_j^T\mathit{Q}-\frac{\overline{\gamma }}{4}{\mathit{QA}}_{j+1}\right)\hfill \\ {\varpi }_{00}^{11}=\frac{\overline{\gamma }}{6}\left(2\mathit{Q}+\frac{{\overline{\gamma }}^2}{16}{\mathit{A}}_0^T{\mathit{QA}}_0+\frac{\overline{\gamma }}{4}{\mathit{A}}_0^T\mathit{Q}+\frac{\overline{\gamma }}{4}{\mathit{QA}}_0\right),{\varpi }_{\mathit{nn}}^{11}=\frac{\overline{\gamma }}{6}\left(2\mathit{Q}+\frac{{\overline{\gamma }}^2}{16}{\mathit{A}}_n^T{\mathit{QA}}_n-\frac{\overline{\gamma }}{4}{\mathit{A}}_n^T\mathit{Q}-\frac{\overline{\gamma }}{4}{\mathit{QA}}_n\right)+{\mathit{S}}_f\hfill \\ {\varpi }_{\mathit{jj}}^{22}=\frac{\overline{\gamma }}{6}\left(2R+\frac{{\overline{\gamma }}^2}{8}{\mathit{B}}_j^T{\mathit{QB}}_j\right),{\varpi }_{\mathit{jj}+1}^{22}=-\overline{\gamma }{\mathit{B}}_j^T{\mathit{QB}}_{j+1}\hfill \\ {\varpi }_{00}^{22}=\frac{\overline{\gamma }}{6}\left(R+\frac{{\overline{\gamma }}^2}{16}{\mathit{B}}_0^T{\mathit{QB}}_0\right),{\varpi }_{\mathit{nn}}^{22}=\frac{\overline{\gamma }}{6}\left(R+\frac{{\overline{\gamma }}^2}{16}{\mathit{B}}_n^T{\mathit{QB}}_n\right)\hfill \\ {\varpi }_{\mathit{jj}}^{12}=\frac{{\overline{\gamma }}^3}{48}{\mathit{A}}_j^T{\mathit{QB}}_j,{\varpi }_{\mathit{jj}+1}^{12}=-\frac{{\overline{\gamma }}^2}{24}\left({\mathit{QB}}_{j+1}+\frac{\overline{\gamma }}{4}{\mathit{A}}_j^T{\mathit{QB}}_{j+1}\right)\hfill \\ {\varpi }_{j+1j}^{12}=\frac{{\overline{\gamma }}^2}{24}\left({\mathit{QB}}_j-\frac{\overline{\gamma }}{4}{\mathit{A}}_{j+1}^T{\mathit{QB}}_j\right),{\varpi }_{00}^{12}=\frac{{\overline{\gamma }}^2}{24}\left({\mathit{QB}}_0+\frac{\overline{\gamma }}{4}{\mathit{A}}_0^T{\mathit{QB}}_0\right),{\varpi }_{\mathit{nn}}^{12}=-\frac{{\overline{\gamma }}^2}{24}\left({\mathit{QB}}_{{\displaystyle \overline{N}}}-\frac{\overline{\gamma }}{4}{\mathit{A}}_n^T{\mathit{QB}}_n\right)\hfill \end{array}$$

where E is the unit matrix which has the same dimension as A_j, and j = 0,1,2,…,n-1.