# Trajectory Optimization for High-Speed and Long-Range Interceptor Based on Improved Adaptive hp Pseudospectral Method

Authors: Changxin Luo; Chijun Zhou; Jiong Li; Humin Lei
Journal: International Journal of Aerospace Engineering
 (2022)

---

## Abstract

The trajectory optimization design of interceptor is very important in the defense combat against a class of high-speed and strong-maneuvering aircraft. At present, an important idea is to use hp pseudospectral method to offline optimizes the trajectory of interceptor, but the solving efficiency of this method needs to be further improved. Aiming at this issue, an improved adaptive hp pseudospectral method is proposed in this paper. In order to shorten the solving time of the algorithm, the proposed method has two main improvements on the basis of the traditional hp pseudospectral method: on one hand, by judging the positions of the control sudden change points, prerefine the mesh around them according to certain rules. On the other hand, the curvature of the system state curve is used as the criterion to segment the original mesh nonuniformly so that more mesh points can be allocated where the curvature is large. These two points together ensure that the collocation point resources can be used more efficiently in the mesh refinement process. The simulation results show that the proposed method can solve the optimized trajectory of interceptor effectively, and it also proves that this method has higher solving efficiency than the traditional adaptive hp method.


---

## Body

## 1. Introduction

Recently, with the rapid development of aerospace science and technology, a class of high-speed and strong-maneuvering vehicle such as near space hypersonic vehicle and high-performance fighter aircraft has gradually entered people’s vision. Compared with traditional aircraft, they have obvious advantages such as faster flight speed and stronger penetration capability, which pose great challenges to the country’s air and space defense operations [1–4]. Therefore, it is urgent to carry out research on intercepting technology against such targets [5].The defense against such high-speed and strong-maneuvering targets belongs to the super-speed and super-range interception issue. In order to intercept such targets better, the range and speed of interceptor are the key breakthrough, and the separation design of trajectory and guidance law is the fundamental way to improve the range and speed of interceptor [6]. In view of this, a kind of high-speed and long-range interceptor should be selected to defend such incoming targets, and the trajectory of interceptor should be rationally optimized. As for the trajectory optimization of interceptor, many scholars have actively carried out research and achieved some results.The trajectory optimization problem is essentially a kind of optimal control problem. However, due to the complexity of the interceptor model and many constraints in the trajectory optimization process, it is difficult to get the analytical solution of the optimized trajectory directly. In order to solve this problem, pseudospectral method [7–14] is an important idea. The method firstly selects the corresponding collocation points (such as Legendre-Gauss (LG) points [7], Legendre-Gauss-Radau (LGR) points [11, 12], or Legendre-Gauss-Lobatto (LGL) points [10]) in the time domain. Then, based on the selected collocation points, Lagrange interpolation is used to transform the original optimal control problem into a nonlinear programming problem (NLP). Finally, NLP solver [15] is used to solve the problem. Unfortunately, for complex optimization problems such as trajectory optimization, a satisfactory solution may not be obtained by using a fixed number of collocation points. Therefore, by changing the number and position of collocation points according to certain rules, many scholars continue to use pseudospectral method to solve the NLP problem iteratively until the expected solution accuracy is obtained. Typical methods include p method [16] and h method [17–19]. p method is mainly used for iterative optimization by increasing the number of collocation points, and it has a better solution effect for optimization problems with smooth solutions [16]. However, if the solution of the problem to be solved is not smooth, h method is needed to obtain higher solving efficiency, which is mainly employed to refine the mesh by constantly dividing the time domain [17–19].Since h method and p method have their own advantages for different types of problems, in order to further improve the efficiency of optimal solution, researchers consider combining h method with p method and propose an adaptive mesh refinement method called hp method [20]-[21]. hp method adaptively adopts h method or p method in mesh refinement through a preset criterion, so that the advantages of both methods can be fully brought into play [20–22]. In recent years, many scholars [23–30] have further enriched and developed adaptive mesh refinement method along the research ideas of Ref. [20-22]. According to the posteriori error analysis, Ref. [23] proposes a new adaptive hp-refinement strategy to solve discrete local problems. Based on Birkhoff interpolation, Ref. [24] proposes a fast mesh refinement strategy to generate well-conditioned PS optimal control discretizations. For nonsmooth control problem, a mesh adaptation method is proposed in Ref. [25]. In addition, Ref. [26] introduces a C++ software program for solving the multiple-phase optimal control problems, which is called CGPOPS. Ref. [27] proposes a new adaptive mesh refinement method, which is based on the decay rate of Legendre polynomial coefficients to determine whether to use h method or p method for mesh refinement, and reduce the mesh size if necessary. In order to shorten the operation time of the algorithm, Ref. [28] refines the mesh with the arc length of the state estimation curve as the criterion, which improves the efficiency of the algorithm. By adopting the mesh density function to generate the grid, Ref. [29] proposes an adaptive control parameterization method to solve the optimal control problem. Ref. [30] designs a simple and effective adaptive mesh refinement algorithm based on interpolation error analysis and improved data compression technology, and verifies the effectiveness of the method in solving optimization problems through simulation.Recently, many different optimization methods have been applied to solve trajectory optimization problems, and there are many meaningful research results [31–38]. Based on a series of LGL collocation points, the trajectory planning problem of interceptor is transformed into NLP, and Ref. [31] combines direct transcription and mesh refinement, and then proposes a new adaptive hp method. Aiming at the trajectory optimization problem of hypersonic reentry vehicle, Refs. [32, 33] effectively solve the trajectory optimization problem based on a multiphase Gauss pseudospectral method and a new adaptive Gauss pseudospectral method, respectively. Refs. [34, 35], respectively, explore the trajectory optimization of solar-powered aircraft and Mars entry by using the adaptive mesh refinement pseudospectral method. Ref. [36] combines Gauss pseudospectral method with Finite Fourier Series and proposes a feasible close-range satellite cooperative trajectory planning scheme. For space manoeuvre vehicles skip trajectory optimization problem, a fuzzy physical programming method, and a new three-layer-hybrid optimal control solver are proposed in Refs. [37, 38], respectively.The above studies on aircraft trajectory optimization have achieved considerable results for their respective research concerns. Considering the specific defense operations background, the high-speed and long-range interceptor will have a period of time for offline nominal trajectory optimization before launching. Based on this, hp pseudospectral method is usually adopted in trajectory design of interceptor. In this paper, the existing hp pseudospectral method will be ameliorated to improve the optimization efficiency of interceptor trajectory.This paper mainly improves the mesh refinement process of hp pseudospectral method, so as to enhance the solving efficiency of the algorithm. Although many scholars have made some achievements in interceptor trajectory optimization based on hp pseudospectral method, further considerations are still needed in the following aspects:(1)
Most of the existing studies [31–38] do not take into account deeply the terminal constraints that the interceptor should meet at the midcourse guidance in the trajectory optimization design of interceptor, and the optimized trajectory cannot ensure that the target can be detected when the interceptor seeker is turned on(2)
When pseudospectral method is used to optimize interceptor trajectory, most studies [23–30] do not pay attention to the sudden changes of optimization control in the process of mesh refinement, which will cause a large number of collocation points to be configured near the control sudden change points, thus increasing the number of mesh iterations and affecting optimization efficiency(3)
In the process of trajectory optimization, the existing researches [23–30] basically insert new mesh points at equal intervals during the piecewise iteration of the mesh, which will cause the waste of collocation points resources and reduce the solving efficiencyIn view of the above considerations, and draw on the ideas of [39], this paper proposes an improved adaptive hp pseudospectral method to solve the trajectory optimization of interceptor based on the combat scenario of high-speed and long-range interceptor defending high-speed strong-maneuvering target. The main innovations are as follows:
(1)
The terminal constraint of interceptor in midcourse guidance are considered and derived. On this basis, the optimized trajectory can ensure that the seeker of interceptor can successfully detect the target and complete the handover between midcourse and terminal guidance smoothly(2)
By judging the position of the control sudden change points, the mesh is refined in advance around it. In this way, the collocation points resources are preferentially configured to reduce the number of mesh iterations(3)
The nonuniform mesh segmentation iteration is carried out based on the curvature of the system state curve, so that more mesh points can be allocated where the curvature is large, so as to further improve the efficiency of solving trajectory optimizationThe rest of the paper is organized as follows. In Section2, the trajectory scheme of interceptor is designed and the mathematical model of high-speed and long-range interceptor is established. Section 3 describes the trajectory optimization problem to be solved in this paper. In Section 4, the improved adaptive hp pseudospectral method proposed in this paper is described in detail. In Section 5, the effectiveness and superiority of the proposed algorithm in solving trajectory optimization problems are verified by simulation. Section 6 summarizes the whole work.


## 2. High-Speed and Long-Range Interceptor Model

### 2.1. Interceptor Trajectory Scheme for High-Speed Strong-Maneuvering Target

Because the velocity ratio between missile and target will directly affect the guidance accuracy of interceptor, in order to effectively defend against a class of high-speed and strong-maneuvering target (hereinafter referred to as “target”), the interceptor needs to have a high interception speed. At the same time, because the target’s great flight speed will further compress the interception window of the defense side, in order to push the interception window forward to increase the interception margin of the defense side, the interceptor needs to have a larger range. As we know, an effective way to improve the range and speed of the interceptor is to adopt a hill-turn-reentry trajectory, as shown in Figure1.Figure 1
High-speed and long-range interceptor trajectory scheme.In order to achieve this type of interceptor trajectory, the high-speed and long-range interceptor adopted in this paper consists of three parts: the first stage booster, the second stage booster, and the third stage interceptor system. The third stage interceptor system mainly includes the following: warhead, fairing, detection system, missile-borne computer, and engine which can switch on and off at any time. As shown in Figure1, at the end of the first and second boost phases, the interceptor separates from the first and second booster, respectively. The first and second stage rockets boost the interceptor rapidly through the lower dense atmosphere, which greatly reduces the velocity loss of the interceptor due to air resistance. At the same time, under the continuous boost of the two stage rocket, the third stage interceptor can be sent to the thin air in the upper atmosphere, and make it have a large flight speed, which is also an important guarantee for the high-speed and long-range interception.Because the target has the characteristics of fast flight speed and strong maneuverability, in order to ensure that the planned trajectory can guide the interceptor to the target, it is necessary to use the correlation prediction algorithm to predict and estimate the target trajectory and get the high probability region (HPR) of the target. At the termination of midcourse guidance, it is necessary to make the detection field of the seeker successfully cover the HPR of the target, that is, to successfully capture the target, so as to realize the handover between midcourse and terminal guidance.


