# Trajectory Optimization Based on Multi-Interval Mesh Refinement Method

Authors: Ningbo Li; Humin Lei; Lei Shao; Tao Liu; Bin Wang
Journal: Mathematical Problems in Engineering
 (2017)

---

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


---

## Body

## 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 Figure1. 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.Figure 1
Air-to-air missile trajectory.


### 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:(1)v˙=Pcos⁡α-CxqSm-gsin⁡θ,θ˙=Psin⁡α+CyqSmv-gcos⁡θv,x˙=vcos⁡θ,y˙=vsin⁡θ,where v, θ, x, and y are velocity, flight path angle, and location in the inertial coordinate respectively; m is the missile mass; g is gravity acceleration; α is angle of attack; P is engine thrust; q is the dynamic pressure; and Cx, Cy denote the drag and lift coefficients, respectively, which are expressed as a function of Ma, the Mach number, and α, the angle of attack:(2)Cy=CyαMaα,Cx=Cx0Ma+KMa·Cy2,where Cx0, K denote zero-lift drag coefficient and induced drag coefficient, respectively and Cyα denotes the partial derivative of Cy with respect to α. q is the dynamic pressure:(3)q=0.5ρv2.Andρ is the atmosphere density expressed as(4)ρ=ρ0exp⁡-yy0,where ρ0=1.2250 kg/m3, y0=7254.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:(5)Xt0=v0θ0x0h0.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(6)vtf≥vfmin,θfmin≤θtf≤θfmax,rtf≤rfmax,where vfmin denotes the minimum velocity that terminal guidance requires to hit the target. θfmin, θfmax denote the minimum and maximum of terminal flight path angle. r(tf) denotes the distance between missile and target at the beginning of terminal guidance. rfmax denotes detection range of the seeker.


#### 2.2.2. Path Constraints

Dynamic pressure constraint:qmin≤q≤qmax, where qmin, qmax 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:n≤nmax, where nmax 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:h≤hmax, where hmax 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:(7)t0i+1=tfi,U0i+1=Ufi,X0i+1=Xfi,where i denotes the ith flight phase and subscripts 0 and fdenote the initial and terminal point, respectively.


#### 2.2.4. Control Constraints

α
≤
α
m
a
x, where αmaxdenotes 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(8)minJ=-xtf.


## 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:(1)v˙=Pcos⁡α-CxqSm-gsin⁡θ,θ˙=Psin⁡α+CyqSmv-gcos⁡θv,x˙=vcos⁡θ,y˙=vsin⁡θ,where v, θ, x, and y are velocity, flight path angle, and location in the inertial coordinate respectively; m is the missile mass; g is gravity acceleration; α is angle of attack; P is engine thrust; q is the dynamic pressure; and Cx, Cy denote the drag and lift coefficients, respectively, which are expressed as a function of Ma, the Mach number, and α, the angle of attack:(2)Cy=CyαMaα,Cx=Cx0Ma+KMa·Cy2,where Cx0, K denote zero-lift drag coefficient and induced drag coefficient, respectively and Cyα denotes the partial derivative of Cy with respect to α. q is the dynamic pressure:(3)q=0.5ρv2.Andρ is the atmosphere density expressed as(4)ρ=ρ0exp⁡-yy0,where ρ0=1.2250 kg/m3, y0=7254.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:(5)Xt0=v0θ0x0h0.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(6)vtf≥vfmin,θfmin≤θtf≤θfmax,rtf≤rfmax,where vfmin denotes the minimum velocity that terminal guidance requires to hit the target. θfmin, θfmax denote the minimum and maximum of terminal flight path angle. r(tf) denotes the distance between missile and target at the beginning of terminal guidance. rfmax denotes detection range of the seeker.


### 2.2.2. Path Constraints

