# Variable-Time-Domain Online Neighboring Optimal Trajectory Modification for Hypersonic Interceptors

Authors: Ningbo Li; Humin Lei; Jin Zhou; Lei Shao; Bin Wang
Journal: International Journal of Aerospace Engineering
 (2017)

---

## Abstract

The predicted impact point (PIP) of hypersonic interception changes continually; therefore the midcourse guidance law must have the ability of online trajectory optimization. In this paper, an online trajectory generation algorithm is designed based on neighboring optimal control (NOC) theory and improved indirect Radau pseudospectral method (IRPM). A trajectory optimization model is designed according to the features of operations in near space. Two-point boundary value problems (TPBVPs) are obtained based on NOC theory. The second-order linear form of transversality conditions is deduced backward to express the modifications of terminal states, costates, and flight time in terms of current state errors and terminal constraints modifications. By treating the current states and the optimal costates modifications as initial constraints and perturbations, the feedback control variables are obtained based on improved IRPM and nominal trajectory information. The simulation results show that when the changes of terminal constraints are not relatively large, this method can generate a modified trajectory effectively with high precision of terminal modifications. The design concept can provide a reference for the design of the online trajectory generation system of hypersonic vehicles.


---

## Body

## 1. Introduction

After decades of development, hypersonic technology has made a lot of progress [1–3]. Hypersonic vehicles will be used in military as precision strike weapons or platforms in the near future. The demand on the research of advanced intercept and defense technology is urgent.In allusion to the high speed and large maneuverability of hypersonic weapons, the interceptor should adopt the compound guidance strategy to improve the success rate of interception. From the end of program control flight to the target acquisition of terminal guidance, the interceptor spends most of the time flying in the midcourse guidance. So the flight performance of midcourse guidance determines the overall performance.If we adopt the traditional guidance method to intercept the high dynamic target directly, it will cause the frequent changes of the trajectory, which is not conducive to the implementation of the interception and leads to unnecessary loss of energy. So it is a reasonable strategy to design the trajectory of midcourse guidance. The main objective of the midcourse guidance is to minimize the energy consumption in the flight process of interceptor and enter the terminal guidance with the best relative geometric relationship [4]. At present, the challenges which the midcourse guidance law design of near space interceptors faces are mainly in the following two aspects: on the one hand, during the flight process of interceptors, problems such as structural strength, thermal protection, normal work of the engine, and control stability give rise to strict requirements on the heat flux, dynamic pressure, overload, control, and other characteristics of the process. Additionally, the predicted impact point (PIP) provides terminal position constraints, and the terminal guidance acquisition state-space becomes the strong constraints conditions of interceptors [5]. On the other hand, in actual flight process, the variation ranges of altitude and velocity are relatively large. Due to the initial conditions errors, atmospheric environment changes, and great uncertainty of aerodynamic model and navigation equipment, it will cause trajectory tracking errors. Furthermore, the nominal trajectory generated before launch is based on remote detection results. In this condition, tracking and prediction of the target have great error. With the approaching of the interceptor and the target, tracking and prediction accuracy will be improved by means of approaching detection with on-board equipment. And the target is also likely to maneuver actively. All those will lead to the change of PIP, and terminal constraints also need updates accordingly. Therefore, the designed midcourse guidance law must have the ability of online trajectory optimization.Taking all kinds of terminal constraints and a series of process characteristics into account, the online optimization problem is actually a complex nonlinear optimization problem [6], and many scholars have studied it. Lin and Tsai [7, 8] obtained a trajectory shaping guidance law by simplifying the model and realized online trajectory generation. But the trajectory was not true optimal due to too much approximation. Indig et al. [6, 9, 10] first linearized the dynamic model of interceptors. According to the location of PIP, the optimal guidance law was designed based on the trajectory shaping guidance. Then, in allusion to the actual nonlinear dynamic model, the terminal angle constraints of interceptor were added according to the midcourse and terminal guidance handover. Finally, the optimal trajectory model was deduced based on Pontryagin minimum principle and solved by Gauss pseudospectral method. Dwivedi et al. [11, 12] studied the problem of midcourse guidance trajectory planning using model predictive static planning method. Based on the state equations of discretized system, the control variables were updated by solving the costate vectors of the whole time period, whose coefficient matrix could be obtained by recursive method. The simulation results show that this method has high accuracy and efficiency. Yakimenko et al. [13–15] proposed a trajectory shaping guidance law based on the inverse dynamic of virtual domain. In this method, the coordinates of aircraft were expressed by higher order polynomials with virtual arc, and the dynamic equations were transformed from the time domain to the virtual domain. Introducing the virtual arc as an argument made it possible to optimize the speed history along the trajectory independently. Additionally, the number of the optimized variables was reduced and the integration process was avoided. By this method, the trajectory can be generated with good convergence robustness. In addition, some scholars have also studied this problem by pseudospectral method [16, 17], intelligent algorithm [18, 19], and so on.The aforementioned scholars have made many accomplishments in the domain of trajectory optimization; however, they mostly focused on the trajectory regeneration in regard to the research of online trajectory optimization problems. They completely abandoned the original trajectory and reoptimized again, which caused large amount of calculation and more time consumption, and put forward higher requirements on the storage capacity and computing ability of on-board computer. So it is not easy to realize in engineering. In allusion to the characteristics of near space targets, such as small overload, long flight time, and large flight range, prediction of impact point will be continually updated and the terminal constraints change correspondingly. When the perturbations of current states and changes of terminal constraints are not very large, according to the neighboring optimal control theory, the modified trajectory can be generated quickly within the neighborhood of the nominal trajectory using nominal trajectory data to realize online trajectory generation.We still need to solve the TPBVPs [20] when adopting above modification algorithm. The time-intensive backward integration has brought difficulties to the online implementation. Recently, Fahroo and Ross [21] and Yan et al. [22] proposed the indirect Legendre pseudospectral method (ILPM) for solving the TPBVPs by matrix inversion which improved solving efficiency a lot. In this method, the modifications of states and costates were solved by initial perturbations of states and the nominal trajectory information, and then the feedback control modifications were obtained to form the closed loop control. In this paper, the idea of ILPM is used for reference and the IRPM is proposed which is improved in two aspects. On the one hand, Fahroo and Hui only solved the trajectory tracking problems with initial perturbations, while the improved IRPM considers not only trajectory tracking errors but also modifications of terminal constraints. If computation speed is fast enough, the trajectory can be generated in real time to realize closed loop guidance without designing the feedback control law. On the other hand, Radau pseudospectral method is superior to Legendre pseudospectral method in approximating precision of control, state, and costate variables. In terms of computational efficiency, the two methods have little difference when solving the same scale problem.The main purpose of this paper is to design an online optimal modification algorithm and realize online trajectory optimization under the condition of changing terminal constraints. The rest of paper is organized as follows. In Section2, the nominal trajectory satisfying terminal and path constraints is generated offline. Next, the trajectory modified model with unspecified terminal time is deduced by NOC in Section 3 and is solved by improved IRPM in Section 4. In Section 5, numerical simulations are presented. Finally, some conclusions are drawn in Section 6.