### 2.2. Mathematical Model of Interceptor

In this paper, it is assumed that the seeker of the third stage interceptor is placed on its axis, that is to say, the detection direction of the seeker depends on the attitude of the interceptor. Therefore, in order to make the interceptor successfully complete the handover between midcourse and terminal guidance, the terminal attitude of the interceptor at the midcourse guidance needs to meet certain constraints. In view of this, in order to facilitate the subsequent analysis, the pitch angleϑ and yaw angle ψ (roll angle is assumed to be 0) and the switch engine uon‐off are taken as the control input. By considering the force applied to the interceptor and the definition of the relevant coordinate system [40], the interceptor model in the ground coordinate system is established as Equation. (1).where x,y,z is the position vector of the interceptor in the ground coordinate system, m represents the mass of the interceptor, g is the gravity constant, T is the engine thrust, and I represents the specific impulse of the engine. Vx, Vy, and Vz are components of velocity vector of interceptor in three axes of ground coordinate system, respectively.
(1)dxdt=VcosθcosψV,dydt=Vsinθ,dzdt=−VcosθsinψV,mdVxdt=Tcosϑcosψuon‐off−DcosθcosψV+LsinγVsinψV−sinθcosγVcosψV+ZsinθsinγVcosψV+cosγVsinψV,mdVydt=Tsinϑuon‐off−Dsinθ+LcosθcosγV−ZcosθsinγV−mg,mdVzdt=−Tcosϑsinψuon‐off+DcosθsinψV+LsinγVcosψV+sinθcosγVinψV+ZcosγVcosψV−sinθsinγVsinψV,dmdt=−TgIuon‐off.Velocity valueV of interceptor can be expressed as V=Vx2+Vy2+Vz2. Flight path angle θ and heading angle ψV can be calculated by Vx, Vy, and Vz.
(2)θ=arcsinVy/V,(3)ψV=−arctanVz/Vx.According to the geometric relation equation of angle, angle of sideslipβ, angle of attack α and bank angle γV can be calculated from pitch angle ϑ, yaw angle ψ, flight path angle θ, and heading angle ψV.
(4)β=arcsincosθsinψ−ψV,(5)α=arccoscosϑcosθcosψ−ψV+sinϑsinθcosβ,(6)γv=arccoscosψ−ψV/cosβ.Based on the solid rocket model in [41] and the background of defense operations, this paper makes reasonable improvements, and designs the relevant performance parameters and aerodynamic parameters of the interceptor as Table 1.Table 1
Overall parameters of interceptor.
First stage boosterSecond stage boosterThird stage interceptorWhole weight/kg122854008900Empty weight/kg41951105400Engine run time/s67.449.3100Specific impulseI/s227.8272.5510.2Engine thrustT/N26793115726125000Reference areaS/m20.70.70.7During the flight of interceptor, drag forceD, lift force L, and side force Z are expressed as
(7)D=ρV2CDS2,L=ρV2CLS2,Z=ρV2CZS2,ρ=ρ0e−y/H.where ρ is the air density, ρ0=1.226kg/m3 is the air density at sea level, H=7254.24m represents the reference height, and y represents the height of the interceptor from sea level. CD, CL, and CZ are drag force coefficient, lift force coefficient and side force coefficient, respectively.
(8)CD=0.29,0≤Ma≤0.8,Ma−0.51,0.8≤Ma≤1.07,0.091+0.5Ma−1,Ma>1.07,(9)CLα=2.8.0≤Ma≤0.25,2.8+0.447Ma−0.25.0.25≤Ma≤1.1,3.18−0.660Ma−1.1.1.1≤Ma≤1.6,2.85+0.35Ma−1.6.1.6≤Ma≤3.6,3.55,Ma>3.6,(10)CZβ=−CLα,(11)CL=CLαα,(12)CZ=CZββ,where Ma stands for Mach number, CLα represents the partial derivatives of CL with respect to α, and CZβ represents the partial derivatives of CZ with respect to β.In the midcourse guidance stage, the aerodynamic parameters of the third stage interceptor are(13)CD=0.1763+0.5655α+3.882α2,(14)CL=4.115α−1.96α2,(15)CZ=−4.115β+1.96β2,In addition,uon‐off in Equation (1) represents the switching amount of the engine. During the first and second boost stages, the rocket engine is continuously turned on, so there is uon‐off≡1. In the midcourse guidance stage, the engine of the third stage interceptor can be switched on and off at any time, and when it is turned on, there is uon‐off=1. When the engine is turned off, there is uon‐off=0.Sinceuon‐off is defined as discrete value (equal to 0 or 1) in the midcourse guidance stage, while it is necessary to ensure that the input control is continuous value in the subsequent trajectory optimization process. Therefore, pulse-width pulse-frequency modulator (PWPF) [42] is introduced here to modulate the continuous control ucont to be optimized into discrete switching instruction uon‐off. Its schematic diagram is shown in Figure 2.Figure 2
The principle diagram of the PWPF.where the gain of the first-order plant isKm=1 and the time constant is Tm=1. The turn-on threshold for relays is Uon=0.2 and the turn-off threshold is Uoff=0.05. ucont is the switching control to be optimized, which is continuous value. ucont is always 1 in the first and second boost stages, and there is 0≤ucont≤1 in the midcourse guidance stage.


## 2.1. Interceptor Trajectory Scheme for High-Speed Strong-Maneuvering Target

Because the velocity ratio between missile and target will directly affect the guidance accuracy of interceptor, in order to effectively defend against a class of high-speed and strong-maneuvering target (hereinafter referred to as “target”), the interceptor needs to have a high interception speed. At the same time, because the target’s great flight speed will further compress the interception window of the defense side, in order to push the interception window forward to increase the interception margin of the defense side, the interceptor needs to have a larger range. As we know, an effective way to improve the range and speed of the interceptor is to adopt a hill-turn-reentry trajectory, as shown in Figure1.Figure 1
High-speed and long-range interceptor trajectory scheme.In order to achieve this type of interceptor trajectory, the high-speed and long-range interceptor adopted in this paper consists of three parts: the first stage booster, the second stage booster, and the third stage interceptor system. The third stage interceptor system mainly includes the following: warhead, fairing, detection system, missile-borne computer, and engine which can switch on and off at any time. As shown in Figure1, at the end of the first and second boost phases, the interceptor separates from the first and second booster, respectively. The first and second stage rockets boost the interceptor rapidly through the lower dense atmosphere, which greatly reduces the velocity loss of the interceptor due to air resistance. At the same time, under the continuous boost of the two stage rocket, the third stage interceptor can be sent to the thin air in the upper atmosphere, and make it have a large flight speed, which is also an important guarantee for the high-speed and long-range interception.Because the target has the characteristics of fast flight speed and strong maneuverability, in order to ensure that the planned trajectory can guide the interceptor to the target, it is necessary to use the correlation prediction algorithm to predict and estimate the target trajectory and get the high probability region (HPR) of the target. At the termination of midcourse guidance, it is necessary to make the detection field of the seeker successfully cover the HPR of the target, that is, to successfully capture the target, so as to realize the handover between midcourse and terminal guidance.


## 2.2. Mathematical Model of Interceptor

In this paper, it is assumed that the seeker of the third stage interceptor is placed on its axis, that is to say, the detection direction of the seeker depends on the attitude of the interceptor. Therefore, in order to make the interceptor successfully complete the handover between midcourse and terminal guidance, the terminal attitude of the interceptor at the midcourse guidance needs to meet certain constraints. In view of this, in order to facilitate the subsequent analysis, the pitch angleϑ and yaw angle ψ (roll angle is assumed to be 0) and the switch engine uon‐off are taken as the control input. By considering the force applied to the interceptor and the definition of the relevant coordinate system [40], the interceptor model in the ground coordinate system is established as Equation. (1).where x,y,z is the position vector of the interceptor in the ground coordinate system, m represents the mass of the interceptor, g is the gravity constant, T is the engine thrust, and I represents the specific impulse of the engine. Vx, Vy, and Vz are components of velocity vector of interceptor in three axes of ground coordinate system, respectively.
(1)dxdt=VcosθcosψV,dydt=Vsinθ,dzdt=−VcosθsinψV,mdVxdt=Tcosϑcosψuon‐off−DcosθcosψV+LsinγVsinψV−sinθcosγVcosψV+ZsinθsinγVcosψV+cosγVsinψV,mdVydt=Tsinϑuon‐off−Dsinθ+LcosθcosγV−ZcosθsinγV−mg,mdVzdt=−Tcosϑsinψuon‐off+DcosθsinψV+LsinγVcosψV+sinθcosγVinψV+ZcosγVcosψV−sinθsinγVsinψV,dmdt=−TgIuon‐off.Velocity valueV of interceptor can be expressed as V=Vx2+Vy2+Vz2. Flight path angle θ and heading angle ψV can be calculated by Vx, Vy, and Vz.
(2)θ=arcsinVy/V,(3)ψV=−arctanVz/Vx.According to the geometric relation equation of angle, angle of sideslipβ, angle of attack α and bank angle γV can be calculated from pitch angle ϑ, yaw angle ψ, flight path angle θ, and heading angle ψV.
(4)β=arcsincosθsinψ−ψV,(5)α=arccoscosϑcosθcosψ−ψV+sinϑsinθcosβ,(6)γv=arccoscosψ−ψV/cosβ.Based on the solid rocket model in [41] and the background of defense operations, this paper makes reasonable improvements, and designs the relevant performance parameters and aerodynamic parameters of the interceptor as Table 1.Table 1
Overall parameters of interceptor.
First stage boosterSecond stage boosterThird stage interceptorWhole weight/kg122854008900Empty weight/kg41951105400Engine run time/s67.449.3100Specific impulseI/s227.8272.5510.2Engine thrustT/N26793115726125000Reference areaS/m20.70.70.7During the flight of interceptor, drag forceD, lift force L, and side force Z are expressed as
(7)D=ρV2CDS2,L=ρV2CLS2,Z=ρV2CZS2,ρ=ρ0e−y/H.where ρ is the air density, ρ0=1.226kg/m3 is the air density at sea level, H=7254.24m represents the reference height, and y represents the height of the interceptor from sea level. CD, CL, and CZ are drag force coefficient, lift force coefficient and side force coefficient, respectively.
(8)CD=0.29,0≤Ma≤0.8,Ma−0.51,0.8≤Ma≤1.07,0.091+0.5Ma−1,Ma>1.07,(9)CLα=2.8.0≤Ma≤0.25,2.8+0.447Ma−0.25.0.25≤Ma≤1.1,3.18−0.660Ma−1.1.1.1≤Ma≤1.6,2.85+0.35Ma−1.6.1.6≤Ma≤3.6,3.55,Ma>3.6,(10)CZβ=−CLα,(11)CL=CLαα,(12)CZ=CZββ,where Ma stands for Mach number, CLα represents the partial derivatives of CL with respect to α, and CZβ represents the partial derivatives of CZ with respect to β.In the midcourse guidance stage, the aerodynamic parameters of the third stage interceptor are(13)CD=0.1763+0.5655α+3.882α2,(14)CL=4.115α−1.96α2,(15)CZ=−4.115β+1.96β2,In addition,uon‐off in Equation (1) represents the switching amount of the engine. During the first and second boost stages, the rocket engine is continuously turned on, so there is uon‐off≡1. In the midcourse guidance stage, the engine of the third stage interceptor can be switched on and off at any time, and when it is turned on, there is uon‐off=1. When the engine is turned off, there is uon‐off=0.Sinceuon‐off is defined as discrete value (equal to 0 or 1) in the midcourse guidance stage, while it is necessary to ensure that the input control is continuous value in the subsequent trajectory optimization process. Therefore, pulse-width pulse-frequency modulator (PWPF) [42] is introduced here to modulate the continuous control ucont to be optimized into discrete switching instruction uon‐off. Its schematic diagram is shown in Figure 2.Figure 2
The principle diagram of the PWPF.where the gain of the first-order plant isKm=1 and the time constant is Tm=1. The turn-on threshold for relays is Uon=0.2 and the turn-off threshold is Uoff=0.05. ucont is the switching control to be optimized, which is continuous value. ucont is always 1 in the first and second boost stages, and there is 0≤ucont≤1 in the midcourse guidance stage.


