Improved successive approximation control for formation flying at libration points of solar-earth system

1.   Introduction

Since the ISEE-3 mission launched in 1978, the detection of translation points has never been interrupted. In the past 30 years, five spacecrafts have been launched, but these spacecrafts are single-satellite missions. Due to the constraints of a single spacecraft (weight, size, fuel, etc.), it is impossible to complete large baseline length or large aperture astronomical observation with high resolution. With the rapid development of small spacecraft in recent years, many spacecrafts work together to complete formation flying tasks, which provides a feasible way to complete high-performance and high-resolution astronomical observation tasks. At present, the international space agency has attached great importance to the formation astronomical observation tasks, such as XUES, MAXIM, SI and Darwin ^[1,2,3].

The deployment of spacecraft formation in translational orbit can achieve interferometry of medium and long baseline, which is of great significance to the study of the origin of the universe. Firstly, the design of formation configuration and the formation preservation based on it are in-volved in the formation of translational spacecraft, and the formation preservation is one of the main difficulties. At present, Hamilton and Folta applied linear Optimal control technology to study formation control near Lissajous orbit in the real ephemeris model. Penin proposed offline planning and online adjustment algorithm for configuration reconstruction of deep space interfer-ometric formation ^[5]. Peng Haijun and Jiang Xin proposed a nonlinear control method for Space-craft Formation Reconfiguration in libration point orbit based on symplectic numerical algo-rithm ^[6,7,8,9]. Chai Runqi, Tsourdos Antonios et al work explores the optimal flight of aero-assisted reentry vehicles during the atmospheric entry flight phase with the consideration of both determin-istic and control chance constraints ^[10,11].

At present, most of the studies have not touched the vicinity of the solar-earth libration points. In this study, the improved successive approximation control method is used to maintain the for-mation flying control at the solar-earth libration points. Compared with the previous control meth-ods (such as LQR control method, SDRE control method, invariant popular control method, etc.), this method can further improve the accuracy of formation flight hold control, but the amount of calculation is large. It is suitable for long-term deep space exploration missions with high for-mation flying. At the same time, with the development of computer science, the computing power will also be enhanced. It is believed that the successive approximation method will get more atten-tion in the field of formation flight maintenance control in the future.

2.   Successive approximation method

Because the dynamic large-scale system has the characteristics of "large scale, complex struc-ture, many factors and function synthesis", the analysis of large-scale system with conventional control theory and method can't achieve good control effect. Aiming at the general nonlinear inter-connected coupled large-scale system and a class of affine nonlinear similar composite large-scale system, a successive approximation method of Optimal control law design is proposed, Firstly, by using this method, the solution of high-order coupled nonlinear two-point boundary value problems is simplified to a series of decoupled linear two-point boundary value problems, and then it is proved that the solutions of the linear two-point boundary value problems uniformly converge to the Optimal control of nonlinear interconnected dynamic large-scale systems.

The successive approximation method is the simplest and systematic method to study the problem of non-linear local analysis based on the solvability of operator equation. The basic idea of this method is: For a given, displayed a Cauchy sequence of defined ele-ments, so that.

Based on the successive approximation method for local analysis of non-linear functions, an approximation method for solving the Optimal control problem of non-linear systems is proposed in reference ^[12]. First, the non-linear term in the system is treated as an external additional dis-turbance to the system. Then, according to the successive approximation method in differential equation theory, the non-linear two-point boundary value problem is transformed into a series of linear non-homogeneous two-point boundary value problems.

Considering Affine Nonlinear Systems

where:

Its cost function is:

Where: $ {Q}_{f} $ and $ Q $ are the weight matrix, R is the error matrix, $ {Q}_{f} $ and $ Q $ are the semi-positive definite constant matrix, R is the positive definite constant matrix.

According to Optimal control theory, the necessary conditions for this Optimal control problem are the boundary value problems of the following two points:

where: $ S = B{R}^{-1}{B}^{T} $, $ F\left(x\right) = \frac{\partial {f}^{T}\left(x\right)}{\partial x} $. The Optimal control law can be described by the following formula:

Among them, the cost function $ J $ is mainly used to find the objective function of the optimal solution. If we want to find the optimal solution of Eq (1), that is to say, that is to say, we need to minimize the value of $ J $. The function of optimal control rate is to describe the functional relationship between controlled state variables and system input signals.

3.   Improved successive approximation method

Construct the following sequence of two-point boundary value problems.

