Fast convex optimization of vehicle trajectories based on improved trust domain

The re-entry process of lift vehicles faces a variety of complex constraints and path limitations, so fast and accurate trajectory optimization for them is an important topic. This paper proposes a fast trajectory generation method based on improved sequential convex optimization, which simplifies the equations of motion with energy as the independent variable, takes the sine and cosine function of the bank angle as the control variable, and convexities multiple constraints to obtain a computational model that satisfies the optimization conditions. On this basis, the "jitter" phenomenon, which is easy to occur in the optimization process, is improved, and the method of variable trust domain and weight function is proposed to reduce the jump problem in the iterative process, and also to improve the terminal accuracy requirement. Finally, numerical simulations show that the algorithm can efficiently and rapidly solve the re-entry trajectory optimization problem under the inclusion of no-fly zone.


Introduction
During the return to the atmosphere, the lift vehicle can make full use of its aerodynamic performance and manoeuvrability to accomplish a wide range, long distance, and high Mach number flight. During flight, it will face complex external aerodynamic environment changes: nonlinear standing heat flow conditions, vehicle overload, and dynamic pressure limitations, as well as flight path constraints (e.g., no-fly zones).Therefore, the trajectory optimization for vehicle re-entry segment has become a hot research topic in recent years, and many effective methods have been developed, such as direct method, indirect method, exact solution method and modern heuristic algorithms. At present, a lot of research and improvement of the above methods are mainly focused on improving the optimization efficiency and shortening the planning time in order to achieve solutions that can be applied to the optimization of real-time trajectories on future bombs [1].
For pseudospectral methods, which mainly include Gauss, Legendre and Redau. And most of the current solutions are adaptations of the above three methods [2].And although modern heuristic algorithms such as genetic algorithms and particle swarm algorithms have unique advantages in terms of iteration speed, the complex constraints on the re-entry process can be less accurate and easily fall into local minima in the later stages of computation.
In contrast, convex optimization, as an optimization algorithm that has gradually emerged in recent years, is insensitive to initial values and has the advantage that the local optimal solution is the global optimal solution when the problem satisfies the model of convex optimization, which can effectively avoid the problem of falling into local minimum value. Moreover, with the improvement of simulation methods, large-scale convex optimization problems can be solved optimally in a short time, and this method is gradually being applied. its pioneering application in astronautics, with Acikmese [3] and others at JPL Laboratory in the United States pioneering the application of convex optimization methods  [4], space orbit transfers [5], extra-terrestrial object landing, and reusable vehicle trajectory optimization scenarios [6][7]. Xiang Zhou [8] et al. proposed the command inverse solution method nested in the convex optimization model, which allows the computational accuracy to be improved and completes the fast planning of re-entrant 3D trajectory profiles. Ben Yang [9] used a bsample curve discrete control volume and a backtracking linear search method to improve the stability and smoothness of the numerical simulation for the "pseudo infeasibility" problem in the re-entry vehicle trajectory optimization process. Jinbo Wang [10] and Lin Ma [11] combined sequential convex optimization and lossless convexification methods to develop a series of studies for application in reusable rocket trajectory optimization and powered soft landing problems.
In this paper, based on the three-degree-of-freedom normalized re-entry model of a lift vehicle, we use a sequential convexity and relaxation method based on the variable-trust domain to convexity the constraints with the energy "e" as the independent variable. The trust region range and the size of the weight function are automatically adjusted to reduce the oscillation problem in the iterations. Finally, the convex problem is discretized and then numerical simulations are performed. The feasibility and computational efficiency of the algorithm are verified by analysing the optimization results.

Re-entry motion model
The re-entry section of the lift vehicle completes its flight under gravity and aerodynamic forces only. In order to simplify the optimization number, energy is chosen as the independent variable, and the influence of the angular velocity of the Earth's rotation is considered, and the three-degree-of-freedom dimensionless motion model of the re-entry section is established as follows.
Where is the radial distance between the centre of mass of the aircraft and the centre of the earth, normalized with the radius of the earth 0 as the unit; and are latitude and longitude respectively, and are the Flight-path angle and heading angle, the units are all in radians; The speed unit is normalized by = ( 0 0 ) 0.5 , where 0 = 9.8 ⋅ −1 is the surface gravitational acceleration. The Colombian acceleration term 1 , 2 ,and the implicated acceleration term 3 , 4 ,caused by the rotation of the earth. and are the dimensionless aerodynamic lift acceleration and drag acceleration.
Let the optimization state quantity be = [ ] , and the control quantity be the bank angle. In order to simplify the post-processing and improve the chattering phenomenon in the optimization process, the control quantity u is decomposed into. = [ ] .

Restrictions and Optimization problem description
The terminal equation is constrained to be altitude, longitude and latitude, and the flight path angle is given within a range. In the same way, the heading angle will also have other flight phases (such as the energy management phase or the terminal guidance phase) after the gliding phase, which also needs to be limited to a certain range. So, it is assumed that the terminal constraints are as follows ( ) = * , ( ) = * , ( ) = * , ( ) = * , 3 ̇= (( 0 0 ) 0.5 ) 3.15 ( ) 0.5 3.15 ≤̇ (4) Heat flow ratėunit is ⋅ −2 ,heat flow coefficient = 9.4369 × 10 −5 ,the unit of dynamic pressure ̄ is ⋅ −2 ,overload is aligned with gravitational acceleration 0 .In addition to the above constraints, the constraints controlling the bank angle of the variables are shown below.