Dynamic pressure constraint:qmin≤q≤qmax, where qmin, qmax 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:n≤nmax, where nmax 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:h≤hmax, where hmax 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:(7)t0i+1=tfi,U0i+1=Ufi,X0i+1=Xfi,where i denotes the ith flight phase and subscripts 0 and fdenote the initial and terminal point, respectively.


### 2.2.4. Control Constraints

α
≤
α
m
a
x, where αmaxdenotes 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.2.1. Boundary Conditions

As the missile is launched from the carrier, the following initial conditions need to be satisfied:(5)Xt0=v0θ0x0h0.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(6)vtf≥vfmin,θfmin≤θtf≤θfmax,rtf≤rfmax,where vfmin denotes the minimum velocity that terminal guidance requires to hit the target. θfmin, θfmax denote the minimum and maximum of terminal flight path angle. r(tf) denotes the distance between missile and target at the beginning of terminal guidance. rfmax denotes detection range of the seeker.


## 2.2.2. Path Constraints

Dynamic pressure constraint:qmin≤q≤qmax, where qmin, qmax 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:n≤nmax, where nmax 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:h≤hmax, where hmax 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:(7)t0i+1=tfi,U0i+1=Ufi,X0i+1=Xfi,where i denotes the ith flight phase and subscripts 0 and fdenote the initial and terminal point, respectively.


## 2.2.4. Control Constraints

α
≤
α
m
a
x, where αmaxdenotes 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(8)minJ=-xtf.


## 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], τ∈[-1,+1] is composed of k mesh intervals Sk=[Tk-1,Tk],  k=1,…,K, where(9)-1=T0<T1<⋯<Tk=+1,⋃k=1KSk=-1,+1.X
(
k
)
(
τ
) and U(k)(τ) denote state and control in interval Sk. They can be discretized at Legendre-Gauss-Radau (LGR) points expressed as follows:(10)Xkτ=∑j=1Nk+1Xjkτljkτ,Ukτ=∑j=1NkUjkτljkτ,ljkτ=∏l=1,l≠jNk+1τ-τlkτjk-τlk,where lj(k)(τ),  j=1,…,Nk+1, denotes the basis of Lagrange polynomials. τ1(k),…,τNk(k) are LGR collocation points of mesh interval [Tk-1,Tk) and τNk+1(k)=Tk is a noncollocated node point.We can get state equations as follows:(11)dXkτdτ=∑j=1Nk+1Xjkdljkτdτ,∑j=1Nk+1DijkXjk=tf-t02fXikτ,Uikτ,tτik,t0,tf,where(12)Dijk=dljkτikdτ,i=1,…,Nk,j=1,…,Nk+1.D
i
j
(
k
) is the LGR differential matrix of dimension Nk×(Nk+1) of mesh interval Sk. After LGR discretization, the nonlinear optimal control problem is transformed into the following NLP problem.Objective function:(13)J=MX11,t0,XNK+1K,tf+tf-t02∑k=1K∑j=1NkωjkLXjk,Ujk,tτjk,t0,tf.Discretized differential equations:(14)∑j=1Nk+1DijkXjk-tf-t02fXik,Uik,tτik,t0,tf=0.Discretized path constraints:(15)cmin≤cXik,Uik,tτik,t0,tf≤cmax.Discretized boundary conditions:(16)bmin≤bX11,t0,XNk+1K,tf≤bmax.The state is continuous at points(T1,…,TK-1), so the condition XNk+1(k)=X1(k+1) must be satisfied.


## 4. Mesh Refinement Method

### 4.1. Relative Error Estimation

