# Robust Trajectory Planning for Hypersonic Glide Vehicle with Parametric Uncertainties

Authors: Chunyun Dong; Zhi Guo; Xiaolong Chen
Journal: Mathematical Problems in Engineering
 (2021)

---

## 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.


---

## Body

## 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 Section2, 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 optimizationPD is as follows:(1)Problem:PDfindUtminJ=MUt0,Xt0,Utf,Xtf,t0,tf+∫t0tfLUt,Xt,tdtas.t.X˙t=FUt,Xt,tbCUt,Xt,t≤0cBXt0,t0,Xtf,tf=0,dwhere Ut 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; t0 and tf denote the initial and final time, respectively; Xt is the state variable vector and X˙t=FUt,Xt,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. 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 problemPR is formulated as follows:(2)Problem:PRfindUtminJA=kμ⋅JμUt,Xt,ξ,t+kσ⋅JσUt,Xt,ξ,tas.t.X˙t,ξ=FUt,Xt,ξ,t,ξbCμUt,Xt,ξ,t≤0BμXt0,ξ,t0,Xtf,ξ,tf=0cCσUt,Xt,ξ,t≤εCBσXt0,ξ,t0,Xtf,ξ,tf≤εB,dwhere ξ is the random variable vector; JA 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σ; Xt,ξ denotes the stochastic state variable vector; X˙t,ξ=FUt,Xt,ξ,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.


