Air-breathing hypersonic vehicle trajectory optimization with uncertain no-fly zones

The air-breathing hypersonic planes are considered to be the future of commercial airlines, which can fly from NYC to London in under an hour. For air-breathing vehicles, 3D trajectory planning will become extremely important due to its significant impact on flight performance. Past research on this issue was not comprehensive enough. Thus, we proposed a hybrid method for generating the optimal trajectories efficiently. No-fly zones are specified for geopolitical restrictions in air-breathing hypersonic vehicle missions. However, previous studies focused on no-fly zone constraints with fixed locations and boundaries. For robust execution, we must take into account no-fly zones’ uncertainties, which arise due to uncertain localization, measurement errors, and no-fly zones’ movement. A chance-constrained approach is presented here to deal with such uncertainties. The key idea of this method is controlling risk when the flight path is close to the uncertain no-fly zones. First, the trajectory planning of air-breathing hypersonic vehicles is modeled as a chance-constrained optimization (CCOPT) problem. With the help of the convexification and linearization techniques, we can approximate the non 3.2.2 CCOPT problem as a convex second-order cone programing under the Gaussian distribution assumption. It can be solved to global optimality by the interior point method. Finally, we give numerical simulation results to prove the effectiveness of this method.


Introduction
Trajectory planning for hypersonic vehicles such as an air-breathing hypersonic vehicle (ABHV) has received a great deal of attention in recent years. [1][2][3][4] An ABHV needs to be able to plan trajectories that take the vehicle from its current location to a goal while avoiding no-fly zones (NFZs). These trajectories should be optimal with respect to a criterion such as time or fuel consumption. This problem has two main challenges: First, the trajectory optimization problem of ABHV is highly complicated and inherently nonconvex. Second, there are a number of sources of uncertainty of no-fly zones, such as uncertain localization, measurement errors, and NFZs' movement. Such uncertainty must be modeled as a form that can be processed by numerical algorithms to achieve robustness.
Many methods can be used to solve the problem of unpowered hypersonic gliding vehicle (HGV) trajectory planning. Model predictive programing, 5 pseudospectral method, 1,4,6 indirect method, 7 and convex optimization 8 all demonstrate their unique advantages in certain scenarios. Considering the need to generate trajectories quickly in some tasks, we proposed an algorithm based on geometric stitching. 9 As an application of heuristics algorithm, particle swarm optimization 10 and gravitational search algorithm 11 can also give appropriate solutions. Furthermore, with the development of machine learning, deep neural networks have been used in HGV's trajectory planning. 12,13 However, ABHV's dynamic is much more complex. Its scramjet engine's performance is highly correlated with the vehicle's altitude, velocity, angle of attack, and throttle (as shown in equation (8)). As a result, the optimization algorithm is required to solve a multi-variable, highly non-linear and complex coupling problem in a limited time. Those methods developed for the unpowered vehicle are unsuitable and computationally expensive in this situation. Some existing work circumvents this problem using relatively simple propulsion system models. 14,15 In this work, we build a hybrid algorithm for ABHV with complex propulsion-pneumatic coupling based on convex optimization, 8 which has shown great application potential in recent years because of its serval superior features: 16 Every local minimum is a global minimum; if the objective function is strictly convex, then the problem has at most one optimal point. Some kinds of convex problems, like the second-order cone programing (SOCP), can be solved efficiently by primal-dual interior-point methods (IPMs) that only have polynomial complexity.
Previous work focuses on the no-fly zone with precisely known locations and boundaries, 8 which is modeled as a cylinder with its center at the specified longitude and latitude and a given radius. However, they do not explicitly take into account uncertainty, which arises for two reasons. First, the location and boundaries are not usually known exactly but are estimated using a statistical model, satellite imagery, or inaccurate intelligence analysis. Second, some NFZs can move even during the flight of ABHV. So in practice, the vehicle may collide with the real NFZ even if the planned trajectory did not. In this paper, we will address this issue.
There are two solutions. A feasible idea is to ensure that the aircraft will not enter NFZ even in the worst case, which leads to a conservative solution. Let NFZ R be the real NFZ, and NFZ P be the areas that NFZ R may exist. The conservative solution asks the ABHV to bypass the NFZ P even if NFZ P is much larger than NFZ R . It is a bad idea because of three reasons. First, it is a huge waste of aircraft performance. Second, ABHV's maneuverability and flight range are limited. Bypassing such a large NFZ P may cause mission failure. Not to mention the possibility that a NFZ P happens to block the way to target. Flying through the NFZ P is inevitable in this case. Third, ABHV's performance is highly influenced by the flight conditions, and bypassing will make the vehicle deviate from the designed optimistic conditions.
Another solution is to tolerate some degree of risk if the probability of failure is below a specified threshold. Therefore, it is reasonable to demand constraint satisfaction with some level of probability of success. This calls for the use of chance constraints. Chance constraints ensure that the decision made by a model leads to a certain probability of complying with constraints. In the chance-constrained optimization (CCOPT), there exist inequality constraints that are formulated as x is the state and Prob 2 ½0, 1 is the probability level. The chance constraints were first introduced by Charnes et al. 17 in 1958 to solve optimization problems under uncertainty. Since then, many studies from Pre´kopa, 18 Charnes and Cooper, 19 Li et al. 20 have looked into more methods of analyzing CCOPT and using CCOPT to solve engineering problems. Blackmore et al. 21 and Ono et al. 22 brought this concept into the field of dynamic trajectory planning and proved its value in a variety of missions, like UAV and Mars landing, etc. Their work inspired us to use a similar technique in uncertain NFZ problems.
This paper focuses on solving ABHV's trajectory planning problem with uncertain NFZ constraints. The rest of this article is organized as follows. Section 2 proposed the major problem of ABHV's trajectory planning. Section 3 gives a unique hybrid method for ABHV's trajectory planning that balances efficiency and accuracy. The convexification of chanceconstrained NFZ is presented in Section 4. Section 5 presents the numerical simulation results of our methods. This paper concludes with Section 6 by providing a summary, some comments, and pointing some future research directions.

