Optimal control of a rectilinear motion of a rocket

In this work, we have modelled the problem of maximizing the velocity of a rocket moving with a rectilinear motion by a linear optimal control problem, where the control represents the action of the pilot on the rocket. In order to solve the obtained model, we applied both analytical and numerical methods. The analytical solution is calculated using the Pontryagin maximum principle while the approximate solution of the problem is found using the shooting method as well as two techniques of discretization: the technique using the Cauchy formula and the one using the Euler formula. In order to compare the different methods, we developed an implementation with MATLAB and presented some simulation results.

The traditional approach to solve optimal control problems is the Pontryagin maximum principle. This principle, proposed by L.S. Pontryagin in 1956 [12], generalizes the Euler-Lagrange equations of the calculus of variations and gives a necessary condition of optimality. The control theory is also applied in different mathematical disciplines: optimal control of partial differential equations, stochastic control theory, game theory, etc. The optimal control problem can be solved analytically using the Pontryagin maximum principle, however when the problem is difficult, we must use numerical methods in order to find an approximate solution.
There exist several numerical methods for solving optimal control problems [13]: indirect methods based on the calculus of variations and direct methods based on discretization techniques and optimization. The indirect shooting method [14] is known for its efficiency and accuracy and it is successfully applied to solve practical problems [4,6]. However, when the problem is difficult, direct discretization methods can be used for finding an approximate solution. In [15], the Cauchy formula for solving systems of linear differential equations is used to transform the linear optimal control problem into an equivalent problem, then discretization is performed to get a linear programming problem which is solved using the adaptive method. Recently in [16], the Cauchy formula is used with discretization to transform the original linear optimal control problem into a linear optimization problem which is solved with the hybrid direction algorithm proposed in [17]. In [18], the Euler discretization formula 282 OPTIMAL CONTROL OF A RECTILINEAR MOTION OF A ROCKET is first applied for the initial linear optimal control problem, then the obtained optimization problem is solved with the fmincon function of Matlab. In [19], the Euler formula is used for the discretization of an optimal control problem corresponding to a viral marketing model and the obtained optimization problem was solved with different optimization solvers.
In this work, we have proposed a linear optimal control model with free final time for maximizing the velocity of a rocket moving with a rectilinear motion from an initial position to a final one. In order to solve the obtained model, we applied both analytical and numerical methods. The analytical solution is calculated using the Pontryagin maximum principle, while the approximate solution of the problem is found using the shooting method and two discretization techniques: the Cauchy discretization technique and the Euler discretization one. The obtained linear programming problems are efficiently solved with an interior-point method implemented in Matlab.
This article is structured as follows. In Section 2, we build the mathematical optimal control model corresponding to the problem of maximizing the velocity of a rocket which moves with a rectilinear motion. In Section 3, we fix the value of the final time by solving an other problem called a minimal time optimal control problem. In Section 4, we solve the problem of maximizing the velocity of the rocket with four methods: the analytical method, the shooting method and two techniques of discretization: the technique using the Cauchy formula and the one using the Euler formula. Then, we compare the obtained results. Finally, we conclude our paper and give some future works.

Problem modelling
Consider a rocket with mass m and velocity v(t) at time t ∈ [0, T ], which moves from an initial height h(0) to a final height h(T ) with a rectilinear motion. Using the laws of physics for bodies with rectilinear motion [1,2], we obtain the following equation: where h(t) and T p (t) represent respectively the travelled distance and the thrust of the rocket at time t ∈ [0, T ], g is the gravitational acceleration. Let us put w(t) = Tp(t) m and dh dt (t) = v(t). Hence, we get If the pilot does not apply any command on the rocket, then the thrust w(t) satisfies the following equation: where α > 0 is the inverse of the response time of the rocket motors to the pilot application. In order to control the movement of the rocket, the pilot applies a force u(t), t ∈ [0, T ]. The ordinary differential equation (2) becomes: The objective is therefore to maximize the velocity of the rocket which moves from the initial point h(0) = h 0 , to a given height h 1 . This leads us to solve the following optimal control problem: This problem is written in the following matrix form: where The symbol " ′ " designates the transposition operation. The Kalman matrix is given by: The rank of matrix K is rank(K) = 3, so the system (5) is controllable [6].

Minimal time problem
Similarly to [6,18], we compute the value of the final time T , which will be used in problem (4), by solving the following minimal time problem: In order to solve Problem (6), we use the Pontryagin maximum principle. Let us denote by , p w (t)) the adjoint vector associated to the problem (6), whose Hamiltonian is given as follows: Since we maximize the Hamiltonian, we set p 0 = −1 [6]. The adjoint vectors are the solutions of the system corresponding to the following Euler-Lagrange equations: From the system (8), we deduce that where The dynamical system is autonomous, so the Hamiltonian is constant for all t ∈ [0, T ]. Therefore, we will have Since the final time T is free, according to the condition of transversality, we get: From relationships (12) and (13), we obtain The Hamiltonian's maximum is given by: Hence, the control maximizing the Hamiltonian is: According to (14) and (15), we will have u * (0) = −1.
If we fix the response time of the motors at 0.7s, we find α = 1.42 and we use the following data to solve our problem: h 0 = 0m, h 1 = 600m and g = 9.80665m.s −2 . The goal is therefore to determine the minimum time, optimum of the problem (6). By using the shooting method, we have obtained the results shown in Figure 1.
The shooting method is based on the Pontryagin Maximum Principle [6]. It consists in finding a zero of the shooting function associated with the minimal time problem. It is a fast, high-precision method that does not require assumptions about the control structure. The shooting method consists in three main steps [3,18]:  From the graphs of Figure 1, we deduce that T min = 3.43s. Knowing this minimal time, our final time should always be chosen slightly larger than T min [18]. So, we choose T = 3, 5 s.

