Parameters of an EDT system.

## 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.

### Keywords

- optimal control
- parallel onboard optimal algorithm
- discretizing Hermite–Simpson method
- nonlinear dynamic system
- aerospace engineering

## 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 *i*-th interval obtained by CUP-1, the superscript and subscript denote the internal number and the states node number. *i* + 1)-th interval, and *i*-1)-th interval, the data used to calculate the error is picked up from the memory. For the first interval, *i* = 1,

### 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 *n* intervals, where

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 *n* subintervals with *n +* 1 nodes at the discretized time

The state vectors and control inputs are discretized at *n* + 1 nodes,

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

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

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,

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

and the constraints are discretized into

where,

The derivation of Eq. (8a) finally leads to a quadratic programming problem to find a programming vector

subject to.

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].

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

where the row vector

It is apparent that the closed-loop control law derived here is a linear proportional feedback control law, and the feedback gain matrix

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 *ti* and *ti* 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

## 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 (

The components of perturbative accelerations are defined in a local frame. The components

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 *y*-axis) followed by an out-of-plane angle *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,

where

The perturbative accelerations (

where the polar radius

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

where *r*_{0} = 6371.2 × 10^{3} km is the reference radius of the Earth, respectively.

The average current in the EDT is defined as

The open-loop optimal control problem for EDT deorbit can be stated as finding a state-control pair

subject to the nonlinear state equations of libration motion

where

Accordingly, the electrodynamic force

#### 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 |

Parameters | Values |
---|---|

0.4 A | |

0 A | |

45 degrees | |

45 degrees |

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.

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.

### 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 *OXY* with its origin at the centre of Earth. Denoting the position of the main satellite (

where

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

subject to the simplified dynamic equations

where *i* denotes the interval number. In (26) the gravitational perturbations and the trigonometric functions are ignored, then they can be simplified following the assumptions:

To ensure the availability and the suppression of the libration angles, following constrains are also required to be subjected *T*_{1max} and *T*_{2max} are the upper bounds of the tension control inputs *T*_{1} and *T*_{2}, respectively. *θ*_{1max}, *θ*_{2max} and *LcLimit* are the magnitudes of libration angles *θ*_{1}, *θ*_{2} and the maximum available length scale of *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_{1} = 500 *kg*, *m*_{2} = 1000 *kg*, *θ*_{1}(0) = *θ*_{2}(0) = 0, *L*_{0} = 20 *km*, *L*_{1}(0) = 19,500 *m*, *Lc* (0) = 0,

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**rad* by the end of the ascending period, see the dashed line in Figure 10. Using the closed-loop control, the tracking trajectory of *rad* by the end of the ascending period. The closed-loop trajectories of *m/s* by 270 s with some slight fluctuations. From 270 s to 355 s, it reduces continuously to −3 m/s. After that, **Figure 11** with frequent changes between its lower bound and upper bound.

## 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,

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