## 2. Nominal Optimal Trajectory

### 2.1. Dynamic Model of Interceptor

When the nominal optimal trajectory is generated offline, the interceptor should be directly guided towards the PIP, so the dynamic model is considered in the longitudinal plane. The states of interceptor are represented by four-dimensional vector(v,θ,x,y) consisting of speed, flight path angle, and location in the inertial coordinate and given as follows [23]:(1)v˙=Pcos⁡α-CxqSm-gsin⁡θ,θ˙=Psin⁡α+CyqSmv-gcos⁡θv,x˙=vcos⁡θ,y˙=vsin⁡θ,where m is interceptor mass; q is the dynamic pressure; S is the reference area; g is gravitational acceleration; P is the engine thrust. In order to achieve a higher flight speed, the two-stage booster program is designed. The thrust calculation formulas are shown in (2)P=Vedmdt,Ve=gIsp,where Ve is the fuel gas velocity and Isp is the specific impulse, whose value is generally 200~300 s for the solid rocket engine [24]. This paper designs a two-stage solid rocket booster. The related parameters of the engine are shown in Table 1.Table 1
The related parameters of booster engine.
Parameter name
Value
I
s
p (stage 1/2)
280/200 (s)
V
e (stage 1/2)
30/8.9 (kg/s)
Burning time (stage 1/2)
20/50 (s)
Mass of booster shell
250 kg
Mass of interceptor
60 kgC
x
,
C
y denote the drag and lift coefficients, respectively, which are expressed as a function of Ma, the Mach number, and α, the angle of attack:(3)Cy=CyαMaα,Cx=Cx0Ma+KMa·Cy2,where Cx0,K denote zero-lift drag coefficient and induced drag coefficient, respectively; Cyα denotes the partial derivative of Cy with respect to α. q denotes the dynamic pressure:(4)q=0.5ρv2.Andρ is the atmosphere density expressed as(5)ρ=ρ0exp⁡-yy0,where ρ0=1.2250kg/m3,y0=7254.3m.


### 2.2. Trajectory Optimization Problem

In order to ensure kinematic kill effect of interception, terminal speed is usually selected as cost function:(6)J=ϕvtf,tf=-vf.In this paper, subscript0 denotes initial conditions and subscript f denotes terminal conditions. As is discussed above, the main objective of the midcourse guidance is to deliver the interceptor to a certain position (usually the PIP) with specified conditions to ensure a successful handover between midcourse and terminal guidance. So the position constraints xf,yf are given by PIP and the flight path angle constraint θf at terminal time is given by the terminal guidance acquisition conditions(7)ψ=θ-θfx-xfy-yfT=0,where 0 denotes zero matrix with the corresponding dimension.The above trajectory optimization problem is solved with the help of the software package GPOPS [25]. The whole trajectory was composed of three phases: boost stages 1 and 2 and nonpropulsive phase. Compared with other direct optimization methods, pseudospectral method can obtain higher accuracy with fewer nodes, which improves efficiency of the algorithm. The nominal optimal trajectory is shown in Figure 1.Figure 1
Nominal trajectory.


## 2.1. Dynamic Model of Interceptor

