# Near-Optimal Guidance with Impact Angle and Velocity Constraints Using Sequential Convex Programming

Authors: Pei Pei; Jiang Wang
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

This paper proposes a near-optimal air-to-ground missile guidance law with impact angle and impact velocity constraints based on sequential convex programming. A realistic aerodynamic model is introduced into the problem formulation, such that traditional optimization theory cannot obtain an analytical solution to the optimization problem under state constraints. The original problem is considered as an optimization problem, and the angle of attack is replaced with the angle of attack rate as a new control variable to reconstruct the problem and simplify the solving process. Next, the independent variable is changed in the differential equations to linearize and discretize the problem such that the reconstructed problem can be solved using sequential convex programming. The results obtained by numerical simulations confirmed that the proposed algorithm is valid and faster than the general-purpose nonlinear optimal control problem solver. Finally, it was verified that different impact angles and impact velocities were achieved.


---

## Body

## 1. Introduction

In modern warfare, the objectives of the guidance law are not only limited to zero miss distance interception but also include target interception at a certain impact angle and impact velocity. The impact angle constraint in a terminal interception engagement is critical for a homing missile attacking modern warships, tanks, and ballistic missiles. Additionally, it can increase the effectiveness and lethality of the missile’s warhead and realize an escape from the limited defense zone of the target [1, 2]. The impact velocity constraint at a terminal homing guidance is crucial to the visibility of seekers and for precise and successful mission execution [3, 4]. The guidance law with impact angle and impact velocity constraints has the advantage of effectively destroying a target.The guidance law considering the impact angle and impact velocity constraints has been widely investigated in previous decades. Li et al. [5] used graph theory to derive the impact angle and time constraint guidance law. Sliding mode theory [6–9] and proportional navigation [10–12] have been introduced to achieve the desired impact angle. Optimal control theory has also been implemented to solve optimal guidance law problems under an impact angle constraint [13–18], but a limit has not been set on the terminal velocity owing to the complicated form of the aerodynamic coefficients. Kim et al. [19] converted the guidance law to a polynomial form and used polynomial coefficients to constrain the impact time and impact angle. Later, Tahk et al. [20] extended the polynomial form of the guidance law to control the impact velocity with another polynomial coefficient, which is clearly not an optimal solution. Moreover, the drag coefficient in the guidance law design was considered as a constant, which is unrealistic and apparently invalid for variable aerodynamic coefficients.Owing to the time-varying and complicated form of the aerodynamic terms in the velocity dynamic equation, it is hard to obtain an analytical form of the optimal guidance law for velocity control. Hence, computational guidance [21] was introduced to handle the optimal real-time control problems or trajectory optimization problems with velocity constraints. Cheng et al. [22] used the Newton–Kantorovich/pseudospectral approach to solve the convex model of ascent trajectory optimization problems, while the drag coefficient is easily considered as constant. Convex optimization [23–29] is increasingly becoming attractive in computational guidance owing to its efficiency and reliability. Liu et al. [30] considered a more realistic form of lift and drag, then used the drag polar curve to turn the drag and lift into a single variable to control the impact velocity, and finally used second-order convex programming (SOCP) to solve the control variable. However, the drag polar curve was not exactly accurate, which caused the real flight trajectory and impact velocity to deviate from the simulations. Moreover, to maintain the equilibrium of the original problem and transformed problem, an additional term must be added to the objective function. The independent variable altitude chosen by Liu makes it impossible to handle a wavy trajectory, which is contrary to the monotonically independent variable.In this study, the realistic form of aerodynamic coefficients is introduced to precisely control the impact angle and impact velocity. The drag and lift coefficients are set as functions of the angle of attack and flight Mach number, which are obtained by curve fitting based on the flight parameters. To overcome the inherent shortcoming of traditional optimization theory, whereby the analytical solution to the realistic velocity control problem cannot be obtained, the sequential convex programming algorithm is proposed to handle the time-varying nonlinear optimization problem. A new control variable and an independent variable are introduced to reformulate the original problem as an optimization problem. Additionally, convexification and discretization render the problem solvable using sequential convex programming. The main novelty of this study is the derivation of a realistic model for velocity control. Moreover, sequential convex programming is proposed to numerically solve the optimization problem in a nonanalytical manner. The numerical solutions confirmed the validity and effectiveness of the proposed approach.The rest of this paper is organized as follows: in Section2, the near-optimal problem is formulated and transformed into a linearized form. Then, in Section 3, the linearized form of the near-optimal problem is discretized and solved using sequential convex programming. The numerical simulations are presented in Section 4. Finally, the conclusions drawn from this study are presented in Section 5.


## 2. Problem Formulation

In this section, we formulate an optimization problem for an air-to-ground aerodynamically controlled missile impacting a target at a desired impact angle and impact velocity with minimal control effort and according to the state inequality constraints on the aerodynamic limitation. Unlike other simplified dynamic and kinematic equations, the drag and lift terms are introduced into the mathematical model.


### 2.1. Control Effort Optimization Interception Problem

