Robust Trajectory Planning for Hypersonic Glide Vehicle with Parametric Uncertainties

. 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 ﬁtness 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 satisﬁed 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 coeﬃcients. The results demonstrate the feasibility and eﬀectiveness of the proposed method.


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][2][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].erefore, the reentry trajectory design, as a core problem of the guidance and control, has become a challenge and hotspot for HGV [5,6].e performance of trajectory design is the essential section of a flight.e 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.
e 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][9][10][11][12].According to Betts' survey on numerical trajectory optimization methods, these algorithms can be divided into two categories [13].e 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].
ese methods can obtain accurate analytical solutions but their initial states are difficult to guess.e 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.
ese 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.e 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].e 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.
e 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.erefore, 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.ese 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.
e 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.
e 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 2 Mathematical Problems in Engineering 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.en, 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.erefore, 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.
e 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.
e number of decision variables as well as the huge amount of computation for uncertainty propagation using NIPC is extremely reduced.
ese above properties of the proposed approach indicate its high efficiency and potential applicability for robust trajectory optimization problems of HGV. e remainder of this paper is organized as follows.In Section 2, the robust trajectory optimization model is established.en, 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.e numerical simulation examples are included in Section 5, and conclusions are drawn in Section 6.

Problem Formulation
e trajectory optimization is a typical nonlinear constrained optimal control problem.A commonly used optimal control formulation of the deterministic trajectory optimization P D is as follows: where U(t) is the control variable vector and t is the time; J is the objective function, in which M is the cost of boundary constraints, and L is the cost of Lagrange (integral) term; t 0 and t f denote the initial and final time, respectively; X(t) is the state variable vector and _ X(t) � F(U(t), X(t), t) is the system ordinary differential equations, i.e., ODEs; C and B 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 t [23,25], whereas the performance objective J, the ODEs, and the constraints C and B are all stochastic functions of the random variables.e goal of robust trajectory optimization is to obtain the optimal control law, which is insensitive to possible disturbances.erefore, 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 P R is formulated as follows: Mathematical Problems in Engineering where ξ is the random variable vector; J A is the augmented objective function; J μ and J σ are the expectation and standard deviation of J, respectively; k μ and k σ are the weight factors of J μ and J σ ; X(t, ξ) denotes the stochastic state variable vector; _ X(t, ξ) � F(U(t), X(t, ξ), t, ξ) denotes the stochastic ODEs; the variables C μ , B μ and C σ , B σ are the expectations and standard deviations of C and B, respectively; ε C and ε B are allowable standard deviation thresholds of C and B.
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.

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 q − dimensional random variable vector ξ, the exact solution X(t, ξ) to the original stochastic ODEs of equation (2b) can be represented as the following NIPC model: where  X(t, ξ) is the approximation of X(t, ξ); p is the number of expansion terms that is determined by p + 1 � (((q + r)!)/(q!r!)) and r is the desired approximation order of orthogonal polynomials;  X i (t) are the PC expansion coefficients; ϕ i (ξ) are the orthogonal PC basis functions belonging to the Askey family [20] and satisfy the following orthogonal properties: where δ ij is the Kronecker delta function, and 〈•, •〉 denotes the inner product, defined by where f(ξ) and g(ξ) are functions of ξ and W(ξ) is a weight function.ere are a corresponding polynomial basis function and weight function for a random variable with a given distribution.Moreover, for the q − dimensional random vector ξ, ϕ i (ξ) and W(ξ) 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  X i (t) 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 where m is the number of quadrature points along each dimension of ξ; ξ k j , k � 1, . . ., m, j � 1, . . ., q, are the kth integration points of the jth component of ξ, and ω k j is the multidimensional quadrature weight 4 Mathematical Problems in Engineering corresponding to the jth dimension of the quadrature point ξ k j , k � 1, . . ., m.To implement the above numerical approximation of equation ( 7), (ξ k 1 , . . ., ξ k q ) and ω k j need to be determined firstly according to a proper quadrature rule, such as Gaussian quadrature.en, compute ϕ i (ξ k 1 , . . ., ξ k q ) and 〈ϕ 2 i 〉.Afterward, calculate X(t, ξ k 1 , . . ., ξ k q ) by substituting (ξ k 1 , . . ., ξ k q ) into equation (2b).According to equation (7), there are  q i�1 m � m q quadrature points of (ξ k 1 , . . ., ξ k q ), i.e., ξ 1 , ξ 2 , . . ., ξ n  , n � m q .erefore, m q deterministic ODEs are generated and need to be solved, i.e., where _ Step 2. Establish the equivalent deterministic trajectory optimization model P E D of P R based on NIPC.
According to the NIPC model of equation ( 3), the original stochastic state X(t, ξ) can be represented as  X(t, ξ).When the PC expansion coefficients  X i (t), i � 0, 1, . . ., p have been calculated, its first two statistical moments can be easily computed as follows: where X μ (t, ξ) and X σ (t, ξ) are the mean and the standard deviations of X(t, ξ). e original performance objective J and the constraints C and B are all functions with respect to the approximate stochastic state  X(t, ξ).erefore, 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 P R is approximately transcribed into a deterministic optimal control problem P E D , i.e., Problem : e 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 P E DZ and use the numerical integration scheme to simplify the calculation of PC expansion coefficients.To solve the above continuous optimization problem P E D 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 P E D , there are m q deterministic ODEs and the number of state variables is quite a lot.us, 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.
e time interval from t 0 to t f is first divided into N subintervals as follows: where Z is the discretized sampling time.
Mathematical Problems in Engineering en, the control variables at the discrete nodes are selected as the decision variables to be optimized, that is, where U(Z) is the discretized control vector.e discrete-time dynamic optimization problem is established by fewer design variables; however, the m q 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 U(Z), during each segment Δt ∈ [t i , t i+1 ], it can be approximated by linear interpolation, i.e., e approximate continuous control variable  U(t) can be represented as follows: For the convenience of description, represent the expanded ODEs of equation ( 8) in a brief form as follows: where X D denotes the deterministic states X(t, ξ 1 ), X(t, ξ 2 ), . . ., X(t, ξ n )   and F D denotes the augmented deterministic ODEs.According to the numerical integration theory, the state variables propagate from t 0 to t f and their values of any time can be obtained.Herein, by providing control inputs  U(t) 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 Determine the discrete control vector U(Z), the swarm size N S , the iterations times N I , and other operating parameters of PSO to reform the problem P ED to P EDZ .n > m q ?Converge?Or j > N I ?
Compare all the results of N S particles and obtain the optimal control law of the j-th iteration.
Calculate the required statistical information using NIPC and obtain the fitness function Obtain the optimal control law for the robust trajectory planning problem.
Update all the particles based on the last iteration and set j = j + 1, i = 1.
For a q-dimensional random variable vector ξ of the problem P R , using NIPC with m quadrature points to transform the stochastic ODEs into deterministic and obtain the problem P ED .
For the i-th particle U i (Z) of the j-th iteration, solve the ODEs with ξ n through numerical interpolation and integration algorithms, i.e.
Initialize all the particles for the first iteration of PSO and set parameters j = 1, i = 1, and n = 1.Mathematical Problems in Engineering