When the nominal optimal trajectory is generated offline, the interceptor should be directly guided towards the PIP, so the dynamic model is considered in the longitudinal plane. The states of interceptor are represented by four-dimensional vector(v,θ,x,y) consisting of speed, flight path angle, and location in the inertial coordinate and given as follows [23]:(1)v˙=Pcos⁡α-CxqSm-gsin⁡θ,θ˙=Psin⁡α+CyqSmv-gcos⁡θv,x˙=vcos⁡θ,y˙=vsin⁡θ,where m is interceptor mass; q is the dynamic pressure; S is the reference area; g is gravitational acceleration; P is the engine thrust. In order to achieve a higher flight speed, the two-stage booster program is designed. The thrust calculation formulas are shown in (2)P=Vedmdt,Ve=gIsp,where Ve is the fuel gas velocity and Isp is the specific impulse, whose value is generally 200~300 s for the solid rocket engine [24]. This paper designs a two-stage solid rocket booster. The related parameters of the engine are shown in Table 1.Table 1
The related parameters of booster engine.
Parameter name
Value
I
s
p (stage 1/2)
280/200 (s)
V
e (stage 1/2)
30/8.9 (kg/s)
Burning time (stage 1/2)
20/50 (s)
Mass of booster shell
250 kg
Mass of interceptor
60 kgC
x
,
C
y denote the drag and lift coefficients, respectively, which are expressed as a function of Ma, the Mach number, and α, the angle of attack:(3)Cy=CyαMaα,Cx=Cx0Ma+KMa·Cy2,where Cx0,K denote zero-lift drag coefficient and induced drag coefficient, respectively; Cyα denotes the partial derivative of Cy with respect to α. q denotes the dynamic pressure:(4)q=0.5ρv2.Andρ is the atmosphere density expressed as(5)ρ=ρ0exp⁡-yy0,where ρ0=1.2250kg/m3,y0=7254.3m.


## 2.2. Trajectory Optimization Problem

In order to ensure kinematic kill effect of interception, terminal speed is usually selected as cost function:(6)J=ϕvtf,tf=-vf.In this paper, subscript0 denotes initial conditions and subscript f denotes terminal conditions. As is discussed above, the main objective of the midcourse guidance is to deliver the interceptor to a certain position (usually the PIP) with specified conditions to ensure a successful handover between midcourse and terminal guidance. So the position constraints xf,yf are given by PIP and the flight path angle constraint θf at terminal time is given by the terminal guidance acquisition conditions(7)ψ=θ-θfx-xfy-yfT=0,where 0 denotes zero matrix with the corresponding dimension.The above trajectory optimization problem is solved with the help of the software package GPOPS [25]. The whole trajectory was composed of three phases: boost stages 1 and 2 and nonpropulsive phase. Compared with other direct optimization methods, pseudospectral method can obtain higher accuracy with fewer nodes, which improves efficiency of the algorithm. The nominal optimal trajectory is shown in Figure 1.Figure 1
Nominal trajectory.


## 3. Trajectory Modification Based on Neighborhood Optimization

### 3.1. Two-Point Boundary Value Problem

In the process of actual flight, when tracking errors of current states and changes of terminal constraints are not very large, according to the neighboring optimal control theory, we can further differentiate the first-order necessary conditions of the optimal trajectory to the second order to obtain perturbation equations, and the modified trajectory satisfying the changed terminal conditions can be obtained using the nominal trajectory data. At the same time, the cost function can still maintain certain optimality.First of all, the first-order optimal necessary conditions are derived based on optimal control. By introducing the costate variableλ=λvλθλxλyT with the same dimension as the state variable x=vθxyT, Hamilton equation can be expressed as(8)H=λTf=λVv˙+λθθ˙+λxx˙+λyy˙,where f=v˙θ˙x˙y˙T denotes system equations.Then the augmented performance index with terminal constraintsJ- can be expressed as(9)J-=ϕ+νTψf+∫t0tfH-λTx˙dt,where ν=νθνxνyT. Canonical equation and coupled equation can be expressed as(10)x˙=∂H∂λ,(11)λ˙=-∂H∂x,(12)∂H∂u=0.Considering the change of terminal conditions, the terminal time of modified trajectory will change accordingly. We can derive the equations of boundary conditions specified at an unspecified terminal time and the terminal constraints of Hamiltonian function along the optimal trajectory based on NOC [26, 27](13)ψxtf=0,(14)λtf=-∂ϕ∂x+νT∂ψ∂xt=tf=-∂Φ∂xt=tf,(15)0=H+∂ϕ∂tf+νT∂ψ∂tft=tf=H+∂Φ∂tft=tf.Equations (10)~(12) can be further differentiated to second order as follows: (16)δx˙=∂2H∂λ∂xδx+∂2H∂λ∂uδu,(17)δλ˙=-∂2H∂x2δx-∂2H∂x∂λδλ-∂2H∂x∂uδu,(18)0=∂2H∂u∂xδx+∂2H∂u∂λδλ+∂2H∂u2δu.If∂2H/∂u2 is nonsingular in the whole flight process, we can get the expression of control modifications δu by (18):(19)δu=-∂2H∂u2-1∂2H∂u∂xδx+∂2H∂u∂λδλ.The dynamic equations of state modificationsδx and costate modifications δλ can be obtained by substituting (19) into (16)~(17):(20)δx˙=Atδx-Btδλ,(21)δλ˙=Ctδx-ATtδλ,where(22)At=∂2H∂λ∂x-∂2H∂λ∂u∂2H∂u2-1∂2H∂u∂x,Bt=∂H∂λ∂u∂2H∂u2-1∂2H∂u∂λ,Ct=∂2H∂x2-∂2H∂x∂u∂2H∂u2-1∂2H∂u∂x.Equations (20)~(21) satisfy initial conditions δxt0=δx0,δλt0=δλ0 and terminal conditions (13)~(15) form the TPBVPs of δx,δλ.


### 3.2. Terminal Linearization and Backward Recursion

