# Optimal Design of Cooperative Penetration Trajectories for Multiaircraft

Authors: Qiangqiang Xu; Jianquan Ge; Tao Yang
Journal: International Journal of Aerospace Engineering
 (2020)

---

## Abstract

At present, two kinds of shortages exist in the research on cooperative combat. One is that radar detection threat (which cannot be ignored) is rarely considered. The other is that limited efforts have been made on the cooperative penetration trajectories under the conditions of long distance, vast airspace, and wide speed range. In order to offset the shortages of the research on cooperative combat, the penetration trajectory optimization method considering the influence of aircraft radar cross-section (RCS) and the cooperative penetration strategy is proposed in this study. Firstly, the RCS data are calculated by the physical optics (PO) method. The radar detection threat model is established considering the influence of the aircraft RCS. Then, a trajectory optimization framework with the dynamic model, constraint conditions, and optimal objectives is formed. Using the hp-adaptive Radau pseudospectral method, the optimal control problem for a single aircraft flight is solved. Finally, a cooperative penetration strategy is proposed to solve the cooperative penetration problem of multiaircraft. The impact time and angle constraints are given, and the virtual target point is introduced for terminal guidance. Two cases are simulated and verified. Simulation results demonstrate that the proposed method is effective. The single aircraft can effectively penetrate, and the multiaircraft can fulfill the requirement of cooperative impact time and angle under the condition of meeting the minimum threat of radar detection.


---

## Body

## 1. Introduction

With the development of military technology, severe challenges on the attack and defense confrontation between the aircraft and the air defense system have become pronounced. The continuous strengthening of the modern defense system has made the traditional single aircraft more difficult to penetrate. Radar is a major threat in the process of aircraft penetration. Therefore, stealth performance plays an important role. In order to improve the stealth performance, on the one hand, some technical means can be adopted to reduce their own characteristic signals. On the other hand, if the target characteristic signals cannot be directly changed, the probability of radar detection can be reduced through trajectory planning. As multiaircraft cooperative combat can effectively improve the intelligent level of the aircraft and the overall operational effectiveness of the formation, it has gradually become a research hotspot in the military field.In the research of aircraft stealth trajectory planning, Xu et al. and Lin et al. [1, 2] only consider the influence of the distance factor between the aircraft and the radar, and the radar cross-section (RCS) of the aircraft is set as an independent parameter of the attitude. However, in the actual battlefield environment, the angle of the aircraft relative to the incident wave of the radar is constantly changing, which leads to the RCS value changing in accordance with the attitude. In response to this, Chen et al. [3] proposed a 3D low-observable trajectory optimization method based on RCS, where the RCS data model is constructed by B-spline fitting. Liu et al. [4, 5] and Chen et al. [6] simplified the RCS model of the unmanned combat aerial vehicle (UCAV) into an ellipsoid model and designed a low-observable trajectory of UCAV. Zhao et al. [7] proposed a trajectory planning algorithm-A∗ algorithm that satisfies the computational speed and unmanned aerial vehicle (UAV) stealth characteristics in a multiradar threat environment.In the research of cooperative penetration trajectory, the focus is on the design of the terminal cooperative guidance law. The main research results include time cooperative guidance law [8–11], angle cooperative guidance law [12], and the impact time/angle cooperative guidance law [13–16]. But the premise is that the guidance law is designed based on the hypothesis of a small angle and the linearization of the kinematic equations. The shortcoming is that the flight path under the conditions of long distance, wide airspace, and speed range is prone to large errors. For the cooperative reentry of hypersonic vehicles, Fang et al. [17] designed a time-controlled reentry guidance law and proposed a two-layer cooperative reentry structure. Wang et al. [18] proposed a predictor-corrector cooperative guidance law based on altitude-velocity profile. But neither of them considered the real situation of radar detection threat with the aircraft.To overcome the shortcomings of the two aspects mentioned above, this study proposes a penetration trajectory optimization method and a cooperative penetration strategy of multiple aircraft, where stealth and cooperative penetrations are combined. The rest of the study is organized as follows: InMultiaircraft Cooperative Penetration Problem Model, the RCS characteristic model of the aircraft is established based on the physical optics method. A radar detection threat model is then built considering the impact of RCS. In The Aircraft Penetration Trajectory Optimization, the trajectory optimization framework with the dynamic model, constraint conditions, and optimal objectives is constructed. The hp-adaptive Radau pseudospectral method is applied to solve the optimal control problem, which is strongly nonlinear and satisfies multiconstraints. In The Cooperative Penetration Strategy of the Multiaircraft, a cooperative penetration strategy is proposed to solve the cooperative penetration problem of the multiaircraft. The impact time and angle constraints are given, and the virtual target point is introduced as the projected target. In Simulation Verification and Analysis, simulation and analysis of the single-aircraft and multiaircraft cooperative penetration trajectory optimizations are carried out.


## 2. Multiaircraft Cooperative Penetration Problem Model

### 2.1. Aircraft Dynamic Model

The unpowered gliding aircraft is taken as the object of study, under the research background that multiple aircraft reach the terminal guidance area cooperatively from the glide starting point. The three-degree-of-freedom dynamic equation for a single aircraft is established without considering the effects of the Earth’s rotation and flattening [19]:
(1)ṙ=Vsinθ,λ̇=Vcosθsinσrcosϕ,ϕ̇=Vcosθcosσr,V̇=−Dm−gsinθ,θ̇=LcosνmV+Vr−gVcosθ,σ̇=LsinνmVcosθ+Vrcosθsinσtanϕ.In this equation,r is the radial distance from the center of the earth to the aircraft and r=Re+h, where Re is earth’s average radius and h is the height of the aircraft from the local level; λ is the longitude; ϕ is the latitude; V is the velocity; θ is the flight path angle; σ is the heading angle; m is the mass of the aircraft; g is the acceleration due to gravity; ν is the bank angle; and the aerodynamic drag and lift forces are D and L, respectively, which can be defined as follows:
(2)L=ρV2ClS/2=ρ0e−βr−r0V2ClS/2,D=ρV2CdS/2=ρ0e−βr−r0V2CdS/2,where S is the reference area, ρ=ρ0e−βr−r0 is the atmospheric density, ρ0=1.225kg/m3, and β=1/7100. CD and CL are the drag and lift coefficients, respectively, which are functions of the angle of attack and Mach number.In the whole process of combat flight, the aircraft should also meet the following constraints:


#### 2.1.1. Terminal Constraints

For the initial state, the aircraft state variables are usually constant, and the terminal state constraints are primarily determined by the mission. If the aircraft is required to reach a specified target point, the position parameters, geocentric distancer, longitude λ, and latitude ϕ are constrained. The terminal constraints can be defined as follows:
(3)rtf=rf,ϕtf=ϕf,λtf=λf,where the subscript f denotes the terminal state.


#### 2.1.2. Control Constraints

The angle of attack and the bank angle are used as control variables. The control constraints can be expressed as follows:(4)αmin≤α≤αmax,νmin≤ν≤νmax,where αmin and αmax are the minimum and maximum values of the angle of attack, respectively. νmin and νmax are the minimum and maximum values of the bank angle, respectively.


#### 2.1.3. Equilibrium Constraints

In order to reduce the pressure of the aircraft control system, the aircraft adopts an improved quasi-equilibrium glide strategy. The equilibrium glide constraints can be expressed as follows:(5)ε=LcosνmV+Vr−gVcosθ,where ε is set to be a small negative number, ensuring the slow descent of the trajectory.


### 2.2. Aircraft RCS Characteristic Model

The aircraft RCS characteristic model is established considering the influence on the radar detection threat from aircraft radar cross-section (RCS). A plane symmetry aircraft is selected as the object of study, whose shape can be seen in Figure1.Figure 1
Shape of the aircraft.According to the shape and parameters of the given aircraft, the omnidirectional RCS data of the aircraft are obtained by fast calculation using the physical optics method (PO) [20], as shown in Figure 2.Figure 2
Omnidirectional RCS data of the aircraft.It can be seen from Figure2 that the raw data are plane symmetry, which is the same with the shape of the aircraft. It can also be seen from Figure 2 that the raw data obtained drastically fluctuate and have poor continuity, which is not conducive to the application of the subsequent algorithm. Therefore, the raw data is first smoothed by Gaussian filtering.The Gaussian filter is a low-pass weighted filter that is very effective for smoothing the data. According to the Gaussian filtering principle, the RCS value of a certain position of the aircraft is not only related to itself, but also is affected by the RCS in its adjacent area. The Gaussian filter function can be expressed as follows:(6)gx,y=12πσ12e−x2+y2/2σ12,where σ1 represents variance and x,y represents the raw data.The aircraft RCS data obtained by Gaussian filtering are shown in Figure3.Figure 3
RCS data obtained by Gaussian filtering.Second, the RCS surface is fitted using a polynomial. The selected polynomial can be expressed as follows:(7)σrcsθr,ϕr=∑i=0n∑j=0mpi,jθriϕrj,where pi,j is the fitting polynomial coefficient, θr is the elevation angle of the aircraft in the radar coordinate system, and ϕr is the azimuth angle in the radar coordinate system.Equation (7) indicates that the higher the degree of the polynomial, the more accurate the result. But the equation will become much more complex. Although higher order can be selected in the fitting equation, this paper takes m=n=5 to ensure a certain fitting equation accuracy with acceptable complexity. The obtained fitting result is shown in Figure 4.Figure 4
RCS data obtained by polynomial fitting.As depicted in Figure4, the fitted surface is not exactly the same as the original surface, but basically reflects the data characteristics of the original surface. The smaller values of the aircraft RCS are concentrated in the nearby intermediate region, where ϕr=90∘, the larger values are distributed around it. The fitting surface preserves the nearby ridge line, where θr=90∘. At the same time, the coefficients of the root mean square error (RMSE) and the square of correlation (R-square) obtained by polynomial fitting are 4.7401 and 0.9107, respectively. This indicates that the fitting accuracy is better than the results obtained in Figure 2.In order to calculate the RCS of an aircraft in a certain flight state, it is necessary to first calculate the elevation and azimuth angles at the current moment.Assuming that the positions of the aircraft and radar at the current moment areλ,ϕ,r and λR,ϕR,rR, respectively. The positions in the geocentric coordinate system E can be expressed as X,Y,Z=rcosϕcosλ,rcosϕsinλ,rsinϕ and XR,YR,ZR=rRcosϕRcosλR,rRcosϕRsinλR,rRsinϕR. In the geocentric coordinate system E, the unit vector of the radar position relative to the aircraft can be expressed as follows:
(8)x→E=XR−X,YR−Y,ZR−ZXR−X2+YR−Y2+YR−Y2.After coordinate transformation, the unit vector can be expressed in the radar coordinate system as follows:(9)x→R=REx→E.From this, parametersθr and ϕr corresponding to the RCS of the aircraft at the current moment can be calculated as follows:
(10)θr=cos−1x→R,3,ϕr=tan−1x→R,2x→R,1,where x→R,i represents the ith component of the unit vector x→R,i.


### 2.3. Radar Detection Threat Model