First, let us consider the planar interception engagement of an aerodynamically controlled missile. The aerodynamic and kinematic equations are expressed as follows:(1)x˙=Vcosθ,y˙=Vsinθ,V˙=−Dm−gsinθ,θ˙=LmV−gcosθV,where x, y, V, and θ denote the missile longitudinal range, altitude, velocity, and flight path angle, respectively; g is the gravitational acceleration; D is the drag force; and L is the lift force, which can be calculated as follows:(2)D=12CDα,MρV2Sref,L=12CLα,MρV2Sref,where ρ=ρ0e−y/hs is the atmospheric density; Sref is the reference area of the missile; and CD and CL are the drag and lift coefficients, respectively, both of which depend on the angle of attack α and Mach number M. In some optimization problems [30], the drag polar curve [31] has been used to simplify the problem, but this approach is not sufficiently precise for guidance law design. For the control effort optimization problem considered in this study, the aerodynamic coefficients are approximated by fitting the data to all flight envelopes, as follows [32]:(3)CD=0.4136−0.008152α+0.0112M+0.00422α2+0.002126Mα,CL=−0.2128+0.2284α+0.07589M+0.04997Mα.In summary, the aerodynamic and kinematic equations can be rearranged as follows:(4)x˙=Vcosθ,y˙=Vsinθ,V˙=−CDα,Mρ0−y/hsV2Sref2m−gsinθ,θ˙=CLα,Mρ0−y/hsVSref2m−gcosθV.In this study, the Mach numberM can be calculated through the missile velocity. Consequently, the angle of attack α is assumed to be the only control variable for the control optimization problem.The physical constraints that must be satisfied in the problem are as follows:(1)
Engagement kinematics: it is expressed by equation (4) with initial conditions xt0, yt0, Vt0, θt0, and α0, which are the current values of these variables from the navigation system:(5)θt0=θ0∗,Vt0=V0∗,xt0=x0,yt0=y0.(2)
Control constraint: owing to the physical constraint of the missile, the angle of attack during the engagement is limited as follows:(6)αt≤αmax,whereαmax is the maximum angle of attack.(3)
Terminal constraints: to impact the target from a specified directionθf∗ and specified velocity Vf∗, we require ttf=0, xtf=xf, and ytf=yf, where the final impact time tf is free and xf and yf are the target coordinates. Therefore, the terminal constraints are expressed as follows:(7)θtf=θf∗,Vtf=Vf∗,xtf=xf,ytf=yf.In addition to satisfying all of these constraints, the optimization objective also consists of minimizing the control energy; then, the following performance index should be minimized, as follows:(8)J=∫t0tfαt2dt.Now, the energy optimization problem can be rewritten as follows:Problem 1.
Minimize equation (8) subject to equations (4)–(7).


### 2.2. Choice of New Control and New Independent Variable

The dynamic and kinematic equations expressed by equation (4) are highly nonlinear in terms of both the state variables and control variable. This causes high-frequency jitters when applying successive linearization to the kinematic and dynamic equations. The control variable α acting on the aerodynamic coefficients is highly nonlinear and coupled with the state variable such that it becomes complicated to obtain a numerical optimization solution to Problem 1. To simplify Problem 1 and eliminate the jitter phenomena, we define the angle of attack rate α˙ as a new control input and consider the missile longitudinal range x as a new independent variable. Note that x˙ is positive-defined, owing to the limits on the flight path angle θ≤π/2, which means that x monotonically increases during the engagement. The following state equation is added to the original equations:(9)α˙=u.It is known that the angle of attack rate cannot be arbitrarily large; thus, it is limited by the physical characteristic of the missile, as follows:(10)u≤α˙max,where α˙max is the maximum allowable angle of attack rate.By adding the angle of attack rate dynamics expressed in equation (9) to equation (4) and changing the independent variable, the augmented equations of the dynamic and kinematic equations can be reformed as follows:(11)x′=fx+Bxu,where x is the state vector given as x=y;V;θ;α and the quotation mark denotes that the differentiations are now with respect to x. The column vectors fx∈ℝ4 and Bx∈ℝ4 are expressed as follows:(12)fx=tanθ−CDρ0e−y/hsVSref2mcosθ−gtanθVCLρ0e−y/hsSref2mcosθ−gV20,Bx=0001Vcosθ.Next, withx as the independent variable, the state, initial, and terminal constraints in equations (5)–(7) can be written as follows:(13)αx≤αmax,(14)θx0=θ0∗,Vx0=V0∗,yx0=y0,(15)θxf=θf∗,Vxf=Vf∗,yxf=yf.Finally, withx as the independent variable, the performance index can be written as follows (using dt=dx/x˙):(16)J=∫x0xfα2Vcosθdx.Remark 1.
The original control variable angle of attack acts on the aerodynamic coefficients, is of high order, and is coupled with the states. This form is complicated for optimal control theory and numerical optimization. With the angle of attack rate as the new control variable, the control variable is now a first-order term, which benefits the numerical optimization.


### 2.3. Reformulated Optimization-Control Problem

With the original dynamic and kinematic equations in equation (4) replaced by the augmented dynamics and new independent variable, Problem 1 can be reformulated into a new optimization-control problem, as follows:Problem 2.
Minimize equation (16) subject to equations (10)–(15).Lemma 1.
If there exists an optimization solution to Problem 2 given byx∗,u∗=y∗,V∗,θ∗,α∗,u∗, then, x∗,y∗,V∗,θ∗,α∗ is a feasible solution to Problem 1.Proof.
Problem 2 is given by introducing a new control variable and new independent variable to the original system in equation (4). The original control constraint on α becomes a boundary condition on the new state variable. The angle of attack rate is introduced as a new control variable into equation (9) and as a constraint on the angle of attack rate in equation (10). Thus, the set of feasible controls for Problem 2 are strictly constrained in the set of feasible controls for Problem 1. Obviously, Problem 2 is merely a relaxation of Problem 1. Additionally, the optimization of Problem 2 satisfies all constraints of Problem 1 and clearly defines a feasible solution to Problem 1. However, a feasible solution to Problem 1 does not necessarily define a solution to the original Problem 1, and the equivalence of the dynamics and kinematics with and without α˙=u is demonstrated by the numerical simulations.Remark 2.
Lemma 1 considers a relationship between the solution of the original Problem 1 and the relaxed Problem 2. The nonlinear and coupled control in equation (4) is avoided by introducing a new control variable and new independent variable. Consequently, a feasible solution to Problem 1 can be obtained by solving the problem. The numerical simulations validate the effectiveness of the relaxed transformation.


## 2.1. Control Effort Optimization Interception Problem

