# Hybrid Trajectory Optimization Method and Tracking Guidance for Variable-Sweep Missiles

Authors: Zhifang Wei; Yang Cheng; Xiangxiang Guo; Senlin Liu
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

In this paper, an offline hybrid trajectory optimization approach is proposed for variable-sweep missiles to explore the superiority in the diving phase. Aiming at the maximal terminal velocity with the impact angle constraint, the trajectory optimization model is formulated under multiple constraints, and the aerodynamic analysis in different sweep angles is discussed. Unlike only the attack angle used for the optimization process traditionally, the two-variable optimization scheme on both the attack angle and sweep angle is investigated for variable-sweep missiles. Then, the trajectory optimization problem is transformed into the nonlinear programming problem via a hybrid optimization strategy combining the Gauss pseudospectral method and direct shooting method to obtain the high precision and fast convergence solution. Finally, to verify the feasibility of the optimal trajectory under uncertainties, the tracking guidance law is designed on basis of the gain scheduled linear quadratic regulator control. Numerical simulation results reveal not only of the proposed hybrid optimization strategy but also of the superiority of variable-sweep missiles compared with traditional missiles.


---

## Body

## 1. Introduction

Missiles are widely used in fiercely confronted environment and increasingly complex modern warfare for high guidance accuracy, strong penetration capability, and good maneuverability [1]. Traditional missiles adopt the fixed aerodynamic configurations, which lead to the limited flight performances. Thus, the morphing wing technology has been brought up in an effort to offer the well-matched aerodynamic characteristics in different flight states by means of changing the wing shape adaptively, such as stretching the wings to obtain the large lift-to-drag ratio in the cruise phase and shrinking the wings to reduce the drag coefficient in the diving phase [2, 3].As the benefits of morphing wings are so promising, many intelligent morphing concepts have been investigated in recent aerospace researches, which can be roughly categorized as folding wings, variable-span wings, variable-chord, variable-sweep wings, and twisting wings [4, 5]. However, the morphing wing technologies are highly relied on the advanced materials, smart sensors, and flexible structures, which are still the primary challenges nowadays. Different from the other morphing ways, the variable-sweep wing technology has been applied for the aircraft design successfully, such as F-14 Fighter and Tu-160 Bomber [6]. Therefore, in order to enhance the missile’s overall performances effectively, the variable-sweep wing is taken into considerations in the terms of the technology feasibility. Through the servomechanism system, the sweep angles can be regulated to reshape the aerodynamic configurations of the missiles, which also increase the aerodynamic complexity and make it harder for the trajectory optimization than the traditional fixed configuration missiles. At present, the trajectory optimization researches on variable-sweep missiles are quite limited. Therefore, to explore the superiority of the variable-sweep missile, the investigation of this paper mainly focuses on its trajectory optimization design.A well-designed trajectory is of great essence to improve the aircraft’s flight performance, and the trajectory optimization problems are always solved as the optimal control problem [7]. Up to now, extensive researches on the trajectory optimization techniques for aircrafts have been presented in the literatures, which can be generally classified as the direct methods and indirect methods [8, 9]. The typical indirect methods, such as the gradient descent method [10], pattern search algorithm [11, 12], and multiple shooting method [13], achieve good optimization effects but are highly sensitive to the initial guesses and need to handle the constraints with skills. In order to avoid the complex mathematic derivations, the direct method is widely used to solve the trajectory optimization problems with numerical solutions. By parameterization, the continuous-time trajectory optimization problem is transcribed to the nonlinear programming problem (NLP), which can be solved the sequence quadratic program algorithm effectively. The direct shooting method (DSM) [14], particle swarm optimization algorithm [15], genetic algorithm [16], and symplectic iterative algorithm [17, 18] achieve satisfying results for the specific missions. Different from other direct methods, the Gauss pseudospectral method (GPM) has obtained great attention for fast convergence, global optimization, and extensive applications in the aerospace field recently by means of approximating the control and state variables with polynomials [19–22].In previous work [23, 24], the GPM is used for the trajectory optimization for the variable-sweep missiles. However, it is noted that the GPM with polynomials for approximation cannot guarantee the exact solution with highly constrained conditions. On the contrary, the DSM can get the high precision results under the condition of the proper initial guesses. Motivated by the requirements of the trajectory optimization for the variable-sweep missiles with complex trajectory characteristics, an offline hybrid sequential trajectory optimization strategy combining the GPM and DSM is proposed to obtain the high precision and fast convergence solution in this paper. Firstly, the trajectory optimization model for variable-sweep missiles is established under multiple constraints, and the aerodynamic analysis in different sweep angles is discussed in detail. Unlike only the attack angle for optimization process in the traditional way, the two-variable optimization scheme on both the attack angle and sweep angle is studied. Secondly, the trajectory of variable-sweep missiles for the maximal terminal velocity is optimized via the hybrid strategy. The linear quadratic regulator (LQR) control has been successfully used for the quadrotor’s trajectory tracking and achieved great effects in engineering [25, 26]. Thus, in order to verify the feasibility of the optimal trajectory in the presence of uncertainties, the tracking guidance law is designed based on the gain scheduled LQR. Numerical simulation results show not only the effectiveness of the proposed hybrid optimization strategy but also the superiority of variable-sweep missiles compared with the traditional missiles.The outline of this paper is organized as follows. Section2 presents the problem formulation part, including the variable-sweep missile flight dynamics, the aerodynamic analysis, and the two-variable optimization scheme. In Section 3, the hybrid strategy combining DSM and GPM is proposed. The gain scheduled LQR controller for trajectory guidance is developed in Section 4. Numerical simulations are performed in Section 5, and conclusions are discussed in Section 6.