The instantaneous detection probability of the radar can be approximated as [21]
(11)Pt=1/1+c2R4/σrcsc1,where R is the distance between the aircraft and the radar. c1 and c2 are constants that are determined by the radar power, signal processing capabilities, frequency parameters, etc. They can be used to adjust the impact of distance and aircraft RCS on the radar detection probability. σrcs is the aircraft RCS. It can be seen from the formula that when the radar performance parameters are determined, the radar instantaneous detection probability is related to the aircraft RCS and the distance between the aircraft and the radar. Setting the constants c1=1.01and c2=1.25×10−24 in Equation (11). The relationship between the instantaneous detection probability of the radar, aircraft RCS, and the distance between the aircraft and the radar are given in the following figures.Figure5(a) shows the variation curve of the radar detection probability with the aircraft RCS when the distance between the aircraft and the radar is 1000 km. As can be seen from this figure, the detection probability increases with the increase of the aircraft RCS. And when the RCS is less than -20 dBsm, the detection probability is close to 0. When the RCS is greater than 20 dBsm, the detection probability is close to 1.Figure 5
Radar detection probability model curve.
(a)
Relationship between discovery probability and RCS whenR=1000km
(b)
Relationship between discovery probability and R whenσrcs=0dBsm
(c)
Detection probability surfaceFigure5(b) shows the variation curve of radar detection probability with the distance between the aircraft and the radar when the aircraft RCS is 0 dBsm. It can be seen from the figure that the detection probability gradually decreases with the increase of the detection distance. When the detection distance is less than 300 km, the detection probability is close to 1. When the detection distance is greater than 1000 km, the detection probability slowly decreases. And when R=3000 km, the detection probability is only 0.9343%. Figure 5(c) shows the radar detection probability surface. It can be seen from this figure that there exists such a region. When the distance is close, the detection probability is close to 1, even if the RCS is very small. When the RCS is large, the detection probability is close to 1, even if the aircraft is far away from the radar. Therefore, if the aircraft is kept at a certain distance from the radar and the aircraft presents a low RCS in an appropriate attitude, it can effectively avoid the threat of radar detection.The radar detection threat to the aircraft during the entire combat flight timet0,tf can be expressed as follows:
(12)PT=∫t0tfPtdt.Assuming that the combat formation consists ofN aircrafts, the total radar detection threat to the entire combat formation can be defined as follows:
(13)PT=∑i=1NwiPT,i,where wi is the threat weight of the ith aircraft and PT,iis the radar detection threat of the ith aircraft.The different models established in this section will be used in the following sections.


## 2.1. Aircraft Dynamic Model

The unpowered gliding aircraft is taken as the object of study, under the research background that multiple aircraft reach the terminal guidance area cooperatively from the glide starting point. The three-degree-of-freedom dynamic equation for a single aircraft is established without considering the effects of the Earth’s rotation and flattening [19]:
(1)ṙ=Vsinθ,λ̇=Vcosθsinσrcosϕ,ϕ̇=Vcosθcosσr,V̇=−Dm−gsinθ,θ̇=LcosνmV+Vr−gVcosθ,σ̇=LsinνmVcosθ+Vrcosθsinσtanϕ.In this equation,r is the radial distance from the center of the earth to the aircraft and r=Re+h, where Re is earth’s average radius and h is the height of the aircraft from the local level; λ is the longitude; ϕ is the latitude; V is the velocity; θ is the flight path angle; σ is the heading angle; m is the mass of the aircraft; g is the acceleration due to gravity; ν is the bank angle; and the aerodynamic drag and lift forces are D and L, respectively, which can be defined as follows:
(2)L=ρV2ClS/2=ρ0e−βr−r0V2ClS/2,D=ρV2CdS/2=ρ0e−βr−r0V2CdS/2,where S is the reference area, ρ=ρ0e−βr−r0 is the atmospheric density, ρ0=1.225kg/m3, and β=1/7100. CD and CL are the drag and lift coefficients, respectively, which are functions of the angle of attack and Mach number.In the whole process of combat flight, the aircraft should also meet the following constraints:


### 2.1.1. Terminal Constraints

For the initial state, the aircraft state variables are usually constant, and the terminal state constraints are primarily determined by the mission. If the aircraft is required to reach a specified target point, the position parameters, geocentric distancer, longitude λ, and latitude ϕ are constrained. The terminal constraints can be defined as follows:
(3)rtf=rf,ϕtf=ϕf,λtf=λf,where the subscript f denotes the terminal state.


### 2.1.2. Control Constraints

The angle of attack and the bank angle are used as control variables. The control constraints can be expressed as follows:(4)αmin≤α≤αmax,νmin≤ν≤νmax,where αmin and αmax are the minimum and maximum values of the angle of attack, respectively. νmin and νmax are the minimum and maximum values of the bank angle, respectively.


### 2.1.3. Equilibrium Constraints

In order to reduce the pressure of the aircraft control system, the aircraft adopts an improved quasi-equilibrium glide strategy. The equilibrium glide constraints can be expressed as follows:(5)ε=LcosνmV+Vr−gVcosθ,where ε is set to be a small negative number, ensuring the slow descent of the trajectory.


## 2.1.1. Terminal Constraints

For the initial state, the aircraft state variables are usually constant, and the terminal state constraints are primarily determined by the mission. If the aircraft is required to reach a specified target point, the position parameters, geocentric distancer, longitude λ, and latitude ϕ are constrained. The terminal constraints can be defined as follows:
(3)rtf=rf,ϕtf=ϕf,λtf=λf,where the subscript f denotes the terminal state.


## 2.1.2. Control Constraints

The angle of attack and the bank angle are used as control variables. The control constraints can be expressed as follows:(4)αmin≤α≤αmax,νmin≤ν≤νmax,where αmin and αmax are the minimum and maximum values of the angle of attack, respectively. νmin and νmax are the minimum and maximum values of the bank angle, respectively.


## 2.1.3. Equilibrium Constraints

In order to reduce the pressure of the aircraft control system, the aircraft adopts an improved quasi-equilibrium glide strategy. The equilibrium glide constraints can be expressed as follows:(5)ε=LcosνmV+Vr−gVcosθ,where ε is set to be a small negative number, ensuring the slow descent of the trajectory.


## 2.2. Aircraft RCS Characteristic Model

