# A Rapid Penetration Trajectory Optimization Method for Hypersonic Vehicles

Authors: Binbin Yan; Ruifan Liu; Pei Dai; Muzeng Xing; Shuangxi Liu
Journal: International Journal of Aerospace Engineering
 (2019)

---

## Abstract

The evasion maneuver problem of hypersonic vehicles differs from those of ballistic missiles and other traditional weapons, showing distinctive properties including expansive maneuver range and weak maneuverability. How to avoid the disadvantage of low available overload and ensure follow-up tasks are the main concerns of the hypersonic penetration. This paper presents a penetration trajectory optimization algorithm for an air-breathing hypersonic vehicle, where the prerequisite penetration condition is analyzed and control costs are chosen as an objective function to minimize the fuel consumption and maneuver range. This paper focuses on how to formulate the complex, highly constrained nonconvex penetration problem to be a sequence of easily solved second-order cone programming through a combination of successive linearization and relaxation techniques. Innovation lies in the raising of the penetration angle and the relaxation technique of nonlinear and nonconvex elements. Various numerical simulations are conducted to verify the validity of penetration condition and to demonstrate that the proposed method is effective and has a good computational performance irrespective of initial guesses.


---

## Body

## 1. Introduction

Since the hypersonic vehicles have arrived, the attack and defense countermeasure game centered on it receives the focal attention by various military powers. The air-breathing hypersonic vehicle cruises at the altitude of the near space, featuring a lasting and wide-range maneuverability. But owing to the usage of the unique integrated propulsion technology and the scramjet engine [1], strict constraints are also imposed on its maneuverability, including the attack angle and bank angle [2, 3]. Therefore, the hypersonic penetration possesses characteristics of a “wide maneuvering range, weak maneuverability, and high flight speed.” Besides, the prominent infrared features of hypersonic vehicles make it easily to be detected and intercepted. Consequently, how to make maximal use of the speed superiority and achieve penetration with minimal energy consumption are the core of hypersonic penetration programming.In the confrontation, the hypersonic vehicle and the interception missile maneuver actively and simultaneously, thus always regarded as a complicated gaming problem. With the game theory, differential gaming models are established. Shaferman and Shima [4, 5] studied the zero-sum differential game theory, the quadratic differential game theory, and various other forms of differential game and analyzed the conditions in which the saddle points of the differential game exist. However, most of studies of differential game theory focus on mathematical modeling and theoretical derivation, while its application to practical engineering remains a challenge [6]. Considering the huge speed gap, the game problem between a hypersonic vehicle and an interception missile is not like the normal one. In other words, only when harsh interception conditions are satisfied can the interceptor succeed. Accordingly, Dwivedi et al. [7] presented a midcourse guidance law against high-speed targets believing that alignment angle constraints in both elevation and azimuth could create a favorable condition for the terminal guidance while Guo et al. [8] tried to guide a hypersonic vehicle to break the alignment angle of the interceptor from an opponent view; thus, they converted the penetration problem into an optimal control problem.In fact, using optimal control techniques for guidance of flight vehicles is not new. Numerical methods, like the shooting method [9], the gradient method [10], and the recently prominent one, the pseudospectral method [11, 12], have been developed to solve two-point boundary value problems (TPBVPs) derived from optimal problems. However, all of them lead to large computations that are infeasible to implement in real time. Other contributions have also been proposed in this regard, such as the model predictive control (MPC) theory [13] combined with approximate dynamic programming (APC) and model predictive static programming (MPSP) [7, 8]. The former one has been successfully used in the industry (especially in slow-varying process control) but is not guaranteed to remain optimal and stable for its sliding time window [13, 14]. As to MPSP which is rather efficient and adaptable, however, when applied to the hypersonic vehicle, it is also not a preferable method since it cannot properly solve the optimal problems with control constraints.This paper is aimed at designing a rapid penetration method to generate an appropriate evasion trajectory of a hypersonic vehicle against a head-on intercepting missile. According to the above content, the main two barriers of this problem are an explicit expression of penetration condition and a rapid, robust optimization algorithm catering both terminal angular constraints and the demanding control constraints of hypersonic vehicles. About the former, this paper introduces the concept of the penetration angle on the basis of interception geometry, giving the clear formula of penetration condition. The validity of the condition is also verified by the numerical simulation in terms of the ultimate miss distance under different penetration angles. As to the feasible optimization method, this paper raises a method based on sequential convex programming, which could meet the control constraints as well as the penetration angle constraints in a short time, and compared with the existing optimal method GPM, the algorithm proposed in this paper has a very competitive convergence rate and does not require good initial guesses.Convex optimization has been widely utilized in recent years for its unique theoretical advantages and the polynomial time complexity [15, 16]. However, applications of convex optimization in aerospace trajectory planning have been limited to problems without complex nonlinear dynamics and aerodynamic forces so far. This paper is mainly devoted to this difficulty in the penetration application. The technique developed in this paper successfully “convexifies” the initial nonconvex penetration problem into a series of second-order cone programming (SOCP) problems. Various convexity strategies are applied to different nonconvex elements. The SOCP-based method proves to be a very robust and rapid tool for hypersonic penetration trajectory optimization. Numerical results are provided to demonstrate the effectiveness and some of the noteworthy aspects of the method.The rest of this paper is organized as follows. In Section2, the countermeasure model is established, including the hypersonic vehicle model, the interceptor model, and the relative motion model, and thus, the original form of penetration problem is described. In Section 3, the convexification process is explained in detail and the converted SOCP problem is presented. Numerical simulation results are presented in Section 4. Finally, some conclusions are drawn in Section 5.


## 2. Mathematical Model for Attack and Defense Countermeasure

The primary objective of this paper is to present an effective penetration guidance logic. Hence, considering solving speed and application feasibility, one only requires point mass mathematical models for both the vehicle and the interception missile. Moreover, with the complex high-fidelity model, the calculation of Jacobi matrices presented in Section3 will be boosted. Therefore, it is crucial to formulate the penetration problem in a certain simplified way conducive to the SOCP approach. This section presents the details of one such formulation of the problem.


