Receding Horizon Trajectory Optimization with Terminal Impact Specifications

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.


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][2][3][4][5][6][7][8][9][10][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][13][14][15] and intelligent optimization algorithms, such as genetic algorithm (GA) [16][17][18], ant colony optimization (ACO) [19,20], and particle swarm optimization (PSO) [21][22][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][25][26][27][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 Section 2. 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.

Problem Formulation
2.1.Kinematics Equations.The missile-target relative kinematics is shown in Figure 1.Without loss of generality, it is assumed that the velocity of the missile, , is a constant.The earth coordinate system is chosen as .The current states of the missile are the position (  ,   ) and the heading angle   .The control variable for the missile is the normal acceleration, , being perpendicular to , which can change the velocity direction.The corresponding states of the target are the position (  ,   ) and the heading angle   .
According to Figure 1, the mass-point dynamics of the missile and the target are [9,29] where   and   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.

Control Constraint.
The acceleration is bounded as which are determined by the airframe and the engine operating conditions of specific missiles.

No-Fly Zone Constraints.
The cruise missile often flies over a wide range of area.There exist several types of nofly 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 where (  ,   ) and   are the center and the radius, respectively, and  is the number of no-fly zones.

Cost Function. The accumulative control energy during the guidance course
is often employed as the cost function to produce a more or less straight trajectory, where  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 As the impact time is fixed,   () can be designed as a monotonically increasing function according to the terminal requirements.

A Brief Introduction to GPM
Consider a dynamic equation where () ∈   and () ∈   are state and control variables, respectively, and the function :   ×   ×  →   .
The traditional optimal control problem can be formulated in a unified form as subject to the boundary equality constraints and the process inequality constraints where  :   ×  ×   ×  →   and  :   →   .A brief review of pseudospectral method (PSM), especially Gauss PM (GPM), is presented here.The detailed description of GPM can be found in [24].
In the GPM, the Legendre-Gauss (LG) points,   (1 ≤  ≤  − 1), distributed on the interval [−1, 1], are defined as the roots of while  0 = −1,   = 1.The continuous state and control variables are approximated by the th polynomials as where the Lagrange interpolation polynomial is Secondly, ( 16) is differentiated at the node points as [28] ẋ where 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 where   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.

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, , the specified impact time   , and terminal heading angle   .Set  0 as the initial time instant and  = 1.
Step 2. According to the guidance measurement information, plan an optimal trajectory from the current position to the target and obtain a control law  * (  ) using GPM.Assume that this computational period is Δ  .During the time interval [  ,   + Δ  ], or when the new optimization results are not available, the corresponding values of  * ( −1 ) on this interval are implemented.And then  * (  ) is used after the time instant of   + Δ  .
Step 3. Let  =  + 1 and repeat Step 2 until  =   .This method can be illustrated in Figure 2.

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 timeweight function   (), 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, , are universally selected as 20.

Selection of Time-Weight Function 𝑓 𝑝 (𝑡).
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 functions   () are investigated and they are (1)   () = 1; (2)   () =  (3)   () = 1/(  −  + 0.5); (4)   () = 1/(  −  + 0.5) 2 ; (5)   () = 1/(  −  + 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   = 50 s and   = 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   () = 1, the terminal control amplitude is the maximal one among all; when   () =  2 , the initial acceleration is large, and the control curve is similar to a bang-bang control.However, this function is free of   , which could be sensitive to the case configuration; when   () = 1/(  −  + 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   () = 1/(  −  + 0.5) 2 is chosen and fixed for the following simulations.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, (Δ, Δ) 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.
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 nofly 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.

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.

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.

Figure 6 :
Figure 6: Location of the no-fly zones.

Table 1 :
Case configurations without inequality constraints.

Table 2 :
Simulation results without inequality constraints.In this scenario, two antiship missiles attack a moving warship without any no-fly zone constraint.The fixed sampling period is Δ = 0.5 s.The velocity of the warship is 15 m/s with   = −135 ∘ .The maximal feasible acceleration is  max = 45 m/s 2 .The case configurations are given in Table

Table 3 :
Case configurations for comparison.