The aircraft RCS characteristic model is established considering the influence on the radar detection threat from aircraft radar cross-section (RCS). A plane symmetry aircraft is selected as the object of study, whose shape can be seen in Figure1.Figure 1
Shape of the aircraft.According to the shape and parameters of the given aircraft, the omnidirectional RCS data of the aircraft are obtained by fast calculation using the physical optics method (PO) [20], as shown in Figure 2.Figure 2
Omnidirectional RCS data of the aircraft.It can be seen from Figure2 that the raw data are plane symmetry, which is the same with the shape of the aircraft. It can also be seen from Figure 2 that the raw data obtained drastically fluctuate and have poor continuity, which is not conducive to the application of the subsequent algorithm. Therefore, the raw data is first smoothed by Gaussian filtering.The Gaussian filter is a low-pass weighted filter that is very effective for smoothing the data. According to the Gaussian filtering principle, the RCS value of a certain position of the aircraft is not only related to itself, but also is affected by the RCS in its adjacent area. The Gaussian filter function can be expressed as follows:(6)gx,y=12πσ12e−x2+y2/2σ12,where σ1 represents variance and x,y represents the raw data.The aircraft RCS data obtained by Gaussian filtering are shown in Figure3.Figure 3
RCS data obtained by Gaussian filtering.Second, the RCS surface is fitted using a polynomial. The selected polynomial can be expressed as follows:(7)σrcsθr,ϕr=∑i=0n∑j=0mpi,jθriϕrj,where pi,j is the fitting polynomial coefficient, θr is the elevation angle of the aircraft in the radar coordinate system, and ϕr is the azimuth angle in the radar coordinate system.Equation (7) indicates that the higher the degree of the polynomial, the more accurate the result. But the equation will become much more complex. Although higher order can be selected in the fitting equation, this paper takes m=n=5 to ensure a certain fitting equation accuracy with acceptable complexity. The obtained fitting result is shown in Figure 4.Figure 4
RCS data obtained by polynomial fitting.As depicted in Figure4, the fitted surface is not exactly the same as the original surface, but basically reflects the data characteristics of the original surface. The smaller values of the aircraft RCS are concentrated in the nearby intermediate region, where ϕr=90∘, the larger values are distributed around it. The fitting surface preserves the nearby ridge line, where θr=90∘. At the same time, the coefficients of the root mean square error (RMSE) and the square of correlation (R-square) obtained by polynomial fitting are 4.7401 and 0.9107, respectively. This indicates that the fitting accuracy is better than the results obtained in Figure 2.In order to calculate the RCS of an aircraft in a certain flight state, it is necessary to first calculate the elevation and azimuth angles at the current moment.Assuming that the positions of the aircraft and radar at the current moment areλ,ϕ,r and λR,ϕR,rR, respectively. The positions in the geocentric coordinate system E can be expressed as X,Y,Z=rcosϕcosλ,rcosϕsinλ,rsinϕ and XR,YR,ZR=rRcosϕRcosλR,rRcosϕRsinλR,rRsinϕR. In the geocentric coordinate system E, the unit vector of the radar position relative to the aircraft can be expressed as follows:
(8)x→E=XR−X,YR−Y,ZR−ZXR−X2+YR−Y2+YR−Y2.After coordinate transformation, the unit vector can be expressed in the radar coordinate system as follows:(9)x→R=REx→E.From this, parametersθr and ϕr corresponding to the RCS of the aircraft at the current moment can be calculated as follows:
(10)θr=cos−1x→R,3,ϕr=tan−1x→R,2x→R,1,where x→R,i represents the ith component of the unit vector x→R,i.


## 2.3. Radar Detection Threat Model

The instantaneous detection probability of the radar can be approximated as [21]
(11)Pt=1/1+c2R4/σrcsc1,where R is the distance between the aircraft and the radar. c1 and c2 are constants that are determined by the radar power, signal processing capabilities, frequency parameters, etc. They can be used to adjust the impact of distance and aircraft RCS on the radar detection probability. σrcs is the aircraft RCS. It can be seen from the formula that when the radar performance parameters are determined, the radar instantaneous detection probability is related to the aircraft RCS and the distance between the aircraft and the radar. Setting the constants c1=1.01and c2=1.25×10−24 in Equation (11). The relationship between the instantaneous detection probability of the radar, aircraft RCS, and the distance between the aircraft and the radar are given in the following figures.Figure5(a) shows the variation curve of the radar detection probability with the aircraft RCS when the distance between the aircraft and the radar is 1000 km. As can be seen from this figure, the detection probability increases with the increase of the aircraft RCS. And when the RCS is less than -20 dBsm, the detection probability is close to 0. When the RCS is greater than 20 dBsm, the detection probability is close to 1.Figure 5
Radar detection probability model curve.
(a)
Relationship between discovery probability and RCS whenR=1000km
(b)
Relationship between discovery probability and R whenσrcs=0dBsm
(c)
Detection probability surfaceFigure5(b) shows the variation curve of radar detection probability with the distance between the aircraft and the radar when the aircraft RCS is 0 dBsm. It can be seen from the figure that the detection probability gradually decreases with the increase of the detection distance. When the detection distance is less than 300 km, the detection probability is close to 1. When the detection distance is greater than 1000 km, the detection probability slowly decreases. And when R=3000 km, the detection probability is only 0.9343%. Figure 5(c) shows the radar detection probability surface. It can be seen from this figure that there exists such a region. When the distance is close, the detection probability is close to 1, even if the RCS is very small. When the RCS is large, the detection probability is close to 1, even if the aircraft is far away from the radar. Therefore, if the aircraft is kept at a certain distance from the radar and the aircraft presents a low RCS in an appropriate attitude, it can effectively avoid the threat of radar detection.The radar detection threat to the aircraft during the entire combat flight timet0,tf can be expressed as follows:
(12)PT=∫t0tfPtdt.Assuming that the combat formation consists ofN aircrafts, the total radar detection threat to the entire combat formation can be defined as follows:
(13)PT=∑i=1NwiPT,i,where wi is the threat weight of the ith aircraft and PT,iis the radar detection threat of the ith aircraft.The different models established in this section will be used in the following sections.