First, let us consider the planar interception engagement of an aerodynamically controlled missile. The aerodynamic and kinematic equations are expressed as follows:(1)x˙=Vcosθ,y˙=Vsinθ,V˙=−Dm−gsinθ,θ˙=LmV−gcosθV,where x, y, V, and θ denote the missile longitudinal range, altitude, velocity, and flight path angle, respectively; g is the gravitational acceleration; D is the drag force; and L is the lift force, which can be calculated as follows:(2)D=12CDα,MρV2Sref,L=12CLα,MρV2Sref,where ρ=ρ0e−y/hs is the atmospheric density; Sref is the reference area of the missile; and CD and CL are the drag and lift coefficients, respectively, both of which depend on the angle of attack α and Mach number M. In some optimization problems [30], the drag polar curve [31] has been used to simplify the problem, but this approach is not sufficiently precise for guidance law design. For the control effort optimization problem considered in this study, the aerodynamic coefficients are approximated by fitting the data to all flight envelopes, as follows [32]:(3)CD=0.4136−0.008152α+0.0112M+0.00422α2+0.002126Mα,CL=−0.2128+0.2284α+0.07589M+0.04997Mα.In summary, the aerodynamic and kinematic equations can be rearranged as follows:(4)x˙=Vcosθ,y˙=Vsinθ,V˙=−CDα,Mρ0−y/hsV2Sref2m−gsinθ,θ˙=CLα,Mρ0−y/hsVSref2m−gcosθV.In this study, the Mach numberM can be calculated through the missile velocity. Consequently, the angle of attack α is assumed to be the only control variable for the control optimization problem.The physical constraints that must be satisfied in the problem are as follows:(1)
Engagement kinematics: it is expressed by equation (4) with initial conditions xt0, yt0, Vt0, θt0, and α0, which are the current values of these variables from the navigation system:(5)θt0=θ0∗,Vt0=V0∗,xt0=x0,yt0=y0.(2)
Control constraint: owing to the physical constraint of the missile, the angle of attack during the engagement is limited as follows:(6)αt≤αmax,whereαmax is the maximum angle of attack.(3)
Terminal constraints: to impact the target from a specified directionθf∗ and specified velocity Vf∗, we require ttf=0, xtf=xf, and ytf=yf, where the final impact time tf is free and xf and yf are the target coordinates. Therefore, the terminal constraints are expressed as follows:(7)θtf=θf∗,Vtf=Vf∗,xtf=xf,ytf=yf.In addition to satisfying all of these constraints, the optimization objective also consists of minimizing the control energy; then, the following performance index should be minimized, as follows:(8)J=∫t0tfαt2dt.Now, the energy optimization problem can be rewritten as follows:Problem 1.
Minimize equation (8) subject to equations (4)–(7).


## 2.2. Choice of New Control and New Independent Variable

The dynamic and kinematic equations expressed by equation (4) are highly nonlinear in terms of both the state variables and control variable. This causes high-frequency jitters when applying successive linearization to the kinematic and dynamic equations. The control variable α acting on the aerodynamic coefficients is highly nonlinear and coupled with the state variable such that it becomes complicated to obtain a numerical optimization solution to Problem 1. To simplify Problem 1 and eliminate the jitter phenomena, we define the angle of attack rate α˙ as a new control input and consider the missile longitudinal range x as a new independent variable. Note that x˙ is positive-defined, owing to the limits on the flight path angle θ≤π/2, which means that x monotonically increases during the engagement. The following state equation is added to the original equations:(9)α˙=u.It is known that the angle of attack rate cannot be arbitrarily large; thus, it is limited by the physical characteristic of the missile, as follows:(10)u≤α˙max,where α˙max is the maximum allowable angle of attack rate.By adding the angle of attack rate dynamics expressed in equation (9) to equation (4) and changing the independent variable, the augmented equations of the dynamic and kinematic equations can be reformed as follows:(11)x′=fx+Bxu,where x is the state vector given as x=y;V;θ;α and the quotation mark denotes that the differentiations are now with respect to x. The column vectors fx∈ℝ4 and Bx∈ℝ4 are expressed as follows:(12)fx=tanθ−CDρ0e−y/hsVSref2mcosθ−gtanθVCLρ0e−y/hsSref2mcosθ−gV20,Bx=0001Vcosθ.Next, withx as the independent variable, the state, initial, and terminal constraints in equations (5)–(7) can be written as follows:(13)αx≤αmax,(14)θx0=θ0∗,Vx0=V0∗,yx0=y0,(15)θxf=θf∗,Vxf=Vf∗,yxf=yf.Finally, withx as the independent variable, the performance index can be written as follows (using dt=dx/x˙):(16)J=∫x0xfα2Vcosθdx.Remark 1.
The original control variable angle of attack acts on the aerodynamic coefficients, is of high order, and is coupled with the states. This form is complicated for optimal control theory and numerical optimization. With the angle of attack rate as the new control variable, the control variable is now a first-order term, which benefits the numerical optimization.


## 2.3. Reformulated Optimization-Control Problem

With the original dynamic and kinematic equations in equation (4) replaced by the augmented dynamics and new independent variable, Problem 1 can be reformulated into a new optimization-control problem, as follows:Problem 2.
Minimize equation (16) subject to equations (10)–(15).Lemma 1.
If there exists an optimization solution to Problem 2 given byx∗,u∗=y∗,V∗,θ∗,α∗,u∗, then, x∗,y∗,V∗,θ∗,α∗ is a feasible solution to Problem 1.Proof.
Problem 2 is given by introducing a new control variable and new independent variable to the original system in equation (4). The original control constraint on α becomes a boundary condition on the new state variable. The angle of attack rate is introduced as a new control variable into equation (9) and as a constraint on the angle of attack rate in equation (10). Thus, the set of feasible controls for Problem 2 are strictly constrained in the set of feasible controls for Problem 1. Obviously, Problem 2 is merely a relaxation of Problem 1. Additionally, the optimization of Problem 2 satisfies all constraints of Problem 1 and clearly defines a feasible solution to Problem 1. However, a feasible solution to Problem 1 does not necessarily define a solution to the original Problem 1, and the equivalence of the dynamics and kinematics with and without α˙=u is demonstrated by the numerical simulations.Remark 2.
Lemma 1 considers a relationship between the solution of the original Problem 1 and the relaxed Problem 2. The nonlinear and coupled control in equation (4) is avoided by introducing a new control variable and new independent variable. Consequently, a feasible solution to Problem 1 can be obtained by solving the problem. The numerical simulations validate the effectiveness of the relaxed transformation.


## 3. Convexification and SCP Method

