# Fast Trajectory Generation Method for Midcourse Guidance Based on Convex Optimization

Authors: Jinlin Zhang; Jiong Li; Chijun Zhou; Humin Lei; Wanli Li
Journal: International Journal of Aerospace Engineering
 (2022)

---

## Abstract

To generate the midcourse guidance trajectory for intercepting the high-speed and high maneuvering target, which is a strongly nonlinear and strongly constrained problem, a two-stage convex optimization method is proposed to solve the optimal trajectory quickly. In the first stage, an initial trajectory generation method is proposed, by which the trajectory’s terminal state is close to the terminal position. And the generated trajectory is used as the initial solution of the convex optimization method. In the second stage, the original nonconvex optimization problem is transformed into a convex optimization problem by linearization and relaxation methods and then solved discretely. In the numerical simulation, the effectiveness of the proposed method is verified, and the robustness of the method is verified in different initial and terminal states. Then, several ablation experiments are operated to verify the advantage of the rapid initial trajectory generation method in the first stage. Finally, compared with the gauss pseudospectral method (GPM), the proposed method is proved to be efficient and has the potential of online trajectory generation.


---

## Body

## 1. Introduction

More and more countries are developing delivery vehicles and weapons to dominate aerospace [1]. The midcourse guidance of interceptor plays an important role in intercepting such targets and can provide a good intercepting situation for terminal guidance with path constraints [2].There are two types of design methods for midcourse guidance. The first is to determine that the guidance law is designed to control the line of sight, when the relative dynamic equation of the interceptor and the target, by which the adverse influence of the target maneuver could be overcome [3]. However, in the midcourse guidance, the accurate information of the high-speed target is difficult, which will lead to the larger calculation error of the sight line [4]. To solve this problem, Wan et al. [5] used the model predictive static programming method to give a suboptimal guidance law with terminal angle constraints. Ebrahimi et al. [6] proposed a new generalized model predictive extended control technology, which smoothed the control signal. But during the intercept, there may be no guarantee that the interceptor meets some path restrictions. Ann et al. [7] combined zero effort trajectory with minimum time trajectory to predict intercept points. Liu et al. [8] proposed a midcourse guidance law considering the influence of random interference, detection range constraint, and target capture probability. However, the guidance accuracy of the first method is easily affected by target maneuver. The ability to intercept high speed and high maneuvering target at high altitude is limited [9].The second is generating the guidance instructions based on the terminal information prediction with several path constrains and terminal constraints [10]. The pseudospectral method has high accuracy in solving the strongly nonlinear and strongly constrained optimal control problems and has been rapidly developed and applied in the aircraft trajectory optimization [11]. Li et al. [12] proposed a multi-interval mesh refinement Radau pseudospectral method, which improved the optimization accuracy and convergence speed of the algorithm. Zhou et al. [13] improved the particle swarm optimization algorithm to solve the trajectory planning problem in the glide phase of hypersonic aircraft. Aiming at target maneuver and interceptor terminal constraints correction, Zhou et al. [14] and Li et al. [15] established a trajectory optimization and correction problem of midcourse guidance with path constraints. Du et al. [16] presented a semianalytical method for solving the exoatmospheric midcourse guidance problem with minimum velocity increment. Li et al. [17] designed an angle-constrained midcourse guidance trajectory according to the offline optimized trajectory information. But the trajectory generation efficiency of the above method is not ideal for intercepting 3d maneuverable high-speed target.In recent years, convex optimization methods have been widely used in aerospace aircraft due to the existence of solutions [18] and high computational efficiency of solving complex polynomials [19]. Liu et al. [20] and Wang and Grant [21] used linearization and relaxation methods to solve the entry trajectory optimization problem with multiple constraints based on second-order cone programming (SOCP). Yan et al. [22] presented the entry problem with multiple constraints as an easily solvable SOCP sequence. Wang and Grant [23] presented an autonomous entry aircraft guidance algorithm based on convex optimization and continuous closed loop method. Wang et al. [24] used convex optimization algorithms to provide high-quality initial guesses for pseudospectral methods. Hong et al. [25] presented an autonomous entry aircraft guidance algorithm based on convex optimization method and continuous closed loop. Liu et al. [26, 27] introduced a nonconvex constraint relaxation technique, which uses regularization technique to constrain the degree of relaxation and proves the effectiveness of the method. In order to improve the convergence rate of convex optimization, Wang and Lu [28] and Zhou et al. [29] proposed an adaptive method to adjust the shrinkage coefficient of trust region. Combining pseudospectral method and convex optimization method, Sagliano and Mooij [30] proposed a new drag-energy scheme for atmospheric entry guidance, which can los-less convexify the formulation. Zhou et al. [31] used an adaptive mesh refinement method to improve the computational efficiency of convex optimization method. The above methods show that convex optimization can solve complex polynomial constraint problems effectively.In midcourse guidance, Yang et al. [32] proposed a midcourse guidance method for boost interception based on an adaptive update strategy of trust region in convex optimization. In multi-interceptor guidance coordination, Jiang et al. [33] used convex optimization technology to solve the multiconstraint optimal proportional guidance problem for each interceptor online to achieve time consistency between interceptors. In this paper, the convex optimization method is applied to the generation of midcourse guidance trajectory for intercepting high-speed and high maneuvering target.In this paper, the trajectory optimization problem with strong nonlinearity and strong constraints is described, and the dynamics model of the interceptor is transformed by constructing affine variables. Then, the free terminal time is affined to the fixed time domain by the time expansion factor. The dynamic and path constraints are linearized when the control variables are relaxed, through which the problem are transformed into a convex optimization problem. Finally, the problem is discretized and the fast trajectory generation method is proposed, and the process of solving the problem is given.The innovations of this paper are as follows:(1)
The midcourse guidance problem with strong nonlinearity and strong constraints is transformed into a SOCP problem that can be solved iteratively by constructing an affine system and using relaxation and linearization methods. The effectiveness of the method is verified by simulation and the robustness of the method is verified under the condition of initial state perturbation and terminal position change(2)
A fast trajectory generation method is proposed as the initial solution of the convex optimization method. The advantages of the fast trajectory generation method are proved by ablation experiments. Compared with the GPM, the proposed method shows its rapidityThe remainder of this paper are introduced as follows. In Section2, the trajectory optimization model of midcourse guidance is constructed. In Section 3, the original problem is transformed into a convex optimization problem by using linearization and relaxation methods. In Section 4, the convex optimization problem is discretized, and the fast trajectory generation method and the process of solving the problem are given. In Section 5, the proposed two-stage convex optimization method is verified.


## 2. Midcourse Guidance Problem

### 2.1. Dynamics Model

