Data-Driven Nonlinear Control for Orbit Transfer: From the Earth to Moon and L4/5

Show more

1. Introduction

During the progresses of linear algebra, mathematicians have designed well-developed methods (e.g. Linear Quadratic Operator, LQR) to figure out solutions of optimal control problems. However, it remains a challenge for scientists to figure out the general optimal control framework for non-linear systems. As the new technology emerges, humans lack simple models for a general optimal control framework for non-linear system of nonlinear systems. This motivates scientists to propose data-driven models for non-linear optimal control problem [1] [2] . Koopman eigenfunction is a promising approach to provide a general solution for non-linear optimal control. Koopman eigenfunction control remodels the strongly non-liner system into a linear framework [1] to use well-developed tools for linear system. Combining with the well-developed linear optimal control solver; e.g. LQR, we obtain the Koopman control technique, which is the core idea of this paper.

Four steps exist in the Koopman control designs: the determination of the system dynamics, finding the Koopman eigenfunction, the incorporation of control, and the design of a controller. One can treat Koopman eigenfunction to be the director for non-linear data-driven control to become linear.

Restricted Body Problem has important non-linear component. As a result, when we consider the problem, we should take priority for its non-linear identity, and we extend the application of Koopman eigenfunction control to this problem. Restricted Three-Body Problem has deep engineering implications: It can help humans to determine the suitable time for the satellite to change its orbit. Restricted Three-Body Problem has five Lagrange points where the net force is zero. The second Lagrange point is the place where the space telescope James Webb [3] works in. The Lagrange point 4 and point 5 are a suitable place to build Spatial VLBI [4] , which guarantees the acceptance of much more information from the outer space. Classical three-body problem has perplexed scientists for hundreds of years. The typical three-body problem involves 18 first order differential equations. Through use of conservation equations and calculus, the order can be reduced to 6. It has still not been solved because there are not enough conservation quantities to allow for further simplification. In this paper, we consider a restricted three-body model neglecting the force from the smallest body to other larger bodies (i.e. the satellite-earth-moon system or the comet coming into the solar system).

Our goal is to simulate the satellite orbits and design the optimal orbit transfer for it. The Koopman method is adopted because it is an advanced idea to obtain linear representations of nonlinear dynamical systems and thus provide us an easier way to design non-linear optimal control strategy.

2. Dynamics and Control

As we all know, it is important to control the three body system with the non-linear dynamics. That is the motivation for us to carry on with next procedures: to figure out Koopman eigenfunction and apply control.

In a restricted Three Body Problem, the system is restricted to an x-y plane where the center of two larger mass is located at the origin. Generally, it should be calculated in a 3-dimensional plane. However, we just leave the z-axis out since motion in z-axis is relatively small compared to other x and y coordinates. We locate the relative location of satellite to the whole system (the system containing an earth and moon), and assume the relative position of earth and moon in rotational orbits (in reality they will move along certain obits, but in this case our subject, the satellite, moves with them. we should not make the situation complex).

Now we assume the mass of two larger objects are ${m}_{1}$ and ${m}_{2}$ ( ${m}_{1}>{m}_{2}$ ). The third mass, in this case the satellite, has a negligent mass. The distance of the two larger mass to the origin is ${r}_{1}$ and ${r}_{2}$ , and they satisfy the equation below:

${r}_{1}=\frac{R{m}_{2}}{{m}_{1}+{m}_{2}}$ (1)

${r}_{2}=\frac{R{m}_{1}}{{m}_{1}+{m}_{2}}$ (2)

The two equations can be transformed into:

${r}_{1}/{r}_{2}={m}_{2}/{m}_{1}$ (3)

The larger two mass both move around the center of mass at the same angular velocity. As a result, the distance between these two bodies will not change. In order to be simple, we normalize all distance by R and thus just the distance of the two larger bodies ${m}_{1}$ and ${m}_{2}$ to be 1 ( ${r}_{1}+{r}_{2}=1$ ). By doing this, we can avoid the sophisticated calculations of the Three Body System. Another parameter is also introduced: $\mu $ to simplify the problem in that only the ratio between two larger mass is of interest.

