Trajectory Optimization Based on Multi-Interval Mesh Refinement Method

In order to improve the optimization accuracy and convergence rate for trajectory optimization of the air-to-air missile, a multi-interval mesh refinement Radau pseudospectral method was introduced. This method made the mesh endpoints converge to the practical nonsmooth points and decreased the overall collocation points to improve convergence rate and computational efficiency. The trajectory was divided into four phases according to the working time of engine and handover of midcourse and terminal guidance, and then the optimization model was built. The multi-interval mesh refinement Radau pseudospectral method with different collocation points in each mesh interval was used to solve the trajectory optimization model. Moreover, this method was compared with traditional ℎ method. Simulation results show that this method can decrease the dimensionality of nonlinear programming (NLP) problem and therefore improve the efficiency of pseudospectral methods for solving trajectory optimization problems.


Introduction
Beyond visual range (BVR) combat is being the main form of future air combat [1][2][3].Developing the missile with BVR combat capability by improving the existing medium range missile has great practical significance.Considering the existing propulsion capability, we can effectively increase the range by trajectory optimization.
Trajectory optimization has been a hot research topic in the field of aircraft guidance and control.Many scholars have studied it in view of the different background.
Considering the long range of the flight, problems such as structural strength, normal work of the engine, and control stability give rise to strict requirements on the dynamic pressure, overload, control, and other characteristics of the process.Therefore, trajectory optimization of the air-to-air missile can be regarded as a complex nonlinear optimization problem.
With the development of computer technology, the direct method based on parameterization has developed remarkably [4][5][6].This method converts the optimal control problem of continuous time into NLP problem and then solves it by NLP solver, such as the sequential quadratic programming (SQP) [7], which is recognized as the best solving method for NLP problem.
During last decades, pseudospectral methods have been a research hotspot in the domain of trajectory optimization [8,9].In this method, the control and state are firstly discretized at certain collocation points in the time interval of interest, which are then approximated by polynomial interpolation.Finally, the optimal control problem is transformed into NLP problem which can be solved by NLP solver.Different pseudospectral methods select different collocation points and base functions of interpolation.In Huntington et al. [10], three kinds of commonly used pseudospectral methods, Legendre, Radau, and Gauss, were compared.The computational efficiency of the three methods is similar, while the Gauss and Radau pseudospectral methods (RPM) are superior to the Legendre pseudospectral method in the aspects of convergence rate and estimation accuracy of the costate.In addition, the RPM produces the most accurate propagated solution among three methods and the discretization of RPM spans the entire interval without overlap.So it is desirable to implement for multi-interval optimization problems.And in this paper, the RPM method is selected for air-to-air missile trajectory optimization.
Many researchers have studied pseudospectral methods for trajectory optimization in different background, such as Mars atmospheric entry trajectory optimization [11], onboard trajectory generation [12], and UAV path planning [13].Though some achievements have been made, there are still some intractable problems to solve.For example, to improve the accuracy of the solution, we can add the number of the collocation points.However, with increase of the points, the NLP problem becomes much more complex and the computation efficiency decreases.Moreover, the high efficiency and accuracy of pseudospectral method depend on the smoothness of the problem.If the control and state have discontinuous points, the convergence rate and accuracy of the solution will decrease.In allusion to above problems, some improved optimization algorithms were introduced [14,15].Zhao and Tsiotras [14] proposed a simple but effective method which employs a probability density function to solve the trajectory optimization problem by generating a fixedorder mesh.Gong et al. [15] designed a spectral algorithm for pseudospectral method which uses the pseudospectral differentiation matrix to locate switches, kinks, corners, and other discontinuities.In this method, the time interval was divided into smaller subintervals and the distribution of the nodes was controlled by knotting method.This paper advances pseudospectral methods in computational optimal control.
In allusion to the problems, such as too many collocation points in mesh optimization, high dimensionality of NLP problem, and the low computational efficiency, a multiinterval mesh refinement method with unfixed collocation points was introduced in this paper.The mesh needs to be further refined if the relative error is larger than the given value, and if nonsmooth point exists in the mesh, the number of mesh intervals will be increased; otherwise the number of collocation points will be increased.

Air-to-Air Missile Trajectory Optimization Problem
The range of air-to-air missile is beyond 100 km.According to the working time of the engine and handover of midcourse and terminal guidance, the missile trajectory can be divided into four phases, as shown in Figure 1.The first phase is from the missile separated from the carrier to the engine starting working, where the missile is in a state of low speed with no thrust.The second phase is the working process of engine, where the missile is speeding up, and the flight altitude, velocity, overload, and dynamic pressure of missile reach the peak value.The third section is from shutdown of the engine to the seeker capturing the target, where the velocity and altitude of missile change slowly.The last phase is the terminal guidance.It is necessary to adopt multi-interval trajectory optimization for the previous three phases because of the sharp change of thrust and the discontinuous control caused by the switch of the engine.