For intercepting the high-speed and high maneuvering target, trajectory optimization of midcourse guidance is a trajectory optimization problem with strong nonlinearity and strong constraints.The model assumes the following(1)
As the earth’s rotation has little influence on the mid-course guidance flight of interceptor, the earth’s rotation is ignored(2)
The interceptor guidance and control errors are not consideredThe interceptor uses unpowered glide in the midcourse guidance phase, and the dynamics model is based on the ground coordinate system. Thex-axis and z-axis point north and east, respectively, and form the right-hand coordinate system with h-axis. The dynamics of the interceptor are modeled as follows:
(1)ḣ=Vsinθ,ż=Vcosθsinψ,ẋ=Vcosθcosψ,V̇=−D−sinθr2,θ̇=LcosσV+V2−1/rcosθVr,ψ̇=LsinσVcosθ,where h,z,x is the position coordinates; r=1+h is the straight-line distance from the center of the earth to the interceptor; V is the velocity of the interceptor relative to the Earth; θ is the flight path angle; ψ is the heading angle of the relative velocity vector; h, z and x are all scaled by re, the radius of the Earth; g0 is the acceleration of gravity at re; and V is scaled by g0re.Dimensionless lift and drag acceleration are as follows:(2)L=0.5ρVg0re2CLSmg0,(3)D=0.5ρVg0re2CDSmg0,whereS is the reference area of the interceptor; CL and CD are the lift coefficient and drag coefficient of the interceptor, respectively, which are related to the angle of attack and Mach number; m is the mass of the interceptor; and ρ is the atmosphere density expressed as
(4)ρh=ρ0e−hre/H,where ρ0=1.225kg/m3 and H=7254.3m.


### 2.2. Control-Affine System

Dynamic equation (1) is directly controlled by the angle of attack and bank angle, and there is no affine relationship between them, so it cannot be directly convex. A new affine system is constructed to relax the control variables [26].First, define a new variable normalized lift coefficientη(5)η=CLCL∗,(6)CD=CD∗1+η22,where CL∗ and CD∗ are the lift coefficient corresponding to the maximum lift-to-drag ratio at a certain Mach, which can be obtained in the case of aerodynamic data sheets.According to equations (5) and (6), equations (2) and (3) can be transformed into
(7)L=L^η,D=D^1+η22,where L^=0.5ρV2reCL∗S/m,D^=0.5ρV2reCD∗S/mSo far, the dynamics equation is a function ofη, the dynamic equation is nonconvex because of η2. To solve this problem, define three new control variables
(8)u1=ηcosσ,u2=ηsinσ,u3=η2,(9)u12+u22−u3=0,where the normalized lift coefficientη is related to the aerodynamic characteristics of the interceptor. Assuming that the lift coefficient in this paper is nonnegative, the upper limit of η is η¯, the value range of A is as follows:
(10)0≤η≤η¯.It can be obtained through Equation (8)
(11)0≤u3≤u¯3,where u¯3=η¯2.According to the defined affine variables, equation (1) can be transformed into an affine system dynamic model
(12)dhdt=Vsinθ,dzdt=Vcosθsinψ,dxdt=Vcosθcosψ,dVdt=−0.51+u3D^−sinθr2,dθdt=u1L^V+V2−1/rcosθVr,dψdt=u2L^Vcosθ.For the convenience of expression, equation (12) is rewritten as the following nonlinear system:
(13)ẋ=fx,t+Βx,tu,wherex=h,z,x,V,θ,ψT;u=u1,u2,u3T and
(14)fx,t=VsinθVcosθsinψVcosθcosψ−0.5D^−sinθr2V2−1/rcosθVr0,(15)Bx,t=03×103×103×100−0.5D^L^V000L^Vcosθ0.


### 2.3. Path Constraints

Interceptors must meet path constraints during interception to ensure normal operation. Interceptor path constraints usually include heat rateQ̇, dynamic pressure q, and load n, which must not exceed the maximum value given by the interceptor.The three path constraints are expressed as follows:(16)Q̇=kQρ0.5Vreg03.15≤Q̇max,(17)q=0.5ρVreg02≤qmax,(18)n=L2+D2=u3L∧2+1+u32D∧24≤nmax.


### 2.4. Boundary Constraints and Objective Function

The starting point of interceptor midcourse guidance is the starting point of trajectory optimization. Initial point constraint is mandatory equality constraint(19)xt0−x0=0,wherex0=h0,z0,x0,V0,θ0,ψ0T is the initial value of the trajectory states.To intercept high-speed target, the interceptor needs a certain speed to ensure the successful destruction of the target in the kinetic energy interception. The direction of the interceptor’s velocity can be enforced as an equality constraint.(20)Vtf≥V¯,θtf=θf∗,ψtf=ψf∗,whereV¯ is the lower limit of interceptor speed and θf∗ and ψf∗ are the expected values of the flight path angle and heading angle at the interceptor terminal moment tf, respectively.The terminal position deviation can be constrained in the objective function. Position constraint relaxation can ensure the existence of viable solutions to the original problem. If the position constraint is directly used as the equality constraint, the original problem may not have a feasible solution, resulting in optimization failure. The result of the objective function is to find the optimal trajectory whose terminal position is closest to the expected value.(21)J0=κ1hf−hf∗+zf−zf∗+κ2xf−xf∗,where hf∗,zf∗,xf∗ is the expected value of the terminal position and κ1 and κ2 indicate the weight coefficient, which is determined by the value and priority of the status in the performance indicator.In summary, the original midcourse guidance trajectory optimization problem is transformed into a continuous time optimal control problem, represented as problem P0 and is represented as follows:(22)P0:minJ0subjectto:9,11,13,16,17,18,19,20.


## 2.1. Dynamics Model

For intercepting the high-speed and high maneuvering target, trajectory optimization of midcourse guidance is a trajectory optimization problem with strong nonlinearity and strong constraints.The model assumes the following(1)
As the earth’s rotation has little influence on the mid-course guidance flight of interceptor, the earth’s rotation is ignored(2)
The interceptor guidance and control errors are not consideredThe interceptor uses unpowered glide in the midcourse guidance phase, and the dynamics model is based on the ground coordinate system. Thex-axis and z-axis point north and east, respectively, and form the right-hand coordinate system with h-axis. The dynamics of the interceptor are modeled as follows:
(1)ḣ=Vsinθ,ż=Vcosθsinψ,ẋ=Vcosθcosψ,V̇=−D−sinθr2,θ̇=LcosσV+V2−1/rcosθVr,ψ̇=LsinσVcosθ,where h,z,x is the position coordinates; r=1+h is the straight-line distance from the center of the earth to the interceptor; V is the velocity of the interceptor relative to the Earth; θ is the flight path angle; ψ is the heading angle of the relative velocity vector; h, z and x are all scaled by re, the radius of the Earth; g0 is the acceleration of gravity at re; and V is scaled by g0re.Dimensionless lift and drag acceleration are as follows:(2)L=0.5ρVg0re2CLSmg0,(3)D=0.5ρVg0re2CDSmg0,whereS is the reference area of the interceptor; CL and CD are the lift coefficient and drag coefficient of the interceptor, respectively, which are related to the angle of attack and Mach number; m is the mass of the interceptor; and ρ is the atmosphere density expressed as
(4)ρh=ρ0e−hre/H,where ρ0=1.225kg/m3 and H=7254.3m.


## 2.2. Control-Affine System