$\mu =\frac{{m}_{2}}{{m}_{1}+{m}_{2}}$ (4)

${r}_{1}=\mu $ (5)

${r}_{2}=1-\mu $ (6)

The term ${p}_{1}$ and ${p}_{2}$ represents the distance vector of the third body relative to the first and second body.

The equations of the motion of the third body must be known. The kinetic energy of the third body equals to:

$T=\frac{1}{2}\times {m}_{3}\times \left({\left(\frac{\text{d}x}{\text{d}t}\right)}^{2}+{\left(\frac{\text{d}y}{\text{d}t}\right)}^{2}\right)$ (7)

The potential equals to:

$V=-\frac{G{m}_{1}{m}_{3}}{{p}_{1}}-\frac{G{m}_{3}{m}_{2}}{{p}_{2}}$ (8)

${p}_{1}$ and ${p}_{2}$ equals to:

${p}_{1}=\sqrt{{\left(x-{x}_{1}\right)}^{2}+{\left(y-{y}_{1}\right)}^{2}}$ (9)

${p}_{2}=\sqrt{\left(x-{x}_{2}\right){}^{2}+{\left(y-{y}_{2}\right)}^{2}}$ (10)

Now we switch to the rotational system to make the two larger bodies stationary in the system. The Langrangian ( $L=T-V$ ) becomes (after switching to the rotational system): (Figure 1)

$L=\frac{1}{2}{m}_{3}\left({x}^{2}+{y}^{2}+2x\omega y-2y\omega x+{\omega}^{2}\left({x}^{2}+{y}^{2}\right)\right)-\frac{G{m}_{1}{m}_{3}}{{p}_{1}}-\frac{G{m}_{3}{m}_{2}}{{p}_{2}}$ (11)

Using Lagrange equation: (This is equivalent to the Newton’s second law in mechanics).

Figure 1. Three body visualization.

$\frac{\text{d}}{\text{d}t}\frac{\partial L}{\partial \stackrel{\dot{}}{q}}-\frac{\partial L}{\partial q}=0$ (12)

where q = x, y is the generalized coordinates. Through normalize the setting $\omega =1$ :

${x}^{\u2033}=2{y}^{\prime}+x-\frac{\left(1-\mu \right)\left(x-{x}_{1}\right)}{{p}_{1}^{3}}-\frac{\mu \left(x-{x}_{2}\right)}{{p}_{2}^{3}}$ (13)

$\text{y \u2033}=2{x}^{\prime}+y-\frac{\left(1-\mu \right)\left(y-{y}_{1}\right)}{{p}_{1}^{3}}-\frac{\mu \left(y-{y}_{2}\right)}{{p}_{2}^{3}}$ (14)

$\frac{{\text{d}}^{2}x}{\text{d}{t}^{2}}\left(\frac{\text{d}x}{\text{d}t}\right)+\frac{{\text{d}}^{2}y}{\text{d}{t}^{2}}\left(\frac{\text{d}y}{\text{d}t}\right)=2\frac{\text{d}U}{\text{d}t}$ (15)

The final equation is:

$\left[\begin{array}{l}{\stackrel{\dot{}}{x}}_{1}\\ {\stackrel{\dot{}}{x}}_{2}\\ {\stackrel{\dot{}}{x}}_{3}\\ {\stackrel{\dot{}}{x}}_{4}\end{array}\right]=\left[\begin{array}{c}{x}_{2}\\ 2{x}_{4}+{x}_{1}-\frac{\left(1-\mu \right)\left({x}_{1}+\mu \right)}{{p}_{1}^{3}}-\frac{\mu \left({x}_{1}-1+\mu \right)}{{p}_{2}^{3}}\\ {x}_{4}\\ -2{x}_{2}+{x}_{3}-\frac{\left(1-\mu \right){x}_{3}}{{p}_{1}^{3}}-\frac{\mu \left({x}_{1}-1+\mu \right)}{{p}_{2}^{3}}\end{array}\right]:=f\left(x\right)$ (16)