### 2.1. Dynamics of the Hypersonic Vehicle

The penetration maneuver of a hypersonic vehicle is a short-range trajectory optimization problem; for this reason, we simplify its kinematic equation sets as follows:(a)
With the rotation of the Earth not taken into account, namely, ignoring the Coriolis acceleration and convective acceleration
(b)
Assuming that the Earth is flat, namely, ignoring the influence of its centrifugal forceAnd as is mentioned above, this paper pays most of its attention to the trajectory programming; therefore, the mass point model is established and some simplifications of the vehicle’s thrust and aerodynamic forces are conducted.Thus, we obtain the following equations of the mass point motion model of the hypersonic vehicle:(1)ḣe=Vesinθe,V̇e=TcosαT−Dm−gsinθe,θ̇e=TsinαTcosσ+LTcosσmVe−gcosθeVe,ψ̇e=TsinαTsinσ+LTsinσmVecosθe,ṁe=−ms,where he is the altitude of the hypersonic vehicle; Ve is the speed; θe and ψe represent the flight path angle and trajectory deflection angle; αT and σ are the angle of attack and bank angle; LT,D, and T are its aerodynamic force and thrust; me is its mass; ms is its mass consumption rate.The thrust model, the aerodynamic force model, and the control constraints are as follows:(2)LT=0.5ρVe2SCL0+CLααT,D=0.5ρVe2SCD0+CDα2αT2,T=0.5ρVe2SCTmsms,αTmin≤αT≤αTmax,σ̇≤σ̇max,0≤ms≤msmax,where αTmin and αTmax represent the boundaries of attack of the angle while σ̇max and msmax represent the upper limits of the bank angle rate and mass consumption rate, respectively.


### 2.2. Dynamics of the Interception Missile

For the midcourse guidance of an interception missile against a high-speed target, it is desirable to shape the missile trajectory that provides a preferable initial state for terminal attack. Instead of the guidance laws carried out according to the nonlinear theory and the MPSP algorithm, whose practicability are always restricted in application, this paper utilizes a more compendious and feasible guidance law called the trajectory shaping guidance law [17, 18].The trajectory shaping guidance law has been developed to travel to the desired trajectory and verified by physical applications, such as landing the Apollo spacecraft on the moon. The guidance law appears to be a form of augmented proportional navigation plus an extra term that is proportional to the difference between the line of sight angle and the desired angle at the end of flight. The application formula for this paper is as follows:(3)θ̇i=N1q̇y+N2tgoqy+θe,ψ̇i=N1q̇z+N3tgoqz−ψe+π,where θi and ψi are the flight path angle and the trajectory deflection angle; qy and qz are the line of sight (LOS) angle in elevation and azimuth, respectively. Figure 1 shows explicit geometric meaning; N1,N2,N3 are the constants of the guidance law, where tgo is the remaining flight time.Figure 1
Relative kinetics in the attack and defense countermeasure.


### 2.3. Attack and Defense Countermeasure Analysis

Under usual circumstances, the speed of a hypersonic vehicle is much higher than that of a defense system. Hence, only the head-on strike by an interception missile can impose a threat to the hypersonic vehicle. Through reasonable assumptions and simplifications, the interception geometry is established, and thus, one can gain the interception missile’s attack area. The attack area is a round spherical surface that includes the interception missile but excludes the hypersonic vehicle, and its boundary is tangent to the hypersonic vehicle, as shown in Figure2. See reference [19] for specific derivation processes.Figure 2
The three-dimensional interception attack area.It is known from interception geometry that if the velocity vector of the hypersonic vehicle is always outside of the interception geometric tangent line, namely,η>ηmax, then it can break through the head-on strike situation of the interception missile, achieving a successful penetration [20].η is defined as the penetration angle of the hypersonic vehicle, namely, the angle between the line-of-sight and the velocity of the hypersonic vehicle. Its cosine value can be obtained with the transformation matrices of the ground-trajectory coordinate system and the ground-LOS coordinate system. The relationship can be formulated as follows:
(4)cosη=−cosqycosθecosqz−ψe−sinθesinqy.


### 2.4. Mathematical Model of Relative Motions

With the penetration angle constraints taken into consideration, the model of relative motion between the hypersonic vehicle and the interception missile is established under the LOS coordinate system. The spatial coordinates and relative kinetics in the attack and defense countermeasure are shown in Figure1.The relative motions are defined by the following equations:(5)ṙ=Vecosθecosqycosψe−qz+sinθesinqy−Vicosθicosqycosψi−qz+sinθisinqy,q̇yr=−Vecosθesinqycosψe−qz−sinθecosqy+Vicosθisinqycosψi−qz−sinθicosqy,q̇zrcosqy=Vecosθesinψe−qz−Vicosθisinψi−qz,where r,qy,qz are the relative distance, the LOS elevation angle, and the LOS azimuth angle between the hypersonic vehicle and the interception missile.


### 2.5. Hypersonic Vehicle’s Penetration Problem

Taking the minimum control as the objective function and the penetration angle constraint into consideration, the hypersonic vehicle’s penetration problem can be described as an optimal control problemP0:
(6)P0:minJ=∫utTRut,s.t.ẋt=fxt,ut,xt0=x0,αTmin≤αTt≤αTmax,σ̇t≤σ̇max,0≤mst≤msmax,ηt>ηmin,where the control vector is u=αT,σ,msT∈R3; x=he,Ve,θe,ψe,me,θi,ψi,r,qy,qzT∈R11 represents the state vector.


## 2.1. Dynamics of the Hypersonic Vehicle

