Abstract
To generate the midcourse guidance trajectory for intercepting the high-speed and high maneuvering target, which is a strongly nonlinear and strongly constrained problem, a two-stage convex optimization method is proposed to solve the optimal trajectory quickly. In the first stage, an initial trajectory generation method is proposed, by which the trajectory’s terminal state is close to the terminal position. And the generated trajectory is used as the initial solution of the convex optimization method. In the second stage, the original nonconvex optimization problem is transformed into a convex optimization problem by linearization and relaxation methods and then solved discretely. In the numerical simulation, the effectiveness of the proposed method is verified, and the robustness of the method is verified in different initial and terminal states. Then, several ablation experiments are operated to verify the advantage of the rapid initial trajectory generation method in the first stage. Finally, compared with the gauss pseudospectral method (GPM), the proposed method is proved to be efficient and has the potential of online trajectory generation.
1. Introduction
More and more countries are developing delivery vehicles and weapons to dominate aerospace [1]. The midcourse guidance of interceptor plays an important role in intercepting such targets and can provide a good intercepting situation for terminal guidance with path constraints [2].
There are two types of design methods for midcourse guidance. The first is to determine that the guidance law is designed to control the line of sight, when the relative dynamic equation of the interceptor and the target, by which the adverse influence of the target maneuver could be overcome [3]. However, in the midcourse guidance, the accurate information of the high-speed target is difficult, which will lead to the larger calculation error of the sight line [4]. To solve this problem, Wan et al. [5] used the model predictive static programming method to give a suboptimal guidance law with terminal angle constraints. Ebrahimi et al. [6] proposed a new generalized model predictive extended control technology, which smoothed the control signal. But during the intercept, there may be no guarantee that the interceptor meets some path restrictions. Ann et al. [7] combined zero effort trajectory with minimum time trajectory to predict intercept points. Liu et al. [8] proposed a midcourse guidance law considering the influence of random interference, detection range constraint, and target capture probability. However, the guidance accuracy of the first method is easily affected by target maneuver. The ability to intercept high speed and high maneuvering target at high altitude is limited [9].
The second is generating the guidance instructions based on the terminal information prediction with several path constrains and terminal constraints [10]. The pseudospectral method has high accuracy in solving the strongly nonlinear and strongly constrained optimal control problems and has been rapidly developed and applied in the aircraft trajectory optimization [11]. Li et al. [12] proposed a multi-interval mesh refinement Radau pseudospectral method, which improved the optimization accuracy and convergence speed of the algorithm. Zhou et al. [13] improved the particle swarm optimization algorithm to solve the trajectory planning problem in the glide phase of hypersonic aircraft. Aiming at target maneuver and interceptor terminal constraints correction, Zhou et al. [14] and Li et al. [15] established a trajectory optimization and correction problem of midcourse guidance with path constraints. Du et al. [16] presented a semianalytical method for solving the exoatmospheric midcourse guidance problem with minimum velocity increment. Li et al. [17] designed an angle-constrained midcourse guidance trajectory according to the offline optimized trajectory information. But the trajectory generation efficiency of the above method is not ideal for intercepting 3d maneuverable high-speed target.
In recent years, convex optimization methods have been widely used in aerospace aircraft due to the existence of solutions [18] and high computational efficiency of solving complex polynomials [19]. Liu et al. [20] and Wang and Grant [21] used linearization and relaxation methods to solve the entry trajectory optimization problem with multiple constraints based on second-order cone programming (SOCP). Yan et al. [22] presented the entry problem with multiple constraints as an easily solvable SOCP sequence. Wang and Grant [23] presented an autonomous entry aircraft guidance algorithm based on convex optimization and continuous closed loop method. Wang et al. [24] used convex optimization algorithms to provide high-quality initial guesses for pseudospectral methods. Hong et al. [25] presented an autonomous entry aircraft guidance algorithm based on convex optimization method and continuous closed loop. Liu et al. [26, 27] introduced a nonconvex constraint relaxation technique, which uses regularization technique to constrain the degree of relaxation and proves the effectiveness of the method. In order to improve the convergence rate of convex optimization, Wang and Lu [28] and Zhou et al. [29] proposed an adaptive method to adjust the shrinkage coefficient of trust region. Combining pseudospectral method and convex optimization method, Sagliano and Mooij [30] proposed a new drag-energy scheme for atmospheric entry guidance, which can los-less convexify the formulation. Zhou et al. [31] used an adaptive mesh refinement method to improve the computational efficiency of convex optimization method. The above methods show that convex optimization can solve complex polynomial constraint problems effectively.
In midcourse guidance, Yang et al. [32] proposed a midcourse guidance method for boost interception based on an adaptive update strategy of trust region in convex optimization. In multi-interceptor guidance coordination, Jiang et al. [33] used convex optimization technology to solve the multiconstraint optimal proportional guidance problem for each interceptor online to achieve time consistency between interceptors. In this paper, the convex optimization method is applied to the generation of midcourse guidance trajectory for intercepting high-speed and high maneuvering target.
In this paper, the trajectory optimization problem with strong nonlinearity and strong constraints is described, and the dynamics model of the interceptor is transformed by constructing affine variables. Then, the free terminal time is affined to the fixed time domain by the time expansion factor. The dynamic and path constraints are linearized when the control variables are relaxed, through which the problem are transformed into a convex optimization problem. Finally, the problem is discretized and the fast trajectory generation method is proposed, and the process of solving the problem is given.
The innovations of this paper are as follows: (1)The midcourse guidance problem with strong nonlinearity and strong constraints is transformed into a SOCP problem that can be solved iteratively by constructing an affine system and using relaxation and linearization methods. The effectiveness of the method is verified by simulation and the robustness of the method is verified under the condition of initial state perturbation and terminal position change(2)A fast trajectory generation method is proposed as the initial solution of the convex optimization method. The advantages of the fast trajectory generation method are proved by ablation experiments. Compared with the GPM, the proposed method shows its rapidity
The remainder of this paper are introduced as follows. In Section 2, the trajectory optimization model of midcourse guidance is constructed. In Section 3, the original problem is transformed into a convex optimization problem by using linearization and relaxation methods. In Section 4, the convex optimization problem is discretized, and the fast trajectory generation method and the process of solving the problem are given. In Section 5, the proposed two-stage convex optimization method is verified.
2. Midcourse Guidance Problem
2.1. Dynamics Model
For intercepting the high-speed and high maneuvering target, trajectory optimization of midcourse guidance is a trajectory optimization problem with strong nonlinearity and strong constraints.
The model assumes the following (1)As the earth’s rotation has little influence on the mid-course guidance flight of interceptor, the earth’s rotation is ignored(2)The interceptor guidance and control errors are not considered
The interceptor uses unpowered glide in the midcourse guidance phase, and the dynamics model is based on the ground coordinate system. The -axis and -axis point north and east, respectively, and form the right-hand coordinate system with -axis. The dynamics of the interceptor are modeled as follows: where is the position coordinates; is the straight-line distance from the center of the earth to the interceptor; is the velocity of the interceptor relative to the Earth; is the flight path angle; is the heading angle of the relative velocity vector; , and are all scaled by , the radius of the Earth; is the acceleration of gravity at ; and is scaled by .
Dimensionless lift and drag acceleration are as follows:
where is the reference area of the interceptor; and are the lift coefficient and drag coefficient of the interceptor, respectively, which are related to the angle of attack and Mach number; is the mass of the interceptor; and is the atmosphere density expressed as where and .
2.2. Control-Affine System
Dynamic equation (1) is directly controlled by the angle of attack and bank angle, and there is no affine relationship between them, so it cannot be directly convex. A new affine system is constructed to relax the control variables [26].
First, define a new variable normalized lift coefficient where and are the lift coefficient corresponding to the maximum lift-to-drag ratio at a certain Mach, which can be obtained in the case of aerodynamic data sheets.
According to equations (5) and (6), equations (2) and (3) can be transformed into where
So far, the dynamics equation is a function of , the dynamic equation is nonconvex because of . To solve this problem, define three new control variables
where the normalized lift coefficient is related to the aerodynamic characteristics of the interceptor. Assuming that the lift coefficient in this paper is nonnegative, the upper limit of is , the value range of is as follows:
It can be obtained through Equation (8) where .
According to the defined affine variables, equation (1) can be transformed into an affine system dynamic model
For the convenience of expression, equation (12) is rewritten as the following nonlinear system:
where and
2.3. Path Constraints
Interceptors must meet path constraints during interception to ensure normal operation. Interceptor path constraints usually include heat rate , dynamic pressure , and load , which must not exceed the maximum value given by the interceptor.
The three path constraints are expressed as follows:
2.4. Boundary Constraints and Objective Function
The starting point of interceptor midcourse guidance is the starting point of trajectory optimization. Initial point constraint is mandatory equality constraint
where is the initial value of the trajectory states.
To intercept high-speed target, the interceptor needs a certain speed to ensure the successful destruction of the target in the kinetic energy interception. The direction of the interceptor’s velocity can be enforced as an equality constraint.
where is the lower limit of interceptor speed and and are the expected values of the flight path angle and heading angle at the interceptor terminal moment , respectively.
The terminal position deviation can be constrained in the objective function. Position constraint relaxation can ensure the existence of viable solutions to the original problem. If the position constraint is directly used as the equality constraint, the original problem may not have a feasible solution, resulting in optimization failure. The result of the objective function is to find the optimal trajectory whose terminal position is closest to the expected value. where is the expected value of the terminal position and and indicate the weight coefficient, which is determined by the value and priority of the status in the performance indicator.
In summary, the original midcourse guidance trajectory optimization problem is transformed into a continuous time optimal control problem, represented as problem P0 and is represented as follows:
3. Problem of Convexification
Problem P0 is a nonconvex problem. Next, it is transformed into a convex optimization problem. Firstly, the terminal time free problem dynamic model is affined to the fixed time domain. Then, the dynamic and path constraints are linearized, and the control variable constraints are relaxed. Finally, the convex optimization problem is summarized.
3.1. Linearization of Dynamic Constraints
Time is in the variable time domain in the problem. and represent the initial and terminal time of the interceptor in problem P0, respectively. is uncertain. It can be transformed into the time domain by affine .
Applying this time domain to equation (13), the following can be obtained: where is the time expansion factor, which scales the time domain .
The reference trajectory is linearized for equation (23), and the result is as follows” where superscript represents the reference trajectory, and
The additional constraints of the trust region are expressed by the following inequality: where and are trust region constraint radius of and , respectively. They are constants large enough to limit the range of deviations from the reference trajectories defined for different systems.
Equations (19) and (20) are transformed into the time domain , as shown below:
3.2. Linearization of Path Constraints
The inequality constraints of Equations (16)–(18) are functions of , , and, which linearizes the path constraints at fixed state ; the form is as follows: where represents the maximum value of the constraint, and
When intercepting the high-speed and high maneuvering target, the interceptor usually intercepts the target in the near space with high speed. Because of the material limitations of the interceptor, it may result in a lower limit of the interceptor’s flight altitude. The lower limit of altitude is mainly determined by the speed of the interceptor, the flight altitude, and the attitude of the aircraft. Assume that the lower limit of the height of an interceptor is , then
3.3. Relaxation of Control Variable Constraints
It can be seen from Equations (8), (9), and (11) that the control variable is a strong equality constraint. The control domain formed by this equation causes the problem of non-convexity. To convexify problem P0, equations (9) and (11) can be relaxed as where is the convex domain containing .
In order to ensure the effectiveness of the control variable after relaxation, the optimal solution of the control variable satisfies equation (9). Add an additional term to the objective function . The small value is to ensure that the original objective function is not affected.
3.4. Convex Optimization Problem
A continuous time optimal control problem P0 can be transformed into a convex optimization problem after linearization of dynamic and path constraints and convexity of control variable constraints.
where equation (39) is the objective function. Equation (24) is the dynamic constraint. Equation (33) is the trust region constraint. Equations (34) and (35) are the boundary condition constraint. Equations (36) and (40) are the path constraints. Equation (41) is the control variable constraints.
4. Problem Solving Process
In this section, convex optimization problem P1 is discretized and transformed into SOCP problem P2. Then, a fast initial trajectory generation method is proposed as the initial reference solution for convex optimization problems. Finally, the process of solving SOCP problem is given.
4.1. Discretization
Problem P1 is a continuous parameter optimization problem in convex domain, which cannot be directly optimized. Therefore, problem P1 should be discretized first to facilitate computer numerical processing.
In this section, variables are discretized in the time domain and all constraints are enforced at discrete points. Assuming that the time domain is divided into discrete intervals , the total number of discrete points is , and the discrete points can be expressed as .
The dynamic constraint equation (24) is discretized by the trapezoidal method: where, , , , , , is the discrete point, .
In order to reduce the optimization space and improve the convergence speed, the variable trust region method is adopted in this paper. The method is to introduce the relaxation coefficient into the trust domain constraint and extend it into the objective function.
After discretization of problem P1, the relaxation coefficient is extended to the objective function. Then, by updating the path constraints and control variable constraints, the problem P1 is transformed into the SOCP problem P2, as follows:
4.2. Fast Trajectory Generation Method
The deviation between the initial trajectory and the optimal trajectory determines the number of optimization iterations and the efficiency of solving convex optimization problems to a certain extent. In this paper, initial trajectories are rapidly generated by interpolation and integration. The steps for initial trajectory generation are as follows.
The maximum normalized lift coefficient is selected to ensure that the initial trajectory generated quickly satisfies the trajectory constraint to the maximum extent. Within the bank angle constraint, the value of equal interval is used as the control variables.
According to the dynamic model, the control variables at different bank angles are integrated, respectively. Generate a trajectory group.
From the trajectory group, two trajectories and closest to the desired position of the terminal are selected. The distances from the desired position of the terminal are and . For the bank angles and corresponding to the two trajectories, the bank angle is selected according to
The final control variables are calculated according to the bank angle . Repeat step two to generate the initial trajectory. Finally, the maximum normalized lift coefficient and the bank angle are transformed into the control variables , , and .
Under the initial conditions in Table 1 and different terminal conditions in Table 2, according to the above fast trajectory generation method, the trajectory is generated as follows.
The flight times of interceptors at three different terminal positions , , and are , , and , respectively. The trajectories are shown in Figure 1. The maximum normalized lift coefficient is taken as the control parameter for all three trajectories, and the corresponding bank angles are shown in Figure 2. The control variables remain unchanged in the flight trajectory. The larger is, the larger the flight time is and the larger the corresponding bank angle is.
(a) Longitudinal trajectory
(b) Lateral trajectory
Compared with the GPM, the rapid trajectory generation method saves the iterative solution time and greatly improves the initial trajectory generation speed. The method takes about 0.1032 s.
The purpose of the fast trajectory generation method is to make the generated trajectory close to the optimal trajectory, reduce the number of convex optimization iterations, and improve the efficiency of problem solving.
4.3. SOCP Method Solving Process
In this section, problem P2 is solved iteratively until the discrete points converge during the solution process. This solution is the highly approximate solution of the original problem P0. The following is the solution flow of problem P2.
Step 1: set . Select and the initial time domain . The initial time expansion factor is obtained. The discrete points are selected on the initial trajectory and the control variable of the initial trajectory.
Step 2: at the iteration, the result of the iteration is used to establish the problem P2. Trajectory state variables, control variables, and time expansion factors are obtained by solving the SOCP problem.
Step 3: verify whether the iterative results meet the following convergence conditions:
If all states and control variables of all discrete points meet the convergence conditions in equation (47), step into Step 4. Otherwise, set and go back to Step 2.
Step 4: iteration stops. The convergent solution is the highly approximate solution of the problem P0.
5. Numerical Simulations
All the calculations in this paper are performed on a laptop equipped with Intel Core I7-10510 2.30 GHz, 8 G RAM, and Windows 10 operating system. ECOS [34] is used to solve SOCP problem P2 on MATLAB, where the number of discrete points was 70.
The interceptor model parameters are and . For the three path constraints, the maximum allowable values are ˙, , and . The value range of is , and the value range of bank angle is . The weight coefficients in the objective function are , , and , respectively.
The terminal constraints are shown in Table 3.
When solving problem P2, the trust region constraint parameters are as follows:
where is given in Equation (10).
The convergence conditions of simulation are given in
In Section 5.1, the feasibility of the two-stage convex optimization method is verified. In Section 5.2, the robustness of the method is verified under the initial angle perturbation condition. In Section 5.3, the robustness of the proposed method is verified by selecting different terminal position states. In Section 5.4, different initial trajectory generation methods are selected to prove the superiority of the fast trajectory generation method. In Section 5.5, compared with the GPM, the efficient solving ability of the two-stage convex optimization method is verified.
5.1. Feasibility Verification
In order to reduce the optimization space of the flight time of interceptors, the range of flight time should be constrained. The lower bound and upper bound of time are generally determined empirically. When the time optimization range is too small, the feasible solution may not be found. However, if the interval is too large, the trajectory may not converge quickly.
In Figure 3, the convergence of the error between the trajectory terminal position and the desired position are shown. The error converges in the terminal position during the first iteration. This is caused by position error which has great influence on objective function and variable trust region method. The fluctuation of position error is caused by the fluctuation of relaxation coefficient and regularization term . The purpose of subsequent iterations is to make the trajectory converge so as to ensure the validity of the iterative solution.
In Figure 4, the position error for each iteration is shown. After each iteration, the trajectory gradually converges. The results show that the third and fourth iterations are convergent.
(a) Longitudinal trajectories during iteration
(b) Lateral trajectories during iteration
In Figure 5, path constraints for convergent trajectories are shown. Heating rate, dynamic pressure, and overload all meet constraints requirements. There are peaks between 150 s and 200 s. This is because the trajectory height of convergence presents a concave curve, and the air density curve presents a convex curve, which leads to the path constraint presenting a convex curve.
Figure 6 shows the control variables of the convergence trajectory. The values of and can be transformed by and bank angle , both of which satisfy the constraints.
The flight time and operation time of four iterations are shown in Table 4. In each iteration, ECOS consumed approximately 0.4072 s of CPU to solve a SOCP problem. The total computing time is the computing time of four iterations plus the generation time of initial trajectory, which is about 1.732 seconds.
5.2. Initial Angle Perturbation
In this section, the robustness of the proposed method is analyzed by using initial angle perturbation, which provides a reference for the initial angle design of the interceptor. The contrastive trajectory is the convergence trajectory of Section 4.1. The initial angle variation condition (1) is the initial flight path angle . Variation condition (2) is heading angle .
In Figure 7, the convergence trajectories under initial angular perturbation are shown. When the initial flight path angle increases, the initial velocity component of the interceptor on the -axis increases, resulting in an increase in the height peak of the trajectory. At the same time, the initial velocity component of the interceptor along the -axis decreases, resulting in a more concave trajectory pattern. When the initial heading angle increases, the initial velocity component of the interceptor is opposite to that when the initial heading angle increases.
(a) Longitudinal trajectories under initial angle perturbation
(b) Lateral trajectories under initial angle perturbation
In Figures 8 and 9, the control variables of the convergent trajectories satisfy the constraints. Flight path angle and heading angle reach the terminal angle constraints from different angles. To some extent, the robustness of the proposed optimization method is proved. In order to obtain the cost time in solving the midcourse guidance problem, 100 Monte Carlo simulations were performed with an average cost time of 1.730 s.
5.3. Terminal Position Change
In the trajectory generation stage, interceptors need to generate optimized trajectory groups to deal with the attacks of target at different positions. This section selects three representative terminal positions ,, and . For other boundary condition constraints, refer to Tables 1 and 3.
Figure 10 shows the convergence trajectories under constraints of different terminal positions. The trajectory lines are concave in Figure 10(b). This is due to a change in heading angle from to and the increasing projection of the interceptor’s velocity on the -axis.
(a) Longitudinal trajectories at different terminal positions
(b) Lateral trajectories at different terminal positions
Figure 11 shows the variation of control variables under the constraints of different terminal positions. When the terminal position is , the flight path angle is from to negative and then to , because the bank angle is first negative and then positive. This results in a negative and then positive projection of the interceptor’s velocity on the -axis. The convergence trajectory goes from 0 km to negative and then to 0 km on the -axis.
Control variables of convergent trajectories under different terminal positions all meet the constraint conditions in Figure 12.
In this section, when the control variables meet the constraints, the terminal positions of the trajectories reach the terminal position constraints. To some extent, the robustness of the proposed optimization method is proved.
5.4. Different Initial Trajectories
To demonstrate the advantages of the fast trajectory generation method in this paper, it is compared with the method in reference [31]. Except terminal position constraint , other boundary condition constraints refer to Table 1 and Table 3. Numerical simulation is as follows.
It can be seen from Table 5 and Figure 13 that the optimization iteration under the initial trajectory one condition reaches the convergence condition after 4 times, and it takes 5 times under the initial trajectory two condition, consuming an extra time of about 0.2974 seconds. The reason is that the initial trajectory two deviates more from the convergence trajectory relative to the initial trajectory one, and more iterations are needed to approach the convergence trajectory.
(a) Longitudinal trajectories below the initial trajectory one
(b) Longitudinal trajectories below the initial trajectory two
(c) Lateral trajectories below the initial trajectory one
(d) Lateral trajectories below the initial trajectory two
The control variables of convergent trajectories under two different initial trajectory conditions, and the control variables all satisfy constraints in Figure 14.
The generation method of initial trajectory should pay attention to two aspects: (a)The rate at which the initial trajectory is generated. If the generation time is long, it is contrary to the high efficiency of convex optimization method(b)The deviation between the generated initial trajectory and the convergent trajectory. If the deviation degree is large, the solving process needs more iterations and takes more operation time
A good initial trajectory generation method should satisfy the need to generate an initial trajectory close to convergence in a short time.
5.5. Different Optimization Methods
In order to demonstrate the high efficiency of the proposed method, the GPOPS software was used to solve the GPM under the condition of the same calculation accuracy of the objective function. The comparison results are as follows.
From Table 6, both the two-order convex optimization method and the GPOPS can obtain the midcourse guidance trajectory satisfying the constraints, and the GPM has higher accuracy. Compared with the GPOPS, the convex optimization method has shorter computing time. The reason is that CVX adopts “warm-start” initial guesses and is insensitive to initial values. Compared with the CVX, GPOPS adopts interpolation “cold-start” initial guesses, which is more sensitive to initial value and with a lower computational efficiency [24].
It can be seen from Figure 15 that the optimization trajectory generated by the two-stage convex optimization method and GPOPS are not exactly the same. This is mainly caused by the differences in the ways of internal matching points and initial trajectory generation.
(a) Longitudinal trajectories under different optimization methods
(b) Lateral trajectories under different optimization methods
Figure 16 shows the changes of control variables under the two methods. All the control variables meet the constraints. The control variable of the convex optimization method fluctuates greatly, which may be caused by the scarcity of matching points and fewer iterations.
In conclusion, for midcourse guidance trajectory optimization, the accuracy of the proposed method is slightly lower than that of GPM. But the running speed of the method is greatly improved, which can provide help for rapid trajectory generation and real-time trajectory correction.
6. Conclusions
In this paper, a two-stage convex optimization method for trajectory generation of high-speed target interception with strong nonlinearity and strong constraints is proposed. By constructing affine system to linearize, relax, and discretize the original nonconvex problem, a trajectory convex optimization problem that is easier to solve is established, which is convenient for engineering application. A fast initial trajectory generation method is designed, which not only ensures the feasibility of convergent solution but also improves the convergence efficiency and solving speed of convex optimization method.
The influence of initial angle and terminal position constraints on the optimization trajectory is analyzed, which provides guidance for trajectory offline design. Compared with the GPM, the proposed method has a higher efficiency in solving the problem of trajectory generation with strong nonlinearity and strong constraints for midcourse guidance. It has potential for engineering applications and online trajectory generation.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This research was supported by the National Natural Science Foundation of China (Grant no. 62173339 and Grant no. 61873278).