## 3. The Aircraft Penetration Trajectory Optimization

The solution framework of penetration trajectory optimization problem proposed in this study is shown in Figure6.Figure 6
Framework for penetration trajectory optimization problem.The framework for penetration trajectory optimization problem is formed through the following steps:Firstly, given the research object, the RCS raw data of the aircraft can be calculated by physical optics method. The aircraft RCS characteristic model is then established through Gaussian filtering and polynomial fitting processing. Based on the aircraft dynamic model and radar detection probability model, the penetration trajectory optimization design model is established, which is introduced inMultiaircraft Cooperative Penetration Problem Model.Secondly, the aircraft dynamic model used as the dynamic constraint, the objective function corresponding to the radar detection probability, the terminal constraint, and the control constraint jointly construct a framework for the optimal control problem.Then, based on the hp-adaptive Radau pseudospectral method, the optimal control problem is transformed into a nonlinear programming problem.Finally, the flight trajectory under multiconstraint conditions is obtained, and the flight time, flight position, and flight altitude data are obtained as the output.In recent years, the pseudospectral method has been widely used in the field of aircraft trajectory optimization [22–24] due to its characteristics such as large convergence radius, insensitivity to initial value, and fast convergence rate.The hp-adaptive Radau pseudospectral method [25] combines the advantages of p method and h method in that both the number of mesh intervals (the sparse calculation of h method) and the degree of the approximating polynomial within each mesh interval (rapid convergence of the p method) can be varied, which can be used to solve the complex optimal control problem. The basic principle of the Radau pseudospectral method is to discretize the state and control variables in each mesh interval on Legendre-Gauss-Radau (LGR) points and take the discrete points as nodes to create a Lagrange interpolation polynomial to fit the approximate state and control variables. By deriving the state variable obtained by the global interpolation polynomial, the derivative of the state variable to time is approximated. Then, the dynamical differential equation constraints of the system can be converted into algebraic constraints, and the integral in the performance index is calculated by Radau integral. Since the terminal state of the global current mesh interval is the initial state of the next mesh interval, the integration process of the terminal state in the global Radau pseudospectral method is avoided. Through the above series of steps, the original optimal control problem can be transformed into an NLP problem.The main solution steps are shown as follows.Step 1.
Time variable transformation.
The time intervalt0,tf is divided into a mesh consisting of K mesh intervalstk−1,tk, k=1,2,⋯,K. The variable t∈tk−1,tk can be transformed into t∈−1,1 by the following transformation:
(14)τ=2t−tk+tk−1tk−tk−1.Step 2.
Discretize the state and control variables.
The number of collocations in thekth interval is Nk, which can be expressed asτ1k,⋯,τNkk. The Nkth degree Legendre polynomial is expressed as follows:
(15)Likτ=∏j=1j≠iNk+1τ−τjτi−τji=1,⋯,Nk+1.
Therefore, the statexkτ is approximated as follows:
(16)xkτ≈Xkτ=∑j=0NkXjkLjkτ.
Similarly, the controlukτ is approximated as follows:
(17)ukτ≈Ukτ=∑j=0NkUjkLjkτ.Step 3.
Discretize dynamical differential equation.
Differentiating the approximation state variable:(18)dXkτdt=∑j=0NkXjkL̇jkτ.
Then,(19)∑j=0NkXjkDijk−tk−tk−12fik=0i=1,⋯,Nk,k=1,⋯,K,where Dijk=L̇jkτik∈Nk×Nk+1i=1,⋯,Nk,j=1,⋯,Nk+1 is the differential approximation matrix in the kth interval.Step 4.
Discretize the path constraints and boundary conditions.
The path constraints can be approximated as follows:(20)CkXik,Uik≤0i=1,⋯,Nk,k=1,⋯,K.
The boundary condition can be approximated as follows:(21)ϕx11,t0,xNk+1K,tK=0.Step 5.
Discretize the objective function.
The optimal problem in this paper is the continuous Bolza problem. And the Bolza objective function can be expressed as the following equation after the above steps:(22)J=ΦX0,t0,Xf,tf+tf−t02∑k=1KωkgXτk,Uτk,τk;t0,tf.
For a single aircraft penetration trajectory optimization in this paper, the objective function can be expressed by Equation (12), which subjects to dynamic (Equation (1)), path (Equation (5)), and boundary constraints (Equations (3) and (4)).


## 4. The Cooperative Penetration Strategy of the Multiaircraft

With the continuous development in air defense and antimissile technology, the modern defense system is improving every day. The difficulty in the single aircraft penetration is increasing. Therefore, it is necessary to adopt a cooperative penetration strategy to increase the penetration capability. The cooperative in this paper refers to the cooperative of the impact time and angle.


### 4.1. The Constraints of Cooperative Impact Time/Angle

For the selection of the cooperative impact time, if the longest flight time of the aircraft is taken, it is necessary to calculate the flight time of each aircraft offline in advance, which is not conducive to online application. For this purpose, this paper implements time coordination by extending state variables. As shown in Equation (1), the motion equation of a single aircraft contains six state and two control variables. The left side of the equation is the differential of state variables to time. Assuming that a combat formation consists of N aircrafts, the number of states and control variables in the equations of motion are 6∗N and 2∗N, respectively. Then, the motion equations are solved as a whole, and time cooperative can be realized accordingly. The motion equation of the multiaircraft can be expressed as follows:
(23)Ẋ1=FX1,Ẋ2=FX2,⋯ẊN=FXN,,where Xi=ri,λi,ϕi,Vi,θi,σi is the state variable of the ith aircraft.The cooperative impact angle is defined as the angle between the terminal velocity of the aircraft and the north direction. The new terminal constraints can be expressed as follows:(24)σtf=σf,co,where σf,co is the specific terminal cooperative impact angle. The cooperative impact angles of different aircrafts can be given the same value or be set different according to typical combat mission requirements.