Finally, we get the original system dynamics that describes the motion and the position of the third body.

Now it is the time to determine the eigenfunction in order to change the dynamics to a linear one to incorporate control.

At first, let me introduce how Koopman eigenfunction works [5] : We have a non-linear ordinary dynamics: (the dynamics are given by ${x}_{k+1}={F}_{k}$ ); F is the transitional pathway from one x to another.

$\frac{\text{d}}{\text{d}t}x\left(t\right)=f\left(x\right)$ (17)

[6] introduced an infinite dimensional linear operator, given by ${K}_{t}$ that acts to advance all measurement functions g with the flow ${K}_{t}g=g\left(F\right)$ . ${K}_{t}$ is the transitional pathway between two y, which is equivalent to g(x). The Koopman operator advances the measurements linearly:

$g\left({x}_{k+1}\right)={K}_{t}g\left({x}_{k}\right)$ (18)

As a result, it is easy for us to deduce: (Figure 2)

$\frac{\text{d}}{\text{d}t}g\left(x\right)={K}_{t}g$ (19)

K is the infinitesimal generator of the Koopman operator. We can rewrite the formula:

${K}_{t}\varphi \left(x\right)=\lambda \varphi \left(x\right)$ (20)

to:

$A\nu =\lambda \nu $ (21)

where $\lambda $ is the eigenvalue and v is the eigenfuction. A Koopman eigenfunction $\varphi \left(x\right)$ corresponding to eigenvalue $\lambda $ satisfies

$\lambda \varphi \left(x\right)=\varphi \left(F\left(x\right)\right)$ (22)

and

$\frac{\text{d}}{\text{d}t}\varphi \left(x\right)=\varphi \left(F\left(x\right)\right)$ (23)

Applying the chain rule to (21) yields:

$\frac{\text{d}}{\text{d}t}\varphi \left(x\right)={\nabla}_{x}\varphi \left(x\right)\times f\left(x\right)$ (24)

Combined with (21) and (19), this yields: (Figure 3)

${\nabla}_{x}\varphi \left(x\right)\times f\left(x\right)=\lambda \varphi \left(x\right)$ (25)

The Koopman eigenfunction can be found by using the above method.

Now, the next step is to find the eigenfunction in the Restricted Body Problem.

Figure 2. Koopman method.

Figure 3. Koopman method.

The Jacobi integral is chosen to be our eigenfunction here because among all the Koopman eigenfunction, it is a constant and thus relatively easy to be controlled. We first consider the possibility of using data-driven as the method to find a Koopman eigenfunction. However, the results are null. Most of the Koopman eigenfunctions have strange shapes and are hard to be characterized by any function in algebra. To our surprise, the failure of using data-driven is not totally of no use. Among all the Koopman eigenfunctions, I find a Koopman eigenfunction which has unique coefficients similar to Jacobi constant, as shown in the expression of the Equation (44):

$C={n}^{2}\left({x}^{2}+{y}^{2}\right)+2\left(\frac{\mu}{{r}_{1}}+\frac{\mu}{{r}_{2}}\right)-\left({\stackrel{\dot{}}{x}}^{2}+{\stackrel{\dot{}}{y}}^{2}+{\stackrel{\dot{}}{z}}^{2}\right)$ (26)

The coefficient for all two variables and their derivatives satisfies a certain ratio, which is similar in the eigenfunction detected by data-driven. The data-driven method [6] is a comprehensive method to find the eigenfunctions. It first sets a library function of candidate function $\theta \left(x\right)$ , and the library function is chosen to be large enough. As a result, the Koopman eigenfunctions can be well approximated, as shown in Equation (45):