## 2. Problem Formulation

### 2.1. Variable-Sweep Missile Dynamics

The variable-sweep missile with different aerodynamic configurations is illustrated in Figure1, where the sweep angles can be regulated via the servomechanism system. Generally, the small sweep angle mode is used to increase the lift-to-drag ratio in the low speed cruise phase; and the large sweep angle mode is utilized to reduce the drag in the high speed diving phase.Figure 1
Missile configurations with different sweep angles. (a)χ=10°. (b) χ=30°. (c) χ=50°. (d) χ=70°.
(a)(b)(c)(d)To be convenient, suppose that the missile is simplified as the point-mass model, and there exists zero lag control system. Under the condition of nonrotating spherical earth, the planar motion of the variable-sweep missile in the unpowered diving phase is described by the following set ordinary differential equations [23]:(1)dxdt=υcosγ,(2)dhdt=υsinγ,(3)dυdt=−1mD+mgsinγ,(4)dγdt=1mυL−mgcosγ,where x is the range, h is the altitude, υ is the velocity, γ is the flight path angle, m denotes the missile mass, g=g01−2h/Re is the gravitational acceleration, g0=9.80665m/s2 at sea level, Re=6378145m is the earth radius, and L and D represent the lift and drag, respectively, given as follows:(5)D=qSCDMa,α,χ,L=qSCLMa,α,χ,where q=0.5ρυ2 is the dynamic pressure, ρ is the air density referring to the 1976 US standard atmospheric model with the exponential fitting equation ρ=ρ01−h/443004.2533, ρ0=1.225kg/m3 at sea level, S is the aerodynamic reference area, CD and CL denote the drag and lift coefficients, respectively, both related with the Mach number Ma, the angle of attack α, and the sweep angle χ.


### 2.2. Aerodynamic Modeling and Analysis

In order to figure out the reliable aerodynamic characteristics of variable-sweep missiles, the Missile DATCOM software (97 version) is used for the aerodynamic calculation [27]. The data are obtained with different Mach numbers, attack angles, and sweep angles. For aerodynamic analysis with an emphasis on sweep angle’s influences, the aerodynamic characteristics are shown in Figure 2 in Ma=1.2, and in Figure 3 in α=4° in terms of the lift coefficient CL, drag coefficient CD, lift-to-drag ratio L/D, and maximum lift-to-drag ratio L/Dmax.Figure 2
Aerodynamic characteristics with different sweep angles inMa=1.2. (a) Lift coefficient CL. (b) Drag coefficient CD. (c) Lift-to-drag ratio L/D. (d) Maximum lift-to-drag ratio L/Dmax.
(a)(b)(c)(d)Figure 3
Aerodynamic characteristics with different sweep angles inα=4°. (a) Lift coefficient CL. (b) Drag coefficient CD. (c) Lift-to-drag ratio L/D. (d) Maximum lift-to-drag ratio L/Dmax.
(a)(b)(c)(d)As seen from Figures2 and 3, both the lift coefficient CL and drag coefficient CD decrease when the sweep angle χ increases. Besides, it also can be observed that the drag coefficient CD reaches the peak value around the sound velocity, which is confirmed in [28].The aerodynamic data calculated by Missile DATCOM are originally arranged in the lookup table, which can be formulated and identified by the least squares estimation. The drag coefficientCD and the lift coefficient CL in different Mach numbers are modeled as follows:(6)CDMa,α,χ=CD0+CDαα+CDχχ+CDα2α2+CDχα2χα2,CLMa,α,χ=CL0+CLαα+CLχαχα,where the aerodynamic derivatives CD0, CDα, CDχ, CDα2, CDχα2, CL0, CLα, and CLχα are given in Appendix A.Obviously, it is indicated in equation (6) that the terms CDχχ and CDχα2χα2 act as the induced drag part on account of the sweep angle χ, which cannot be neglected in the variable-sweep process. Though the term CLχαχα is contributed to the lift L by the sweep angle χ, it is noted that the variable-sweep in χ=10° has the largest maximum lift-to-drag ratio L/Dmax shown in Figures 2(d) and 3(d). Consequently, the increasing sweep angle χ can induce more drag than lift.


### 2.3. Constraints and Two-Variable Optimization Scheme

To maximize the warhead’s effective damages, the maximal terminal velocity is selected as the optimization index under the constraints of the boundary conditions and path constraints including the control constraints. The cost function is described as(7)minJ=−υf,subject to(8)υ0−υt0=0,x0−xt0=0,xf−xtf=0,h0−ht0=0,hf−htf=0,γ0−γt0=0,γf−γtf=0,ny=Lcosα+Dsinα≤nymax,q≤qmax,αmin≤α≤αmax,χmin≤χ≤χmax,where the subscriptd “0” and “f” under the states υ, x, h, and γ denote the actual initial and final states, respectively, and the subscripts “0” and “f” under the time t mean the given constraints; ny is the normal overload, nymax is the allowed maximum overload, qmax represents the maximum dynamic pressure; αmin and αmax are the minimum and maximum attack angles, and χmin and χmax are the minimum and maximum sweep angles.Since the final timetf is not specified, this optimization problem is classified as the free-final-time problem. To solve this problem, the final time tf is generally taken considerations as being known with an initial guess.Unlike only the attack angle used for optimization process in the traditional way, the two-variable optimization scheme on both the attack angle and sweep angle is presented for variable-sweep missiles. As exhibited in Figure4, the traditional way of missile’s trajectory optimization is intended to search the optimal attack angle curve to increase the range during the whole flight phase, which can be regarded as the one-dimensional optimization, whereas the two-variable optimization scheme for variable-sweep missiles extends the curve into the surface, which can be taken as the two-dimensional optimization to dig out the potential performance greatly.Figure 4
Comparison of different optimization process between (a) traditional missiles and (b) variable-sweep missiles.
(a)(b)