A highly coupled multi constraints aircraft model for integrated flight/propulsion
In this section, we describe the generic trajectory optimization problem for the cruising phase of the airbreathing hypersonic vehicle. Due to the flight/ propulsion are highly coupled, the original problem must be formulated in a certain way conducive to the SOCP approach.
Regardless of the Earth rotation, a 3-DOF pointmass translation equations of motion of the ABHV with respect to non-dimensional time t can be described as _ g = L cos s + T sin a cos s + F Ty cos a cos s + (V 2 À 1=r) cos g=r _ c = (L sin s + T sin a sin s + F Ty cos a sin s)= cos g + V 2 cos g sin c tan f=r where r is the radial distance from the Earth center to the vehicle; u and f are the longitude and latitude, respectively; g is the flight-path angle of the Earth-relative velocity vector; c is the heading angle, which is measured clockwise in the local horizontal plane from the north. V donates the velocity of ABHV. All these variables are non-dimensional: the velocity is normalized by ffiffiffiffiffiffiffiffiffi ffi g 0 R 0 p , where g 0 = 9:81m=s 2 and the Earth radius is R 0 = 6378:135km. The length is scaled by R 0 and time t is normalized by ffiffiffiffiffiffiffiffiffiffiffiffi R 0 =g 0 p . F Ty is the additional normal force caused by the scramjet engine. 23 Let L be the non-dimensional lift acceleration in g 0 , D be the non-dimensional drag acceleration, and T be the non-dimensional thrust acceleration. They are given by where r is the local atmosphere density. C L , C D , and C T correspond to the lift, drag, and thrust coefficient, respectively. They are functions of Mach number, angle of attack a, non-dimensional throttle h, and altitude r. S ref and m are the reference area and mass of the vehicle, respectively. S Engine is the reference area of the scramjet engine. C f (Mach, a, r) is the fuel consumption rate.