Equations (6), (15), and (16) were differentiated to second order with unspecified terminal time. Since Φ=ϕ+νTψ, we can get(23)dψ=∂ψ∂xδx+∂ψ∂t+∂ψ∂xdxdtdtft=tf,(24)dλtf=δλtf+dλdtt=tfdtf=∂2Φ∂x2δx+∂ψ∂xTdν+∂2Φ∂x∂tdtft=tf,(25)0=∂H∂x+∂2Φ∂t∂xδx+∂ψ∂tTdν+∂H∂λδλ+∂H∂t+∂2Φ∂t2dtft=tf.Equations (23)~(26) consist of n+q+1 linear equations with 2n+q+1 unknown variables. They are state modifications δx of dimension n, costate modifications δλ of dimension n, variables dν of dimension q, and variable dtf of dimension one. Variables dν and state modifications of dimension n-q unspecified at terminal time are defined as free variables (26)dμ=dν1,dν2,…,dνq,δxq+1,δxq+2,…,δxnt=tf.Since the changes of terminal constraints are not very large, (23)~(25) can be linearized based on small perturbation hypothesis, so that we get(27)δxδλdtft=tf=Pdμdψf=P11P12P21P22P31P32dμdψf,where(28)P11=1V˙λ˙V000V˙λ˙Vθ˙2θ˙x˙θ˙y˙θ˙λ˙Vθ˙x˙x˙2x˙y˙x˙λ˙Vθ˙y˙x˙y˙y˙2y˙λ˙V,P12=1V˙λ˙V000V˙λ˙V+θ˙λ˙θθ˙λ˙xθ˙λ˙yx˙λ˙θV˙λ˙V+x˙λ˙xx˙λ˙yy˙λ˙θy˙λ˙xV˙λ˙V+y˙λ˙y,P21=1V˙λ˙Vθ˙λ˙Vx˙λ˙Vy˙λ˙Vλ˙V2-V˙λ˙V+θ˙λ˙θx˙λ˙θy˙λ˙θλ˙Vλ˙θθ˙λ˙x-V˙λ˙V+x˙λ˙xy˙λ˙xλ˙Vλ˙xθ˙λ˙yx˙λ˙y-V˙λ˙V+y˙λ˙yλ˙Vλ˙y,P22=1V˙λ˙Vλ˙Vλ˙θλ˙Vλ˙xλ˙Vλ˙yλ˙θ2λ˙θλ˙xλ˙θλ˙yλ˙θλ˙xλ˙x2λ˙xλ˙yλ˙θλ˙yλ˙xλ˙yλ˙y2,P31=1V˙λ˙V-θ˙-x˙-y˙-λ˙V,P32=1V˙λ˙V-λ˙θ-λ˙x-λ˙y.Equations (20)~(21) can be expressed as the following form:(29)δΘ˙=AtBtCt-AtδΘ=GδΘ,where δΘ=δxδλT and the variables of above differential equation are discretized as M+1 points by Δt. We can get by backward recursion(30)δΘM-1=δΘM+GMδΘM-Δt=I-GMΔtδΘM,δΘk=δΘk+1+Gk+1δΘk+1-Δt=I-Gk+1ΔtδΘk+1=I-Gk+1ΔtI-Gk+2ΔtδΘk+2.Passing to the first point in turn we can get(31)δΘ0==GδΘM,where G=I-G1ΔtI-G2Δt⋯I-GMΔt, and according to (27) we have(32)δxδλt=tf=P11P12P21P22dμdψf.Substituting (32) into (31)(33)δxt0δλt0=G1G2G3G4P11P12P21P22dμdψf=Π1Π2Π3Π4dμdψf,where(34)Π1=G1P11+G2P21,Π2=G1P12+G2P22,Π3=G3P11+G4P21,Π4=G3P12+G4P22,and if Π1 is reversible, dμ can be expressed as(35)dμ=Π1-1δx0-Π2dψf.Substituting (35) into (27) we can get(36)δxf=P11Π1-1δx0-Π2dψf+P12dψf,δλf=P21Π1-1δx0-Π2dψf+P22dψf,dtf=P31Π1-1δx0-Π2dψf+P32dψf.


### 3.3. The Optimality Theory

Theorem 1.
The modified trajectories obtained in the nominal trajectory neighborhood based on the neighborhood optimal control theory meet the changed terminal constraints, and the performance index has second-order optimality.Proof.
According to the division integral theorem, the performance index can be expressed as(37)J-=ϕ+νTψf-λTxt0tf+∫t0tfH+λ˙Txdt.
Considering the first-order variation of the performance index and taking canonical equations (10)-(11) and coupled equation (12) into it, we can get by ignoring the high order terms:(38)δJ-=λTδx0+νTδψf+12δxT∂2Φ∂x2δxf+∫t0tf12δxTδuT∂2H∂x2∂2H∂x∂u∂2H∂u∂x∂2H∂u2δxδudt.
Multiplying (16) and (23) by the δλ and δν, respectively, we can get by taking the resulting equations into (38):(39)δJ-=λTδx0+νTδψf+12δxT∂2Φ∂x2δxf+dνT∂ψ∂xδxf+∫t0tfδλTδx˙-∂2H∂λ∂xδx+∂2H∂λ∂uδudt+∫t0tf12δxTδuT∂2H∂x2∂2H∂x∂u∂2H∂u∂x∂2H∂u2δxδudt.
Considering the second-order variation of (39), we can get according to the division integral theorem:(40)δ2J-=δxT∂2Φ∂x2+dνT∂ψ∂x-δλTδ2xf+∫t0tfδxT∂2H∂x∂u+δλT∂2H∂λ∂u+δuT∂2H∂u2δ2udt+∫t0tfddtδλT+δxT∂2H∂x2+δλT∂2H∂λ∂x+δuT∂2H∂u∂xδ2xdt.
Taking (17)-(18) and (24) into (40), we can get δ2J-=0 for any δ2x and δ2u. Therefore, the second-order optimality of the performance index can be guaranteed with satisfying the changed terminal constraints.