## 2.1. Variable-Sweep Missile Dynamics

The variable-sweep missile with different aerodynamic configurations is illustrated in Figure1, where the sweep angles can be regulated via the servomechanism system. Generally, the small sweep angle mode is used to increase the lift-to-drag ratio in the low speed cruise phase; and the large sweep angle mode is utilized to reduce the drag in the high speed diving phase.Figure 1
Missile configurations with different sweep angles. (a)χ=10°. (b) χ=30°. (c) χ=50°. (d) χ=70°.
(a)(b)(c)(d)To be convenient, suppose that the missile is simplified as the point-mass model, and there exists zero lag control system. Under the condition of nonrotating spherical earth, the planar motion of the variable-sweep missile in the unpowered diving phase is described by the following set ordinary differential equations [23]:(1)dxdt=υcosγ,(2)dhdt=υsinγ,(3)dυdt=−1mD+mgsinγ,(4)dγdt=1mυL−mgcosγ,where x is the range, h is the altitude, υ is the velocity, γ is the flight path angle, m denotes the missile mass, g=g01−2h/Re is the gravitational acceleration, g0=9.80665m/s2 at sea level, Re=6378145m is the earth radius, and L and D represent the lift and drag, respectively, given as follows:(5)D=qSCDMa,α,χ,L=qSCLMa,α,χ,where q=0.5ρυ2 is the dynamic pressure, ρ is the air density referring to the 1976 US standard atmospheric model with the exponential fitting equation ρ=ρ01−h/443004.2533, ρ0=1.225kg/m3 at sea level, S is the aerodynamic reference area, CD and CL denote the drag and lift coefficients, respectively, both related with the Mach number Ma, the angle of attack α, and the sweep angle χ.


## 2.2. Aerodynamic Modeling and Analysis

In order to figure out the reliable aerodynamic characteristics of variable-sweep missiles, the Missile DATCOM software (97 version) is used for the aerodynamic calculation [27]. The data are obtained with different Mach numbers, attack angles, and sweep angles. For aerodynamic analysis with an emphasis on sweep angle’s influences, the aerodynamic characteristics are shown in Figure 2 in Ma=1.2, and in Figure 3 in α=4° in terms of the lift coefficient CL, drag coefficient CD, lift-to-drag ratio L/D, and maximum lift-to-drag ratio L/Dmax.Figure 2
Aerodynamic characteristics with different sweep angles inMa=1.2. (a) Lift coefficient CL. (b) Drag coefficient CD. (c) Lift-to-drag ratio L/D. (d) Maximum lift-to-drag ratio L/Dmax.
(a)(b)(c)(d)Figure 3
Aerodynamic characteristics with different sweep angles inα=4°. (a) Lift coefficient CL. (b) Drag coefficient CD. (c) Lift-to-drag ratio L/D. (d) Maximum lift-to-drag ratio L/Dmax.
(a)(b)(c)(d)As seen from Figures2 and 3, both the lift coefficient CL and drag coefficient CD decrease when the sweep angle χ increases. Besides, it also can be observed that the drag coefficient CD reaches the peak value around the sound velocity, which is confirmed in [28].The aerodynamic data calculated by Missile DATCOM are originally arranged in the lookup table, which can be formulated and identified by the least squares estimation. The drag coefficientCD and the lift coefficient CL in different Mach numbers are modeled as follows:(6)CDMa,α,χ=CD0+CDαα+CDχχ+CDα2α2+CDχα2χα2,CLMa,α,χ=CL0+CLαα+CLχαχα,where the aerodynamic derivatives CD0, CDα, CDχ, CDα2, CDχα2, CL0, CLα, and CLχα are given in Appendix A.Obviously, it is indicated in equation (6) that the terms CDχχ and CDχα2χα2 act as the induced drag part on account of the sweep angle χ, which cannot be neglected in the variable-sweep process. Though the term CLχαχα is contributed to the lift L by the sweep angle χ, it is noted that the variable-sweep in χ=10° has the largest maximum lift-to-drag ratio L/Dmax shown in Figures 2(d) and 3(d). Consequently, the increasing sweep angle χ can induce more drag than lift.


## 2.3. Constraints and Two-Variable Optimization Scheme

To maximize the warhead’s effective damages, the maximal terminal velocity is selected as the optimization index under the constraints of the boundary conditions and path constraints including the control constraints. The cost function is described as(7)minJ=−υf,subject to(8)υ0−υt0=0,x0−xt0=0,xf−xtf=0,h0−ht0=0,hf−htf=0,γ0−γt0=0,γf−γtf=0,ny=Lcosα+Dsinα≤nymax,q≤qmax,αmin≤α≤αmax,χmin≤χ≤χmax,where the subscriptd “0” and “f” under the states υ, x, h, and γ denote the actual initial and final states, respectively, and the subscripts “0” and “f” under the time t mean the given constraints; ny is the normal overload, nymax is the allowed maximum overload, qmax represents the maximum dynamic pressure; αmin and αmax are the minimum and maximum attack angles, and χmin and χmax are the minimum and maximum sweep angles.Since the final timetf is not specified, this optimization problem is classified as the free-final-time problem. To solve this problem, the final time tf is generally taken considerations as being known with an initial guess.Unlike only the attack angle used for optimization process in the traditional way, the two-variable optimization scheme on both the attack angle and sweep angle is presented for variable-sweep missiles. As exhibited in Figure4, the traditional way of missile’s trajectory optimization is intended to search the optimal attack angle curve to increase the range during the whole flight phase, which can be regarded as the one-dimensional optimization, whereas the two-variable optimization scheme for variable-sweep missiles extends the curve into the surface, which can be taken as the two-dimensional optimization to dig out the potential performance greatly.Figure 4
Comparison of different optimization process between (a) traditional missiles and (b) variable-sweep missiles.
(a)(b)


