#### Abstract

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.

#### 1. Introduction

Beyond visual range (BVR) combat is being the main form of future air combat [1–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–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 fixed-order 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 multi-interval 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.

#### 2. 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.

##### 2.1. 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 , , , 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 , 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 aswhere kg/m^{3}, m.

##### 2.2. Constraint Conditions

###### 2.2.1. 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 aswhere denotes the minimum velocity that terminal guidance requires to hit the target. , denote the minimum and maximum of terminal flight path angle. denotes the distance between missile and target at the beginning of terminal guidance. denotes detection range of the seeker.

###### 2.2.2. Path Constraints

Dynamic pressure constraint: , where , 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: , where 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: , where denotes the maximum allowable altitude. To ensure the normal work of engine, the altitude cannot be too high.

###### 2.2.3. 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 flight phase and subscripts and denote the initial and terminal point, respectively.

###### 2.2.4. Control Constraints

, where 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.

##### 2.3. Objective Function

The objective is to increase the range with predetermined thrust profile. The performance index can be expressed as

#### 3. 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.

In the mesh refinement Radau pseudospectral method [16], is composed of mesh intervals , , where

and denote state and control in interval . They can be discretized at Legendre-Gauss-Radau (LGR) points expressed as follows:where , , denotes the basis of Lagrange polynomials. are LGR collocation points of mesh interval and is a noncollocated node point.

We can get state equations as follows:where

is the LGR differential matrix of dimension of mesh interval . After LGR discretization, the nonlinear optimal control problem is transformed into the following NLP problem.

Objective function:

Discretized differential equations:

Discretized path constraints:

Discretized boundary conditions:

The state is continuous at points , so the condition must be satisfied.

#### 4. Mesh Refinement Method

##### 4.1. Relative Error Estimation

It is assumed that there are LGR collocation points at mesh interval . When estimating whether to add the collocation points, we suppose that there are LGR collocation points , where , . The state is approximated as at LGR collocation points . The basis of Lagrange interpolation of control is expressed as follows:

And can be expressed as follows:where , , are integral weight matrices of dimension of LGR points . The absolute error and relative error between and are expressed, respectively, as follows:

The maximum relative error estimation at mesh interval is expressed as follows:

##### 4.2. Hp-Adaptive Mesh Update

###### 4.2.1. Nonsmooth Point Location

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 , is the local maximum of . Similarly, is the local maximum of , . denotes the current number of iterations.

Ifwhere denotes the given ratio, mesh interval is nonsmooth [17].

###### 4.2.2. Mesh Interval Division

If 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 on the number of submesh intervals can be calculated as follows.

If (such as ), can be 15–25 [18], while if , is close to 0

###### 4.2.3. Collocation Points Increasing

If and at mesh interval , the mesh interval is smooth. In order to maintain the relative error estimation under , the error must be reduced 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 iteration and at iteration. And can be expressed as follows:where is relevant to . Because 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 at all mesh intervals, then terminate, or go to Step 3.

*Step 3. *If 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.

#### 5. Simulation and Analysis

The working time of engine is 15 s. Simulation initial constraints are as follows: m, km, m/s, and ; simulation process constraints: km, g, , and ; simulation terminal constraints: m/s, °.

The number of collocation points at each mesh interval in multi-interval method which is expressed as is unfixed. , 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 , . 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 method is adopted, the distribution of collocation points is shown in Figure 6. And when , the collocation points in method are similar to those in multi-interval method. But when , 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 multi-interval 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, is the smallest of the multi-interval method and is the smallest of method. The total number of collocation points of method () is less than method (). In fact, as shown in Table 1, for and , if and 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 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.

#### 6. 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.

#### Conflicts of Interest

The authors declare no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

#### Acknowledgments

This study was cosupported by the National Natural Science Foundation of China (Grant no. 61573374 and no. 61503408) and Aeronautical Science Foundation of China (Grant no. 20140196004 and no. 20150196006).