Motion Model.
When we design the nominal optimal trajectory, the air-to-air missile should be directly guided towards the target, so the dynamic model is considered in the longitudinal plane.The mass point motion equations are given as follows: where V, , , and  are velocity, flight path angle, and location in the inertial coordinate respectively;  is the missile mass;  is gravity acceleration;  is angle of attack;  is engine thrust;  is the dynamic pressure; and   ,   denote the drag and lift coefficients, respectively, which are expressed as a function of Ma, the Mach number, and , the angle of attack: where  0 ,  denote zero-lift drag coefficient and induced drag coefficient, respectively and    denotes the partial derivative of   with respect to .  is the dynamic pressure: And  is the atmosphere density expressed as where  0 = 1.2250 kg/m 3 ,  0 = 7254.3m.

Boundary Conditions.
As the missile is launched from the carrier, the following initial conditions need to be satisfied: The initial point of the trajectory optimization is the starting control point after launching from the carrier, and the end point is the initial point of terminal guidance.In order to successfully capture the target, the seeker must satisfy certain constraints which can be expressed as where V min denotes the minimum velocity that terminal guidance requires to hit the target. min ,  max denote the minimum and maximum of terminal flight path angle.(  ) denotes the distance between missile and target at the beginning of terminal guidance. max denotes detection range of the seeker.

Path Constraints
Dynamic pressure constraint:  min ≤  ≤  max , where  min ,  max denote allowable minimum and maximum dynamic pressure.In order to maintain the basic aerodynamic flight, the dynamic pressure cannot be too small or too large.
Overload constraint:  ≤  max , where  max denotes the maximum allowable overload.Due to the stability limitation of the missile body structure, it is necessary to limit overload during the flight.
Altitude constraint: ℎ ≤ ℎ max , where ℎ max denotes the maximum allowable altitude.To ensure the normal work of engine, the altitude cannot be too high.

Connection Point Constraints.
In order to ensure smooth transition between two phases, the state and control should be the same at each phase connection point: where  denotes the th flight phase and subscripts 0 and denote the initial and terminal point, respectively.

Control Constraints. |𝛼| ≤ 𝛼 max
, where  max denotes the maximum angle of attack.Due to the stability limitation of the body structure and attitude control system, the angle of attack should not be too large.

Objective Function.
The objective is to increase the range with predetermined thrust profile.The performance index can be expressed as min  = − (  ) . (8)

Multi-Interval Radau Pseudospectral Method
For concision and without loss of generality, let us begin with discussing the standard form of nonlinear optimal control problem.
()  is the LGR differential matrix of dimension   × (  + 1) of mesh interval   .After LGR discretization, the nonlinear optimal control problem is transformed into the following NLP problem.

Mesh Interval Division.
If  () max >  and   ≥  at mesh interval   , where  denotes maximum allowable error, the nonsmooth points exist in this mesh interval.It is necessary to divide this mesh interval into  submesh intervals.And  can be calculated as follows: The upper bound  max on the number of submesh intervals can be calculated as follows.

Collocation Points Increasing. If 𝑒 (𝑘)
max >  and   <  at mesh interval   , the mesh interval   is smooth.In order to maintain the relative error estimation under , the error must be reduced  ()  max / times as its current value.If the requirements are satisfied by increasing the number of collocation points, the collocation points of mesh interval   can be  ()  at th iteration and  (+1)  at ( + 1)th iteration.And  (+1)  can be expressed as follows: where  is relevant to  ()  .Because  (+1)  is an integer, (19) can be transformed into the following form: 4.3.Algorithm Flow.The calculation steps of multi-interval mesh refinement pseudospectral method are as follows and the flow chart is shown in Figure 2.
Step 1. Initialize the mesh and discretize it by the method proposed in this paper.Transform the optimal control problem into a NLP problem which can be solved by SNOPT then.
Step 2. If  () max <  at all mesh intervals, then terminate, or go to Step 3.
Step 3. If  ()  max <  at mesh interval   , go to Step 5, or go to Step 4.
Step 4. If   ≤  at grid interval   , increase the number of collocation points, or increase the number of the mesh intervals; then go to Step 5.
Step 5.According to the result of Step 2 to Step 4, build new mesh intervals and collocation points, and solve the current NLP problem; then return to Step 2.