where: $ 0 < t\le {t}_{f} $, $ k = \mathrm{1, 2}, 3\cdots $.

Definition: $ {G}^{\left(k\right)}\left(t\right) = G\left({x}^{\left(k\right)}\right) $

The Eq (8) is a linear homogeneous two-boundary value problem, and its endpoint condition is also linear. It is known that the relationship between $ {\lambda }^{\left(1\right)}\left(t\right) $ and $ {x}^{\left(1\right)}\left(t\right) $ must be linear.

By deriving t on both sides of equation (10) and by the first and second equations in equation group (8), we can obtain:

Considering that $ {x}^{\left(1\right)}\left(t\right)\equiv 0 $ is not valid, Then the matrix differential equation satisfied by $ {P}^{\left(1\right)}\left(t\right) $ can be obtained as follows:

By using the endpoint condition, we can deduce that $ {P}^{\left(k\right)}\left(t\right) $ satisfies the following conditions:

The corresponding Optimal control sequence is:

The closed-loop system sequence of affine nonlinear systems is constructed as follows:

Algorithm 1 Improved successive approximation algorithm for Optimal control approximation design of affine nonlinear systems

Step 1: Given the normal number $ \varepsilon $. Let $ k = 1 $, $ N = 1 $, from Eq (8), we can see that $ {x}^{\left(0\right)}\left(t\right) = 0 $. According to definition (9), $ {G}^{\left(0\right)}\left(t\right) $ and $ {F}^{\left(0\right)}\left(t\right) $ can be obtained.

Step 2: $ {P}^{\left(1\right)}\left(t\right) $ can be obtained from matrix differential equation (12) and its endpoint condition (13). Then, $ {u}^{\left(1\right)}\left(t\right) $ can be obtained from Eq (14). According to Eq (9), $ {G}^{\left(1\right)}\left(t\right) $ and $ {F}^{\left(1\right)}\left(t\right) $ are obtained. Substituting $ {x}^{\left(1\right)}\left(t\right) $ and $ {u}^{\left(1\right)}\left(t\right) $ into cost function (3) to calculate $ {J}^{\left(1\right)} $. Then, let $ k = 2 $;

Step 3: Matrix $ {P}^{\left(k\right)}\left(t\right) $ is obtained by matrix differential equation (12) and its endpoint condition (13). Then, $ {u}^{\left(k\right)}\left(t\right) $ is obtained from (14). $ {x}^{\left(k\right)}\left(t\right) $ is obtained from (15).

Step 4: Substituting $ {u}^{\left(k\right)}\left(t\right) $ and $ {x}^{\left(k\right)}\left(t\right) $ into cost function (3) to obtain $ {J}^{\left(k\right)} $, and the convergence index $ {J}_{e}^{\left(k\right)} $ of performance index $ {J}^{\left(k\right)} $ is calculated according to Eq (16).

Step 5: If $ \left|{J}_{e}^{\left(k\right)}\right| < \varepsilon $, let $ N = k $, then output is $ {u}_{N}\left(t\right) $; Now termination.

Step 6: $ {G}^{\left(k\right)}\left(t\right) $ and $ {F}^{\left(k\right)}\left(t\right) $ can be calculated from definition (9). Let $ k = k+1 $, then go to step 3.

4.   The proof of convergence

Firstly, the following radiation nonlinear systems are considered.

Where: $ f\left(x\right) = G\left(x\right)x\left(t\right) $, $ x\in {\mathbb{R}}^{n} $ is the state vector of the system, $ u\in {\mathbb{R}}^{n} $ is the control vector of the system, $ A\left(t\right) $ is a time-varying matrix of suitable dimension, $ {x}_{0} $ is the initial state vector of the state vector $ x\left(t\right) $.

Assuming that $ {t}_{f} $ is the terminal time of the system and based on the improved successive approximation method, the sequence of state equations for the two-point boundary value problem is:

The vectors sequence of Eq (13) is defined as follows:

Where: $ \varPhi (t, 0) $ is the state transition matrix for $ A\left(t\right) $.

Since $ G\left(x\right) $ is Lipschitz on $ {\mathbb{R}}^{n} $, for any $ {x}^{\left(i\right)}, {x}^{\left(j\right)}\in U $, there is:

The Lipschitz constant of $ f\left(x\right) $ is assumed to be $ \beta $.

Where: $ M $, $ \gamma $ is a positive constant. $ ‖\varPhi \left(\mathrm{0, 0}\right)‖ = ‖I‖ = 1 $, the moral is that, $ M\ge 1 $.