Problem 2 is still a highly nonlinear optimization problem with nonlinear differential state equation constraints and an objective function, which adds complexity when applying a convex-optimization solver to obtain the solution. However, except for the differential state equations and objective function, all other constraints on the state variables and control variables are linear and convex. Therefore, the nonlinear part must be converted to tractable formulations for convex programming. A small-disturbance-based linearization method is used to partially linearize the terms that are unrelated to the control [30]. Subsequently, a sequence of convex-optimization problem, specifically to a sequence of SOCP problem, is introduced to approximate Problem 2. The numerical simulations in Section 4 demonstrate that the accuracy of the successive linearization approximation of Problem 2 is acceptable.


### 3.1. Successive Linear Approximation

First, letxk,uk be the kth successive solution, which has already been obtained. For the fast convergence of the successive solution, we linearize the term fx at xk in equation (11) but approximate Bx with Bxk to avoid the presence of uk in the resulting linearized system. Then, we convert equation (11) to the following partially linearized system:(17)x′=Axkx+Bxku+Cxk,where xk=ykVkθkαkT and(18)Axk=dfdxxk=00a130a21a22a23a24a31a32a33a340000,Bxk=000b4,cxk=fxk−Axkxk=c1c2c30.The coefficients are defined as follows, and subscriptk of the state variables is ignored for simplicity:(19)a13=1cos2θ,a21=CDρ0e−y/hsVSref2mcosθhs,a22=−CDρ0e−y/hsSref2mcosθ−CDVρ0e−y/hsVSref2mcosθ+gtanθV2,a23=−CDρ0e−y/hsVSrefsinθ2mcos2θ−gVcos2θ,a24=−CDαρ0e−y/hsVSref2mcosθ,a31=−CLρ0e−y/hsSref2mcosθhs,a32=2gV3+CLVρ0e−y/hsSref2mcosθ,a33=CLρ0e−y/hsSrefsinθ2mcos2θ,a34=CLαρ0e−y/hsSref2mcosθ,b4=1Vcosθ,c1=tanθ−θcos2θ,c2=−2gtanθV−yCDρ0e−y/hsVSref2mcosθhs+CDρ0e−y/hsVSrefsinθθ2mcos2θ+gθVcos2θ+CDVρ0e−y/hsV2Sref2mcosθ+CDαρ0e−y/hsVSrefα2mcosθ,c3=CLρ0e−y/hsSref2mcosθ−3gV2+yCLρ0e−y/hsSref2mcosθhs−CLρ0e−y/hsSrefsinθθ2mcos2θ−CLVρ0e−y/hsVSref2mcosθ−CLαρ0e−y/hsSrefα2mcosθ,CDV=dCDdV,CDα=dCDdα,CLV=dCLdV,CLα=dCLdα.Additionally, a trust region must be added to maintain the validity of the proposed linearization, as follows:(20)x−xk≤δ,where δ∈ℝ4 is a constant vector, and the inequality is applied componentwise.Next, minimizing the preceding performance index expressed in equation (16) is equivalent to minimizing the following objective function:(21)J=∫x0xfηdx.With an additional constraint expressed as follows:(22)α2Vcosθ≤η,where η is a slack variable. Obviously, the performance index in equation (16) is a linear integrand. The inequality constraint in equation (22) can be approximated using a “lagging” technique, as follows:(23)α2≤wkη,where wk=Vkcosθk.In summary, the optimization-control problem approximating Problem 2 is presented as follows.Problem 3.
Minimize equation (21) subject to equations (13), (15), (17), (20), and (23).
The derived Problem 3 can be very efficiently solved using a convex-optimization algorithm because the discretization problem will yield a standard SOCP problem, which has been mathematically defined in the literature [33], as follows:
Minimize(24)fTx.
According to(25)Aix+bi2≤ciTx+di,i=1,…,m,(26)Fx=g,where x∈ℝn is the optimization variable, and the problem parameters include f∈ℝn, F∈ℝp×n, g∈ℝp, Ai∈ℝni×n, bi∈ℝn, and di∈ℝn. Note that, when Ai=0, the second-order cone constraints are reduced to general linear constraints.
For Problem 3, in the fixed intervalx0,xf, the performance index in equation (21) can be discretized into a linear form of equation (24). The inequality constraints in equations (13), (20), and (23) are converted to a sequence of second-order cone constraints. The dynamic and kinematic equations expressed by equation (17) and the terminal constraints in equation (15) are transformed to linear constraints.Remark 2.
The trapezoidal rule is implemented to discretize the performance index in equation (21) and the differential equations in equation (17) at the n+1 uniformly distributed discretized point in x0,xf; then, the step length is e=xf−x0/n. Therefore, the performance index can be discretized to a linear form, as follows:(27)J=∫x0xfηdx=∑i=2ne2ηi−1+ηi.
Similarly, the differential equations in equation (17) can be discretized and rearranged as follows:(28)Hi−1xi−1−Hixi+Gi−1ui−1+Giui+Ci−1+Ci=0,where the coefficients are expressed as follows:(29)Hi−1=I+e2Ai−1xk;Hi=I−e2Aixk;Gi−1=e2Bi−1xk;Gi=e2Bixk;Ci−1=−e2Ci−1xk;Ci=−e2Cixk.
In summary, the nonlinear integrand performance index and differential equations can be transformed into the convex constraints in equation (26).


### 3.2. Sequential Convex Programming