## 3. Hybrid Optimization Strategy for Trajectory Optimization

### 3.1. Description of the Trajectory Optimization Problem

The trajectory optimization of the variable-sweep missile with multiple constraints can be classified as the Bolza optimal control problem, which is described by the cost function:(9)J=Φxt0,t0,xtf,tf+∫t0tfLxt,ut,tdt,where xt∈Rn is the state variable and ut∈Rm is the control variable; t0 and tf are the initial and final time, respectively; Φxt0,t0,xtf,tf is the final weight function; and Lxt,ut,t is the path integral function.The optimization problem is subject to the dynamic constraint, boundary condition, and path constraint as(10)x˙t=fxt,ut,t,t∈t0,tf,(11)gxt0,t0,xtf,tf=0,(12)hxt,ut,t≤0,t∈t0,tf.


### 3.2. Direct Shooting Method

The DSM is the simple and common technique of the direct method to solve the trajectory optimization problem. The way to transform the trajectory optimization problem into the NLP problem for the DSM is to discretize the control variable.The time intervalt0,tf is divided into a set of subintervals as follows:(13)t0<t1<⋯<tk<⋯<tN=tf,k=0,1,…,N.Assume that the guessed control variableut in the discrete points isuk=u0,u1,…,uN. Hence, the guessed control variable ut is obtained by the linear interpolation method:(14)ut=uk−1+t−tk−1tk−tk−1uk+1−uk,tk≤t≤tk+1.Since the trajectory optimization of the variable-sweep missile is the free-final-time problem, denote the design variableU as U=u0T,u1T,…,uNT,tfT. Then, the trajectory optimization problem is transformed into the NLP problem described as follows:(15)minJ=JU,subject to(16)gU=0,hU≤0.


### 3.3. Gauss Pseudospectral Method

The GPM has been proven to solve the trajectory optimization problem effectively. To implement the GPM, both the control and state variables are discretized to transform the trajectory optimization problem into the NLP problem on a set of LegendreGauss (LG) collocation points.The time intervalt0,tf is converted into −1,1 by the affine transformation:(17)τ=2ttf−t0−tf+t0tf−t0.DefineK as the collocation point number, which also indicates the rank of the Legendre polynomials. The state and control variables are approximated by Lagrange interpolating basis polynomials Liτ=∏j=0,j≠iKτ−τj/τi−τj,i=0,1,…,K with the initial point τ0=−1:(18)xτ≈Xτ=∑i=0KLiτxτi,(19)uτ≈Uτ=∑i=0KLiτuτi.Differentiating equation (10) and combining equation (18), the dynamic constraint is approximated on each LG point in the algebraic form(20)x˙τk≈X˙τk=∑i=0KL˙iτkXτi,k=0,1,…,K,under the condition of the constraint expressed by(21)∑i=0KL˙iτkXτi−tf−t02fXτk,Uτk,τ;t0,tf=0.For the terminal constraint, as the terminal pointτf=1 is not included in the collocation points, it should also satisfy the dynamic constraint equation (10), which is expressed by the Gauss quadrature in the algebraic form:(22)Xτf=Xτ0+τf−τ02∑k=1KωkfXτk,Uτk,τ;t0,tf,where ωk is the Gauss weight function.In the same way, the cost function equation (9), the boundary conditions equation (11) and the path constraint equation (12) are approximated in the following algebraic form:(23)J=ΦX0,t0,Xf,tf+tf−t02∑k=1KωkLXk,Uk,τk;t0,tf,gX0,t0,Xf,tf=0,hXk,Uk,τk;t0,tf≤0.


### 3.4. Hybrid Optimization Strategy

It should be pointed out that the trajectory optimization of the variable-sweep missiles studied in this paper is used for the offline design. In previous researches on variable-sweep missiles [23, 24], the GPM is used for the trajectory design with fast convergence and global optimization. However, as the GPM is approximated by polynomials, it cannot guarantee the exact solution with highly constrained conditions, especially for the variable-sweep missiles with complex trajectory characteristics. By contrast, the DSM provides the accurate solution without approximation under the condition of a satisfying initial guess. In order to obtain the high precision and fast convergence optimized trajectory of variable-sweep missiles, a hybrid sequential optimization strategy combining the GPM and DSM is presented.By choosing the LG collocation points, the initial solution is generated via the GPM. As this optimization problem is the free-final-time problem, the control variables and the final timetf are both included in the initial solution. Then, taking the initial solution into the DSM, the accurate solution can be obtained with the proper discrete points.It should be pointed out that it is hard for the DSM to deal with the equality constraints, such as the boundary conditions. Consequently, the terminal boundary conditions in equation (8) are converted into the inequalities by adding the thresholds εx, εh, and εγ expressed as(24)xf−xtf≤εx,hf−htf≤εh,γf−γtf≤εγ.


## 3.1. Description of the Trajectory Optimization Problem

The trajectory optimization of the variable-sweep missile with multiple constraints can be classified as the Bolza optimal control problem, which is described by the cost function:(9)J=Φxt0,t0,xtf,tf+∫t0tfLxt,ut,tdt,where xt∈Rn is the state variable and ut∈Rm is the control variable; t0 and tf are the initial and final time, respectively; Φxt0,t0,xtf,tf is the final weight function; and Lxt,ut,t is the path integral function.The optimization problem is subject to the dynamic constraint, boundary condition, and path constraint as(10)x˙t=fxt,ut,t,t∈t0,tf,(11)gxt0,t0,xtf,tf=0,(12)hxt,ut,t≤0,t∈t0,tf.