Dynamic equation (1) is directly controlled by the angle of attack and bank angle, and there is no affine relationship between them, so it cannot be directly convex. A new affine system is constructed to relax the control variables [26].First, define a new variable normalized lift coefficientη(5)η=CLCL∗,(6)CD=CD∗1+η22,where CL∗ and CD∗ are the lift coefficient corresponding to the maximum lift-to-drag ratio at a certain Mach, which can be obtained in the case of aerodynamic data sheets.According to equations (5) and (6), equations (2) and (3) can be transformed into
(7)L=L^η,D=D^1+η22,where L^=0.5ρV2reCL∗S/m,D^=0.5ρV2reCD∗S/mSo far, the dynamics equation is a function ofη, the dynamic equation is nonconvex because of η2. To solve this problem, define three new control variables
(8)u1=ηcosσ,u2=ηsinσ,u3=η2,(9)u12+u22−u3=0,where the normalized lift coefficientη is related to the aerodynamic characteristics of the interceptor. Assuming that the lift coefficient in this paper is nonnegative, the upper limit of η is η¯, the value range of A is as follows:
(10)0≤η≤η¯.It can be obtained through Equation (8)
(11)0≤u3≤u¯3,where u¯3=η¯2.According to the defined affine variables, equation (1) can be transformed into an affine system dynamic model
(12)dhdt=Vsinθ,dzdt=Vcosθsinψ,dxdt=Vcosθcosψ,dVdt=−0.51+u3D^−sinθr2,dθdt=u1L^V+V2−1/rcosθVr,dψdt=u2L^Vcosθ.For the convenience of expression, equation (12) is rewritten as the following nonlinear system:
(13)ẋ=fx,t+Βx,tu,wherex=h,z,x,V,θ,ψT;u=u1,u2,u3T and
(14)fx,t=VsinθVcosθsinψVcosθcosψ−0.5D^−sinθr2V2−1/rcosθVr0,(15)Bx,t=03×103×103×100−0.5D^L^V000L^Vcosθ0.


## 2.3. Path Constraints

Interceptors must meet path constraints during interception to ensure normal operation. Interceptor path constraints usually include heat rateQ̇, dynamic pressure q, and load n, which must not exceed the maximum value given by the interceptor.The three path constraints are expressed as follows:(16)Q̇=kQρ0.5Vreg03.15≤Q̇max,(17)q=0.5ρVreg02≤qmax,(18)n=L2+D2=u3L∧2+1+u32D∧24≤nmax.


## 2.4. Boundary Constraints and Objective Function

The starting point of interceptor midcourse guidance is the starting point of trajectory optimization. Initial point constraint is mandatory equality constraint(19)xt0−x0=0,wherex0=h0,z0,x0,V0,θ0,ψ0T is the initial value of the trajectory states.To intercept high-speed target, the interceptor needs a certain speed to ensure the successful destruction of the target in the kinetic energy interception. The direction of the interceptor’s velocity can be enforced as an equality constraint.(20)Vtf≥V¯,θtf=θf∗,ψtf=ψf∗,whereV¯ is the lower limit of interceptor speed and θf∗ and ψf∗ are the expected values of the flight path angle and heading angle at the interceptor terminal moment tf, respectively.The terminal position deviation can be constrained in the objective function. Position constraint relaxation can ensure the existence of viable solutions to the original problem. If the position constraint is directly used as the equality constraint, the original problem may not have a feasible solution, resulting in optimization failure. The result of the objective function is to find the optimal trajectory whose terminal position is closest to the expected value.(21)J0=κ1hf−hf∗+zf−zf∗+κ2xf−xf∗,where hf∗,zf∗,xf∗ is the expected value of the terminal position and κ1 and κ2 indicate the weight coefficient, which is determined by the value and priority of the status in the performance indicator.In summary, the original midcourse guidance trajectory optimization problem is transformed into a continuous time optimal control problem, represented as problem P0 and is represented as follows:(22)P0:minJ0subjectto:9,11,13,16,17,18,19,20.


## 3. Problem of Convexification

Problem P0 is a nonconvex problem. Next, it is transformed into a convex optimization problem. Firstly, the terminal time free problem dynamic model is affined to the fixed time domain. Then, the dynamic and path constraints are linearized, and the control variable constraints are relaxed. Finally, the convex optimization problem is summarized.


### 3.1. Linearization of Dynamic Constraints

Timet is in the variable time domain t0,tf in the problem. t0 and tf represent the initial and terminal time of the interceptor in problem P0, respectively. tf is uncertain. It can be transformed into the time domain τ∈0,1 by affine τ=t−t0/tf−t0.Applying this time domainτ to equation (13), the following can be obtained:
(23)ẋ=λfxk,τ+Βxk,τuτ,where λ is the time expansion factor, which scales the time domain t0,tf.The reference trajectoryxk,uk,λk is linearized for equation (23), and the result is as follows”
(24)ẋ=λkAxk,τxτ+λkBxk,τuτ+λFxk,uk,τ+λkCk,where superscript k represents the kth reference trajectory, and
(25)Fxk,uk,τ=fxk,τ+Βxk,τuk,(26)Ck=−Axk+Buk,(27)Axk,τ=000a14a150000a24a25a26000a34a35a36a4100a44a450a5100a54a550000000x=xk,(28)a14=sinθ,a15=Vcosθ,(29)a24=cosθsinψ,a25=−Vsinθsinψ,a26=Vcosθcosψ,(30)a34=cosθcosψ,a35=−Vsinθcosψ,a36=−Vcosθsinψ,(31)a41=reD^2H+2sinθr3,a44=−D^V,a45=−cosθr2,(32)a51=2cosθVr3−Vcosθr2,a54=cosθr+cosθVr2,a55=sinθVr2−Vsinθr.The additional constraints of the trust region are expressed by the following inequality:(33)x−xk≤εx,u−uk≤εu,where εx∈R6 and εu∈R6 are trust region constraint radius of xk and uk, respectively. They are constants large enough to limit the range of deviations from the reference trajectories defined for different systems.Equations (19) and (20) are transformed into the time domain τ, as shown below:
(34)xτ=0−x0=0,(35)Vτ=1≥V¯,θτ=1=θf,ψτ=1=ψf.


### 3.2. Linearization of Path Constraints

The inequality constraints of Equations (16)–(18) are functions of h, V, andu3, which linearizes the path constraints at fixed state hk,Vk,u3k; the form is as follows:
(36)fih,V,u3≈fihk,Vk,u3k+fi′hk,Vk,u3kh−hk;V−Vk;u3−u3k,≤fi,max,i=1,2,3,where fi,max represents the maximum value of the ith constraint, and
(37)f1′hk,Vk,u3k=δf1δh,δf1δV,δf1δu3=−kQreρ0.5Vreg03.15H,3.15kQρ0.5reg03.15V2.15,0hk,Vk,uk,(38)f2′hk,Vk,u3k=δf2δh,δf2δV,δf2δu3=−ρreVreg022H,ρreg02V,0hk,Vk,uk,(39)f3′hk,Vk,u3k=δf3δh,δf3δV,δf3δu3=−L2+D2reH,2L2+D2V,L2+D2C^D2+2C^L2+C^D2u3C^D2+2C^D2+4C^L2u3+C^D2u32hk,Vk,uk.When intercepting the high-speed and high maneuvering target, the interceptor usually intercepts the target in the near space with high speed. Because of the material limitations of the interceptor, it may result in a lower limit of the interceptor’s flight altitude. The lower limit of altitude is mainly determined by the speed of the interceptor, the flight altitude, and the attitude of the aircraft. Assume that the lower limit of the height of an interceptor ishmin, then
(40)h≥hmin.