The penetration maneuver of a hypersonic vehicle is a short-range trajectory optimization problem; for this reason, we simplify its kinematic equation sets as follows:(a)
With the rotation of the Earth not taken into account, namely, ignoring the Coriolis acceleration and convective acceleration
(b)
Assuming that the Earth is flat, namely, ignoring the influence of its centrifugal forceAnd as is mentioned above, this paper pays most of its attention to the trajectory programming; therefore, the mass point model is established and some simplifications of the vehicle’s thrust and aerodynamic forces are conducted.Thus, we obtain the following equations of the mass point motion model of the hypersonic vehicle:(1)ḣe=Vesinθe,V̇e=TcosαT−Dm−gsinθe,θ̇e=TsinαTcosσ+LTcosσmVe−gcosθeVe,ψ̇e=TsinαTsinσ+LTsinσmVecosθe,ṁe=−ms,where he is the altitude of the hypersonic vehicle; Ve is the speed; θe and ψe represent the flight path angle and trajectory deflection angle; αT and σ are the angle of attack and bank angle; LT,D, and T are its aerodynamic force and thrust; me is its mass; ms is its mass consumption rate.The thrust model, the aerodynamic force model, and the control constraints are as follows:(2)LT=0.5ρVe2SCL0+CLααT,D=0.5ρVe2SCD0+CDα2αT2,T=0.5ρVe2SCTmsms,αTmin≤αT≤αTmax,σ̇≤σ̇max,0≤ms≤msmax,where αTmin and αTmax represent the boundaries of attack of the angle while σ̇max and msmax represent the upper limits of the bank angle rate and mass consumption rate, respectively.


## 2.2. Dynamics of the Interception Missile

For the midcourse guidance of an interception missile against a high-speed target, it is desirable to shape the missile trajectory that provides a preferable initial state for terminal attack. Instead of the guidance laws carried out according to the nonlinear theory and the MPSP algorithm, whose practicability are always restricted in application, this paper utilizes a more compendious and feasible guidance law called the trajectory shaping guidance law [17, 18].The trajectory shaping guidance law has been developed to travel to the desired trajectory and verified by physical applications, such as landing the Apollo spacecraft on the moon. The guidance law appears to be a form of augmented proportional navigation plus an extra term that is proportional to the difference between the line of sight angle and the desired angle at the end of flight. The application formula for this paper is as follows:(3)θ̇i=N1q̇y+N2tgoqy+θe,ψ̇i=N1q̇z+N3tgoqz−ψe+π,where θi and ψi are the flight path angle and the trajectory deflection angle; qy and qz are the line of sight (LOS) angle in elevation and azimuth, respectively. Figure 1 shows explicit geometric meaning; N1,N2,N3 are the constants of the guidance law, where tgo is the remaining flight time.Figure 1
Relative kinetics in the attack and defense countermeasure.


## 2.3. Attack and Defense Countermeasure Analysis

Under usual circumstances, the speed of a hypersonic vehicle is much higher than that of a defense system. Hence, only the head-on strike by an interception missile can impose a threat to the hypersonic vehicle. Through reasonable assumptions and simplifications, the interception geometry is established, and thus, one can gain the interception missile’s attack area. The attack area is a round spherical surface that includes the interception missile but excludes the hypersonic vehicle, and its boundary is tangent to the hypersonic vehicle, as shown in Figure2. See reference [19] for specific derivation processes.Figure 2
The three-dimensional interception attack area.It is known from interception geometry that if the velocity vector of the hypersonic vehicle is always outside of the interception geometric tangent line, namely,η>ηmax, then it can break through the head-on strike situation of the interception missile, achieving a successful penetration [20].η is defined as the penetration angle of the hypersonic vehicle, namely, the angle between the line-of-sight and the velocity of the hypersonic vehicle. Its cosine value can be obtained with the transformation matrices of the ground-trajectory coordinate system and the ground-LOS coordinate system. The relationship can be formulated as follows:
(4)cosη=−cosqycosθecosqz−ψe−sinθesinqy.


## 2.4. Mathematical Model of Relative Motions

With the penetration angle constraints taken into consideration, the model of relative motion between the hypersonic vehicle and the interception missile is established under the LOS coordinate system. The spatial coordinates and relative kinetics in the attack and defense countermeasure are shown in Figure1.The relative motions are defined by the following equations:(5)ṙ=Vecosθecosqycosψe−qz+sinθesinqy−Vicosθicosqycosψi−qz+sinθisinqy,q̇yr=−Vecosθesinqycosψe−qz−sinθecosqy+Vicosθisinqycosψi−qz−sinθicosqy,q̇zrcosqy=Vecosθesinψe−qz−Vicosθisinψi−qz,where r,qy,qz are the relative distance, the LOS elevation angle, and the LOS azimuth angle between the hypersonic vehicle and the interception missile.


## 2.5. Hypersonic Vehicle’s Penetration Problem

Taking the minimum control as the objective function and the penetration angle constraint into consideration, the hypersonic vehicle’s penetration problem can be described as an optimal control problemP0:
(6)P0:minJ=∫utTRut,s.t.ẋt=fxt,ut,xt0=x0,αTmin≤αTt≤αTmax,σ̇t≤σ̇max,0≤mst≤msmax,ηt>ηmin,where the control vector is u=αT,σ,msT∈R3; x=he,Ve,θe,ψe,me,θi,ψi,r,qy,qzT∈R11 represents the state vector.


## 3. Penetration Problem Convexification

The hypersonic penetration problemP0 has nonlinear inequation constraints and other nonconvex factors brought about by the quadratic integral objective function, nonlinear differential equation constraints, and penetration angle constraints. Therefore, how to transform the penetration problem into a convex problem that can be easily solved is crucial for applying a convex optimization algorithm. The following section will explain in some detail the processes.


### 3.1. Linearization and Discretization of Dynamics Equations