## 3.2. Direct Shooting Method

The DSM is the simple and common technique of the direct method to solve the trajectory optimization problem. The way to transform the trajectory optimization problem into the NLP problem for the DSM is to discretize the control variable.The time intervalt0,tf is divided into a set of subintervals as follows:(13)t0<t1<⋯<tk<⋯<tN=tf,k=0,1,…,N.Assume that the guessed control variableut in the discrete points isuk=u0,u1,…,uN. Hence, the guessed control variable ut is obtained by the linear interpolation method:(14)ut=uk−1+t−tk−1tk−tk−1uk+1−uk,tk≤t≤tk+1.Since the trajectory optimization of the variable-sweep missile is the free-final-time problem, denote the design variableU as U=u0T,u1T,…,uNT,tfT. Then, the trajectory optimization problem is transformed into the NLP problem described as follows:(15)minJ=JU,subject to(16)gU=0,hU≤0.


## 3.3. Gauss Pseudospectral Method

The GPM has been proven to solve the trajectory optimization problem effectively. To implement the GPM, both the control and state variables are discretized to transform the trajectory optimization problem into the NLP problem on a set of LegendreGauss (LG) collocation points.The time intervalt0,tf is converted into −1,1 by the affine transformation:(17)τ=2ttf−t0−tf+t0tf−t0.DefineK as the collocation point number, which also indicates the rank of the Legendre polynomials. The state and control variables are approximated by Lagrange interpolating basis polynomials Liτ=∏j=0,j≠iKτ−τj/τi−τj,i=0,1,…,K with the initial point τ0=−1:(18)xτ≈Xτ=∑i=0KLiτxτi,(19)uτ≈Uτ=∑i=0KLiτuτi.Differentiating equation (10) and combining equation (18), the dynamic constraint is approximated on each LG point in the algebraic form(20)x˙τk≈X˙τk=∑i=0KL˙iτkXτi,k=0,1,…,K,under the condition of the constraint expressed by(21)∑i=0KL˙iτkXτi−tf−t02fXτk,Uτk,τ;t0,tf=0.For the terminal constraint, as the terminal pointτf=1 is not included in the collocation points, it should also satisfy the dynamic constraint equation (10), which is expressed by the Gauss quadrature in the algebraic form:(22)Xτf=Xτ0+τf−τ02∑k=1KωkfXτk,Uτk,τ;t0,tf,where ωk is the Gauss weight function.In the same way, the cost function equation (9), the boundary conditions equation (11) and the path constraint equation (12) are approximated in the following algebraic form:(23)J=ΦX0,t0,Xf,tf+tf−t02∑k=1KωkLXk,Uk,τk;t0,tf,gX0,t0,Xf,tf=0,hXk,Uk,τk;t0,tf≤0.


## 3.4. Hybrid Optimization Strategy

It should be pointed out that the trajectory optimization of the variable-sweep missiles studied in this paper is used for the offline design. In previous researches on variable-sweep missiles [23, 24], the GPM is used for the trajectory design with fast convergence and global optimization. However, as the GPM is approximated by polynomials, it cannot guarantee the exact solution with highly constrained conditions, especially for the variable-sweep missiles with complex trajectory characteristics. By contrast, the DSM provides the accurate solution without approximation under the condition of a satisfying initial guess. In order to obtain the high precision and fast convergence optimized trajectory of variable-sweep missiles, a hybrid sequential optimization strategy combining the GPM and DSM is presented.By choosing the LG collocation points, the initial solution is generated via the GPM. As this optimization problem is the free-final-time problem, the control variables and the final timetf are both included in the initial solution. Then, taking the initial solution into the DSM, the accurate solution can be obtained with the proper discrete points.It should be pointed out that it is hard for the DSM to deal with the equality constraints, such as the boundary conditions. Consequently, the terminal boundary conditions in equation (8) are converted into the inequalities by adding the thresholds εx, εh, and εγ expressed as(24)xf−xtf≤εx,hf−htf≤εh,γf−γtf≤εγ.


## 4. Gain Scheduled LQR Controller for Tracking Guidance

After the optimal trajectory of the variable-sweep missile is obtained via the hybrid optimization strategy, it appears to be feasible and realizable in the open-loop tracking guidance. However, there exist many uncertainty factors in reality, such as the wind and initial state biases. Compared with the model predictive control for the reference tracking in [29], the LQR control has been successfully used for the quadrotor’s trajectory tracking and achieved great effects in engineering. In order to demonstrate the feasibility of the optimal trajectory of the variable-sweep missile under uncertainties, the gain-scheduled LQR controller is designed for tracking guidance in this paper.


### 4.1. Dynamics Linearization

The optimal trajectory generated by the hybrid optimization algorithm is taken as the reference trajectory. As the missile velocityυ changes fast in the diving phase and cannot be controlled directly, the altitude h and the flight path angle γ are chosen as the feedback control terms to track the reference trajectory in this paper. For the better tracking robustness, the range is chosen as the independent variable. The independent variable transformation between the range and time is expressed as follows:(25)d∗dx=d∗dx⋅dtdt=d∗dt1υcosγ.Thus, the dynamics of variable-sweep missiles in equations (2) and (4) can be rewritten as(26)dhdx=f1h,γ,α,χ,dγdx=f2h,γ,α,χ,where f1h,γ,α,χ=tanγ and f2h,γ,α,χ=1/mυ2L/cosγ−mg.Define the tracking error vectorz˜=h˜γ˜T and control vector u˜=α˜χ˜T, where h˜=h−href is the tracking error with respect to the reference altitude href and γ˜=γ−γref is the tracking error with respect to the reference flight path angle γref. Based on the small-perturbation theory, the linearized dynamics equation of the variable-sweep missile in the state space form is written as(27)dz˜dx=∂f1∂h∂f1∂γ∂f2∂h∂f2∂γz+∂f1∂α∂f1∂χ∂f2∂α∂f2∂χu˜=Az˜+Bu˜,where the detailed expressions of the Jacobian matrices A and B are given in the Appendix B.