We already known that Problem 3 is an approximation to Problem 2. For the sake of obtaining an accurate solution to Problem 2, we propose the following sequential convex programming by solving a sequence of convex programming problems defined by Problem 3 to approximate the exact solution to Problem 2. The detail procedures are as follows:(1)
Setk=0; initialize the states yx0=y0, Vx0=V0, θx0=θ0, αx0=0, yxf=yf, Vxf=Vf∗, θxf=θf∗, and αxf=0. Note that the initial states in the convex-optimization algorithm do not need to satisfy the constraints; therefore, we can generate x0 by linearly interpolating between the initial condition and the final condition.(2)
Fork>1, compute the xk−1-dependent parameters in equations (17) and (28) using xk−1. Then, Problem 3 can be transformed into a SOCP problem and solved with the initial states xk−1 and preceding calculated parameters. Therefore, we can obtain the solution xk,uk.(3)
Check the following convergence condition:(30)supx0<x<xfxk−xk−1≤ξ,whereξ∈ℝ4 is the prescribed tolerance value for convergence. If the preceding convergence condition is satisfied, go to Step (4); otherwise, set k=k+1 and go to Step (2).(4)
The solution to Problem 2 is determined asx∗=xk and u∗=uk.Remark 4.
In equations (17) and (23), and during the discretization process, the parameters xk and wk depend on the solution but are unknown before the problem is solved. However, we can approximate their values based on the solution obtained in the previous solution procedure. The first solution is given as the initial condition in Step (1) and does not need to satisfy all of the constraints in Problem 3, but rather the different initial states will influence the iteration times to converge in the sequential convex programming.Remark 5.
The theoretical proof for the convergence of the sequential convex programming is very challenging owing to the presence of various constraints, nonlinear dynamics, and kinematics. Lu and Liu [34] and Wang and Grant [35, 36] provided several examples of sequential convex programming convergence. The numerical simulations presented in the next part demonstrate the convergence and validity of the proposed algorithm.


## 3.1. Successive Linear Approximation

First, letxk,uk be the kth successive solution, which has already been obtained. For the fast convergence of the successive solution, we linearize the term fx at xk in equation (11) but approximate Bx with Bxk to avoid the presence of uk in the resulting linearized system. Then, we convert equation (11) to the following partially linearized system:(17)x′=Axkx+Bxku+Cxk,where xk=ykVkθkαkT and(18)Axk=dfdxxk=00a130a21a22a23a24a31a32a33a340000,Bxk=000b4,cxk=fxk−Axkxk=c1c2c30.The coefficients are defined as follows, and subscriptk of the state variables is ignored for simplicity:(19)a13=1cos2θ,a21=CDρ0e−y/hsVSref2mcosθhs,a22=−CDρ0e−y/hsSref2mcosθ−CDVρ0e−y/hsVSref2mcosθ+gtanθV2,a23=−CDρ0e−y/hsVSrefsinθ2mcos2θ−gVcos2θ,a24=−CDαρ0e−y/hsVSref2mcosθ,a31=−CLρ0e−y/hsSref2mcosθhs,a32=2gV3+CLVρ0e−y/hsSref2mcosθ,a33=CLρ0e−y/hsSrefsinθ2mcos2θ,a34=CLαρ0e−y/hsSref2mcosθ,b4=1Vcosθ,c1=tanθ−θcos2θ,c2=−2gtanθV−yCDρ0e−y/hsVSref2mcosθhs+CDρ0e−y/hsVSrefsinθθ2mcos2θ+gθVcos2θ+CDVρ0e−y/hsV2Sref2mcosθ+CDαρ0e−y/hsVSrefα2mcosθ,c3=CLρ0e−y/hsSref2mcosθ−3gV2+yCLρ0e−y/hsSref2mcosθhs−CLρ0e−y/hsSrefsinθθ2mcos2θ−CLVρ0e−y/hsVSref2mcosθ−CLαρ0e−y/hsSrefα2mcosθ,CDV=dCDdV,CDα=dCDdα,CLV=dCLdV,CLα=dCLdα.Additionally, a trust region must be added to maintain the validity of the proposed linearization, as follows:(20)x−xk≤δ,where δ∈ℝ4 is a constant vector, and the inequality is applied componentwise.Next, minimizing the preceding performance index expressed in equation (16) is equivalent to minimizing the following objective function:(21)J=∫x0xfηdx.With an additional constraint expressed as follows:(22)α2Vcosθ≤η,where η is a slack variable. Obviously, the performance index in equation (16) is a linear integrand. The inequality constraint in equation (22) can be approximated using a “lagging” technique, as follows:(23)α2≤wkη,where wk=Vkcosθk.In summary, the optimization-control problem approximating Problem 2 is presented as follows.Problem 3.
Minimize equation (21) subject to equations (13), (15), (17), (20), and (23).
The derived Problem 3 can be very efficiently solved using a convex-optimization algorithm because the discretization problem will yield a standard SOCP problem, which has been mathematically defined in the literature [33], as follows:
Minimize(24)fTx.
According to(25)Aix+bi2≤ciTx+di,i=1,…,m,(26)Fx=g,where x∈ℝn is the optimization variable, and the problem parameters include f∈ℝn, F∈ℝp×n, g∈ℝp, Ai∈ℝni×n, bi∈ℝn, and di∈ℝn. Note that, when Ai=0, the second-order cone constraints are reduced to general linear constraints.
For Problem 3, in the fixed intervalx0,xf, the performance index in equation (21) can be discretized into a linear form of equation (24). The inequality constraints in equations (13), (20), and (23) are converted to a sequence of second-order cone constraints. The dynamic and kinematic equations expressed by equation (17) and the terminal constraints in equation (15) are transformed to linear constraints.Remark 2.
The trapezoidal rule is implemented to discretize the performance index in equation (21) and the differential equations in equation (17) at the n+1 uniformly distributed discretized point in x0,xf; then, the step length is e=xf−x0/n. Therefore, the performance index can be discretized to a linear form, as follows:(27)J=∫x0xfηdx=∑i=2ne2ηi−1+ηi.
Similarly, the differential equations in equation (17) can be discretized and rearranged as follows:(28)Hi−1xi−1−Hixi+Gi−1ui−1+Giui+Ci−1+Ci=0,where the coefficients are expressed as follows:(29)Hi−1=I+e2Ai−1xk;Hi=I−e2Aixk;Gi−1=e2Bi−1xk;Gi=e2Bixk;Ci−1=−e2Ci−1xk;Ci=−e2Cixk.
In summary, the nonlinear integrand performance index and differential equations can be transformed into the convex constraints in equation (26).


## 3.2. Sequential Convex Programming

