Open access peer-reviewed chapter

Piecewise Parallel Optimal Algorithm

By Zheng Hong Zhu and Gefei Shi

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

DOI: 10.5772/intechopen.76625

Downloaded: 217

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 Si 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 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 n intervals, 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

Ji=titi+1ΠxudtE1

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

γk=τk+1τk,k=1nγk=ti+1tiE2

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,

xk+0.5=12xk+xk+1+γ8ΓxkυkτkΓ(xk+1υk+1τk+1)E3

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

Jγ6k=0n1Πxkυkτk+4Π(xk+0.5υk+0.5τk+0.5)+Π(xk+1υk+1τk+1)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

γ6Πxkυktk+4Π(xk+0.5υk+0.5τk+0.5)+Π(xk+1υk+1τk+1)+xkxk+1=0E5

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,

x0=xstart,xminxkxmax,υminυkυmax,υminυk+0.5υmaxE6

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

G=12δxi+1TSδxi+1+γ12k=0n1δxkTQδxk+4δxk+0.5TQδxk+0.5+δxk+1TQδxk+1+Rδυk2+4δυk+0.52+δυk+12E7

and the constraints are discretized into

δxkδxk+1+γ6Akδxk+Bkδυk+4Ak+0.5δxk+0.5+4Bk+0.5δυk+0.5+Ak+1δxk+1+Bk+1δυk+1=0E8
δxk=xτkxoptτkE9
δxk+0.5=12δxk+δxk+1+γ8Akδxk+BkδυkAk+1δxk+1Bk+1δυk+1E10

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:

G=12ZTMZE11

subject to.

CZ=X    X=δx0T000TE12

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

Z=M1CTCM1CT1XE13

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

δυti=VZ=VM1CTCM1CT1XKtinthxtixopttiE14

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

υti=υoptti+δυti=υoptti+KtinthxtixopttiE15

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 M and C are both formulated by the influence matrices at certain discretization notes (Ak, Bk, Ak + 0.5, Bk + 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 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 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, Ω, i, ex, ey, φ)

dadt=2n1e2σzesinν+σxprE16
dΩdt=σyrsinuna21e2siniE17
didt=σyrcosuna21e2E18
dexdt=1e2naσzsinu+σx1+rpcosu+rpex+dΩdteycosiE19
deydt=1e2naσzcosu+σx1+rpsinu+rpeydΩdtexcosiE20
dt=n1naσz2ra+1e21+1e2ecosνσx1+rp1e21+1e2esinνσyrcosisinuna21e2siniE21

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,

α¨+ν¨2α̇+ν̇β̇tanβ+3μr3sinαcosα=Qαm˜L2cos2βE22
β¨+α̇+ν̇2sinβcosβ+3μr3cos2αsinβcosβ=Qβm˜L2E23

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,

hg=rrpo1eE2cos2θ1/2E24

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

Bϕ=1sinθn=1r0rn+2m=0nmgnmsinhnmcosPnmθcBθ=n=1r0rn+2m=0ngnmcos+hnmsinPnmθcθcBr=n=1r0rn+2n+1m=0ngnmcos+hnmsinPnmθcE25

where r0 = 6371.2 × 103 km is the reference radius of the Earth, respectively.

The average current in the EDT is defined as

Iave=1L0LIsdsE26

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

J=titi+1Fe·vdtttti+1ΠxυtdtE27

subject to the nonlinear state equations of libration motion

x1=x2x2=2η3esinν+2x2+η2x4tanx33η3sinx1cosx1+sinitanx32sinucosx1cosusinx1cosiςλIaveμmμm˜η3E28
x3=x4x4=x2+η22sinx3cosx33η3cos2x1sinx3cosx3sini2sinusinx1+cosucosx1ςλIaveμmμm˜η3E29

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,

B=μmr3cosusinieox+μmr3cosieoy2μmr3sinusinieozE30

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

Fe=0LB×Ilds=IaveLB×lE31

3.1.2. Results and discussion

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

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

ParametersValues
Imax(equatorial)0.4 A
Imax(inclined)0.1sinisinβ¯BsinΩG+α¯BΩ+cosicosβ¯Bcos1iβ¯BA
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 OXY with 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.