It is assumed that there areNk LGR collocation points at mesh interval Sk=[Tk-1,Tk]. When estimating whether to add the collocation points, we suppose that there are Mk=Nk+1 LGR collocation points τ^1(k),…,τ^Mk(k), where τ^1(k)=τ1(k)=Tk-1, τ^Mk+1(k)=Tk. The state is approximated as X(τ^1(k)),…,X(τ^Mk(k)) at LGR collocation points τ^1(k),…,τ^Mk(k). The basis of Lagrange interpolation of control is expressed as follows:(17)Ukτ=∑j=1NkUjkτl^jkτ,l^jkτ=∏l=1l≠jNkτ-τlkτjk-τlk.AndX^(k)(τ^j(k)) can be expressed as follows:(18)X^kτ^jk=Xkτk-1+tf-t02∑l=1MkI^jlkfXkτ^lk,Ukτ^lk,tτ^lk,t0,tf,j=2,…,Mk+1,where I^jl(k), j,l=1,…,Mk, are integral weight matrices of dimension Mk×Mk of LGR points τ^1(k),…,τ^Mk(k). The absolute error and relative error between X(τ^j(k)) and X^(τ^j(k)) are expressed, respectively, as follows:(19)Eikτ^lk=X^ikτ^lk-Xikτ^lk,(20)eikτ^lk=Eikτ^lk1+maxj∈1,…,Nk+1,k∈1,…,K⁡Xikτlk,i=1,…,nx,l=1,…,Mk+1.The maximum relative error estimation at mesh intervalSk is expressed as follows:(21)emaxk=maxi∈1,…,nx,j∈1,…,Mk+1eikτ^jk.


### 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 intervalSk, Pij(M) is the local maximum of X¨i(M)(τ). Similarly, Pij(M-1) is the local maximum of X¨i(M-1)(τ), τij∈Sk  (i=1,…,nx,j=1,…,Mk). M denotes the current number of iterations.If(22)Rij=PijMPijM-1≥R-,where R- denotes the given ratio, mesh interval Sk is nonsmooth [17].


#### 4.2.2. Mesh Interval Division

Ifemax(k)>ε and Rij≥R- at mesh interval Sk, where ε denotes maximum allowable error, the nonsmooth points exist in this mesh interval. It is necessary to divide this mesh interval into S submesh intervals. And S can be calculated as follows:(23)S=logNkekε.The upper boundHmax on the number of submesh intervals can be calculated as follows.Ife(k)≫ε (such as 106), Hmax can be 15–25 [18], while if e(k)→ε, Hmax is close to 0(24)S=min⁡Hk+1,Hmax.


#### 4.2.3. Collocation Points Increasing

Ifemax(k)>ε and Rij<R- at mesh interval Sk, the mesh interval Sk is smooth. In order to maintain the relative error estimation under ε, the error must be reduced emax(k)/ε times as its current value. If the requirements are satisfied by increasing the number of collocation points, the collocation points of mesh interval Sk can be Nk(M) at Mth iteration and Nk(M+1) at M+1th iteration. And Nk(M+1) can be expressed as follows:(25)NkM+1=NkMekMε1/q-2.5,where q is relevant to Nk(M). Because Nk(M+1) is an integer, (19) can be transformed into the following form:(26)NkM+1=NkMekMε1/q-2.5.


### 4.3. Algorithm Flow

The calculation steps of multi-interval mesh refinement pseudospectral method are as follows and the flow chart is shown in Figure2.Figure 2
Multi-interval mesh refinement method.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.
Ifemax(k)<ε at all mesh intervals, then terminate, or go to Step 3.Step 3.
Ifemax(k)<ε at mesh interval Sk, go to Step 5, or go to Step 4.Step 4.
IfRij≤R- at grid interval Sk, 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 Step2 to Step 4, build new mesh intervals and collocation points, and solve the current NLP problem; then return to Step 2.


## 4.1. Relative Error Estimation