## 3.1. Two-Point Boundary Value Problem

In the process of actual flight, when tracking errors of current states and changes of terminal constraints are not very large, according to the neighboring optimal control theory, we can further differentiate the first-order necessary conditions of the optimal trajectory to the second order to obtain perturbation equations, and the modified trajectory satisfying the changed terminal conditions can be obtained using the nominal trajectory data. At the same time, the cost function can still maintain certain optimality.First of all, the first-order optimal necessary conditions are derived based on optimal control. By introducing the costate variableλ=λvλθλxλyT with the same dimension as the state variable x=vθxyT, Hamilton equation can be expressed as(8)H=λTf=λVv˙+λθθ˙+λxx˙+λyy˙,where f=v˙θ˙x˙y˙T denotes system equations.Then the augmented performance index with terminal constraintsJ- can be expressed as(9)J-=ϕ+νTψf+∫t0tfH-λTx˙dt,where ν=νθνxνyT. Canonical equation and coupled equation can be expressed as(10)x˙=∂H∂λ,(11)λ˙=-∂H∂x,(12)∂H∂u=0.Considering the change of terminal conditions, the terminal time of modified trajectory will change accordingly. We can derive the equations of boundary conditions specified at an unspecified terminal time and the terminal constraints of Hamiltonian function along the optimal trajectory based on NOC [26, 27](13)ψxtf=0,(14)λtf=-∂ϕ∂x+νT∂ψ∂xt=tf=-∂Φ∂xt=tf,(15)0=H+∂ϕ∂tf+νT∂ψ∂tft=tf=H+∂Φ∂tft=tf.Equations (10)~(12) can be further differentiated to second order as follows: (16)δx˙=∂2H∂λ∂xδx+∂2H∂λ∂uδu,(17)δλ˙=-∂2H∂x2δx-∂2H∂x∂λδλ-∂2H∂x∂uδu,(18)0=∂2H∂u∂xδx+∂2H∂u∂λδλ+∂2H∂u2δu.If∂2H/∂u2 is nonsingular in the whole flight process, we can get the expression of control modifications δu by (18):(19)δu=-∂2H∂u2-1∂2H∂u∂xδx+∂2H∂u∂λδλ.The dynamic equations of state modificationsδx and costate modifications δλ can be obtained by substituting (19) into (16)~(17):(20)δx˙=Atδx-Btδλ,(21)δλ˙=Ctδx-ATtδλ,where(22)At=∂2H∂λ∂x-∂2H∂λ∂u∂2H∂u2-1∂2H∂u∂x,Bt=∂H∂λ∂u∂2H∂u2-1∂2H∂u∂λ,Ct=∂2H∂x2-∂2H∂x∂u∂2H∂u2-1∂2H∂u∂x.Equations (20)~(21) satisfy initial conditions δxt0=δx0,δλt0=δλ0 and terminal conditions (13)~(15) form the TPBVPs of δx,δλ.


## 3.2. Terminal Linearization and Backward Recursion

Equations (6), (15), and (16) were differentiated to second order with unspecified terminal time. Since Φ=ϕ+νTψ, we can get(23)dψ=∂ψ∂xδx+∂ψ∂t+∂ψ∂xdxdtdtft=tf,(24)dλtf=δλtf+dλdtt=tfdtf=∂2Φ∂x2δx+∂ψ∂xTdν+∂2Φ∂x∂tdtft=tf,(25)0=∂H∂x+∂2Φ∂t∂xδx+∂ψ∂tTdν+∂H∂λδλ+∂H∂t+∂2Φ∂t2dtft=tf.Equations (23)~(26) consist of n+q+1 linear equations with 2n+q+1 unknown variables. They are state modifications δx of dimension n, costate modifications δλ of dimension n, variables dν of dimension q, and variable dtf of dimension one. Variables dν and state modifications of dimension n-q unspecified at terminal time are defined as free variables (26)dμ=dν1,dν2,…,dνq,δxq+1,δxq+2,…,δxnt=tf.Since the changes of terminal constraints are not very large, (23)~(25) can be linearized based on small perturbation hypothesis, so that we get(27)δxδλdtft=tf=Pdμdψf=P11P12P21P22P31P32dμdψf,where(28)P11=1V˙λ˙V000V˙λ˙Vθ˙2θ˙x˙θ˙y˙θ˙λ˙Vθ˙x˙x˙2x˙y˙x˙λ˙Vθ˙y˙x˙y˙y˙2y˙λ˙V,P12=1V˙λ˙V000V˙λ˙V+θ˙λ˙θθ˙λ˙xθ˙λ˙yx˙λ˙θV˙λ˙V+x˙λ˙xx˙λ˙yy˙λ˙θy˙λ˙xV˙λ˙V+y˙λ˙y,P21=1V˙λ˙Vθ˙λ˙Vx˙λ˙Vy˙λ˙Vλ˙V2-V˙λ˙V+θ˙λ˙θx˙λ˙θy˙λ˙θλ˙Vλ˙θθ˙λ˙x-V˙λ˙V+x˙λ˙xy˙λ˙xλ˙Vλ˙xθ˙λ˙yx˙λ˙y-V˙λ˙V+y˙λ˙yλ˙Vλ˙y,P22=1V˙λ˙Vλ˙Vλ˙θλ˙Vλ˙xλ˙Vλ˙yλ˙θ2λ˙θλ˙xλ˙θλ˙yλ˙θλ˙xλ˙x2λ˙xλ˙yλ˙θλ˙yλ˙xλ˙yλ˙y2,P31=1V˙λ˙V-θ˙-x˙-y˙-λ˙V,P32=1V˙λ˙V-λ˙θ-λ˙x-λ˙y.Equations (20)~(21) can be expressed as the following form:(29)δΘ˙=AtBtCt-AtδΘ=GδΘ,where δΘ=δxδλT and the variables of above differential equation are discretized as M+1 points by Δt. We can get by backward recursion(30)δΘM-1=δΘM+GMδΘM-Δt=I-GMΔtδΘM,δΘk=δΘk+1+Gk+1δΘk+1-Δt=I-Gk+1ΔtδΘk+1=I-Gk+1ΔtI-Gk+2ΔtδΘk+2.Passing to the first point in turn we can get(31)δΘ0==GδΘM,where G=I-G1ΔtI-G2Δt⋯I-GMΔt, and according to (27) we have(32)δxδλt=tf=P11P12P21P22dμdψf.Substituting (32) into (31)(33)δxt0δλt0=G1G2G3G4P11P12P21P22dμdψf=Π1Π2Π3Π4dμdψf,where(34)Π1=G1P11+G2P21,Π2=G1P12+G2P22,Π3=G3P11+G4P21,Π4=G3P12+G4P22,and if Π1 is reversible, dμ can be expressed as(35)dμ=Π1-1δx0-Π2dψf.Substituting (35) into (27) we can get(36)δxf=P11Π1-1δx0-Π2dψf+P12dψf,δλf=P21Π1-1δx0-Π2dψf+P22dψf,dtf=P31Π1-1δx0-Π2dψf+P32dψf.


