Abstract
The predicted impact point (PIP) of hypersonic interception changes continually; therefore the midcourse guidance law must have the ability of online trajectory optimization. In this paper, an online trajectory generation algorithm is designed based on neighboring optimal control (NOC) theory and improved indirect Radau pseudospectral method (IRPM). A trajectory optimization model is designed according to the features of operations in near space. Two-point boundary value problems (TPBVPs) are obtained based on NOC theory. The second-order linear form of transversality conditions is deduced backward to express the modifications of terminal states, costates, and flight time in terms of current state errors and terminal constraints modifications. By treating the current states and the optimal costates modifications as initial constraints and perturbations, the feedback control variables are obtained based on improved IRPM and nominal trajectory information. The simulation results show that when the changes of terminal constraints are not relatively large, this method can generate a modified trajectory effectively with high precision of terminal modifications. The design concept can provide a reference for the design of the online trajectory generation system of hypersonic vehicles.
1. Introduction
After decades of development, hypersonic technology has made a lot of progress [1–3]. Hypersonic vehicles will be used in military as precision strike weapons or platforms in the near future. The demand on the research of advanced intercept and defense technology is urgent.
In allusion to the high speed and large maneuverability of hypersonic weapons, the interceptor should adopt the compound guidance strategy to improve the success rate of interception. From the end of program control flight to the target acquisition of terminal guidance, the interceptor spends most of the time flying in the midcourse guidance. So the flight performance of midcourse guidance determines the overall performance.
If we adopt the traditional guidance method to intercept the high dynamic target directly, it will cause the frequent changes of the trajectory, which is not conducive to the implementation of the interception and leads to unnecessary loss of energy. So it is a reasonable strategy to design the trajectory of midcourse guidance. The main objective of the midcourse guidance is to minimize the energy consumption in the flight process of interceptor and enter the terminal guidance with the best relative geometric relationship [4]. At present, the challenges which the midcourse guidance law design of near space interceptors faces are mainly in the following two aspects: on the one hand, during the flight process of interceptors, problems such as structural strength, thermal protection, normal work of the engine, and control stability give rise to strict requirements on the heat flux, dynamic pressure, overload, control, and other characteristics of the process. Additionally, the predicted impact point (PIP) provides terminal position constraints, and the terminal guidance acquisition state-space becomes the strong constraints conditions of interceptors [5]. On the other hand, in actual flight process, the variation ranges of altitude and velocity are relatively large. Due to the initial conditions errors, atmospheric environment changes, and great uncertainty of aerodynamic model and navigation equipment, it will cause trajectory tracking errors. Furthermore, the nominal trajectory generated before launch is based on remote detection results. In this condition, tracking and prediction of the target have great error. With the approaching of the interceptor and the target, tracking and prediction accuracy will be improved by means of approaching detection with on-board equipment. And the target is also likely to maneuver actively. All those will lead to the change of PIP, and terminal constraints also need updates accordingly. Therefore, the designed midcourse guidance law must have the ability of online trajectory optimization.
Taking all kinds of terminal constraints and a series of process characteristics into account, the online optimization problem is actually a complex nonlinear optimization problem [6], and many scholars have studied it. Lin and Tsai [7, 8] obtained a trajectory shaping guidance law by simplifying the model and realized online trajectory generation. But the trajectory was not true optimal due to too much approximation. Indig et al. [6, 9, 10] first linearized the dynamic model of interceptors. According to the location of PIP, the optimal guidance law was designed based on the trajectory shaping guidance. Then, in allusion to the actual nonlinear dynamic model, the terminal angle constraints of interceptor were added according to the midcourse and terminal guidance handover. Finally, the optimal trajectory model was deduced based on Pontryagin minimum principle and solved by Gauss pseudospectral method. Dwivedi et al. [11, 12] studied the problem of midcourse guidance trajectory planning using model predictive static planning method. Based on the state equations of discretized system, the control variables were updated by solving the costate vectors of the whole time period, whose coefficient matrix could be obtained by recursive method. The simulation results show that this method has high accuracy and efficiency. Yakimenko et al. [13–15] proposed a trajectory shaping guidance law based on the inverse dynamic of virtual domain. In this method, the coordinates of aircraft were expressed by higher order polynomials with virtual arc, and the dynamic equations were transformed from the time domain to the virtual domain. Introducing the virtual arc as an argument made it possible to optimize the speed history along the trajectory independently. Additionally, the number of the optimized variables was reduced and the integration process was avoided. By this method, the trajectory can be generated with good convergence robustness. In addition, some scholars have also studied this problem by pseudospectral method [16, 17], intelligent algorithm [18, 19], and so on.
The aforementioned scholars have made many accomplishments in the domain of trajectory optimization; however, they mostly focused on the trajectory regeneration in regard to the research of online trajectory optimization problems. They completely abandoned the original trajectory and reoptimized again, which caused large amount of calculation and more time consumption, and put forward higher requirements on the storage capacity and computing ability of on-board computer. So it is not easy to realize in engineering. In allusion to the characteristics of near space targets, such as small overload, long flight time, and large flight range, prediction of impact point will be continually updated and the terminal constraints change correspondingly. When the perturbations of current states and changes of terminal constraints are not very large, according to the neighboring optimal control theory, the modified trajectory can be generated quickly within the neighborhood of the nominal trajectory using nominal trajectory data to realize online trajectory generation.
We still need to solve the TPBVPs [20] when adopting above modification algorithm. The time-intensive backward integration has brought difficulties to the online implementation. Recently, Fahroo and Ross [21] and Yan et al. [22] proposed the indirect Legendre pseudospectral method (ILPM) for solving the TPBVPs by matrix inversion which improved solving efficiency a lot. In this method, the modifications of states and costates were solved by initial perturbations of states and the nominal trajectory information, and then the feedback control modifications were obtained to form the closed loop control. In this paper, the idea of ILPM is used for reference and the IRPM is proposed which is improved in two aspects. On the one hand, Fahroo and Hui only solved the trajectory tracking problems with initial perturbations, while the improved IRPM considers not only trajectory tracking errors but also modifications of terminal constraints. If computation speed is fast enough, the trajectory can be generated in real time to realize closed loop guidance without designing the feedback control law. On the other hand, Radau pseudospectral method is superior to Legendre pseudospectral method in approximating precision of control, state, and costate variables. In terms of computational efficiency, the two methods have little difference when solving the same scale problem.
The main purpose of this paper is to design an online optimal modification algorithm and realize online trajectory optimization under the condition of changing terminal constraints. The rest of paper is organized as follows. In Section 2, the nominal trajectory satisfying terminal and path constraints is generated offline. Next, the trajectory modified model with unspecified terminal time is deduced by NOC in Section 3 and is solved by improved IRPM in Section 4. In Section 5, numerical simulations are presented. Finally, some conclusions are drawn in Section 6.
2. Nominal Optimal Trajectory
2.1. Dynamic Model of Interceptor
When the nominal optimal trajectory is generated offline, the interceptor should be directly guided towards the PIP, so the dynamic model is considered in the longitudinal plane. The states of interceptor are represented by four-dimensional vector consisting of speed, flight path angle, and location in the inertial coordinate and given as follows [23]:where is interceptor mass; is the dynamic pressure; is the reference area; is gravitational acceleration; is the engine thrust. In order to achieve a higher flight speed, the two-stage booster program is designed. The thrust calculation formulas are shown in where is the fuel gas velocity and is the specific impulse, whose value is generally 200~300 s for the solid rocket engine [24]. This paper designs a two-stage solid rocket booster. The related parameters of the engine are shown in Table 1.
denote the drag and lift coefficients, respectively, which are expressed as a function of , the Mach number, and , the angle of attack:where denote zero-lift drag coefficient and induced drag coefficient, respectively; denotes the partial derivative of with respect to . denotes the dynamic pressure:
And is the atmosphere density expressed aswhere .
2.2. Trajectory Optimization Problem
In order to ensure kinematic kill effect of interception, terminal speed is usually selected as cost function:
In this paper, subscript denotes initial conditions and subscript denotes terminal conditions. As is discussed above, the main objective of the midcourse guidance is to deliver the interceptor to a certain position (usually the PIP) with specified conditions to ensure a successful handover between midcourse and terminal guidance. So the position constraints are given by PIP and the flight path angle constraint at terminal time is given by the terminal guidance acquisition conditionswhere denotes zero matrix with the corresponding dimension.
The above trajectory optimization problem is solved with the help of the software package GPOPS [25]. The whole trajectory was composed of three phases: boost stages 1 and 2 and nonpropulsive phase. Compared with other direct optimization methods, pseudospectral method can obtain higher accuracy with fewer nodes, which improves efficiency of the algorithm. The nominal optimal trajectory is shown in Figure 1.
3. Trajectory Modification Based on Neighborhood Optimization
3.1. Two-Point Boundary Value Problem
In the process of actual flight, when tracking errors of current states and changes of terminal constraints are not very large, according to the neighboring optimal control theory, we can further differentiate the first-order necessary conditions of the optimal trajectory to the second order to obtain perturbation equations, and the modified trajectory satisfying the changed terminal conditions can be obtained using the nominal trajectory data. At the same time, the cost function can still maintain certain optimality.
First of all, the first-order optimal necessary conditions are derived based on optimal control. By introducing the costate variable with the same dimension as the state variable , Hamilton equation can be expressed aswhere denotes system equations.
Then the augmented performance index with terminal constraints can be expressed aswhere . Canonical equation and coupled equation can be expressed as
Considering the change of terminal conditions, the terminal time of modified trajectory will change accordingly. We can derive the equations of boundary conditions specified at an unspecified terminal time and the terminal constraints of Hamiltonian function along the optimal trajectory based on NOC [26, 27]
Equations (10)~(12) can be further differentiated to second order as follows:
If is nonsingular in the whole flight process, we can get the expression of control modifications by (18):
The dynamic equations of state modifications and costate modifications can be obtained by substituting (19) into (16)~(17):where
Equations (20)~(21) satisfy initial conditions and terminal conditions (13)~(15) form the TPBVPs of .
3.2. Terminal Linearization and Backward Recursion
Equations (6), (15), and (16) were differentiated to second order with unspecified terminal time. Since , we can get
Equations (23)~(26) consist of linear equations with unknown variables. They are state modifications of dimension , costate modifications of dimension , variables of dimension , and variable of dimension one. Variables and state modifications of dimension unspecified at terminal time are defined as free variables
Since the changes of terminal constraints are not very large, (23)~(25) can be linearized based on small perturbation hypothesis, so that we getwhere
Equations (20)~(21) can be expressed as the following form:where and the variables of above differential equation are discretized as points by . We can get by backward recursion
Passing to the first point in turn we can getwhere , and according to (27) we haveSubstituting (32) into (31)whereand if is reversible, can be expressed asSubstituting (35) into (27) we can get
3.3. The Optimality Theory
Theorem 1. The modified trajectories obtained in the nominal trajectory neighborhood based on the neighborhood optimal control theory meet the changed terminal constraints, and the performance index has second-order optimality.
Proof. According to the division integral theorem, the performance index can be expressed asConsidering the first-order variation of the performance index and taking canonical equations (10)-(11) and coupled equation (12) into it, we can get by ignoring the high order terms:Multiplying (16) and (23) by the and , respectively, we can get by taking the resulting equations into (38):Considering the second-order variation of (39), we can get according to the division integral theorem:Taking (17)-(18) and (24) into (40), we can get for any and . Therefore, the second-order optimality of the performance index can be guaranteed with satisfying the changed terminal constraints.
4. Solving Algorithm Based on Improved IRPM
In the last section, considering the changes of terminal constraints of some state variables, the optimal terminal modification of costates variables and optimal modification of terminal time are derived at an unspecified terminal time based on the NOC. Improved IRPM will be used to solve the modified optimal trajectory according to and in this section.
Since has been calculated in the last section, the nominal trajectory information which will be used to solve the modified trajectory changes correspondingly. The new nominal trajectory information can be obtained by interpolating the original nominal trajectory. Because the trajectory tracking errors have been taken into account when we solved in the last section, we can regard current states of the interceptor as initial constraints and as the perturbations when solving the modified trajectory. Then we can use the new nominal trajectory information to design the solving algorithm to generate a trajectory from the current states to the modified terminal states, whose performance index can still maintain certain optimality.
The boundary conditions in last section change to the following form:
Equations (41)~(42) can be expressed as follows:
Variation of (43) is solved as
Now we use improved IRPM to solve the optimal trajectory modified model. Firstly, the variables are discretized at Legendre-Gauss-Radau (LGR) points. Since the time domain for the above problem is , while Radau pseudospectral method is used in the time domain , it is necessary to introduce a new time variable :
States and costates variables are approximated as follows:where and are LGR points distributed at the interval , which are defined as zero points of . is Legendre polynomial of order . And is Lagrange interpolation polynomial basis function:
The derivatives of and at the point of are obtained by differentiating (46) as follows:where is the difference matrix of dimension :
The TPBVPs of are transformed into the following algebraic equations:
represents the zero matrix with the corresponding dimension. For simplicity, substituting into (50)~(51) we can getwhere are matrices of dimension :where and are zero and identity matrices of dimension . The purpose of this algorithm is to solve (53) subjected to the boundary conditions (52). These equations can be summarized in matrix form aswhere , and are matrices of dimension . . Dividing and as and we getwhere of dimension and of dimension are block matrices of . is a vector of dimension , which can be obtained by solving
Since , we can getwhere and of dimension are the block matrices of . Then we havewhere of dimension is block matrix of and of dimension is block matrix of . Subscript represents the th LGR point.
Substituting (59) into (19), we can get the optimal control modifications:
Since has been solved in the last section, the modifications of states, costates, and controls at the LGR points can be obtained by (59)-(60) and the values at instants of time between the nodes can be obtained by interpolation. The above solutions are obtained without a wide range of iterative optimization process, so it can ensure the high computational efficiency.
The modified algorithm can be concluded as Figure 2 and specific procedures are as follows:(1)Considering various process constraints and terminal constraints and selecting the maximum velocity as the performance index, we can establish the optimal model according to the dynamic equations of interceptor.(2)After further differentiating the canonical equations and coupled equations to second order, we can get the expression of feedback control and the TPBVP equations of .(3)After further differentiating the boundary conditions at unspecified terminal time to second order, we can get the expressions of in terms of free variables and terminal constraints ; the relationship between and is established by backward recursion, which can be used to deduce the expressions of in terms of , and then we get the expressions of in terms of .(4)The flight time of original nominal trajectory plus can be used as new flight time of the new nominal trajectory which can be obtained by interpolating the original nominal trajectory.(5)Regarding current states of the interceptor as initial constraint and as perturbation, we solve the feedback control using the improved IRPM and finally realize online trajectory optimization of midcourse guidance.
5. Simulations and Discussion
In order to verify the validity of the proposed algorithm, the following assumption of intercept operations is designed. Because the nonpropulsive phase is more difficult for trajectory adjustment, the terminal PIP is modified at the initial time of the nonpropulsive phase with altitude increasing 2 km and flight path angle increasing 5 deg. And the current tracking errors of flight path angle and altitude are 0.1 deg and 10 m.
The GPM has high precision in solving optimization problems. So when the tracking error generates or the terminal constraints change, the reoptimized trajectory by GPM, which is solved by MATLAB programs GPOPS [25], will be used as a standard optimal modified trajectory for comparison. The comparison of modified trajectory is shown in Figure 3. NOC represents the modified trajectory generated based on NOC theory and the improved IRPM solving algorithm designed in this paper. Nominal represents the nominal trajectory with the unchanged terminal constraints; GPOPS represents the modified trajectory reoptimized by GPM.
All simulation results of this paper are obtained on a Lenovo laptop with an Intel 2.50 GHz quad-core processor, using windows 7 operating system. All the codes are run under the MATLAB® R2014a environment.
From Figure 3 we know that the modified trajectory (NOC), which is substantially overlapping the reoptimized trajectory (GPOPS), can satisfy the changed terminal constraints caused by the changed PIP. And the flight time of modified trajectory is approximatively equal to that of the reoptimized trajectory by GPM. These results prove the high accuracy of the algorithm designed in this paper. In the initial time after modification, the error of angle of attack between the modified trajectory and the nominal trajectory is not large, no more than 0.5 deg, which represents that the angle of attack does not have a large jump and therefore ensures the stability of the control.
In order to verify the higher precision and efficiency of improved IRPM, the contrast simulation is designed as follows: in the fourth section, according to the tracking error and terminal constraints modifications, and are obtained by backward recursive algorithm. And then we solve feedback control variables by inversely integrating differential equations (19)–(21) of trajectory modification model rather than improved IRPM. In the integration process, 50, 500, and 1000 nodes are, respectively, selected at the same time interval. The curves of flight path angle modification and altitude modification are shown in Figures 4–6; the comparison of time consumption and terminal constraints modifications is shown in Table 2.
From Figures 4 and 5 and Table 2 we know that compared with the modified trajectory solved by improved IGPM with that reoptimized by GPM, the most part of curves is close except for initial several seconds after modification with flight path angle error of no more than 0.15 degrees and altitude error of no more than 120 m. Furthermore, the two methods have similar precision in terminal modifications, but the time consumption of the former is much shorter than the latter. The method of solving differential equations directly has relatively large errors compared to reoptimization by GPM. With the increase of the number of nodes, the errors decrease gradually and the precision of terminal modifications accordingly increases. But the time consumption increases as well. And even if 1000 nodes are selected, the precision of terminal modifications still has a certain gap compared with improved IGPM while time consumption increases a lot.
In summary, the concept of modified trajectory design using nominal trajectory information can guarantee the modified precision with improving efficiency a lot; and the solving method based on improved IGPM method is superior to the method by solving differential equations directly in accuracy and efficiency.
In order to verify the adaptivity of modified algorithm, the simulation is designed as follows: the value of is set randomly within , and the value of is set from the −5 deg to 5 deg at the interval of 1 deg, which adds up to 11 sets of different changed terminal constraints. The results are shown in Figure 7. Similarly, the value of is set randomly within , the value of is set from the −3000 m to 3000 at the interval of 600 m which adds up to 11 sets of different changed terminal constraints. The results are shown in Figure 7. We can know from Figures 6 and 7 that the algorithm has the modified ability for different changes of terminal constraints. All initial errors of altitude are close to 0 and the maximal initial error of flight path angle is no more than 0.1 deg, which can be adjusted by the trajectory tracking control process.
is set as 10 deg and is set to 0; the curve of modifications is shown in Figure 8. Similarly, is set to 8000 m and is set to 0; the curve of modifications is shown in Figure 9. According to Figures 8 and 9 we can know that when terminal constraints change largely, the error between the modified trajectory generated by the algorithm of this paper and the reoptimized trajectory by GPM is relatively large, which means the optimality of modified trajectory is difficult to guarantee. Because the modified algorithm of this paper is based on the NOC theory, when the changes of terminal constraints are beyond the neighborhood of nominal trajectory, the algorithm is no longer applicable.
6. Conclusion
An online trajectory optimization algorithm of midcourse guidance for hypersonic interception is designed based on the NOC theory and the improved IGPM. The main conclusions are as follows:(1)The algorithm considers not only the tracking error but also the change of terminal constraints. Therefore, it has the disturbance rejection ability and the ability to adjust terminal constraints online.(2)By using the original nominal trajectory information, a modified trajectory with higher accuracy can be rapidly generated based on neighboring optimal control theory. The iterative optimization process is avoided, and the trajectory generation is realized online.(3)When the changes of terminal constraints are relatively large, beyond the neighborhood of the original nominal trajectory, the online modified algorithm is no longer applicable. It needs to reoptimize online to get the modified trajectory.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This study was cosupported by the National Natural Science Foundation of China (Grant no. 61573374 and no. 61503408) and Aeronautical Science Foundation of China (Grant no. 20140196004 and no. 20150196006).