### 4.2. Gain-Scheduled LQR Controller

For the given linear time-varying system in equation (27), the gain-scheduled LQR controller provides the optimal solution. The control law takes the following form:(28)u=−Hz˜+uref,where uref=αrefχrefT is the reference control vector and H is the feedback gain matrix.In order to determine the optimal gain matrixH, the quadratic cost function is defined as(29)Jz˜,u˜=∫t0tfz˜TQz˜+u˜TRu˜dt,where Q and R are the positive-definite matrices. Assume that the optimal quadratic cost function is J∗z˜=z˜TPz˜, and the feedback gain matrix H=−R−1BTP is obtained by solving the Ricatti equation:(30)ATP+PA−PBR−1BTP+Q=0.It is noted that both the reference trajectory and LQR gain matrices are obtained offline. Like the gain-scheduled technique for PID controller in engineering, the feedback gain matrices can be generated via the interpolation in real time.


## 4.1. Dynamics Linearization

The optimal trajectory generated by the hybrid optimization algorithm is taken as the reference trajectory. As the missile velocityυ changes fast in the diving phase and cannot be controlled directly, the altitude h and the flight path angle γ are chosen as the feedback control terms to track the reference trajectory in this paper. For the better tracking robustness, the range is chosen as the independent variable. The independent variable transformation between the range and time is expressed as follows:(25)d∗dx=d∗dx⋅dtdt=d∗dt1υcosγ.Thus, the dynamics of variable-sweep missiles in equations (2) and (4) can be rewritten as(26)dhdx=f1h,γ,α,χ,dγdx=f2h,γ,α,χ,where f1h,γ,α,χ=tanγ and f2h,γ,α,χ=1/mυ2L/cosγ−mg.Define the tracking error vectorz˜=h˜γ˜T and control vector u˜=α˜χ˜T, where h˜=h−href is the tracking error with respect to the reference altitude href and γ˜=γ−γref is the tracking error with respect to the reference flight path angle γref. Based on the small-perturbation theory, the linearized dynamics equation of the variable-sweep missile in the state space form is written as(27)dz˜dx=∂f1∂h∂f1∂γ∂f2∂h∂f2∂γz+∂f1∂α∂f1∂χ∂f2∂α∂f2∂χu˜=Az˜+Bu˜,where the detailed expressions of the Jacobian matrices A and B are given in the Appendix B.


## 4.2. Gain-Scheduled LQR Controller

For the given linear time-varying system in equation (27), the gain-scheduled LQR controller provides the optimal solution. The control law takes the following form:(28)u=−Hz˜+uref,where uref=αrefχrefT is the reference control vector and H is the feedback gain matrix.In order to determine the optimal gain matrixH, the quadratic cost function is defined as(29)Jz˜,u˜=∫t0tfz˜TQz˜+u˜TRu˜dt,where Q and R are the positive-definite matrices. Assume that the optimal quadratic cost function is J∗z˜=z˜TPz˜, and the feedback gain matrix H=−R−1BTP is obtained by solving the Ricatti equation:(30)ATP+PA−PBR−1BTP+Q=0.It is noted that both the reference trajectory and LQR gain matrices are obtained offline. Like the gain-scheduled technique for PID controller in engineering, the feedback gain matrices can be generated via the interpolation in real time.


## 5. Numerical Simulations

In this section, numerical simulations are implemented to validate the performances of both the proposed hybrid trajectory optimization strategy and the gain scheduled LQR controller for tracking guidance. The setup parameters including the initial states and final states of the missile are given in Table1, and the path constraints are nymax=4 and qmax=100kPa. The control constraints are α∈−10°10° and χ∈10°80°. Besides, the missile mass is set as 1200 kg.Table 1
Parameter setup of variable missile.
Parameterxt0ht0υt0γt0Value0 m5000 m440 m/s0°ParameterxtfhtfυtfγtfValue8500 m0 mFree−60°


### 5.1. Verification of the Proposed Hybrid Optimization Strategy

The GPM is used for comparison with the proposed hybrid optimization strategy under the same simulation conditions, and 30 LG collocation points are employed. For the hybrid optimization strategy, the results generated by the GPM are used as the initial guesses for the DSM. Besides, 50 discrete points are chosen for the DSM, and the fourth-order RungeKutta method is used for numerical quadrature with the fixed step 5 ms. The thresholds for the DSM are set asεx=1m, εh=1m, and εγ=1°.During this simulation, the sweep angle of the missile keeps atχ=10°, and only the attack angle is employed for the optimization process. The results are shown in Table 2 and Figure 5.Table 2
Final states for GPM and hybrid optimization strategy.
Methodxfmhfmυfm/sγf°GPM8503.520324.48−59.97Hybrid8500.000324.50−60.00Figure 5
Altitude along the range for GPM and hybrid optimization strategy.As displayed from Table2, in contrast with the GPM, the hybrid optimization strategy can meet all the given constraints and guarantee the high precision solution.


### 5.2. Comparison of Fixed Configurations with the Variable-Sweep Configuration