### 3.3. Relaxation of Control Variable Constraints

It can be seen from Equations (8), (9), and (11) that the control variable is a strong equality constraint. The control domain formed by this equation causes the problem of non-convexity. To convexify problem P0, equations (9) and (11) can be relaxed as
(41)u∈convU=u12+u22≤u3,0≤u3≤u¯3,where convU is the convex domain containing u12+u22−u3=0.In order to ensure the effectiveness of the control variable after relaxation, the optimal solution of the control variable satisfies equation (9). Add an additional term κ3∫01ψτdτ to the objective function J0. The small value κ3 is to ensure that the original objective function J0 is not affected.


### 3.4. Convex Optimization Problem

A continuous time optimal control problem P0 can be transformed into a convex optimization problem after linearization of dynamic and path constraints and convexity of control variable constraints.(42)P1:minJ0+κ3∫01ψτdτsubjectto:24,33,34,35,36,40,41,where equation (39) is the objective function. Equation (24) is the dynamic constraint. Equation (33) is the trust region constraint. Equations (34) and (35) are the boundary condition constraint. Equations (36) and (40) are the path constraints. Equation (41) is the control variable constraints.


## 3.1. Linearization of Dynamic Constraints

Timet is in the variable time domain t0,tf in the problem. t0 and tf represent the initial and terminal time of the interceptor in problem P0, respectively. tf is uncertain. It can be transformed into the time domain τ∈0,1 by affine τ=t−t0/tf−t0.Applying this time domainτ to equation (13), the following can be obtained:
(23)ẋ=λfxk,τ+Βxk,τuτ,where λ is the time expansion factor, which scales the time domain t0,tf.The reference trajectoryxk,uk,λk is linearized for equation (23), and the result is as follows”
(24)ẋ=λkAxk,τxτ+λkBxk,τuτ+λFxk,uk,τ+λkCk,where superscript k represents the kth reference trajectory, and
(25)Fxk,uk,τ=fxk,τ+Βxk,τuk,(26)Ck=−Axk+Buk,(27)Axk,τ=000a14a150000a24a25a26000a34a35a36a4100a44a450a5100a54a550000000x=xk,(28)a14=sinθ,a15=Vcosθ,(29)a24=cosθsinψ,a25=−Vsinθsinψ,a26=Vcosθcosψ,(30)a34=cosθcosψ,a35=−Vsinθcosψ,a36=−Vcosθsinψ,(31)a41=reD^2H+2sinθr3,a44=−D^V,a45=−cosθr2,(32)a51=2cosθVr3−Vcosθr2,a54=cosθr+cosθVr2,a55=sinθVr2−Vsinθr.The additional constraints of the trust region are expressed by the following inequality:(33)x−xk≤εx,u−uk≤εu,where εx∈R6 and εu∈R6 are trust region constraint radius of xk and uk, respectively. They are constants large enough to limit the range of deviations from the reference trajectories defined for different systems.Equations (19) and (20) are transformed into the time domain τ, as shown below:
(34)xτ=0−x0=0,(35)Vτ=1≥V¯,θτ=1=θf,ψτ=1=ψf.


## 3.2. Linearization of Path Constraints

The inequality constraints of Equations (16)–(18) are functions of h, V, andu3, which linearizes the path constraints at fixed state hk,Vk,u3k; the form is as follows:
(36)fih,V,u3≈fihk,Vk,u3k+fi′hk,Vk,u3kh−hk;V−Vk;u3−u3k,≤fi,max,i=1,2,3,where fi,max represents the maximum value of the ith constraint, and
(37)f1′hk,Vk,u3k=δf1δh,δf1δV,δf1δu3=−kQreρ0.5Vreg03.15H,3.15kQρ0.5reg03.15V2.15,0hk,Vk,uk,(38)f2′hk,Vk,u3k=δf2δh,δf2δV,δf2δu3=−ρreVreg022H,ρreg02V,0hk,Vk,uk,(39)f3′hk,Vk,u3k=δf3δh,δf3δV,δf3δu3=−L2+D2reH,2L2+D2V,L2+D2C^D2+2C^L2+C^D2u3C^D2+2C^D2+4C^L2u3+C^D2u32hk,Vk,uk.When intercepting the high-speed and high maneuvering target, the interceptor usually intercepts the target in the near space with high speed. Because of the material limitations of the interceptor, it may result in a lower limit of the interceptor’s flight altitude. The lower limit of altitude is mainly determined by the speed of the interceptor, the flight altitude, and the attitude of the aircraft. Assume that the lower limit of the height of an interceptor ishmin, then
(40)h≥hmin.


## 3.3. Relaxation of Control Variable Constraints

It can be seen from Equations (8), (9), and (11) that the control variable is a strong equality constraint. The control domain formed by this equation causes the problem of non-convexity. To convexify problem P0, equations (9) and (11) can be relaxed as
(41)u∈convU=u12+u22≤u3,0≤u3≤u¯3,where convU is the convex domain containing u12+u22−u3=0.In order to ensure the effectiveness of the control variable after relaxation, the optimal solution of the control variable satisfies equation (9). Add an additional term κ3∫01ψτdτ to the objective function J0. The small value κ3 is to ensure that the original objective function J0 is not affected.


## 3.4. Convex Optimization Problem

A continuous time optimal control problem P0 can be transformed into a convex optimization problem after linearization of dynamic and path constraints and convexity of control variable constraints.(42)P1:minJ0+κ3∫01ψτdτsubjectto:24,33,34,35,36,40,41,where equation (39) is the objective function. Equation (24) is the dynamic constraint. Equation (33) is the trust region constraint. Equations (34) and (35) are the boundary condition constraint. Equations (36) and (40) are the path constraints. Equation (41) is the control variable constraints.


## 4. Problem Solving Process

In this section, convex optimization problem P1 is discretized and transformed into SOCP problem P2. Then, a fast initial trajectory generation method is proposed as the initial reference solution for convex optimization problems. Finally, the process of solving SOCP problem is given.


### 4.1. Discretization

Problem P1 is a continuous parameter optimization problem in convex domain, which cannot be directly optimized. Therefore, problem P1 should be discretized first to facilitate computer numerical processing.In this section, variables are discretized in the time domainτ∈0,1 and all constraints are enforced at discrete points. Assuming that the time domain τ∈0,1 is divided into N discrete intervals Δτ=1/N, the total number of discrete points is N+1, and the discrete points can be expressed as τj,j=0,1,2,⋯,N−1,N.The dynamic constraint equation (24) is discretized by the trapezoidal method:
(43)xj=xj−1+λkAj−1kxj−1+λkBj−1kuj−1+λFj−1k+λkCj−1k+λkAjkxj+λkBjkuj+λFjk+λkCjkΔτ2,wherexj=xτj, uj=uτj, Ajk=Axk,τj, Bjk=Bxk,τj, Fjk=Fxk,uk,τj, Cjk=Cxk,uk,τj, j is the jth discrete point, j=0,1,2,⋯,N−1,N.In order to reduce the optimization space and improve the convergence speed, the variable trust region method is adopted in this paper. The method is to introduce the relaxation coefficientξ into the trust domain constraint and extend it into the objective function.
(44)xj−xjk≤ξεx,uj−ujk≤ξεu.After discretization of problem P1, the relaxation coefficient is extended to the objective function. Then, by updating the path constraints and control variable constraints, the problem P1 is transformed into the SOCP problem P2, as follows:(45)P2:minJ0+κ3∑j=0NψjΔτ+κ4ξsubjectto:34,35,43,44,fihjk,Vjk,u3jk+fi′hjk,Vjk,u3jk⋅hj−hjk;Vj−Vjk;u3j−u3jk≤fi,max,i=1,2,3,hj≥hmin,u1j2+u2j2≤u3j,0≤u3j≤u¯3.