### 4.2. Virtual Target Point Setting

In order to facilitate the multiaircraft to strike the target from different places at different impact angles, a cylinder to describe the boundary of the terminal guidance area is defined, whose center isλT,ϕT. It has a radius of Rter and an height of hter. The upper surface circle of the cylinder is defined as the midterminal guidance transition area. It is assumed that when the multiple aircraft reach the specified position on the arc at the end of middle guidance, the effective multidirection penetration strikes are triggered and can then be carried out. For this purpose, a virtual target point is introduced, and its position is set as follows:
(25)λdT=λT+Rtersinψ/Re,ϕdT=ϕT+Rtercosψ/Re,where λdT,ϕdT is the longitude and latitude of the virtual target, and ψ is the angle between the north direction and the projection of the line connecting the virtual target point and the target point on the ground. The value range of this angle is 0o,360o.In the rest of the paper, the virtual target is calculated first to meet the requirements of the combat mission for each of the three targets.


## 4.1. The Constraints of Cooperative Impact Time/Angle

For the selection of the cooperative impact time, if the longest flight time of the aircraft is taken, it is necessary to calculate the flight time of each aircraft offline in advance, which is not conducive to online application. For this purpose, this paper implements time coordination by extending state variables. As shown in Equation (1), the motion equation of a single aircraft contains six state and two control variables. The left side of the equation is the differential of state variables to time. Assuming that a combat formation consists of N aircrafts, the number of states and control variables in the equations of motion are 6∗N and 2∗N, respectively. Then, the motion equations are solved as a whole, and time cooperative can be realized accordingly. The motion equation of the multiaircraft can be expressed as follows:
(23)Ẋ1=FX1,Ẋ2=FX2,⋯ẊN=FXN,,where Xi=ri,λi,ϕi,Vi,θi,σi is the state variable of the ith aircraft.The cooperative impact angle is defined as the angle between the terminal velocity of the aircraft and the north direction. The new terminal constraints can be expressed as follows:(24)σtf=σf,co,where σf,co is the specific terminal cooperative impact angle. The cooperative impact angles of different aircrafts can be given the same value or be set different according to typical combat mission requirements.


## 4.2. Virtual Target Point Setting

In order to facilitate the multiaircraft to strike the target from different places at different impact angles, a cylinder to describe the boundary of the terminal guidance area is defined, whose center isλT,ϕT. It has a radius of Rter and an height of hter. The upper surface circle of the cylinder is defined as the midterminal guidance transition area. It is assumed that when the multiple aircraft reach the specified position on the arc at the end of middle guidance, the effective multidirection penetration strikes are triggered and can then be carried out. For this purpose, a virtual target point is introduced, and its position is set as follows:
(25)λdT=λT+Rtersinψ/Re,ϕdT=ϕT+Rtercosψ/Re,where λdT,ϕdT is the longitude and latitude of the virtual target, and ψ is the angle between the north direction and the projection of the line connecting the virtual target point and the target point on the ground. The value range of this angle is 0o,360o.In the rest of the paper, the virtual target is calculated first to meet the requirements of the combat mission for each of the three targets.


## 5. Simulation Verification and Analysis

In order to fully verify the effectiveness and applicability of the method presented in this paper, the following combat scenario is proposed. At the initial moment, three aircrafts with the same performance parameters are launched at different points and are required to strike on the same target from different directions approaching the midterminal guidance transition area. During the flight, they are threatened by a noncooperative radar at the position (35°E, 2°N). Parameters of the terminal guidance area are set asRter=50km and hter=30km. The simulation parameters of the combat mission of the three aircrafts are given in Table 1.Table 1
Simulation parameter setting.
Aircraft
Initial position
Target position
Height (km)
Velocity (m/s)
Flight path angle (°)
Heading angle (°)
Impact angle (°)
A
(0° E, 5° N)
(78° E, 0° N)
60
7000
0
90
180
B
(0° E, -5°N)
(78° E, 0° N)
60
7000
0
90
0
C
(0° E, 0° N)
(78° E, 0° N)
60
7000
0
90
90


### 5.1. Single Aircraft Penetration