It is assumed that there areNk LGR collocation points at mesh interval Sk=[Tk-1,Tk]. When estimating whether to add the collocation points, we suppose that there are Mk=Nk+1 LGR collocation points τ^1(k),…,τ^Mk(k), where τ^1(k)=τ1(k)=Tk-1, τ^Mk+1(k)=Tk. The state is approximated as X(τ^1(k)),…,X(τ^Mk(k)) at LGR collocation points τ^1(k),…,τ^Mk(k). The basis of Lagrange interpolation of control is expressed as follows:(17)Ukτ=∑j=1NkUjkτl^jkτ,l^jkτ=∏l=1l≠jNkτ-τlkτjk-τlk.AndX^(k)(τ^j(k)) can be expressed as follows:(18)X^kτ^jk=Xkτk-1+tf-t02∑l=1MkI^jlkfXkτ^lk,Ukτ^lk,tτ^lk,t0,tf,j=2,…,Mk+1,where I^jl(k), j,l=1,…,Mk, are integral weight matrices of dimension Mk×Mk of LGR points τ^1(k),…,τ^Mk(k). The absolute error and relative error between X(τ^j(k)) and X^(τ^j(k)) are expressed, respectively, as follows:(19)Eikτ^lk=X^ikτ^lk-Xikτ^lk,(20)eikτ^lk=Eikτ^lk1+maxj∈1,…,Nk+1,k∈1,…,K⁡Xikτlk,i=1,…,nx,l=1,…,Mk+1.The maximum relative error estimation at mesh intervalSk is expressed as follows:(21)emaxk=maxi∈1,…,nx,j∈1,…,Mk+1eikτ^jk.


## 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 intervalSk, Pij(M) is the local maximum of X¨i(M)(τ). Similarly, Pij(M-1) is the local maximum of X¨i(M-1)(τ), τij∈Sk  (i=1,…,nx,j=1,…,Mk). M denotes the current number of iterations.If(22)Rij=PijMPijM-1≥R-,where R- denotes the given ratio, mesh interval Sk is nonsmooth [17].


### 4.2.2. Mesh Interval Division

Ifemax(k)>ε and Rij≥R- at mesh interval Sk, where ε denotes maximum allowable error, the nonsmooth points exist in this mesh interval. It is necessary to divide this mesh interval into S submesh intervals. And S can be calculated as follows:(23)S=logNkekε.The upper boundHmax on the number of submesh intervals can be calculated as follows.Ife(k)≫ε (such as 106), Hmax can be 15–25 [18], while if e(k)→ε, Hmax is close to 0(24)S=min⁡Hk+1,Hmax.


### 4.2.3. Collocation Points Increasing

Ifemax(k)>ε and Rij<R- at mesh interval Sk, the mesh interval Sk is smooth. In order to maintain the relative error estimation under ε, the error must be reduced emax(k)/ε times as its current value. If the requirements are satisfied by increasing the number of collocation points, the collocation points of mesh interval Sk can be Nk(M) at Mth iteration and Nk(M+1) at M+1th iteration. And Nk(M+1) can be expressed as follows:(25)NkM+1=NkMekMε1/q-2.5,where q is relevant to Nk(M). Because Nk(M+1) is an integer, (19) can be transformed into the following form:(26)NkM+1=NkMekMε1/q-2.5.


## 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 intervalSk, Pij(M) is the local maximum of X¨i(M)(τ). Similarly, Pij(M-1) is the local maximum of X¨i(M-1)(τ), τij∈Sk  (i=1,…,nx,j=1,…,Mk). M denotes the current number of iterations.If(22)Rij=PijMPijM-1≥R-,where R- denotes the given ratio, mesh interval Sk is nonsmooth [17].


## 4.2.2. Mesh Interval Division

Ifemax(k)>ε and Rij≥R- at mesh interval Sk, where ε denotes maximum allowable error, the nonsmooth points exist in this mesh interval. It is necessary to divide this mesh interval into S submesh intervals. And S can be calculated as follows:(23)S=logNkekε.The upper boundHmax on the number of submesh intervals can be calculated as follows.Ife(k)≫ε (such as 106), Hmax can be 15–25 [18], while if e(k)→ε, Hmax is close to 0(24)S=min⁡Hk+1,Hmax.