State and control constraints
The control inputs s and a respectively donates the bank angle and angle of attack. Their saturation limitation can be presented by following equations a min ł a ł a max ð9Þ The throttle is also limited We use u 0 = ½a coss sin s h to replace the ½a s h in equations (1)- (8). It is because Liu et al., 27 and Zhao and Song 24 pointed out that directly using s as a control variable introduces significant high-frequency chatters, which is generally detrimental to the convergence of the successive solution process. Obviously, the new control input u must satisfy equation (12).
For an HGV, the total potential energy is monotonically decreasing during the gliding process because of drag, so an energy-like variable is used as the independent variable of the programing. 27 However, when it comes to an ABHV with scramjet, it is not suitable. As a result, we need to find a new independent variable. In all states in the equations (1)-(7), time t and mass m are monotonical during the flight. In this paper, we prefer mass to time as the independent variable due to the following reason: 1. The dynamic equations (1)-(7) can be reduced to 6 dimensions from 7 dimensions. 2. If we choose time as the independent variable and the terminal time is unavailable, it is hard to define the terminal constraints of the problem. In contrast, terminal mass is easy to be given in a maximum-downrange problem.
m is decreasing during the flight, so let the increasing variablem = À m be the new independent variable.
Substitute equation (13) into equations (1)-(7), we obtain dc=dm = (L sin s + T sin a sin s + F Ty cos a sin s)= cos g + V 2 cos g sin c tan f=r Thus the inequality constraints on controls that derived from equations (9)- (12), are given as The initial and terminal constraints are described respectively by x wherem 0 andm f correspond to the initial and terminal m of the vehicle, respectively.

Path and NFZ chance constraints
The heating rate _ Q, dynamic pressure q, and normal load n are also constrained: Equations (23)- (25) can also be expressed as altitude constraints because all of them are functions of r and V .
which can also be written as corresponding to differentm. Meanwhile, to avoid flying too high, the upper boundary is given by Thus, Besides, the velocity of ABHV is also constrained Equations (26), (29), and (30) yields the normalized path constraints: An NFZ is an area over which aircraft are not permitted to fly. NFZs may be established for safety, environmental protection, or other reasons. In previous works, the determined NFZs are modeled as no-fly cylinders with given radii and center coordinates: where u NFZ and f NFZ are the center coordinates of the NFZs and r NFZ is the radius.
When it comes to uncertainty, things have changed. Due to the initial estimation error or the movement of NFZ R , the real NFZ R may appear anywhere in a certain area NFZ P (As shown in Figure 1). When it becomes inevitable to cross this NFZ P , a very natural idea is to turn the original hard NFZ constraint into a soft one: Reducing the risk of passing through the NFZ R (which means mission failed) to a certain level. Then we have the following optimization goal: completing the original task based on satisfying the specific probability of success requirement. In this paper, a chance constraint is used to describe the probability of success requirement, which can be derived from the preceding hard NFZ constraint (32): where PrfAg is the probability of an event A, PÃ is the probability of success threshold, which means that the probability of not-colliding with the NFZ R must be greater than or equal to PÃ. Thus, the original problem turns into a CCOPT problem. In past research, the CCOPT problem can be solved in a variety of ways, such as the Monte Carlo method, particle swarm method, etc. But most of them are too slow to be used for an online algorithm. For the sake of computation efficiency, we propose an analytic approximation strategy as described in Section 3, which has been proven to be effective and highly efficient.

ABHV trajectory planning problem with uncertain NFZs
The optimal trajectory planning problem of ABHV with uncertain NFZs is formulated as the following optimal control problem.
The details of the performance index J will be discussed in Section 3.2.

Trajectory optimization for air-breathing hypersonic vehicle
Considering the lack of a widely recognized ABHV trajectory optimization method so far, this article first proposes a SOCP-based hybrid trajectory planning algorithm before introducing uncertain NFZ.
A typical SOCP has the following form x + e i is called secondorder cone constraint and its dimension is n i .
P0 in Section 2.4 contains highly nonlinear dynamics, nonlinear path constraints, and nonconvex control constraints that cannot be directly solved with the framework of SOCP. Blackmore and Liu [25][26][27] introduced lossless convexification techniques and relaxation methodology that convert the original aerospace problem into the specific form required in SOCP. In this chapter, similar techniques are used. Of course, the existing framework 27 needs to be further modified to accommodate the specificity of the air-breathing vehicle problem.
Furthermore, a hybrid method that combines online and offline computing is proposed to balance computational complexity and accuracy.
A hybrid method P0 contains a four-dimensional control input vector u 0 = ½a coss sin s h, and the elements in this vector are coupled with each other, as shown in equations (15)- (18).
Solving such a complex optimal control problem in a given time is challenging. To reduce the difficulty of computation and improve efficiency, a hybrid method is introduced: We hope ABHV can reach the maximum range if needed, which can be approximated in a discretized form where s is flight range, and t togo is time-to-go. The entire flight process is divided into n-1 parts and i represents the i-th discrete point. Dm i is the mass change between each discrete point. Equation (40) reveals that for a given fuel consumption, maximizing the total range is equivalent to maximizing the indicator and h directly affects this indicator in the cruising phase, it is reasonable to give them higher calculation priority than coss and sin s. Besides, if we solve a and Table 1. The procedure of the hybrid method. . This is a static optimization problem that can be easily solved by mature algorithms such as the trust-region reflection method. The result is the optimal combination (a i , h i )Ã that guarantees the maximum range capability. h first, the control input of the original problem will be decoupled, and the quantity of computation will be significantly reduced. According to the above ideas, a hybrid method is obtained (Table 1). First, in the offline part, the optimal a and h combination in the whole flight envelope is solved offline. Then, the actual optimal control problem is solved online.