## 3. Trajectory Optimization Problem Description

Combined with the above model, the first and second boost stages are programmed guidance stage, and the principle of trajectory optimization in this stage is to send the interceptor to a certain altitude and make it have the maximum flight speed. While for midcourse guidance, the principle of trajectory optimization is to make the third stage interceptor smoothly complete the handover between midcourse and terminal guidance. At the same time, the fuel consumption of interceptor engine in the midcourse guidance should be reduced as much as possible to ensure that the third stage interceptor can have a large interception speed and strong maneuverability in the terminal guidance stage. According to the above principles, it is necessary to solve the trajectory optimization in these two stages, respectively. Firstly, the trajectory optimization problem is described as follows.(16)minJ=Exf,tfs.t.ẋt=fxt,utexf,tf=0hxt,ut≤0,where J=Exf,tf is the cost function, ẋt=fxt,ut is the state equation of the system, exf,tf=0 represents the terminal constraint, and hxt,ut≤0 is the process constraint. xt=x,y,z,Vx,Vy,Vz,m is the system state variable, ut=ϑ,ψ,ucont is the control variable to be optimized, and tf and xf represent terminal time and terminal state, respectively.


### 3.1. Trajectory Constraints in the Programmed Guidance Stage

For programmed guidance stage (first and second boost stages), the cost function is to maximize the velocity at the terminal of the second boost stage, i.e.,(17)J=−Vf.The state equation is shown in Equation (1). The terminal constraint is that the flight path angle and heading angle of the interceptor reach the expected value at the terminal of the second boost stage, i.e.,
(18)θf−θref=0,ψVf−ψref=0,where θf and ψVf represent the flight path angle and heading angle at the terminal moment, respectively, θref and ψref represent the desired flight path angle and heading angle, respectively.In addition, during the programmed guidance stage, the interceptor needs to meet the process constraints of dynamic pressureq, heat flux density Q, overload n, i.e.,
(19)0≤q≤qmax,0≤Q≤Qmax,0≤n≤nmax,where qmax, Qmax, and nmax represent upper limits of q, Q and n, respectively, and
(20)q=12ρV2,Q=Cρ0.5V3.05,n=T+D+L+Zmg,where C=1.291×10−4 is the heat flux constant and T+D+L+Z is the vector sum of thrust, drag, lift, and side forces.


### 3.2. Trajectory Constraints in the Midcourse Guidance Stage

For the midcourse guidance stage, the cost function is to minimize the fuel consumption of the third stage interceptor during flight, that is, to maximize the remaining mass of the third stage interceptor at the end of the midcourse guidance, i.e.(21)J=−mf.The state equation of the interceptor in this flight phase is shown in Equation (1). In order to achieve the handover between midcourse and terminal guidance, it is necessary to ensure that the third stage interceptor can detect the target at the end of the midcourse guidance stage, that is, to ensure that the detection range of the third stage interceptor’s seeker can successfully cover the HPR of the target, which is displayed in Figure 3.Figure 3
The handover between midcourse and terminal guidance.As shown in Figure3, it is assumed that the HPR of the target is a spherical region with radius r in space, the coordinate of its spherical center Tp is xp,yp,zp, the position of the interceptor’s centroid M is x,y,z, and the distance between M and Tp is d. In order to ensure that the seeker’s detection field can cover the target’s HPR more effectively, it is necessary to make the axis direction of the seeker’s detection field consistent with the line direction of MTp, that is, to make the pitch angle ϑ and yaw angle ψ of the interceptor meet the following geometric relations.
(22)ϑ=arcsinyp−yd,ψ′=arcsinz−zpd2−yp−y2,ψ=ψ′,xp≥x,sgnz−zp·π−ψ′,xp<x.Furthermore, assuming that the detection field angle of the seeker isϕ and the maximum detection range is Rmax, the following constraints should be met to ensure that the detection field can completely cover the HPR of the target.
(23)d≥rtanϕ,d≤Rmax−r.Equations (22) and (23) are the terminal constraints of interceptor in the midcourse guidance stage. In addition, the third interceptor also needs to meet various process constraints in the midcourse guidance stage, as shown in Equations (19) and (20).


## 3.1. Trajectory Constraints in the Programmed Guidance Stage

For programmed guidance stage (first and second boost stages), the cost function is to maximize the velocity at the terminal of the second boost stage, i.e.,(17)J=−Vf.The state equation is shown in Equation (1). The terminal constraint is that the flight path angle and heading angle of the interceptor reach the expected value at the terminal of the second boost stage, i.e.,
(18)θf−θref=0,ψVf−ψref=0,where θf and ψVf represent the flight path angle and heading angle at the terminal moment, respectively, θref and ψref represent the desired flight path angle and heading angle, respectively.In addition, during the programmed guidance stage, the interceptor needs to meet the process constraints of dynamic pressureq, heat flux density Q, overload n, i.e.,
(19)0≤q≤qmax,0≤Q≤Qmax,0≤n≤nmax,where qmax, Qmax, and nmax represent upper limits of q, Q and n, respectively, and
(20)q=12ρV2,Q=Cρ0.5V3.05,n=T+D+L+Zmg,where C=1.291×10−4 is the heat flux constant and T+D+L+Z is the vector sum of thrust, drag, lift, and side forces.


## 3.2. Trajectory Constraints in the Midcourse Guidance Stage

For the midcourse guidance stage, the cost function is to minimize the fuel consumption of the third stage interceptor during flight, that is, to maximize the remaining mass of the third stage interceptor at the end of the midcourse guidance, i.e.(21)J=−mf.The state equation of the interceptor in this flight phase is shown in Equation (1). In order to achieve the handover between midcourse and terminal guidance, it is necessary to ensure that the third stage interceptor can detect the target at the end of the midcourse guidance stage, that is, to ensure that the detection range of the third stage interceptor’s seeker can successfully cover the HPR of the target, which is displayed in Figure 3.Figure 3
The handover between midcourse and terminal guidance.As shown in Figure3, it is assumed that the HPR of the target is a spherical region with radius r in space, the coordinate of its spherical center Tp is xp,yp,zp, the position of the interceptor’s centroid M is x,y,z, and the distance between M and Tp is d. In order to ensure that the seeker’s detection field can cover the target’s HPR more effectively, it is necessary to make the axis direction of the seeker’s detection field consistent with the line direction of MTp, that is, to make the pitch angle ϑ and yaw angle ψ of the interceptor meet the following geometric relations.
(22)ϑ=arcsinyp−yd,ψ′=arcsinz−zpd2−yp−y2,ψ=ψ′,xp≥x,sgnz−zp·π−ψ′,xp<x.Furthermore, assuming that the detection field angle of the seeker isϕ and the maximum detection range is Rmax, the following constraints should be met to ensure that the detection field can completely cover the HPR of the target.
(23)d≥rtanϕ,d≤Rmax−r.Equations (22) and (23) are the terminal constraints of interceptor in the midcourse guidance stage. In addition, the third interceptor also needs to meet various process constraints in the midcourse guidance stage, as shown in Equations (19) and (20).


## 4. Trajectory Optimization Based on Improved Adaptive Hp Algorithm

In order to intercept the target successfully, the current interceptor mostly adopts the compound guidance mode. Under the guidance system, the trajectories of programmed guidance stage and midcourse guidance stage are optimized offline in advance combining with the prediction of target position. Since the trajectories of interceptor in these two stages do not need to be acquired online, the trajectories can be solved by offline optimization algorithm, and hp pseudospectral method is commonly used to solve this problem. In order to improve the solving efficiency of trajectory, this paper proposes an improved adaptive hp pseudospectral method based on the traditional hp method, which mainly improves the mesh refinement process in two aspects: (1) By evaluating the change rate of the control, the grid is prerefined near the sudden change points of the control according to certain rules. (2) According to the curvature value of the state curve, the mesh interval is divided nonuniformly to ensure that more mesh points are allocated to the position with large curvature value.


### 4.1. Optimization Problem Transformation Based on LGR Collocation Method

