# Receding Horizon Trajectory Optimization with Terminal Impact Specifications

Authors: Limin Zhang; Mingwei Sun; Zengqiang Chen; Zenghui Wang; Yongkun Wang
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

The trajectory optimization problem subject to terminal impact time and angle specifications can be reformulated as a nonlinear programming problem using the Gauss pseudospectral method. The cost function of the trajectory optimization problem is modified to reduce the terminal control energy. A receding horizon optimization strategy is implemented to reject the errors caused by the motion of a surface target. Several simulations were performed to validate the proposed method via the C programming language. The simulation results demonstrate the effectiveness of the proposed algorithm and that the real-time requirement can be easily achieved if the C programming language is used to realize it.


---

## Body

## 1. Introduction

To enhance the penetration probability of a cruise missile, controlling both the impact angle and the impact time is of special significance in practice. In the past three decades, several guidance laws with terminal impact angle specifications were proposed within different frameworks [1–11]. It should be noted that the impact time is a crucial factor in the salvo attack against a highly valued surface target equipped with advanced air defense installations such as close-in weapon systems (CIWS). However the impact time specification is difficult to achieve and there were few reports in this field. Jeon et al. [7] proposed a 2-dimensional salvo attack strategy with terminal impact angle and time specifications, which was obtained by a linearized model based optimal control. An analytical guidance law was derived explicitly for this strategy. There are two control loops in this approach: the outer and inner loops are designed to regulate the time-to-go and the angle, respectively. When the time-to-go prediction is larger than the specified time, a trajectory curve will be planned online to consume the excessive time; otherwise, a comparatively straight line will be generated to meet the impact time requirement. Despite its simplicity, two severe pitfalls accompany this method: (
1
) the small heading angle assumption could cause large terminal errors when the practical angle is large and (
2
) the process constraints with the inequality forms, for example, the familiar no-fly zone constraints, are unable to be incorporated. To cope with the no-fly zone constraints, many trajectory planning methodologies based on A
* algorithm [12–15] and intelligent optimization algorithms, such as genetic algorithm (GA) [16–18], ant colony optimization (ACO) [19, 20], and particle swarm optimization (PSO) [21–23], have been extensively investigated. Nevertheless, these approaches can only generate an open-loop optimal trajectory offline due to the computational complexity and it is not applicable to use them to impact a moving surface target. Therefore, it is urgent to solve the trajectory optimization problem with terminal impact angle and time specifications as well as process inequality constraints simultaneously in real-time.Recently, the Gauss pseudospectral method (GPM) has attracted a wide attention within aerospace industry, especially in the fields of guidance and trajectory design [24–28]. GPM is a direct optimal control solving approach for the general nonlinear systems with various constraints. Nowadays, a mature GPM solving software GPOPS (Gauss Pseudospectral OPtimization Software) is available, which makes GPM a technology to be more than an art. In this paper, the GPM is employed in a receding horizon way to obtain an adaptive trajectory to impact a moving surface target, which satisfies the terminal and process constraints and can be realized online. In the design, a time-to-go dependent penalty function is utilized to ensure a smooth trajectory at the terminal end and reduce the miss distance.The remaining parts of this paper are organized as follows. The problem is formulated in Section2. Section 3 gives a brief introduction about GPM. The receding horizon strategy to implement GPM is presented in Section 4. In Section 5, the simulation results are provided and analyzed. The concluding remarks are summarized in Section 6.


## 2. Problem Formulation

### 2.1. Kinematics Equations