## 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 aq−dimensional random variable vector ξ, the exact solution Xt,ξ to the original stochastic ODEs of equation (2b) can be represented as the following NIPC model:(3)Xt,ξ≈X˜t,ξ=∑i=0pX˜itϕiξ,whereX˜t,ξ is the approximation of Xt,ξ; 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˜it are the PC expansion coefficients; ϕiξ are the orthogonal PC basis functions belonging to the Askey family [20] and satisfy the following orthogonal properties:(4)<ϕiξ,ϕjx>=δij<ϕiξ2>,whereδij is the Kronecker delta function, and ⋅,⋅ denotes the inner product, defined by(5)<fξ,gξ>=∫fξgξWξdξ,wherefξ and gξ are functions of ξ and Wξ 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 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 coefficientsX˜it of equation (3) can be calculated by the following integral [26]:(6)X˜it=1ϕiξ2∫Xt,ξϕiξWξdξ,i=0,1,2,…,p.Based on the above formulation, the full tensor-product numerical quadrature rule can be used to approximateX˜it, i.e.,(7)X˜it=∑k1=1m⋯∑kq=1mXt,ξk1,…,ξkqϕiξk1,…,ξkqϕi2∏j=1qωkj,wherem is the number of quadrature points along each dimension of ξ; ξkj,k=1,…,m,j=1,…,q, are the kth integration points of the jth component of ξ, and ωkj is the multidimensional quadrature weight corresponding to the jth dimension of the quadrature point ξkj,k=1,…,m.To implement the above numerical approximation of equation (7), ξk1,…,ξkq and ωkj need to be determined firstly according to a proper quadrature rule, such as Gaussian quadrature. Then, compute ϕiξk1,…,ξkq and ϕi2. Afterward, calculate Xt,ξk1,…,ξkq by substituting ξk1,…,ξkq into equation (2b). According to equation (7), there are ∏i=1qm=mq quadrature points of ξk1,…,ξkq, i.e., ξ1,ξ2,…,ξn,n=mq. Therefore, mq deterministic ODEs are generated and need to be solved, i.e.,(8)X˙t,ξ1=FUt,Xt,ξ1,t,ξ1X˙t,ξ2=FUt,Xt,ξ2,t,ξ2⋮X˙t,ξn=FUt,Xt,ξn,t,ξn,n=mq,whereX˙t,ξn=FUt,Xt,ξn,t,ξn is the augmented deterministic ODEs with ξn.Step 2. Establish the equivalent deterministic trajectory optimization model PED of PR based on NIPC.According to the NIPC model of equation (3), the original stochastic state Xt,ξ can be represented as X˜t,ξ. When the PC expansion coefficients X˜it,i=0,1,…,p have been calculated, its first two statistical moments can be easily computed as follows:(9)Xμt,ξ=X˜0t,Xσt,ξ=∑i=1pX˜it2<ϕi2>,whereXμt,ξ and Xσt,ξ are the mean and the standard deviations of Xt,ξ.The original performance objectiveJ and the constraints C and B are all functions with respect to the approximate stochastic state X˜t,ξ. 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 problemPR is approximately transcribed into a deterministic optimal control problem PED, i.e.,(10)Problem:PEDfindUtminJA=kμ⋅JμUt,X˜t,ξ,t+kσ⋅JσUt,X˜t,ξ,tas.t.X˙t,ξ1=FUt,Xt,ξ1,t,ξ1X˙t,ξ2=FUt,Xt,ξ2,t,ξ2,n=mq⋮X˙t,ξn=FUt,Xt,ξn,t,ξnbCμUt,X˜t,ξ,t≤0BμX˜t0,ξ,t0,X˜tf,ξ,tf=0cCσUt,X˜t,ξ,t≤εC,BσX˜t0,ξ,t0,X˜tf,ξ,tf≤εB.dThe 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 PEDZ and use the numerical integration scheme to simplify the calculation of PC expansion coefficients.To solve the above continuous optimization problemPED 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 PED, there are mq 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 fromt0 to tf is first divided into N subintervals as follows:(11)Z=t0,t1,…,tN,t0<t1<⋯<tN−1<tN=tf,whereZ is the discretized sampling time.Then, the control variables at the discrete nodes are selected as the decision variables to be optimized, that is,(12)UZ=Ut0T,Ut1T,…,UtNTT,whereUZ is the discretized control vector.The discrete-time dynamic optimization problem is established by fewer design variables; however, themq 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 vectorUZ, during each segment Δt∈ti,ti+1, it can be approximated by linear interpolation, i.e.,(13)UΔti≈Uti+Uti+1−Utiti+1−tiΔti−ti,Δti∈ti,ti+1,i=0,1,…,N−1.The approximate continuous control variableU˜t can be represented as follows:(14)U˜t=∑i=0N−1UΔti.For the convenience of description, represent the expanded ODEs of equation (8) in a brief form as follows:(15)X˙Dt=FDU˜t,XDt,t,whereXD denotes the deterministic states Xt,ξ1,Xt,ξ2,…,Xt,ξn and FD denotes the augmented deterministic ODEs. According to the numerical integration theory, the state variables propagate from t0 to tf 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 XDtk is already obtained, XDtk+1=XDtk+h is iteratively computed by(16)XDtk+1=XDtk+h6k1+2k2+2k3+k4,where(17)k1=FDU˜tk,XDtk,tk,k2=FDU˜tk+h2,XDtk+h2k1,tk+h2,k3=FDU˜tk+h2,XDtk+h2k2,tk+h2k4=FDU˜tk+h,XDtk+hk3,tk+h,whereh=tk+1−tk 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 UZ. Moreover, the solution to the mq deterministic ODEs of the PC expansion coefficients is simplified through numerical integrations. The discrete-time dynamic optimization problem PEDZ is as follows:(18)Problem:PEDZfindUZminJ=kμ⋅JμU˜t,X˜t,ξ,t+kσ⋅JσU˜t,X˜t,ξ,tas.t.X˙Dt=FDU˜t,XDt,tbCμU˜t,X˜t,ξ,t≤0BμX˜t0,ξ,t0,X˜tf,ξ,tf=0cCσU˜t,X˜t,ξ,t≤εC,BσX˜t0,ξ,t0,X˜tf,ξ,tf≤εB.dStep 4. Incorporate the NIPC with PSO and solve the problem PEDZ 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.,(19)Problem:PEDZ−PSOfindUZminP=pJA⋅JA+pCμ⋅Cμ+pBμ⋅Bμ+pCσ⋅Cσ+pBσ⋅Bσas.t.X˙Dt=FDU˜t,XDt,t,bwhereP is the integrated penalty function and p∗ 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 Bμ 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 ofNS particles, i.e., Ωp=uk,vk,k=1,2,…,NS, and the maximum number of iterations is set to be NI. Each particle represents a potential solution to the optimization problem PEDZ−PSO and consists of a position vector and a velocity vector as follows:(20)uk=uk,1,uk,2,…,uk,nZT,vk=vk,1,vk,2,…,vk,nZT,whereuk and vk denote the position vector and velocity vector of particle k, respectively; nZ is the number of discrete points of decision parameters in equation (12); namely, nZ=dimUZ.For the first iteration, a swarm ofNS particles Ωp0=uk0,vk0,k=1,2,…,NS is generated randomly, i.e.,(21)uk0=uL+ru⋅uU−uLvk0=vL+rv⋅vU−vL,k=1,2,⋯,NS,whereuL and uU are, respectively, the lower and upper bounds of the design parameters UZ; ru and rv are random coefficients between 0 and 1; vL and vU are the lower and upper bounds for the velocity vector.For thes−th iteration, the performance of the swarm of NS particles is represented by their fitness functions Pks,k=1,2,…,NS, 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 pks of each particle k is recorded and the global best position pgs of the whole swarm up to the current iteration is updated, i.e.,(22)pks=uklm,lm=argminp1=1,2,…,sPkp1,k=1,2,…,NS,pgs=uqmzm,zm,qm=argminz1=1,2,…,sq1=1,2,…,NSPq1z1.Afterward, the velocity and position vector of each particle are updated in the outer loop as follows:(23)vks+1=ω⋅vks+c1⋅r1⋅pks−uks+c2⋅r2⋅pgs−uksuks+1=uks+vks+1,k=1,2,…,NS,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/NI,s=1,2,…,NI; c1 and c2 are two trust parameters, which represent confidence the particle has in itself and in the swarm, respectively; r1 and r2 are two independent random variables between 0 and 1. To keep the position vector of the next generation within the feasible region, vL and vU are chosen as follows:(24)vL=uL−uU,vU=uU−uL.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 Figure1.Figure 1
Flow chart of the algorithm program.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 variables(25)U=αcσcT,which minimize the objective function,(26)J=tf,and subject to the dynamic constraints(27)dydt=Vsinγ,dθdt=Vcosγsinψycosφ,dφdt=Vcosγcosψy,dVdt=−D¯−sinγy2+Ω2ycosφsinγcosφ−cosγsinφcosψ,dγdt=1VL¯cosσc+cosγyV2−1y+2ΩVcosφsinψ+Ω2ycosφcosγcosφ+sinγsinφcosψ,dψdt=1VL¯sinσccosγ+V2ycosγsinψtanφ−2ΩVtanγcosφcosψ−sinφ+Ω2ycosγsinφcosφsinψ.The path constraints are(28)C=Q˙=C1Rρρs0.5V3.15≤Q˙max,qr=12ρV⋅Vc2≤qmax,nr=L¯2+D¯2≤nmax,and the boundary constraints are(29)B=ytf−yf≤εfy,θtf−θf≤εfθ,φtf−φf≤εfϕ,Vtf−Vf≤εfV,where αc is the angle of attack and σc is the bank angle; θ, φ, γ, 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 R0, Vc=g0R00.5, R0/g00.5, and g0/R00.5, respectively; R0 is the average radius of the Earth and g0 is the gravitational acceleration at sea level; L¯ and D¯ represent the dimensionless lift and drag acceleration; Q˙, qr, and nr 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 C1 is a constant; εfy, εfθ, εfϕ, and εfV are error bounds for corresponding terminal states.In this paper, the path constraints are given asQ˙max=1000kW/m2, qmax=400kPa, and nmax=6, 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 εfy=0, εfθ=0, εfϕ=0, and εfV=0, respectively.Table 1
Initial conditions and terminal constraints of the states.
Constraintsykmθ°φ°Vm/sγ°ψ°Initial conditions100160257200-258Terminal constraints20232371000Aerodynamic 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:(30)L¯=ρ⋅V⋅Vc2⋅SV⋅CL2mV⋅g0,D¯=ρ⋅V⋅Vc2⋅SV⋅CD2mV⋅g0,where mV=907.2kg is the vehicle mass; SV=0.4839m2 is the pneumatic reference area; CL and CD 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 coefficientsCL and CD. 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 CL and CD 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:(31)ξ=ξC,ξIT,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:(32)C˜L=CL1+ξC,C˜D=CD1+ξC,X˜0=X01+ξI,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]. 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:(33)J=tf,Q˙,qr,nr,ytf−yf,θtf−θf,φtf−φf,Vtf−VfT.A maximum normalization method is adopted to map magnitudes of different penalty components into comparable scales, i.e.,(34)J¯j=JjJjU×100,j=1,2,…,h0,where Jj=Jj,j=1,2,…,h0, indicate the values of diffident components in equation (33) and h0 is the number of all the terms; JjU is the corresponding upper limit of Jj and J¯j is its normalized value. After normalization, all the objective functions are mapped to the interval of [0, 100]. The objective function turns out to be(35)minP=kμ∑j=1h0μJ¯jkμj+kσ∑j=1h0σJ¯jkσj,where 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 c1=c2=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. The Gauss points and the corresponding weights are shown in Table 2.Table 2
Gauss Legendre quadrature points and weights.
PointsWeights−0.932470.17132−0.661210.36076−0.238620.467910.238620.467910.661210.360760.932470.17132All 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 Table3.Table 3
Complexity analysis of the three cases.
MethodsParticlesIterationsNumber of ODEsNumber of integrationsRunning time (s)DO5010015000213.9TRO501005025000010550.2NIPC-PSO50100361800007670.4Since 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 Table4, and the more intuitive comparisons are shown in Figures 2–8 distinctly.Table 4
Statistical results of the three cases.
Mean valuesStandard deviationsDOTRONIPC-PSODOTRONIPC-PSOP80.5792.1941.1542.6123.7813.99Q˙maxkW/m2854.13776.84815.9727.3328.5525.26qmaxkPa77.3647.6841.1014.8648.625.12nmax6.083.803.100.680.440.57ytfm1972622833200882405.611803.71780.09θtf°231.75231.05231.533.032.611.83φtf°37.0137.1536.890.270.290.18Figure 2
Boundaries of the trajectory states in MC tests of the three cases.Figure 3
Comparison of the trajectory state standard deviations in MC tests of the three cases.Figure 4
Comparison of the normalized penalty functions in MC tests of the three cases.Figure 5
Mean values of path constraints in MC tests of the three cases.Figure 6
Distributions of terminal altitude in MC tests of the three cases.Figure 7
Distributions of terminal longitude in MC tests of the three cases.Figure 8
Distributions of terminal latitude in MC tests of the three cases.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 Table4, 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 Figure2. 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, Figure3 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.Figure4 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 Figure5. 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 Figures6–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 Section4, 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.Figure 9
Mean values of states of MC tests and NIPC method in NIPC-PSO case.Figure 10
Standard deviations of states of MC tests and NIPC method in NIPC-PSO case.As shown in Figure9, 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.



---

*Source: 3676810-2021-01-16.xml*