Yes
where where h � t k+1 − t k is the fixed-step length of integrations.is approximation makes it possible to represent the PC expansion coefficients as a function of the discrete control variables U(Z).Moreover, the solution to the m q deterministic ODEs of the PC expansion coefficients is simplified through numerical integrations.e discrete-time dynamic optimization problem P E DZ is as follows: 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., Problem: P E DZ− PSO find U(Z) where P is the integrated penalty function and p * denote the coefficients of different penalty components.ese coefficients can be determined according to the importance of different penalty components using some evaluation methods, such as the analytic hierarchy process (AHP) [34].e original equality constraint B μ can be expressed by a pair of inequality constraints.It should be noted that significant Mathematical Problems in Engineering 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 N S particles, i.e., Ω p � (u k , v k )  , k � 1, 2, . . ., N S , and the maximum number of iterations is set to be N I .Each particle represents a potential solution to the optimization problem P E DZ− PSO and consists of a position vector and a velocity vector as follows: where u k and v k denote the position vector and velocity vector of particle k, respectively; n Z is the number of discrete points of decision parameters in equation ( 12); namely, n Z � dim(U(Z)).
For the first iteration, a swarm of is generated randomly, i.e., where u L and u U are, respectively, the lower and upper bounds of the design parameters U(Z); r u and r v are random coefficients between 0 and 1; v L and v U are the lower and upper bounds for the velocity vector.For thes − th iteration, the performance of the swarm of N S particles is represented by their fitness functions P s k , k � 1, 2, . . ., N S , 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 p s k of each particle k is recorded and the global best position p s g of the whole swarm up to the current iteration is updated, i.e., z m , q m  � arg min z 1 �1,2,...,s  Mathematical Problems in Engineering 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., ω s+1 � ω max − s((ω max − ω min )/N I ), s� 1, 2, . . ., N I ; c 1 and c 2 are two trust parameters, which represent confidence the particle has in itself and in the swarm, respectively; r 1 and r 2 are two independent random variables between 0 and 1.To keep the position vector of the next generation within the feasible region, v L and v U 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

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 variables which minimize the objective function, and subject to the dynamic constraints e path constraints are and the boundary constraints are where α c is the angle of attack and σ c is the bank angle; θ, φ, c, and ψ are the longitude, the latitude, the path angle, and the azimuth, respectively; y, V, t, and Ω are dimensionless variables of the geocentric distance, the velocity, the time, and the angular rate of Earth rotation, and their dimensionless parameters are R 0 , V c � (g 0 R 0 ) 0.5 , (R 0 /g 0 ) 0.5 , and (g 0 /R 0 ) 0.5 , respectively; R 0 is the average radius of the Earth and g 0 is the gravitational acceleration at sea level; L and D represent the dimensionless lift and drag acceleration; _ Q, q r , and n r are the heat flux, the dynamic pressure, and the overload, respectively, in which R is the vehicle head radius of curvature, ρ is the atmospheric density, ρ s is the atmospheric density at sea level, and C 1 is a constant; ε f y , ε f θ , ε f ϕ , and ε f V are error bounds for corresponding terminal states.
In this paper, the path constraints are given as _ Q max � 1000 kW/m 2 , q max � 400 kPa, and n max � 6, respectively.
e 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 ε f y � 0, ε f θ � 0, ε f ϕ � 0, and ε f V � 0, respectively.