### 4.2. Fast Trajectory Generation Method

The deviation between the initial trajectory and the optimal trajectory determines the number of optimization iterations and the efficiency of solving convex optimization problems to a certain extent. In this paper, initial trajectories are rapidly generated by interpolation and integration. The steps for initial trajectory generation are as follows.The maximum normalized lift coefficientη¯ is selected to ensure that the initial trajectory generated quickly satisfies the trajectory constraint to the maximum extent. Within the bank angle constraint, the value of equal interval is used as the control variables.According to the dynamic model, the control variables at different bank angles are integrated, respectively. Generate a trajectory group.From the trajectory group, two trajectoriesl1 and l2 closest to the desired position of the terminal are selected. The distances from the desired position of the terminal are Δl1 and Δl2. For the bank angles σ1 and σ2 corresponding to the two trajectories, the bank angle σ3 is selected according to
(46)σ3=σ2Δl1+σ1Δl2Δl1+Δl2.The final control variables are calculated according to the bank angleσ3. Repeat step two to generate the initial trajectory. Finally, the maximum normalized lift coefficient η¯ and the bank angle σ3 are transformed into the control variables u1, u2, and u3.Under the initial conditions in Table1 and different terminal conditions in Table 2, according to the above fast trajectory generation method, the trajectory is generated as follows.Table 1
Fast trajectories generate initial conditions.
h0kmz0kmx0kmV0m/sθ0°ψ0°40003000100Table 2
Fast trajectories generate different terminal conditions.
xf∗kmzf1∗kmzf2∗kmzf3∗km75002040The flight times of interceptors at three different terminal positionszf1∗, zf2∗, and zf3∗ are tf1=325s, tf2=326s, and tf3=329s, respectively. The trajectories are shown in Figure 1. The maximum normalized lift coefficient η¯ is taken as the control parameter for all three trajectories, and the corresponding bank angles are shown in Figure 2. The control variables remain unchanged in the flight trajectory. The larger zf∗ is, the larger the flight time is and the larger the corresponding bank angle is.Figure 1
Three initial trajectories of different terminals.⊳ represents the starting position of the trajectory. ○ represents the expected end position of the trajectory. The meaning is the same in the figure below.
(a)
Longitudinal trajectory(b)
Lateral trajectoryFigure 2
Three initial trajectory bank angles of different terminals.Compared with the GPM, the rapid trajectory generation method saves the iterative solution time and greatly improves the initial trajectory generation speed. The method takes about 0.1032 s.The purpose of the fast trajectory generation method is to make the generated trajectory close to the optimal trajectory, reduce the number of convex optimization iterations, and improve the efficiency of problem solving.


### 4.3. SOCP Method Solving Process

In this section, problem P2 is solved iteratively until the discrete points converge during the solution process. This solution is the highly approximate solution of the original problem P0. The following is the solution flow of problem P2.Step 1: setk=0. Select N and the initial time domain t0,tf. The initial time expansion factor λ0 is obtained. The discrete points τj,j=0,1,2,⋯,N−1,N are selected on the initial trajectory x0=h0,z0,x0,V0,θ0,ψ0T and the control variable u0=u10,u20,u30T of the initial trajectory.Step 2: at thekth iteration, the result of the k−1th iteration xk−1,uk−1,λk−1 is used to establish the problem P2. Trajectory state variables, control variables, and time expansion factors xk,uk,λk are obtained by solving the SOCP problem.Step 3: verify whether the iterative results meet the following convergence conditions:(47)xk−xk−1≤δx,uk−uk−1≤δu.If all states and control variables of all discrete points meet the convergence conditions in equation (47), step into Step 4. Otherwise, set k=k+1 and go back to Step 2.Step 4: iteration stops. The convergent solutionxk,uk,λk is the highly approximate solution of the problem P0.


## 4.1. Discretization

Problem P1 is a continuous parameter optimization problem in convex domain, which cannot be directly optimized. Therefore, problem P1 should be discretized first to facilitate computer numerical processing.In this section, variables are discretized in the time domainτ∈0,1 and all constraints are enforced at discrete points. Assuming that the time domain τ∈0,1 is divided into N discrete intervals Δτ=1/N, the total number of discrete points is N+1, and the discrete points can be expressed as τj,j=0,1,2,⋯,N−1,N.The dynamic constraint equation (24) is discretized by the trapezoidal method:
(43)xj=xj−1+λkAj−1kxj−1+λkBj−1kuj−1+λFj−1k+λkCj−1k+λkAjkxj+λkBjkuj+λFjk+λkCjkΔτ2,wherexj=xτj, uj=uτj, Ajk=Axk,τj, Bjk=Bxk,τj, Fjk=Fxk,uk,τj, Cjk=Cxk,uk,τj, j is the jth discrete point, j=0,1,2,⋯,N−1,N.In order to reduce the optimization space and improve the convergence speed, the variable trust region method is adopted in this paper. The method is to introduce the relaxation coefficientξ into the trust domain constraint and extend it into the objective function.
(44)xj−xjk≤ξεx,uj−ujk≤ξεu.After discretization of problem P1, the relaxation coefficient is extended to the objective function. Then, by updating the path constraints and control variable constraints, the problem P1 is transformed into the SOCP problem P2, as follows:(45)P2:minJ0+κ3∑j=0NψjΔτ+κ4ξsubjectto:34,35,43,44,fihjk,Vjk,u3jk+fi′hjk,Vjk,u3jk⋅hj−hjk;Vj−Vjk;u3j−u3jk≤fi,max,i=1,2,3,hj≥hmin,u1j2+u2j2≤u3j,0≤u3j≤u¯3.


## 4.2. Fast Trajectory Generation Method