The nonlinear differential equations are approximated by a series of linear time-varying systems, and the solution sequence of linear systems, through several iterations, will be convergent to the optimal solution of the original system [21].For simplicity of notation, the nonlinear state equation is expressed as follows:(7)ẋ=fx,u,xt0=x0.In the form of state correlation coefficient, the nonlinear state equations can be presented as follows:(8)ẋ=Fxx,ux+Fux,uu+bx,u,bx,u=fx,u−Fxx,ux−Fux,uu,where Fx=∂f/∂x and Fu=∂f/∂u are the Jacobi matrices of the state function f in terms of x and u. If the above equations meet the local Lipschitz condition, then the solution of the nonlinear state equation can be approximated with the following linear time-varying equation:
(9)ẋk+1=Fxxk,ukxk+1+Fuxk,ukuk+1+bxk,uk,k=0,1,⋯,where xk and uk denote the values of the state variable and control variable in the k times of iteration.The sequences of solutions obtained through solving the above linear time-varying equations converge to the solutions of original nonlinear equations [22], namely,
(10)limk→∞xkt→xt.To ensure the effectiveness of state equation linearization, we add thrust domain constraints to state variables.(11)x−xk≤δ,where xk denotes the solution of state variables in the k times of iteration; δ∈R11 denotes the trust threshold value and is a constant vector.What follows is the discretization of the linearized state equation, but before discretization, it is imperative to confirm the optimization range and discrete points. We choose the relative distancer as the independent variable of the hypersonic vehicle’s penetration problem under investigation in this paper. The work has the following advantages: (1) compared with the separate state of the hypersonic vehicle and the interception missile at a certain moment, we are more concerned about the attack and defense countermeasure in terms of a certain relative distance and (2) if we use time as an independent variable, the appearance of uncertain terminal time tf will make the optimal control problem more complicated. With the relative distance served as the independent variable, the final value of r can be assumed to be the ending point of midcourse guidance of the interception missile. Thus, the terminal penetration angle constraint can also be determined.Accordingly,r0 and rf are chosen as the optimization range of the penetration problem and the problem is discretized by N+1 numbers of discrete points whose intervals are equal. Thus, the discrete step length is Δr=rf−r0/N. Firstly, we discretize the state variable and the control variable on discrete points, namely, xi=xri,ui=uri,i=0,1,⋯,N. Then according to the trapezoidal formula, we perform the mathematical integration of the sequentially linearized state equation (namely, equation (9)) with the following formula:
(12)xi=xi−1+Δr2Ai−1kxi−1+Bi−1kui−1+bi−1k+Aikxi+Bikui+bik,where Aik≜Fxxkri,ukri, Bik≜Fuxkri,ukri, and bik=bxkri,ukri.Equation (12) is further represented as follows:
(13)Hi−1xi−1+Hixi+Gi−1ui−1+Giui=−Δr2bi−1k+bik,where Hi−1=I+Δr/2Ai−1k, Hi=−I+Δr/2Aik, Gi−1=Δr/2Bi−1k, and Gi=Δr/2Bik; the number of dimensions of unit matrix I is in agreement with that of coefficient matrix A.All the variables to be optimized are unitarily represented as optimized vectorz=x0T⋯xNTu0T⋯uNTT; the linearly discretized state equations are written into the expressions related to z:
(14)Mz=F,where
(15)M=I00⋯00000⋯00H0H10⋯00G0G10⋯00⋮⋮⋮⋱⋮⋮⋮⋮⋮⋱⋮⋮000⋯HN−1HN000⋯GN−1GN,F=−Δe2−2Δex0b0k+b1k⋮bN−1k+bNk.It should be pointed out that the first line of equation in the equality constraints (namely, equation (14)) obtained through linearizing the initial state constraint xr0=x0.


### 3.2. Control Constraint Convexification

The original problemP0 has control constraints with the angle of attack, bank angle, and mass consumption rate. The constraints on both the angle of attack aT and the mass consumption rate ms are all linear constraints while the constraint on the bank angle is nonlinear and handled in the following ways: the bank angle rate σ̇ is introduced as a control variable, namely, u=αT,σ̇,msT∈R3, and the state variable σ is also added to the differential function.


### 3.3. Penetration Angle Constraint Convexification

The interception missile is assumed to use the head-on interception method, and the anglesqy+θe and π+qz−ψe formed by the missile’s line of sight and the hypersonic vehicle’s velocity vector are less than 45 degrees. Therefore, we perform the following slackening on formula (4):
(16)cosη=−cosqycosθecosqz−ψe−sinθesinqy<cosqy+θecosπ+qz−ψe.The nonlinear inequality constraint is then transformed into the concave inequality constraint:(17)cosqy+θecosπ+qz−ψe<cosηmax.To show that equation (17) is a concave inequation constraint, we further transform the left hand-side term of the inequation into the following:
(18)cosqy+θecosπ+qz−ψe=cosqy+θe+qz−ψe+π+cosqy+θe−qz+ψe−π2.According to the abovementioned assumptionsqy+θe∈−π/4,π/4 and π+qz−ψe∈−π/4,π/4, the value ranges of the two cosine functions at the left side of equation (18) are all in −π/2,π/2; within this scope, the cosine function cosx is a concave function. Refer to the compound affine mapping theorem and the vector compound theorem [23] in the convex optimization theory. Equation (17) is deduced to be a concave inequation constraint.Similar to the differential equation convexification method, the concave inequation constraint can also be approximated with the sequential linearization; the sequential linearization of a concave inequation has already been proved to be able to converge to its optimal solutions [24]. The formula for sequentially linearizing the penetration angle constraint is as follows:
(19)∂g∂qyxkqy+∂g∂qzxkqz+∂g∂θexkθe+∂g∂ψexkψe<cosηmax+b¯,where
(20)∂g∂qyxk=∂g∂θexk=−sinqyk+θekcosπ+qzk−ψek,∂g∂qzxk=−sinqyk+θekcosπ+qzk−ψek,∂g∂ψexk=sinqyk+θekcosπ+qzk−ψek.qyk, qzk, θek, and ψek denote the solutions of the k times of iteration; the linearized residual error is as follows:
(21)b¯=−cosqyk+θekcosπ+qzk−ψek+∂g∂qyxkqyk+∂g∂qzxkqzk+∂g∂θexkθek+∂g∂ψexkψek.Thus, the penetration angle constraint is represented in the form of terminal state constraint qualification, namely,(22)Cxfkxf+dx<0,where the coefficient matrix is obtained with sequential linearization equation (19).


### 3.4. Objective Convexification