12
Mathematical Problems in Engineering Aerodynamic forces play an important role in the unpowered reentry process of HGV. e dimensionless lift and drag acceleration shown in equation ( 27) is expressed as follows: where m V � 907.2 kg is the vehicle mass; S V � 0.4839 m 2 is the pneumatic reference area; C L and C D 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 C L and C D .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.erefore, it is assumed that the variation in C L and C D 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 ξ C indicates the aerodynamic uncertainty in the lift and drag coefficients, and ξ I indicates the state deviation uncertainty in the initial reentry conditions.It is assumed that both ξ C and ξ I are uniformly distributed, and they are independent of each other.In summary, the corresponding stochastic model is as follows: where  C L and  C D are the stochastic lift and drag coefficients, respectively;  X 0 is the stochastic initial reentry states.It is assumed that the aerodynamic uncertainty ξ C lies in the interval of [−0.2, 0.2] and the state deviation uncertainty ξ I lies in the interval of [−0.01, 0.01].
is means that the 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.
e 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 J j � J(j), j � 1, 2, . . ., h 0 , indicate the values of diffident components in equation ( 33) and h 0 is the number of all the terms; J jU is the corresponding upper limit of J j and J j is its normalized value.After normalization, all the objective functions are mapped to the interval of [0, 100].e objective function turns out to be where all the penalty coefficients are user-defined parameters determined through the AHP method.

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 c 1 � c 2 � 2; the angle of attack is limited within [10 °, 20 °] and the bank angle is limited within [−80 °, 80 °].Since the twice-order PC expansion is accurate enough for common applications [37], herein, the NIPC expansion order is set as r � 2 and the number of Gauss points is set as m � 6 in order to balance the calculation burden and accuracy.e 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.e 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.us, 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.
e 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 timeconsuming 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  Mathematical Problems in Engineering 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.e 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.is 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.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. is also suggests that the optimal control law of NIPC-PSO is better in tolerability and stability to the uncertainties.e 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.is 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. is conclusion is in accordance with the performance of their penalty functions.
e 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.e 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.

Mathematical Problems in Engineering
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.e 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.is is because the PC expansions adopted in this paper are only twice-order and the altitude changes sharply during the skip reentry process.e 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.

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.e stochastic optimal control problem is transformed into a deterministic vision with higher-dimensional ODEs and robustness considerations using the NIPC method.e 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.e 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.
e proposed procedure improves the feasibility and effectiveness of NIPC-based optimization methods.

Data Availability
e CAV data used to support the findings of this study are included within the article.
e more detailed algorithm parameters are available from the corresponding author upon request.

Figure 1 :
Figure 1: Flow chart of the algorithm program.

Figure 2 :
Figure 2: Boundaries of the trajectory states in MC tests of the three cases.

Figure 3 :
Figure 3: Comparison of the trajectory state standard deviations in MC tests of the three cases.

Figure 4 :
Figure 4: Comparison of the normalized penalty functions in MC tests of the three cases.

Figure 5 :
Figure 5: Mean values of path constraints in MC tests of the three cases.

Figure 6 :
Figure 6: Distributions of terminal altitude in MC tests of the three cases.

Figure 7 :
Figure 7: Distributions of terminal longitude in MC tests of the three cases.

Figure 8 :
Figure 8: Distributions of terminal latitude in MC tests of the three cases.

Figure 9 :
Figure 9: Mean values of states of MC tests and NIPC method in NIPC-PSO case.

Figure 10 :
Figure 10: Standard deviations of states of MC tests and NIPC method in NIPC-PSO case.

Table 1 :
Initial conditions and terminal constraints of the states.

Table 2 :
Gauss Legendre quadrature points and weights.

Table 3 :
Complexity analysis of the three cases.
Mathematical Problems in EngineeringStep 4. Incorporate the NIPC with PSO and solve the problem P E DZ by the double-loop NIPC-PSO framework.

Table 4 :
Statistical results of the three cases.