Next, we prove that $ \left\{{x}^{\left(k\right)}\right(t\left)\right\} $ and $ \left\{{u}^{\left(k\right)}\right(t\left)\right\} $ respectively, converge uniformly to the optimal state trajectories $ {x}^{*}\left(t\right) $ and $ {u}^{*}\left(t\right) $ of the system (8) for quadratic performance indices.

The following lemma needs to be proved first.

Lemma 1 Assume that a and B are sufficiently small and satisfy $ M\alpha \gamma {t}_{f} < 1 $, $ \beta {t}_{f} < 1 $.

Then:

1) the system (8) has a unique solution in interval $ [0, {t}_{f}] $;

2) Sequence $ \left\{{x}^{\left(k\right)}\right(t\left)\right\} $ is the only solution that converges uniformly to the system (8).

Prove: (1) If $ f\left(x\right) $ is Lipschitz on $ {\mathbb{R}}^{n} $, and $ \beta {t}_{f} < 1 $, then the affine system (8) has a unique solution in interval $ [0, {t}_{f}] $.

(2) Consider $ \left\{{x}^{\left(k\right)}\left(t\right)\right\} $ as a sequence on $ {C}^{N}[0, {t}_{f}] $. From Eq (14) the following equation can be obtained.

Then,

For the same reason, the following equation can be obtained.

According to Eq (15), the following equation can be obtained.

According to hypothetical conditions, the following equation can be obtained.

Using Mathematical Induction, the following equation can be obtained.

As Eq (16) is known, the right side of the inequality is enlarged.

According to the trigonometric inequality, for any number j, the following equation can be obtained.

When $ k\to \infty $, $ \frac{{M}^{k+1}{\alpha }^{k+1}{\gamma }^{k+2}{t}_{f}^{k+1}}{(k+1)!{\left(1-M\alpha \gamma {t}_{f}\right)}^{k}} $ is a power series tending to zero and $ \mathit{exp}\left(\frac{M\alpha \gamma {t}_{f}}{1-M\alpha \gamma {t}_{f}}\right) $ is a constant. Then:

That is to say, for any number j, the sequence $ \left\{{x}^{\left(k\right)}\right(t\left)\right\} $ is a Cauchy sequence, which is the only solution uniformly convergent to system (8). The proof is complete.

Theorem 2 Assuming that the solution sequence of the constructed two-point boundary value problem is $ \left\{{x}^{\left(k\right)}\right(t\left)\right\} $ and the control sequence is $ \left\{{u}^{\left(k\right)}\right(t\left)\right\} $, they converge uniformly to the optimal state trajectory $ {x}^{*}\left(t\right) $ and the Optimal control theory $ {u}^{*}\left(t\right) $ of the system (8) with respect to the quadratic performance index, respectively.

Prove: According to lemma 1, the sequence $ \left\{{x}^{\left(k\right)}\right(t\left)\right\} $ is uniformly convergent. When $ k\to \infty $, the sequence $ \left\{{x}^{\left(k\right)}\right(t\left)\right\} $ converges uniformly to the optimal trajectory $ {x}^{*}\left(t\right) $ of the system.

It can be seen from the Eqs (8) and (9) that $ \left\{{u}^{\left(k\right)}\right(t\left)\right\} $ is only related to $ \left\{{x}^{\left(k\right)}\right(t\left)\right\} $ and $ Q $. Obviously, $ Q $ is a semi-positive constant matrix, so $ \left\{{u}^{\left(k\right)}\right(t\left)\right\} $ is uniformly convergent. When $ k\to \infty $, $ \left\{{u}^{\left(k\right)}\right(t\left)\right\} $ converges uniformly to the Optimal control theory $ {u}^{*}\left(t\right) $ of the system for the quadratic performance index. The proof is complete.

5.   Establishment of formation flight motion equation at solar-earth libration points

Assume that In the convergent coordinate system, the coordinate origins O are located at the centroid of $ {M}_{1} $, $ {M}_{2} $ and $ {M}_{3} $. Their masses are $ {m}_{1} $, $ {m}_{2} $ and $ {m}_{3} $. The coordinate system schematic diagram of the three-body problem is shown in Figure 1.

According to the dynamic equation of the circular restricted three-body problem in reference 1, it is obtained that:

Let $ \mu = \frac{{m}_{2}+{m}_{3}}{{m}_{1}+{m}_{2}+{m}_{3}} $, $ {\mu }_{0} = \frac{{m}_{3}}{{m}_{1}+{m}_{2}+{m}_{3}} $, $ k = \frac{d}{R} $, where d is the average distance between the earth and the moon, and R is the average distance between the sun and the earth. Then:

Where: $ {\varOmega }_{1 } $ is the classical restricted trisomy problem and $ {\varOmega }_{2} $ is the lunar perturbation force.

Where: $ {n}_{1} $ is the average angular velocity of the moon relative to the earth. From Eq (31), it is known that the conditions for the existence of equilibrium solutions are as follows:

Five libration points can be solved from Eq (3), of which $ {L}_{4} $ and $ {L}_{5} $ are stable points and $ {L}_{1} $, $ {L}_{2} $ and $ {L}_{3} $ are collinear libration points, which belong to unstable points. This paper chooses the libration point $ {L}_{2} $ of the solar-terrestrial system as the research object.

Assuming the control variable of the main spacecraft is $ {u}^{A} = {\left[\begin{array}{ccc}{u}_{x}^{A}& {u}_{y}^{A}& {u}_{z}^{A}\end{array}\right]}^{T} $, the state vector of the master spacecraft is defined as:

Then, the equation of state of the system is:

Similarly, the equation of state of the slave spacecraft system is as follows:

The relative motion diagram of spacecraft in the libration point orbit is shown in Figure 2. The relative motion equation of the spacecraft can be obtained by subtracting equation (36) from Eq (37) as follows:

Where: $ {u}_{x} = {u}_{x}^{B}-{u}_{x}^{A} $, $ {u}_{y} = {u}_{y}^{B}-{u}_{y}^{A} $, $ {u}_{z} = {u}_{z}^{B}-{u}_{z}^{A} $, $ {x}_{i} = {x}_{i}^{B}-{x}_{i}^{A}(i = \mathrm{1, 2}, \cdots, 6) $.

6.   Design of formation flight controller

The above relative motion equation is written in a general form:

The master spacecraft:

The slave spacecraft:

The Optimal control method can be used to design the controller of the master spacecraft separately. The improved successive approximation method is used to design the controller of the slave spacecraft. In order to use the improved successive approximation method, condition $ f\left(0\right) = 0 $ needs to be satisfied. But the following terms in equation of motion (33) will deviate.

As well as the following two terms in the above formula, the Eq (33) is deviated and the system state cannot be zero. For this reason, these effects need to be eliminated.

And,

In this paper, $ {\dot{s}}_{1} = -{\lambda }_{1}{s}_{1} $ and $ {\dot{s}}_{2} = -{\lambda }_{2}{s}_{2} $ are used to eliminate the influence of deviation term. The following equations can be obtained from the equation of relative motion (33):

$ A = \left[\begin{array}{cccccccc}0& 1& 0& 0& 0& 0& 0& 0\\ 1& 0& 0& 2& 0& 0& {a}_{27}& 0\\ 0& 0& 0& 1& 0& 0& 0& 0\\ 0& -2& 1& 0& 0& 0& 0& {a}_{48}\\ 0& 0& 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& {a}_{65}& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& {\lambda }_{1}& 0\\ 0& 0& 0& 0& 0& 0& 0& {\lambda }_{2}\end{array}\right] $, $ B = \left[\begin{array}{ccc}0& 0& 0\\ 1& 0& 0\\ 0& 0& 0\\ 0& 1& 0\\ 0& 0& 0\\ 0& 0& 1\\ 0& 0& 0\\ 0& 0& 0\end{array}\right] $,

$ f\left(x\right) = \left[\begin{array}{c}0\\ {b}_{27}\\ 0\\ {b}_{48}\\ 0\\ {b}_{65}\end{array}\right] $, $ G\left(x\right) = \left[\begin{array}{c}\mathrm{0, 0}, \mathrm{0, 0}, \mathrm{0, 0}\\ {b}_{27}/{x}_{1}, \mathrm{0, 0}, \mathrm{0, 0}, 0\\ \mathrm{0, 0}, \mathrm{0, 0}, \mathrm{0, 0}\\ \mathrm{0, 0}, {b}_{48}/{x}_{3}, \mathrm{0, 0}, 0\\ \mathrm{0, 0}, \mathrm{0, 0}, \mathrm{0, 0}\\ \mathrm{0, 0}, \mathrm{0, 0}, {b}_{65}/{x}_{5}, 0\end{array}\right] $.