Through introducing the intermediate variableξ, the quadratic performance index in problem P0 can be equivalently discretized into the following:
(23)minJ=∑i=1Nζi,R1/2ui≤ζi,i=1,⋯,N,where R1/2 is obtained with the square root matrix decomposition formula, namely, R=R1/2TR1/2.So far, through slackening, sequentially linearizing, and discretizing objective functions, differential equations, and nonlinear inequations, the original penetration problemP0 is transformed into a convex problem in its discrete form, noted as problem P1:
(24)P1:minJ=∑i=1Nζi,s.t.Mz=F,αTmin≤αT_i≤αTmax,σ̇i≤σ̇max,0≤ms_i≤msmax,Cxfkxf+dx<0,R1/2ui≤ζi,xi−xik≤δ,where · is a 2 norm; i=1,2,⋯,N.ProblemP1 is a sort of convex problem named second-order cone problem (SOCP), with a linear objective function and all convex constraints including linear constraints and second-order cone constraints. The optimal solution of the original penetration problem P0 can be obtained through iteratively solving problem P1. Each iteration process uses solution xk of the last iteration to update problem P1, constantly solves problem P1 within a thrust domain, and gradually approximates the optimal solution of problem P0 till maxixik+1−xik≤ε. Then the optimization algorithm converges and obtains optimal solutions.


## 3.1. Linearization and Discretization of Dynamics Equations

The nonlinear differential equations are approximated by a series of linear time-varying systems, and the solution sequence of linear systems, through several iterations, will be convergent to the optimal solution of the original system [21].For simplicity of notation, the nonlinear state equation is expressed as follows:(7)ẋ=fx,u,xt0=x0.In the form of state correlation coefficient, the nonlinear state equations can be presented as follows:(8)ẋ=Fxx,ux+Fux,uu+bx,u,bx,u=fx,u−Fxx,ux−Fux,uu,where Fx=∂f/∂x and Fu=∂f/∂u are the Jacobi matrices of the state function f in terms of x and u. If the above equations meet the local Lipschitz condition, then the solution of the nonlinear state equation can be approximated with the following linear time-varying equation:
(9)ẋk+1=Fxxk,ukxk+1+Fuxk,ukuk+1+bxk,uk,k=0,1,⋯,where xk and uk denote the values of the state variable and control variable in the k times of iteration.The sequences of solutions obtained through solving the above linear time-varying equations converge to the solutions of original nonlinear equations [22], namely,
(10)limk→∞xkt→xt.To ensure the effectiveness of state equation linearization, we add thrust domain constraints to state variables.(11)x−xk≤δ,where xk denotes the solution of state variables in the k times of iteration; δ∈R11 denotes the trust threshold value and is a constant vector.What follows is the discretization of the linearized state equation, but before discretization, it is imperative to confirm the optimization range and discrete points. We choose the relative distancer as the independent variable of the hypersonic vehicle’s penetration problem under investigation in this paper. The work has the following advantages: (1) compared with the separate state of the hypersonic vehicle and the interception missile at a certain moment, we are more concerned about the attack and defense countermeasure in terms of a certain relative distance and (2) if we use time as an independent variable, the appearance of uncertain terminal time tf will make the optimal control problem more complicated. With the relative distance served as the independent variable, the final value of r can be assumed to be the ending point of midcourse guidance of the interception missile. Thus, the terminal penetration angle constraint can also be determined.Accordingly,r0 and rf are chosen as the optimization range of the penetration problem and the problem is discretized by N+1 numbers of discrete points whose intervals are equal. Thus, the discrete step length is Δr=rf−r0/N. Firstly, we discretize the state variable and the control variable on discrete points, namely, xi=xri,ui=uri,i=0,1,⋯,N. Then according to the trapezoidal formula, we perform the mathematical integration of the sequentially linearized state equation (namely, equation (9)) with the following formula:
(12)xi=xi−1+Δr2Ai−1kxi−1+Bi−1kui−1+bi−1k+Aikxi+Bikui+bik,where Aik≜Fxxkri,ukri, Bik≜Fuxkri,ukri, and bik=bxkri,ukri.Equation (12) is further represented as follows:
(13)Hi−1xi−1+Hixi+Gi−1ui−1+Giui=−Δr2bi−1k+bik,where Hi−1=I+Δr/2Ai−1k, Hi=−I+Δr/2Aik, Gi−1=Δr/2Bi−1k, and Gi=Δr/2Bik; the number of dimensions of unit matrix I is in agreement with that of coefficient matrix A.All the variables to be optimized are unitarily represented as optimized vectorz=x0T⋯xNTu0T⋯uNTT; the linearly discretized state equations are written into the expressions related to z:
(14)Mz=F,where
(15)M=I00⋯00000⋯00H0H10⋯00G0G10⋯00⋮⋮⋮⋱⋮⋮⋮⋮⋮⋱⋮⋮000⋯HN−1HN000⋯GN−1GN,F=−Δe2−2Δex0b0k+b1k⋮bN−1k+bNk.It should be pointed out that the first line of equation in the equality constraints (namely, equation (14)) obtained through linearizing the initial state constraint xr0=x0.


## 3.2. Control Constraint Convexification

The original problemP0 has control constraints with the angle of attack, bank angle, and mass consumption rate. The constraints on both the angle of attack aT and the mass consumption rate ms are all linear constraints while the constraint on the bank angle is nonlinear and handled in the following ways: the bank angle rate σ̇ is introduced as a control variable, namely, u=αT,σ̇,msT∈R3, and the state variable σ is also added to the differential function.


## 3.3. Penetration Angle Constraint Convexification