Firstly, the penetration trajectory of the aircraft C is analyzed to verify the penetration trajectory optimization method. It is noted that the deployment position of the radar is set as (35° E, 0° N) in the following two examples.Example 1.
In order to strike the time-sensitive target, the optimal trajectory from the initial to target point is designed, whose objective function is the minimum flight time. The trajectory satisfies all the constraints mentioned above in the flight process and is denoted as MT.Example 2.
In order to achieve effective stealth and reduce the threat of radar detection, an optimal trajectory from the initial to target point is designed with the minimum radar detection probability integral as the objective function during the entire flight process. The trajectory satisfies all the constraints mentioned above in the flight process and is denoted as MP.Simulation results of these two examples are shown and compared with each other in Figure7.Figure 7
Simulation results of MT and MP trajectories.
(a)
Longitude-latitude
(b)
Time-detection probability
(c)
Time-RCS
(d)
Attitude angle trajectory
(e)
Time-angle of attack
(f)
Time-bank angleIt can be seen from Figure7(a) that when the radar detection threat is not considered, the projection of trajectory MT on the ground directly passes through the radar detection center, while trajectory MP passes close to the radar detection center to reduce the radar detection threat to the aircraft. The maneuvering is realized with a certain evasive action. Both the trajectories (MT and MP) reach the virtual target point with a specified attack angle.It can also be seen from Figure7(b) that since the trajectory MT directly passes through the radar detection center, the radar detection probability has a peak value of 0.9921 at 648.6 s. In comparison to the results of trajectory MT, the radar detection probability of the trajectory MP is smaller during the entire process, and its peak value is 0.003145 at 564.4 s.Moreover, it can be found from Figure7(c) that the RCS of the trajectory MT is higher than that of trajectory MP. This is because when considering the radar threat, the RCS is reduced by the flight attitude adjustment to reduce the probability of radar detection. The corresponding flight attitude angle trajectory is demonstrated in Figure 7(d). At the same time, it should be noted that from 1630 s to the end, there is a sudden change in the RCS of the trajectory MP. This is because, after the aircraft crosses the radar detection area, the attitude of the aircraft is adjusted to meet the specified attack angle, which causes the change in RCS. The corresponding control variable curves are shown in Figures 7(e) and 7(f). However, since the aircraft is far away from the radar currently, the RCS has little influence on the radar detection probability.To sum up, when the radar detection threat is not considered, the aircraft can quickly pass through the radar detection area without lateral maneuvering and reach the target point in the shortest time to meet the specified attack angle. When the radar detection threat is considered, the aircraft can reduce its own RCS and perform lateral maneuvering by adjusting its attitude. Besides, the aircraft is able to reach the target point according to the specified attack angle while satisfying the minimum threat.To explore the possibility of cooperative strike formatted by multiple aircraft, a scenario of independent penetrations from three aircrafts is established.The radar position is readjusted to (35° E, 2° N) and the single penetration trajectory of three aircrafts A, B, and C is analyzed and compared. All the radar detection threat in the whole combat process is taken as the optimization objective. The simulation results are shown in Figures8 and 9.Figure 8
Longitude-latitude.Figure 9
Time-azimuth angle.As can be seen from Figures8 and 9, the flight time of the three aircrafts are 1766, 1952, and 2212 s, respectively. In order to minimize the threat of radar detection during the entire process, the flight around strategy is adopted when passing through the radar detection area. And the aircraft can reach the target point according to the specified attack angle at the terminal moment, which lays a foundation for further cooperative strike.


### 5.2. Multiple Aircraft Cooperative Penetration

It can be found from Figures8 and 9 that the minimum flight time interval is 186 s when single aircraft sequentially strikes the target. When the time interval is longer, the penetration is harder to be conducted. The multiaircraft cooperative penetration strategy is adopted to improve the overall combat effectiveness by locating multiaircraft at the boundary of the terminal guidance area at exactly the same time. Two different scenarios are considered below.When only aircrafts A and B’s cooperative penetration is considered, the minimum of total radar detection threat received by A and B is taken as the optimization target. The objective function can be expressed as follows:(26)J=PT=wAPT,A+wBPT,B.The weights of aircrafts A and B are 0.5 and 0.5, respectively. The parameters of the combat missions of the two aircrafts are shown in Table1. The simulation results are shown in Figures 9 and 10.Figure 10
Longitude-latitude.It can be seen from Figures10 and 11 that when aircrafts A and B adopt a cooperative strategy, they can reach their respective virtual target points at 1783 s. The heading angles of A and B are 180° and 0° , respectively, which meet the requirements of the cooperative impact angles. Compared with Figures 8 and 9, the flight time of aircraft A increases, and the flight time of aircraft B decreases to satisfy the cooperative impact time.Figure 11
Time-azimuth angle.When cooperative penetration of the three aircrafts A, B, and C is adopted, the minimum of total radar detection threat received by aircrafts A, B, and C is taken as the optimization target. The objective function can be expressed as follows:(27)J=PT=wAPT,A+wBPT,B+wCPT,C.The weights of aircrafts A, B, and C are 1/3, 1/3, and 1/3, respectively. The parameters of the combat mission for the three aircrafts are shown in Table1. The simulation results are shown in Figures 12 and 13.Figure 12
Longitude-latitude.Figure 13
Time-azimuth angle.As can be seen from Figures12 and 13, when aircrafts A, B, and C’s formation adopt the cooperative penetration strategy, they reach their respective virtual target points at 1880 s. The terminal heading angles of A, B, and C are 180°, 0°, and 90°, respectively, which meet the requirements of the cooperative impact angles.Comparing Figures8–13, it is evident that the aircrafts can reach their respective virtual target points from designated cooperative impact angles when the cooperative penetration strategy is adopted, which lays a foundation for further cooperative strike in the terminal phase. It also should be noted that the maneuvering is realized with a certain evasive action for every aircraft to reduce the total radar detection threat.


## 5.1. Single Aircraft Penetration