In order to solve the trajectory optimization problem (16), the idea of hp pseudospectral method [20–22] is adopted here. Firstly, K+1 mesh points T0,T1,⋯,TK are selected on the whole time interval 0,tf, and then 0,tf is divided into K mesh intervals Sk=Tk−1,Tk,k=1,2,⋯,K, where T0=0, TK=tf, and ∪k=1KSk=0,tf. At the same time, let xkt and ukt represent the state variable and control variable on the mesh interval Sk, so the original optimization problem (16) can be transformed into the following equivalent problem B.
(24)BminJ=ExKTKs.t.dxktdt=fxkt,uktexKTK=0hxkt,ukt≤0k=1,⋯,K,It is worth noting that since the state variable is continuous in time interval0,tf, each mesh point in Equation (24) needs to meet the following conditions.
(25)xkTk−=xk+1Tk+,k=1,2,⋯,K−1.Then problem B is transformed into NLP by LGR collocation method. First, we need to convert the time variablet in each mesh interval Sk to interval −1,1. Here, we introduce a new time variable τ to perform the following transformation.
(26)τ=2tTk−Tk−1−Tk+Tk−1Tk−Tk−1k=1,2,⋯,K.Lethk represent the length of each mesh interval Sk, i.e., hk=Tk−Tk−1. Then, combining Equations (24)–(26), problem B can be equivalent to the following problem.
(27)BminJ=ExK+1s.t.dxktdt=hk2fxkt,uktexK+1=0hxkt,ukt≤0k=1,⋯,Kxk+1=xk+1−1,k=1,⋯,K−1.Next,Nk LGR collocation points τ1k,⋯,τNkk∈−1,1 are selected on each mesh interval Sk. In order to discretize Equation (27), the selected LGR collocation points τ1k,⋯,τNkk and noncollocation points τNk+1k=+1 are used as sampling points to perform Lagrange interpolation for state variable xkτ on interval Sk. The Lagrange interpolation polynomial based on the Nk+1 sampling points is expressed as follows
(28)ljkτ=∏i=1i≠jNk+1τ−τikτjk−τik,j=1,2,⋯,Nk+1.Then, continuous state variablexkτ on mesh interval Sk can be expressed as
(29)xkτ≈Xkτ=∑j=1Nk+1Xjkljkτ,where Xjk,j=1,2,⋯,Nk+1 is the approximate value of xktjk. Take the derivative of Equation (29) at the collocation points τ1k,⋯,τNkk, then
(30)Ẋktik=∑j=1Nk+1Xjkl̇jktik=∑j=1Nk+1XjkDijk,i=1,2,⋯,Nk.where Dijk=l̇jkτik. Next, let Ujk,j=1,⋯,Nk represent the approximate value of the control variable uktjk. Combining Equations. (27)–(30), problem B can be transformed into the following NLP, denoted as problem BN(31)BNminJ=EXNK+1Ks.t.∑j=1Nk+1XjkDijk=hk2fXik,UikeXNK+1K=0hXik,Uik≤0i=1,⋯,Nkk=1,⋯,KXNk+1k=X1k+1k=1,⋯,K−1.For problemBN, the NLP solver SNOPT [15] can be used to solve it. However, solving it only once may not achieve the solution accuracy required by the original problem B. Therefore, it is necessary to continuously adjust the position of mesh points T0,T1,⋯,TK (h method) or the number of LGR collocation points in each mesh interval (p method) to solve it iteratively until the required solution accuracy is achieved. In view of this, this paper proposes an improved adaptive hp mesh refinement method.


### 4.2. Improved Adaptive Hp Mesh Refinement Method

When solving the trajectory optimization problem in this paper, the obtained control variables may suffer from sudden change. However, the traditional hp method [20–22] will allocate a large number of LGR collocation points near the sudden change points or segment the grid frequently, which will increase the solving dimension of NLP and affect the solving efficiency. In addition, the traditional hp method often performs average mesh segmentation in the mesh refinement process, and the excessive curvature of state curve will further increase the solving burden of hp method. To solve the above problems, an improved adaptive hp mesh refinement method is proposed in this paper, and the detailed steps are as follows.


#### 4.2.1. Evaluation of Solution Error

When solving problemBN with NLP solver, it is necessary to determine whether the obtained solution can make the error between Xkτ and xkτ in Equation (29) within accuracy tolerance range ε. When it exceeds ε, the original mesh needs to be refined and the optimization problem needs to be solved again. Here we adopt the following method to evaluate the solution error in each mesh interval.First, an LGR collocation point is added in each mesh interval. In the mesh intervalSk, the new Γk collocation points are denoted as τ^1k,⋯,τ^Γkk, where Γk=Nk+1 and τ^1k=−1. At the same time, take the noncollocation point τ^Γk+1k=+1.Then, the state variableXkτ^ik at τ^1k,⋯,τ^Γk+1k and the control variable Ukτ^ik at τ^1k,⋯,τ^Γkk are obtained by Lagrange interpolation combined with the obtained numerical solution Xjk,Ujk and the corresponding sampling point τ1k,⋯,τNk+1k. These two can be expressed as
(32)Xkτ^ik=∑j=1Nk+1Xjkljkτ^ik,i=1,⋯,Γk+1,ljkτ=∏i=1i≠jNk+1τ−τikτjk−τik,j=1,2,⋯,Nk+1,(33)Ukτ^ik=∑j=1NkUjkljkτ^ik,i=1,⋯,Γk,ljkτ=∏i=1i≠jNkτ−τikτjk−τik,j=1,2,⋯,Nk.Finally, combiningXkτ^ik and Ukτ^ik obtained above, the state variable X∧kt^ik at τ^1k,⋯,τ^Γk+1k is reintegrated by using the integral form of the discrete state equation.So far, we can useXkτ^ik−X∧kτ^ik to approximate the error between numerical solution Xkτ and analytical solution xkτ of problem B, as shown below
(34)elkτ^ik=Xlkτ^ik−X^lkτ^ik,emaxk=maxl=1,⋯,nxi=1,⋯,Γk+1elkτ^ik.where Xlk and X^lk represent the l th element of vectors Xk and X∧k, respectively.


#### 4.2.2. The Mesh Prerefinement around Control Sudden Change Points

When hp pseudospectral method is used to solve problemBN iteratively, the sudden change of control variables will greatly increase the number of iterations and affect the efficiency of solving. In order to effectively solve the above problems in the improved adaptive hp method in this paper, we prerefine the mesh based on the control sudden change points. The basic idea is to accurately locate the control sudden change points, and then prioritize the mesh refinement near these points to improve the solving efficiency.First, in order to locate the sudden change points, the first derivative of time is taken attii=1,2,⋯,Nk with respect to the control variable ukt.
(35)u̇kti=ukti+1−uktiti+1−ti.where ti=τi+1hk/2+Tk−1. Equation (35) can be further written as
(36)u̇kti=2hk·Ui+1k−Uikτi+1−τi.Then, according to the derivative of the control variable and the set thresholdη, the interval where the control variable changes too fast is selected. Assume that the derivatives of control variables at each collocation point in interval Sk are shown in the figure below.As shown in Figure4, when the derivative of the control variable exceeds η (such as tα and tβ−1,tβ,tβ+1 in the figure), we believe that there may be control sudden change points in this area, so we need to refine the mesh in advance. We will add new mesh points tα−1,tα+1,tβ−2,tβ+2 to further refine region tα−1,tα+1 and tβ−2,tβ+2. In order to describe the subsequent mesh refinement process, the new mesh points combination obtained is denoted as T0′,T1′,⋯,TK′′ and the new mesh interval is denoted as Sk′=Tk−1′,Tk′,k=1,2,⋯,K′.Figure 4
The derivative of control variable.


#### 4.2.3. Mesh Refinement Classification Criterion Based On State Variable Curvature

When hp method is used to solve the optimal control problem, different methods should be used to refine the mesh according to the characteristics of different mesh [31]. In some mesh intervals, the state variables change dramatically, which is suitable for segmentation, that is, h method; For those mesh regions with relatively stable state changes, it is appropriate to refine them by adding collocation points, namely p method. In the improved adaptive hp method in this paper, we take the curvature of state variables as the criterion of classification.The set of collocation points in the new mesh intervalSk′ is denoted as Ωtk, and the state variables at each collocation points are expressed as Ωxtk. Based on Ωtk and Ωxtk, N points with equal time interval in interval Tk−1′,Tk′ are selected for Lagrange interpolation, and the interpolation points obtained are marked as χik,Yik, where χik=Tk−1′+i⋅Tk′−Tk−1′/N−1. The curvature of each point can be calculated by using the coordinates of interpolation points
(37)κik=Y¨ik1+Ẏik23/2i=0,1,⋯,N−1.Combined with the curvature obtained, the classification refinement criterion is set as(38)Crik=maxi=0,1,⋯,N−1κik∑i=0N−1κik/N.If the curvature threshold parameter is set asρ, and if there is Crik>ρ, it can be considered that the state variable changes too fast in the mesh interval and needs to be refined by segmenting. On the contrary, if there is Crik≤ρ, it is necessary to refine the mesh by adding collocation points.


#### 4.2.4. Increasing Collocation Points Based on Maximum Approximation Error

According to the above judgment, if mesh intervalSk′ needs to further increase the number of collocation points, based on the maximum approximation error emaxk in Equation (34), we increase the number of collocation points according to the following equation
(39)Nk′=Nk+log10emaxkε.where Nk is the original number of collocation points in mesh interval Sk′, Nk′ represents the number of collocation points after the increase, and ε is the set accuracy tolerance. · is the least integer function.


#### 4.2.5. Mesh Segmentation Method Based on Cumulative Sum of Curvature Values

In the traditional adaptive hp method, the mesh interval that needs to be segmented is usually evenly divided into several subintervals, which do not maximize the efficiency of segmentation. In order to further improve the efficiency of mesh refinement, we design a segmentation method based on the accumulative sum of curvature value, which can allocate more mesh points at the position with larger curvature of state curve.First, estimate the number of subintervals that need to be added. The number of subintervals segmented within mesh intervalSk′ is determined by the following equation.
(40)Mk=max2,σ·log10emaxkε,where max∘,∗ represents the maximum value of ∘ and ∗, and σ is the set product factor.Then, based on the curvature of state curve in mesh intervalSk′, the position of new mesh points is determined. According to Equation (40), interval Sk′ needs to be divided into Mk subintervals, that is, Mk−1 mesh points should be added between Tk−1′ and Tk′. In order to determine the positions of these mesh points, the curvature value κik at each interpolation point χik,Yik in Equation (37) needs to be taken into account to ensure that more dense mesh points are configured at positions where the state variables change dramatically.According to the above principles, the mesh point configuration rule of the improved adaptive hp method in this paper is as follows:Mk−1 points are selected from N interpolation points χik in interval Sk′, denoted as mik,i=1,⋯,Mk−1, so that the curvature sum of interpolation points in each subinterval mi−1k,mik is approximately equal. According to the rule, the selection algorithm of newly added mesh points is as follows:
(41)mik=χjk∑ι=0jκιk≥i·S¯,∑ι=0j−1κιk<i·S¯,where S¯=∑i=0N−1κik/Mk.