## 3.3. The Optimality Theory

Theorem 1.
The modified trajectories obtained in the nominal trajectory neighborhood based on the neighborhood optimal control theory meet the changed terminal constraints, and the performance index has second-order optimality.Proof.
According to the division integral theorem, the performance index can be expressed as(37)J-=ϕ+νTψf-λTxt0tf+∫t0tfH+λ˙Txdt.
Considering the first-order variation of the performance index and taking canonical equations (10)-(11) and coupled equation (12) into it, we can get by ignoring the high order terms:(38)δJ-=λTδx0+νTδψf+12δxT∂2Φ∂x2δxf+∫t0tf12δxTδuT∂2H∂x2∂2H∂x∂u∂2H∂u∂x∂2H∂u2δxδudt.
Multiplying (16) and (23) by the δλ and δν, respectively, we can get by taking the resulting equations into (38):(39)δJ-=λTδx0+νTδψf+12δxT∂2Φ∂x2δxf+dνT∂ψ∂xδxf+∫t0tfδλTδx˙-∂2H∂λ∂xδx+∂2H∂λ∂uδudt+∫t0tf12δxTδuT∂2H∂x2∂2H∂x∂u∂2H∂u∂x∂2H∂u2δxδudt.
Considering the second-order variation of (39), we can get according to the division integral theorem:(40)δ2J-=δxT∂2Φ∂x2+dνT∂ψ∂x-δλTδ2xf+∫t0tfδxT∂2H∂x∂u+δλT∂2H∂λ∂u+δuT∂2H∂u2δ2udt+∫t0tfddtδλT+δxT∂2H∂x2+δλT∂2H∂λ∂x+δuT∂2H∂u∂xδ2xdt.
Taking (17)-(18) and (24) into (40), we can get δ2J-=0 for any δ2x and δ2u. Therefore, the second-order optimality of the performance index can be guaranteed with satisfying the changed terminal constraints.


## 4. Solving Algorithm Based on Improved IRPM