Multi-objective optimization performance index
Without loss of generality, the performance index or cost function of ABHV is given in following form in this paper: where k is weight factor. In order to satisfy the form of SOCP, the performance index must be linear. Thus, the above equation must be reorganized before application.
Minimum terminal error. In the general trajectory planning problem, we hope that the vehicle will eventually reach a certain designated position. However, direct use of terminal constraint (22) is not good for algorithm convergence.
Thus, the terminal position constraints are transformed into a performance index where k is weight factor.
Maximum flight range. In a real mission, straight&level flight (s = 0) in the cruising phase probably does not keep the scramjet engine in the most efficient interval. In some cases, flying with a certain bank angle s (s 6 ¼ 0) may be a more sensible strategy. In addition, when ABHV has to circumvent obstacles like NFZs, we hope that the performance loss caused by bank angle s is minimal. So, use the indicator V i C f i h i mentioned in section 3.1 to minimize the range loss Equation (45) is noncovex, so it can be further expanded into a second-order Taylor series to fit the form requirement of SOCP as shown in Liu et al. 27 Minimum time. The performance index with time is Equation (46) also can be expanded into a secondorder Taylor series at each discrete point.

Convexification
Obviously, the problem P0 does not conform to the general SOCP form (39), so it still needs to be convexified and simplified. The process is given below.
Convexification of control constraints. As shown in Table 1, the only two control inputs that need to be handled are u 0 (2) = coss and u 0 (3) = sin s. Let us redefine the control input u = ½coss sins. Liu et al. 27 introduced an effective control constraint relaxation method for typical HGV, which is also suitable for the problem of this article. As the active assurances of the relaxed control constraints are already proved independently in Liu et al., 27 the detailed proof is not presented here. As a result, the control constraint (12) is relaxed into the following form: Liu et al. 27 assures that this condition will hold in the solution of the original problem. Obviously, the restructured constraint (47) is second-order cone and is suitable for SOCP.
Convexification of path constraints. In Section 2, the heating rate, dynamic pressure, and overload constraints, (a) Discretize the state equation of the optimal control problem. Generate an initial trajectory guess x 0 at the beginning of the iteration. (b) At each discrete point we calculate the (V, r, m) along the initial trajectory, then the (a, h) command can be achieved by interpolation from the mapping table. Substitute (a, h) into state equation. As a result, the only unknown variable that exists in the state equation is ½coss sin s. Solve the problem using the convex optimization method like the SOCP described in Section 3.5. (c) Update the initial trajectory with the result x 1 , go back to b), repeat the above process until we find an x k (k = 0, 1, 2, . . . ) which fully meets the optimization goal.
which are regarded as path constraints, were restructured into altitude constraint (26). The max function in r(m) ø max (r _ Q (m), r q (m), r n (m)) is convex, so the discrete form (aboutm i ) is linear and convex wherem i corresponds to the i-th discrete point.