In order to verify the superiority of variable-sweep missiles compared with traditional missiles, the hybrid optimization strategy is applied for the trajectory optimization. The fixed configurations, whose sweep angles are set as10°, 30°, 60°, and 80°, respectively, only use the attack angle for optimization process. The variable-sweep configuration can change the sweep angle from 10° and 80°, and the two-variable optimization scheme on the attack angle and sweep angle is adopted. The simulation results are shown in Figures 6 and 7, and the terminal velocities of different configurations are displayed in Table 3.Figure 6
State variables for different configurations. (a) Curve of velocity. (b) Curve of final velocity. (c) Curve of altitude along range. (d) Curve of flight path angle.
(a)(b)(c)(d)Figure 7
Control variables for different configurations. (a) Curve of attack angle. (b) Curve of sweep angle.
(a)(b)Table 3
Terminal velocities of different configurations.
Configurationsχ=10°χ=30°χ=60°χ=80°Variable-sweepυfm/s324.50328.82335.95339.34340.85Figure6(b) and Table 3 reveal that the terminal velocity is highly related with the sweep angle, or rather, the missile aerodynamic configuration. As analyzed in Section 2.2, the drag coefficient decreases when the sweep angle increases. Thus, the fixed configuration with χ=80° has the largest terminal velocity reaching 339.34 m/s compared with the others. However, on basis of the two-variable optimization scheme for the variable-sweep configuration, the terminal velocity can reach 340.85 m/s. Figure 7 presents the control variables of different configurations, and it is found in Figure 7(b) that the sweep angle of the variable-sweep configuration changes fast during 10 s–15 s, when the velocity in Figure 6(a) is around the sound speed and the drag coefficient is higher than the other phases.The results indicate that, in contrast with the fixed configurations, the variable-sweep missile can adaptively change the sweep angles to provide the optimal aerodynamic characteristics in the light of the current flight states.


### 5.3. Tracking Guidance under Uncertainties

The gain-scheduled LQR controller is used as the closed-loop scheme to track the optimal trajectory generated by the hybrid optimization strategy under uncertainties accurately. To demonstrate the effectiveness, the open-loop scheme is taken for comparison.In simulation, the parameter matrices for the LQR controller are set asQ=diag1000,1000 and R=diag1000,1000. The wind and initial state biases are considered as the uncertainties. The wind model used in [30] shows that the wind speed varies with the altitude. The wind speed on the ground is set as 5m/s, and initial state biases are set as Δυ0=−10m/s, Δh0=10m, and Δγ0=3°. The simulation results are shown in Figure 8, and the terminal state biases compared with the reference trajectory are displayed in Table 4.Figure 8
State variables for different tracking guidance schemes. (a) Curve of velocity. (b) Curve of altitude along range. (c) Curve of flight path angle. (d) Curve of tracking error.
(a)(b)(c)(d)Table 4
Terminal state biases compared with reference trajectory.
GuidanceΔxfmΔhfmΔυfm/sΔγf°Open-loop224.510−0.48−1.07LQR−0.340−0.780.00Apparently, the open-loop scheme cannot track the reference trajectory accurately with the terminal range bias achieving 224.51 m in the presence of uncertainties. By contrast, the LQR controller keeps the tracking errors bounded, such as the tracking errorγ˜ of the flight path angle depicted in Figure 8(d).To verify the robustness the LQR controller for tracking guidance under uncertainties, the wind speed on the ground is set as within±5m/s randomly, and initial state biases are all set as within Δυ0=±10m/s, Δh0=±10m, and Δγ0=±3° randomly. The Monte-Carlo simulations of 100 times are carried out. The dispersions of the terminal states are illustrated in Figure 9, and the statistical results of the terminal states in terms of the mean values and standard deviations are exhibited in Table 5.Figure 9
Monte-Carlo simulation results for different tracking guidance schemes. Dispersion of (a) terminal velocity, and (b) terminal range, (b) terminal flight path angle.
(a)(b)(c)Table 5
Statistical results of Monte-Carlo simulations for tracking guidance.
Statistical resultsOpen-loopLQRmeanxfm8602.458499.90stdxfm327.340.26meanυfm/s340.69340.70stdυfm/s0.530.43meanγf°−59.54−59.97stdγf°1.590.00According to the Monte-Carlo simulation results, it can be observed that the impact points for the open-loop scheme deviate from the reference trajectory with the mean value 8602.45 m and standard deviation 327.34 m. The LQR controller for tracking guidance has a distinct advantage over the open-loop scheme in terms of the impact accuracy and constrained angle under uncertainties. To conclude, the effectiveness and robustness of the gained schemed LQR controller for tracking guidance are verified under uncertainties in the numerical simulations.


## 5.1. Verification of the Proposed Hybrid Optimization Strategy

The GPM is used for comparison with the proposed hybrid optimization strategy under the same simulation conditions, and 30 LG collocation points are employed. For the hybrid optimization strategy, the results generated by the GPM are used as the initial guesses for the DSM. Besides, 50 discrete points are chosen for the DSM, and the fourth-order RungeKutta method is used for numerical quadrature with the fixed step 5 ms. The thresholds for the DSM are set asεx=1m, εh=1m, and εγ=1°.During this simulation, the sweep angle of the missile keeps atχ=10°, and only the attack angle is employed for the optimization process. The results are shown in Table 2 and Figure 5.Table 2
Final states for GPM and hybrid optimization strategy.
Methodxfmhfmυfm/sγf°GPM8503.520324.48−59.97Hybrid8500.000324.50−60.00Figure 5
Altitude along the range for GPM and hybrid optimization strategy.As displayed from Table2, in contrast with the GPM, the hybrid optimization strategy can meet all the given constraints and guarantee the high precision solution.


## 5.2. Comparison of Fixed Configurations with the Variable-Sweep Configuration