Where:

The quadratic performance index function is selected as follows:

Set the control precision to: $ \varepsilon = 0.01 $.

7.   Numerical examples and simulation results

In the numerical example, the master spacecraft in halo orbit of L2 point of solar-earth system is adopted. Halo orbit is a kind of special orbit, which can maintain spacecraft in space without consuming too much fuel. It is an approximate circular orbit formed on the equilibrium point (translation point) when multiple gravitational sources act on spacecraft. The first application of halo orbit was the international comet probe launched in 1978. It marched to L1 point of solar-earth system and stayed there for several years. Because of its particularity, it is usually used as an important observation point in space missions. The specific values are as follows:

$ {x}_{0} = 87028.508409273 $ km, $ {y}_{0} = -24739.51269980 $ km, $ {z}_{0} = -229952.974656271 $ km, $ {\dot{x}}_{0} = -8.985877859 $ m/s, $ {\dot{y}}_{0} = -121.605674977 $ m/s, $ {\dot{z}}_{0} = 9.457953755 $ m/s.

According to the calculation method of the improved successive approximation method, the calculation is carried out according to algorithm 1. After seven iterations, their performance index $ {J}^{\left(k\right)} $ and its convergence index $ {J}_{e}^{\left(k\right)} $ are calculated. The calculation results are shown in Table 1.

It can be seen from Table 1 that performance indicators $ {J}^{\left(k\right)} $ and $ {J}_{e}^{\left(k\right)} $ are convergent, which indicates that performance index $ J $ is uniformly convergent. The tracking error e(t) of the system tends to zero in the optimal way. As for the selection of parameter K, when the performance index $ {J}^{\left(k\right)} $ is basically unchanged, that is, when the difference between the performance index $ {J}^{\left(k\right)} $ of this iteration and the performance index $ {J}^{(k-1)} $ of the previous iteration is close to zero, the value of K is the last iteration.

The relative initial conditions for formation are:

Figure 3 shows the tracking error $ {e}_{x} $, $ {e}_{y} $, $ {e}_{z} $. Where: $ ({e}_{x}, {e}_{y}, {e}_{z}) = ({x}_{r}-x, {y}_{r}-y, {z}_{r}-z) $, and $ {e}_{r} = \sqrt{{e}_{x}^{2}+{e}_{y}^{2}+{e}_{z}^{2}} $. In order to clearly represent the instantaneous response of the error, Figure 4 provides the error history of the first day. The error dropped quickly to near zero in one day. As can be seen from Figure 4, the error after two days fluctuated within a range of 10mm. Figures 5 and 6 show the control history of formation flight. As can be seen from Figure 5, the control force drops rapidly from 4 $ N $ to near zero. From Figure 6, it can be seen that the average control force to keep formation flying is about 250 $ \mu N $.

8.   Conclusions

In this paper, a non-linear dynamic model of formation flying is established by derivation. After proving the convergence of the improved successive approximation, a formation flight controller is designed. According to the improved successive approximation algorithm, the performance index function is iterated for seven times, and good convergence is obtained. Finally, the error and control simulations of the formation flying of the spacecraft with the libration points of the solar-earth system are carried out for two days. The simulation results show that the error decreases to near zero in one day and fluctuates in the range of 10 mm after two days. The control power on the first day dropped quickly from 4 $ N $ to near zero, and the average control force after two days was about 250 $ \mu N $. The results show that the method has the advantages of high control accuracy and short control time in formation keeping control of solar-earth libration points. The disadvantage is that the calculation is complicated. But, with the development of computer technology, it will have a certain practical value in the future research of team maintenance control.

Recently, convex optimization method is becoming more and more popular ^[12]. Convex optimization method is mostly used in artificial intelligence technology. Recently, the author also tries to use convex optimization method to study formation flight keeping control. In order to achieve better control results. Of course, the control effect of the improved successive approximation method is the highest in terms of accuracy. However, the convex optimization method can greatly reduce the amount of computation and obtain better computational efficiency. The research of convex optimization in formation flight control technology will not be described here, but will be described in another paper of the author.

Acknowledgments

The authors acknowledge the National Natural Science Foundation of China (Grant: 61174202), the National Natural Science Foundation of China (Grant: 61502391), the 13th Five-Year Plan of Shaanxi Province in 2020 (Gant: GH20Y649), School level scientific research pro-jects(Gant:33/115040393, Gant:33/115040441).

Conflict of interest

The authors declare that they have no competing interests.