Linearization of dynamics
The nonlinear dynamics must be approximately linearized before applying the optimization algorithm. To improve the convergence and robustness of the algorithm, the computing process is designed as an iteration: A pre-designed linear trajectory is used as the initial guess reference trajectory at the first iteration, then the reference trajectory will be updated with the result of k-th iteration x(m) k , u(m) k È É at the beginning of the k + 1th iteration. The pre-designed trajectory can be a linear trajectory and its feasibility does not affect the result.
Decouple the right side of equation (35) and rewrite it in a concise form 25 where f 0 = V sin g=C f h (V 2 À 1=r) cos g=(rVC f h) V cos g sin c=(rC f h cos f) V cos g sin c tan f=(rC f h) cos g sin c tan f=(rD) ½T cos a À D À sin g In the offline part, the optimal a and h combination have been already obtained. Thus, L, T, F Ty , and D in equations (51) and (52) are determined only by the flight status x.
To solve the dynamic problem, the equation (51) must be linearized. Let x k = x k (m) and u k be the solution pair in k-th iteration. So at the (k + 1)-th iteration, linearizing of equation (51) yields 27 where A(x k ) = ∂f 0 (x k ) ∂x k and b(x k ) = f 0 (x k ) À A(x k )x k . In order to ensure the reliability of linearization, a new constraint (54) including a trust-region column vector d 2 R 6 is added.
As a result, P0 is transformed into P1 with given (a, h) Ã that obtained in offline part. s:t: where k is the weight factor.

Discretization and the solution procedure of online part
In order to be solved by numerical algorithms, continuous systems in optimization problems are always discretized in practice. Of course, all the constraints must be satisfied at these discrete points. In this case, the independent variable ism. So if we discretize equation (53) into N + 1 points equally, then the step size ofm is Dm =m i Àm iÀ1 =m N Àm 0 ð Þ =N . Hence, the states in problem P1 can be denoted as x(m 0 ), x(m 1 ), x(m 2 ), Á Á Á , x(m N ) f g , and similarly, the control can be denoted as u(m 0 ), u(m 1 ), u(m 2 ), Á Á Á , u(m N ) f g . Let x i = x(m i ) and u i = u(m i ), the trapezoidal rule is applied to obtain the integral result from equation (53), which is given by where where I is a unit matrix with the same dimension as that of A.
By combining the state variable x i and control variable u i , together with the accessorial vector q and u (introduced in section 3.2.1), we can obtain the optimization variable vector z: As a result, the whole discrete equation is given by Obviously, the discretized constraints and performance index respectively illuminated from Section 3.2 through 3.3 are a linear function of z. The performance index is denoted by c T z where c is a constant vector. Thus, let a discretized version of problem P1 is given in the following: s:t: r min ł r i ł r max ð75Þ u All of the constraints in P2 are either linear or second-order, and the performance index is linear, too. Therefore, the problem P2 is SOCP, which can be solved with primal-dual interior point method. Note that the dimension of M is pretty large, directly depending on N. So superabundant discrete points will significantly depress the computational efficiency of the solution procedure. Although primal-dual interior point method shows its unique advantage in dealing with a large-scale sparse matrix like M, N still needs to be chosen properly. The details of how to solve P2 are illustrated in Table 2.
According to several test results and Liu et al., 27 the d in equation (73) Table 2. Details of the iteration procedure of P2 in the online part.
(a) Set k = 0 and build an initial profile Just integrate the motion equation without control until the value of the flight-path angle changed from negative to positive for the first time. Let state at this time be x 0 initial and set zero bank angle, perform linear interpolation aboutm between x 0 initial and x Ã f . The result will be the initial profile. (b) At the (k + 1)th iteration (k ø 0), use x k to update those corresponding parameters in problem P2 (e.g. M and F in equation (72), x k i in the trust-region constraint and r i min in the path constraint). Then, solve problem P2 to get a solution z k + 1 . (c) Check if the solution meets the stopping criterion: max )<e (81) If so, output the result. Otherwise, update z k with z k + 1 , set k = k + 1, and go back to (b).

Convexification of NFZ chanceconstraints
Single NFZ chance-constraint In Section 2.3, the uncertain NFZ is converted into a chance constraint optimization problem. These NFZ constraints are obviously concave so they cannot be substituted into the solving process directly. Whereas as shown in Liu et al., 27 if the concave constraints are discounted, the concave constraints can be approximated successively by its linearized form, for example: The hard NFZ constraint (32) is equal to where u (k) and f (k) are from kth iteration and Similarly, the chance constraint is used to describe the probability of success requirement, which can be derived from the preceding linearized hard NFZ constraint (82): Without loss of generality, the exact location of NFZ R is unknown, so let the center coordinates be individual Gaussian-distribution uncertain parameters.
where u NFZ and f NFZ corresponds to the mean of longitude and latitude, respectively. S represents the corresponding variances. Their value can be approximated using prior knowledge.
For the convenience of expression, let w = ½0 2(u NFZ À u (k) ) 2(f NFZ À f (k) ) 0 0 T and 1 = À (d 2 NFZ + d NFZ ) be the generalized uncertain parameter and determined parameter in equation (84), respectively. Obviously, w obey Gaussian distribution where the corresponding variance S NFZ = ½2S u NFZ 2S f NFZ , w is the mean of the new distribution.
Then equation (84) can be rewritten as then hence where F is the Gauss cumulative distribution function and thus which happens to be a second-order cone constraint for any 0:5 ł P Ã ł 1 (i.e. F À1 (P Ã ) ø 0). Hence, the chance constraint is converted into a hard boundary constraint. The inverse function F À1 (P Ã ) has no analytical form but can be calculated using a look-up method, or approximated as 28 Note that if those uncertain parameters are not individual, it will become a much more challenging problem and requires more undiscovered skills to be solved.