In the last section, considering the changes of terminal constraints of some state variables, the optimal terminal modification of costates variablesδλf and optimal modification of terminal time dtf are derived at an unspecified terminal time based on the NOC. Improved IRPM will be used to solve the modified optimal trajectory according to δλf and dtf in this section.Sincedtf has been calculated in the last section, the nominal trajectory information which will be used to solve the modified trajectory changes correspondingly. The new nominal trajectory information can be obtained by interpolating the original nominal trajectory. Because the trajectory tracking errors have been taken into account when we solved δλf in the last section, we can regard current states of the interceptor as initial constraints and δλf as the perturbations when solving the modified trajectory. Then we can use the new nominal trajectory information to design the solving algorithm to generate a trajectory from the current states to the modified terminal states, whose performance index can still maintain certain optimality.The boundary conditions in last section change to the following form:(41)ψ′=θ-θ0x-x0y-y0T=0,(42)λt0=-νT∂ψ′∂xt=t0.Equations (41)~(42) can be expressed as follows:(43)Ωxt0,λt0=0.Variation of (43) is solved as(44)∂Ω∂xδx+∂Ω∂λδλt=t0=0.Now we use improved IRPM to solve the optimal trajectory modified model. Firstly, the variables are discretized at Legendre-Gauss-Radau (LGR) points. Since the time domain for the above problem is[t0,tf], while Radau pseudospectral method is used in the time domain [-1,1], it is necessary to introduce a new time variable τ:(45)τ=2ttf-t0-tf+t0tf-t0.States and costates variables are approximated as follows:(46)δx≈δxNτ=∑l=0N-1δxτlLlτ,δλ≈δλNτ=∑l=0N-1δλτlLlτ,where τ0=-1 and τl(l=1,2,…,N) are N LGR points distributed at the interval (-1,1], which are defined as zero points of P˙Nτ-P˙N-1τ. P˙N is Legendre polynomial of order N. And Li is Lagrange interpolation polynomial basis function:(47)Li=∏j=0,j≠iKτ-τjτi-τj.The derivatives ofδxN(τ) and δλN(τ) at the point of τl are obtained by differentiating (46) as follows:(48)δx˙τk≈δx˙Nτk=∑l=0N-1δxτlL˙lτ=∑l=0N-1Dklδxτl,δλτk≈δλNτk=∑l=0N-1δλτlL˙lτ=∑l=0N-1Dklδλτl,where Dkl is the difference matrix of dimension (N+1)×N:(49)Dkl=q˙τkτk-τlq˙τl,l≠k,q¨τl2q˙τl,l=k,qτl=1+τlPNτl-PN-1τl.The TPBVPs ofδx,δλ are transformed into the following algebraic equations:(50)∑j=1NDijδxjτ-tf-t02Aiτδxiτ-Biτδλiτ=0,(51)∑j=1NDijδλjτ+tf-t02Ciτδxiτ+AiTτδλiτ=0,(52)∂Ωxδx+∂Ωλδλτ=τ0=0.0 represents the zero matrix with the corresponding dimension. For simplicity, substituting X=[δx0T,δx1T,…,δxN-1T]T,Λ=[δλ0T,δλ1T,…,δλN-1T]T into (50)~(51) we can get(53)EX-tf-t02FΛ=0,tf-t02GX+HΛ=0,where E,F,G,H are matrices of dimension [n(N+1)×nN]:(54)Eij=DijIni≠jDijIn-tf-t02Aii=j,Fij=0ni≠jBii=j,Gij=0ni≠jCii=j,Hij=DijIni≠jDijIn+tf-t02AiTi=j,where 0n and In are zero and identity matrices of dimension n×n. The purpose of this algorithm is to solve (53) subjected to the boundary conditions (52). These equations can be summarized in matrix form as(55)Γ1Γ2tf-t02GHE-tf-t02FXΛ=VZ=000,where ZT=XT,ΛT, Γ1 and Γ2 are matrices of dimension n×n(N+1). Γ1=0n,…,0n,∂Ω/∂x,Γ2=0n,…,0n,∂Ω/∂λ. Dividing V and Z as V=VeVf and Z=Zeδλf we get(56)VeZe+VfδλN=0,where Ve of dimension n(2N+1)×n(2N+1) and Vf of dimension n(2N+1)×n are block matrices of V. Ze=[δx0T,δx1T,…,δxNT,δλ0T,…,δλN-1T]T is a vector of dimension n(2N+1), which can be obtained by solving (57)Ze=-Ve∖VfδλN=WδλN.SinceZ=[δx0T,…,δxNT,δλ0T,…,δλNT]T, we can get(58)Z=ZeδλN=WInδλN=WxWλδλN,where Wx and Wλ of dimension n(N+1)×n are the block matrices of WInT. Then we have(59)δxk=WxkδλN,δλk=WλkδλN,where Wxk of dimension n×n is block matrix of Wx and Wλk of dimension n×n is block matrix of Wλ. Subscript k represents the kth LGR point.Substituting (59) into (19), we can get the optimal control modifications:(60)δuτi=-Huu-1HuxWxi+fuTWλiδλN.SinceδλN=δλf has been solved in the last section, the modifications of states, costates, and controls at the LGR points can be obtained by (59)-(60) and the values at instants of time between the nodes can be obtained by interpolation. The above solutions are obtained without a wide range of iterative optimization process, so it can ensure the high computational efficiency.The modified algorithm can be concluded as Figure2 and specific procedures are as follows:(1)
Considering various process constraints and terminal constraints and selecting the maximum velocity as the performance index, we can establish the optimal model according to the dynamic equations of interceptor.(2)
After further differentiating the canonical equations and coupled equations to second order, we can get the expression of feedback controlδu and the TPBVP equations of δx,δλ.(3)
After further differentiating the boundary conditions at unspecified terminal time to second order, we can get the expressions ofδx,δλ in terms of free variables dμ and terminal constraints dψf; the relationship between δx0,δλ0 and δxf,δλf is established by backward recursion, which can be used to deduce the expressions of δx0,δλ0 in terms of dμ,dψf, and then we get the expressions of δλf,dtf in terms of δx0,dψf.(4)
The flight time of original nominal trajectorytf plus dtf can be used as new flight time of the new nominal trajectory which can be obtained by interpolating the original nominal trajectory.(5)
Regarding current states of the interceptor as initial constraint andδλf as perturbation, we solve the feedback control δu using the improved IRPM and finally realize online trajectory optimization of midcourse guidance.Figure 2
Online neighboring optimal modified algorithm.


## 5. Simulations and Discussion