The deviation between the initial trajectory and the optimal trajectory determines the number of optimization iterations and the efficiency of solving convex optimization problems to a certain extent. In this paper, initial trajectories are rapidly generated by interpolation and integration. The steps for initial trajectory generation are as follows.The maximum normalized lift coefficientη¯ is selected to ensure that the initial trajectory generated quickly satisfies the trajectory constraint to the maximum extent. Within the bank angle constraint, the value of equal interval is used as the control variables.According to the dynamic model, the control variables at different bank angles are integrated, respectively. Generate a trajectory group.From the trajectory group, two trajectoriesl1 and l2 closest to the desired position of the terminal are selected. The distances from the desired position of the terminal are Δl1 and Δl2. For the bank angles σ1 and σ2 corresponding to the two trajectories, the bank angle σ3 is selected according to
(46)σ3=σ2Δl1+σ1Δl2Δl1+Δl2.The final control variables are calculated according to the bank angleσ3. Repeat step two to generate the initial trajectory. Finally, the maximum normalized lift coefficient η¯ and the bank angle σ3 are transformed into the control variables u1, u2, and u3.Under the initial conditions in Table1 and different terminal conditions in Table 2, according to the above fast trajectory generation method, the trajectory is generated as follows.Table 1
Fast trajectories generate initial conditions.
h0kmz0kmx0kmV0m/sθ0°ψ0°40003000100Table 2
Fast trajectories generate different terminal conditions.
xf∗kmzf1∗kmzf2∗kmzf3∗km75002040The flight times of interceptors at three different terminal positionszf1∗, zf2∗, and zf3∗ are tf1=325s, tf2=326s, and tf3=329s, respectively. The trajectories are shown in Figure 1. The maximum normalized lift coefficient η¯ is taken as the control parameter for all three trajectories, and the corresponding bank angles are shown in Figure 2. The control variables remain unchanged in the flight trajectory. The larger zf∗ is, the larger the flight time is and the larger the corresponding bank angle is.Figure 1
Three initial trajectories of different terminals.⊳ represents the starting position of the trajectory. ○ represents the expected end position of the trajectory. The meaning is the same in the figure below.
(a)
Longitudinal trajectory(b)
Lateral trajectoryFigure 2
Three initial trajectory bank angles of different terminals.Compared with the GPM, the rapid trajectory generation method saves the iterative solution time and greatly improves the initial trajectory generation speed. The method takes about 0.1032 s.The purpose of the fast trajectory generation method is to make the generated trajectory close to the optimal trajectory, reduce the number of convex optimization iterations, and improve the efficiency of problem solving.


## 4.3. SOCP Method Solving Process

In this section, problem P2 is solved iteratively until the discrete points converge during the solution process. This solution is the highly approximate solution of the original problem P0. The following is the solution flow of problem P2.Step 1: setk=0. Select N and the initial time domain t0,tf. The initial time expansion factor λ0 is obtained. The discrete points τj,j=0,1,2,⋯,N−1,N are selected on the initial trajectory x0=h0,z0,x0,V0,θ0,ψ0T and the control variable u0=u10,u20,u30T of the initial trajectory.Step 2: at thekth iteration, the result of the k−1th iteration xk−1,uk−1,λk−1 is used to establish the problem P2. Trajectory state variables, control variables, and time expansion factors xk,uk,λk are obtained by solving the SOCP problem.Step 3: verify whether the iterative results meet the following convergence conditions:(47)xk−xk−1≤δx,uk−uk−1≤δu.If all states and control variables of all discrete points meet the convergence conditions in equation (47), step into Step 4. Otherwise, set k=k+1 and go back to Step 2.Step 4: iteration stops. The convergent solutionxk,uk,λk is the highly approximate solution of the problem P0.


## 5. Numerical Simulations

All the calculations in this paper are performed on a laptop equipped with Intel Core I7-10510 2.30 GHz, 8 G RAM, and Windows 10 operating system. ECOS [34] is used to solve SOCP problem P2 on MATLAB, where the number of discrete points was 70.The interceptor model parameters arem0=900kg and S=0.4839m2. For the three path constraints, the maximum allowable values are ˙Q̇max=3800kw/m2, qmax=370kPa, and nmax=3.5g0. The value range of η is 0,5.5869, and the value range of bank angle σ is −60ο,60ο. The weight coefficients in the objective function are κ1=100, κ2=10, κ3=0.01 and κ4=0.01, respectively.The terminal constraints are shown in Table3.Table 3
Terminal constraints.
hf∗kmzf3∗kmxf∗kmV¯fm/sθf∗°ψf∗°30407501800100When solving problem P2, the trust region constraint parameters are as follows:(48)εx=5e3re,5e3re,5e4re,500g0re,10π180,10π180,εu=η¯2,η¯2,η¯22,whereη¯ is given in Equation (10).The convergence conditions of simulation are given in(49)δx=100re,100re,500re,10g0re,0.5π180,0.5π180,(50)δu=η¯10,η¯10,η¯210.In Section5.1, the feasibility of the two-stage convex optimization method is verified. In Section 5.2, the robustness of the method is verified under the initial angle perturbation condition. In Section 5.3, the robustness of the proposed method is verified by selecting different terminal position states. In Section 5.4, different initial trajectory generation methods are selected to prove the superiority of the fast trajectory generation method. In Section 5.5, compared with the GPM, the efficient solving ability of the two-stage convex optimization method is verified.


### 5.1. Feasibility Verification

In order to reduce the optimization space of the flight time of interceptors, the range of flight time should be constrained. The lower boundtmin and upper bound tmax of time are generally determined empirically. When the time optimization range is too small, the feasible solution may not be found. However, if the interval is too large, the trajectory may not converge quickly.In Figure3, the convergence of the error between the trajectory terminal position and the desired position are shown. The error converges in the terminal position during the first iteration. This is caused by position error which has great influence on objective function and variable trust region method. The fluctuation of position error is caused by the fluctuation of relaxation coefficient ξ and regularization term ψ. The purpose of subsequent iterations is to make the trajectory converge so as to ensure the validity of the iterative solution.Figure 3
The error between the trajectory terminal position and the terminal position constraint.In Figure4, the position error for each iteration is shown. After each iteration, the trajectory gradually converges. The results show that the third and fourth iterations are convergent.Figure 4
Trajectories during iteration.
(a)
Longitudinal trajectories during iteration(b)
Lateral trajectories during iterationIn Figure5, path constraints for convergent trajectories are shown. Heating rate, dynamic pressure, and overload all meet constraints requirements. There are peaks between 150 s and 200 s. This is because the trajectory height of convergence presents a concave curve, and the air density curve presents a convex curve, which leads to the path constraint presenting a convex curve.Figure 5
Path constraints for convergent solutions.Figure6 shows the control variables of the convergence trajectory. The values of u1 and u2 can be transformed by u3 and bank angle σ, both of which satisfy the constraints.Figure 6
Control variables for convergent solutions.The flight time and operation time of four iterations are shown in Table4. In each iteration, ECOS consumed approximately 0.4072 s of CPU to solve a SOCP problem. The total computing time is the computing time of four iterations plus the generation time of initial trajectory, which is about 1.732 seconds.Table 4
The flight time and operation time of each iteration trajectory.
Iterative sequenceFlight timeOperation time1291.40.46402287.70.42103287.40.36424287.30.3797


### 5.2. Initial Angle Perturbation