Multiple NFZ chance-constraints
Let Q be the number of NFZs. The multiple NFZ chance-constraint is given as where j = 1, 2, . . . , Q.
In order to apply the proposed analytic approximation method, we give the following relation In fact, equation (94) follows from this generic inequality max (x, y) ł ln(e x + e y ) Consequently, Pr ln implies that So inequality (96) is stronger than constraint (93). If (96) holds, the trajectory must satisfy the multiple NFZ chance-constraints. Note that the log-sum-exp function in equation (96) is convex and the detailed proof can be found in book. 16 As a result, the multiple NFZ chanceconstraints do not affect the convexity of the whole problem. To handle with (96), the log-sum-exp can be firstly linearized, then convert the result into a hard boundary with the technique in Section 3.6.1.

Numerical experiment
In this section, a series of numerical simulation demonstrations will be given to prove the effectiveness of the method described in this paper.
All these experiments are based on a generic research ABHV model called BABHV-1 (Figure 2), which is developed by Beihang University. It is a concept design for a 5-ton class reusable single-stage spaceplane, using a scramjet engine propulsion system. It has a tailless delta wing similar to SR-71. The purpose of this aircraft is to study the flight characteristics and control strategies of the air-breathing hypersonic vehicle at different speeds. Currently, a scaled model has been produced for the low-speed test flight. This design is capable of cruising at an altitude ranging from 24 to 32 km (Relative sea level) with a cruise Mach number of 5.0-7.5. As shown in Figures 3 and 4, the aerodynamic and thrust coefficients vary with flight conditions and throttle.
In all applications in this section, we select a value N = 200 in the discretization as described in Section 3.5. The Trust Region Reflective method is called in the offline part (as shown in Table 1). We also use the primal-dual interior-point methods (IPMs) in commercial optimization software for solving the SOCP problem in the online part (as shown in Table 2). All the results are obtained by running the proposed algorithm on a laptop with Intel Core i7-4720HQ 2.60 GHz.

Verify the accuracy of the algorithm
Before introducing the uncertain NFZ, a test is conducted to assure that the trajectory provided by the proposed algorithm is accurate and flyable. This verification is based on the following assumptions: We first simulate an open-loop trajectory randomly and then use the terminal conditions of the open-loop trajectory as the terminal constraint to test whether the task can   be reproduced by SOCP. If it can, it proves that the proposed method is accurate and can be trusted.
For this mission, the initial and the terminal condition are given in Table 3. The optimization index J is as follows The values of the limits in the path constraints in equations (24)-(26) are as follows: _ Q max = 500kW=m 2 , q max = 80000N=m 2 , and n max = 2:0g 0 . In addition, the upper bounds on the magnitude of the bank angle in equation (10)

