Open access peer-reviewed chapter

Piecewise Parallel Optimal Algorithm

Written By

Zheng Hong Zhu and Gefei Shi

Submitted: October 14th, 2017 Reviewed: March 20th, 2018 Published: September 5th, 2018

DOI: 10.5772/intechopen.76625

Chapter metrics overview

901 Chapter Downloads

View Full Metrics


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.


  • 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 Siwill 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 S˜ii+1=S˜ii+12ei1, where S˜iidenotes 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. S˜ii+1is the estimated initial state of the (i + 1)-th interval, and ei1=Si1S˜i1i1is 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, S0=S˜01and e0=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 titi+1is discretized into nintervals, where tiand ti+1are the initial and final time, respectively. ti+1can 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


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 titi+1is discretized into nsubintervals with n + 1 nodes at the discretized time τk(k=0,1,.,n).


The state vectors and control inputs are discretized at n + 1 nodes, x0, x1, x2, …, xnand υ0, υ1, υ2, …, υn. Further, denote the state vectors and the control inputs at mid-points between adjacent nodes by x0.5, x1.5, x2.5, …, xn0.5and υ0.5, υ1.5, υ2.5, …, υn0.5and the mid-point state vectors, xk+0.5, can be derived by the Hermite interpolation scheme,


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, Ak=Aτk, Ak+0.5=Aτk+0.5, Bk=Bτk, Bk+0.5=Bτk+0.5.

The derivation of Eq. (8a) finally leads to a quadratic programming problem to find a programming vector Z=δx0Tδx1TδxnTδυ0δυ1δυnδυ0.5δυ1.5δυn0.5T, which minimizes the cost function:


subject to.

CZ=X    X=δx0T000TE12

where the matrices Cand Mare 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 Vis defined to “choose” the target value from the optimal solution, and the position of “1” in the row vector Vis the same as the position of δυ0in the column vector Z. Finally, the control input of the closed-loop control, υti, is


It is apparent that the closed-loop control law derived here is a linear proportional feedback control law, and the feedback gain matrix Kis a function of time. Without any explicit integration of differential equations, Kcan 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 Mand Care both formulated by the influence matrices at certain discretization notes (Ak, Bk, Ak + 0.5, Bk + 0.5). If tiand tiare 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 Mis unchanged if treating the terminal horizon time ti+thas the time-to-go and keep the future horizon interval titi+1unchanged. This means the inverse of Mcould 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 Oat 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, Ω, i, ex, ey, φ)


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 σxand σzare in the orbital plane, and σzis the radial component pointing outwards. The out-of-plane component σycompletes 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 eox, eoy, and eoz, respectively. Then, the instantaneous attitudes of the EDT system are described by an in-plane angle α(pitch angle, rotating about the y-axis) followed by an out-of-plane angle β(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,


where m˜=m1m2+m1+m2mt/3+mt2/12mEDT1is the equivalent mass, (Qα, Qβ) are the corresponding perturbative torques by the perturbative forces to be discussed below.

The perturbative accelerations (σx, σy, σz) and torques (Qα, Qβ) 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,


where the polar radius rpoand the Earth’s eccentricity eEare 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


where r0 = 6371.2 × 103 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 xtυtover each time interval titi+1to minimize a cost function of the negative work done by the electrodynamic force


subject to the nonlinear state equations of libration motion


where x1x2x3x4=ααββ, η=1+ecosν, ν̇=μ/p30.51+ecosν2, r=p1+ecosν1, λ=m1+0.5mt/mEDTis determined by the mass ratio between the end-bodies, and ςis determined by the distribution of current along the EDT, such that, ς=Iave1L20LsIsds. Accordingly, ς=0.5is used for the assumption of a constant current in the EDT. The initial conditions xti=xstartand the box constraint ααmax,ββmax,IminIaveImax. The environmental perturbations are simplified by considering only the electrodynamic force with a simple non-tilted dipole model of geomagnetic field, such that,


Accordingly, the electrodynamic force Feexerting on the EDT can be obtained as,


3.1.2. Results and discussion

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

Mass of the main satellite5 kg
Mass of subsatellite1.75 kg
Mass of the tether0.25 kg
Dimensions of main satellite0.2 × 0.2 × 0.2 m
Dimensions of subsatellite0.1 × 0.17 × 0.1 m
Tether length500 m
Tether diameter0.0005 m
Tether conductivity (aluminum)3.4014 × 107 Ω−1 m−1
Tether current lower/upper limits0 ~ 0.8 A for the equatorial orbit
Orbital altitudes700 ~ 800 km

Table 1.

Parameters of an EDT system.

Imax(equatorial)0.4 A
Imin(equatorial & inclined)0 A
αmax(equatorial & inclined)45 degrees
β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 L1and L2, 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, m1and m2) 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 OXYwith its origin at the centre of Earth. Denoting the position of the main satellite (M) by a vector rmeasuring from the centre of Earth. The climber m1is connected to the main satellite Mby a tether 1 with the length of L1and a libration angle θ1measured from the vector r. The distance between them is controlled by reeling in/out tether 1 at main satellite. The end body m2is connected to m1by a tether 2 with the length of L2and a libration angle θ2measured from the vector r. The length of tether 2 L2is 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.