The missile-target relative kinematics is shown in Figure1. Without loss of generality, it is assumed that the velocity of the missile, V, is a constant. The earth coordinate system is chosen as X
O
Y. The current states of the missile are the position (
X
m
,
Y
m
) and the heading angle γ
m. The control variable for the missile is the normal acceleration, A, being perpendicular to V, which can change the velocity direction. The corresponding states of the target are the position (
X
t
,
Y
t
) and the heading angle γ
t.Figure 1
Missile-target relative kinematics.According to Figure1, the mass-point dynamics of the missile and the target are [9, 29]
(1)
d
X
m
d
t
=
V
cos
⁡
(
γ
m
)
,
d
Y
m
d
t
=
V
sin
(
γ
m
)
,
d
γ
m
d
t
=
u
m
=
A
V
,
d
X
t
d
t
=
V
t
cos
⁡
⁡
(
γ
t
)
,
d
Y
t
d
t
=
V
t
sin
(
γ
t
)
.


### 2.2. Optimal Trajectory Requirement

#### 2.2.1. Terminal Constraints

The terminal impact conditions are(2)
X
m
(
t
f
)
=
X
t
,
Y
m
(
t
f
)
=
Y
t
,
γ
m
=
γ
f
,
where γ
f and t
f are the specified impact angle and time, respectively. When multiple missiles are synchronized by their onboard clocks separately, a salvo attack can be achieved automatically without any communication. Therefore, there is less risk when encountering sophisticated electronic jamming and deception environment.


#### 2.2.2. Control Constraint

The acceleration is bounded as(3)
|
A
|
≤
A
max
⁡
,
A
max
⁡
>
0
,
which are determined by the airframe and the engine operating conditions of specific missiles.


#### 2.2.3. No-Fly Zone Constraints

The cruise missile often flies over a wide range of area. There exist several types of no-fly zones, such as dangerous terrains and adverse air defense areas, which should be evaded completely. Without loss of generality, these no-fly zones can be represented by disks as(4)
(
X
m
(
t
)
-
x
i
)
2
+
(
Y
m
(
t
)
-
y
i
)
2
≤
R
i
2
,
i
=
1,2
,
…
n
,
where (
x
i
,
y
i
) and R
i are the center and the radius, respectively, and n is the number of no-fly zones.


#### 2.2.4. Cost Function

The accumulative control energy during the guidance course(5)
J
(
t
0
)
=
1
2
min
⁡
∫
t
0
t
f
u
m
2
d
t
is often employed as the cost function to produce a more or less straight trajectory, where t
0 is the initial time instant.For most guidance laws, the maximal control energy often occurs at the terminal end, which could cause large miss distance. To solve this problem, a time-weight function is introduced in the cost objective as(6)
J
(
t
0
)
=
1
2
min
⁡
∫
t
0
t
f
(
f
p
(
t
)
·
u
m
2
)
d
t
.
As the impact time is fixed, f
p
(
t
) can be designed as a monotonically increasing function according to the terminal requirements.


## 2.1. Kinematics Equations

The missile-target relative kinematics is shown in Figure1. Without loss of generality, it is assumed that the velocity of the missile, V, is a constant. The earth coordinate system is chosen as X
O
Y. The current states of the missile are the position (
X
m
,
Y
m
) and the heading angle γ
m. The control variable for the missile is the normal acceleration, A, being perpendicular to V, which can change the velocity direction. The corresponding states of the target are the position (
X
t
,
Y
t
) and the heading angle γ
t.Figure 1
Missile-target relative kinematics.According to Figure1, the mass-point dynamics of the missile and the target are [9, 29]
(1)
d
X
m
d
t
=
V
cos
⁡
(
γ
m
)
,
d
Y
m
d
t
=
V
sin
(
γ
m
)
,
d
γ
m
d
t
=
u
m
=
A
V
,
d
X
t
d
t
=
V
t
cos
⁡
⁡
(
γ
t
)
,
d
Y
t
d
t
=
V
t
sin
(
γ
t
)
.


## 2.2. Optimal Trajectory Requirement

### 2.2.1. Terminal Constraints

The terminal impact conditions are(2)
X
m
(
t
f
)
=
X
t
,
Y
m
(
t
f
)
=
Y
t
,
γ
m
=
γ
f
,
where γ
f and t
f are the specified impact angle and time, respectively. When multiple missiles are synchronized by their onboard clocks separately, a salvo attack can be achieved automatically without any communication. Therefore, there is less risk when encountering sophisticated electronic jamming and deception environment.