Single uncertain NFZ
In this mission, we assume that there is an uncertain cylindrical NFZ R , the center coordinates of which are difficult to determine. The source of uncertainty can be    vague intelligence, modeling errors, or movement of the NFZ R itself. However, an initial guess of the center coordinates ( u NFZ , f NFZ ) can be given. Without loss of generality, NFZ R 's real center coordinates satisfy individual Gaussian distribution. By measuring the speed and heading of the object, the standard deviation (S 1=2 u , S 1=2 f ) of its position after a period of time can also be estimated. Unfortunately, the area where NFZ R may exist (NFZ P ) just blocks the way to target, and bypassing the NFZ P exceeds the vehicle's capacity. So it is inevitable to bear the risk and fly through the NFZ P . The probability of success threshold P* is set to 90%, which means that the probability of the vehicle not flying through the NFZ R and reaching the target should be greater than or equal to 90%. The performance index is which helps ABHV reach the target in the shortest time.
The initial and terminal conditions of the mission are given in Table 4. The target has a downrange of 792.2 km and a left crossrange of 107.6 km. The values of the limits in the path constraints are as follows: _ Q max = 500kW=m 2 , q max = 80000N=m 2 , and n max = 2:0g 0 . In addition, the upper bounds on the magnitude of the bank angle are set at s max = 60 deg. The parameters of uncertain NFZ are given in Table 5. Figure 8 shows the ground track. For comparison, the solution without uncertain NFZ is also obtained as the baseline. Both trajectories do fly through the NFZ P , but the solution that considering the uncertain NFZ is farther away from the ( u NFZ , f NFZ ) compares with the baseline one. To verify the correctness of the trajectory, the 68-95-99.7 rule is introduced. The baseline trajectory is quite close to the 1 ffiffiffi ffi S p boundary, which means the probability of success is only about 68%. Thus, it  does not meet the mission requirements. Solution 2 steers away from the NFZ center, grazes 2 ffiffiffi ffi S p the boundary, which means that the probability of success is slightly less than 95%, larger than the threshold P*. This figure illustrates the proposed algorithm can successfully handle the uncertain NFZ problem.
More comparison is given in Figures 9 to 13. In Figure 9 we can see the vehicle turns right then left to steer away from the uncertain NFZ in solution 2. Its trajectory is more curved than the baseline. Figure 12 shows that angle of attack and throttle changes with flight conditions. As the altitude decreases, the air density gradually increases, so the required trim angle of     attack decreases. This will result in a decrease in drag and the throttle required. The heating rate and dynamic pressure, and normal load histories are plotted in Figure 13. All two solution satisfies path constraints.

Multiple uncertain NFZs
In this mission, we assume that there exist two uncertain NFZs ( Table 6). The areas where they may exist happens to be on the path to the target. The vehicle's path constraints, initial and terminal conditions are the same as in section 4.3. Use the method described in section 4.2 to handle multiple chance constraints. The value of the probability of success threshold P* is still 90%.
The ground track of the solution subject to the multiple NFZs chance constraints is shown in Figure 14.
The trajectory without NFZ constraints flies through both NFZ P , whereas the trajectory with the uncertain NFZ constraints grazes the 90% boundary. In order to verify the correctness of the results, a Monte Carlo simulation was carried out. First, 1000 NFZs combinations that follow the distribution in Table 6 are randomly generated. Next, we check whether the trajectory violates the constraints formed by the these NFZs combinations. Table 7 shows that the trajectory has a 91.7% probability of not colliding with the real no-fly zone, which is higher than the set threshold. This result strongly supports the method in this paper.  The corresponding control input and states are shown in Figures 15 and 16.

Capability test
The previous section demonstrates that the algorithm can handle multiple uncertain no-fly zones. Then, further tests were performed to show how the shape and location of uncertain NFZs influenced the trajectory, in other words, how robust the new algorithm is.
The position of NFZ1 is changed to compress the vehicle's activity space (Table 8). We want to test whether the vehicle can still satisfy the chance constraint under more extreme conditions.
Compared to the previous simulation results, the vehicle performs a more violent maneuver when approaching NFZ1 to achieve a given success rate ( Figure 17). Similarly, a Monte Carlo simulation is performed for the new trajectory (Table 9), which guarantees the satisfaction of chance constraints (92.3% . p*).
In addition, we change the variance of the coordinate position of the NFZ's center (Table 10). Table 11 and Figure 18 demonstrate the same conclusion: the algorithm proposed in this paper has a strong robustness.

Conclusion
According to the flight characteristics of ABHV, this paper proposes a hybrid optimization method first. The numerical demonstrations show that this method can effectively solve the trajectory optimization problem of ABHV. Then we discuss the situation with uncertain NFZ. As a result, the uncertain NFZs that may exist are modeled as chance constraints. On this basis, a convexification method is proposed, which can convert the concave chance constraint into a deterministic constraint that can be solved by convex optimization.   Furthermore, it is further extended to situations where there are multiple uncertain NFZs.
In the past research, environmental uncertainty in hypersonic vehicle trajectory optimization was rarely discussed, let alone integrated into the optimization algorithm. The preliminary results obtained expand our understanding of convex optimization in complex aerospace engineering.