$\varphi \left(x\right)\approx {\displaystyle \underset{k=1}{\overset{P}{\sum}}{\theta}_{k}\left(x\right){\xi}_{k}}=\theta \left(x\right)\xi $ (27)

For a specific value of $\lambda $ , the Koopman may be based on data, yielding:

$\left(\lambda \theta \left(x\right)-r\left(x,dot\left(x\right)\right)\right)\xi =0$ (28)

Therefore, it not only proves that our original method (as indicated in Equation (33) and Equation (22)) is correct, but also guarantees that the Koopman eigenfunction we find is the simplest and best among all. The function U, corresponding to the negative effective potential, is useful in further analysis of the system.

$U=\frac{1}{2}\left({x}^{2}+{y}^{2}\right)+\frac{{\mu}_{1}}{{p}_{1}}+\frac{{\mu}_{2}}{{p}_{2}}$ (29)

When taking the partial derivatives with respect to x and y respectively, it is apparent that Equations (25) can be rewritten as:

$\frac{{\text{d}}^{2}x}{\text{d}t}-2\stackrel{\dot{}}{y}=\frac{\partial U}{\partial x}$ (30)

$\frac{{\text{d}}^{2}y}{\text{d}t}+2\stackrel{\dot{}}{x}=\frac{\partial U}{\partial y}$ (31)

Multiplying Equation (26) by $\frac{\text{d}x}{\text{d}t}$ and Equation (28) by $\frac{\text{d}y}{\text{d}t}$ and adding the two, terms cancel to give:

$\frac{{\text{d}}^{2}x}{\text{d}{t}^{2}}\left(\frac{\text{d}x}{\text{d}t}\right)+\frac{{\text{d}}^{2}y}{\text{d}{t}^{2}}\left(\frac{\text{d}y}{\text{d}t}\right)=\frac{\text{d}x}{\text{d}t}\frac{\partial U}{\partial x}+\frac{\text{d}y}{\text{d}t}\frac{\partial U}{\partial y}$ (32)

U depends on x and y, and after taking integration respect to time, the Equation (29) becomes:

${\left(\frac{\text{d}x}{\text{d}t}\right)}^{2}+{\left(\frac{\text{d}y}{\text{d}t}\right)}^{2}=2U-C$ (33)

C is the Jacobi’s constant. The left side of the equation is the velocity of the third body squared, so it can be rewritten as:

${V}^{2}=2U-C$ (34)

Due to C is a constant, we can treat it as an eigenfunction with an eigenvalue $\lambda $ equals to 0:

$\frac{\text{d}C}{\text{d}t}=0$ (35)

Now we can incorporate control, but let me first introduce how to incorporate control. Consider a control-affine system:

$\frac{\text{d}}{\text{d}t}x\left(t\right)=f\left(x\right)+Bu$ (36)

The linear system in Koopman eigenfunction coordinates becomes:

$\frac{\text{d}}{\text{d}t}\varphi \left(x\right)={\nabla}_{x}\varphi \left(x\right)\cdot \left(f\left(x\right)+Bu\right)$ (37)

$=\lambda \varphi \left(x\right)+{\nabla}_{x}\varphi \left(x\right)\cdot Bu$ (38)

The control input enters the dynamics of $\varphi $ via an additional term leading to a control-affine system, which is linear in $\varphi $ .

$\frac{\text{d}x}{\text{d}t}=A\left(x\right)+Bu$ (39)

That Equation (36) is non-linear. A(x) is the dynamics of Restricted Three Body Problem above. u is added to the acceleration function of the satellite. At the end, control (the symbol u) is added as a force in the eigenfunction to ascertain the best orbit, since we add it in the function to determine accelerations in both axes. Consequently, the obit is designed in a way that guarantees the least thrust to drive the satellite. However, the control equation is linear for eigenfunction C, Jacobi’s constant.