In this section, the robustness of the proposed method is analyzed by using initial angle perturbation, which provides a reference for the initial angle design of the interceptor. The contrastive trajectory is the convergence trajectory of Section4.1. The initial angle variation condition (1) is the initial flight path angle θ0=11ο. Variation condition (2) is heading angle ψ0=1ο.In Figure7, the convergence trajectories under initial angular perturbation are shown. When the initial flight path angle increases, the initial velocity component of the interceptor on the h-axis increases, resulting in an increase in the height peak of the trajectory. At the same time, the initial velocity component of the interceptor along the z-axis decreases, resulting in a more concave trajectory pattern. When the initial heading angle increases, the initial velocity component of the interceptor is opposite to that when the initial heading angle increases.Figure 7
Trajectories under initial angle perturbation.
(a)
Longitudinal trajectories under initial angle perturbation(b)
Lateral trajectories under initial angle perturbationIn Figures8 and 9, the control variables of the convergent trajectories satisfy the constraints. Flight path angle and heading angle reach the terminal angle constraints from different angles. To some extent, the robustness of the proposed optimization method is proved. In order to obtain the cost time in solving the midcourse guidance problem, 100 Monte Carlo simulations were performed with an average cost time of 1.730 s.Figure 8
Comparison of flight path angle and heading angleFigure 9
Comparison of control variables


### 5.3. Terminal Position Change

In the trajectory generation stage, interceptors need to generate optimized trajectory groups to deal with the attacks of target at different positions. This section selects three representative terminal positionszf1∗=0km,zf2∗=20km, and zf3∗=40km. For other boundary condition constraints, refer to Tables 1 and 3.Figure10 shows the convergence trajectories under constraints of different terminal positions. The trajectory lines are concave in Figure 10(b). This is due to a change in heading angle from 0° to 10° and the increasing projection of the interceptor’s velocity on the z-axis.Figure 10
Trajectories at different terminal positions.
(a)
Longitudinal trajectories at different terminal positions(b)
Lateral trajectories at different terminal positionsFigure11 shows the variation of control variables under the constraints of different terminal positions. When the terminal position is zf1∗=0km, the flight path angle is from 0° to negative and then to 10°, because the bank angle is first negative and then positive. This results in a negative and then positive projection of the interceptor’s velocity on the z-axis. The convergence trajectory goes from 0 km to negative and then to 0 km on the z-axis.Figure 11
Comparison of path constraints.Control variables of convergent trajectories under different terminal positions all meet the constraint conditions in Figure12.Figure 12
Comparison of control variables.In this section, when the control variables meet the constraints, the terminal positions of the trajectories reach the terminal position constraints. To some extent, the robustness of the proposed optimization method is proved.


### 5.4. Different Initial Trajectories

To demonstrate the advantages of the fast trajectory generation method in this paper, it is compared with the method in reference [31]. Except terminal position constraint zf2∗=20km, other boundary condition constraints refer to Table 1 and Table 3. Numerical simulation is as follows.It can be seen from Table5 and Figure 13 that the optimization iteration under the initial trajectory one condition reaches the convergence condition after 4 times, and it takes 5 times under the initial trajectory two condition, consuming an extra time of about 0.2974 seconds. The reason is that the initial trajectory two deviates more from the convergence trajectory relative to the initial trajectory one, and more iterations are needed to approach the convergence trajectory.Table 5
Comparison of the results of different initial trajectories.
Position errorFlight timeTotal solution timeInitial trajectory one5.8e−4279.21.6456Initial trajectory two1.6e−3278.11.9430Figure 13
Iterative trajectories of different initial reference trajectories. Initial trajectory one is generated by the fast trajectory generation method in this paper, and initial trajectory two is generated by the method in reference [31].
(a)
Longitudinal trajectories below the initial trajectory one(b)
Longitudinal trajectories below the initial trajectory two(c)
Lateral trajectories below the initial trajectory one(d)
Lateral trajectories below the initial trajectory twoThe control variables of convergent trajectories under two different initial trajectory conditions, and the control variables all satisfy constraints in Figure14.Figure 14
Control variables under different initial trajectories.The generation method of initial trajectory should pay attention to two aspects:(a)
The rate at which the initial trajectory is generated. If the generation time is long, it is contrary to the high efficiency of convex optimization method(b)
The deviation between the generated initial trajectory and the convergent trajectory. If the deviation degree is large, the solving process needs more iterations and takes more operation timeA good initial trajectory generation method should satisfy the need to generate an initial trajectory close to convergence in a short time.


### 5.5. Different Optimization Methods

In order to demonstrate the high efficiency of the proposed method, the GPOPS software was used to solve the GPM under the condition of the same calculation accuracy of the objective function. The comparison results are as follows.From Table6, both the two-order convex optimization method and the GPOPS can obtain the midcourse guidance trajectory satisfying the constraints, and the GPM has higher accuracy. Compared with the GPOPS, the convex optimization method has shorter computing time. The reason is that CVX adopts “warm-start” initial guesses and is insensitive to initial values. Compared with the CVX, GPOPS adopts interpolation “cold-start” initial guesses, which is more sensitive to initial value and with a lower computational efficiency [24].Table 6
Comparison of results between two methods.
Position errorFlight timeIterations numberTotal solution timeCVX5.8e−5279.241.732GPOPS1.2e−7282.01536.497It can be seen from Figure15 that the optimization trajectory generated by the two-stage convex optimization method and GPOPS are not exactly the same. This is mainly caused by the differences in the ways of internal matching points and initial trajectory generation.Figure 15
Trajectories under different optimization methods.
(a)
Longitudinal trajectories under different optimization methods(b)
Lateral trajectories under different optimization methodsFigure16 shows the changes of control variables under the two methods. All the control variables meet the constraints. The control variable of the convex optimization method fluctuates greatly, which may be caused by the scarcity of matching points and fewer iterations.Figure 16
Control variables under different optimization methods.In conclusion, for midcourse guidance trajectory optimization, the accuracy of the proposed method is slightly lower than that of GPM. But the running speed of the method is greatly improved, which can provide help for rapid trajectory generation and real-time trajectory correction.


## 5.1. Feasibility Verification

In order to reduce the optimization space of the flight time of interceptors, the range of flight time should be constrained. The lower boundtmin and upper bound tmax of time are generally determined empirically. When the time optimization range is too small, the feasible solution may not be found. However, if the interval is too large, the trajectory may not converge quickly.In Figure3, the convergence of the error between the trajectory terminal position and the desired position are shown. The error converges in the terminal position during the first iteration. This is caused by position error which has great influence on objective function and variable trust region method. The fluctuation of position error is caused by the fluctuation of relaxation coefficient ξ and regularization term ψ. The purpose of subsequent iterations is to make the trajectory converge so as to ensure the validity of the iterative solution.Figure 3
The error between the trajectory terminal position and the terminal position constraint.In Figure4, the position error for each iteration is shown. After each iteration, the trajectory gradually converges. The results show that the third and fourth iterations are convergent.Figure 4
Trajectories during iteration.
(a)
Longitudinal trajectories during iteration(b)
Lateral trajectories during iterationIn Figure5, path constraints for convergent trajectories are shown. Heating rate, dynamic pressure, and overload all meet constraints requirements. There are peaks between 150 s and 200 s. This is because the trajectory height of convergence presents a concave curve, and the air density curve presents a convex curve, which leads to the path constraint presenting a convex curve.Figure 5
Path constraints for convergent solutions.Figure6 shows the control variables of the convergence trajectory. The values of u1 and u2 can be transformed by u3 and bank angle σ, both of which satisfy the constraints.Figure 6
Control variables for convergent solutions.The flight time and operation time of four iterations are shown in Table4. In each iteration, ECOS consumed approximately 0.4072 s of CPU to solve a SOCP problem. The total computing time is the computing time of four iterations plus the generation time of initial trajectory, which is about 1.732 seconds.Table 4
The flight time and operation time of each iteration trajectory.
Iterative sequenceFlight timeOperation time1291.40.46402287.70.42103287.40.36424287.30.3797