Description of the optimization algorithm
The performance metrics, equations of state and constraints of this vehicle are difficult to solve directly by convex optimization. Therefore, it needs to be convexities. In this paper, we mainly use sequential convex optimization methods and relaxation methods to deal with the model. By solving a series of convex optimization sub-problems, the solution of the original problem is gradually approached.

Convexation of motion equations and control variables
According to the state quantity and control quantity selected above, linearize the motion equation group and rewrite it into the following form: The coefficients in can be obtained by derivation and calculation. Since the linearization expansion occurs at , trust region constraints need to ensure the effectiveness of its linearization.
| − | ≤ (8) In the constraint condition (5) for the control variable u, the first inequality constraint satisfies the convexity requirement, but the second constraint is a non-convex quadratic equality constraint. It is planned to use a relaxation method to transform it into an inequality constraint: ‖ ‖ 2 = 1 ⇒ ‖ ‖ 2 ≤ 1 (9) For the process constraint, it is transformed into a logarithmic form with respect to the vector diameter and is integrated to obtain the following form.

Improved trust region algorithm
For the sequential convex optimization algorithm, it transforms the non-convex problem by linearization, but the process needs to satisfy the trust domain index required for linearization. A fixed reliance domain can lead to large oscillation problems in the preliminary optimization. To improve this problem, an exponential decreasing method based on variable trust domain and weight function is designed as follows to make the iterative process approach the optimized value faster. = − 0 = − − (15) k is the number of iterations, a and b are the setting parameters respectively, and suitable values can be selected according to requirements. In this article, a=1.19 and b=-0.36 are set. Similarly, the penalty function factors c and d in the performance index are also set according to the above method.

Problem description
The paper proposes to use a flight model of reusable aircraft (RLV) in the United States [12]. Its mass is 104305kg, the reference area is 391.22m 2 , and its aerodynamic model fitting function is as follows: = −0.041065 + 0.016292 + 0.0002602 2 = 0.080505 − 0.03026 + 0.86495 2 (16) The relevant setting parameters are shown in Table 1 below. The angle of attack scheme uses a given piecewise function angle of attack profile as literature [9]. ) 1500 In order to ensure the validity of the sequence linearization and the accuracy of the final result in the simulation, the trust region and convergence range are given below.
The method solves convex optimization problems insensitive to initial values, and the initial iterative trajectory state quantities can be determined by linear interpolation.

Analysis of calculation results
The simulation in this paper uses the Mosek solver under the MATLAB platform to solve the convex optimization problem. According to the flight results, the flight time is 1374.05 seconds and the flight distance is 5824.06 kilometres. The results of some flight parameters will be displayed below.  Table 2, it can be seen that the simulation flight process satisfies the constraint requirements. Due to the penalty factor added to the performance index at the early stage of optimization, the terminal latitude and longitude have a little error, but meet the index requirements. The changes in the values of the control volume are concentrated in the 10-degree interval, and the entire trajectory is arc-shaped due to the flight overload limitation, and its form corresponds exactly to the change in the bank angle. When approaching the terminal, the corresponding state amount needs to be adjusted to meet the terminal requirements, so the bank angle appears a large change.
(a) Three-dimensional trajectory (b) Flight trajectory Figure 2. Three-dimensional flight trajectory Consider the path constraints of the no-fly zone as ( − ) 2 + ( − ) 2 ≥ 2 .Among them, and are the longitude and latitude of the circle where the no-fly zone is located, and is the radius. It is easy to know that the constraint does not meet the requirements of convex optimization, and then it is relaxed and convex.̄can be obtained by deducing.
2( − ) + 2( − ) ≥ 2 +̄2 (18)  This article establishes two no-fly zone restrictions, the specific locations of which are shown in Table 3 above. The trajectory under the no-fly zone is shown in Figure 2. It seen that the trajectory obtained by iterative calculation can effectively avoid the no-fly zone and still meet the relevant constraints.
Through calculation, as shown in Figure 3 below, by adding the variable trust region, the oscillation amplitude of the curve becomes significantly smaller, which can effectively avoid the simulation optimization failure caused by large-scale chattering during the iteration process. And the optimized endpoint parameters converge to the terminal values more quickly under the weighting function.
(a) Fixed trust region (b) Optimized trust region Figure 3. Partial view of flight trajectory

Conclusion
In this paper, taking the lift aircraft as the research object, the fast trajectory optimization calculation is carried out for its re-entry section, and the improved sequence convex optimization algorithm of variable trust region is used to effectively reduce the chattering phenomenon in the process and make the optimization result more accurate. Compared with other methods, this algorithm has a shorter calculation time and does not require high initial value. It can be roughly estimated under unknown flight trajectories without affecting the final result, so it will have good application value in the future.