#### 4.2.6. Algorithm Flow

At this point, the key steps of the improved adaptive hp method proposed in this paper have been explained, and the description of the algorithm can be summarized as Algorithm1.Algorithm 1: The description of the proposed algorithm.
The improved adaptive hp pseudospectral methodStep 1: Select initial mesh points T0,T1,⋯,TK and collocation points τ1k,⋯,τNkk.Step 2: Solve problem BN by NLP solver, get the preliminary solution Xjk,Ujk.Step 3: Based on the solution obtained in Step 2, evaluate whether the solution error emaxk is greater than ε .Step 4: If the result of Step 3 is no, end the iteration. Otherwise, jump to Step 5.Step 5: Refined mesh in advance based on change rate of control, get new mesh points T0′,T1′,⋯TK′′.Step 6: Evaluate whether the curvature Crik of state Xik is greater than ρ.Step 7: If the result of Step 6 is yes, divide mesh interval based on the accumulation of the curvature. Otherwise, increase collocation points.Step 8: Based on the new mesh points and collocation points, go back to Step 2 and solve again.


## 4.1. Optimization Problem Transformation Based on LGR Collocation Method

In order to solve the trajectory optimization problem (16), the idea of hp pseudospectral method [20–22] is adopted here. Firstly, K+1 mesh points T0,T1,⋯,TK are selected on the whole time interval 0,tf, and then 0,tf is divided into K mesh intervals Sk=Tk−1,Tk,k=1,2,⋯,K, where T0=0, TK=tf, and ∪k=1KSk=0,tf. At the same time, let xkt and ukt represent the state variable and control variable on the mesh interval Sk, so the original optimization problem (16) can be transformed into the following equivalent problem B.
(24)BminJ=ExKTKs.t.dxktdt=fxkt,uktexKTK=0hxkt,ukt≤0k=1,⋯,K,It is worth noting that since the state variable is continuous in time interval0,tf, each mesh point in Equation (24) needs to meet the following conditions.
(25)xkTk−=xk+1Tk+,k=1,2,⋯,K−1.Then problem B is transformed into NLP by LGR collocation method. First, we need to convert the time variablet in each mesh interval Sk to interval −1,1. Here, we introduce a new time variable τ to perform the following transformation.
(26)τ=2tTk−Tk−1−Tk+Tk−1Tk−Tk−1k=1,2,⋯,K.Lethk represent the length of each mesh interval Sk, i.e., hk=Tk−Tk−1. Then, combining Equations (24)–(26), problem B can be equivalent to the following problem.
(27)BminJ=ExK+1s.t.dxktdt=hk2fxkt,uktexK+1=0hxkt,ukt≤0k=1,⋯,Kxk+1=xk+1−1,k=1,⋯,K−1.Next,Nk LGR collocation points τ1k,⋯,τNkk∈−1,1 are selected on each mesh interval Sk. In order to discretize Equation (27), the selected LGR collocation points τ1k,⋯,τNkk and noncollocation points τNk+1k=+1 are used as sampling points to perform Lagrange interpolation for state variable xkτ on interval Sk. The Lagrange interpolation polynomial based on the Nk+1 sampling points is expressed as follows
(28)ljkτ=∏i=1i≠jNk+1τ−τikτjk−τik,j=1,2,⋯,Nk+1.Then, continuous state variablexkτ on mesh interval Sk can be expressed as
(29)xkτ≈Xkτ=∑j=1Nk+1Xjkljkτ,where Xjk,j=1,2,⋯,Nk+1 is the approximate value of xktjk. Take the derivative of Equation (29) at the collocation points τ1k,⋯,τNkk, then
(30)Ẋktik=∑j=1Nk+1Xjkl̇jktik=∑j=1Nk+1XjkDijk,i=1,2,⋯,Nk.where Dijk=l̇jkτik. Next, let Ujk,j=1,⋯,Nk represent the approximate value of the control variable uktjk. Combining Equations. (27)–(30), problem B can be transformed into the following NLP, denoted as problem BN(31)BNminJ=EXNK+1Ks.t.∑j=1Nk+1XjkDijk=hk2fXik,UikeXNK+1K=0hXik,Uik≤0i=1,⋯,Nkk=1,⋯,KXNk+1k=X1k+1k=1,⋯,K−1.For problemBN, the NLP solver SNOPT [15] can be used to solve it. However, solving it only once may not achieve the solution accuracy required by the original problem B. Therefore, it is necessary to continuously adjust the position of mesh points T0,T1,⋯,TK (h method) or the number of LGR collocation points in each mesh interval (p method) to solve it iteratively until the required solution accuracy is achieved. In view of this, this paper proposes an improved adaptive hp mesh refinement method.


## 4.2. Improved Adaptive Hp Mesh Refinement Method

When solving the trajectory optimization problem in this paper, the obtained control variables may suffer from sudden change. However, the traditional hp method [20–22] will allocate a large number of LGR collocation points near the sudden change points or segment the grid frequently, which will increase the solving dimension of NLP and affect the solving efficiency. In addition, the traditional hp method often performs average mesh segmentation in the mesh refinement process, and the excessive curvature of state curve will further increase the solving burden of hp method. To solve the above problems, an improved adaptive hp mesh refinement method is proposed in this paper, and the detailed steps are as follows.


### 4.2.1. Evaluation of Solution Error

When solving problemBN with NLP solver, it is necessary to determine whether the obtained solution can make the error between Xkτ and xkτ in Equation (29) within accuracy tolerance range ε. When it exceeds ε, the original mesh needs to be refined and the optimization problem needs to be solved again. Here we adopt the following method to evaluate the solution error in each mesh interval.First, an LGR collocation point is added in each mesh interval. In the mesh intervalSk, the new Γk collocation points are denoted as τ^1k,⋯,τ^Γkk, where Γk=Nk+1 and τ^1k=−1. At the same time, take the noncollocation point τ^Γk+1k=+1.Then, the state variableXkτ^ik at τ^1k,⋯,τ^Γk+1k and the control variable Ukτ^ik at τ^1k,⋯,τ^Γkk are obtained by Lagrange interpolation combined with the obtained numerical solution Xjk,Ujk and the corresponding sampling point τ1k,⋯,τNk+1k. These two can be expressed as
(32)Xkτ^ik=∑j=1Nk+1Xjkljkτ^ik,i=1,⋯,Γk+1,ljkτ=∏i=1i≠jNk+1τ−τikτjk−τik,j=1,2,⋯,Nk+1,(33)Ukτ^ik=∑j=1NkUjkljkτ^ik,i=1,⋯,Γk,ljkτ=∏i=1i≠jNkτ−τikτjk−τik,j=1,2,⋯,Nk.Finally, combiningXkτ^ik and Ukτ^ik obtained above, the state variable X∧kt^ik at τ^1k,⋯,τ^Γk+1k is reintegrated by using the integral form of the discrete state equation.So far, we can useXkτ^ik−X∧kτ^ik to approximate the error between numerical solution Xkτ and analytical solution xkτ of problem B, as shown below
(34)elkτ^ik=Xlkτ^ik−X^lkτ^ik,emaxk=maxl=1,⋯,nxi=1,⋯,Γk+1elkτ^ik.where Xlk and X^lk represent the l th element of vectors Xk and X∧k, respectively.


### 4.2.2. The Mesh Prerefinement around Control Sudden Change Points

When hp pseudospectral method is used to solve problemBN iteratively, the sudden change of control variables will greatly increase the number of iterations and affect the efficiency of solving. In order to effectively solve the above problems in the improved adaptive hp method in this paper, we prerefine the mesh based on the control sudden change points. The basic idea is to accurately locate the control sudden change points, and then prioritize the mesh refinement near these points to improve the solving efficiency.First, in order to locate the sudden change points, the first derivative of time is taken attii=1,2,⋯,Nk with respect to the control variable ukt.
(35)u̇kti=ukti+1−uktiti+1−ti.where ti=τi+1hk/2+Tk−1. Equation (35) can be further written as
(36)u̇kti=2hk·Ui+1k−Uikτi+1−τi.Then, according to the derivative of the control variable and the set thresholdη, the interval where the control variable changes too fast is selected. Assume that the derivatives of control variables at each collocation point in interval Sk are shown in the figure below.As shown in Figure4, when the derivative of the control variable exceeds η (such as tα and tβ−1,tβ,tβ+1 in the figure), we believe that there may be control sudden change points in this area, so we need to refine the mesh in advance. We will add new mesh points tα−1,tα+1,tβ−2,tβ+2 to further refine region tα−1,tα+1 and tβ−2,tβ+2. In order to describe the subsequent mesh refinement process, the new mesh points combination obtained is denoted as T0′,T1′,⋯,TK′′ and the new mesh interval is denoted as Sk′=Tk−1′,Tk′,k=1,2,⋯,K′.Figure 4
The derivative of control variable.


### 4.2.3. Mesh Refinement Classification Criterion Based On State Variable Curvature

When hp method is used to solve the optimal control problem, different methods should be used to refine the mesh according to the characteristics of different mesh [31]. In some mesh intervals, the state variables change dramatically, which is suitable for segmentation, that is, h method; For those mesh regions with relatively stable state changes, it is appropriate to refine them by adding collocation points, namely p method. In the improved adaptive hp method in this paper, we take the curvature of state variables as the criterion of classification.The set of collocation points in the new mesh intervalSk′ is denoted as Ωtk, and the state variables at each collocation points are expressed as Ωxtk. Based on Ωtk and Ωxtk, N points with equal time interval in interval Tk−1′,Tk′ are selected for Lagrange interpolation, and the interpolation points obtained are marked as χik,Yik, where χik=Tk−1′+i⋅Tk′−Tk−1′/N−1. The curvature of each point can be calculated by using the coordinates of interpolation points
(37)κik=Y¨ik1+Ẏik23/2i=0,1,⋯,N−1.Combined with the curvature obtained, the classification refinement criterion is set as(38)Crik=maxi=0,1,⋯,N−1κik∑i=0N−1κik/N.If the curvature threshold parameter is set asρ, and if there is Crik>ρ, it can be considered that the state variable changes too fast in the mesh interval and needs to be refined by segmenting. On the contrary, if there is Crik≤ρ, it is necessary to refine the mesh by adding collocation points.


### 4.2.4. Increasing Collocation Points Based on Maximum Approximation Error

