Abstract
A hybrid double-loop optimization algorithm combing particle swarm optimization (PSO) and nonintrusive polynomial chaos (NIPC) is proposed for solving the robust trajectory optimization of hypersonic glide vehicle (HGV) under uncertainties. In the outer loop, the PSO method searches globally for the robust optimal control law according to a penalized fitness function that contains the system robustness considerations. In the inner loop, uncertainty propagation of the stochastic system is performed using the NIPC method, to provide statistical moments for the iterative scheme of the PSO method in the outer loop. Only control variables are discretized, and the state constraints are satisfied implicitly through the numerical integration process, which reduces the number of decision variables as well as the huge amount of computation increased by NIPC. In the end, the robust optimal control law is achieved conveniently. Numerical simulations are carried out considering a classical time-optimal trajectory optimization problem of HGV with uncertainties in both initial states and aerodynamic coefficients. The results demonstrate the feasibility and effectiveness of the proposed method.
1. Introduction
Hypersonic glide vehicle (HGV) is generally released from solid rocket boosters and then reentry glides through the atmosphere at hypersonic speed without power. Depending on the high lift-to-drag (L/D) ratio shape design and proper aerodynamic control techniques, it can achieve long flight distance and strong maneuverability. In recent years, HGV has attracted worldwide attention for its broad application prospects in both military and civilian fields [1–3].
However, the reentry region of HGV is quite narrow and always suffers from complex uncertainties, due to large space span, long flight time, changeable aerodynamic environment, and shortage of flight experience [4]. Therefore, the reentry trajectory design, as a core problem of the guidance and control, has become a challenge and hotspot for HGV [5, 6]. The performance of trajectory design is the essential section of a flight. The trajectory design framework of offline trajectory optimization and online tracking guidance is usually adopted in engineering, because of the limitation of the current vehicle-borne computing capacity and the high real-time calculation requirements [6, 7]. Under this framework, the reference trajectory is optimized offline and then stored in the onboard computer for the guidance and control system to track.
The reference trajectory optimization process affords an overall analysis of multidisciplinary designs and provides a theoretical basis for many other disciplines, such as guidance, control, defense penetration, and thermal protection. It is essentially a nonlinear optimal control problem with multiple constraints, including the control limitation, dynamic pressure, heating rate, and aerodynamic load. It is difficult to solve this problem analytically; thus, many numerical optimization algorithms have been developed and applied by now [8–12]. According to Betts’ survey on numerical trajectory optimization methods, these algorithms can be divided into two categories [13]. The first one is the indirect method, in which the trajectory optimization problem is converted to a Hamiltonian boundary value problem using the variation method or Pontryagin’s maximum principle [11]. These methods can obtain accurate analytical solutions but their initial states are difficult to guess. The other one is the direct methods, in which the optimal control problem is transformed into a nonlinear programming problem and then solved by some numerical methods. These methods are weakly dependent on their initial states and require no system derivative information. In the field of direct methods for trajectory optimization, the metaheuristic algorithms, especially the swarm intelligence algorithms, have received widespread attention [10]. However, the deterministic optimization methods lack robustness considerations for uncertainties and, therefore, will increase the design burden of the online tracking guidance system for the case of unpredictable fluctuations or deviations from the original reference trajectory [14].
To solve this problem, trajectories that are statistically less sensitive to uncertainties are preferred, and robust reference trajectory optimization is needed. The improvement of robust trajectory optimization over the conventional trajectory design is to integrate the independent design steps into an interactive whole, by means of building a reverse information flow between the offline trajectory optimization and the online reentry environment. Compared with some online trajectory optimization and tracking guidance techniques [15, 16], robust trajectory optimization displays obvious advantages in reducing the burden of the guidance and control system effectively.
Robust trajectory optimization of HGV is a complex stochastic process with multidisciplinary attributes such as flight system modeling, uncertainty analysis, and trajectory numerical solution. Since the deterministic optimization technology is relatively mature, the commonly used robust strategy is to transform the stochastic system into an equivalent deterministic system. Generally, this is achieved by adding appropriate statistical information (such as the mean and the standard deviation) of the user-defined qualities to the original objective function as well as the constraints. To this end, the efficient and accurate uncertainty propagation (UP) techniques [17] are necessary, which usually include Monte Carlo (MC) simulation, linear methods, and nonlinear methods. MC simulation relies on a large random sample space to approach the true value and is usually used to verify the effectiveness of the target algorithm. Janson et al. proposed an MC motion planning method for the robot trajectory optimization problem under uncertainty and proved its feasibility and effectiveness [18]. Linear methods simplify the uncertainty propagation analytically but are inapplicable for highly nonlinear hypersonic vehicles. As one of the popular nonlinear methods, polynomial chaos (PC) uses the sum of orthogonal polynomial chaos expansions to fit the probabilistic model output response. It is based on the spectral representation of uncertainties and can obtain the system’s high-order statistical moments conveniently [19, 20]. The common implementation forms of PC include the generalized polynomial chaos (gPC) and the nonintrusive polynomial chaos (NIPC) that produce almost the same accuracy. However, the NIPC is more promising for it treats the system as a black box and requires no modification on the existing code [21].
In recent years, the PC method has been implemented to some aerospace engineering problems. Cottrill and Harmon [22] proved that the gPC algorithm displayed improved calculation efficiency while it maintained the same accuracy with Monte Carlo simulation. Subsequently, they proposed a robust trajectory optimization model that used the polynomial chaos method and the Gaussian pseudospectral method to transform the stochastic optimization problem into a set of similar deterministic optimization problems. Furtherly, Okamoto and Tsuchiya [23] combined the gPC method with the Legendre-Gaussian pseudospectral method to study the robust optimization of aircraft landing trajectories in a microburst environment. Boutselis et al. [24] combined the gPC method with the dynamic programming algorithm for the trajectory optimization of quadcopters with uncertain model parameters, and they adopted variational integrators to improve the performance of the proposed framework.
The robust optimization methods reviewed above are concise. However, when the target system is complicated, the repeated solution of the deterministic optimization problem will be unacceptable. Therefore, Fisher and Bhattacharya [25] transformed the stochastic optimization problem into an expanded deterministic optimization problem in the higher-dimensional state space and ensured the robustness considerations by adding optimization constraints. Afterward, Li et al. [26] conducted a preliminary study on the robust dynamic trajectory optimization of the aircraft for a small number of random variables. Huang et al. [27, 28] combined the NIPC with the most probable point search algorithm and the evidence theory, respectively, to study the robust trajectory optimization of Mars entry. Wang et al. [29] also addressed the Mars entry problem by combining the gPC with convex optimization techniques, which improved the efficiency and accuracy of the robust trajectory optimization. These methods avoid the repeated solution of corresponding deterministic optimization problems, but a larger number of constraints, especially the expanded ordinary differential equations, increase the difficulty of optimization convergence.
Consequently, how to imply the PC method to solve the complex HGV trajectory optimization problem with uncertainties efficiently is still a crucial problem. The key difficulty in the application of polynomial chaos is that it is easy to cause dimension disaster with the increase of random variables. As a result, the computational burden is unacceptable. The sparse tensor-product theory is a potential solution to this problem [30, 31]. However, if proper mathematical techniques can be used to simplify the application process of polynomial chaos in trajectory optimization, these difficulties will be alleviated in a wider scope.
In this paper, a stochastic trajectory optimization model of HGV is studied in the aspect of uncertainties in both initial states and aerodynamic coefficients. Specifically, it is assumed that the exact values of these uncertainties are unknown, but their probability density functions can be accessed. On this basis, the objective function and the constraint function are augmented to include the statistical moments of the system, and the dynamic robust optimization model is established. Then, the NIPC method is utilized to perform the uncertainty propagation over time for the statistical information and transform the original stochastic ordinary differential equations into deterministic ones. Afterward, an expanded deterministic optimal control problem with a stochastic flavor is obtained depending on the NIPC representation of the dynamics. Since the dimension of deterministic states is expanded, the control parameterization method is preferred in this paper for reducing the number of decision variables. Therefore, a novel version of the PSO method [32, 33] is incorporated to solve the obtained problem, and a hybrid double-loop optimization algorithm called NIPC-PSO is developed. In the proposed algorithm, the PSO particles which represent different control laws are iteratively optimized in the outer loop, based on their probabilistic behaviors that are calculated via NIPC in the inner loop. Moreover, the computing efficiency of the proposed framework is further improved by incorporating the concept of a numerical integration scheme for solving the complex dynamic equations.
The main contribution of this work lies in developing a concise formulation for robust reference trajectory optimization of HGV with uncertainty in dynamics. Only control variables are discretized, and the expanded state constraints are satisfied implicitly utilizing the benefits of the numerical integration scheme. The number of decision variables as well as the huge amount of computation for uncertainty propagation using NIPC is extremely reduced. These above properties of the proposed approach indicate its high efficiency and potential applicability for robust trajectory optimization problems of HGV.
The remainder of this paper is organized as follows. In Section 2, the robust trajectory optimization model is established. Then, the double-loop NIPC-PSO framework is described step by step in Section 3. In Section 4, a classical trajectory optimization problem of HGV with uncertainties is introduced. The numerical simulation examples are included in Section 5, and conclusions are drawn in Section 6.
2. Problem Formulation
The trajectory optimization is a typical nonlinear constrained optimal control problem. A commonly used optimal control formulation of the deterministic trajectory optimization is as follows:where is the control variable vector and is the time; is the objective function, in which is the cost of boundary constraints, and is the cost of Lagrange (integral) term; and denote the initial and final time, respectively; is the state variable vector and is the system ordinary differential equations, i.e., ODEs; and are the path constraints and boundary constraints, respectively.
With the consideration of uncertainties, random variables obeying some probability distribution functions (PDFs) should be included [23]. In the present study, the control variables are assumed as deterministic functions of only [23, 25], whereas the performance objective , the ODEs, and the constraints and are all stochastic functions of the random variables. The goal of robust trajectory optimization is to obtain the optimal control law, which is insensitive to possible disturbances. Therefore, when uncertainty occurs, the expected trajectory dispersion would be minimized statistically.
To achieve the goal of robust trajectory optimization, the objective and constraint functions are often augmented to include robustness considerations, usually the statistical moments of the stochastic system, such as the mean and the standard deviation. In this study, the robust trajectory optimization problem is formulated as follows:where is the random variable vector; is the augmented objective function; and are the expectation and standard deviation of , respectively; and are the weight factors of and ; denotes the stochastic state variable vector; denotes the stochastic ODEs; the variables , and , are the expectations and standard deviations of and , respectively; and are allowable standard deviation thresholds of and .
For the above robust optimization model, the deterministic demands are satisfied using mean values and the stability requirements are guaranteed with standard deviations. Robustness of both the objective function and the constraint functions is carried out simultaneously. When this model is performed for trajectory planning problems of HGV, it is necessary to solve the complex stochastic ODEs as well as system robustness considerations.
3. The Proposed Robust Trajectory Optimization Scheme
To solve the above robust trajectory optimization problem, a step-by-step description of the proposed hybrid procedure using NIPC and PSO is given as follows. Step 1. Transform the stochastic ODEs into augmented higher-dimensional deterministic ODEs using the NIPC method. For a random variable vector , the exact solution to the original stochastic ODEs of equation (2b) can be represented as the following NIPC model: where is the approximation of ; is the number of expansion terms that is determined by and is the desired approximation order of orthogonal polynomials; are the PC expansion coefficients; are the orthogonal PC basis functions belonging to the Askey family [20] and satisfy the following orthogonal properties: where is the Kronecker delta function, and denotes the inner product, defined by where and are functions of and is a weight function. There are a corresponding polynomial basis function and weight function for a random variable with a given distribution. Moreover, for the random vector , and can be obtained based on the product of corresponding functions related to random variables in each dimension. According to the NIPC expansion scheme, the unknown PC expansion coefficients of equation (3) can be calculated by the following integral [26]: Based on the above formulation, the full tensor-product numerical quadrature rule can be used to approximate , i.e., where is the number of quadrature points along each dimension of ; , are the integration points of the component of , and is the multidimensional quadrature weight corresponding to the dimension of the quadrature point . To implement the above numerical approximation of equation (7), and need to be determined firstly according to a proper quadrature rule, such as Gaussian quadrature. Then, compute and . Afterward, calculate by substituting into equation (2b). According to equation (7), there are quadrature points of , i.e., . Therefore, deterministic ODEs are generated and need to be solved, i.e., where is the augmented deterministic ODEs with . Step 2. Establish the equivalent deterministic trajectory optimization model of based on NIPC. According to the NIPC model of equation (3), the original stochastic state can be represented as . When the PC expansion coefficients have been calculated, its first two statistical moments can be easily computed as follows: where and are the mean and the standard deviations of . The original performance objective and the constraints and are all functions with respect to the approximate stochastic state . Therefore, their statistical moments can also be calculated deterministically based on the PC expansion coefficients and their specific function forms using corresponding computation rules, such as addition, subtraction, multiplication, and division [29]. Afterward, the original stochastic optimal control problem is approximately transcribed into a deterministic optimal control problem , i.e., The main difficulty in the application of the above model lies in the repeated solutions of the augmented deterministic ODEs to obtain the PC expansion coefficients. In particular, for the complex flight systems of HGV, this process requires a lot of calculation and time. Some mathematical techniques are adopted to improve this problem in this study. Step 3. Discretize the continuous optimal control problem to and use the numerical integration scheme to simplify the calculation of PC expansion coefficients. To solve the above continuous optimization problem using direct optimization algorithms, discretization strategies are necessary and the number of decision variables always affects the calculation amount of the discrete optimization process significantly. In the problem of , there are deterministic ODEs and the number of state variables is quite a lot. Thus, for reducing the number of decision variables, the control parameterization method [33] instead of the state and control parameterization method is adopted here to discretize the continuous optimization problem. The time interval from to is first divided into subintervals as follows: where is the discretized sampling time. Then, the control variables at the discrete nodes are selected as the decision variables to be optimized, that is, where is the discretized control vector. The discrete-time dynamic optimization problem is established by fewer design variables; however, the expanded ODEs are still not resolved. In this study, the numerical integration scheme is incorporated to simplify the calculation of PC expansion coefficients. For the discretized control vector , during each segment , it can be approximated by linear interpolation, i.e., The approximate continuous control variable can be represented as follows: For the convenience of description, represent the expanded ODEs of equation (8) in a brief form as follows: where denotes the deterministic states and denotes the augmented deterministic ODEs. According to the numerical integration theory, the state variables propagate from to and their values of any time can be obtained. Herein, by providing control inputs and necessary initial conditions, a brief fixed-step fourth-order Runge-Kutta method is applied to effectively reduce the amount of calculation; i.e., when is already obtained, is iteratively computed by where where is the fixed-step length of integrations. This approximation makes it possible to represent the PC expansion coefficients as a function of the discrete control variables . Moreover, the solution to the deterministic ODEs of the PC expansion coefficients is simplified through numerical integrations. The discrete-time dynamic optimization problem is as follows: Step 4. Incorporate the NIPC with PSO and solve the problem by the double-loop NIPC-PSO framework. Since the iteration optimization process of PSO depends on the performance of its fitness function, it is not directly available for the constrained optimization problem of equation (18). To solve this problem, a PSO application model is established here using the penalty function method, i.e., where is the integrated penalty function and denote the coefficients of different penalty components. These coefficients can be determined according to the importance of different penalty components using some evaluation methods, such as the analytic hierarchy process (AHP) [34]. The original equality constraint can be expressed by a pair of inequality constraints. It should be noted that significant differences between the magnitudes of the penalty components will result in calculation difficulties. For fairness, the normalization process should be applied to map different objectives into comparable scales in some cases. For PSO, its population is represented by a swarm of particles, i.e., , and the maximum number of iterations is set to be . Each particle represents a potential solution to the optimization problem and consists of a position vector and a velocity vector as follows: where and denote the position vector and velocity vector of particle , respectively; is the number of discrete points of decision parameters in equation (12); namely, . For the first iteration, a swarm of particles is generated randomly, i.e., where and are, respectively, the lower and upper bounds of the design parameters ; and are random coefficients between 0 and 1; and are the lower and upper bounds for the velocity vector. For the iteration, the performance of the swarm of particles is represented by their fitness functions , which are calculated, respectively, in the inner loop based on the NIPC method and the numerical integration scheme. When the calculation of the inner loop is all finished, the best position of each particle is recorded and the global best position of the whole swarm up to the current iteration is updated, i.e.,
Afterward, the velocity and position vector of each particle are updated in the outer loop as follows:where denotes the inertia parameter and can be linearly decreased with iteration for better optimization effect [35], i.e., ; and are two trust parameters, which represent confidence the particle has in itself and in the swarm, respectively; and are two independent random variables between 0 and 1. To keep the position vector of the next generation within the feasible region, and are chosen as follows:
Moreover, if any element in the updated velocity vector and position vector violates its limits, it is set as the closest boundary value before the next iteration starts.
To show the double-loop optimization algorithm more clearly, a flow chart of its implementation program is shown in Figure 1.
According to the above double-loop optimization framework, the PSO iterates for the global optimal solution in the outer loop, while the corresponding fitness functions with robustness considerations are calculated in the inner loop by the NIPC method and numerical integrations. Finally, the best particle of PSO which represents the optimal solution to the original robust optimization problem is obtained conveniently.
4. The Robust Trajectory Optimization Problem of HGV
In this section, a typical robust trajectory optimization problem of HGV is introduced, and the reference aircraft is chosen as the CAV-H which has been widely used in the research of hypersonic vehicles [9, 36].
Considering the Earth rotation model, the specific trajectory optimization formulation is as follows. Find the design variableswhich minimize the objective function,and subject to the dynamic constraints
The path constraints areand the boundary constraints arewhere is the angle of attack and is the bank angle; , , , and are the longitude, the latitude, the path angle, and the azimuth, respectively; , , , and are dimensionless variables of the geocentric distance, the velocity, the time, and the angular rate of Earth rotation, and their dimensionless parameters are , , , and , respectively; is the average radius of the Earth and is the gravitational acceleration at sea level; and represent the dimensionless lift and drag acceleration; , , and are the heat flux, the dynamic pressure, and the overload, respectively, in which is the vehicle head radius of curvature, is the atmospheric density, is the atmospheric density at sea level, and is a constant; , , , and are error bounds for corresponding terminal states.
In this paper, the path constraints are given as , , and , respectively. The initial conditions and desired terminal states are shown in Table 1. In the proposed method, the system state variables are obtained using numerical integrations based on given initial conditions; thus, only terminal conditions are required as boundary constraints. As described in equation (29), the corresponding error bounds of the terminal constraints are , , , and , respectively.
Aerodynamic forces play an important role in the unpowered reentry process of HGV. The dimensionless lift and drag acceleration shown in equation (27) is expressed as follows:where is the vehicle mass; is the pneumatic reference area; and are the lift and drag coefficients, which are functions of the angle of attack and Mach number.
Generally, there are no analytical expressions for aerodynamic coefficients and . Instead, they are calculated through numerical fitting based on experimental aerodynamic data. However, due to the complex atmospheric environment, these experimental data cannot be as accurate as anticipated. Therefore, it is assumed that the variation in and is a kind of uncertainty that can be represented by a single random variable [14, 26]. Otherwise, for the hypersonic vehicle, it is always not easy to achieve the designed reentry point accurately [36]. As a result, random deviations are also unavoidable at the beginning of the reentry process, i.e., the uncertain initial conditions. Based on the above discussions, two kinds of uncertainties are included in this paper, and the uncertain parameters are as follows:where indicates the aerodynamic uncertainty in the lift and drag coefficients, and indicates the state deviation uncertainty in the initial reentry conditions. It is assumed that both and are uniformly distributed, and they are independent of each other. In summary, the corresponding stochastic model is as follows:where and are the stochastic lift and drag coefficients, respectively; is the stochastic initial reentry states. It is assumed that the aerodynamic uncertainty lies in the interval of [−0.2, 0.2] and the state deviation uncertainty lies in the interval of [−0.01, 0.01]. This means that the deviation of the aerodynamic coefficients around their nominal values is about [−20%, +20%], while the deviation of the initial reentry conditions around their nominal values is about [−1%, +1%].
It should be noted that although the state vector includes six variables, we only introduce uncertainty to the altitude, the longitude, the latitude, and the speed due to their more demanding requirements. When equation (32) is incorporated, the deterministic trajectory optimization problem is now reformulated to a robust trajectory optimization problem.
The penalty components of the PSO application model in equation (19) are as follows:
A maximum normalization method is adopted to map magnitudes of different penalty components into comparable scales, i.e.,where , indicate the values of diffident components in equation (33) and is the number of all the terms; is the corresponding upper limit of and is its normalized value. After normalization, all the objective functions are mapped to the interval of [0, 100]. The objective function turns out to bewhere all the penalty coefficients are user-defined parameters determined through the AHP method.
5. Simulation and Results
In this section, three simulation cases are carried out to illustrate the proposed algorithm for the above HGV reentry problem, i.e., the deterministic trajectory optimization without consideration of uncertainties (denoted as DO), the traditional robust trajectory optimization using MC (denoted as TRO) [18], and the proposed robust trajectory optimization in this paper (denoted as NIPC-PSO). Herein, TRO is set to a similar structure of NIPC-PSO whereas its uncertainty propagation process is performed using 50 runs of MC simulations instead of the NIPC.
During optimization, the number of control discretization nodes is set as 100; the number of numerical integration points for the states is set as 500. Moreover, the swarm size of the PSO particles is 50; the maximum number of iterations is 100; is linearly decreased from 0.9 to 0.4, and ; the angle of attack is limited within and the bank angle is limited within . Since the twice-order PC expansion is accurate enough for common applications [37], herein, the NIPC expansion order is set as and the number of Gauss points is set as in order to balance the calculation burden and accuracy. The Gauss points and the corresponding weights are shown in Table 2.
All the simulations are performed with Matlab2014 on a personal computer with the Intel Core i5-8400 2.81 GHz processor and 16.0 G of RAM. After abundant numerical simulations, the convergence results of DO, NIPC-PSO, and TRO are achieved. The complexity analysis of the three cases is shown in Table 3.
Since the control parameterization method is adopted in this paper, only control variables are discretized and the number of design variables in the three cases is the same. As introduced in equation (27), HGV’s dynamic constraints are complex and time-consuming to solve. Thus, the complexity of the three cases is mainly reflected in their augmented ODEs. As shown in Table 3, DO only needs to deal with the original six-dimensional ODEs, while the number of ODEs in TRO and NIPC-PSO increases to 50 and 36, respectively. The number of numerical integrations as well as the corresponding running time of the three cases is directly proportional to their number of ODEs. It is indicated that both TRO and NIPC-PSO are more complex and time-consuming than DO, but NIPC-PSO is obviously simpler than that TRO of 50 MC simulations. Since the trajectory optimization is carried out in the offline phase of trajectory design, it is worth sacrificing some design time to obtain better robustness performance. Furthermore, under the double-loop optimization framework of this paper, the main cost of robust optimization is the repeated numerical integrations which is much easier to solve than adding design variables.
To compare the robustness of the designed trajectories in the three cases, the more detailed simulation verification and analysis are carried based on the optimal control laws of DO, TRO, and NIPC-PSO. Specifically, 10,000 times MC tests are implemented through numerical integrations under the same stochastic model of robust optimization. The obtained statistical trajectory data are summarized in Table 4, and the more intuitive comparisons are shown in Figures 2–8 distinctly.
To the statistical data of MC tests, the mean value of different parameters represents corresponding trajectory constraint satisfaction, while the standard deviation is usually regarded as a direct measure of robustness. As shown in Table 4, the mean values of path constraints and boundary constraints for NIPC-PSO all meet the design requirements well, whereas the maximum dynamic pressure of DO exceeds the upper limit of 6 and the terminal altitude of TRO has a large deviation from the reference value of 20,000 m. All the standard deviations of NIPC-PSO are significantly better than those of DO. Moreover, except for the maximum dynamic pressure, they are also better than the standard deviations of TRO. It should be noted that, in some indicators, the standard deviations of TRO are improved, but in others, they become worse. This is because there are only 50 MC simulations in TRO, and the uncertainty propagation process is not accurate enough. To improve this problem, more MC simulations are needed, but the complexity of TRO will be greatly increased.
In summary, the trajectories of DO are more likely to become infeasible with disturbances, and this is because no uncertainty of the online flight is considered by the traditional deterministic optimization. Although there is a little sacrifice in the original objective function of flight time for the inclusion of system robustness considerations, all the trajectory indicators of NIPC-PSO are less sensitive to uncertainties which is exactly what this article prefers. Furthermore, compared with TRO, NIPC-PSO can achieve better performance with much less computation.
Based on the MC tests of the three cases, the evolution of different states along the reentry process over time is shown in Figure 2. As can be seen, all the states fluctuate in a certain range around their nominal values and the states of NIPC-PSO have smaller fluctuations.
More clearly, Figure 3 shows the evolution standard deviations of the system states in the three cases. We can see that the standard deviations of NIPC-PSO are always smaller than those of DO and TRO, and it becomes more obvious when approaching the terminal time. This also suggests that the optimal control law of NIPC-PSO is better in tolerability and stability to the uncertainties. The conclusion is consistent with Table 4.
Figure 4 shows the normalized optimal penalty function values of MC tests in the three cases. As the figure reveals, both the stability and optimal values of the penalty functions in NIPC-PSO are the best in the three cases. According to Table 4, the standard deviations of these normalized penalty functions in DO, TRO, and NIPC-PSO are 42.61, 23.78, and 13.99, respectively. Compared with DO, these standard deviations of TRO and NIPC-PSO have been improved by 44% and 67%, respectively. This indicates that both TRO and NIPC-PSO have better robustness performance than DO, whereas the optimal penalty function value of TRO needs further improvement.
Since path constraints are significant to the trajectory performance, the mean values of the heat flux, the dynamic pressure, and the overload of MC tests in the three cases are shown in Figure 5. As we can see, although under uncertainties, path constraints in TRO and NIPC-PSO still show better capability than DO. This conclusion is in accordance with the performance of their penalty functions.
The terminal constraint distributions of altitude, longitude, and latitude are shown in Figures 6–8, respectively. Compared with DO and TRO, the accuracy of terminal state distributions of NIPC-PSO is improved obviously. The accuracy of TRO is not stable enough, and for better performance, it needs more simulations. Simulation results indicate that the proposed algorithm can achieve better robustness with less computation.
In the proposed algorithm, the NIPC method plays an important role in providing statistical information for the double-loop optimization framework. To verify the efficiency and accuracy of the NIPC method adopted in this paper, the statistical information obtained by the NIPC method is compared with that of MC tests. Specifically, based on the control law of NIPC-PSO and the same stochastic HGV model of Section 4, the uncertainty propagations are carried out using the NIPC method and 1000 MC tests, respectively. Statistical analysis is then conducted to obtain the mean and standard diversions of these six trajectory states. The comparison results of the NIPC method and MC tests are shown in Figures 9 and 10, respectively.
As shown in Figure 9, the mean values of the two methods are almost the same for all the states. In Figure 10, their corresponding standard deviations basically match each other well. As expected, the NIPC method is suitable and accurate enough for our robust trajectory planning problem and it has significantly fewer samples than the MC method. However, it should be noted that there are some errors between the two methods in Figure 10, especially in the altitude. This is because the PC expansions adopted in this paper are only twice-order and the altitude changes sharply during the skip reentry process. The approximation accuracy of NIPC can be improved by using the higher-order PC expansions, but the computational amount will be increased somewhat. Considering the trajectory performance achieved above, these errors are acceptable in this paper and the NIPC has good application potential in such complex robust optimization problems.
6. Conclusion
To improve the robustness of the reference trajectory for HGV with parametric uncertainties, a hybrid double-loop NIPC-PSO methodology is developed in this paper. The stochastic optimal control problem is transformed into a deterministic vision with higher-dimensional ODEs and robustness considerations using the NIPC method. The large number of calculations of the PC expansion coefficients is reduced through the control parameterization method and the numerical integration scheme. Based on an unconstrained model with a penalized fitness function, the PSO method searches globally for the optimal control law in the outer loop, while the penalized fitness function with robustness considerations is calculated repeatedly in the inner loop. The HGV reentry simulation of a classical time-optimal trajectory optimization problem with uncertainties in both initial states and aerodynamic coefficients shows that the NIPC-PSO of our proposed method has better robustness and performance in constraint satisfaction compared with the DO and TRO. The proposed procedure improves the feasibility and effectiveness of NIPC-based optimization methods.
Data Availability
The CAV data used to support the findings of this study are included within the article. The more detailed algorithm parameters are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there are no conflicts of interest.
Acknowledgments
This work was financially supported by the Natural Science Foundation of Shaanxi Province (No. 2020JQ-318) and the Fundamental Research Funds for the Central Universities (No. XJS190404 and JB190411). The authors gratefully acknowledge all the supports.