2.1. Control scheme

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/nanosatellites 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

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

elevator (PSE) system are considered and tested.

94 Optimization Algorithms - Examples

Assume two CUPs are used to process the calculation. CUP-1 is used to determine the openloop 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~<sup>i</sup>þ<sup>1</sup> <sup>i</sup> <sup>¼</sup> <sup>S</sup>~<sup>i</sup> <sup>i</sup> <sup>þ</sup> <sup>1</sup> <sup>2</sup> <sup>e</sup><sup>i</sup>�1, where <sup>S</sup>~<sup>i</sup> <sup>i</sup> denotes the estimated finial state of the i-th interval obtained by CUP-1, the superscript and subscript denote the internal number and the states node number. S~<sup>i</sup>þ<sup>1</sup> <sup>i</sup> is the estimated initial state of the (i + 1)-th interval, and <sup>e</sup><sup>i</sup>�<sup>1</sup> <sup>¼</sup> <sup>S</sup><sup>i</sup>�<sup>1</sup> � <sup>S</sup>~<sup>i</sup>�<sup>1</sup> <sup>i</sup>�<sup>1</sup> is the error between the real and the estimated final state of the (i-1)-th interval, the data used to calculate the error is picked up from the memory. For the first interval, <sup>i</sup> = 1, <sup>S</sup><sup>0</sup> <sup>¼</sup> <sup>S</sup>~<sup>1</sup> <sup>0</sup> and e<sup>0</sup> ¼ 0. The computational diagram of the entire control strategy is given as shown in Figure 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 ti ½ � ; tiþ<sup>1</sup> is discretized into n intervals, where ti and tiþ<sup>1</sup> are the initial and final time, respectively. tiþ<sup>1</sup> can be obtained in terms of the transfer length of the climber and its speed. To make calculation convenient to be realized in practical condition, the transfer process is divided evenly. The optimal trajectory should be found to satisfy the desired cost functions for each time intervals as

$$J\_i = \int\_{t\_i}^{t\_{i+1}} \Pi(\mathfrak{x}, \mathfrak{u}) \, dt \tag{1}$$

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

Assume that the time interval ti ½ � ; tiþ<sup>1</sup> is discretized into n subintervals with n + 1 nodes at the

Xn k¼1

The state vectors and control inputs are discretized at n + 1 nodes, x0, x1, x2, …, x<sup>n</sup> and υ0, υ1, υ2, …, υn. Further, denote the state vectors and the control inputs at mid-points between adjacent nodes by x<sup>0</sup>:5, x<sup>1</sup>:5, x<sup>2</sup>:5, …, x<sup>n</sup>�0:<sup>5</sup> and υ<sup>0</sup>:5, υ<sup>1</sup>:5, υ<sup>2</sup>:5, …, υ<sup>n</sup>�0:<sup>5</sup> and the mid-point state

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

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

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

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

<sup>6</sup> <sup>Π</sup> <sup>x</sup>k; <sup>υ</sup><sup>k</sup> <sup>½</sup> ð Þþ ; tk <sup>4</sup>Πðx<sup>k</sup>þ0:<sup>5</sup>; <sup>υ</sup><sup>k</sup>þ0:<sup>5</sup>; <sup>τ</sup><sup>k</sup>þ0:<sup>5</sup>Þ þ <sup>Π</sup>ð Þ <sup>x</sup><sup>k</sup>þ<sup>1</sup>; <sup>υ</sup><sup>k</sup>þ<sup>1</sup>; <sup>τ</sup><sup>k</sup>þ<sup>1</sup> � þ <sup>x</sup><sup>k</sup> � <sup>x</sup><sup>k</sup>þ<sup>1</sup> <sup>¼</sup> <sup>0</sup> (5)

x<sup>0</sup> ¼ xstart, xmin ≤ x<sup>k</sup> ≤ xmax, υmin ≤ υ<sup>k</sup> ≤ υmax, υmin ≤ υ<sup>k</sup>þ0:<sup>5</sup> ≤ υmax (6)

be denoted by discretized equations using the Simpson integration formula, such that

initial states and the continuous box constraints with the discretized constraints,

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

Π xk; υ<sup>k</sup> ½ � ð Þþ ; τ<sup>k</sup> 4Πðx<sup>k</sup>þ0:<sup>5</sup>; υ<sup>k</sup>þ0:<sup>5</sup>; τ<sup>k</sup>þ0:<sup>5</sup>Þ þ Πð Þ x<sup>k</sup>þ<sup>1</sup>; υ<sup>k</sup>þ<sup>1</sup>; τ<sup>k</sup>þ<sup>1</sup> (4)

γ 8

γ<sup>k</sup> ¼ tiþ<sup>1</sup> � ti (2)

Piecewise Parallel Optimal Algorithm http://dx.doi.org/10.5772/intechopen.76625 97

Γ xk; υ<sup>k</sup> ½ � ð Þ� ; τ<sup>k</sup> Γð Þ x<sup>k</sup>þ<sup>1</sup>; υ<sup>k</sup>þ<sup>1</sup>; τ<sup>k</sup>þ<sup>1</sup> (3)

scheme [14, 17], is adopted because of its simplicity and accuracy.

vectors, x<sup>k</sup>þ0:5, can be derived by the Hermite interpolation scheme,

ð Þþ x<sup>k</sup> þ x<sup>k</sup>þ<sup>1</sup>

<sup>x</sup><sup>k</sup>þ0:<sup>5</sup> <sup>¼</sup> <sup>1</sup> 2 γ<sup>k</sup> ¼ τ<sup>k</sup>þ<sup>1</sup> � τk,

discretized time τ<sup>k</sup> (k ¼ 0, 1, :…, n).

<sup>J</sup> ffi <sup>γ</sup> 6 Xn�1 k¼0

γ

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

Figure 1. Control scheme.

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.

The calculation of the optimal trajectory in CPU-1 is a prediction state trajectory whose initial

interval obtained by CUP-1, the superscript and subscript denote the internal number and the

interval, the data used to calculate the error is picked up from the memory. For the first

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 ti ½ � ; tiþ<sup>1</sup> is discretized into n intervals, where ti and tiþ<sup>1</sup> are the initial and final time, respectively. tiþ<sup>1</sup> can be obtained in terms of the transfer length of the climber and its speed. To make calculation convenient to be realized in practical condition, the transfer process is divided evenly. The optimal trajectory should be found to satisfy the desired cost functions

Ji ¼

ðtiþ<sup>1</sup> ti

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

<sup>i</sup> denotes the estimated finial state of the i-th

Πð Þ x; u dt (1)

<sup>i</sup> is the estimated initial state of the (i + 1)-th interval, and

<sup>0</sup> and e<sup>0</sup> ¼ 0. The computational diagram of the entire control strategy is

<sup>i</sup>�<sup>1</sup> is the error between the real and the estimated final state of the (i-1)-th

<sup>2</sup> <sup>e</sup><sup>i</sup>�1, where <sup>S</sup>~<sup>i</sup>

state is estimated as S~<sup>i</sup>þ<sup>1</sup>

96 Optimization Algorithms - Examples

states node number. S~<sup>i</sup>þ<sup>1</sup>

given as shown in Figure 1.

for each time intervals as

Figure 1. Control scheme.

2.2. Open-loop control trajectory

<sup>e</sup><sup>i</sup>�<sup>1</sup> <sup>¼</sup> <sup>S</sup><sup>i</sup>�<sup>1</sup> � <sup>S</sup>~<sup>i</sup>�<sup>1</sup>

interval, <sup>i</sup> = 1, <sup>S</sup><sup>0</sup> <sup>¼</sup> <sup>S</sup>~<sup>1</sup>

<sup>i</sup> <sup>¼</sup> <sup>S</sup>~<sup>i</sup>

<sup>i</sup> <sup>þ</sup> <sup>1</sup>

Assume that the time interval ti ½ � ; tiþ<sup>1</sup> is discretized into n subintervals with n + 1 nodes at the discretized time τ<sup>k</sup> (k ¼ 0, 1, :…, n).

$$
\gamma\_k = \tau\_{k+1} - \tau\_{k\prime} \sum\_{k=1}^n \gamma\_k = t\_{i+1} - t\_i \tag{2}
$$

The state vectors and control inputs are discretized at n + 1 nodes, x0, x1, x2, …, x<sup>n</sup> and υ0, υ1, υ2, …, υn. Further, denote the state vectors and the control inputs at mid-points between adjacent nodes by x<sup>0</sup>:5, x<sup>1</sup>:5, x<sup>2</sup>:5, …, x<sup>n</sup>�0:<sup>5</sup> and υ<sup>0</sup>:5, υ<sup>1</sup>:5, υ<sup>2</sup>:5, …, υ<sup>n</sup>�0:<sup>5</sup> and the mid-point state vectors, x<sup>k</sup>þ0:5, can be derived by the Hermite interpolation scheme,

$$\mathbf{x}\_{k+0.5} = \frac{1}{2}(\mathbf{x}\_k + \mathbf{x}\_{k+1}) + \frac{\mathcal{V}}{8}[\Gamma(\mathbf{x}\_k, \boldsymbol{\nu}\_k, \boldsymbol{\tau}\_k) - \Gamma(\mathbf{x}\_{k+1}, \boldsymbol{\nu}\_{k+1}, \boldsymbol{\tau}\_{k+1})] \tag{3}$$

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