Simulation and Analysis
The working time of engine is 15 s.Simulation initial constraints are as follows:  0 = 0 m, ℎ 0 = 10 km,  0 = 400 m/s, and  0 = 0 ∘ ; simulation process constraints: ℎ ≤ 30 km, || ≤ 10 g, || ≤ 20 ∘ , and 10 kPa ≤  ≤ 200 kPa; simulation terminal constraints: The number of collocation points at each mesh interval in multi-interval method which is expressed as ℎ − ( min ,  max ) is unfixed. min ,  max denote the minimum and maximum allowable number of collocation points.ℎ −   denotes the ℎ method, and the number of collocation points at each mesh interval is a fixed value,   .The simulation results of this paper are obtained from Matlab simulation on Lenovo CPU 3.4 GHz Core i7 Intel computer, and NLP problem is solved by SNOPT. denotes the iteration times.And the maximum allowable error of mesh optimization accuracy is  = 10 −6 ,  = 1.2.During the initialization, the trajectory is divided into 3 phases, each with 10 mesh intervals.There are 2 initial collocation points in each mesh interval.
The trajectory curve is shown in Figure 4.In the simulation, the solution converges after 4 iterations.The time consumption of simulation is 9.7 s.In order to analyze its convergence, the local part of Figure 4 is enlarged, as shown in Figure 5, from which we can see that the effect of initiation is good, but the first iteration error is large, then the mesh converges quickly, and the fourth iteration satisfies the accuracy requirement, which means the error is under .
We can see from Figure 3 that the altitude changes rapidly at the beginning and the end, which indicates that the states change rapidly at the beginning and end of the mesh.The distribution of the mesh collocation points is shown in Figure 5, from which we can see that collocation points at the beginning and the end are obviously much more than that at the middle section, which verifies the adaptation of the method.
When the ℎ − 4 method is adopted, the distribution of collocation points is shown in Figure 6.And when  = 0, 1, the collocation points in ℎ method are similar to those in multi-interval method.But when  = 2, 3, 4, the collocation points of ℎ method are more evenly distributed in each mesh interval, and the total collocation points are significantly more than those in multi-interval method, because the multiinterval method has fewer collocation points in the middle section with small change of state and control.
Then the multi-interval method and ℎ method are compared in terms of computational efficiency.The running time, the total number of collocation points , and the mesh iteration times  are shown in Table 1.
Table 1 shows that, in the aspect of time consumption, ℎ− (3,8) is the smallest of the multi-interval method and ℎ−4 is the smallest of ℎ method.The total number of collocation points of ℎ − (3, 8) method ( = 175) is less than ℎ − 4 method ( = 180).In fact, as shown in Table 1, for ℎ−  and ℎ−( min ,  max ), if   and  min are the same, the colocation points of multi-interval method are fewer than those of ℎ method, which indicates that the multi-interval method is more efficient than ℎ method.In addition, the performance index of these two methods is consistent.The simulation results of ℎ − (3, 12) method are shown in Figures 7-12.In Figure 7, the trajectory is smooth and the maximum altitude is 30 km.In order to maximize the range, the trajectory is optimized as high altitude trajectory.In practical application, the boost-glide trajectory has the characteristics of increasing range.In this paper, the maximized range is considered as the performance index and the trajectory has the boost-glide process, which is consistent with the practical situation and verifies the effectiveness of the method.From Figures 8-12, we can know that the change of angle of attack is smooth, and dynamic pressure and overload are satisfying the path constraints, which ensures the stability of flight.

Conclusion
In order to increase the range of existing air-to-air missile, the trajectory was divided into four phases according to the working time of the engine and handover of midcourse and terminal guidance, and then the optimization model was built.In allusion to the trajectory optimization problems with parametric discontinuity and poor fitting effect of state and control at the nonsmooth points, multi-interval mesh refinement method was introduced in this paper.Compared with ℎ method, this method can effectively decrease the dimensionality of NLP problems and improve the convergence rate with the same precision.Its effectiveness was verified by simulation.
],∈[1,...,  +1] When the mesh interval needs to be further refined, it is necessary to decide whether to add the number of collocation points or the number of mesh intervals.In this paper, if the mesh interval is nonsmooth, we increase the number of mesh intervals; otherwise we increase the number of collocation points.For simplicity, in the mesh interval   , () ()|.Similarly,  (−1)  is the local maximum of | Ẍ(−1)  ()|,   ∈   ( = 1, . . .,   ,  = 1, . . .,   ). denotes the current number of iterations.

Table 1 :
Comparison of two methods.