According to the above judgment, if mesh intervalSk′ needs to further increase the number of collocation points, based on the maximum approximation error emaxk in Equation (34), we increase the number of collocation points according to the following equation
(39)Nk′=Nk+log10emaxkε.where Nk is the original number of collocation points in mesh interval Sk′, Nk′ represents the number of collocation points after the increase, and ε is the set accuracy tolerance. · is the least integer function.


### 4.2.5. Mesh Segmentation Method Based on Cumulative Sum of Curvature Values

In the traditional adaptive hp method, the mesh interval that needs to be segmented is usually evenly divided into several subintervals, which do not maximize the efficiency of segmentation. In order to further improve the efficiency of mesh refinement, we design a segmentation method based on the accumulative sum of curvature value, which can allocate more mesh points at the position with larger curvature of state curve.First, estimate the number of subintervals that need to be added. The number of subintervals segmented within mesh intervalSk′ is determined by the following equation.
(40)Mk=max2,σ·log10emaxkε,where max∘,∗ represents the maximum value of ∘ and ∗, and σ is the set product factor.Then, based on the curvature of state curve in mesh intervalSk′, the position of new mesh points is determined. According to Equation (40), interval Sk′ needs to be divided into Mk subintervals, that is, Mk−1 mesh points should be added between Tk−1′ and Tk′. In order to determine the positions of these mesh points, the curvature value κik at each interpolation point χik,Yik in Equation (37) needs to be taken into account to ensure that more dense mesh points are configured at positions where the state variables change dramatically.According to the above principles, the mesh point configuration rule of the improved adaptive hp method in this paper is as follows:Mk−1 points are selected from N interpolation points χik in interval Sk′, denoted as mik,i=1,⋯,Mk−1, so that the curvature sum of interpolation points in each subinterval mi−1k,mik is approximately equal. According to the rule, the selection algorithm of newly added mesh points is as follows:
(41)mik=χjk∑ι=0jκιk≥i·S¯,∑ι=0j−1κιk<i·S¯,where S¯=∑i=0N−1κik/Mk.


### 4.2.6. Algorithm Flow

At this point, the key steps of the improved adaptive hp method proposed in this paper have been explained, and the description of the algorithm can be summarized as Algorithm1.Algorithm 1: The description of the proposed algorithm.
The improved adaptive hp pseudospectral methodStep 1: Select initial mesh points T0,T1,⋯,TK and collocation points τ1k,⋯,τNkk.Step 2: Solve problem BN by NLP solver, get the preliminary solution Xjk,Ujk.Step 3: Based on the solution obtained in Step 2, evaluate whether the solution error emaxk is greater than ε .Step 4: If the result of Step 3 is no, end the iteration. Otherwise, jump to Step 5.Step 5: Refined mesh in advance based on change rate of control, get new mesh points T0′,T1′,⋯TK′′.Step 6: Evaluate whether the curvature Crik of state Xik is greater than ρ.Step 7: If the result of Step 6 is yes, divide mesh interval based on the accumulation of the curvature. Otherwise, increase collocation points.Step 8: Based on the new mesh points and collocation points, go back to Step 2 and solve again.


## 4.2.1. Evaluation of Solution Error

When solving problemBN with NLP solver, it is necessary to determine whether the obtained solution can make the error between Xkτ and xkτ in Equation (29) within accuracy tolerance range ε. When it exceeds ε, the original mesh needs to be refined and the optimization problem needs to be solved again. Here we adopt the following method to evaluate the solution error in each mesh interval.First, an LGR collocation point is added in each mesh interval. In the mesh intervalSk, the new Γk collocation points are denoted as τ^1k,⋯,τ^Γkk, where Γk=Nk+1 and τ^1k=−1. At the same time, take the noncollocation point τ^Γk+1k=+1.Then, the state variableXkτ^ik at τ^1k,⋯,τ^Γk+1k and the control variable Ukτ^ik at τ^1k,⋯,τ^Γkk are obtained by Lagrange interpolation combined with the obtained numerical solution Xjk,Ujk and the corresponding sampling point τ1k,⋯,τNk+1k. These two can be expressed as
(32)Xkτ^ik=∑j=1Nk+1Xjkljkτ^ik,i=1,⋯,Γk+1,ljkτ=∏i=1i≠jNk+1τ−τikτjk−τik,j=1,2,⋯,Nk+1,(33)Ukτ^ik=∑j=1NkUjkljkτ^ik,i=1,⋯,Γk,ljkτ=∏i=1i≠jNkτ−τikτjk−τik,j=1,2,⋯,Nk.Finally, combiningXkτ^ik and Ukτ^ik obtained above, the state variable X∧kt^ik at τ^1k,⋯,τ^Γk+1k is reintegrated by using the integral form of the discrete state equation.So far, we can useXkτ^ik−X∧kτ^ik to approximate the error between numerical solution Xkτ and analytical solution xkτ of problem B, as shown below
(34)elkτ^ik=Xlkτ^ik−X^lkτ^ik,emaxk=maxl=1,⋯,nxi=1,⋯,Γk+1elkτ^ik.where Xlk and X^lk represent the l th element of vectors Xk and X∧k, respectively.


## 4.2.2. The Mesh Prerefinement around Control Sudden Change Points

When hp pseudospectral method is used to solve problemBN iteratively, the sudden change of control variables will greatly increase the number of iterations and affect the efficiency of solving. In order to effectively solve the above problems in the improved adaptive hp method in this paper, we prerefine the mesh based on the control sudden change points. The basic idea is to accurately locate the control sudden change points, and then prioritize the mesh refinement near these points to improve the solving efficiency.First, in order to locate the sudden change points, the first derivative of time is taken attii=1,2,⋯,Nk with respect to the control variable ukt.
(35)u̇kti=ukti+1−uktiti+1−ti.where ti=τi+1hk/2+Tk−1. Equation (35) can be further written as
(36)u̇kti=2hk·Ui+1k−Uikτi+1−τi.Then, according to the derivative of the control variable and the set thresholdη, the interval where the control variable changes too fast is selected. Assume that the derivatives of control variables at each collocation point in interval Sk are shown in the figure below.As shown in Figure4, when the derivative of the control variable exceeds η (such as tα and tβ−1,tβ,tβ+1 in the figure), we believe that there may be control sudden change points in this area, so we need to refine the mesh in advance. We will add new mesh points tα−1,tα+1,tβ−2,tβ+2 to further refine region tα−1,tα+1 and tβ−2,tβ+2. In order to describe the subsequent mesh refinement process, the new mesh points combination obtained is denoted as T0′,T1′,⋯,TK′′ and the new mesh interval is denoted as Sk′=Tk−1′,Tk′,k=1,2,⋯,K′.Figure 4
The derivative of control variable.


## 4.2.3. Mesh Refinement Classification Criterion Based On State Variable Curvature

When hp method is used to solve the optimal control problem, different methods should be used to refine the mesh according to the characteristics of different mesh [31]. In some mesh intervals, the state variables change dramatically, which is suitable for segmentation, that is, h method; For those mesh regions with relatively stable state changes, it is appropriate to refine them by adding collocation points, namely p method. In the improved adaptive hp method in this paper, we take the curvature of state variables as the criterion of classification.The set of collocation points in the new mesh intervalSk′ is denoted as Ωtk, and the state variables at each collocation points are expressed as Ωxtk. Based on Ωtk and Ωxtk, N points with equal time interval in interval Tk−1′,Tk′ are selected for Lagrange interpolation, and the interpolation points obtained are marked as χik,Yik, where χik=Tk−1′+i⋅Tk′−Tk−1′/N−1. The curvature of each point can be calculated by using the coordinates of interpolation points
(37)κik=Y¨ik1+Ẏik23/2i=0,1,⋯,N−1.Combined with the curvature obtained, the classification refinement criterion is set as(38)Crik=maxi=0,1,⋯,N−1κik∑i=0N−1κik/N.If the curvature threshold parameter is set asρ, and if there is Crik>ρ, it can be considered that the state variable changes too fast in the mesh interval and needs to be refined by segmenting. On the contrary, if there is Crik≤ρ, it is necessary to refine the mesh by adding collocation points.


## 4.2.4. Increasing Collocation Points Based on Maximum Approximation Error

According to the above judgment, if mesh intervalSk′ needs to further increase the number of collocation points, based on the maximum approximation error emaxk in Equation (34), we increase the number of collocation points according to the following equation
(39)Nk′=Nk+log10emaxkε.where Nk is the original number of collocation points in mesh interval Sk′, Nk′ represents the number of collocation points after the increase, and ε is the set accuracy tolerance. · is the least integer function.


## 4.2.5. Mesh Segmentation Method Based on Cumulative Sum of Curvature Values

In the traditional adaptive hp method, the mesh interval that needs to be segmented is usually evenly divided into several subintervals, which do not maximize the efficiency of segmentation. In order to further improve the efficiency of mesh refinement, we design a segmentation method based on the accumulative sum of curvature value, which can allocate more mesh points at the position with larger curvature of state curve.First, estimate the number of subintervals that need to be added. The number of subintervals segmented within mesh intervalSk′ is determined by the following equation.
(40)Mk=max2,σ·log10emaxkε,where max∘,∗ represents the maximum value of ∘ and ∗, and σ is the set product factor.Then, based on the curvature of state curve in mesh intervalSk′, the position of new mesh points is determined. According to Equation (40), interval Sk′ needs to be divided into Mk subintervals, that is, Mk−1 mesh points should be added between Tk−1′ and Tk′. In order to determine the positions of these mesh points, the curvature value κik at each interpolation point χik,Yik in Equation (37) needs to be taken into account to ensure that more dense mesh points are configured at positions where the state variables change dramatically.According to the above principles, the mesh point configuration rule of the improved adaptive hp method in this paper is as follows:Mk−1 points are selected from N interpolation points χik in interval Sk′, denoted as mik,i=1,⋯,Mk−1, so that the curvature sum of interpolation points in each subinterval mi−1k,mik is approximately equal. According to the rule, the selection algorithm of newly added mesh points is as follows:
(41)mik=χjk∑ι=0jκιk≥i·S¯,∑ι=0j−1κιk<i·S¯,where S¯=∑i=0N−1κik/Mk.


## 4.2.6. Algorithm Flow