$$J \cong \frac{\mathcal{V}}{6} \sum\_{k=0}^{n-1} \left[ \Pi(\mathbf{x}\_k, \boldsymbol{\nu}\_k, \boldsymbol{\tau}\_k) + 4 \Pi(\mathbf{x}\_{k+0.5}, \boldsymbol{\nu}\_{k+0.5}, \boldsymbol{\tau}\_{k+0.5}) + \Pi(\mathbf{x}\_{k+1}, \boldsymbol{\nu}\_{k+1}, \boldsymbol{\tau}\_{k+1}) \right] \tag{4}$$

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

$$\frac{\mathcal{V}}{6} \left[ \Pi(\mathbf{x}\_k, \boldsymbol{\nu}\_k, t\_k) + 4 \Pi(\mathbf{x}\_{k+0.5}, \boldsymbol{\nu}\_{k+0.5}, \boldsymbol{\tau}\_{k+0.5}) + \Pi(\mathbf{x}\_{k+1}, \boldsymbol{\nu}\_{k+1}, \boldsymbol{\tau}\_{k+1}) \right] + \mathbf{x}\_k - \mathbf{x}\_{k+1} = \mathbf{0} \tag{5}$$

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,

$$\mathbf{x}\_0 = \mathbf{x}\_{\text{start}} \\ \mathbf{x}\_{\text{min}} \le \mathbf{x}\_k \le \mathbf{x}\_{\text{max}} \\ \boldsymbol{\nu}\_{\text{min}} \le \boldsymbol{\nu}\_k \le \boldsymbol{\nu}\_{\text{max}} \\ \boldsymbol{\nu}\_{\text{min}} \le \boldsymbol{\nu}\_{k+0.5} \le \boldsymbol{\nu}\_{\text{max}} \tag{6}$$

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

$$\mathbf{G} = \frac{1}{2} \delta \mathbf{x}\_{i+1}^T \mathbf{S} \delta \mathbf{x}\_{i+1}$$

$$+ \frac{\mathcal{V}}{12} \sum\_{k=0}^{n-1} \left[ \delta \mathbf{x}\_k^T \mathbf{Q} \delta \mathbf{x}\_k + 4 \delta \mathbf{x}\_{k+0.5}^T \mathbf{Q} \delta \mathbf{x}\_{k+0.5} + \delta \mathbf{x}\_{k+1}^T \mathbf{Q} \delta \mathbf{x}\_{k+1} + \mathbf{R} \left( \delta \mathbf{u}\_k^2 + 4 \delta \mathbf{u}\_{k+0.5}^2 + \delta \mathbf{u}\_{k+1}^2 \right) \right] \tag{7}$$

and the constraints are discretized into

$$\delta \mathbf{x}\_k - \delta \mathbf{x}\_{k+1} + \frac{\mathcal{V}}{6} [\mathbf{A}\_k \delta \mathbf{x}\_k + \mathbf{B}\_k \delta \nu\_k + \mathbf{4} \mathbf{A}\_{k+0.5} \delta \mathbf{x}\_{k+0.5} + \mathbf{4} \mathbf{B}\_{k+0.5} \delta \nu\_{k+0.5} + \mathbf{A}\_{k+1} \delta \mathbf{x}\_{k+1} + \mathbf{B}\_{k+1} \delta \nu\_{k+1}] = \mathbf{0} \tag{8}$$

$$
\delta \mathfrak{x}\_k = \mathfrak{x}(\tau\_k) - \mathfrak{x}\_{opt}(\tau\_k) \tag{9}
$$

υð Þ¼ ti υoptð Þþ ti δυð Þ¼ ti υoptð Þþ ti K ti ð Þ ; n; th xð Þ� ti xoptð Þti

It is apparent that the closed-loop control law derived here is a linear proportional feedback control law, and the feedback gain matrix K is a function of time. Without any explicit integration of differential equations, K can either be determined offline or online depending

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 M is unchanged if treating the terminal horizon time ti þ th as the time-to-go and keep the future horizon interval ti ½ � ; tiþ<sup>1</sup> unchanged. This means the inverse of M could be calculated only once in the same interval. It is attractive for the online implementation of HRC, where the computational effort is critical. As the entire interval can be discretized into small intervals, if these intervals are sufficiently small relative to the computing power of the satellite, then the calculation process can be carried by the onboard computer. Furthermore, for small intervals, the computation for the open-loop optimal trajectory of the next interval can be done by CPU-1 while the tracking is still in process, see Figure 1. This makes of the proposed optimal suppression control a parallel online implementation, which is another advantage of this control scheme.

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

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

<sup>1</sup> � <sup>e</sup><sup>2</sup> <sup>p</sup> <sup>σ</sup>ze sin <sup>ν</sup> <sup>þ</sup> <sup>σ</sup><sup>x</sup>

dt <sup>¼</sup> <sup>σ</sup>yr sin <sup>u</sup> na<sup>2</sup> ffiffiffiffiffiffiffiffiffiffiffiffi

p r

� � (16)

<sup>1</sup> � <sup>e</sup><sup>2</sup> <sup>p</sup> sin <sup>i</sup> (17)

on the computation and restoration ability onboard the satellite.

3.1. Parallel optimal algorithm in attitudes control of EDT system

equations of six independent orbital elements (a, Ω, i, ex, ey, φ)

da dt <sup>¼</sup> <sup>2</sup> n ffiffiffiffiffiffiffiffiffiffiffiffi

dΩ

3. Cases study

with some existing control methods.

3.1.1. Problem formulation

� � (15)

Piecewise Parallel Optimal Algorithm http://dx.doi.org/10.5772/intechopen.76625 99

$$\delta \mathbf{x}\_{k+0.5} = \frac{1}{2} (\delta \mathbf{x}\_k + \delta \mathbf{x}\_{k+1}) + \frac{\gamma}{8} (\mathbf{A}\_k \delta \mathbf{x}\_k + \mathbf{B}\_k \delta \upsilon\_k - \mathbf{A}\_{k+1} \delta \mathbf{x}\_{k+1} - \mathbf{B}\_{k+1} \delta \upsilon\_{k+1}) \tag{10}$$

where, A<sup>k</sup> ¼ Að Þ τ<sup>k</sup> , A<sup>k</sup>þ0:<sup>5</sup> ¼ Að Þ τ<sup>k</sup>þ0:<sup>5</sup> , B<sup>k</sup> ¼ Bð Þ τ<sup>k</sup> , B<sup>k</sup>þ0:<sup>5</sup> ¼ Bð Þ τ<sup>k</sup>þ0:<sup>5</sup> .

The derivation of Eq. (8a) finally leads to a quadratic programming problem to find a programming vector <sup>Z</sup> <sup>¼</sup> <sup>δ</sup>x<sup>T</sup> <sup>0</sup> δx<sup>T</sup> <sup>1</sup> … δx<sup>T</sup> <sup>n</sup> δυ<sup>0</sup> δυ<sup>1</sup> … δυ<sup>n</sup> δυ<sup>0</sup>:<sup>5</sup> δυ<sup>1</sup>:<sup>5</sup> … δυ<sup>n</sup>�0:<sup>5</sup> � �<sup>T</sup> , which minimizes the cost function:

$$\mathbf{G} = \frac{1}{2} \mathbf{Z}^T \mathbf{M} \mathbf{Z} \tag{11}$$

subject to.

$$\mathbf{CZ} = \mathbf{X} \quad \mathbf{X} = \begin{bmatrix} \delta \mathbf{x}\_0^T \ \mathbf{0} \ \mathbf{0} \ \dots \ \mathbf{0} \end{bmatrix}^T \tag{12}$$

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

$$\mathbf{Z}^\* = \mathbf{M}^{-1} \mathbf{C}^T \left(\mathbf{C} \mathbf{M}^{-1} \mathbf{C}^T\right)^{-1} \mathbf{X} \tag{13}$$

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

$$\delta\boldsymbol{\nu}(t\_i) = \mathbf{V}\mathbf{Z}^\* = \mathbf{V}\mathbf{M}^{-1}\mathbf{C}^T \left(\mathbf{C}\mathbf{M}^{-1}\mathbf{C}^T\right)^{-1} \mathbf{X} \triangleq \mathbf{K}(t\_i, \boldsymbol{\nu}, t\_h) \left[\mathbf{x}(t\_i) - \mathbf{x}\_{\text{opt}}(t\_i)\right] \tag{14}$$

where the row vector V is defined to "choose" the target value from the optimal solution, and the position of "1" in the row vector V is the same as the position of δυ<sup>0</sup> in the column vector Z. Finally, the control input of the closed-loop control, υð Þti , is

$$\boldsymbol{\nu}(t\_i) = \boldsymbol{\nu}\_{opt}(t\_i) + \delta\boldsymbol{\nu}(t\_i) = \boldsymbol{\nu}\_{opt}(t\_i) + \mathbf{K}(t\_i, \boldsymbol{n}, t\_h) \left[ \mathbf{x}(t\_i) - \mathbf{x}\_{opt}(t\_i) \right] \tag{15}$$

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