In order to verify the superiority of variable-sweep missiles compared with traditional missiles, the hybrid optimization strategy is applied for the trajectory optimization. The fixed configurations, whose sweep angles are set as10°, 30°, 60°, and 80°, respectively, only use the attack angle for optimization process. The variable-sweep configuration can change the sweep angle from 10° and 80°, and the two-variable optimization scheme on the attack angle and sweep angle is adopted. The simulation results are shown in Figures 6 and 7, and the terminal velocities of different configurations are displayed in Table 3.Figure 6
State variables for different configurations. (a) Curve of velocity. (b) Curve of final velocity. (c) Curve of altitude along range. (d) Curve of flight path angle.
(a)(b)(c)(d)Figure 7
Control variables for different configurations. (a) Curve of attack angle. (b) Curve of sweep angle.
(a)(b)Table 3
Terminal velocities of different configurations.
Configurationsχ=10°χ=30°χ=60°χ=80°Variable-sweepυfm/s324.50328.82335.95339.34340.85Figure6(b) and Table 3 reveal that the terminal velocity is highly related with the sweep angle, or rather, the missile aerodynamic configuration. As analyzed in Section 2.2, the drag coefficient decreases when the sweep angle increases. Thus, the fixed configuration with χ=80° has the largest terminal velocity reaching 339.34 m/s compared with the others. However, on basis of the two-variable optimization scheme for the variable-sweep configuration, the terminal velocity can reach 340.85 m/s. Figure 7 presents the control variables of different configurations, and it is found in Figure 7(b) that the sweep angle of the variable-sweep configuration changes fast during 10 s–15 s, when the velocity in Figure 6(a) is around the sound speed and the drag coefficient is higher than the other phases.The results indicate that, in contrast with the fixed configurations, the variable-sweep missile can adaptively change the sweep angles to provide the optimal aerodynamic characteristics in the light of the current flight states.


## 5.3. Tracking Guidance under Uncertainties

The gain-scheduled LQR controller is used as the closed-loop scheme to track the optimal trajectory generated by the hybrid optimization strategy under uncertainties accurately. To demonstrate the effectiveness, the open-loop scheme is taken for comparison.In simulation, the parameter matrices for the LQR controller are set asQ=diag1000,1000 and R=diag1000,1000. The wind and initial state biases are considered as the uncertainties. The wind model used in [30] shows that the wind speed varies with the altitude. The wind speed on the ground is set as 5m/s, and initial state biases are set as Δυ0=−10m/s, Δh0=10m, and Δγ0=3°. The simulation results are shown in Figure 8, and the terminal state biases compared with the reference trajectory are displayed in Table 4.Figure 8
State variables for different tracking guidance schemes. (a) Curve of velocity. (b) Curve of altitude along range. (c) Curve of flight path angle. (d) Curve of tracking error.
(a)(b)(c)(d)Table 4
Terminal state biases compared with reference trajectory.
GuidanceΔxfmΔhfmΔυfm/sΔγf°Open-loop224.510−0.48−1.07LQR−0.340−0.780.00Apparently, the open-loop scheme cannot track the reference trajectory accurately with the terminal range bias achieving 224.51 m in the presence of uncertainties. By contrast, the LQR controller keeps the tracking errors bounded, such as the tracking errorγ˜ of the flight path angle depicted in Figure 8(d).To verify the robustness the LQR controller for tracking guidance under uncertainties, the wind speed on the ground is set as within±5m/s randomly, and initial state biases are all set as within Δυ0=±10m/s, Δh0=±10m, and Δγ0=±3° randomly. The Monte-Carlo simulations of 100 times are carried out. The dispersions of the terminal states are illustrated in Figure 9, and the statistical results of the terminal states in terms of the mean values and standard deviations are exhibited in Table 5.Figure 9
Monte-Carlo simulation results for different tracking guidance schemes. Dispersion of (a) terminal velocity, and (b) terminal range, (b) terminal flight path angle.
(a)(b)(c)Table 5
Statistical results of Monte-Carlo simulations for tracking guidance.
Statistical resultsOpen-loopLQRmeanxfm8602.458499.90stdxfm327.340.26meanυfm/s340.69340.70stdυfm/s0.530.43meanγf°−59.54−59.97stdγf°1.590.00According to the Monte-Carlo simulation results, it can be observed that the impact points for the open-loop scheme deviate from the reference trajectory with the mean value 8602.45 m and standard deviation 327.34 m. The LQR controller for tracking guidance has a distinct advantage over the open-loop scheme in terms of the impact accuracy and constrained angle under uncertainties. To conclude, the effectiveness and robustness of the gained schemed LQR controller for tracking guidance are verified under uncertainties in the numerical simulations.


## 6. Conclusions

In this paper, an offline hybrid trajectory optimization approach is proposed for variable-sweep missiles to explore the superiority in the diving phase. The trajectory optimization model is firstly established under multiple constraints, and the aerodynamic characteristics in different sweep angles are analyzed, which indicate that the lift coefficient and drag coefficient decrease when the sweep angle increases. In addition, different from the traditional attack angle for optimization process, the two-variable optimization scheme on both the attack angle and sweep angle is presented. Then, a hybrid optimization strategy combining the GPM and DSM is presented to obtain the high precision solution, and the gain-scheduled LQR controller is designed for the tracking guidance. Numerical simulations show the effectiveness of the proposed hybrid optimization strategy and the superiority of variable-sweep missiles in contrast with traditional missiles.In the future study, the main work about the variable-sweep missiles is to design the stable and reliable control system in order to track the attack angle and sweep angle generated by the LQR controller for guidance tracking.



---

*Source: 6669618-2021-03-31.xml*