We already known that Problem 3 is an approximation to Problem 2. For the sake of obtaining an accurate solution to Problem 2, we propose the following sequential convex programming by solving a sequence of convex programming problems defined by Problem 3 to approximate the exact solution to Problem 2. The detail procedures are as follows:(1)
Setk=0; initialize the states yx0=y0, Vx0=V0, θx0=θ0, αx0=0, yxf=yf, Vxf=Vf∗, θxf=θf∗, and αxf=0. Note that the initial states in the convex-optimization algorithm do not need to satisfy the constraints; therefore, we can generate x0 by linearly interpolating between the initial condition and the final condition.(2)
Fork>1, compute the xk−1-dependent parameters in equations (17) and (28) using xk−1. Then, Problem 3 can be transformed into a SOCP problem and solved with the initial states xk−1 and preceding calculated parameters. Therefore, we can obtain the solution xk,uk.(3)
Check the following convergence condition:(30)supx0<x<xfxk−xk−1≤ξ,whereξ∈ℝ4 is the prescribed tolerance value for convergence. If the preceding convergence condition is satisfied, go to Step (4); otherwise, set k=k+1 and go to Step (2).(4)
The solution to Problem 2 is determined asx∗=xk and u∗=uk.Remark 4.
In equations (17) and (23), and during the discretization process, the parameters xk and wk depend on the solution but are unknown before the problem is solved. However, we can approximate their values based on the solution obtained in the previous solution procedure. The first solution is given as the initial condition in Step (1) and does not need to satisfy all of the constraints in Problem 3, but rather the different initial states will influence the iteration times to converge in the sequential convex programming.Remark 5.
The theoretical proof for the convergence of the sequential convex programming is very challenging owing to the presence of various constraints, nonlinear dynamics, and kinematics. Lu and Liu [34] and Wang and Grant [35, 36] provided several examples of sequential convex programming convergence. The numerical simulations presented in the next part demonstrate the convergence and validity of the proposed algorithm.


## 4. Numerical Simulations

In this section, first, we verify the equivalence of the original problem and new problem with the new control variable by comparison to the nonlinear solver. Simultaneously, the convergence and validity of the proposed sequential convex programming algorithm are also be verified. Next, the effectiveness of the proposed algorithm is presented for different impact angles and impact velocities, respectively.To begin with, we set the trust region asδ=500010050π/18020π/180 in equation (19) and the convergence condition as ξ=510.1π/1800.1π/180 in equation (30). The simulations are implemented in MATLAB R2017b on a laptop computer equipped with Intel Core i5-8250U 1.6 GHz and 8 GB RAM. The convex programming problems are solved using MOSEK [37] in CVX [38].The vertical interception engagement of an air-to-ground missile is considered as follows. The missile initial positionx0,y0=0,8000; the target position xf,yf=18000,0; the missile initial velocity Vt0=200m/s. The initial flight path angle θt0=0∘; the initial angle of attack αt0=0∘. Additionally, the limits on the angle of attack and angle of attack rate are αmax=20∘ and α˙max=4∘/s, respectively.


### 4.1. Validity and Convergence of Proposed Algorithm

In this part, we set the impact angleθf∗=−60∘ and impact velocity Vf∗=250m/s. To validate the equivalence and validity of the proposed algorithm, we used the independent Gauss Pseudospectral Optimization Software (GPOPS) [39] to directly solve the original Problem 1.The performance index values of the proposed algorithm in each iteration are presented in Figure1(a). This indicates that, for the considered problem, the sequential algorithm converged in 13 steps without an initial guess. To further illustrate the converge process, more quantitative results are presented in Table 1 and reveal that the difference between each iteration can decrease below the tolerance value. If looser tolerance is used or an appropriate initial guess is made, fewer steps are required for convergence. Each step requires approximately 0.3–0.5 s of CPU time. In summary, the proposed algorithm requires 5.73 s to solve Problem 3, while the GPOPS requires 43.32 s to solve the original Problem 1. The calculation time substantially decreases if a high-performance CPU is used; therefore, the proposed algorithm can be effectively used online in practical situations.Figure 1
(a) Performance index values of each iteration; (b) trajectory profiles; (c) flight path angle profiles; (d) velocity profiles; (e) angle of attack profiles obtained by proposed algorithm and GPOPS; (f) angle of attack rate profile obtained by the proposed algorithm.
(a)
(b)
(c)
(d)
(e)
(f)Table 1
The difference of the states between each iteration for Case A.
Iteration number
Δ
y
m
Δ
V
m
/
s
Δ
θ
deg
Δ
α
deg
1
1609.9
99.583
34.035
9.9714
2
1153.1
22.390
20.542
4.9468
3
392.14
26.228
5.8535
4.6356
4
109.04
7.3673
2.0949
0.7755
5
172.03
7.4144
3.1740
0.9927
6
24.083
1.4514
0.6722
1.1575
7
12.165
0.4772
0.1826
0.0942
8
6.8327
0.3021
0.1312
0.0663
9
7.9768
0.5660
0.2726
0.3389
10
6.1521
0.5061
0.2743
0.4331
11
6.3616
0.5203
0.2565
0.4382
12
6.4974
0.5345
0.2853
0.4701
13
0.7938
0.0694
0.0382
0.0506Figures1(b)–1(e) shows the trajectory profiles, flight path angle profiles, velocity profiles, and angle of attack profiles obtained by the proposed algorithm and GPOPS, respectively. As can be seen, both the proposed sequential convex programming and GPOPS successfully satisfied the impact angle and impact velocity constraints and the additional constraints. Note that the solution obtained by the GPOPS is approximately identical to that obtained by the proposed algorithm, which proves that our approach of handling and transforming the problem is valid. Therefore, Problem 3 is equivalent to the original problem, and Lemma 1 is verified. The small difference between the proposed algorithm and the GPOPS may have occurred for several reasons, such as the linearization and approximation of differential equations or different discretization strategies. Figure 1(f) shows the value of the angle of attack rate obtained by the proposed method. However, the chattering phenomena of the new selected control variable, which were caused by the linearization strategy, did not significantly influence the angle of attack obtained in this study. In other words, the actual control variable angle of attack did not exhibit chattering phenomena and can thus be used in practical interception engagement. Generally, the proposed algorithm can obtain a successful solution to the optimization problem without an initial guess and faster than the GPOPS.


