Abstract
The trajectory optimization design of interceptor is very important in the defense combat against a class of high-speed and strong-maneuvering aircraft. At present, an important idea is to use hp pseudospectral method to offline optimizes the trajectory of interceptor, but the solving efficiency of this method needs to be further improved. Aiming at this issue, an improved adaptive hp pseudospectral method is proposed in this paper. In order to shorten the solving time of the algorithm, the proposed method has two main improvements on the basis of the traditional hp pseudospectral method: on one hand, by judging the positions of the control sudden change points, prerefine the mesh around them according to certain rules. On the other hand, the curvature of the system state curve is used as the criterion to segment the original mesh nonuniformly so that more mesh points can be allocated where the curvature is large. These two points together ensure that the collocation point resources can be used more efficiently in the mesh refinement process. The simulation results show that the proposed method can solve the optimized trajectory of interceptor effectively, and it also proves that this method has higher solving efficiency than the traditional adaptive hp method.
1. Introduction
Recently, with the rapid development of aerospace science and technology, a class of high-speed and strong-maneuvering vehicle such as near space hypersonic vehicle and high-performance fighter aircraft has gradually entered people’s vision. Compared with traditional aircraft, they have obvious advantages such as faster flight speed and stronger penetration capability, which pose great challenges to the country’s air and space defense operations [1–4]. Therefore, it is urgent to carry out research on intercepting technology against such targets [5].
The defense against such high-speed and strong-maneuvering targets belongs to the super-speed and super-range interception issue. In order to intercept such targets better, the range and speed of interceptor are the key breakthrough, and the separation design of trajectory and guidance law is the fundamental way to improve the range and speed of interceptor [6]. In view of this, a kind of high-speed and long-range interceptor should be selected to defend such incoming targets, and the trajectory of interceptor should be rationally optimized. As for the trajectory optimization of interceptor, many scholars have actively carried out research and achieved some results.
The trajectory optimization problem is essentially a kind of optimal control problem. However, due to the complexity of the interceptor model and many constraints in the trajectory optimization process, it is difficult to get the analytical solution of the optimized trajectory directly. In order to solve this problem, pseudospectral method [7–14] is an important idea. The method firstly selects the corresponding collocation points (such as Legendre-Gauss (LG) points [7], Legendre-Gauss-Radau (LGR) points [11, 12], or Legendre-Gauss-Lobatto (LGL) points [10]) in the time domain. Then, based on the selected collocation points, Lagrange interpolation is used to transform the original optimal control problem into a nonlinear programming problem (NLP). Finally, NLP solver [15] is used to solve the problem. Unfortunately, for complex optimization problems such as trajectory optimization, a satisfactory solution may not be obtained by using a fixed number of collocation points. Therefore, by changing the number and position of collocation points according to certain rules, many scholars continue to use pseudospectral method to solve the NLP problem iteratively until the expected solution accuracy is obtained. Typical methods include p method [16] and h method [17–19]. p method is mainly used for iterative optimization by increasing the number of collocation points, and it has a better solution effect for optimization problems with smooth solutions [16]. However, if the solution of the problem to be solved is not smooth, h method is needed to obtain higher solving efficiency, which is mainly employed to refine the mesh by constantly dividing the time domain [17–19].
Since h method and p method have their own advantages for different types of problems, in order to further improve the efficiency of optimal solution, researchers consider combining h method with p method and propose an adaptive mesh refinement method called hp method [20]-[21]. hp method adaptively adopts h method or p method in mesh refinement through a preset criterion, so that the advantages of both methods can be fully brought into play [20–22]. In recent years, many scholars [23–30] have further enriched and developed adaptive mesh refinement method along the research ideas of Ref. [20-22]. According to the posteriori error analysis, Ref. [23] proposes a new adaptive hp-refinement strategy to solve discrete local problems. Based on Birkhoff interpolation, Ref. [24] proposes a fast mesh refinement strategy to generate well-conditioned PS optimal control discretizations. For nonsmooth control problem, a mesh adaptation method is proposed in Ref. [25]. In addition, Ref. [26] introduces a C++ software program for solving the multiple-phase optimal control problems, which is called CGPOPS. Ref. [27] proposes a new adaptive mesh refinement method, which is based on the decay rate of Legendre polynomial coefficients to determine whether to use h method or p method for mesh refinement, and reduce the mesh size if necessary. In order to shorten the operation time of the algorithm, Ref. [28] refines the mesh with the arc length of the state estimation curve as the criterion, which improves the efficiency of the algorithm. By adopting the mesh density function to generate the grid, Ref. [29] proposes an adaptive control parameterization method to solve the optimal control problem. Ref. [30] designs a simple and effective adaptive mesh refinement algorithm based on interpolation error analysis and improved data compression technology, and verifies the effectiveness of the method in solving optimization problems through simulation.
Recently, many different optimization methods have been applied to solve trajectory optimization problems, and there are many meaningful research results [31–38]. Based on a series of LGL collocation points, the trajectory planning problem of interceptor is transformed into NLP, and Ref. [31] combines direct transcription and mesh refinement, and then proposes a new adaptive hp method. Aiming at the trajectory optimization problem of hypersonic reentry vehicle, Refs. [32, 33] effectively solve the trajectory optimization problem based on a multiphase Gauss pseudospectral method and a new adaptive Gauss pseudospectral method, respectively. Refs. [34, 35], respectively, explore the trajectory optimization of solar-powered aircraft and Mars entry by using the adaptive mesh refinement pseudospectral method. Ref. [36] combines Gauss pseudospectral method with Finite Fourier Series and proposes a feasible close-range satellite cooperative trajectory planning scheme. For space manoeuvre vehicles skip trajectory optimization problem, a fuzzy physical programming method, and a new three-layer-hybrid optimal control solver are proposed in Refs. [37, 38], respectively.
The above studies on aircraft trajectory optimization have achieved considerable results for their respective research concerns. Considering the specific defense operations background, the high-speed and long-range interceptor will have a period of time for offline nominal trajectory optimization before launching. Based on this, hp pseudospectral method is usually adopted in trajectory design of interceptor. In this paper, the existing hp pseudospectral method will be ameliorated to improve the optimization efficiency of interceptor trajectory.
This paper mainly improves the mesh refinement process of hp pseudospectral method, so as to enhance the solving efficiency of the algorithm. Although many scholars have made some achievements in interceptor trajectory optimization based on hp pseudospectral method, further considerations are still needed in the following aspects: (1)Most of the existing studies [31–38] do not take into account deeply the terminal constraints that the interceptor should meet at the midcourse guidance in the trajectory optimization design of interceptor, and the optimized trajectory cannot ensure that the target can be detected when the interceptor seeker is turned on(2)When pseudospectral method is used to optimize interceptor trajectory, most studies [23–30] do not pay attention to the sudden changes of optimization control in the process of mesh refinement, which will cause a large number of collocation points to be configured near the control sudden change points, thus increasing the number of mesh iterations and affecting optimization efficiency(3)In the process of trajectory optimization, the existing researches [23–30] basically insert new mesh points at equal intervals during the piecewise iteration of the mesh, which will cause the waste of collocation points resources and reduce the solving efficiency
In view of the above considerations, and draw on the ideas of [39], this paper proposes an improved adaptive hp pseudospectral method to solve the trajectory optimization of interceptor based on the combat scenario of high-speed and long-range interceptor defending high-speed strong-maneuvering target. The main innovations are as follows: (1)The terminal constraint of interceptor in midcourse guidance are considered and derived. On this basis, the optimized trajectory can ensure that the seeker of interceptor can successfully detect the target and complete the handover between midcourse and terminal guidance smoothly(2)By judging the position of the control sudden change points, the mesh is refined in advance around it. In this way, the collocation points resources are preferentially configured to reduce the number of mesh iterations(3)The nonuniform mesh segmentation iteration is carried out based on the curvature of the system state curve, so that more mesh points can be allocated where the curvature is large, so as to further improve the efficiency of solving trajectory optimization
The rest of the paper is organized as follows. In Section 2, the trajectory scheme of interceptor is designed and the mathematical model of high-speed and long-range interceptor is established. Section 3 describes the trajectory optimization problem to be solved in this paper. In Section 4, the improved adaptive hp pseudospectral method proposed in this paper is described in detail. In Section 5, the effectiveness and superiority of the proposed algorithm in solving trajectory optimization problems are verified by simulation. Section 6 summarizes the whole work.
2. High-Speed and Long-Range Interceptor Model
2.1. Interceptor Trajectory Scheme for High-Speed Strong-Maneuvering Target
Because the velocity ratio between missile and target will directly affect the guidance accuracy of interceptor, in order to effectively defend against a class of high-speed and strong-maneuvering target (hereinafter referred to as “target”), the interceptor needs to have a high interception speed. At the same time, because the target’s great flight speed will further compress the interception window of the defense side, in order to push the interception window forward to increase the interception margin of the defense side, the interceptor needs to have a larger range. As we know, an effective way to improve the range and speed of the interceptor is to adopt a hill-turn-reentry trajectory, as shown in Figure 1.
In order to achieve this type of interceptor trajectory, the high-speed and long-range interceptor adopted in this paper consists of three parts: the first stage booster, the second stage booster, and the third stage interceptor system. The third stage interceptor system mainly includes the following: warhead, fairing, detection system, missile-borne computer, and engine which can switch on and off at any time. As shown in Figure 1, at the end of the first and second boost phases, the interceptor separates from the first and second booster, respectively. The first and second stage rockets boost the interceptor rapidly through the lower dense atmosphere, which greatly reduces the velocity loss of the interceptor due to air resistance. At the same time, under the continuous boost of the two stage rocket, the third stage interceptor can be sent to the thin air in the upper atmosphere, and make it have a large flight speed, which is also an important guarantee for the high-speed and long-range interception.
Because the target has the characteristics of fast flight speed and strong maneuverability, in order to ensure that the planned trajectory can guide the interceptor to the target, it is necessary to use the correlation prediction algorithm to predict and estimate the target trajectory and get the high probability region (HPR) of the target. At the termination of midcourse guidance, it is necessary to make the detection field of the seeker successfully cover the HPR of the target, that is, to successfully capture the target, so as to realize the handover between midcourse and terminal guidance.
2.2. Mathematical Model of Interceptor
In this paper, it is assumed that the seeker of the third stage interceptor is placed on its axis, that is to say, the detection direction of the seeker depends on the attitude of the interceptor. Therefore, in order to make the interceptor successfully complete the handover between midcourse and terminal guidance, the terminal attitude of the interceptor at the midcourse guidance needs to meet certain constraints. In view of this, in order to facilitate the subsequent analysis, the pitch angle and yaw angle (roll angle is assumed to be 0) and the switch engine are taken as the control input. By considering the force applied to the interceptor and the definition of the relevant coordinate system [40], the interceptor model in the ground coordinate system is established as Equation. (1).where is the position vector of the interceptor in the ground coordinate system, represents the mass of the interceptor, is the gravity constant, is the engine thrust, and represents the specific impulse of the engine. , , and are components of velocity vector of interceptor in three axes of ground coordinate system, respectively.
Velocity value of interceptor can be expressed as . Flight path angle and heading angle can be calculated by , , and .
According to the geometric relation equation of angle, angle of sideslip , angle of attack and bank angle can be calculated from pitch angle , yaw angle , flight path angle , and heading angle .
Based on the solid rocket model in [41] and the background of defense operations, this paper makes reasonable improvements, and designs the relevant performance parameters and aerodynamic parameters of the interceptor as Table 1.
During the flight of interceptor, drag force , lift force , and side force are expressed as where is the air density, is the air density at sea level, represents the reference height, and represents the height of the interceptor from sea level. , , and are drag force coefficient, lift force coefficient and side force coefficient, respectively. where Ma stands for Mach number, represents the partial derivatives of with respect to , and represents the partial derivatives of with respect to .
In the midcourse guidance stage, the aerodynamic parameters of the third stage interceptor are
In addition, in Equation (1) represents the switching amount of the engine. During the first and second boost stages, the rocket engine is continuously turned on, so there is . In the midcourse guidance stage, the engine of the third stage interceptor can be switched on and off at any time, and when it is turned on, there is . When the engine is turned off, there is .
Since is defined as discrete value (equal to 0 or 1) in the midcourse guidance stage, while it is necessary to ensure that the input control is continuous value in the subsequent trajectory optimization process. Therefore, pulse-width pulse-frequency modulator (PWPF) [42] is introduced here to modulate the continuous control to be optimized into discrete switching instruction . Its schematic diagram is shown in Figure 2.
where the gain of the first-order plant is and the time constant is . The turn-on threshold for relays is and the turn-off threshold is . is the switching control to be optimized, which is continuous value. is always 1 in the first and second boost stages, and there is in the midcourse guidance stage.
3. Trajectory Optimization Problem Description
Combined with the above model, the first and second boost stages are programmed guidance stage, and the principle of trajectory optimization in this stage is to send the interceptor to a certain altitude and make it have the maximum flight speed. While for midcourse guidance, the principle of trajectory optimization is to make the third stage interceptor smoothly complete the handover between midcourse and terminal guidance. At the same time, the fuel consumption of interceptor engine in the midcourse guidance should be reduced as much as possible to ensure that the third stage interceptor can have a large interception speed and strong maneuverability in the terminal guidance stage. According to the above principles, it is necessary to solve the trajectory optimization in these two stages, respectively. Firstly, the trajectory optimization problem is described as follows. where is the cost function, is the state equation of the system, represents the terminal constraint, and is the process constraint. is the system state variable, is the control variable to be optimized, and and represent terminal time and terminal state, respectively.
3.1. Trajectory Constraints in the Programmed Guidance Stage
For programmed guidance stage (first and second boost stages), the cost function is to maximize the velocity at the terminal of the second boost stage, i.e.,
The state equation is shown in Equation (1). The terminal constraint is that the flight path angle and heading angle of the interceptor reach the expected value at the terminal of the second boost stage, i.e., where and represent the flight path angle and heading angle at the terminal moment, respectively, and represent the desired flight path angle and heading angle, respectively.
In addition, during the programmed guidance stage, the interceptor needs to meet the process constraints of dynamic pressure , heat flux density , overload , i.e., where , , and represent upper limits of , and , respectively, and where is the heat flux constant and is the vector sum of thrust, drag, lift, and side forces.
3.2. Trajectory Constraints in the Midcourse Guidance Stage
For the midcourse guidance stage, the cost function is to minimize the fuel consumption of the third stage interceptor during flight, that is, to maximize the remaining mass of the third stage interceptor at the end of the midcourse guidance, i.e.
The state equation of the interceptor in this flight phase is shown in Equation (1). In order to achieve the handover between midcourse and terminal guidance, it is necessary to ensure that the third stage interceptor can detect the target at the end of the midcourse guidance stage, that is, to ensure that the detection range of the third stage interceptor’s seeker can successfully cover the HPR of the target, which is displayed in Figure 3.
As shown in Figure 3, it is assumed that the HPR of the target is a spherical region with radius in space, the coordinate of its spherical center is , the position of the interceptor’s centroid is , and the distance between and is . In order to ensure that the seeker’s detection field can cover the target’s HPR more effectively, it is necessary to make the axis direction of the seeker’s detection field consistent with the line direction of , that is, to make the pitch angle and yaw angle of the interceptor meet the following geometric relations.
Furthermore, assuming that the detection field angle of the seeker is and the maximum detection range is , the following constraints should be met to ensure that the detection field can completely cover the HPR of the target.
Equations (22) and (23) are the terminal constraints of interceptor in the midcourse guidance stage. In addition, the third interceptor also needs to meet various process constraints in the midcourse guidance stage, as shown in Equations (19) and (20).
4. Trajectory Optimization Based on Improved Adaptive Hp Algorithm
In order to intercept the target successfully, the current interceptor mostly adopts the compound guidance mode. Under the guidance system, the trajectories of programmed guidance stage and midcourse guidance stage are optimized offline in advance combining with the prediction of target position. Since the trajectories of interceptor in these two stages do not need to be acquired online, the trajectories can be solved by offline optimization algorithm, and hp pseudospectral method is commonly used to solve this problem. In order to improve the solving efficiency of trajectory, this paper proposes an improved adaptive hp pseudospectral method based on the traditional hp method, which mainly improves the mesh refinement process in two aspects: (1) By evaluating the change rate of the control, the grid is prerefined near the sudden change points of the control according to certain rules. (2) According to the curvature value of the state curve, the mesh interval is divided nonuniformly to ensure that more mesh points are allocated to the position with large curvature value.
4.1. Optimization Problem Transformation Based on LGR Collocation Method
In order to solve the trajectory optimization problem (16), the idea of hp pseudospectral method [20–22] is adopted here. Firstly, mesh points are selected on the whole time interval , and then is divided into mesh intervals , where , , and . At the same time, let and represent the state variable and control variable on the mesh interval , so the original optimization problem (16) can be transformed into the following equivalent problem B.
It is worth noting that since the state variable is continuous in time interval , each mesh point in Equation (24) needs to meet the following conditions.
Then problem B is transformed into NLP by LGR collocation method. First, we need to convert the time variable in each mesh interval to interval . Here, we introduce a new time variable to perform the following transformation.
Let represent the length of each mesh interval , i.e., . Then, combining Equations (24)–(26), problem B can be equivalent to the following problem.
Next, LGR collocation points are selected on each mesh interval . In order to discretize Equation (27), the selected LGR collocation points and noncollocation points are used as sampling points to perform Lagrange interpolation for state variable on interval . The Lagrange interpolation polynomial based on the sampling points is expressed as follows
Then, continuous state variable on mesh interval can be expressed as where is the approximate value of . Take the derivative of Equation (29) at the collocation points , then where . Next, let represent the approximate value of the control variable . Combining Equations. (27)–(30), problem B can be transformed into the following NLP, denoted as problem
For problem , the NLP solver SNOPT [15] can be used to solve it. However, solving it only once may not achieve the solution accuracy required by the original problem B. Therefore, it is necessary to continuously adjust the position of mesh points (h method) or the number of LGR collocation points in each mesh interval (p method) to solve it iteratively until the required solution accuracy is achieved. In view of this, this paper proposes an improved adaptive hp mesh refinement method.
4.2. Improved Adaptive Hp Mesh Refinement Method
When solving the trajectory optimization problem in this paper, the obtained control variables may suffer from sudden change. However, the traditional hp method [20–22] will allocate a large number of LGR collocation points near the sudden change points or segment the grid frequently, which will increase the solving dimension of NLP and affect the solving efficiency. In addition, the traditional hp method often performs average mesh segmentation in the mesh refinement process, and the excessive curvature of state curve will further increase the solving burden of hp method. To solve the above problems, an improved adaptive hp mesh refinement method is proposed in this paper, and the detailed steps are as follows.
4.2.1. Evaluation of Solution Error
When solving problem with NLP solver, it is necessary to determine whether the obtained solution can make the error between and in Equation (29) within accuracy tolerance range . When it exceeds , the original mesh needs to be refined and the optimization problem needs to be solved again. Here we adopt the following method to evaluate the solution error in each mesh interval.
First, an LGR collocation point is added in each mesh interval. In the mesh interval , the new collocation points are denoted as , where and . At the same time, take the noncollocation point .
Then, the state variable at and the control variable at are obtained by Lagrange interpolation combined with the obtained numerical solution and the corresponding sampling point . These two can be expressed as
Finally, combining and obtained above, the state variable at is reintegrated by using the integral form of the discrete state equation.
So far, we can use to approximate the error between numerical solution and analytical solution of problem B, as shown below where and represent the th element of vectors and , respectively.
4.2.2. The Mesh Prerefinement around Control Sudden Change Points
When hp pseudospectral method is used to solve problem iteratively, the sudden change of control variables will greatly increase the number of iterations and affect the efficiency of solving. In order to effectively solve the above problems in the improved adaptive hp method in this paper, we prerefine the mesh based on the control sudden change points. The basic idea is to accurately locate the control sudden change points, and then prioritize the mesh refinement near these points to improve the solving efficiency.
First, in order to locate the sudden change points, the first derivative of time is taken at with respect to the control variable . where . Equation (35) can be further written as
Then, according to the derivative of the control variable and the set threshold , the interval where the control variable changes too fast is selected. Assume that the derivatives of control variables at each collocation point in interval are shown in the figure below.
As shown in Figure 4, when the derivative of the control variable exceeds (such as and in the figure), we believe that there may be control sudden change points in this area, so we need to refine the mesh in advance. We will add new mesh points to further refine region and . In order to describe the subsequent mesh refinement process, the new mesh points combination obtained is denoted as and the new mesh interval is denoted as .
4.2.3. Mesh Refinement Classification Criterion Based On State Variable Curvature
When hp method is used to solve the optimal control problem, different methods should be used to refine the mesh according to the characteristics of different mesh [31]. In some mesh intervals, the state variables change dramatically, which is suitable for segmentation, that is, h method; For those mesh regions with relatively stable state changes, it is appropriate to refine them by adding collocation points, namely p method. In the improved adaptive hp method in this paper, we take the curvature of state variables as the criterion of classification.
The set of collocation points in the new mesh interval is denoted as , and the state variables at each collocation points are expressed as . Based on and , points with equal time interval in interval are selected for Lagrange interpolation, and the interpolation points obtained are marked as , where . The curvature of each point can be calculated by using the coordinates of interpolation points
Combined with the curvature obtained, the classification refinement criterion is set as
If the curvature threshold parameter is set as , and if there is , it can be considered that the state variable changes too fast in the mesh interval and needs to be refined by segmenting. On the contrary, if there is , it is necessary to refine the mesh by adding collocation points.
4.2.4. Increasing Collocation Points Based on Maximum Approximation Error
According to the above judgment, if mesh interval needs to further increase the number of collocation points, based on the maximum approximation error in Equation (34), we increase the number of collocation points according to the following equation where is the original number of collocation points in mesh interval , represents the number of collocation points after the increase, and is the set accuracy tolerance. is the least integer function.
4.2.5. Mesh Segmentation Method Based on Cumulative Sum of Curvature Values
In the traditional adaptive hp method, the mesh interval that needs to be segmented is usually evenly divided into several subintervals, which do not maximize the efficiency of segmentation. In order to further improve the efficiency of mesh refinement, we design a segmentation method based on the accumulative sum of curvature value, which can allocate more mesh points at the position with larger curvature of state curve.
First, estimate the number of subintervals that need to be added. The number of subintervals segmented within mesh interval is determined by the following equation. where represents the maximum value of and , and is the set product factor.
Then, based on the curvature of state curve in mesh interval , the position of new mesh points is determined. According to Equation (40), interval needs to be divided into subintervals, that is, mesh points should be added between and . In order to determine the positions of these mesh points, the curvature value at each interpolation point in Equation (37) needs to be taken into account to ensure that more dense mesh points are configured at positions where the state variables change dramatically.
According to the above principles, the mesh point configuration rule of the improved adaptive hp method in this paper is as follows: points are selected from interpolation points in interval , denoted as , so that the curvature sum of interpolation points in each subinterval is approximately equal. According to the rule, the selection algorithm of newly added mesh points is as follows: where .
4.2.6. Algorithm Flow
At this point, the key steps of the improved adaptive hp method proposed in this paper have been explained, and the description of the algorithm can be summarized as Algorithm 1.
| ||||||||||||||||||
5. Simulation Results
Next, based on the interceptor model in Section 2 and the trajectory optimization problem in Section 3, the improved adaptive hp method proposed in this paper is used to conduct trajectory optimization simulation. In this part, all the simulation relies on the computer hardware of Lenovo desktop computer, whose processor is Intel(R) Core(TM) i5-4570, main frequency is 3.20GHz, memory is 4.00G, and the software platform of simulation is Matlab2014a. In addition, the algorithm proposed in this paper is developed based on the open-source hp pseudospectral optimization software GPOPS [43], and the NLP solver used is SNOPT [15].
5.1. Trajectory Optimization Simulation of Programmed Guidance Stage
Firstly, we use the improved adaptive hp algorithm proposed in this paper to optimize the trajectory of the programmed guidance stage (first and second boost stages). The initial state of interceptor in programmed guidance stage is , , and . Its initial control variable is . It is worth noting that since it is stipulated that the engine of the interceptor is always on in the programmed guidance stage, there is . Meanwhile, the upper limit of process constraint of interceptor in programmed guidance stage is set as , , and .
In the improved adaptive hp algorithm proposed in this paper, there are four parameters to be designed:, , , and . The smaller the accuracy tolerance range is, the higher the accuracy is required. Generally, satisfactory accuracy is achieved when is less than 10-3. The larger the product factor is, the more mesh points will be added during grid iteration, which will improve the accuracy of numerical solution to a certain extent. However, too large will lead to a large number of new mesh points, thus greatly increasing the dimension of NLP problem and affecting the solving speed. Generally, the value of can be selected between 1 and 2. The smaller threshold is, the more sensitive the proposed algorithm is to the change of control, that is, the smaller the change of control will make the grid to be prerefined, and the value between 0.5 and 2 is appropriate. The smaller the curvature threshold is, the easier it is to add mesh points instead of collocation points during grid iteration. Generally, it can be set at 1 ~ 3.
According to the above parameter selection principle, relevant parameters in the improved adaptive hp algorithm are selected as follows: , , ,and .
Usually, when setting the initial grid, the grid interval is not segmented, that is, usually only two initial mesh points are selected, one at the initial moment and one at the terminal moment. In addition, although the approximation accuracy of numerical solution can be improved if there are too many initial collocation points in the entire initial grid interval, the solving speed of the algorithm will be seriously affected. Generally, it is appropriate to choose the number of LGR collocation points between 8 and 15. According to this principle, the initial mesh point is set as , and the initial LGR collocation points are set as 12.
In addition, under the same mesh accuracy tolerance range , hp method in [21] is used to solve the problem, and the simulation is compared with the method in this paper. In the simulation results, P1 represents the first boost stage, and P2 represents the second boost stage. As described in Section 3.1, the interceptor also needs to meet the constraints of flight path angle and heading angle at the end of the programmed guidance stage. Here, the expected value of the two is set as . Simulation results are shown in Figure 5 and Table 2.
(a) Trajectory of the missile
(b) Velocity vs. time
(c) Mass of missile
(d) Flight path angle and heading angle
(e) Control vs. time
(f) Process constraint
(g) Various grids of improved hp
(h) Various grids of hp
As can be seen from Figures 5(a)–5(f), when the improved adaptive hp method proposed in this paper and the hp method in [21] are used to optimize the trajectory of the programmed guidance stage, respectively, the optimization results obtained by the two methods are basically the same. As shown in Figure 5(d), at the end of the programmed guidance stage, the interceptor’s flight path angle converges to 15 degrees and the heading angle is stable at 0 degree, that is, the optimized trajectory obtained by the two methods can meet the preset terminal constraints. Meanwhile, Figure 5(f) shows that both methods can make the interceptor meet the preset process constraints in the flight process of the programmed guidance stage.
However, it is worth noting that in this simulation, in order to achieve the same optimization solution accuracy, the hp method needs 4 mesh iterations, while the improved adaptive hp method only needs 3 iterations (see Figures 5(g) and 5(h)). By observing Figures 5(g) and 5(h), it can be found that the original grid is segmented uniformly by the traditional hp pseudospectral method during grid iteration. However, the improved hp pseudospectral method proposed in this paper is different from this method. According to the curvature of the state curve, this method divides the original grid into nonuniform segments, so as to ensure that the position with large curvature can be allocated to more LGR points more quickly, and then ensure that the obtained solution can meet the set accuracy more quickly. This is more concretely shown in Table 2. As shown in Table 2, the terminal velocity obtained by trajectory optimization using these methods is basically the same (obtained by the improved adaptive hp method slightly higher). Furthermore, in order to achieve the optimization result, the hp method needs to carry out 4 mesh iterations and configure 111 LGR points, which takes 25.48 s. However, the improved adaptive hp method only needs to carry out 3 iterations and configure 87 LGR collocation points, which takes 15.66 s. The results show that the optimization results obtained by the two methods are basically the same, while the improved hp adaptive method takes less time and has higher efficiency.
In addition, the following explanations are needed. The defensive operation background considered in this paper is ultra-long-range interception. When low-orbit satellites or high-power radars find the target, it is still far away from the defensive position and the remaining attack time is long. In order to intercept such targets effectively, the trajectory of the target should be predicted first, and then the possible position of the target in the future can be obtained. Before launching interceptor, it is necessary to offline optimize the trajectories of programmed guidance stage and midcourse guidance stage based on the predicted target position, and then get the nominal trajectory of interceptor, and the algorithm proposed in this paper is mainly for this purpose. Therefore, in this combat scenario, the time for trajectory optimization before launching of interceptor is relatively long, and the CPU time in simulation is within the acceptable range.
5.2. Trajectory Optimization Simulation of Midcourse Guidance Stage
The terminal state of the interceptor in Section 5.1 is taken as the initial state of the interceptor in the midcourse guidance stage, and the improved adaptive hp method proposed in this paper is also used for the midcourse guidance trajectory optimization simulation. Set the center position of the HPR of target as (800 km, 30 km, and 35 km) and radius as 6 km. The detection field angle of the seeker is set as , and its maximum detection range is set as 100 km. In the simulation, the relevant parameters of the improved adaptive hp algorithm are set as: , , , and . The initial mesh point is , and the number of initial LGR collocation points is 10. In addition, the process constraint upper limit of interceptor is set in accordance with Section 5.1. Similarly, the hp method in [21] is used for trajectory optimization, and the simulation is compared with the method in this paper.
Meanwhile, the position and attitude data of the interceptor at the terminal of midcourse guidance are summarized as follows:
According to the terminal state of interceptor in Table 3, combined with the HPR of target set in the simulation and related parameters of seeker, the simulation is carried out as follows to verify whether the optimized trajectory obtained by the two methods can ensure that the seeker of interceptor successfully detect the target.
As can be seen from the above simulation results, when the two methods are used for trajectory optimization simulation of midcourse guidance stage, the obtained optimization results are basically the same (see Figures 6(a), 6(b), and 6(d)). As shown in Figure 6(f), when the interceptor flies along the optimized trajectory obtained by the two methods, the dynamic pressure, heat flux density, and overload can all meet the process constraints. In addition, at the end of midcourse guidance stage, the state of the interceptor under the two methods can ensure that the detectable region of the seeker effectively covers the HPR of the target, thus successfully completing the handover between midcourse and terminal guidance (see Figures 7 and 8). Furthermore, the control variables of the two methods are slightly different, especially the switch engine (see Figure 6(e)), which leads to a difference in the fuel consumption of the third stage interceptor in the midcourse guidance stage. The fuel consumption of the improved hp method is less than that of the hp method (see Figure 6(c)).
(a) Trajectory of the missile
(b) Velocity vs. time
(c) Mass of missile
(d) Flight path angle and heading angle
(e) Control vs. time
(f) Process constraint
(g) Various grids of improved hp
(h) Various grids of hp
Although the optimization effects of the two methods are almost the same, it is worth noting that the efficiency of the two methods is different when solving trajectory optimization. As shown in Figures 6(g) and 6(h), in order to achieve the same optimization accuracy, the hp method requires 3 mesh iterations, while the improved adaptive hp method requires only 2 iterations and requires significantly fewer LGR collocation points. This is more concretely reflected in the comparison in Table 4. As shown in Table 4, 131 LGR collocation points need to be configured for the optimization by hp method, which takes 27.38 s, while the improved adaptive hp method only needs 61 collocation points, which takes 6.97 s. In general, in the trajectory optimization simulation of midcourse guidance stage, the two methods can get the optimized trajectory satisfying the constraints. Meanwhile, the optimization effect of the improved adaptive hp method is slightly better than that of the hp method, and its solving efficiency is significantly higher than that of the hp method.
6. Conclusion
In this paper, an improved adaptive hp pseudospectral method is proposed to solve the trajectory optimization problem of high-speed and long-range interceptor. Different from the mesh refinement process of the traditional adaptive hp method, the proposed method refines the mesh in advance around the control sudden change points. In addition, the method takes the curvature of the state variable curve as the criterion for the nonuniform mesh segmentation iteration; so that more mesh points can be allocated where the curvature is large. These two methods make the allocation of collocation points resources more reasonable and improve the optimization efficiency. The simulation results show that the optimized trajectory of the interceptor can be effectively solved by using the algorithm proposed in this paper, so that the interceptor can successfully realize the handover between midcourse and terminal guidance. And the efficiency of this method in trajectory optimization is higher than that of the traditional adaptive hp method.
Data Availability
All data generated or analysed during this study are included in this published article.
Conflicts of Interest
The authors declare no conflict of interest.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (Grant no. 61873278 and no. 61773398).