where L0is the initial total length of two pieces of the tethers and Lcis the length increment relates to L0.

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


subject to the simplified dynamic equations


where u=T1T2, x=θ1θ̇1θ2θ̇2L1L̇1and idenotes the interval number. In (26) the gravitational perturbations and the trigonometric functions are ignored, then they can be simplified following the assumptions: sinθjθj,cosθj1j=12. 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 0T1T1max,0T2T2max,θ1θ1max,θ2θ2max,LcLcLimitL̇1mL̇1L̇1M,L̇cmL̇cL̇cM, where T1maxand T2maxare the upper bounds of the tension control inputs T1 and T2, respectively. θ1max, θ2maxand LcLimitare the magnitudes of libration angles θ1, θ2 and the maximum available length scale of Lc, respectively. L̇1mand L̇1Mare the lower and upper bounds of climber’s moving speed L̇1, respectively. L̇cmand L̇cMare the lower and upper bounds of end-bodies’ moving speed 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 titi+1evenly into nsubintervals. 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, m1 = 500 kg, m2 = 1000 kg, θ1(0) = θ2(0) = 0, L0 = 20 km, L1(0) = 19,500 m, Lc(0) = 0, θ̇10=θ̇20=0, L̇10=20m/s, and L̇c0=0for 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 T1max=T2max=200N, θ1max=θ2max=0.3rad, L̇1m=15m/s, L̇1M=25m/s, L̇cm=10m/sand L̇cM=10m/s.

The simulation results of this case are shown in Figures 911. The climber’s open-loop libration angle approaches its upper bound at 850 s. After 850 sθ1is kept at 0.3 radby the end of the ascending period, see the dashed line in Figure 10. Using the closed-loop control, the tracking trajectory of θ1matches 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 θ2.The trajectory of θ2obtained by closed-loop control tracks the open-loop trajectory well and reaches 0.1 radby the end of the ascending period. The closed-loop trajectories of L1and Lcare shown in Figure 10. They are the reflections of the control inputs. Both L1and Lcshow smooth fluctuations between 40s and 350 s. Figure 10 shows the time history of trajectories of L̇1and L̇c, respectively. In the first 140 s the trajectory of L̇cincreases continuously until reaches its upper bound. Then, it keeps at 10 m/sby 270 s with some slight fluctuations. From 270 s to 355 s, it reduces continuously to −3 m/s. After that, L̇cfluctuates around −3 m/s by the end of the transfer period. As a reflection of the control input, L̇1also 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 11with frequent changes between its lower bound and upper bound.

Figure 9.

Libration angle ofθ1andθ2with variable climber speed.

Figure 10.