It is worth to point out some advantages of this approach. Firstly, the matrices M and C are both formulated by the influence matrices at certain discretization notes (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 M is unchanged if treating the terminal horizon time ti þ th as the time-to-go and keep the future horizon interval ti ½ � ; tiþ<sup>1</sup> unchanged. This means the inverse of M could be calculated only once in the same interval. It is attractive for the online implementation of HRC, where the computational effort is critical. As the entire interval can be discretized into small intervals, if these intervals are sufficiently small relative to the computing power of the satellite, then the calculation process can be carried by the onboard computer. Furthermore, for small intervals, the computation for the open-loop optimal trajectory of the next interval can be done by CPU-1 while the tracking is still in process, see Figure 1. This makes of the proposed optimal suppression control a parallel online implementation, which is another advantage of this control scheme.

## 3. Cases study

optimization problem, the same direct collocation based on the Hermite-Simpson method is

By using the similar notations of discretization above, the cost function is discretized using the

<sup>i</sup>þ<sup>1</sup>Sδx<sup>i</sup>þ<sup>1</sup>

� � � �

<sup>6</sup> <sup>½</sup>AkδxkþBkδυ<sup>k</sup> <sup>þ</sup> <sup>4</sup>A<sup>k</sup>þ0:<sup>5</sup>δx<sup>k</sup>þ0:<sup>5</sup> <sup>þ</sup> <sup>4</sup>B<sup>k</sup>þ0:<sup>5</sup>δυ<sup>k</sup>þ0:<sup>5</sup> <sup>þ</sup> <sup>A</sup><sup>k</sup>þ<sup>1</sup>δx<sup>k</sup>þ<sup>1</sup>þB<sup>k</sup>þ<sup>1</sup>δυ<sup>k</sup>þ<sup>1</sup>� ¼ <sup>0</sup>

<sup>k</sup>þ<sup>1</sup>Qδx<sup>k</sup>þ<sup>1</sup> <sup>þ</sup> <sup>R</sup> <sup>δ</sup>υ<sup>2</sup>

<sup>n</sup> δυ<sup>0</sup> δυ<sup>1</sup> … δυ<sup>n</sup> δυ<sup>0</sup>:<sup>5</sup> δυ<sup>1</sup>:<sup>5</sup> … δυ<sup>n</sup>�0:<sup>5</sup>

<sup>k</sup> <sup>þ</sup> <sup>4</sup>δυ<sup>2</sup>

δx<sup>k</sup> ¼ xð Þ� τ<sup>k</sup> xoptð Þ τ<sup>k</sup> (9)

<sup>8</sup> ð Þ <sup>A</sup>kδx<sup>k</sup> <sup>þ</sup> <sup>B</sup>kδυ<sup>k</sup> � <sup>A</sup><sup>k</sup>þ<sup>1</sup>δx<sup>k</sup>þ<sup>1</sup> � <sup>B</sup><sup>k</sup>þ<sup>1</sup>δυ<sup>k</sup>þ<sup>1</sup> (10)

Z<sup>T</sup>MZ (11)

<sup>0</sup> 0 0 … <sup>0</sup> � �<sup>T</sup> (12)

X ≜ K ti ð Þ ; n; th xð Þ� ti xoptð Þti

X (13)

� � (14)

<sup>k</sup>þ0:<sup>5</sup> <sup>þ</sup> <sup>δ</sup>υ<sup>2</sup>

kþ1

(7)

(8)

, which mini-

<sup>G</sup> <sup>¼</sup> <sup>1</sup> 2 δx<sup>T</sup>

<sup>k</sup>þ0:<sup>5</sup>Qδx<sup>k</sup>þ0:<sup>5</sup> <sup>þ</sup> <sup>δ</sup>x<sup>T</sup>

γ

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

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

C<sup>T</sup> CM�<sup>1</sup> C<sup>T</sup> � ��<sup>1</sup>

where the row vector V is defined to "choose" the target value from the optimal solution, and the position of "1" in the row vector V is the same as the position of δυ<sup>0</sup> in the column vector Z.

<sup>G</sup> <sup>¼</sup> <sup>1</sup> 2

CZ¼X X <sup>¼</sup> <sup>δ</sup>x<sup>T</sup>

� �<sup>T</sup>

used to discretize the RHC problem.

Simpson integration formula as

98 Optimization Algorithms - Examples

δx<sup>T</sup>

and the constraints are discretized into

<sup>δ</sup>x<sup>k</sup>þ0:<sup>5</sup> <sup>¼</sup> <sup>1</sup>

gramming vector <sup>Z</sup> <sup>¼</sup> <sup>δ</sup>x<sup>T</sup>

mizes the cost function:

subject to.

2

ð Þþ δx<sup>k</sup> þ δx<sup>k</sup>þ<sup>1</sup>

<sup>0</sup> δx<sup>T</sup>

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

δυð Þ¼ ti VZ<sup>∗</sup>

where, A<sup>k</sup> ¼ Að Þ τ<sup>k</sup> , A<sup>k</sup>þ0:<sup>5</sup> ¼ Að Þ τ<sup>k</sup>þ0:<sup>5</sup> , B<sup>k</sup> ¼ Bð Þ τ<sup>k</sup> , B<sup>k</sup>þ0:<sup>5</sup> ¼ Bð Þ τ<sup>k</sup>þ0:<sup>5</sup> .

<sup>1</sup> … δx<sup>T</sup>

Z∗ <sup>¼</sup>M�<sup>1</sup>

C<sup>T</sup> CM�<sup>1</sup> C<sup>T</sup> � ��<sup>1</sup>

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

<sup>¼</sup>VM�<sup>1</sup>

Finally, the control input of the closed-loop control, υð Þti , is

<sup>k</sup> <sup>Q</sup>δx<sup>k</sup> <sup>þ</sup> <sup>4</sup>δx<sup>T</sup>

þ γ 12 Xn�1 k¼0

<sup>δ</sup>x<sup>k</sup> � <sup>δ</sup>x<sup>k</sup>þ<sup>1</sup> <sup>þ</sup> <sup>γ</sup>

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

$$\frac{da}{dt} = \frac{2}{n\sqrt{1 - e^2}} \left( \sigma\_z e \sin \nu + \sigma\_x \frac{p}{r} \right) \tag{16}$$

$$\frac{d\Omega}{dt} = \frac{\sigma\_y r \sin u}{na^2 \sqrt{1 - e^2} \sin i} \tag{17}$$

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

$$\frac{di}{dt} = \frac{\sigma\_y r \cos u}{na^2 \sqrt{1 - e^2}}\tag{18}$$

<sup>β</sup>€ <sup>þ</sup> ð Þ <sup>α</sup>\_ <sup>þ</sup> <sup>ν</sup>\_ <sup>2</sup> sin <sup>β</sup> cos <sup>β</sup> <sup>þ</sup> <sup>3</sup>μ<sup>r</sup>

corresponding perturbative torques by the perturbative forces to be discussed below.

hg ¼ r � rpo 1 � e

m¼0

Xn m¼0

gm

where <sup>r</sup><sup>0</sup> = 6371.2 � <sup>10</sup><sup>3</sup> km is the reference radius of the Earth, respectively.

where the polar radius rpo and the Earth's eccentricity eE are provided by NASA [19].

m gm

<sup>n</sup> cos <sup>m</sup><sup>ϕ</sup> � � <sup>þ</sup> <sup>h</sup><sup>m</sup>

gm

Iave <sup>¼</sup> <sup>1</sup> L ðL 0

The open-loop optimal control problem for EDT deorbit can be stated as finding a state-control pair f g xð Þt ; υð Þt over each time interval ti ½ � ; tiþ<sup>1</sup> to minimize a cost function of the negative work

> ðtiþ<sup>1</sup> tt

2 <sup>E</sup> cos <sup>2</sup>

<sup>n</sup> sin <sup>m</sup><sup>ϕ</sup> � � � hm

<sup>n</sup> sin <sup>m</sup><sup>ϕ</sup> � � � � <sup>∂</sup>Pm

<sup>n</sup> cos <sup>m</sup><sup>ϕ</sup> � � <sup>þ</sup> hm

<sup>n</sup> cos <sup>m</sup><sup>ϕ</sup> � � � � <sup>P</sup><sup>m</sup>

<sup>n</sup> sin <sup>m</sup><sup>ϕ</sup> � � � � <sup>P</sup><sup>m</sup>

Moreover, the local strength of geomagnetic field is described by the IGRF2000 model [20–22]

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

<sup>t</sup> <sup>=</sup><sup>12</sup> � �m�<sup>1</sup>

in a body-fixed spherical coordinates of the Earth, such that

X∞ n¼1

r0 r � �<sup>n</sup>þ<sup>2</sup>X<sup>n</sup>

r0 r � �<sup>n</sup>þ<sup>2</sup>

r0 r � �<sup>n</sup>þ<sup>2</sup>X<sup>n</sup>

m¼0

J ¼ ðtiþ<sup>1</sup> ti Fe ! � v ! dt≜

subject to the nonlinear state equations of libration motion

ð Þ n þ 1

<sup>B</sup><sup>ϕ</sup> <sup>¼</sup> <sup>1</sup> sin θ

<sup>B</sup><sup>θ</sup> <sup>¼</sup> <sup>X</sup><sup>∞</sup> n¼1

Br <sup>¼</sup> <sup>X</sup><sup>∞</sup> n¼1

The average current in the EDT is defined as

done by the electrodynamic force

where <sup>m</sup><sup>~</sup> <sup>¼</sup> <sup>m</sup>1m<sup>2</sup> <sup>þ</sup> ð Þ <sup>m</sup><sup>1</sup> <sup>þ</sup> <sup>m</sup><sup>2</sup> mt=<sup>3</sup> <sup>þ</sup> <sup>m</sup><sup>2</sup>

ellipsoidal surface, such that,

�<sup>3</sup> cos <sup>2</sup>

<sup>α</sup> sin <sup>β</sup> cos <sup>β</sup> <sup>¼</sup> <sup>Q</sup><sup>β</sup>

EDT is the equivalent mass, (Qα, Qβ) are the

θ � ��1=<sup>2</sup> (24)

<sup>n</sup> ð Þ θ<sup>c</sup> ∂θ<sup>c</sup>

<sup>n</sup> ð Þ θ<sup>c</sup>

(25)

<sup>n</sup> ð Þ θ<sup>c</sup>

I sð Þds (26)

Πð Þ x; υ; t dt (27)

mL<sup>~</sup> <sup>2</sup> (23)

101

Piecewise Parallel Optimal Algorithm http://dx.doi.org/10.5772/intechopen.76625

$$\frac{d\mathbf{e}\_x}{dt} = \frac{\sqrt{1-\dot{\mathbf{c}}^2}}{na} \left\{ \sigma\_z \sin u + \sigma\_x \left[ \left( 1 + \frac{r}{p} \right) \cos u + \frac{r}{p} \mathbf{e}\_x \right] \right\} + \frac{d\Omega}{dt} \mathbf{e}\_y \cos i \tag{19}$$

$$\frac{d\sigma\_y}{dt} = \frac{\sqrt{1-\varepsilon^2}}{na} \left\{-\sigma\_z \cos u + \sigma\_x \left[ \left(1 + \frac{r}{p}\right) \sin u + \frac{r}{p} e\_y \right] \right\} - \frac{d\Omega}{dt} e\_x \cos i \tag{20}$$

$$\frac{d\phi}{dt} = n - \frac{1}{na} \left[ \sigma\_z \left( \frac{2r}{a} + \frac{\sqrt{1-e^2}}{1+\sqrt{1-e^2}} e \cos \nu \right) - \sigma\_x \left( 1 + \frac{r}{p} \right) \frac{\sqrt{1-e^2}}{1+\sqrt{1-e^2}} e \sin \nu \right] - \frac{\sigma\_y r \cos i \sin \nu}{na^2 \sqrt{1-e^2} \sin i} \tag{21}$$

The components of perturbative accelerations are defined in a local frame. The components σ<sup>x</sup> and σ<sup>z</sup> are in the orbital plane, and σ<sup>z</sup> is the radial component pointing outwards. The out-ofplane component σ<sup>y</sup> completes a right-hand coordinate system. The components of perturbative accelerations depend on the tether attitude and the EDT' orbital dynamics is coupled with the tether libration motion.

The libration motion of a rigid EDT system is described in an orbital coordinate system shown in Figure 2(b). The z-axis of the orbital coordinate system points from the Earth's center to the CM of the EDT system, the x-axis lies in the orbital plane and points to the direction of the EDT orbital motion, perpendicular to the z-axis. The y-axis completes a right-hand coordinate system. The unit vectors along each axis are expressed as 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,

$$
\ddot{a} + \ddot{\nu} - 2(\dot{a} + \dot{\nu})\dot{\beta}\tan\beta + 3\mu r^{-3}\sin a \cos a = \frac{Q\_a}{\ddot{m}L^2\cos^2\beta} \tag{22}
$$

Piecewise Parallel Optimal Algorithm http://dx.doi.org/10.5772/intechopen.76625 101

$$\ddot{\beta} + (\dot{a} + \dot{\nu})^2 \sin \beta \cos \beta + 3 \mu r^{-3} \cos^2 \alpha \sin \beta \cos \beta = \frac{Q\_{\beta}}{\tilde{m}L^2} \tag{23}$$

where <sup>m</sup><sup>~</sup> <sup>¼</sup> <sup>m</sup>1m<sup>2</sup> <sup>þ</sup> ð Þ <sup>m</sup><sup>1</sup> <sup>þ</sup> <sup>m</sup><sup>2</sup> mt=<sup>3</sup> <sup>þ</sup> <sup>m</sup><sup>2</sup> <sup>t</sup> <sup>=</sup><sup>12</sup> � �m�<sup>1</sup> EDT is 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,

$$h\_{\rm g} = r - r\_{\rm pv} \left(1 - e\_E^2 \cos^2 \theta\right)^{-1/2} \tag{24}$$

where the polar radius rpo and the Earth's eccentricity eE are provided by NASA [19].

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

$$\begin{aligned} B\_{\phi} &= \frac{1}{\sin\Theta} \sum\_{n=1}^{\infty} \left( \frac{r\_0}{r} \right)^{n+2} \sum\_{m=0}^{n} m \left[ g\_n^m \sin \left( m\phi \right) - h\_n^m \cos \left( m\phi \right) \right] P\_n^m(\theta\_c) \\ B\_{\theta} &= \sum\_{n=1}^{\infty} \left( \frac{r\_0}{r} \right)^{n+2} \sum\_{m=0}^{n} \left[ g\_n^m \cos \left( m\phi \right) + h\_n^m \sin \left( m\phi \right) \right] \frac{\partial P\_n^m(\theta\_c)}{\partial \theta\_c} \\ B\_r &= \sum\_{n=1}^{\infty} \left( \frac{r\_0}{r} \right)^{n+2} \left( n+1 \right) \sum\_{m=0}^{n} \left[ g\_n^m \cos \left( m\phi \right) + h\_n^m \sin \left( m\phi \right) \right] P\_n^m(\theta\_c) \end{aligned} \tag{25}$$

where <sup>r</sup><sup>0</sup> = 6371.2 � <sup>10</sup><sup>3</sup> km is the reference radius of the Earth, respectively.

The average current in the EDT is defined as

di

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

σ<sup>z</sup> sin u þ σ<sup>x</sup> 1 þ

na �σ<sup>z</sup> cos <sup>u</sup> <sup>þ</sup> <sup>σ</sup><sup>x</sup> <sup>1</sup> <sup>þ</sup>

<sup>1</sup> � <sup>e</sup><sup>2</sup> <sup>p</sup> <sup>e</sup> cos <sup>ν</sup>

dex dt ¼

dey dt ¼

dφ

dt <sup>¼</sup> <sup>n</sup> � <sup>1</sup>

na σz 2r a þ

100 Optimization Algorithms - Examples

tether libration motion.

ffiffiffiffiffiffiffiffiffiffiffiffi <sup>1</sup> � <sup>e</sup><sup>2</sup> <sup>p</sup> na

ffiffiffiffiffiffiffiffiffiffiffiffi <sup>1</sup> � <sup>e</sup><sup>2</sup> <sup>p</sup>

> ffiffiffiffiffiffiffiffiffiffiffiffi <sup>1</sup> � <sup>e</sup><sup>2</sup> <sup>p</sup> <sup>1</sup> <sup>þ</sup> ffiffiffiffiffiffiffiffiffiffiffiffi

system. The unit vectors along each axis are expressed as eox

<sup>α</sup>€ <sup>þ</sup> <sup>ν</sup>€ � <sup>2</sup>ð Þ <sup>α</sup>\_ <sup>þ</sup> <sup>ν</sup>\_ <sup>β</sup>\_ tan <sup>β</sup> <sup>þ</sup> <sup>3</sup>μ<sup>r</sup>

motion of the EDT system can be derived as,

!

dt <sup>¼</sup> <sup>σ</sup>yr cos <sup>u</sup> na<sup>2</sup> ffiffiffiffiffiffiffiffiffiffiffiffi

> r p � �

> > r p � �

> > > r p

� � � �

� � � �

� σ<sup>x</sup> 1 þ

The components of perturbative accelerations are defined in a local frame. The components σ<sup>x</sup> and σ<sup>z</sup> are in the orbital plane, and σ<sup>z</sup> is the radial component pointing outwards. The out-ofplane component σ<sup>y</sup> completes a right-hand coordinate system. The components of perturbative accelerations depend on the tether attitude and the EDT' orbital dynamics is coupled with the

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

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

" #

cos u þ

sin u þ

� � ffiffiffiffiffiffiffiffiffiffiffiffi

r p ex

> r p ey

<sup>1</sup> � <sup>e</sup><sup>2</sup> <sup>p</sup> <sup>1</sup> <sup>þ</sup> ffiffiffiffiffiffiffiffiffiffiffiffi

! , eoy

�<sup>3</sup> sin <sup>α</sup> cos <sup>α</sup> <sup>¼</sup> <sup>Q</sup><sup>α</sup>

! , and eoz

<sup>1</sup> � <sup>e</sup><sup>2</sup> <sup>p</sup> (18)

� <sup>d</sup><sup>Ω</sup>

dt ey cos <sup>i</sup> (19)

dt ex cos <sup>i</sup> (20)

� <sup>σ</sup>yr cos <sup>i</sup> sin <sup>u</sup> na<sup>2</sup> ffiffiffiffiffiffiffiffiffiffiffiffi <sup>1</sup> � <sup>e</sup><sup>2</sup> <sup>p</sup> sin <sup>i</sup>

! , respectively. Then,

mL<sup>~</sup> <sup>2</sup> cos <sup>2</sup><sup>β</sup> (22)

(21)

þ dΩ

<sup>1</sup> � <sup>e</sup><sup>2</sup> <sup>p</sup> <sup>e</sup> sin <sup>ν</sup>

$$I\_{\rm ave} = \frac{1}{L} \int\_0^L I(\mathbf{s}) d\mathbf{s} \tag{26}$$

The open-loop optimal control problem for EDT deorbit can be stated as finding a state-control pair f g xð Þt ; υð Þt over each time interval ti ½ � ; tiþ<sup>1</sup> to minimize a cost function of the negative work done by the electrodynamic force

$$J = \int\_{t\_l}^{t\_{l+1}} \overrightarrow{F}\_{\varepsilon} \cdot \overrightarrow{v} \, dt \triangleq \int\_{t\_l}^{t\_{l+1}} \Pi(\mathfrak{x}, \upsilon, t) dt \tag{27}$$

subject to the nonlinear state equations of libration motion

$$\begin{aligned} \mathbf{x}'\_1 &= \mathbf{x}\_2 \\ \mathbf{x}'\_2 &= 2\eta^3 \varepsilon \sin \upsilon + 2 \left( \mathbf{x}\_2 + \eta^2 \right) \mathbf{x}\_4 \tan \mathbf{x}\_3 - 3\eta^3 \sin \mathbf{x}\_1 \cos \mathbf{x}\_1 \\ &+ \left[ \sin i \tan \mathbf{x}\_3 (2 \sin u \cos \mathbf{x}\_1 - \cos u \sin \mathbf{x}\_1) - \cos i \right] (\varepsilon - \lambda) I\_{\text{arc}} \frac{\mu\_m}{\mu^m} \eta^3 \\\\ \mathbf{x}'\_3 &= \mathbf{x}\_4 \end{aligned} \tag{28}$$
 
$$\begin{aligned} \mathbf{x}'\_4 &= - \left( \mathbf{x}\_2 + \eta^2 \right)^2 \sin \mathbf{x}\_3 \cos \mathbf{x}\_3 - 3\eta^3 \cos^2 \mathbf{x}\_1 \sin \mathbf{x}\_3 \cos \mathbf{x}\_3 \\ &- \sin i \left( 2 \sin u \sin \mathbf{x}\_1 + \cos u \cos \mathbf{x}\_1 \right) (\varepsilon - \lambda) I\_{\text{arc}} \frac{\mu\_m}{\mu m} \eta^3 \end{aligned} \tag{29}$$

where ð Þ¼ x1; x2; x3; x<sup>4</sup> α; α<sup>0</sup> ; <sup>β</sup>; ; <sup>β</sup><sup>0</sup> � �, <sup>η</sup> <sup>¼</sup> <sup>1</sup> <sup>þ</sup> <sup>e</sup> cos <sup>ν</sup>, <sup>ν</sup>\_ <sup>¼</sup> <sup>μ</sup>=p<sup>3</sup> � �<sup>0</sup>:<sup>5</sup> ð Þ 1 þ e cos ν 2 , r ¼ pð1 þ e cos νÞ �1 , λ ¼ ð Þ m<sup>1</sup> þ 0:5mt =mEDT is determined by the mass ratio between the end-bodies, and ς is determined by the distribution of current along the EDT, such that, ς ¼ I �1 aveL�<sup>2</sup> <sup>Ð</sup> <sup>L</sup> <sup>0</sup> sI sð Þds. Accordingly, ς ¼ 0:5 is used for the assumption of a constant current in the EDT. The initial conditions xð Þ¼ ti xstart and the box constraint j j α ≤ αmax, β � � � � ≤ βmax, Imin ≤ Iave ≤ Imax. The environmental perturbations are simplified by considering only the electrodynamic force with a simple non-tilted dipole model of geomagnetic field, such that,

$$\overrightarrow{\mathbf{B}} = \frac{\mu\_m}{r^3} \cos u \sin i \ \overrightarrow{\mathbf{e}}\_{ox} + \frac{\mu\_m}{r^3} \cos i \ \overrightarrow{\mathbf{e}}\_{oy} - \frac{2\mu\_m}{r^3} \sin u \sin i \ \overrightarrow{\mathbf{e}}\_{az} \tag{30}$$

Accordingly, the electrodynamic force F<sup>e</sup> ! exerting on the EDT can be obtained as,

$$\stackrel{\rightarrow}{\mathbf{F}\_{\epsilon}} = -\int\_{0}^{L} \stackrel{\rightarrow}{\mathbf{B}} \times I \stackrel{\rightarrow}{\mathbf{I}} \, ds = -I\_{dve}L\stackrel{\rightarrow}{\mathbf{B}} \times \stackrel{\rightarrow}{\mathbf{I}}\tag{31}$$

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

cos �<sup>1</sup> <sup>i</sup> � <sup>β</sup><sup>B</sup>

Piecewise Parallel Optimal Algorithm http://dx.doi.org/10.5772/intechopen.76625

A

103

<sup>I</sup>max (inclined) <sup>0</sup>:<sup>1</sup> � sin <sup>i</sup> sin <sup>β</sup><sup>B</sup> sin <sup>ð</sup>Ω<sup>G</sup> <sup>þ</sup> <sup>α</sup><sup>B</sup> � <sup>Ω</sup>Þ þ cos <sup>i</sup> cos <sup>β</sup><sup>B</sup>

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

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

Figure 3. Time history of average current in the equatorial orbit (non-tilted dipole geomagnetic model). Solid line: Open-

geomagnetic fields are conducted respectively.

loop state trajectory. Dashed line: Close-loop tracking trajectory.

Parameters Values Imax (equatorial) 0.4 A

Imin (equatorial & inclined) 0 A αmax (equatorial & inclined) 45 degrees βmax (equatorial & inclined) 45 degrees

#### 3.1.2. Results and discussion

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


Table 1. Parameters of an EDT system.


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

x0 <sup>1</sup> ¼ x<sup>2</sup>

102 Optimization Algorithms - Examples

x0 <sup>2</sup> <sup>¼</sup> <sup>2</sup>η<sup>3</sup>

where ð Þ¼ x1; x2; x3; x<sup>4</sup> α; α<sup>0</sup>

cos νÞ �1 x0 <sup>3</sup> ¼ x<sup>4</sup>

x0

B ! <sup>¼</sup> <sup>μ</sup><sup>m</sup>

Accordingly, the electrodynamic force F<sup>e</sup>

3.1.2. Results and discussion

Table 1. Parameters of an EDT system.

<sup>4</sup> ¼ � <sup>x</sup><sup>2</sup> <sup>þ</sup> <sup>η</sup><sup>2</sup> � �<sup>2</sup>

conditions xð Þ¼ ti xstart and the box constraint j j α ≤ αmax, β

simple non-tilted dipole model of geomagnetic field, such that,

<sup>r</sup><sup>3</sup> cos <sup>u</sup> sin <sup>i</sup> <sup>e</sup>

Fe ! ¼ � ðL 0 B ⇀ �I l ⇀

<sup>e</sup> sin <sup>ν</sup> <sup>þ</sup> <sup>2</sup> <sup>x</sup><sup>2</sup> <sup>þ</sup> <sup>η</sup><sup>2</sup> � �x<sup>4</sup> tan <sup>x</sup><sup>3</sup> � <sup>3</sup>η<sup>3</sup> sin <sup>x</sup><sup>1</sup> cos <sup>x</sup><sup>1</sup>

þ ½ � sin i tan x3ð2 sin u cos x<sup>1</sup> � cos u sin x1Þ � cos i ð Þ ς � λ Iave

� sin ið Þ 2 sin u sin x<sup>1</sup> þ cos u cos x<sup>1</sup> ð Þ ς � λ Iave

; <sup>β</sup>; ; <sup>β</sup><sup>0</sup> � �, <sup>η</sup> <sup>¼</sup> <sup>1</sup> <sup>þ</sup> <sup>e</sup> cos <sup>ν</sup>, <sup>ν</sup>\_ <sup>¼</sup> <sup>μ</sup>=p<sup>3</sup> � �<sup>0</sup>:<sup>5</sup>

Accordingly, ς ¼ 0:5 is used for the assumption of a constant current in the EDT. The initial

ronmental perturbations are simplified by considering only the electrodynamic force with a

<sup>r</sup><sup>3</sup> cos <sup>i</sup> <sup>e</sup> ⇀

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

ς is determined by the distribution of current along the EDT, such that, ς ¼ I

⇀ ox <sup>þ</sup> <sup>μ</sup><sup>m</sup>

Parameters Values Mass of the main satellite 5 kg Mass of subsatellite 1.75 kg Mass of the tether 0.25 kg

Tether length 500 m Tether diameter 0.0005 m

Dimensions of main satellite 0.2 � 0.2 � 0.2 m Dimensions of subsatellite 0.1 � 0.17 � 0.1 m

Tether conductivity (aluminum) 3.4014 � <sup>10</sup><sup>7</sup> <sup>Ω</sup>�<sup>1</sup> <sup>m</sup>�<sup>1</sup>

Orbital altitudes 700 ~ 800 km

Tether current lower/upper limits 0 ~ 0.8 A for the equatorial orbit

!

, λ ¼ ð Þ m<sup>1</sup> þ 0:5mt =mEDT is determined by the mass ratio between the end-bodies, and

� � �

oy � <sup>2</sup>μ<sup>m</sup>

ds ¼ �IaveL B

<sup>r</sup><sup>3</sup> sin <sup>u</sup> sin <sup>i</sup> <sup>e</sup>

exerting on the EDT can be obtained as,

⇀ � l ⇀

sin <sup>x</sup><sup>3</sup> cos <sup>x</sup><sup>3</sup> � <sup>3</sup>η<sup>3</sup> cos <sup>2</sup>

μm <sup>μ</sup>m<sup>~</sup> <sup>η</sup><sup>3</sup>

x<sup>1</sup> sin x<sup>3</sup> cos x<sup>3</sup>

μm <sup>μ</sup>m<sup>~</sup> <sup>η</sup><sup>3</sup>

ð Þ 1 þ e cos ν

2

�1 aveL�<sup>2</sup> <sup>Ð</sup> <sup>L</sup>

oz (30)

� ≤ βmax, Imin ≤ Iave ≤ Imax. The envi-

⇀

, r ¼ pð1 þ e

<sup>0</sup> sI sð Þds.

(31)

(28)

(29)

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

Figure 3. Time history of average current in the equatorial orbit (non-tilted dipole geomagnetic model). Solid line: Openloop 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.

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

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

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

Piecewise Parallel Optimal Algorithm http://dx.doi.org/10.5772/intechopen.76625 105

not lead to a tumbling of the EDT system.

trajectory. Dashed line: Close-loop tracking trajectory.

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

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.

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

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.

104 Optimization Algorithms - Examples

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.

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

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

mainly concern obtaining the local time optimization.

3.2.1. Problem formulation

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

Consider an in-plane PSE system in a circular orbit is shown in Figure 8, where the main satellite, climber and the end body are connected by two inelastic tethers L<sup>1</sup> and 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, m<sup>1</sup> and 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 r measuring from the centre of Earth. The climber m<sup>1</sup> is connected to the main satellite M by a tether 1 with the length of L<sup>1</sup> and a libration angle θ<sup>1</sup> measured from the vector r. The distance between them is controlled by reeling in/out tether 1 at main satellite. The end body m<sup>2</sup> is connected to m<sup>1</sup> by a tether 2 with the length of L<sup>2</sup> and a libration angle θ<sup>2</sup> measured from the vector r. The length of tether 2 L<sup>2</sup> is controlled by reeling in or out tether 2 at end body. The mass of the main satellite is assumed much greater than the masses of the climber and the end body. Therefore, the CM of the PSE system can be assumed residing in the main satellite that moves in a circular orbit.

Based on the aforementioned assumptions, the dynamic equations can be written as.

<sup>2</sup> � <sup>2</sup> <sup>ω</sup> <sup>þ</sup> <sup>θ</sup>\_

<sup>θ</sup><sup>1</sup> <sup>þ</sup> <sup>2</sup>ωL1θ\_ <sup>1</sup> <sup>þ</sup> <sup>L</sup>1θ\_

m<sup>1</sup>

where L<sup>0</sup> is the initial total length of two pieces of the tethers and Lc is the length increment

The libration angles are required to be kept between the desired upper/lower bounds in the climber transfer process. The accessing process is divided into a series of intervals. In this case, we modified the aforementioned parallel optimal algorithm. The total transfer length L<sup>10</sup> � L1<sup>f</sup>

L\_

<sup>2</sup> � <sup>2</sup> <sup>ω</sup> <sup>þ</sup> <sup>θ</sup>\_ <sup>2</sup>

1 L\_

L1

L<sup>0</sup> � L<sup>1</sup> þ Lc

1

<sup>c</sup> � <sup>L</sup>\_ 1 

> 1 2 � <sup>T</sup><sup>1</sup> m<sup>1</sup> þ

þ

� sin ð Þ <sup>θ</sup><sup>1</sup> � <sup>θ</sup><sup>2</sup> <sup>T</sup><sup>2</sup> L1m<sup>1</sup>

> sin ð Þ θ<sup>1</sup> � θ<sup>2</sup> T<sup>1</sup> ð Þ L<sup>0</sup> � L<sup>1</sup> þ Lc m<sup>1</sup>

Piecewise Parallel Optimal Algorithm http://dx.doi.org/10.5772/intechopen.76625 107

cos ð Þ θ<sup>1</sup> � θ<sup>2</sup> T<sup>2</sup> m<sup>1</sup>

� ½ � <sup>m</sup><sup>1</sup> � <sup>m</sup><sup>2</sup> cos ð Þþ <sup>θ</sup><sup>1</sup> � <sup>θ</sup><sup>2</sup> <sup>m</sup><sup>2</sup> <sup>T</sup><sup>2</sup> m1m<sup>2</sup>

ð Þ <sup>L</sup><sup>0</sup> � <sup>L</sup><sup>1</sup> <sup>þ</sup> Lc <sup>θ</sup>\_ <sup>2</sup>

(32)

(33)

(34)

(35)

<sup>θ</sup>€<sup>1</sup> ¼ � <sup>3</sup>ω<sup>2</sup> sin 2θ<sup>1</sup>

L<sup>1</sup> cos <sup>2</sup>

<sup>L</sup>€<sup>c</sup> <sup>¼</sup> <sup>3</sup>ω<sup>2</sup>ð Þ <sup>L</sup><sup>0</sup> � <sup>L</sup><sup>1</sup> <sup>þ</sup> Lc cos <sup>2</sup>θ<sup>2</sup> <sup>þ</sup> <sup>3</sup>ω<sup>2</sup>L<sup>1</sup> cos <sup>2</sup>θ<sup>1</sup> <sup>þ</sup> <sup>2</sup><sup>ω</sup> <sup>þ</sup> <sup>θ</sup>\_ <sup>2</sup>

<sup>L</sup>1θ\_ <sup>1</sup> <sup>þ</sup> ½ � cos ð Þ� <sup>θ</sup><sup>1</sup> � <sup>θ</sup><sup>2</sup> <sup>1</sup> <sup>T</sup><sup>1</sup>

<sup>θ</sup>€<sup>2</sup> <sup>¼</sup> �3ω<sup>2</sup> sin 2θ<sup>2</sup>

<sup>L</sup>€<sup>1</sup> <sup>¼</sup> <sup>3</sup>ω<sup>2</sup>

1

<sup>þ</sup> <sup>2</sup><sup>ω</sup> <sup>þ</sup> <sup>θ</sup>\_

relates to L0.

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

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

ing control.

106 Optimization Algorithms - Examples

Figure 8. PSE system.

Consider an in-plane PSE system in a circular orbit is shown in Figure 8, where the main satellite, climber and the end body are connected by two inelastic tethers L<sup>1</sup> and 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, m<sup>1</sup> and 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 r measuring from the centre of Earth. The climber m<sup>1</sup> is connected to the main satellite M by a tether 1 with the length of L<sup>1</sup> and a libration angle θ<sup>1</sup> measured from the vector r. The distance between them is controlled by reeling in/out tether 1 at main satellite. The end body m<sup>2</sup> is connected to m<sup>1</sup> by a tether 2 with the length of L<sup>2</sup> and a libration angle θ<sup>2</sup> measured from the vector r. The length of tether 2 L<sup>2</sup> is controlled by reeling in or out tether 2 at end body. The mass of the main satellite is assumed much greater than the masses of the climber and the end body. Therefore, the CM of the PSE system can be assumed residing in the main satellite that moves in a circular orbit. Based on the aforementioned assumptions, the dynamic equations can be written as.

$$\ddot{\theta}\_1 = -\frac{3\omega^2 \sin 2\theta\_1}{2} - \frac{2\left(\omega + \dot{\theta}\_1\right)\dot{L}\_1}{L\_1} - \frac{\sin\left(\theta\_1 - \theta\_2\right)T\_2}{L\_1 m\_1} \tag{32}$$

$$\ddot{\theta}\_2 = \frac{-3\omega^2 \sin 2\theta\_2}{2} - \frac{2\left(\omega + \dot{\theta}\_2\right)\left(\dot{L}\_c - \dot{L}\_1\right)}{L\_0 - L\_1 + L\_c} + \frac{\sin\left(\theta\_1 - \theta\_2\right)T\_1}{(L\_0 - L\_1 + L\_c)m\_1} \tag{33}$$

$$\ddot{L}\_1 = 3\omega^2 L\_1 \cos^2 \theta\_1 + 2\omega L\_1 \dot{\theta}\_1 + L\_1 \dot{\theta}\_1^2 - \frac{T\_1}{m\_1} + \frac{\cos\left(\theta\_1 - \theta\_2\right)T\_2}{m\_1} \tag{34}$$

$$\begin{aligned} \ddot{L}\_c &= 3\omega^2 (L\_0 - L\_1 + L\_c) \cos^2 \theta\_2 + 3\omega^2 L\_1 \cos^2 \theta\_1 + (2\omega + \dot{\theta}\_2)(L\_0 - L\_1 + L\_c) \dot{\theta}\_2 \\ &+ (2\omega + \dot{\theta}\_1) L\_1 \dot{\theta}\_1 + \frac{[\cos\left(\theta\_1 - \theta\_2\right) - 1]T\_1}{m\_1} - \frac{[m\_1 - m\_2 \cos\left(\theta\_1 - \theta\_2\right) + m\_2]T\_2}{m\_1 m\_2} \end{aligned} \tag{35}$$

where L<sup>0</sup> is the initial total length of two pieces of the tethers and Lc is 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 L<sup>10</sup> � L1<sup>f</sup>

of the climber is discretized evenly. The optimal trajectory should be found to make the transfer time minimize in each equal tether transferring length. Then the cost function can be rewritten as

$$J\_i = \mathbf{t}\_f \tag{36}$$

m1 = 500 kg, m<sup>2</sup> = 1000 kg, θ1(0) = θ2(0) = 0, L<sup>0</sup> = 20 km, L1(0) = 19,500 m, Lc (0) = 0,

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:3 rad,

cm ¼ �<sup>10</sup> <sup>m</sup>=<sup>s</sup> and <sup>L</sup>\_

between 40s and 350 s. Figure 10 shows the time history of trajectories of L\_

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 θ<sup>1</sup> is kept at 0.3 rad by the end of the ascending period, see the dashed line in Figure 10. Using the closed-loop control, the tracking trajectory of θ<sup>1</sup> matches the open-loop trajectory very well overall, see solid line in Figure 9. A short gap appears between 875 s – 880 s, this is caused by the errors of the model and computation. Figure 9 also shows the changes of the trajectories of θ2.The trajectory of θ<sup>2</sup> obtained by closed-loop control tracks the open-loop trajectory well and reaches 0.1 rad by the end of the ascending period. The closed-loop trajectories of L<sup>1</sup> and Lc are shown in Figure 10. They are the reflections of the control inputs. Both L<sup>1</sup> and Lc show smooth fluctuations

<sup>c</sup>ð Þ¼ 0 0 for the ascending process. The whole

Piecewise Parallel Optimal Algorithm http://dx.doi.org/10.5772/intechopen.76625

<sup>1</sup> and L\_

c,

109

cM ¼ 10 m=s.

<sup>1</sup>ð Þ¼� <sup>0</sup> <sup>20</sup> <sup>m</sup>=s, and <sup>L</sup>\_

<sup>1</sup><sup>M</sup> ¼ �<sup>25</sup> <sup>m</sup>=s, <sup>L</sup>\_

Figure 9. Libration angle of θ<sup>1</sup> and θ<sup>2</sup> with variable climber speed.

Figure 10. Length and its changing ratio of L1, Lc with variable climber speed.

θ\_

L\_

<sup>1</sup>ð Þ¼ <sup>0</sup> <sup>θ</sup>\_ <sup>2</sup>ð Þ¼ <sup>0</sup> 0, <sup>L</sup>\_

<sup>1</sup><sup>m</sup> ¼ �<sup>15</sup> <sup>m</sup>=s, <sup>L</sup>\_

subject to the simplified dynamic equations

$$\ddot{\theta}\_1 = -3\omega^2 \theta\_1 - \frac{2(\omega + \dot{\theta}\_1)\dot{L}\_1}{L\_1} - \frac{(\theta\_1 - \theta\_2)T\_2}{L\_1 m\_1} \tag{37}$$

$$\ddot{\theta}\_2 = -3\omega^2 \theta\_2 - \frac{2\left(\omega + \dot{\theta}\_2\right)\left(\dot{L}\_c - \dot{L}\_1\right)}{L\_0 - L\_1 + L\_c} + \frac{(\theta\_1 - \theta\_2)T\_1}{(L\_0 - L\_1 + L\_c)m\_1} \tag{38}$$

$$\ddot{L}\_1 = 3\omega^2 L\_1 + 2\omega L\_1 \dot{\theta}\_1 + L\_1 \dot{\theta}\_1^{\;\;2} + \frac{T\_2 - T\_1}{m\_1} \tag{39}$$

$$\ddot{L}\_c = 3\omega^2 (L\_0 + L\_c) + \left(2\omega + \dot{\theta}\_2\right)(L\_0 - L\_1 + L\_c)\dot{\theta}\_2 + \left(2\omega + \dot{\theta}\_1\right)L\_1\dot{\theta}\_1 - \frac{T\_2}{m\_2} \tag{40}$$

where <sup>u</sup> <sup>¼</sup> ð Þ <sup>T</sup>1; <sup>T</sup><sup>2</sup> , <sup>x</sup> <sup>¼</sup> <sup>θ</sup>1; <sup>θ</sup>\_ <sup>1</sup>; θ2; θ\_ <sup>2</sup>; L1; L\_ 1 and i denotes the interval number. In (26) the gravitational perturbations and the trigonometric functions are ignored, then they can be simplified following the assumptions: sin θ<sup>j</sup> � θj, cos θ<sup>j</sup> � 1 ð Þ j ¼ 1; 2 . All the errors between simple model and the entire model are regarded as perturbations.

To ensure the availability and the suppression of the libration angles, following constrains are also required to be subjected 0 ≤ T<sup>1</sup> ≤ T1max, 0 ≤ T<sup>2</sup> ≤ T2max, j j θ<sup>1</sup> ≤ θ1max, j j θ<sup>2</sup> ≤ θ2max, Lc j j ≤ LcLimit L\_ <sup>1</sup><sup>m</sup> ≤ L\_ <sup>1</sup> ≤ L\_ <sup>1</sup>M, L\_ cm ≤ L\_ <sup>c</sup> ≤ L\_ cM, where T1max and T2max are the upper bounds of the tension control inputs T<sup>1</sup> and T2, respectively. θ1max, θ2max and LcLimit are the magnitudes of libration angles θ1, θ<sup>2</sup> and the maximum available length scale of Lc, respectively. L\_ <sup>1</sup><sup>m</sup> and L\_ <sup>1</sup><sup>M</sup> are the lower and upper bounds of climber's moving speed L\_ 1, respectively. L\_ cm and L\_ cM are the lower and upper bounds of end-bodies' moving speed L\_ <sup>c</sup>, respectively. It should be noting that, to avoid the tether slacking, the control tensions are not allowed smaller than zero. Dividing the time interval ti ½ � ; tiþ<sup>1</sup> evenly into n subintervals. The cost function minimization problem for each time interval can be solved by Hermite–Simpson method, due its simplicity and accuracy [17]. Then nonlinear programming problem is to search optimal values for the programming variables that minimize the cost function for each interval shown in (25). The closed-loop optimal tracking control method, is same as that in case 1. Direct transcription methods are routinely implemented with standard nonlinear programming (NLP) software. The sparse sequential quadratic programming software SNOPT is used via a MATLAB-executable file interface.

#### 3.2.2. Results and discussion

The proposed control scheme is used to suppress the libration angles of the PSE system in the ascending process with following system parameters and initial conditions: r = 7100 km, m1 = 500 kg, m<sup>2</sup> = 1000 kg, θ1(0) = θ2(0) = 0, L<sup>0</sup> = 20 km, L1(0) = 19,500 m, Lc (0) = 0, θ\_ <sup>1</sup>ð Þ¼ <sup>0</sup> <sup>θ</sup>\_ <sup>2</sup>ð Þ¼ <sup>0</sup> 0, <sup>L</sup>\_ <sup>1</sup>ð Þ¼� <sup>0</sup> <sup>20</sup> <sup>m</sup>=s, and <sup>L</sup>\_ <sup>c</sup>ð Þ¼ 0 0 for the ascending process. The whole transfer trajectory is divided into 50 intervals. The climber's ascending speed along tether 1 is allowed to be controlled to help suppress the libration angles and keep the states of the system in an acceptable area. The constrains are set as T1max ¼ T2max ¼ 200N, θ1max ¼ θ2max ¼ 0:3 rad, L\_ <sup>1</sup><sup>m</sup> ¼ �<sup>15</sup> <sup>m</sup>=s, <sup>L</sup>\_ <sup>1</sup><sup>M</sup> ¼ �<sup>25</sup> <sup>m</sup>=s, <sup>L</sup>\_ cm ¼ �<sup>10</sup> <sup>m</sup>=<sup>s</sup> and <sup>L</sup>\_ cM ¼ 10 m=s.

The simulation results of this case are shown in Figures 9–11. The climber's open-loop libration angle approaches its upper bound at 850 s. After 850 s θ<sup>1</sup> is kept at 0.3 rad by the end of the ascending period, see the dashed line in Figure 10. Using the closed-loop control, the tracking trajectory of θ<sup>1</sup> matches the open-loop trajectory very well overall, see solid line in Figure 9. A short gap appears between 875 s – 880 s, this is caused by the errors of the model and computation. Figure 9 also shows the changes of the trajectories of θ2.The trajectory of θ<sup>2</sup> obtained by closed-loop control tracks the open-loop trajectory well and reaches 0.1 rad by the end of the ascending period. The closed-loop trajectories of L<sup>1</sup> and Lc are shown in Figure 10. They are the reflections of the control inputs. Both L<sup>1</sup> and Lc show smooth fluctuations between 40s and 350 s. Figure 10 shows the time history of trajectories of L\_ <sup>1</sup> and L\_ c,

Figure 9. Libration angle of θ<sup>1</sup> and θ<sup>2</sup> with variable climber speed.

of the climber is discretized evenly. The optimal trajectory should be found to make the transfer time minimize in each equal tether transferring length. Then the cost function can be

> <sup>θ</sup><sup>1</sup> � <sup>2</sup> <sup>ω</sup> <sup>þ</sup> <sup>θ</sup>\_ <sup>1</sup> L\_

L\_

<sup>L</sup><sup>1</sup> <sup>þ</sup> <sup>2</sup>ωL1θ\_

2 ð Þ <sup>L</sup><sup>0</sup> � <sup>L</sup><sup>1</sup> <sup>þ</sup> Lc <sup>θ</sup>\_

> <sup>2</sup>; L1; L\_ 1

gravitational perturbations and the trigonometric functions are ignored, then they can be simplified following the assumptions: sin θ<sup>j</sup> � θj, cos θ<sup>j</sup> � 1 ð Þ j ¼ 1; 2 . All the errors

To ensure the availability and the suppression of the libration angles, following constrains are also required to be subjected 0 ≤ T<sup>1</sup> ≤ T1max, 0 ≤ T<sup>2</sup> ≤ T2max, j j θ<sup>1</sup> ≤ θ1max, j j θ<sup>2</sup> ≤ θ2max, Lc j j ≤ LcLimit

inputs T<sup>1</sup> and T2, respectively. θ1max, θ2max and LcLimit are the magnitudes of libration angles θ1,

tether slacking, the control tensions are not allowed smaller than zero. Dividing the time interval ti ½ � ; tiþ<sup>1</sup> evenly into n subintervals. The cost function minimization problem for each time interval can be solved by Hermite–Simpson method, due its simplicity and accuracy [17]. Then nonlinear programming problem is to search optimal values for the programming variables that minimize the cost function for each interval shown in (25). The closed-loop optimal tracking control method, is same as that in case 1. Direct transcription methods are routinely implemented with standard nonlinear programming (NLP) software. The sparse sequential quadratic programming software SNOPT is used via a MATLAB-executable file interface.

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,

1, respectively. L\_

L<sup>0</sup> � L<sup>1</sup> þ Lc

<sup>θ</sup><sup>2</sup> � <sup>2</sup> <sup>ω</sup> <sup>þ</sup> <sup>θ</sup>\_ <sup>2</sup>

L1

<sup>c</sup> � <sup>L</sup>\_ 1 

> <sup>1</sup> <sup>þ</sup> <sup>L</sup>1θ\_ 1 2 þ

1

� ð Þ <sup>θ</sup><sup>1</sup> � <sup>θ</sup><sup>2</sup> <sup>T</sup><sup>2</sup> L1m<sup>1</sup>

<sup>þ</sup> ð Þ <sup>θ</sup><sup>1</sup> � <sup>θ</sup><sup>2</sup> <sup>T</sup><sup>1</sup> ð Þ L<sup>0</sup> � L<sup>1</sup> þ Lc m<sup>1</sup>

> T<sup>2</sup> � T<sup>1</sup> m<sup>1</sup>

<sup>2</sup> <sup>þ</sup> <sup>2</sup><sup>ω</sup> <sup>þ</sup> <sup>θ</sup>\_

and i denotes the interval number. In (26) the

cM, where T1max and T2max are the upper bounds of the tension control

cm and L\_

1 L1θ\_

<sup>1</sup><sup>m</sup> and L\_

<sup>c</sup>, respectively. It should be noting that, to avoid the

<sup>1</sup> � <sup>T</sup><sup>2</sup> m<sup>2</sup>

<sup>1</sup><sup>M</sup> are the lower and

cM are the lower and upper

Ji ¼ tf (36)

(37)

(38)

(39)

(40)

rewritten as

108 Optimization Algorithms - Examples

subject to the simplified dynamic equations

<sup>L</sup>€<sup>c</sup> <sup>¼</sup> <sup>3</sup>ω<sup>2</sup>

where <sup>u</sup> <sup>¼</sup> ð Þ <sup>T</sup>1; <sup>T</sup><sup>2</sup> , <sup>x</sup> <sup>¼</sup> <sup>θ</sup>1; <sup>θ</sup>\_

cm ≤ L\_ <sup>c</sup> ≤ L\_

upper bounds of climber's moving speed L\_

bounds of end-bodies' moving speed L\_

3.2.2. Results and discussion

L\_ <sup>1</sup><sup>m</sup> ≤ L\_

<sup>1</sup> ≤ L\_ <sup>1</sup>M, L\_ <sup>θ</sup>€<sup>2</sup> ¼ �3ω<sup>2</sup>

<sup>θ</sup>€<sup>1</sup> ¼ �3ω<sup>2</sup>

<sup>L</sup>€<sup>1</sup> <sup>¼</sup> <sup>3</sup>ω<sup>2</sup>

ð Þþ <sup>L</sup><sup>0</sup> <sup>þ</sup> Lc <sup>2</sup><sup>ω</sup> <sup>þ</sup> <sup>θ</sup>\_

<sup>1</sup>; θ2; θ\_

θ<sup>2</sup> and the maximum available length scale of Lc, respectively. L\_

between simple model and the entire model are regarded as perturbations.

Figure 10. Length and its changing ratio of L1, Lc with variable climber speed.

(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

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

M<sup>11</sup> M<sup>12</sup> 0

0 0

⋱

0 0

<sup>A</sup><sup>j</sup>þ0:<sup>5</sup>A<sup>j</sup>þ<sup>1</sup>

B<sup>n</sup>�0:<sup>5</sup>

χ

B<sup>0</sup>:<sup>5</sup>

Piecewise Parallel Optimal Algorithm http://dx.doi.org/10.5772/intechopen.76625 M<sup>22</sup> 0

0 0 M<sup>33</sup>

, C<sup>3</sup> <sup>¼</sup> <sup>2</sup> γ

<sup>A</sup><sup>j</sup>þ0:<sup>5</sup> � <sup>γ</sup><sup>2</sup>

M<sup>T</sup>

0 0

⋱ ⋱

0 0

<sup>j</sup> ¼ �<sup>E</sup> <sup>þ</sup> <sup>γ</sup>

χ3 <sup>n</sup>�<sup>1</sup> <sup>χ</sup><sup>4</sup> n�1

<sup>A</sup><sup>j</sup>þ0:<sup>5</sup>B<sup>j</sup>þ<sup>1</sup> <sup>þ</sup> <sup>γ</sup>

<sup>A</sup><sup>j</sup>þ<sup>1</sup> <sup>þ</sup> <sup>γ</sup>

<sup>B</sup><sup>j</sup>þ<sup>1</sup>

and the Chinese Scholarship Council Scholarship No. 201606290135.

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

C ¼ ½ � C<sup>1</sup> C<sup>2</sup> C<sup>3</sup> , M ¼

, C<sup>2</sup> ¼

<sup>A</sup><sup>j</sup>þ0:<sup>5</sup>Aj, <sup>χ</sup><sup>2</sup>

<sup>B</sup>j, <sup>χ</sup><sup>4</sup>

χ3 χ<sup>4</sup> 

χ3 χ<sup>4</sup> 

<sup>j</sup> ¼ � <sup>γ</sup><sup>2</sup>

within an allowed limit.

Acknowledgements

A. Appendix

χ1 χ<sup>2</sup> 

χ1 <sup>j</sup> <sup>¼</sup> <sup>E</sup> <sup>þ</sup> <sup>γ</sup>

χ3 <sup>j</sup> <sup>¼</sup> <sup>γ</sup><sup>2</sup>

χ1 χ<sup>2</sup> 

E 0

⋱ ⋱

E

<sup>A</sup><sup>j</sup> <sup>þ</sup> <sup>γ</sup>

<sup>A</sup><sup>j</sup>þ0:<sup>5</sup>B<sup>j</sup> <sup>þ</sup> <sup>γ</sup>

χ1 <sup>n</sup>�<sup>1</sup> <sup>χ</sup><sup>2</sup> n�1

<sup>A</sup><sup>j</sup>þ0:<sup>5</sup> <sup>þ</sup> <sup>γ</sup><sup>2</sup>

C<sup>1</sup> ¼

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

respectively. In the first 140 s the trajectory of L\_ <sup>c</sup> increases continuously until reaches its upper bound. Then, it keeps at 10 m/s by 270 s with some slight fluctuations. From 270 s to 355 s, it reduces continuously to 3 m/s. After that, <sup>L</sup>\_ <sup>c</sup> fluctuates around 3 m/s by the end of the transfer period. As a reflection of the control input, L\_ also shows fluctuation during the transfer phase with obvious small-scale fluctuations appear in the period of 120 s – 260 s and 360 s – 750 s. This impacts during the whole transfer period, the changeable speed of the climber has the ability to help the suppression of the libration angles and states trajectory tracking. This time history of the control inputs is shown in Figure 11 with frequent changes between its lower bound and upper bound.