## 4.2.3. Collocation Points Increasing

Ifemax(k)>ε and Rij<R- at mesh interval Sk, the mesh interval Sk is smooth. In order to maintain the relative error estimation under ε, the error must be reduced emax(k)/ε times as its current value. If the requirements are satisfied by increasing the number of collocation points, the collocation points of mesh interval Sk can be Nk(M) at Mth iteration and Nk(M+1) at M+1th iteration. And Nk(M+1) can be expressed as follows:(25)NkM+1=NkMekMε1/q-2.5,where q is relevant to Nk(M). Because Nk(M+1) is an integer, (19) can be transformed into the following form:(26)NkM+1=NkMekMε1/q-2.5.


## 4.3. Algorithm Flow

The calculation steps of multi-interval mesh refinement pseudospectral method are as follows and the flow chart is shown in Figure2.Figure 2
Multi-interval mesh refinement method.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.
Ifemax(k)<ε at all mesh intervals, then terminate, or go to Step 3.Step 3.
Ifemax(k)<ε at mesh interval Sk, go to Step 5, or go to Step 4.Step 4.
IfRij≤R- at grid interval Sk, 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 Step2 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:x0=0 m, h0=10 km, V0=400 m/s, and θ0=0°; simulation process constraints: h≤30 km, n≤10 g, α≤20°, and 10kPa≤q≤200kPa; simulation terminal constraints: Vf≥400 m/s, θf=-45°.The number of collocation points at each mesh interval in multi-interval method which is expressed asph-(Nmin,Nmax) is unfixed. Nmin, Nmax denote the minimum and maximum allowable number of collocation points. h-Nf denotes the h method, and the number of collocation points at each mesh interval is a fixed value, Nf. 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. M denotes the iteration times. And the maximum allowable error of mesh optimization accuracy is ε=10-6, R-=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 Figure4. 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 Figure3 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.Figure 3
Trajectory curves of iterative process.Figure 4
Local trajectory curves of iterative process.Figure 5
Collocation points distribution ofph-(3,12).When theh-4 method is adopted, the distribution of collocation points is shown in Figure 6. And when M=0,1, the collocation points in h method are similar to those in multi-interval method. But when M=2,3,4, the collocation points of h 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.Figure 6
Collocation points distribution ofh-4.Then the multi-interval method andh method are compared in terms of computational efficiency. The running time, the total number of collocation points N, and the mesh iteration times M are shown in Table 1.Table 1
Comparison of two methods.
N
m
i
n
N
m
a
x
Times/s
N
M
h
2
2
34.5
346
6
h
3
3
14.5
200
5
h
4
4
13.1
180
5
h
5
5
14.5
190
5
ph
3
4
15.0
183
5
ph
3
5
16.7
177
5
ph
3
6
13.9
175
5
ph
3
8
9.69
175
4
ph
3
12
14.0
177
5
ph
4
5
12.4
158
5
ph
4
8
12.6
158
5
ph
4
12
12.7
158
5
ph
5
6
12.0
181
5
ph
5
8
11.7
181
5
ph
5
12
12.5
181
5Table1 shows that, in the aspect of time consumption, ph-(3,8) is the smallest of the multi-interval method and h-4 is the smallest of h method. The total number of collocation points of ph-(3,8) method (N=175) is less than h-4 method (N=180). In fact, as shown in Table 1, for h-Nf and ph-(Nmin,Nmax), if Nf and Nmin are the same, the colocation points of multi-interval method are fewer than those of h method, which indicates that the multi-interval method is more efficient than h method. In addition, the performance index of these two methods is consistent.The simulation results ofph-(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.Figure 7
Curve of trajectory.Figure 8
Curve of angle of attack.Figure 9
Curve of flight path angle.Figure 10
Curve of velocity.Figure 11
Curve of dynamic pressure.Figure 12
Curve of overload.


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



---

*Source: 8521368-2017-08-23.xml*