Length and its changing ratio ofL1,Lcwith 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.



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 Cand Mare shown as followed,


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


  1. 1. Zhong R, Zhu ZH. Optimal control of Nanosatellite fast Deorbit using Electrodynamic tether. Journal of Guidance Control and Dynamics. 2014;37:1182-1194. DOI: 10.2514/1.62154
  2. 2. Jablonski AM, Scott R. Deorbiting of microsatellites in low earth orbit (LEO) – An introduction. Canadian Aeronautics and Space Journal. 2009;55:55-67
  3. 3. Stevens RE, Baker WP. Optimal control of a librating Electrodynamic tether performing a multirevolution orbit change. Journal of Guidance Control and Dynamics. 2009;32:1497-1507. DOI: 10.2514/1.42679
  4. 4. Williams P. Optimal control of Electrodynamic tether orbit transfers using timescale separation. Journal of Guidance Control and Dynamics. 2010;33:88-98. DOI: 10.2514/1.45250
  5. 5. Aslanov V, Ledkov A. Dynamics of Tethered Satellite Systems. Cambridge, UK: Elsevier; 2012
  6. 6. Wen H, Zhu ZH, Jin DP, Hu HY. Space tether deployment control with explicit tension constraint and saturation function. Journal of Guidance, Control, and Dynamics. 2016;39:916-921. DOI: 10.2514/1.G001356
  7. 7. Ma ZQ, Sun GH, Li ZK. Dynamic adaptive saturated sliding mode control for deployment of tethered satellite system. Aerospace Science and Technology. 2017;66:355-365. DOI: 10.1016/j.ast.2017.01.030
  8. 8. Jin DP, Hu HY. Optimal control of a tethered subsatellite of three degrees of freedom. Nonlinear Dynamics. 2006;46:161-178. DOI: 10.1007/s11071-006-9021-4
  9. 9. Williams P. Deployment/retrieval optimization for flexible tethered satellite systems. Nonlinear Dynamics. 2008;52:159-179. DOI: 10.1007%2Fs11071-007-9269-3
  10. 10. Williams P, Ockels W. Climber motion optimization for the tethered space elevator. Acta Astronautica. 2010;66:1458-1467. DOI: 10.1016/j.actaastro.2009.11.003
  11. 11. Kojima H, Fukatsu K, Trivailo PM. Mission-function control of tethered satellite/climber system. Acta Astronautica. 2015;106:24-32. DOI: 10.1016/j.actaastro.2014.10.024
  12. 12. Williams P. Optimal control of Electrodynamic tether orbit transfers using timescale separation. Journal of Guidance Control and Dynamics. 2010;33:88-98. DOI: 10.2514/1.45250
  13. 13. Betts JT. Practical methods for optimal control using nonlinear programming. Society for Industrial and Applied Mathematics: Advances in Control and Design Series, Philadelphia, PA. 2001:65-72
  14. 14. Zhong R, Xu S. Orbit-transfer control for TSS using direct collocation method. Acta Aeronautica et Astronautica Sinica. 2010;31:572-578
  15. 15. Bryson AE, Ho YC. Chapter 2. In: Applied Optimal Control. New York: Hemisphere; 1975
  16. 16. Herman AL, Conway BA. Direct optimization using collocation based on high-order gauss-Lobatto Quadrature rules. Journal of Guidance Control and Dynamics. 1996;19:592-599. DOI: 10.2514/3.21662
  17. 17. Hargraves CR, Paris SW. Direct trajectory optimization using nonlinear programming and collocation. Journal of Guidance, Control, and Dynamics. 1987;10:338-342. DOI: 10.2514/3.20223
  18. 18. Vallado DA. Fundamentals of Astrodynamics and Applications. New York, Springer: Microcosm Press; 2007
  19. 19. NASA, Earth Fact Sheet,, [retrieved March 16, 2012]
  20. 20. Davis J. Technical note. In: Mathematical Modeling of Earth's Magnetic Field. Blacksburg: Virginia Tech; 2004
  21. 21. Zhong R, Zhu ZH. Long-term libration dynamics and stability analysis of Electrodynamic tethers in spacecraft Deorbit. Journal of Aerospace Engineering. 2012. DOI: 10.1061/(ASCE)AS.1943-5525.0000310
  22. 22. Shi G, Zhu Z, Zhu ZH. Libration suppression of tethered space system with a moving climber in circular orbit. Nonlinear Dynamics. 2017:1-15. DOI: 10.1007/s11071-017-3919-x
  23. 23. Wen H, Zhu ZH, Jin DP, Hu H. Constrained tension control of a tethered space-tug system with only length measurement. Acta Astronautica. 2016;119:110-117. DOI:
  24. 24. Fletcher R. Chapter 10. In: Practical Methods of Optimization. 2nd ed. New York: Wiley; 1989

Written By

Zheng Hong Zhu and Gefei Shi

Submitted: October 14th, 2017 Reviewed: March 20th, 2018 Published: September 5th, 2018