The interception missile is assumed to use the head-on interception method, and the anglesqy+θe and π+qz−ψe formed by the missile’s line of sight and the hypersonic vehicle’s velocity vector are less than 45 degrees. Therefore, we perform the following slackening on formula (4):
(16)cosη=−cosqycosθecosqz−ψe−sinθesinqy<cosqy+θecosπ+qz−ψe.The nonlinear inequality constraint is then transformed into the concave inequality constraint:(17)cosqy+θecosπ+qz−ψe<cosηmax.To show that equation (17) is a concave inequation constraint, we further transform the left hand-side term of the inequation into the following:
(18)cosqy+θecosπ+qz−ψe=cosqy+θe+qz−ψe+π+cosqy+θe−qz+ψe−π2.According to the abovementioned assumptionsqy+θe∈−π/4,π/4 and π+qz−ψe∈−π/4,π/4, the value ranges of the two cosine functions at the left side of equation (18) are all in −π/2,π/2; within this scope, the cosine function cosx is a concave function. Refer to the compound affine mapping theorem and the vector compound theorem [23] in the convex optimization theory. Equation (17) is deduced to be a concave inequation constraint.Similar to the differential equation convexification method, the concave inequation constraint can also be approximated with the sequential linearization; the sequential linearization of a concave inequation has already been proved to be able to converge to its optimal solutions [24]. The formula for sequentially linearizing the penetration angle constraint is as follows:
(19)∂g∂qyxkqy+∂g∂qzxkqz+∂g∂θexkθe+∂g∂ψexkψe<cosηmax+b¯,where
(20)∂g∂qyxk=∂g∂θexk=−sinqyk+θekcosπ+qzk−ψek,∂g∂qzxk=−sinqyk+θekcosπ+qzk−ψek,∂g∂ψexk=sinqyk+θekcosπ+qzk−ψek.qyk, qzk, θek, and ψek denote the solutions of the k times of iteration; the linearized residual error is as follows:
(21)b¯=−cosqyk+θekcosπ+qzk−ψek+∂g∂qyxkqyk+∂g∂qzxkqzk+∂g∂θexkθek+∂g∂ψexkψek.Thus, the penetration angle constraint is represented in the form of terminal state constraint qualification, namely,(22)Cxfkxf+dx<0,where the coefficient matrix is obtained with sequential linearization equation (19).


## 3.4. Objective Convexification

Through introducing the intermediate variableξ, the quadratic performance index in problem P0 can be equivalently discretized into the following:
(23)minJ=∑i=1Nζi,R1/2ui≤ζi,i=1,⋯,N,where R1/2 is obtained with the square root matrix decomposition formula, namely, R=R1/2TR1/2.So far, through slackening, sequentially linearizing, and discretizing objective functions, differential equations, and nonlinear inequations, the original penetration problemP0 is transformed into a convex problem in its discrete form, noted as problem P1:
(24)P1:minJ=∑i=1Nζi,s.t.Mz=F,αTmin≤αT_i≤αTmax,σ̇i≤σ̇max,0≤ms_i≤msmax,Cxfkxf+dx<0,R1/2ui≤ζi,xi−xik≤δ,where · is a 2 norm; i=1,2,⋯,N.ProblemP1 is a sort of convex problem named second-order cone problem (SOCP), with a linear objective function and all convex constraints including linear constraints and second-order cone constraints. The optimal solution of the original penetration problem P0 can be obtained through iteratively solving problem P1. Each iteration process uses solution xk of the last iteration to update problem P1, constantly solves problem P1 within a thrust domain, and gradually approximates the optimal solution of problem P0 till maxixik+1−xik≤ε. Then the optimization algorithm converges and obtains optimal solutions.


## 4. Numerical Demonstrations

### 4.1. Settings of Simulation Parameters