At this point, the key steps of the improved adaptive hp method proposed in this paper have been explained, and the description of the algorithm can be summarized as Algorithm1.Algorithm 1: The description of the proposed algorithm.
The improved adaptive hp pseudospectral methodStep 1: Select initial mesh points T0,T1,⋯,TK and collocation points τ1k,⋯,τNkk.Step 2: Solve problem BN by NLP solver, get the preliminary solution Xjk,Ujk.Step 3: Based on the solution obtained in Step 2, evaluate whether the solution error emaxk is greater than ε .Step 4: If the result of Step 3 is no, end the iteration. Otherwise, jump to Step 5.Step 5: Refined mesh in advance based on change rate of control, get new mesh points T0′,T1′,⋯TK′′.Step 6: Evaluate whether the curvature Crik of state Xik is greater than ρ.Step 7: If the result of Step 6 is yes, divide mesh interval based on the accumulation of the curvature. Otherwise, increase collocation points.Step 8: Based on the new mesh points and collocation points, go back to Step 2 and solve again.


## 5. Simulation Results

Next, based on the interceptor model in Section2 and the trajectory optimization problem in Section 3, the improved adaptive hp method proposed in this paper is used to conduct trajectory optimization simulation. In this part, all the simulation relies on the computer hardware of Lenovo desktop computer, whose processor is Intel(R) Core(TM) i5-4570, main frequency is 3.20GHz, memory is 4.00G, and the software platform of simulation is Matlab2014a. In addition, the algorithm proposed in this paper is developed based on the open-source hp pseudospectral optimization software GPOPS [43], and the NLP solver used is SNOPT [15].


### 5.1. Trajectory Optimization Simulation of Programmed Guidance Stage

Firstly, we use the improved adaptive hp algorithm proposed in this paper to optimize the trajectory of the programmed guidance stage (first and second boost stages). The initial state of interceptor in programmed guidance stage isx0,y0,z0=0,0,0, Vx0,Vy0,Vz0=0,50m/s,0, and m0=17193kg. Its initial control variable is ϑ,ψ,uon−off=π/2,0,1. It is worth noting that since it is stipulated that the engine of the interceptor is always on in the programmed guidance stage, there is uon‐off≡1. Meanwhile, the upper limit of process constraint of interceptor in programmed guidance stage is set as qmax=100kPa, Qmax=5×105J/m2⋅s, and nmax=6.In the improved adaptive hp algorithm proposed in this paper, there are four parameters to be designed:ε, σ, η, and ρ. The smaller the accuracy tolerance range ε is, the higher the accuracy is required. Generally, satisfactory accuracy is achieved when ε is less than 10-3. The larger the product factor σ is, the more mesh points will be added during grid iteration, which will improve the accuracy of numerical solution to a certain extent. However, too large σ will lead to a large number of new mesh points, thus greatly increasing the dimension of NLP problem and affecting the solving speed. Generally, the value of σ can be selected between 1 and 2. The smaller threshold η is, the more sensitive the proposed algorithm is to the change of control, that is, the smaller the change of control will make the grid to be prerefined, and the value between 0.5 and 2 is appropriate. The smaller the curvature threshold ρ is, the easier it is to add mesh points instead of collocation points during grid iteration. Generally, it can be set at 1 ~ 3.According to the above parameter selection principle, relevant parameters in the improved adaptive hp algorithm are selected as follows:ε=10−5, σ=1.2, η=1,and ρ=2.Usually, when setting the initial grid, the grid interval is not segmented, that is, usually only two initial mesh points are selected, one at the initial moment and one at the terminal moment. In addition, although the approximation accuracy of numerical solution can be improved if there are too many initial collocation points in the entire initial grid interval, the solving speed of the algorithm will be seriously affected. Generally, it is appropriate to choose the number of LGR collocation points between 8 and 15. According to this principle, the initial mesh point is set as−1,1, and the initial LGR collocation points are set as 12.In addition, under the same mesh accuracy tolerance rangeε, hp method in [21] is used to solve the problem, and the simulation is compared with the method in this paper. In the simulation results, P1 represents the first boost stage, and P2 represents the second boost stage. As described in Section 3.1, the interceptor also needs to meet the constraints of flight path angle and heading angle at the end of the programmed guidance stage. Here, the expected value θref,ψref of the two is set as 15°,0°. Simulation results are shown in Figure 5 and Table 2.Figure 5
Programmed guidance stage trajectory optimization results.
(a)
Trajectory of the missile(b)
Velocity vs. time(c)
Mass of missile(d)
Flight path angle and heading angle(e)
Control vs. time(f)
Process constraint(g)
Various grids of improved hp(h)
Various grids of hpTable 2
Comparison of mesh refinement results.
MethodCPU time(s)Number of grid iterationsTotal number of collocation pointsTerminal velocity (m/s)Improved hp15.663872968.1hp25.4841112967.2As can be seen from Figures5(a)–5(f), when the improved adaptive hp method proposed in this paper and the hp method in [21] are used to optimize the trajectory of the programmed guidance stage, respectively, the optimization results obtained by the two methods are basically the same. As shown in Figure 5(d), at the end of the programmed guidance stage, the interceptor’s flight path angle converges to 15 degrees and the heading angle is stable at 0 degree, that is, the optimized trajectory obtained by the two methods can meet the preset terminal constraints. Meanwhile, Figure 5(f) shows that both methods can make the interceptor meet the preset process constraints in the flight process of the programmed guidance stage.However, it is worth noting that in this simulation, in order to achieve the same optimization solution accuracy, the hp method needs 4 mesh iterations, while the improved adaptive hp method only needs 3 iterations (see Figures5(g) and 5(h)). By observing Figures 5(g) and 5(h), it can be found that the original grid is segmented uniformly by the traditional hp pseudospectral method during grid iteration. However, the improved hp pseudospectral method proposed in this paper is different from this method. According to the curvature of the state curve, this method divides the original grid into nonuniform segments, so as to ensure that the position with large curvature can be allocated to more LGR points more quickly, and then ensure that the obtained solution can meet the set accuracy more quickly. This is more concretely shown in Table 2. As shown in Table 2, the terminal velocity obtained by trajectory optimization using these methods is basically the same (obtained by the improved adaptive hp method slightly higher). Furthermore, in order to achieve the optimization result, the hp method needs to carry out 4 mesh iterations and configure 111 LGR points, which takes 25.48 s. However, the improved adaptive hp method only needs to carry out 3 iterations and configure 87 LGR collocation points, which takes 15.66 s. The results show that the optimization results obtained by the two methods are basically the same, while the improved hp adaptive method takes less time and has higher efficiency.In addition, the following explanations are needed. The defensive operation background considered in this paper is ultra-long-range interception. When low-orbit satellites or high-power radars find the target, it is still far away from the defensive position and the remaining attack time is long. In order to intercept such targets effectively, the trajectory of the target should be predicted first, and then the possible position of the target in the future can be obtained. Before launching interceptor, it is necessary to offline optimize the trajectories of programmed guidance stage and midcourse guidance stage based on the predicted target position, and then get the nominal trajectory of interceptor, and the algorithm proposed in this paper is mainly for this purpose. Therefore, in this combat scenario, the time for trajectory optimization before launching of interceptor is relatively long, and the CPU time in simulation is within the acceptable range.


### 5.2. Trajectory Optimization Simulation of Midcourse Guidance Stage

The terminal state of the interceptor in Section5.1 is taken as the initial state of the interceptor in the midcourse guidance stage, and the improved adaptive hp method proposed in this paper is also used for the midcourse guidance trajectory optimization simulation. Set the center position of the HPR of target as (800 km, 30 km, and 35 km) and radius as 6 km. The detection field angle of the seeker is set as 4∘, and its maximum detection range is set as 100 km. In the simulation, the relevant parameters of the improved adaptive hp algorithm are set as: ε=10−7, σ=1.2, η=1, and ρ=2. The initial mesh point is −1,1, and the number of initial LGR collocation points is 10. In addition, the process constraint upper limit of interceptor is set in accordance with Section 5.1. Similarly, the hp method in [21] is used for trajectory optimization, and the simulation is compared with the method in this paper.Meanwhile, the position and attitude data of the interceptor at the terminal of midcourse guidance are summarized as follows:According to the terminal state of interceptor in Table3, combined with the HPR of target set in the simulation and related parameters of seeker, the simulation is carried out as follows to verify whether the optimized trajectory obtained by the two methods can ensure that the seeker of interceptor successfully detect the target.Table 3
Terminal state of missile in the midcourse flight.
MethodTerminal position (km)Pitch angle (deg)Yaw angle (deg)Improved hp(708.90, 27.47, 11.98)1.59-14.18hp(708.85, 27.53, 12.18)1.55-14.05As can be seen from the above simulation results, when the two methods are used for trajectory optimization simulation of midcourse guidance stage, the obtained optimization results are basically the same (see Figures6(a), 6(b), and 6(d)). As shown in Figure 6(f), when the interceptor flies along the optimized trajectory obtained by the two methods, the dynamic pressure, heat flux density, and overload can all meet the process constraints. In addition, at the end of midcourse guidance stage, the state of the interceptor under the two methods can ensure that the detectable region of the seeker effectively covers the HPR of the target, thus successfully completing the handover between midcourse and terminal guidance (see Figures 7 and 8). Furthermore, the control variables of the two methods are slightly different, especially the switch engine (see Figure 6(e)), which leads to a difference in the fuel consumption of the third stage interceptor in the midcourse guidance stage. The fuel consumption of the improved hp method is less than that of the hp method (see Figure 6(c)).Figure 6
Midcourse guidance stage trajectory optimization results.
(a)
Trajectory of the missile(b)
Velocity vs. time(c)
Mass of missile(d)
Flight path angle and heading angle(e)
Control vs. time(f)
Process constraint(g)
Various grids of improved hp(h)
Various grids of hpFigure 7
Handover of improved hp.Figure 8
Handover of hp.Although the optimization effects of the two methods are almost the same, it is worth noting that the efficiency of the two methods is different when solving trajectory optimization. As shown in Figures6(g) and 6(h), in order to achieve the same optimization accuracy, the hp method requires 3 mesh iterations, while the improved adaptive hp method requires only 2 iterations and requires significantly fewer LGR collocation points. This is more concretely reflected in the comparison in Table 4. As shown in Table 4, 131 LGR collocation points need to be configured for the optimization by hp method, which takes 27.38 s, while the improved adaptive hp method only needs 61 collocation points, which takes 6.97 s. In general, in the trajectory optimization simulation of midcourse guidance stage, the two methods can get the optimized trajectory satisfying the constraints. Meanwhile, the optimization effect of the improved adaptive hp method is slightly better than that of the hp method, and its solving efficiency is significantly higher than that of the hp method.Table 4
Comparison of mesh refinement results.
MethodCPU time(s)Number of grid iterationsTotal number of collocation pointsTerminal mass of missile (kg)Improved hp6.97261875hp27.383131870