In order to verify the validity of the proposed algorithm, the following assumption of intercept operations is designed. Because the nonpropulsive phase is more difficult for trajectory adjustment, the terminal PIP is modified at the initial time of the nonpropulsive phase with altitude increasing 2 km and flight path angle increasing 5 deg. And the current tracking errors of flight path angle and altitude are 0.1 deg and 10 m.The GPM has high precision in solving optimization problems. So when the tracking error generates or the terminal constraints change, the reoptimized trajectory by GPM, which is solved by MATLAB programs GPOPS [25], will be used as a standard optimal modified trajectory for comparison. The comparison of modified trajectory is shown in Figure 3. NOC represents the modified trajectory generated based on NOC theory and the improved IRPM solving algorithm designed in this paper. Nominal represents the nominal trajectory with the unchanged terminal constraints; GPOPS represents the modified trajectory reoptimized by GPM.Figure 3
Curves of altitude and flight path angle.All simulation results of this paper are obtained on a Lenovo laptop with an Intel 2.50 GHz quad-core processor, using windows 7 operating system. All the codes are run under the MATLAB® R2014a environment.From Figure3 we know that the modified trajectory (NOC), which is substantially overlapping the reoptimized trajectory (GPOPS), can satisfy the changed terminal constraints caused by the changed PIP. And the flight time of modified trajectory is approximatively equal to that of the reoptimized trajectory by GPM. These results prove the high accuracy of the algorithm designed in this paper. In the initial time after modification, the error of angle of attack between the modified trajectory and the nominal trajectory is not large, no more than 0.5 deg, which represents that the angle of attack does not have a large jump and therefore ensures the stability of the control.In order to verify the higher precision and efficiency of improved IRPM, the contrast simulation is designed as follows: in the fourth section, according to the tracking error and terminal constraints modifications,δxf and δλf are obtained by backward recursive algorithm. And then we solve feedback control variables by inversely integrating differential equations (19)–(21) of trajectory modification model rather than improved IRPM. In the integration process, 50, 500, and 1000 nodes are, respectively, selected at the same time interval. The curves of flight path angle modification and altitude modification are shown in Figures 4–6; the comparison of time consumption and terminal constraints modifications is shown in Table 2.Table 2
Comparison of time consumption terminal modification.
Method
GPOPS
NOC
50
500
1000
t/s
7.116
0.049
0.144
1.248
2.532
δ
θ
f
/
d
e
g
5.00
4.994
4.081
4.667
4.707
δ
y
f
/
m
2000
1924
1391
2376
2144Figure 4
Comparison of flight path angle modification.Figure 5
Comparison of altitude modification.Figure 6
Flight path angle modifications with different terminals.From Figures4 and 5 and Table 2 we know that compared with the modified trajectory solved by improved IGPM with that reoptimized by GPM, the most part of curves is close except for initial several seconds after modification with flight path angle error of no more than 0.15 degrees and altitude error of no more than 120 m. Furthermore, the two methods have similar precision in terminal modifications, but the time consumption of the former is much shorter than the latter. The method of solving differential equations directly has relatively large errors compared to reoptimization by GPM. With the increase of the number of nodes, the errors decrease gradually and the precision of terminal modifications accordingly increases. But the time consumption increases as well. And even if 1000 nodes are selected, the precision of terminal modifications still has a certain gap compared with improved IGPM while time consumption increases a lot.In summary, the concept of modified trajectory design using nominal trajectory information can guarantee the modified precision with improving efficiency a lot; and the solving method based on improved IGPM method is superior to the method by solving differential equations directly in accuracy and efficiency.In order to verify the adaptivity of modified algorithm, the simulation is designed as follows: the value ofδy is set randomly within -30003000, and the value of δθ is set from the −5 deg to 5 deg at the interval of 1 deg, which adds up to 11 sets of different changed terminal constraints. The results are shown in Figure 7. Similarly, the value of δθ is set randomly within -55, the value of δy is set from the −3000 m to 3000 at the interval of 600 m which adds up to 11 sets of different changed terminal constraints. The results are shown in Figure 7. We can know from Figures 6 and 7 that the algorithm has the modified ability for different changes of terminal constraints. All initial errors of altitude are close to 0 and the maximal initial error of flight path angle is no more than 0.1 deg, which can be adjusted by the trajectory tracking control process.Figure 7
Altitude modifications with different terminals.δ
θ is set as 10 deg and δy is set to 0; the curve of modifications is shown in Figure 8. Similarly, δy is set to 8000 m and δθ is set to 0; the curve of modifications is shown in Figure 9. According to Figures 8 and 9 we can know that when terminal constraints change largely, the error between the modified trajectory generated by the algorithm of this paper and the reoptimized trajectory by GPM is relatively large, which means the optimality of modified trajectory is difficult to guarantee. Because the modified algorithm of this paper is based on the NOC theory, when the changes of terminal constraints are beyond the neighborhood of nominal trajectory, the algorithm is no longer applicable.Figure 8
Curves of modifications with flight path angle increasing 10 deg.Figure 9
Curves of modifications with altitude increasing 8000 m.


## 6. Conclusion

An online trajectory optimization algorithm of midcourse guidance for hypersonic interception is designed based on the NOC theory and the improved IGPM. The main conclusions are as follows:(1)
The algorithm considers not only the tracking error but also the change of terminal constraints. Therefore, it has the disturbance rejection ability and the ability to adjust terminal constraints online.(2)
By using the original nominal trajectory information, a modified trajectory with higher accuracy can be rapidly generated based on neighboring optimal control theory. The iterative optimization process is avoided, and the trajectory generation is realized online.(3)
When the changes of terminal constraints are relatively large, beyond the neighborhood of the original nominal trajectory, the online modified algorithm is no longer applicable. It needs to reoptimize online to get the modified trajectory.



---

*Source: 9456179-2017-07-12.xml*