The basic parameters and the relevant simulation parameter settings in the attack and defense countermeasure simulation model of the hypersonic vehicle are given in Tables1 and 2. The boundary values in the control constraints are set as αTmin=−4°, αTmax=6°, σ̇max=5°, and msmax=1; the boundary value in the penetration angle constraints is set as ηmax=20°.Table 1
The vehicle’s basic parameters [8].
Parameters
CL0
CLα (rad−1)
CD0
CDα2 (rad−2)
CTms (s·kg−1)
Parameter values
-0.013
1.833
0.015
4.596
0.162Table 2
The initial values of conditions of the attack side and the defense side.
Simulation parameters
Parameter values
xe0,ye0,ze0 (km)
100,25,10
Ve0 (m·s−1)
1500
θe0,ψe0 (rad)
0,π
xi0,yi0,zi0 (km)
0,0,0
Vi (m·s−1)
1200
θi0,ψi0 (rad)
(0.2438, -0.0997)
N1,N2,N3
5,0.1,0.1i/N
R
diag5,0.1,1
rf (km)
20Simulation conditions are set as follows: in the first iteration, control guesses are set tou=2°,0,0.2T, the number of discrete points is set as N=200, and the thrust domain constraint values in the formula are set as follows:
(25)δ=200020π18020π180100020π18020π1801000120π18020π18020π180T.The error threshold values in the simulation iteration termination condition are set as follows:(26)ε=1000.5π1800.5π1801000.5π1800.5π18050010.5π1800.5π1800.5π180T.The above simulation conditions are used to perform numerical simulation with the YALMIP modeling environment of MATLAB; the primal-dual interior point algorithm in the MOSEK [25] software package is used to solve the SOCP problem. After five times of sequential iteration, the trajectory converges. The time required by each iteration is around 0.4 s; the resolving time required by the interior point algorithm is around 0.2 s. The penetration angle and the hypersonic vehicle’s angle of attack calculated by the first four iterations are shown in Figures 3 and 4. The figures show that the sequential convex optimization algorithm converges quickly; the angle-of-attack value solved in the second iteration is already close to the ultimate optimal solution.Figure 3
Penetration angle in the first four successive iterations of the SOCP-based algorithm.Figure 4
Angle-of-attack profile in the first four successive iterations of the SOCP-based algorithm.In addition, in order to validate the effectiveness of the second-order cone programming (SOCP) method proposed in the paper, numerical results are compared with those obtained by the Gaussian pseudospectral method (GPM). The GPM solves the problem with the sequential quadratic programming (SQP), and we program it by the GPOPS5.1 software package [26].The GPM has already been proved to obtain its globally optimal solution [27].The trajectory profile by the SOCP-based method and the GPM (the mesh tolerance is set to 10-3) are shown in Figures 5–7, while the solving times required and optimal performance indexes are given in Table 3. Figure 8 presents the required calculating time of the two methods using different random control initial guesses within control constraints; the horizontal axis only represents the number of initial guesses without any specific physical meaning.Figure 5
The penetrator’s attack angle profiles, using the two methods.Figure 6
The penetrator’s banking angle profiles, using the two methods.Figure 7
The penetrator’s fuel consumption rate profiles, using the two methods.Table 3
Comparison of the time consumed with performance indexes, using the SOCP and GPM methods.
Method
Discrete point
Iteration
Time (s)
Performance index
GPM
40
—
7.76
751
SOCP
200
5
1.42
770Figure 8
Calculating time for different initial guesses.Based on the figures, the SOCP and GPM solutions show very similar trends; the control profiles in Figures5–7 match well in general. The relatively small differences between the two methods are attributable to the different discretization strategies used in the two solutions, which also occur in reference [28]. On the other hand, the performance of the GPM is only slightly optimal than that of SOCP in terms of the values of 751 and 770, respectively. Therefore, we can at least reckon that the SOCP method is near optimal. Despite the compromise of the optimum objective, the SOCP method performs much better in terms of computational complexity based on Table 3 and Figures 4 and 8. The SOCP-based method in this paper takes an average of 1.3 seconds to generate trajectories, and the peak calculating time for different initial guesses is only 2.45 seconds, while the GPM method takes an average of 7.98 seconds and the longest time can reach 13.1 seconds. To tackle penetration problems, a method with competitive convergence speed and robust solution process is what we always seek for, though it is not rigorously optimal.Further simulation is conducted to verify whether the penetration angle could function as an index for evaluating the penetration effects when the distance between the missile and the hypersonic vehicle is fixed. In other words, the hypersonic vehicle can successfully escape once it meets certain penetration angle constraints in the attack and defense countermeasure.Therefore, we design different penetration angle constraintsη=10°,15°,20°,25°,30° and conduct the simulation in which the continuing terminal-course relative motion is considered. In the simulation, when the distance between the attack side and the defense side equals to rf, the countermeasure enters the terminal phase where the two sides continue to approach each other till the vehicle is intercepted or successfully escapes. During this period, the guidance law of the interception missile is the pure proportional guidance law with the proportional constant of 3, while the vehicle adopts two kinds of maneuver in order to illustrate the impact of angle constraints. It does not make any maneuver and continues flying at the current speed and direction or maneuvers with maximum overload in the opposite direction of the attacked missile. The maximum overload, about 2.8g, is obtained by the maximum available attack angle and maximum thrust for a hypersonic vehicle. The formulas for calculating its overload and maneuver direction are as follows:
(27)n¯cmax=ncmaxI¯n,ncmax=TmaxsinαTmax+LTαTmaxm,I¯n=V¯T×Los¯×V¯TV¯T×Los¯×V¯T,where I¯n represents the maneuver direction; V¯T and Los¯ represent the vehicle’s velocity vector and line-of-sight vector, respectively. With various factors being considered, the time constant of an interception missile is set as 0.05 s and the miss distance is designed as 1 m.Under different penetration angle constraints and two scenarios of the hypersonic vehicle’s continuing motion method, the miss distances of the interception missile in the two scenarios are shown in Table4. The whole trajectory profiles for the attack side and the defense side are given in Figures 9 and 10, while the figure for comparing the terminal trajectory (20°penetration angle) in the horizonal plane between different maneuver methods is also shown in Figure 11.Table 4
Miss distances in the two scenarios of the terminal maximum overload maneuver or no maneuver.
Penetration angle
Miss distance (no maneuver)
Miss distance (maneuver)
10
0.15
0.26
15
0.30
0.78
20
2.89
74.36
25
253.25
757.9
30
1705.34
1827.27Figure 9
The attack and defense trajectories in the scenario of no terminal maneuver.Figure 10
The attack and defense trajectories in the scenario of the maximum overload terminal maneuver.Figure 11
The attack and defense trajectories in the scenarios of no maneuver and maximum overload maneuver (horizonal plane).Table4 shows that if the penetration angle in the countermeasure is larger than 20 degrees, no matter whether or not the vehicle maneuvers in its terminal countermeasure, it can successfully penetrate. Otherwise, the interception missile will create a head-on strike situation; then even the high-overload terminal maneuver of the hypersonic vehicle cannot evade the attack by an interception missile. The penetration angle successfully transforms the complicated countermeasure process of the hypersonic vehicle’s evasion into a simple penetration angle constraint. The penetration can be accomplished without any energy waste.


## 4.1. Settings of Simulation Parameters