## 5.2. Initial Angle Perturbation

In this section, the robustness of the proposed method is analyzed by using initial angle perturbation, which provides a reference for the initial angle design of the interceptor. The contrastive trajectory is the convergence trajectory of Section4.1. The initial angle variation condition (1) is the initial flight path angle θ0=11ο. Variation condition (2) is heading angle ψ0=1ο.In Figure7, the convergence trajectories under initial angular perturbation are shown. When the initial flight path angle increases, the initial velocity component of the interceptor on the h-axis increases, resulting in an increase in the height peak of the trajectory. At the same time, the initial velocity component of the interceptor along the z-axis decreases, resulting in a more concave trajectory pattern. When the initial heading angle increases, the initial velocity component of the interceptor is opposite to that when the initial heading angle increases.Figure 7
Trajectories under initial angle perturbation.
(a)
Longitudinal trajectories under initial angle perturbation(b)
Lateral trajectories under initial angle perturbationIn Figures8 and 9, the control variables of the convergent trajectories satisfy the constraints. Flight path angle and heading angle reach the terminal angle constraints from different angles. To some extent, the robustness of the proposed optimization method is proved. In order to obtain the cost time in solving the midcourse guidance problem, 100 Monte Carlo simulations were performed with an average cost time of 1.730 s.Figure 8
Comparison of flight path angle and heading angleFigure 9
Comparison of control variables


## 5.3. Terminal Position Change

In the trajectory generation stage, interceptors need to generate optimized trajectory groups to deal with the attacks of target at different positions. This section selects three representative terminal positionszf1∗=0km,zf2∗=20km, and zf3∗=40km. For other boundary condition constraints, refer to Tables 1 and 3.Figure10 shows the convergence trajectories under constraints of different terminal positions. The trajectory lines are concave in Figure 10(b). This is due to a change in heading angle from 0° to 10° and the increasing projection of the interceptor’s velocity on the z-axis.Figure 10
Trajectories at different terminal positions.
(a)
Longitudinal trajectories at different terminal positions(b)
Lateral trajectories at different terminal positionsFigure11 shows the variation of control variables under the constraints of different terminal positions. When the terminal position is zf1∗=0km, the flight path angle is from 0° to negative and then to 10°, because the bank angle is first negative and then positive. This results in a negative and then positive projection of the interceptor’s velocity on the z-axis. The convergence trajectory goes from 0 km to negative and then to 0 km on the z-axis.Figure 11
Comparison of path constraints.Control variables of convergent trajectories under different terminal positions all meet the constraint conditions in Figure12.Figure 12
Comparison of control variables.In this section, when the control variables meet the constraints, the terminal positions of the trajectories reach the terminal position constraints. To some extent, the robustness of the proposed optimization method is proved.


## 5.4. Different Initial Trajectories

To demonstrate the advantages of the fast trajectory generation method in this paper, it is compared with the method in reference [31]. Except terminal position constraint zf2∗=20km, other boundary condition constraints refer to Table 1 and Table 3. Numerical simulation is as follows.It can be seen from Table5 and Figure 13 that the optimization iteration under the initial trajectory one condition reaches the convergence condition after 4 times, and it takes 5 times under the initial trajectory two condition, consuming an extra time of about 0.2974 seconds. The reason is that the initial trajectory two deviates more from the convergence trajectory relative to the initial trajectory one, and more iterations are needed to approach the convergence trajectory.Table 5
Comparison of the results of different initial trajectories.
Position errorFlight timeTotal solution timeInitial trajectory one5.8e−4279.21.6456Initial trajectory two1.6e−3278.11.9430Figure 13
Iterative trajectories of different initial reference trajectories. Initial trajectory one is generated by the fast trajectory generation method in this paper, and initial trajectory two is generated by the method in reference [31].
(a)
Longitudinal trajectories below the initial trajectory one(b)
Longitudinal trajectories below the initial trajectory two(c)
Lateral trajectories below the initial trajectory one(d)
Lateral trajectories below the initial trajectory twoThe control variables of convergent trajectories under two different initial trajectory conditions, and the control variables all satisfy constraints in Figure14.Figure 14
Control variables under different initial trajectories.The generation method of initial trajectory should pay attention to two aspects:(a)
The rate at which the initial trajectory is generated. If the generation time is long, it is contrary to the high efficiency of convex optimization method(b)
The deviation between the generated initial trajectory and the convergent trajectory. If the deviation degree is large, the solving process needs more iterations and takes more operation timeA good initial trajectory generation method should satisfy the need to generate an initial trajectory close to convergence in a short time.


## 5.5. Different Optimization Methods

In order to demonstrate the high efficiency of the proposed method, the GPOPS software was used to solve the GPM under the condition of the same calculation accuracy of the objective function. The comparison results are as follows.From Table6, both the two-order convex optimization method and the GPOPS can obtain the midcourse guidance trajectory satisfying the constraints, and the GPM has higher accuracy. Compared with the GPOPS, the convex optimization method has shorter computing time. The reason is that CVX adopts “warm-start” initial guesses and is insensitive to initial values. Compared with the CVX, GPOPS adopts interpolation “cold-start” initial guesses, which is more sensitive to initial value and with a lower computational efficiency [24].Table 6
Comparison of results between two methods.
Position errorFlight timeIterations numberTotal solution timeCVX5.8e−5279.241.732GPOPS1.2e−7282.01536.497It can be seen from Figure15 that the optimization trajectory generated by the two-stage convex optimization method and GPOPS are not exactly the same. This is mainly caused by the differences in the ways of internal matching points and initial trajectory generation.Figure 15
Trajectories under different optimization methods.
(a)
Longitudinal trajectories under different optimization methods(b)
Lateral trajectories under different optimization methodsFigure16 shows the changes of control variables under the two methods. All the control variables meet the constraints. The control variable of the convex optimization method fluctuates greatly, which may be caused by the scarcity of matching points and fewer iterations.Figure 16
Control variables under different optimization methods.In conclusion, for midcourse guidance trajectory optimization, the accuracy of the proposed method is slightly lower than that of GPM. But the running speed of the method is greatly improved, which can provide help for rapid trajectory generation and real-time trajectory correction.


## 6. Conclusions

In this paper, a two-stage convex optimization method for trajectory generation of high-speed target interception with strong nonlinearity and strong constraints is proposed. By constructing affine system to linearize, relax, and discretize the original nonconvex problem, a trajectory convex optimization problem that is easier to solve is established, which is convenient for engineering application. A fast initial trajectory generation method is designed, which not only ensures the feasibility of convergent solution but also improves the convergence efficiency and solving speed of convex optimization method.The influence of initial angle and terminal position constraints on the optimization trajectory is analyzed, which provides guidance for trajectory offline design. Compared with the GPM, the proposed method has a higher efficiency in solving the problem of trajectory generation with strong nonlinearity and strong constraints for midcourse guidance. It has potential for engineering applications and online trajectory generation.



---

*Source: 7188718-2022-06-25.xml*