## 5.1. Trajectory Optimization Simulation of Programmed Guidance Stage

Firstly, we use the improved adaptive hp algorithm proposed in this paper to optimize the trajectory of the programmed guidance stage (first and second boost stages). The initial state of interceptor in programmed guidance stage isx0,y0,z0=0,0,0, Vx0,Vy0,Vz0=0,50m/s,0, and m0=17193kg. Its initial control variable is ϑ,ψ,uon−off=π/2,0,1. It is worth noting that since it is stipulated that the engine of the interceptor is always on in the programmed guidance stage, there is uon‐off≡1. Meanwhile, the upper limit of process constraint of interceptor in programmed guidance stage is set as qmax=100kPa, Qmax=5×105J/m2⋅s, and nmax=6.In the improved adaptive hp algorithm proposed in this paper, there are four parameters to be designed:ε, σ, η, and ρ. The smaller the accuracy tolerance range ε is, the higher the accuracy is required. Generally, satisfactory accuracy is achieved when ε is less than 10-3. The larger the product factor σ is, the more mesh points will be added during grid iteration, which will improve the accuracy of numerical solution to a certain extent. However, too large σ will lead to a large number of new mesh points, thus greatly increasing the dimension of NLP problem and affecting the solving speed. Generally, the value of σ can be selected between 1 and 2. The smaller threshold η is, the more sensitive the proposed algorithm is to the change of control, that is, the smaller the change of control will make the grid to be prerefined, and the value between 0.5 and 2 is appropriate. The smaller the curvature threshold ρ is, the easier it is to add mesh points instead of collocation points during grid iteration. Generally, it can be set at 1 ~ 3.According to the above parameter selection principle, relevant parameters in the improved adaptive hp algorithm are selected as follows:ε=10−5, σ=1.2, η=1,and ρ=2.Usually, when setting the initial grid, the grid interval is not segmented, that is, usually only two initial mesh points are selected, one at the initial moment and one at the terminal moment. In addition, although the approximation accuracy of numerical solution can be improved if there are too many initial collocation points in the entire initial grid interval, the solving speed of the algorithm will be seriously affected. Generally, it is appropriate to choose the number of LGR collocation points between 8 and 15. According to this principle, the initial mesh point is set as−1,1, and the initial LGR collocation points are set as 12.In addition, under the same mesh accuracy tolerance rangeε, hp method in [21] is used to solve the problem, and the simulation is compared with the method in this paper. In the simulation results, P1 represents the first boost stage, and P2 represents the second boost stage. As described in Section 3.1, the interceptor also needs to meet the constraints of flight path angle and heading angle at the end of the programmed guidance stage. Here, the expected value θref,ψref of the two is set as 15°,0°. Simulation results are shown in Figure 5 and Table 2.Figure 5
Programmed guidance stage trajectory optimization results.
(a)
Trajectory of the missile(b)
Velocity vs. time(c)
Mass of missile(d)
Flight path angle and heading angle(e)
Control vs. time(f)
Process constraint(g)
Various grids of improved hp(h)
Various grids of hpTable 2
Comparison of mesh refinement results.
MethodCPU time(s)Number of grid iterationsTotal number of collocation pointsTerminal velocity (m/s)Improved hp15.663872968.1hp25.4841112967.2As can be seen from Figures5(a)–5(f), when the improved adaptive hp method proposed in this paper and the hp method in [21] are used to optimize the trajectory of the programmed guidance stage, respectively, the optimization results obtained by the two methods are basically the same. As shown in Figure 5(d), at the end of the programmed guidance stage, the interceptor’s flight path angle converges to 15 degrees and the heading angle is stable at 0 degree, that is, the optimized trajectory obtained by the two methods can meet the preset terminal constraints. Meanwhile, Figure 5(f) shows that both methods can make the interceptor meet the preset process constraints in the flight process of the programmed guidance stage.However, it is worth noting that in this simulation, in order to achieve the same optimization solution accuracy, the hp method needs 4 mesh iterations, while the improved adaptive hp method only needs 3 iterations (see Figures5(g) and 5(h)). By observing Figures 5(g) and 5(h), it can be found that the original grid is segmented uniformly by the traditional hp pseudospectral method during grid iteration. However, the improved hp pseudospectral method proposed in this paper is different from this method. According to the curvature of the state curve, this method divides the original grid into nonuniform segments, so as to ensure that the position with large curvature can be allocated to more LGR points more quickly, and then ensure that the obtained solution can meet the set accuracy more quickly. This is more concretely shown in Table 2. As shown in Table 2, the terminal velocity obtained by trajectory optimization using these methods is basically the same (obtained by the improved adaptive hp method slightly higher). Furthermore, in order to achieve the optimization result, the hp method needs to carry out 4 mesh iterations and configure 111 LGR points, which takes 25.48 s. However, the improved adaptive hp method only needs to carry out 3 iterations and configure 87 LGR collocation points, which takes 15.66 s. The results show that the optimization results obtained by the two methods are basically the same, while the improved hp adaptive method takes less time and has higher efficiency.In addition, the following explanations are needed. The defensive operation background considered in this paper is ultra-long-range interception. When low-orbit satellites or high-power radars find the target, it is still far away from the defensive position and the remaining attack time is long. In order to intercept such targets effectively, the trajectory of the target should be predicted first, and then the possible position of the target in the future can be obtained. Before launching interceptor, it is necessary to offline optimize the trajectories of programmed guidance stage and midcourse guidance stage based on the predicted target position, and then get the nominal trajectory of interceptor, and the algorithm proposed in this paper is mainly for this purpose. Therefore, in this combat scenario, the time for trajectory optimization before launching of interceptor is relatively long, and the CPU time in simulation is within the acceptable range.


## 5.2. Trajectory Optimization Simulation of Midcourse Guidance Stage

The terminal state of the interceptor in Section5.1 is taken as the initial state of the interceptor in the midcourse guidance stage, and the improved adaptive hp method proposed in this paper is also used for the midcourse guidance trajectory optimization simulation. Set the center position of the HPR of target as (800 km, 30 km, and 35 km) and radius as 6 km. The detection field angle of the seeker is set as 4∘, and its maximum detection range is set as 100 km. In the simulation, the relevant parameters of the improved adaptive hp algorithm are set as: ε=10−7, σ=1.2, η=1, and ρ=2. The initial mesh point is −1,1, and the number of initial LGR collocation points is 10. In addition, the process constraint upper limit of interceptor is set in accordance with Section 5.1. Similarly, the hp method in [21] is used for trajectory optimization, and the simulation is compared with the method in this paper.Meanwhile, the position and attitude data of the interceptor at the terminal of midcourse guidance are summarized as follows:According to the terminal state of interceptor in Table3, combined with the HPR of target set in the simulation and related parameters of seeker, the simulation is carried out as follows to verify whether the optimized trajectory obtained by the two methods can ensure that the seeker of interceptor successfully detect the target.Table 3
Terminal state of missile in the midcourse flight.
MethodTerminal position (km)Pitch angle (deg)Yaw angle (deg)Improved hp(708.90, 27.47, 11.98)1.59-14.18hp(708.85, 27.53, 12.18)1.55-14.05As can be seen from the above simulation results, when the two methods are used for trajectory optimization simulation of midcourse guidance stage, the obtained optimization results are basically the same (see Figures6(a), 6(b), and 6(d)). As shown in Figure 6(f), when the interceptor flies along the optimized trajectory obtained by the two methods, the dynamic pressure, heat flux density, and overload can all meet the process constraints. In addition, at the end of midcourse guidance stage, the state of the interceptor under the two methods can ensure that the detectable region of the seeker effectively covers the HPR of the target, thus successfully completing the handover between midcourse and terminal guidance (see Figures 7 and 8). Furthermore, the control variables of the two methods are slightly different, especially the switch engine (see Figure 6(e)), which leads to a difference in the fuel consumption of the third stage interceptor in the midcourse guidance stage. The fuel consumption of the improved hp method is less than that of the hp method (see Figure 6(c)).Figure 6
Midcourse guidance stage trajectory optimization results.
(a)
Trajectory of the missile(b)
Velocity vs. time(c)
Mass of missile(d)
Flight path angle and heading angle(e)
Control vs. time(f)
Process constraint(g)
Various grids of improved hp(h)
Various grids of hpFigure 7
Handover of improved hp.Figure 8
Handover of hp.Although the optimization effects of the two methods are almost the same, it is worth noting that the efficiency of the two methods is different when solving trajectory optimization. As shown in Figures6(g) and 6(h), in order to achieve the same optimization accuracy, the hp method requires 3 mesh iterations, while the improved adaptive hp method requires only 2 iterations and requires significantly fewer LGR collocation points. This is more concretely reflected in the comparison in Table 4. As shown in Table 4, 131 LGR collocation points need to be configured for the optimization by hp method, which takes 27.38 s, while the improved adaptive hp method only needs 61 collocation points, which takes 6.97 s. In general, in the trajectory optimization simulation of midcourse guidance stage, the two methods can get the optimized trajectory satisfying the constraints. Meanwhile, the optimization effect of the improved adaptive hp method is slightly better than that of the hp method, and its solving efficiency is significantly higher than that of the hp method.Table 4
Comparison of mesh refinement results.
MethodCPU time(s)Number of grid iterationsTotal number of collocation pointsTerminal mass of missile (kg)Improved hp6.97261875hp27.383131870


## 6. Conclusion

In this paper, an improved adaptive hp pseudospectral method is proposed to solve the trajectory optimization problem of high-speed and long-range interceptor. Different from the mesh refinement process of the traditional adaptive hp method, the proposed method refines the mesh in advance around the control sudden change points. In addition, the method takes the curvature of the state variable curve as the criterion for the nonuniform mesh segmentation iteration; so that more mesh points can be allocated where the curvature is large. These two methods make the allocation of collocation points resources more reasonable and improve the optimization efficiency. The simulation results show that the optimized trajectory of the interceptor can be effectively solved by using the algorithm proposed in this paper, so that the interceptor can successfully realize the handover between midcourse and terminal guidance. And the efficiency of this method in trajectory optimization is higher than that of the traditional adaptive hp method.



---

*Source: 7958272-2022-12-22.xml*