### 4.2. Effectiveness of Proposed Algorithm for Different Impact Angles

To validate that the proposed algorithm can successfully achieve the desired impact angles, we set different impact angle parameters, as follows:θf∗=−45∘,−60∘,−75∘; the impact velocity is the same as in the previous situation. The simulation results, including the trajectory profiles, velocity profiles, flight path angle profiles, angle of attack profiles, and angle of attack rate profiles of the three different impact angles, are presented in Figure 2. The flight path angle in Figure 2(b) confirms that the proposed algorithm can satisfy different impact angle constraints. Figure 2(c) shows that the air-to-ground missile can intercept a target with the desired impact velocity. Additionally, all other constraints are satisfied. Figure 2(e) shows that the biggest and smallest angles of attack are proportional to the impact angles. The angle of attack rate in Figure 2(e) confirms that the variation of the attack angle, large angle of attack at the beginning, and large negative angle of attack at the end of the engagement can achieve a large impact angle, which conforms to the physical law.Figure 2
(a) Trajectory profiles; (b) flight path angle profiles; (c) velocity profiles; (d) attack angle; (e) angle of attack rate profiles of different impact angles.
(a)
(b)
(c)
(d)
(e)


### 4.3. Effectiveness of Proposed Algorithm for Different Impact Velocities

The velocity constraint is another important part of the proposed algorithm. In this section, we set the different impact velocities asVf∗=220m/s,250m/s,280m/s; the impact angle is set as θf∗=−60∘. Figure 3 presents the trajectory profiles, flight path angle profiles, velocity profiles, attack angles profiles, and angle of attack rate profiles for different impact velocities. Figures 3(b)–3(e) show that the proposed algorithm can obtain a maximum impact velocity and achieve the desired near-vertical impact angle. Moreover, it satisfies the angle of attack and angle of attack rate constraints. The trajectory shown in Figure 3(a) indicates that the missile requires more flight distance to slow down. The flight path angle profiles shown in Figure 3(b) result in a large flight distance for a lower impact velocity. Figure 3(d) shows that the high impact velocity requires a big angle of attack at the beginning of flight engagement, which results in a small flight path angle. According to equation (4), a big angle of attack is required because the effect of gravity is countered with the lift generated by the attack angle. Figure 3(e) shows the angle of attack rate profiles coinciding with the angle of attack variation and satisfying the maximum angle of attack rate constraint.Figure 3
(a) Trajectory profiles; (b) flight path angle profiles; (c) velocity profiles; (d) attack angle; (e) angle of attack rate profiles of different impact velocities.
(a)
(b)
(c)
(d)
(e)


### 4.4. Maximum Velocity for Near-Vertical Final Impact Angle

Let us consider a more realistic practical application, wherein the vertical attack with a maximum velocity scenario is attractive. In this part, the objective function is selected as follows:(31)J=−Vxf.The maximum velocity is achieved when the objective function is minimized. The final impact angle is set asθf∗=−89∘. Figure 4 shows the trajectory profiles, flight path angle profiles, velocity profiles, attack angle profiles, and angle of attack rate profiles for maximum impact velocity. In this case, the maximum impact velocity shown in Figures 4(b)–4(e) indicates that the proposed algorithm can obtain the maximum impact velocity with the desired impact angle and can also satisfy the angle of attack and angle of attack rate constraints. As can be seen in Figure 4(c), the maximum impact velocity is 273.1 m/s. The trajectory shown in Figure 4(a) demonstrates that the missile requires a larger flight distance at the end of the trajectory to speed up. The flight path angle profiles in Figure 4(b) show that the flight path angle is approximately 90° at the end of the trajectory. Figure 4(d) shows that the angle of attack for achieving the maximum impact velocity is similar to that in the previous case. As can be seen in Figure 4(e), the limits on the angle of the attack rate are satisfied.Figure 4
(a) Trajectory profiles; (b) flight path angle profiles; (c) velocity profiles; (d) attack angle; (e) angle of attack rate profiles of different impact velocities.
(a)
(b)
(c)
(d)
(e)


## 4.1. Validity and Convergence of Proposed Algorithm

In this part, we set the impact angleθf∗=−60∘ and impact velocity Vf∗=250m/s. To validate the equivalence and validity of the proposed algorithm, we used the independent Gauss Pseudospectral Optimization Software (GPOPS) [39] to directly solve the original Problem 1.The performance index values of the proposed algorithm in each iteration are presented in Figure1(a). This indicates that, for the considered problem, the sequential algorithm converged in 13 steps without an initial guess. To further illustrate the converge process, more quantitative results are presented in Table 1 and reveal that the difference between each iteration can decrease below the tolerance value. If looser tolerance is used or an appropriate initial guess is made, fewer steps are required for convergence. Each step requires approximately 0.3–0.5 s of CPU time. In summary, the proposed algorithm requires 5.73 s to solve Problem 3, while the GPOPS requires 43.32 s to solve the original Problem 1. The calculation time substantially decreases if a high-performance CPU is used; therefore, the proposed algorithm can be effectively used online in practical situations.Figure 1
(a) Performance index values of each iteration; (b) trajectory profiles; (c) flight path angle profiles; (d) velocity profiles; (e) angle of attack profiles obtained by proposed algorithm and GPOPS; (f) angle of attack rate profile obtained by the proposed algorithm.
(a)
(b)
(c)
(d)
(e)
(f)Table 1
The difference of the states between each iteration for Case A.
Iteration number
Δ
y
m
Δ
V
m
/
s
Δ
θ
deg
Δ
α
deg
1
1609.9
99.583
34.035
9.9714
2
1153.1
22.390
20.542
4.9468
3
392.14
26.228
5.8535
4.6356
4
109.04
7.3673
2.0949
0.7755
5
172.03
7.4144
3.1740
0.9927
6
24.083
1.4514
0.6722
1.1575
7
12.165
0.4772
0.1826
0.0942
8
6.8327
0.3021
0.1312
0.0663
9
7.9768
0.5660
0.2726
0.3389
10
6.1521
0.5061
0.2743
0.4331
11
6.3616
0.5203
0.2565
0.4382
12
6.4974
0.5345
0.2853
0.4701
13
0.7938
0.0694
0.0382
0.0506Figures1(b)–1(e) shows the trajectory profiles, flight path angle profiles, velocity profiles, and angle of attack profiles obtained by the proposed algorithm and GPOPS, respectively. As can be seen, both the proposed sequential convex programming and GPOPS successfully satisfied the impact angle and impact velocity constraints and the additional constraints. Note that the solution obtained by the GPOPS is approximately identical to that obtained by the proposed algorithm, which proves that our approach of handling and transforming the problem is valid. Therefore, Problem 3 is equivalent to the original problem, and Lemma 1 is verified. The small difference between the proposed algorithm and the GPOPS may have occurred for several reasons, such as the linearization and approximation of differential equations or different discretization strategies. Figure 1(f) shows the value of the angle of attack rate obtained by the proposed method. However, the chattering phenomena of the new selected control variable, which were caused by the linearization strategy, did not significantly influence the angle of attack obtained in this study. In other words, the actual control variable angle of attack did not exhibit chattering phenomena and can thus be used in practical interception engagement. Generally, the proposed algorithm can obtain a successful solution to the optimization problem without an initial guess and faster than the GPOPS.