$\frac{\text{d}\varphi \left(x\right)}{\text{d}t}=0\times \varphi \left(x\right)+{\nabla}_{x}\varphi \left(x\right)\cdot Bu$ (40)

In this equation, u is affected by eigenfunction C, and at the same time C is always changed by u according to u = −KC:

$\frac{\text{d}x}{\text{d}t}=A\left(x\right)-BKC$ (41)

Now we introduce a method to calculate the solution of linear optimal control problem, LQR. The LQR optimal design refers to the design of the state feedback controller K to make the quadratic objective function J take the minimum value.

For a continuous-time linear system in the Restricted Three Body problem, described by:

$\stackrel{\dot{}}{x}=A\left(x\right)+Bu$ (42)

has the cost function J:

$J={x}^{\text{T}}\left({t}_{1}\right)F\left({t}_{1}\right)x\left({t}_{1}\right)+{\displaystyle {\int}_{{t}_{0}}^{{t}_{1}}\left({x}^{\text{T}}Qx+{u}^{\text{T}}Ru+2{x}^{\text{T}}Nu\right)\text{d}t}$ (43)

The feedback control law that minimizes the value of the cost is:

$u=-KC$ (44)

where K is determined by the weight matrix Q and R uniquely, so the choice of Q and R is particularly important, given by

$K={R}^{-1}\left({B}^{\text{T}}P\left(t\right)+{N}^{\text{T}}\right)$ (45)

Now we are not trying to study Lagrange points individually. We plot the zero velocity curve of satellites which go through these those points, and these curves limit possible orbits of satellite. The energy level of these curves marks the barrier of satellite to go across on curve. Such curves are the boundaries for satellites in their closed field. We plot the zero velocity curve of satellites which go through these those points, and these curves limit possible orbits of satellite. The energy level of these curves marks the barrier of satellite to go across on curve. Such curves are indications about the moment that satellite can change its orbit (Figure 4).

Figure 4. Orbit transfer Lagrange point 4 and 5.

3. Results

The first important result is the orbit change. The initial condition of the simulation is shown (Figure 5): the satellite starts from a point left of the Lagrange point 1. The coordinate is (0.83689291982029 - 0.01, 0) following the coordinates shown in the Three Body Problem visualization Figure 1. T he satellite starts from left of the first Lagrange point. The blue line is the orbit when no thrusts are added to the satellite. H is kept at −3.19 as shown in Figure 6. The green line is the orbit when an acceleration is added, and the satellite possess higher potential energy. H is increased to −3.17 corresponding to the period between 10 and 15 seconds as indicated in Figure 6. It starts to go out of the limitation of the zero velocity line that crosses the first Lagrange point. Thus, it can let the satellite to enter into the moon’s orbit from the earth’s. Then the red line indicates a decrease in speed, and the satellite starts to enter into the sphere limited by zero velocity curve that crosses the first Lagrange point again. H is returned back to −3.19 corresponding to the period between 15 seconds and 30 seconds as indicated in Figure 6.

Figure 5. Orbit transfer: blue dashed line: the satellite is kept moving around the earth; green solid line: the satellite is transferring from the earth to the moon; red dashed dot line: the satellite is kept moving around the moon.

Figure 6. Time history of Jacobean integral.

The thrust is changed by alternating the value of H:

H is the Jacobi constant C, as mentioned above. By changing H, we are actually changing the total energy of the three body system; i.e. kinetic energy and potential energy. Although our designed controller is linear with respect to Koopman eigenfunction $\varphi $ Equation (39), the controller is nonlinear for original state variable x; i.e., the velocity and position.

We are going to extend the current framework to let the satellite pass through ${L}_{4}$ Lagrange points. As far as the scientists concerned, the satellites which enter around the Lagrange point 4 and point 5 has quite deep scientific applications. Scientists have long concerned to build Spatial VLBI [4] at these two points, which guarantees the acceptance of much more information from the outer space.