θ¨1=3ω2sin2θ122ω+θ̇1L̇1L1sinθ1θ2T2L1m1E32
θ¨2=3ω2sin2θ222ω+θ̇2L̇cL̇1L0L1+Lc+sinθ1θ2T1L0L1+Lcm1E33
L¨1=3ω2L1cos2θ1+2ωL1θ̇1+L1θ̇12T1m1+cosθ1θ2T2m1E34
L¨c=3ω2L0L1+Lccos2θ2+3ω2L1cos2θ1+2ω+θ̇2L0L1+Lcθ̇2+2ω+θ̇1L1θ̇1+cosθ1θ21T1m1m1m2cosθ1θ2+m2T2m1m2E35

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

Ji=tfE36

subject to the simplified dynamic equations

θ¨1=3ω2θ12ω+θ̇1L̇1L1θ1θ2T2L1m1E37
θ¨2=3ω2θ22ω+θ̇2L̇cL̇1L0L1+Lc+θ1θ2T1L0L1+Lcm1E38
L¨1=3ω2L1+2ωL1θ̇1+L1θ̇12+T2T1m1E39
L¨c=3ω2L0+Lc+2ω+θ̇2L0L1+Lcθ̇2+2ω+θ̇1L1θ̇1T2m2E40

where u=T1T2, x=θ1θ̇1θ2θ̇2L1L̇1and 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θ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 T1max and T2max are the upper bounds of the tension control inputs T1 and T2, respectively. θ1max, θ2max and LcLimit are 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 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, 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 rad by 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 rad by 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/s by 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 11 with frequent changes between its lower bound and upper bound.

Figure 9.

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

Figure 10.

Length and its changing ratio of L1, Lc 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,

C=C1C2C3,M=M11M120M12TM22000M33
C1=[E0χ01χ02χ11χ12χn11χn120E],C2=[00χ03χ04χ13χ14χn13χn1400],C3=23γ¯[00B0.5Bn0.500]χj1=E+γ¯6Aj+γ¯3Aj+0.5+γ¯212Aj+0.5Aj,χj2=E+γ¯6Aj+1+γ¯3Aj+0.5γ¯212Aj+0.5Aj+1χj3=γ¯212Aj+0.5Bj+γ¯6Bj,χj4=γ¯212Aj+0.5Bj+1+γ¯6Bj+1
M11=ϖ0011ϖ01110ϖ0111Tϖ1111ϖn1n111ϖn1n110ϖn1n11Tϖnn11M22=ϖ0022ϖ01220ϖ0122Tϖ1122ϖn1n122ϖn1n220ϖn1n22Tϖnn22
M12=ϖ0012ϖ01120ϖ1012ϖ1112ϖn1n112ϖn1n120ϖnn112ϖnn12M33=23γ¯R00R
ϖjj11=γ¯64Q+γ¯28AjTQAj,ϖjj+111=γ¯6Qγ¯216AjTQAj+1+γ¯4AjTQγ¯4QAj+1ϖ0011=γ¯62Q+γ¯216A0TQA0+γ¯4A0TQ+γ¯4QA0,ϖnn11=γ¯62Q+γ¯216AnTQAnγ¯4AnTQγ¯4QAn+Sfϖjj22=γ¯62R+γ¯28BjTQBj,ϖjj+122=γ¯BjTQBj+1ϖ0022=γ¯6R+γ¯216B0TQB0,ϖnn22=γ¯6R+γ¯216BnTQBnϖjj12=γ¯348AjTQBj,ϖjj+112=γ¯224QBj+1+γ¯4AjTQBj+1ϖj+1j12=γ¯224QBjγ¯4Aj+1TQBj,ϖ0012=γ¯224QB0+γ¯4A0TQB0,ϖnn12=γ¯224QBN¯γ¯4AnTQBn

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

© 2018 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

How to cite and reference

Link to this chapter Copy to clipboard

Cite this chapter Copy to clipboard

Zheng Hong Zhu and Gefei Shi (September 5th 2018). Piecewise Parallel Optimal Algorithm, Optimization Algorithms - Examples, Jan Valdman, IntechOpen, DOI: 10.5772/intechopen.76625. Available from:

chapter statistics

217total chapter downloads

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

Related Content

This Book

Next chapter

Bilevel Disjunctive Optimization on Affine Manifolds

By Constantin Udriste, Henri Bonnel, Ionel Tevy and Ali Sapeeh Rasheed

Related Book

First chapter

Digital Image Processing with MATLAB

By Mahmut Sinecen

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More about us