### 2.2.2. Control Constraint

The acceleration is bounded as(3)
|
A
|
≤
A
max
⁡
,
A
max
⁡
>
0
,
which are determined by the airframe and the engine operating conditions of specific missiles.


### 2.2.3. No-Fly Zone Constraints

The cruise missile often flies over a wide range of area. There exist several types of no-fly zones, such as dangerous terrains and adverse air defense areas, which should be evaded completely. Without loss of generality, these no-fly zones can be represented by disks as(4)
(
X
m
(
t
)
-
x
i
)
2
+
(
Y
m
(
t
)
-
y
i
)
2
≤
R
i
2
,
i
=
1,2
,
…
n
,
where (
x
i
,
y
i
) and R
i are the center and the radius, respectively, and n is the number of no-fly zones.


### 2.2.4. Cost Function

The accumulative control energy during the guidance course(5)
J
(
t
0
)
=
1
2
min
⁡
∫
t
0
t
f
u
m
2
d
t
is often employed as the cost function to produce a more or less straight trajectory, where t
0 is the initial time instant.For most guidance laws, the maximal control energy often occurs at the terminal end, which could cause large miss distance. To solve this problem, a time-weight function is introduced in the cost objective as(6)
J
(
t
0
)
=
1
2
min
⁡
∫
t
0
t
f
(
f
p
(
t
)
·
u
m
2
)
d
t
.
As the impact time is fixed, f
p
(
t
) can be designed as a monotonically increasing function according to the terminal requirements.


## 2.2.1. Terminal Constraints

The terminal impact conditions are(2)
X
m
(
t
f
)
=
X
t
,
Y
m
(
t
f
)
=
Y
t
,
γ
m
=
γ
f
,
where γ
f and t
f are the specified impact angle and time, respectively. When multiple missiles are synchronized by their onboard clocks separately, a salvo attack can be achieved automatically without any communication. Therefore, there is less risk when encountering sophisticated electronic jamming and deception environment.


## 2.2.2. Control Constraint

The acceleration is bounded as(3)
|
A
|
≤
A
max
⁡
,
A
max
⁡
>
0
,
which are determined by the airframe and the engine operating conditions of specific missiles.


## 2.2.3. No-Fly Zone Constraints

The cruise missile often flies over a wide range of area. There exist several types of no-fly zones, such as dangerous terrains and adverse air defense areas, which should be evaded completely. Without loss of generality, these no-fly zones can be represented by disks as(4)
(
X
m
(
t
)
-
x
i
)
2
+
(
Y
m
(
t
)
-
y
i
)
2
≤
R
i
2
,
i
=
1,2
,
…
n
,
where (
x
i
,
y
i
) and R
i are the center and the radius, respectively, and n is the number of no-fly zones.


## 2.2.4. Cost Function

The accumulative control energy during the guidance course(5)
J
(
t
0
)
=
1
2
min
⁡
∫
t
0
t
f
u
m
2
d
t
is often employed as the cost function to produce a more or less straight trajectory, where t
0 is the initial time instant.For most guidance laws, the maximal control energy often occurs at the terminal end, which could cause large miss distance. To solve this problem, a time-weight function is introduced in the cost objective as(6)
J
(
t
0
)
=
1
2
min
⁡
∫
t
0
t
f
(
f
p
(
t
)
·
u
m
2
)
d
t
.
As the impact time is fixed, f
p
(
t
) can be designed as a monotonically increasing function according to the terminal requirements.


## 3. A Brief Introduction to GPM