The process of transmitting a satellite to these two points is the same. By changing the total energy H, we increase the potential energy and let the satellite enter into a high potential energy orbit in Figure 7. However, we face a big trouble: we cannot control the kinetic energy of the satellite in Figure 8. Therefore, the speed of the satellite remains very high when it approaches the two points, which make it hard to stop. We solve this problem by changing the R in the Equation (43) to decrease the rate of acceleration. Therefore, the kinetic energy can be smaller than that in the previous time. The potential energy will be greater when the total energy is constant. Also, the H is increased to reach the potential energy level around Lagrange point 4 and 5. These results indicate that Koopman eigenfunction control is a promising framework to perform the orbit transfer. We design a non-linear controller using Jacobi integral to control the total energy of the system and obtain reference from zero velocity curve to design our control strategy.

The satellite also starts from the same point, a bit left from the Lagrange point 1 (0.83689291982029 - 0.01, 0).

Figure 7. Time history of Jacobean integral to Lagrange point 4 and 5.

Figure 8. Orbit Transfer Lagrange point 4 and 5; the satellite is sent around to the orbit close to Lagrange 4 and 5; Orbit transfer: blue dashed line: the satellite is kept moving around the earth; green solid line: the satellite is transferring from the earth to the moon; red dashed dot line: the satellite is kept moving around the moon; black dot line: the satellite is moving around Lagrange point 4 and 5.

4. Conclusion

As the sum of all, we figure out the best way to change orbit from the earth to the moon. During the orbit change, the zero velocity curve serves as a indication of the potential energy level that the satellite must process to travel across. By controlling the total energy of the satellite, we guarantee that it processes the amount of required potential energy. For further extension of the results, we also use the same strategies to design the optimal path to travel to L4 and L5 Lagrange points. The whole system takes places in Restricted Three-Body Problem, and the governing equation is highly non-linear. Due to the non-linear character of this system, the data-driven method proposed by Koopman is used to design the eigenfunction of the whole non-linear Three-Body System that can help the system become linear. A library of candidate function is created, and then the sparse regression is adopted to select the optimal Koopman eigenfunction. Then the LQR, linear operator, is adopted to design the optimal path by solving the transformed linear system of Three-Body Problem. The objective of LQR is to make the cost function J take the minimum value.

References

[1] Kaiser, E., Kutz, J.N. and Brunton, S.L. (2018) Data-Driven Discovery of Koopman Eigen Functions for Control. arXiv preprint arXiv:1707.01146

[2] Brunton, S.L., Proctor, J.L. and Kutz, J.N. (2016) Discovering Governing Equations from Data by Sparse Identification of Nonlinear Dynamical Systems. Proceedings of the National Academy of Sciences of the USA, 113, 3932-3937.

https://doi.org/10.1073/pnas.1517384113

[3] Gardner, J.P., Mather, J.C., Clampin, M., Doyon, R., Greenhouse, M.A., Hammel, H.B., Hutchings, J.B., Jakobsen, P., Lilly, S.J., Long, K.S., et al. (2006) The James Webb Space Telescope. Space Science Reviews, 123, 485-606.

https://doi.org/10.1007/s11214-006-8315-7

[4] Benson, J.M. and Mutel, R. (1979) Multibaseline vlbi Observations of the 1612 mhz oh Masers toward nml cygni and vy canis majoris. The Astrophysical Journal, 233, 119-126.

https://doi.org/10.1086/157372

[5] Brunton, S.L., Brunton, B.W., Proctor, J.L. and Kutz, J.N. (2016) Koopman Invariant Subspaces and Finite Linear Representations of Nonlinear Dynamical Systems for Control. PLoS ONE, 11, e0150171.

https://doi.org/10.1371/journal.pone.0150171

[6] Koopman, B.O. (1931) Hamiltonian Systems and Transformation in Hilbert Space. Proceedings of the National Academy of Sciences of the USA, 17, 315-318.

https://doi.org/10.1073/pnas.17.5.315