The basic parameters and the relevant simulation parameter settings in the attack and defense countermeasure simulation model of the hypersonic vehicle are given in Tables1 and 2. The boundary values in the control constraints are set as αTmin=−4°, αTmax=6°, σ̇max=5°, and msmax=1; the boundary value in the penetration angle constraints is set as ηmax=20°.Table 1
The vehicle’s basic parameters [8].
Parameters
CL0
CLα (rad−1)
CD0
CDα2 (rad−2)
CTms (s·kg−1)
Parameter values
-0.013
1.833
0.015
4.596
0.162Table 2
The initial values of conditions of the attack side and the defense side.
Simulation parameters
Parameter values
xe0,ye0,ze0 (km)
100,25,10
Ve0 (m·s−1)
1500
θe0,ψe0 (rad)
0,π
xi0,yi0,zi0 (km)
0,0,0
Vi (m·s−1)
1200
θi0,ψi0 (rad)
(0.2438, -0.0997)
N1,N2,N3
5,0.1,0.1i/N
R
diag5,0.1,1
rf (km)
20Simulation conditions are set as follows: in the first iteration, control guesses are set tou=2°,0,0.2T, the number of discrete points is set as N=200, and the thrust domain constraint values in the formula are set as follows:
(25)δ=200020π18020π180100020π18020π1801000120π18020π18020π180T.The error threshold values in the simulation iteration termination condition are set as follows:(26)ε=1000.5π1800.5π1801000.5π1800.5π18050010.5π1800.5π1800.5π180T.The above simulation conditions are used to perform numerical simulation with the YALMIP modeling environment of MATLAB; the primal-dual interior point algorithm in the MOSEK [25] software package is used to solve the SOCP problem. After five times of sequential iteration, the trajectory converges. The time required by each iteration is around 0.4 s; the resolving time required by the interior point algorithm is around 0.2 s. The penetration angle and the hypersonic vehicle’s angle of attack calculated by the first four iterations are shown in Figures 3 and 4. The figures show that the sequential convex optimization algorithm converges quickly; the angle-of-attack value solved in the second iteration is already close to the ultimate optimal solution.Figure 3
Penetration angle in the first four successive iterations of the SOCP-based algorithm.Figure 4
Angle-of-attack profile in the first four successive iterations of the SOCP-based algorithm.In addition, in order to validate the effectiveness of the second-order cone programming (SOCP) method proposed in the paper, numerical results are compared with those obtained by the Gaussian pseudospectral method (GPM). The GPM solves the problem with the sequential quadratic programming (SQP), and we program it by the GPOPS5.1 software package [26].The GPM has already been proved to obtain its globally optimal solution [27].The trajectory profile by the SOCP-based method and the GPM (the mesh tolerance is set to 10-3) are shown in Figures 5–7, while the solving times required and optimal performance indexes are given in Table 3. Figure 8 presents the required calculating time of the two methods using different random control initial guesses within control constraints; the horizontal axis only represents the number of initial guesses without any specific physical meaning.Figure 5
The penetrator’s attack angle profiles, using the two methods.Figure 6
The penetrator’s banking angle profiles, using the two methods.Figure 7
The penetrator’s fuel consumption rate profiles, using the two methods.Table 3
Comparison of the time consumed with performance indexes, using the SOCP and GPM methods.
Method
Discrete point
Iteration
Time (s)
Performance index
GPM
40
—
7.76
751
SOCP
200
5
1.42
770Figure 8
Calculating time for different initial guesses.Based on the figures, the SOCP and GPM solutions show very similar trends; the control profiles in Figures5–7 match well in general. The relatively small differences between the two methods are attributable to the different discretization strategies used in the two solutions, which also occur in reference [28]. On the other hand, the performance of the GPM is only slightly optimal than that of SOCP in terms of the values of 751 and 770, respectively. Therefore, we can at least reckon that the SOCP method is near optimal. Despite the compromise of the optimum objective, the SOCP method performs much better in terms of computational complexity based on Table 3 and Figures 4 and 8. The SOCP-based method in this paper takes an average of 1.3 seconds to generate trajectories, and the peak calculating time for different initial guesses is only 2.45 seconds, while the GPM method takes an average of 7.98 seconds and the longest time can reach 13.1 seconds. To tackle penetration problems, a method with competitive convergence speed and robust solution process is what we always seek for, though it is not rigorously optimal.Further simulation is conducted to verify whether the penetration angle could function as an index for evaluating the penetration effects when the distance between the missile and the hypersonic vehicle is fixed. In other words, the hypersonic vehicle can successfully escape once it meets certain penetration angle constraints in the attack and defense countermeasure.Therefore, we design different penetration angle constraintsη=10°,15°,20°,25°,30° and conduct the simulation in which the continuing terminal-course relative motion is considered. In the simulation, when the distance between the attack side and the defense side equals to rf, the countermeasure enters the terminal phase where the two sides continue to approach each other till the vehicle is intercepted or successfully escapes. During this period, the guidance law of the interception missile is the pure proportional guidance law with the proportional constant of 3, while the vehicle adopts two kinds of maneuver in order to illustrate the impact of angle constraints. It does not make any maneuver and continues flying at the current speed and direction or maneuvers with maximum overload in the opposite direction of the attacked missile. The maximum overload, about 2.8g, is obtained by the maximum available attack angle and maximum thrust for a hypersonic vehicle. The formulas for calculating its overload and maneuver direction are as follows:
(27)n¯cmax=ncmaxI¯n,ncmax=TmaxsinαTmax+LTαTmaxm,I¯n=V¯T×Los¯×V¯TV¯T×Los¯×V¯T,where I¯n represents the maneuver direction; V¯T and Los¯ represent the vehicle’s velocity vector and line-of-sight vector, respectively. With various factors being considered, the time constant of an interception missile is set as 0.05 s and the miss distance is designed as 1 m.Under different penetration angle constraints and two scenarios of the hypersonic vehicle’s continuing motion method, the miss distances of the interception missile in the two scenarios are shown in Table4. The whole trajectory profiles for the attack side and the defense side are given in Figures 9 and 10, while the figure for comparing the terminal trajectory (20°penetration angle) in the horizonal plane between different maneuver methods is also shown in Figure 11.Table 4
Miss distances in the two scenarios of the terminal maximum overload maneuver or no maneuver.
Penetration angle
Miss distance (no maneuver)
Miss distance (maneuver)
10
0.15
0.26
15
0.30
0.78
20
2.89
74.36
25
253.25
757.9
30
1705.34
1827.27Figure 9
The attack and defense trajectories in the scenario of no terminal maneuver.Figure 10
The attack and defense trajectories in the scenario of the maximum overload terminal maneuver.Figure 11
The attack and defense trajectories in the scenarios of no maneuver and maximum overload maneuver (horizonal plane).Table4 shows that if the penetration angle in the countermeasure is larger than 20 degrees, no matter whether or not the vehicle maneuvers in its terminal countermeasure, it can successfully penetrate. Otherwise, the interception missile will create a head-on strike situation; then even the high-overload terminal maneuver of the hypersonic vehicle cannot evade the attack by an interception missile. The penetration angle successfully transforms the complicated countermeasure process of the hypersonic vehicle’s evasion into a simple penetration angle constraint. The penetration can be accomplished without any energy waste.


## 5. Conclusions

Confronted with the serious interception threat, hypersonic vehicles need to be empowered with the ability to generate a penetration trajectory, cope with air defense missiles, and attack targets simultaneously. In this paper, a rapid penetration trajectory programming method is proposed to guide the hypersonic vehicle to evade against a head-on interception missile with the minimum energy consumption. It has been demonstrated that the method raised in this paper is effective and has competitive computational performances including a quick convergence rate, rapid calculation speed, and robustness of initial guesses. Moreover, the penetration angle summarized in this paper proves to be a valid index of the penetration effect, facilitating the application of optimization methods to hypersonic penetration problems. As a whole, these research results contribute to the future online application of hypersonic penetration programming.



---

*Source: 1490342-2019-07-09.xml*