## 4.2. Effectiveness of Proposed Algorithm for Different Impact Angles

To validate that the proposed algorithm can successfully achieve the desired impact angles, we set different impact angle parameters, as follows:θf∗=−45∘,−60∘,−75∘; the impact velocity is the same as in the previous situation. The simulation results, including the trajectory profiles, velocity profiles, flight path angle profiles, angle of attack profiles, and angle of attack rate profiles of the three different impact angles, are presented in Figure 2. The flight path angle in Figure 2(b) confirms that the proposed algorithm can satisfy different impact angle constraints. Figure 2(c) shows that the air-to-ground missile can intercept a target with the desired impact velocity. Additionally, all other constraints are satisfied. Figure 2(e) shows that the biggest and smallest angles of attack are proportional to the impact angles. The angle of attack rate in Figure 2(e) confirms that the variation of the attack angle, large angle of attack at the beginning, and large negative angle of attack at the end of the engagement can achieve a large impact angle, which conforms to the physical law.Figure 2
(a) Trajectory profiles; (b) flight path angle profiles; (c) velocity profiles; (d) attack angle; (e) angle of attack rate profiles of different impact angles.
(a)
(b)
(c)
(d)
(e)


## 4.3. Effectiveness of Proposed Algorithm for Different Impact Velocities

The velocity constraint is another important part of the proposed algorithm. In this section, we set the different impact velocities asVf∗=220m/s,250m/s,280m/s; the impact angle is set as θf∗=−60∘. Figure 3 presents the trajectory profiles, flight path angle profiles, velocity profiles, attack angles profiles, and angle of attack rate profiles for different impact velocities. Figures 3(b)–3(e) show that the proposed algorithm can obtain a maximum impact velocity and achieve the desired near-vertical impact angle. Moreover, it satisfies the angle of attack and angle of attack rate constraints. The trajectory shown in Figure 3(a) indicates that the missile requires more flight distance to slow down. The flight path angle profiles shown in Figure 3(b) result in a large flight distance for a lower impact velocity. Figure 3(d) shows that the high impact velocity requires a big angle of attack at the beginning of flight engagement, which results in a small flight path angle. According to equation (4), a big angle of attack is required because the effect of gravity is countered with the lift generated by the attack angle. Figure 3(e) shows the angle of attack rate profiles coinciding with the angle of attack variation and satisfying the maximum angle of attack rate constraint.Figure 3
(a) Trajectory profiles; (b) flight path angle profiles; (c) velocity profiles; (d) attack angle; (e) angle of attack rate profiles of different impact velocities.
(a)
(b)
(c)
(d)
(e)


## 4.4. Maximum Velocity for Near-Vertical Final Impact Angle

Let us consider a more realistic practical application, wherein the vertical attack with a maximum velocity scenario is attractive. In this part, the objective function is selected as follows:(31)J=−Vxf.The maximum velocity is achieved when the objective function is minimized. The final impact angle is set asθf∗=−89∘. Figure 4 shows the trajectory profiles, flight path angle profiles, velocity profiles, attack angle profiles, and angle of attack rate profiles for maximum impact velocity. In this case, the maximum impact velocity shown in Figures 4(b)–4(e) indicates that the proposed algorithm can obtain the maximum impact velocity with the desired impact angle and can also satisfy the angle of attack and angle of attack rate constraints. As can be seen in Figure 4(c), the maximum impact velocity is 273.1 m/s. The trajectory shown in Figure 4(a) demonstrates that the missile requires a larger flight distance at the end of the trajectory to speed up. The flight path angle profiles in Figure 4(b) show that the flight path angle is approximately 90° at the end of the trajectory. Figure 4(d) shows that the angle of attack for achieving the maximum impact velocity is similar to that in the previous case. As can be seen in Figure 4(e), the limits on the angle of the attack rate are satisfied.Figure 4
(a) Trajectory profiles; (b) flight path angle profiles; (c) velocity profiles; (d) attack angle; (e) angle of attack rate profiles of different impact velocities.
(a)
(b)
(c)
(d)
(e)


## 5. Conclusion

This paper proposed a sequential convex programming algorithm for air-to-ground missile optimization-energy guidance under impact angle and impact velocity constraints and maximum angle of attack limits. This problem is formulated as an optimization problem, and the angle of attack rate is set as a new control variable. Next, the problem is linearized and discretized to be solved using sequential convex programming. The numerical simulation results are compared with those obtained by the GPOPS to confirm the validity of the transformation and the optimality of the proposed sequential convex programming approach. Finally, the simulation results confirmed that the proposed algorithm can intercept a target at different impact angles and impact velocities.



---

*Source: 2065730-2019-10-23.xml*