Consider a dynamic equation(7)
x
˙
=
f
(
x
(
t
)
,
u
(
t
)
,
t
)
,
where x
(
t
)
∈
R
n and u
(
t
)
∈
R
m are state and control variables, respectively, and the function f: R
n
×
R
m
×
R
→
R
n. The traditional optimal control problem can be formulated in a unified form as
(8)
J
=
Φ
(
x
(
t
f
)
,
t
f
)
+
∫
t
0
t
f
g
(
x
(
t
)
,
u
(
t
)
,
t
)
d
t
,
subject to the boundary equality constraints
(9)
ϕ
(
x
(
t
0
)
,
t
0
,
x
(
t
f
)
,
t
f
)
=
0
,
and the process inequality constraints
(10)
C
(
x
(
t
)
,
u
(
t
)
,
t
0
,
t
f
)
≤
0
,
where ϕ
:
R
n
×
R
×
R
n
×
R
→
R
q and C
:
R
n
→
R
r.A brief review of pseudospectral method (PSM), especially Gauss PM (GPM), is presented here. The detailed description of GPM can be found in [24].Firstly, the original time interval[
t
0
,
t
f
] is transformed to the interval [
-
1,1
] as
(11)
τ
=
2
t
t
f
-
t
0
-
t
f
+
t
0
t
f
-
t
0
.
Using (11), (8) can be reformulated in terms of τ as
(12)
J
=
Φ
(
x
(
1
)
,
t
f
)
+
t
f
-
t
0
2
∫
-
1
1
g
(
x
(
τ
)
,
u
(
τ
)
,
τ
)
d
τ
.Similarly, the dynamic (7) and the boundary constraints (9) are replaced by
(13)
2
t
f
-
t
0
d
x
d
τ
=
f
(
x
(
τ
)
,
u
(
τ
)
,
τ
)
,
(14)
ϕ
(
x
(
-
1
)
,
t
0
,
x
(
1
)
,
t
f
)
=
0
,
respectively.In the GPM, the Legendre-Gauss (LG) points,τ
k (1
≤
k
≤
N
-
1), distributed on the interval [
-
1,1
], are defined as the roots of
(15)
P
N
(
τ
)
=
1
2
N
N
!
d
N
d
τ
N
[
(
τ
2
-
1
)
N
]
,
while τ
0
=
-
1, τ
N
=
1. The continuous state and control variables are approximated by the Nth polynomials as
(16)
x
(
τ
)
≈
X
(
τ
)
=
∑
i
=
0
N
L
i
(
τ
)
X
(
τ
i
)
,
(17)
u
(
τ
)
≈
U
(
τ
)
=
∑
i
=
1
N
L
i
(
τ
)
U
(
τ
i
)
,
where the Lagrange interpolation polynomial is
(18)
L
i
(
τ
)
=
∏
j
=
0
,
j
≠
i
N
τ
-
τ
i
τ
i
-
τ
j
.Secondly, (16) is differentiated at the node points as [28]
(19)
x
˙
(
τ
)
≈
X
˙
(
τ
)
=
∑
i
=
0
N
L
˙
i
(
τ
)
X
(
τ
i
)
,
where
(20)
D
k
i
=
L
˙
i
(
τ
k
)
=
{
(
1
+
τ
k
)
P
˙
K
(
τ
k
)
+
P
K
(
τ
k
)
(
τ
k
-
τ
i
)
[
(
1
+
τ
i
)
P
˙
K
(
τ
i
)
+
P
K
(
τ
i
)
]
,
i
≠
k
,
(
1
+
τ
i
)
P
¨
K
(
τ
i
)
+
2
P
˙
K
(
τ
i
)
2
[
(
1
+
τ
i
)
P
˙
K
(
τ
i
)
+
P
K
(
τ
i
)
]
,
i
=
k
.
The static form of (7) can be obtained by equating the right side of both (19) and (7) at the discrete node points.Thirdly, the continuous-time cost function (8) is approximated by using Gauss quadrature formula as
(21)
J
=
Φ
(
X
0
,
t
0
,
X
f
,
t
f
)
+
t
f
-
t
0
2
∑
k
=
1
N
w
k
g
(
X
k
,
U
k
,
τ
k
;
t
0
,
t
f
)
,
where w
k are the Gauss weights.Similarly, the terminal and process constraints can also be reformulated.According to the procedures described above, the original optimal control problem can be approximated as a static nonlinear programing (NLP) one. The essence of GPM is to replace the original infinite-dimension dynamic optimal control problem with a finite-dimension static NLP by eliminating differential and integral equations. There are many effective methods to solve NLP and among them the sequential quadratic programing (SQP) is a famous one which has been widely used because of its maturity. Nowadays, a kind of reliable software, SNOPT (Sparse Nonlinear OPTimizer), is available for solving SQP problems in a unified framework.


## 4. Receding Horizon GPM

The GPM provides us with an effective way to deal with nonlinear optimal control problem subject to various constraints directly, avoiding the shortcomings existing in the linearization approximation approach [7]. Moreover, a free GPOPS is available to seek GPM solution in real-time and reduce the huge computational complexity existing in the traditional trajectory planning [14, 19, 21]. However, the optimal solution obtained by using GPOPS is an open-loop one which is unable to cope with moving surface targets such as warships. Therefore, a closed-loop optimal solution should be generated instead. Enlightened by the mechanism of model predictive control [30], this obstacle can be eliminated by using receding horizon scheme which can be realized online.The receding horizon GPM can be implemented as follows.Step 1.
Input the number of LG points,N, the specified impact time t
f, and terminal heading angle γ
f. Set t
0 as the initial time instant and i
=
1.Step 2.
According to the guidance measurement information, plan an optimal trajectory from the current position to the target and obtain a control lawu
*
(
t
i
) using GPM. Assume that this computational period is Δ
T
i. During the time interval [
t
i
,
t
i
+
Δ
T
i
], or when the new optimization results are not available, the corresponding values of u
*
(
t
i
-
1
) on this interval are implemented. And then u
*
(
t
i
) is used after the time instant of t
i
+
Δ
T
i.Step 3.
Leti
=
i
+
1 and repeat Step 2 until t
=
t
f.This method can be illustrated in Figure2.Figure 2
Receding horizon optimization schematic diagram.


## 5. Numerical Simulations

In this section, the effectiveness of the proposed method is illustrated by several numerical simulations. The first simulation is to illustrate the effectiveness of the time-weight functionf
p
(
t
), wherein a stationary target is employed without using receding horizon optimization. The second and the third simulations combined are designed to test the capacity of the proposed method to deal with no-fly zone constraints. Here a moving target is used to evaluate the receding horizon optimization strategy. In all the simulations, the velocities of missiles keep a constant of 250 m/s and the initial numbers of LG points, N, are universally selected as 20.


### 5.1. Selection of Time-Weight Functionf
p
(
t
)

For most guidance laws, the control energy almost reaches its peak at the terminal end due to the divergent line-of-sight (LOS) rate, which would lead to a large miss distance. Therefore, it would be preferred to ensure the terminal trajectory to be straight for a period of time. In order to find a better time-weight function, five types of time-dependent functionsf
p
(
t
) are investigated and they are(1)
f
p
(
t
)
=
1;
(2)
f
p
(
t
)
=
t
2;
(3)
f
p
(
t
)
=
1
/
(
t
f
-
t
+
0.5
);
(4)
f
p
(
t
)
=
1
/
(
t
f
-
t
+
0.5
)
2;
(5)
f
p
(
t
)
=
1
/
(
t
f
-
t
+
0.5
)
3.Because of the terminal impact time specification, these functions can be easily realized. Assume a scenario that the missile impacts a stationary target subject to t
f
=
50 s and γ
f
=
0
°. The simulation results are shown in Figures 3 and 4 when using a one-phase GPM. According to Figure 3, it can be seen that all the terminal constraints are not violated for each function. From Figure 4, it is clear that when f
p
(
t
)
=
1, the terminal control amplitude is the maximal one among all; when f
p
(
t
)
=
t
2, the initial acceleration is large, and the control curve is similar to a bang-bang control. However, this function is free of t
f, which could be sensitive to the case configuration; when f
p
(
t
)
=
1
/
(
t
f
-
t
+
0.5
)
2, the terminal flight trajectory is straight. But with the increment of power in the denominator, both the initial acceleration and the total energy consumption are increased. Therefore, a comprehensive tradeoff is made and f
p
(
t
)
=
1
/
(
t
f
-
t
+
0.5
)
2 is chosen and fixed for the following simulations.Figure 3
Flight trajectories.Figure 4
Acceleration curves.


### 5.2. Without Process Inequality Constraint

In this scenario, two antiship missiles attack a moving warship without any no-fly zone constraint. The fixed sampling period isΔ
T
=
0.5 s. The velocity of the warship is 15 m/s with γ
t
=
-
13
5
°. The maximal feasible acceleration is A
max
⁡
=
45 m/s2. The case configurations are given in Table 1. The simulation results are shown in Figure 5 and Table 2. In Figure 5(a), the real lines are the actual trajectories and the point lines are the planned optimal trajectories at different sampling time instants. It can be found that both missiles fly along circular arcs during their early stages in order to consume the excessive time. According to Figures 5(b) and 5(c), the terminal trajectories are becoming straight gradually towards specified impact direction prior to impact instant; therefore tiny accelerations are needed in the final phase. In Table 2, (
Δ
X
,
Δ
Y
) and Δ
γ are the terminal errors of position and impact angle, respectively. It can be found that they are negligible from Table 2. This fact suggests that the proposed method can effectively impact a moving surface target subject to terminal constraints. It should also be noted that it is not necessary to pay a lot of attention to the selection of the initial values. In other words, this algorithm is insensitive to the quality of initial solution. The setup of the computer is i5 processor and 2G memory. The GPOPS is running on Matlab platform. When the algorithm is implemented in C code, the convergent time can be reduced by a factor of up to 100. For this case, each trajectory can be obtained within 10 s by using Matlab; therefore it meets real-time requirement when using C programming language onboard. In fact, we rewrote the corresponding GPOPS functions used in the abovementioned path planning by the C programming language on the basis of the commercial SNOPT software for practitioners. When implemented on the practical DOS platform, the running time can be kept below 0.3 s for almost all situations.Table 1
Case configurations without inequality constraints.
Missile
(
X
0
,
Y
0
) (km)
γ
0 (deg)
γ
f (deg)
t
f (s)
1
(−3.5, 8.5)
0
−105
50
2
(−9, −4.5)
45
0
50Table 2
Simulation results without inequality constraints.
Missile
Δ
X (m)
Δ
Y (m)
Δ
γ (deg)
1
8.8
e
-
08
4.1
e
-
8
0
2
3.7
e
-
08
-
8.5
e
-
8
0Salvo attack without inequality constraints.
(a)
Flight trajectories
(b)
Acceleration curves of the first missile
(c)
Heading angle curves of the second missile


### 5.3. With No-Fly Zone Constraints

In this case, we assume that there are two no-fly zones along the trajectories of the previous simulation and they are(1)
(
X
m
+
1300
)
2
+
(
Y
m
-
8000
)
2
≤
122500;
(2)
(
X
m
+
8000
)
2
+
(
Y
m
-
0
)
2
≤
122500.They are shown as two disks in Figure 6, which can be considered as the process constraints. It is obvious that the no-fly zones are just located along the current trajectories; hence the trajectories need to be replanned to evade the no-fly zones. The new trajectories are shown in Figure 7 when adding the no-fly zone constraints. According to Figure 7, the two missiles evade the no-fly zones successfully in the simulations. On the other hand, the running time is almost unchanged compared with the one in the above section; that is, the process inequality constraints have little effect on the optimization time for GPM.Figure 6
Location of the no-fly zones.Figure 7
Salvo attack with no-fly zones.


### 5.4. Comparison with the Analytical Guidance Law

In [7], an analytical guidance law with impact angle and time specifications was derived based on the linearization approximation using small heading angle assumption. In this section, the comparative simulation is performed. The case configurations are shown in Table 3. The first and second missiles are guided by the method in [7] while the trajectory of the third missile is obtained by using the proposed receding horizon GPM. The simulation results are shown in Figure 8 and Table 4. It is evident that poor miss distance is caused for the second missile due to the linearization approximation errors. However, the third missile can still guarantee impact accuracy.Table 3
Case configurations for comparison.
Missile
(
X
0
,
Y
0
) (km)
γ
0 (deg)
γ
f (deg)
t
f (s)
1
(−3, −3)
105
0
50
2
(−3, −3)
105
−25
50
3
(−3, −3)
105
−25
50Table 4
Comparative precisions.
Missile
Δ
X (m)
Δ
Y (m)
Δ
γ (deg)
1
3.7
e
-
08
-
8.5
e
-
8
−0.0002
2
−503.85
533.83
0.0699
3
1.2
e
-
09
-
2.1
e
-
9
0Figure 8
Comparative simulation results.


## 5.1. Selection of Time-Weight Functionf
p
(
t
)

For most guidance laws, the control energy almost reaches its peak at the terminal end due to the divergent line-of-sight (LOS) rate, which would lead to a large miss distance. Therefore, it would be preferred to ensure the terminal trajectory to be straight for a period of time. In order to find a better time-weight function, five types of time-dependent functionsf
p
(
t
) are investigated and they are(1)
f
p
(
t
)
=
1;
(2)
f
p
(
t
)
=
t
2;
(3)
f
p
(
t
)
=
1
/
(
t
f
-
t
+
0.5
);
(4)
f
p
(
t
)
=
1
/
(
t
f
-
t
+
0.5
)
2;
(5)
f
p
(
t
)
=
1
/
(
t
f
-
t
+
0.5
)
3.Because of the terminal impact time specification, these functions can be easily realized. Assume a scenario that the missile impacts a stationary target subject to t
f
=
50 s and γ
f
=
0
°. The simulation results are shown in Figures 3 and 4 when using a one-phase GPM. According to Figure 3, it can be seen that all the terminal constraints are not violated for each function. From Figure 4, it is clear that when f
p
(
t
)
=
1, the terminal control amplitude is the maximal one among all; when f
p
(
t
)
=
t
2, the initial acceleration is large, and the control curve is similar to a bang-bang control. However, this function is free of t
f, which could be sensitive to the case configuration; when f
p
(
t
)
=
1
/
(
t
f
-
t
+
0.5
)
2, the terminal flight trajectory is straight. But with the increment of power in the denominator, both the initial acceleration and the total energy consumption are increased. Therefore, a comprehensive tradeoff is made and f
p
(
t
)
=
1
/
(
t
f
-
t
+
0.5
)
2 is chosen and fixed for the following simulations.Figure 3
Flight trajectories.Figure 4
Acceleration curves.


## 5.2. Without Process Inequality Constraint

In this scenario, two antiship missiles attack a moving warship without any no-fly zone constraint. The fixed sampling period isΔ
T
=
0.5 s. The velocity of the warship is 15 m/s with γ
t
=
-
13
5
°. The maximal feasible acceleration is A
max
⁡
=
45 m/s2. The case configurations are given in Table 1. The simulation results are shown in Figure 5 and Table 2. In Figure 5(a), the real lines are the actual trajectories and the point lines are the planned optimal trajectories at different sampling time instants. It can be found that both missiles fly along circular arcs during their early stages in order to consume the excessive time. According to Figures 5(b) and 5(c), the terminal trajectories are becoming straight gradually towards specified impact direction prior to impact instant; therefore tiny accelerations are needed in the final phase. In Table 2, (
Δ
X
,
Δ
Y
) and Δ
γ are the terminal errors of position and impact angle, respectively. It can be found that they are negligible from Table 2. This fact suggests that the proposed method can effectively impact a moving surface target subject to terminal constraints. It should also be noted that it is not necessary to pay a lot of attention to the selection of the initial values. In other words, this algorithm is insensitive to the quality of initial solution. The setup of the computer is i5 processor and 2G memory. The GPOPS is running on Matlab platform. When the algorithm is implemented in C code, the convergent time can be reduced by a factor of up to 100. For this case, each trajectory can be obtained within 10 s by using Matlab; therefore it meets real-time requirement when using C programming language onboard. In fact, we rewrote the corresponding GPOPS functions used in the abovementioned path planning by the C programming language on the basis of the commercial SNOPT software for practitioners. When implemented on the practical DOS platform, the running time can be kept below 0.3 s for almost all situations.Table 1
Case configurations without inequality constraints.
Missile
(
X
0
,
Y
0
) (km)
γ
0 (deg)
γ
f (deg)
t
f (s)
1
(−3.5, 8.5)
0
−105
50
2
(−9, −4.5)
45
0
50Table 2
Simulation results without inequality constraints.
Missile
Δ
X (m)
Δ
Y (m)
Δ
γ (deg)
1
8.8
e
-
08
4.1
e
-
8
0
2
3.7
e
-
08
-
8.5
e
-
8
0Salvo attack without inequality constraints.
(a)
Flight trajectories
(b)
Acceleration curves of the first missile
(c)
Heading angle curves of the second missile


## 5.3. With No-Fly Zone Constraints

In this case, we assume that there are two no-fly zones along the trajectories of the previous simulation and they are(1)
(
X
m
+
1300
)
2
+
(
Y
m
-
8000
)
2
≤
122500;
(2)
(
X
m
+
8000
)
2
+
(
Y
m
-
0
)
2
≤
122500.They are shown as two disks in Figure 6, which can be considered as the process constraints. It is obvious that the no-fly zones are just located along the current trajectories; hence the trajectories need to be replanned to evade the no-fly zones. The new trajectories are shown in Figure 7 when adding the no-fly zone constraints. According to Figure 7, the two missiles evade the no-fly zones successfully in the simulations. On the other hand, the running time is almost unchanged compared with the one in the above section; that is, the process inequality constraints have little effect on the optimization time for GPM.Figure 6
Location of the no-fly zones.Figure 7
Salvo attack with no-fly zones.


## 5.4. Comparison with the Analytical Guidance Law

In [7], an analytical guidance law with impact angle and time specifications was derived based on the linearization approximation using small heading angle assumption. In this section, the comparative simulation is performed. The case configurations are shown in Table 3. The first and second missiles are guided by the method in [7] while the trajectory of the third missile is obtained by using the proposed receding horizon GPM. The simulation results are shown in Figure 8 and Table 4. It is evident that poor miss distance is caused for the second missile due to the linearization approximation errors. However, the third missile can still guarantee impact accuracy.Table 3
Case configurations for comparison.
Missile
(
X
0
,
Y
0
) (km)
γ
0 (deg)
γ
f (deg)
t
f (s)
1
(−3, −3)
105
0
50
2
(−3, −3)
105
−25
50
3
(−3, −3)
105
−25
50Table 4
Comparative precisions.
Missile
Δ
X (m)
Δ
Y (m)
Δ
γ (deg)
1
3.7
e
-
08
-
8.5
e
-
8
−0.0002
2
−503.85
533.83
0.0699
3
1.2
e
-
09
-
2.1
e
-
9
0Figure 8
Comparative simulation results.


## 6. Conclusions

This paper presented a trajectory optimization approach based on the receding horizon GPM to impact a moving surface target subject to the terminal impact angle and the time specifications. A time-dependent weighting function was introduced to reduce the control energy at the terminal end. The receding horizon strategy could be applied to the target in motion. When using GPOPS, this comprehensive scheme can handle the process inequality constraints without causing the increment of computational complexity and no linearization needs to be done. The optimization algorithm could be implemented real-time with an onboard computer. The effectiveness of the proposed method was illustrated by several numerical examples.



---

*Source: 604705-2014-05-14.xml*