Firstly, the penetration trajectory of the aircraft C is analyzed to verify the penetration trajectory optimization method. It is noted that the deployment position of the radar is set as (35° E, 0° N) in the following two examples.Example 1.
In order to strike the time-sensitive target, the optimal trajectory from the initial to target point is designed, whose objective function is the minimum flight time. The trajectory satisfies all the constraints mentioned above in the flight process and is denoted as MT.Example 2.
In order to achieve effective stealth and reduce the threat of radar detection, an optimal trajectory from the initial to target point is designed with the minimum radar detection probability integral as the objective function during the entire flight process. The trajectory satisfies all the constraints mentioned above in the flight process and is denoted as MP.Simulation results of these two examples are shown and compared with each other in Figure7.Figure 7
Simulation results of MT and MP trajectories.
(a)
Longitude-latitude
(b)
Time-detection probability
(c)
Time-RCS
(d)
Attitude angle trajectory
(e)
Time-angle of attack
(f)
Time-bank angleIt can be seen from Figure7(a) that when the radar detection threat is not considered, the projection of trajectory MT on the ground directly passes through the radar detection center, while trajectory MP passes close to the radar detection center to reduce the radar detection threat to the aircraft. The maneuvering is realized with a certain evasive action. Both the trajectories (MT and MP) reach the virtual target point with a specified attack angle.It can also be seen from Figure7(b) that since the trajectory MT directly passes through the radar detection center, the radar detection probability has a peak value of 0.9921 at 648.6 s. In comparison to the results of trajectory MT, the radar detection probability of the trajectory MP is smaller during the entire process, and its peak value is 0.003145 at 564.4 s.Moreover, it can be found from Figure7(c) that the RCS of the trajectory MT is higher than that of trajectory MP. This is because when considering the radar threat, the RCS is reduced by the flight attitude adjustment to reduce the probability of radar detection. The corresponding flight attitude angle trajectory is demonstrated in Figure 7(d). At the same time, it should be noted that from 1630 s to the end, there is a sudden change in the RCS of the trajectory MP. This is because, after the aircraft crosses the radar detection area, the attitude of the aircraft is adjusted to meet the specified attack angle, which causes the change in RCS. The corresponding control variable curves are shown in Figures 7(e) and 7(f). However, since the aircraft is far away from the radar currently, the RCS has little influence on the radar detection probability.To sum up, when the radar detection threat is not considered, the aircraft can quickly pass through the radar detection area without lateral maneuvering and reach the target point in the shortest time to meet the specified attack angle. When the radar detection threat is considered, the aircraft can reduce its own RCS and perform lateral maneuvering by adjusting its attitude. Besides, the aircraft is able to reach the target point according to the specified attack angle while satisfying the minimum threat.To explore the possibility of cooperative strike formatted by multiple aircraft, a scenario of independent penetrations from three aircrafts is established.The radar position is readjusted to (35° E, 2° N) and the single penetration trajectory of three aircrafts A, B, and C is analyzed and compared. All the radar detection threat in the whole combat process is taken as the optimization objective. The simulation results are shown in Figures8 and 9.Figure 8
Longitude-latitude.Figure 9
Time-azimuth angle.As can be seen from Figures8 and 9, the flight time of the three aircrafts are 1766, 1952, and 2212 s, respectively. In order to minimize the threat of radar detection during the entire process, the flight around strategy is adopted when passing through the radar detection area. And the aircraft can reach the target point according to the specified attack angle at the terminal moment, which lays a foundation for further cooperative strike.


## 5.2. Multiple Aircraft Cooperative Penetration

It can be found from Figures8 and 9 that the minimum flight time interval is 186 s when single aircraft sequentially strikes the target. When the time interval is longer, the penetration is harder to be conducted. The multiaircraft cooperative penetration strategy is adopted to improve the overall combat effectiveness by locating multiaircraft at the boundary of the terminal guidance area at exactly the same time. Two different scenarios are considered below.When only aircrafts A and B’s cooperative penetration is considered, the minimum of total radar detection threat received by A and B is taken as the optimization target. The objective function can be expressed as follows:(26)J=PT=wAPT,A+wBPT,B.The weights of aircrafts A and B are 0.5 and 0.5, respectively. The parameters of the combat missions of the two aircrafts are shown in Table1. The simulation results are shown in Figures 9 and 10.Figure 10
Longitude-latitude.It can be seen from Figures10 and 11 that when aircrafts A and B adopt a cooperative strategy, they can reach their respective virtual target points at 1783 s. The heading angles of A and B are 180° and 0° , respectively, which meet the requirements of the cooperative impact angles. Compared with Figures 8 and 9, the flight time of aircraft A increases, and the flight time of aircraft B decreases to satisfy the cooperative impact time.Figure 11
Time-azimuth angle.When cooperative penetration of the three aircrafts A, B, and C is adopted, the minimum of total radar detection threat received by aircrafts A, B, and C is taken as the optimization target. The objective function can be expressed as follows:(27)J=PT=wAPT,A+wBPT,B+wCPT,C.The weights of aircrafts A, B, and C are 1/3, 1/3, and 1/3, respectively. The parameters of the combat mission for the three aircrafts are shown in Table1. The simulation results are shown in Figures 12 and 13.Figure 12
Longitude-latitude.Figure 13
Time-azimuth angle.As can be seen from Figures12 and 13, when aircrafts A, B, and C’s formation adopt the cooperative penetration strategy, they reach their respective virtual target points at 1880 s. The terminal heading angles of A, B, and C are 180°, 0°, and 90°, respectively, which meet the requirements of the cooperative impact angles.Comparing Figures8–13, it is evident that the aircrafts can reach their respective virtual target points from designated cooperative impact angles when the cooperative penetration strategy is adopted, which lays a foundation for further cooperative strike in the terminal phase. It also should be noted that the maneuvering is realized with a certain evasive action for every aircraft to reduce the total radar detection threat.


## 6. Conclusions

In order to cope with the severe challenges brought about by the perfect modern defense system and to improve the aircraft penetration capability, this study proposes a penetration trajectory optimization method and a cooperative penetration strategy of multiple aircraft, considering the practical need to reduce the radar detection threat. The multiaircraft cooperative penetration problem models are established. A trajectory optimization framework is formed, which includes the RCS characteristic model and the radar detection threat model. Moreover, the hp-adaptive Radau pseudospectral method is proposed to solve the optimal control problem. In order to improve the penetration capability, the cooperative penetration strategy of multiaircraft is designed. The method and strategy are simulated and verified. Simulation results show that (1) when a single aircraft penetrates, it can adjust its attitude to reduce its own RCS and implement maneuvering. Under the condition of meeting the minimum threat of radar detection, it can carry out penetration strike according to the designated combat mission; (2) the multiaircraft can conduct penetration attack at a cooperative impact time and angle according to the cooperative penetration strategy with the condition of meeting the minimum radar detection threat of the entire formation. Thus, the method and strategy can be applicable to many-to-one engagement scenarios of multiple aircraft, and the penetration probability can be increased.



---

*Source: 8490531-2020-01-25.xml*