Solving the problem of the maximum velocity of a rocket
In this section, taking as final time, the time T found in the previous section, i.e., T = 3.5s. We first solve the problem (4) analytically; then we proceed to the numerical resolution with three methods: shooting method, Cauchy discretization method and Euler discretization method.
In order to solve the problem (4), we first apply the Pontryagin maximum principle. The Hamiltonian of the problem (4) is given for t ∈ [0, T ] by: The adjoint vectors are solutions of the following system: From the system (17), we deduce that The Hamiltonian's maximum is given by: The control which maximizes the Hamiltonian is The vector x(T ) must satisfy the constraint Since the final time T is fixed in problem (4), we get the following transversality conditions: where g(x(T )) = −v(T ) and b ∈ R.
From which, we obtain the following boundary conditions: By using the boundary conditions (24) and the relationships (18), we get:

Analytical solution
Since the optimal control is equal to the opposite sign of the adjoint vector p w (t), then the commutation times are given by the roots of the equation p w (t) = 0. Since A is a matrix of order 3 and all its eigenvalues are real, the problem has at most a commutation time t c < T [6].
In order to choose the optimal strategy, consider the following possible strategies:

Strategy 1:
We have, Using the initial conditions, we obtain: For t = T , we will have h(T ) = 549.0243 < 600. Therefore Strategy 1 is not feasible. Strategy 2:

OPTIMAL CONTROL OF A RECTILINEAR MOTION OF A ROCKET
We have · w(t) = α[w(t) + 1], which implies: Using the initial conditions, we obtain: For t = T , we will have: h(T ) = 673.7083 > 600, so Strategy 2 is not feasible.

Strategy 3:
The equations of the first set of trajectories, when The equations of the second set of trajectories, when u(t) = −1, t ∈ [t c , T ], are given by: At the intersection point of the two sets, we obtain: This last condition leads to the following equation: The The equations of the second set of trajectories, when u(t) = 1, t ∈ [t c , T ], are: (34) At the intersection point of the two sets, we obtain: So let's calculate t c such that h(T ) = 600m. This last condition yields the following equation: The numerical solution of this last equation is given by: t c = 0.331 s. According to the system (34), the objective value of the problem (4) is equal to: Note that the objective value given by Strategy 3 is higher than the one found by Strategy 4, which is of the bang-bang type. The theoretical results are plotted in Figures 2 and 3.

Figure 2. t −→ u(t) and t −→ pw(t)
We have therefore determined the optimal trajectory that maximizes v(T ) and satisfies the boundary condition h(T ) = 600m, which is given by

Numerical resolution by the shooting method
The Pontryagin maximum principle leads us to the following boundary value problem: x 1 (0) = 0, x 2 (0) = 0, x 3 (0) = g, where and We construct the following shooting function: Hence, the problem (37) is equivalent to the following problem: We determine the numerical solution of the nonlinear system G(p(0), p(T )) = 0 using Newton method. By implementing this method with Matlab, we find the results plotted in Figure 4. These results show that the commutation time is t c = 0.557s, the maximum velocity is v * (T ) = 940.8949 m.s −1 . The execution time of the shooting method is t = 1.9781s. Let us remark that the shooting method have found the same results as those of the analytical method in a short CPU time. This shows that the shooting method is fast and gives accurate results for the studied problem.

Numerical resolution by the Cauchy discretization method
For a number of discretization subintervals N chosen in advance, the discretization step is θ = T N . The solution of the dynamical system of the problem (4) is given by: where F (t) is the resolvent, which is the solution of the following system: Using this last solution, the initial problem (4) takes the following form:  where We compute the numbers C j and X j as follows: Hence, we obtain the following linear programming problem: We have solved this linear program for different values of N with the interior-point method implemented in MATLAB2009b executed with a PC having a Core i5 microprocessor, 2.40Ghz and a RAM of 6 GO. The execution

Numerical resolution by the Euler discretization method
For a number of subintervals N chosen in advance, we will have the step of discretization θ = T N and the following times: The application of the Euler discretization scheme for solving boundary value problems gives us the following linear programming problem: We also solved the linear program (43) with the MATLAB2009b's interior-point solver. The obtained results (execution time of the interior-point method T , number of iterations N it and the maximum speed v(T )) are presented in Table 2.
From Table 2, we can remark that the Euler discretization method is very fast, however even with a large discretization step, it can not reach the desired accuracy.

Numerical comparison
The optimal speeds found by the four methods are almost similar, however the execution times of these methods are quite different. Note that the shooting method gave an optimal speed with a great accuracy and required a very short execution time (1.97 s). For N = 900, the Cauchy's discretization method gave the maximum velocity with an accuracy of 4 decimals with an execution time of 506.38 s, while the Euler discretization method gave for N = 3000 an optimal speed with an error equal to 0.34 and an execution time of 1.55 s. This shows that:   • The shooting method is highly accurate and fast. • The Cauchy discretization method is highly accurate but it is slow. • The Euler discretization method is very fast but less accurate than the other methods.

Conclusion
In this work, we have solved a practical problem arising in the aerospace field by formulating it as a linear optimal control problem. For the theoretical resolution, the principle of Pontryagin maximum has been used which gives a necessary condition of optimality. Numerically, the considered problem is solved with the shooting method and two techniques of discretization (technique using the Cauchy formula and the one using the Euler formula). The obtained optimal solutions are almost similar, but a large difference in the execution time of the three numerical methods has been observed. In the future, we will introduce a new parameter (heel angle) and other constraints, and then we will model the motion of the rocket as a nonlinear optimal control problem.