Abstract
The uncertainties during the return trajectory of vertical takeoff and vertical landing reusable launch vehicle weaken the ability of precision landing and make the return process more challenging. This paper is devoted to quantifying the probability uncertainty of return trajectory with uncertain parameters. The uncertainty model of return multi-flight-phase under the uncertainties of initial flight path angle, axial aerodynamic coefficient, and atmospheric density is established using the generalized polynomial chaos expansion method. By parameterizing random uncertainties and introducing random parameters into the uncertainty model, the uncertainty analysis problem of return trajectory is transformed into stochastic trajectory approximation problem. The coefficients of the polynomial basis function are solved by the stochastic collocation method. Then state solutions, statistical properties, and global sensitivity with Sobol index are established based on coefficients. The simulation results show the efficiency and accuracy of this method compared with the Monte Carlo method, the evolution process of main output parameters under random parameters, and relative importance for random parameters. Through the uncertainty analysis of the return trajectory, the robustness of return trajectory can be quantified, which is contributed to improving the safety, reliability, and robustness of recovery and landing mission.
1. Introduction
The technologies of retropropulsive return, descent, and pinpoint landing of vertical takeoff and vertical landing (VTVL) reusable launch vehicle (RLV) have received significant research attention in the past years [1] because of their advantages in reducing launch cost and improving safety and realizing reusability [2, 3]. So far, VTVL has been studied by many research institutions, especially industrial companies [4], and successfully demonstrated by SpaceX. However, due to the uncertainties in the return process and the failure of some actuators, some recovery missions of Falcon 9 still failed to land [5]. Guiding a rocket booster to return to and land on Earth with atmosphere is not a trivial mission [6], especially in the case of uncertainties, which weaken the ability of trajectory deviation correction and precision landing and make the return process more challenging [7]. Therefore, it is necessary to analyze the uncertainties of return trajectory to reduce flight risk, improve landing accuracy, and enhance emergency response.
In recent years, many effective strategies have been put forward for spacecraft, RLV, and VTVL trajectory planning, including trajectory optimization [8–13] and guidance and control [14–16] methods. Chai’s team have researched many trajectory optimization algorithms for space vehicle, including an initial guess generator using a violation learning differential evolution algorithm [8], a computational framework addressing the problem of stochastic trajectory optimization [9], a two-step strategy incorporating fuzzy multiobjective transcription and deep neural network [10], and a bilevel structure incorporating desensitized trajectory optimization and deep neural network [11]. Moreover, the convex optimization methods have been very popular in VTVL rocket in recent years. Liu used convex optimization to analyze the cooperative effect of thrust and aerodynamic forces on the fuel cost of return phases [12]. Szmuk et al. studied the fuel-optimal landing problem by combining lossless and convexification techniques [13]. Furthermore, convexification has also been used in guidance algorithm for VTVL trajectory planning. To solve the highly nonlinear fuel-optimal rocket landing problem, Wang et al. proposed a novel guidance algorithm based on convex optimization, pseudospectral discretization, and a model predictive control framework [6]. Scharf et al. reported a successful airborne implementation experiment using lossless-convexification guidance algorithm [14]. Ma et al. proposed a finite-element collocation based convexification method for powered landing guidance of RLV in the presence of aerodynamic drag forces [15]. Besides, Zhang et al. developed an extended state observer based nonsingular fast terminal sliding mode control with fixed-time convergence [16]. The optimization, guidance, and control methods mentioned above do not consider uncertainties. However, using reentry optimization or guidance methods alone is not robust enough to overcome the uncertainties and disturbances encountered during the flight, which may cause the rocket to fail to land at the scheduled landing point.
For return mission analysis and design, it is necessary to know the potential state of each flight phase and landing point state under the given separation state and control law, because it will help the designer to evaluate whether the scheduled engine start-up and shutdown points and landing point are appropriate and whether the mission needs to be replanned or updated [17]. Usually, the actual flight performance of launch vehicle is not as expected because of the uncertainty of some important return dynamics parameters, such as separation point state, aerodynamic parameters, and atmospheric density. Before we further develop future VTVL RLV retropropulsive return, descent, and pinpoint landing technologies, the challenges of uncertainty quantification must be addressed. Therefore, knowing the evolution law of return trajectory under uncertainty is beneficial to mission design and analysis, especially for emergency return mission.
In fact, the uncertainty propagation analysis of nonlinear dynamics has always been a research focus in scientific and engineering fields. In this paper, the uncertainty propagation problem of return trajectory is considered to deal with the uncertainty of random parameters with known probability distribution. Generally speaking, there are mainly three kinds of uncertainty analysis methods: Monte Carlo (MC) method, linear method, and nonlinear method [18]. MC method is relatively easy to implement, but it requires a large number of samples to accurately obtain statistical information, which makes the calculation cost expensive. So it is usually used as a reference benchmark to verify the accuracy of linear and nonlinear propagation tools [19, 20]. The linear uncertainty propagation method can perform well in the treatment of linear system problems [21] but not well for the high-precision entry dynamics system due to the existence of small disturbance hypothesis. The method based on generalized polynomial chaos expansion (GPCE), as a nonlinear uncertainty propagation analysis method, has been successfully applied to many space missions because of its convergence, accuracy, and computational efficiency in orbit and flight dynamics [22–24]. The current return trajectory uncertainty analysis is mainly applied to the vehicles during Earth reentry [25–27] and Mars entry [23, 24, 28, 29]. In order to deliver an RLV to its landing site safely and reliably, trajectory robustness theorem with uncertainties is applied to the trajectory generation in the reentry to Earth [25]. Uncertainties analyses of the initial state and atmospheric density were also introduced to improve error modeling and estimation in short-term reentry predictions [26]. Frayssinet et al. researched drag coefficient stochastic uncertainties in design of atmospheric hypersonic reentry vehicle and drag coefficient was evaluated for each random shape [27]. Many methods of uncertainty analysis have been also used for Mars entry dynamics to improve the accuracy and efficiency. For example, Jiang and Li developed an innovative approach for planning the Mars entry trajectory under uncertainties considering both dynamics parameter uncertainty and initial state uncertainty [28]. Bhattacharya’s team developed some excellent uncertainty analysis approaches for hypersonic flight dynamics [23, 29]. Their team employed the GPC method, in which the vehicle’s states were approximated, and the resulting stochastic dynamical system was converted into a deterministic dynamical system in a high-dimensional space via the Galerkin projection. As mentioned above, although there have been many researches on the design of VTVL RLV return trajectory, there are few references on the evolution of VTVL RLV return trajectory under uncertainties. In this paper, the uncertainty model is established based on GPCE method and the parameter uncertainty is analyzed. MC method is also used for comparison with GPCE and analyzing the state probability distribution of some main output parameters at the connection points between flight phases.
The main contribution of this paper lies in the uncertainty analysis of the recovery and landing mission of VTVL RLV, in which the uncertainty propagation problem is converted into the stochastic trajectory approximation problem. Stochastic collocation (SC) method, which is one of the nonintrusive methods [24], is used to solve the coefficient of polynomial basis function. The uncertainty analysis based on GPCE will be further developed in this paper to analyze the propagation process of multi-flight-phase trajectory under the uncertainties of initial flight path angle, axial aerodynamic coefficient, and atmospheric density. In the framework of GPCE, the state equation of stochastic space is expressed as a function of random variable and approximated by a set of orthogonal polynomials. Due to the connection points state between flight phases, the start-up and shutdown conditions of the engine, and the landing point state playing an important role in the return trajectory [30, 31], this paper focuses on the influence of the uncertain parameters on them. Because the constraints, including path and control constraints, also have a certain impact on the return trajectory, they are also considered inuncertainty analysis. With the evolution results of the return trajectory under parameter uncertainties, some important mission design parameters, such as height, velocity, longitude, and latitude, can be evaluated and determined. Moreover, the relative importance for uncertain parameters can be quantified by introducing global sensitivity with Sobol index based on the GPCE. Therefore, the uncertainty analysis of return trajectory can provide an assessment tool for mission planning.
The structure of this paper is as follows. The multi-flight-phase division and model of return trajectory are introduced in Section 2. The uncertainty model and detailed uncertainty analysis of return trajectory based on the GPCE method are presented in Section 3. To show the application of uncertainty model and analysis, the numerical simulation results are illustrated in Section 4. Finally, the conclusions are given in Section 5.
2. Return Multi-Flight-Phase Problem Formulation
2.1. Division of Return Multi-Flight-Phase
In this paper, the trajectory analysis under uncertainties is only for the return process of the first stage of VTVL RLV, so the ascending stage is not considered. Concerning the recovery pattern of launch vehicle, two distinct mission profiles are considered: downrange landing (DRL), in which the first stage often lands close to the drone ship on the sea site, and return to launch site (RTLS), where the first stage uses additional firing to return to its launch site [32, 33]. The DRL is considered in this paper. A typical VTVL rocket, which does not return to the original field, usually experiences multiple flight phases in the return process, such as attitude-adjusting phase, high-height-sliding phase, powered-descending phase, aero-braking phase, and vertical-landing phase [34]. Considering that the attitude-adjusting phase and high-height-sliding phase are mainly outside the dense atmosphere, the aerodynamic influence can be ignored and they can be combined into one flight phase. So the model of the return phases division is simplified in this paper, including four flight phases: attitude-adjusting phase, powered-descending phase, aero-braking phase, and vertical-landing phase. As shown in Figure 1, the recovery mission profile from separation point to touchdown point is given.
After the rocket’s ascent, the first stage cuts off its engine and then separates (Point 1 of Figure 1). The second stage then ignites the engine and flies toward the payload’s destination orbit. Meanwhile, during the attitude-adjusting phase, the first stage adjusts its attitude by the reaction control system (RCS) to make the engine face the flight direction, and velocity actually decreases when approaching the apogee. At a prespecified height or time, the first stage reignites its engine (Point 2) for the deceleration; then the engine is shut down at a specified height. During aero-braking phase (Point 3 to Point 4), the grid fins are mainly used as the actuator to guide the vehicle to fly to a predetermined position above the landing site at a certain angle. Finally, the engine is reignited (Point 4) at a prespecified height to enter the vertical-landing phase, which brings the vehicle from its current position and velocity to a soft touchdown at the drone ship. During powered-descending phase and vertical-landing phase, the vehicle mainly relies on RCS and grid fins to adjust and control its attitude. In addition, for a more effective management of aerothermal requirements, recovery guidance can be explicitly divided into three parts: (1) a return burn aimed at decelerating the vehicle, (2) a second engine cutoff, and (3) a landing burn to ensure a precise touchdown [7].
The slender first stage must be equipped with a variety of actuators such as RCS and grid fins to ensure the safe recovery of the vehicle. Since the first stage should be controllable not only in the exoatmospheric part (mainly the attitude-adjusting phase) of the return trajectory (with RCS) but also in the aerodynamic environment, the aerodynamic control devices such as grid fins are necessary. According to the characteristics of the actuators, the grid fins work continuously, while the RCS system works discretely. So it is necessary to distribute the torque control command from the control system to the actuators and then calculate the control signals of each actuator [35]. Grid fin is a kind of honeycomb structure composed of many grids. When the engine is shut down, the grid fins can also operate independently for attitude control. Therefore, it is necessary to install not only the RCS but also the aerodynamics control devices such as grid fins on the first stage, so that the aerodynamic forces can be controlled well and the precise control can be achieved during descent. Furthermore, the pinpoint landing of the first stage requires some kind of landing gear or landing legs [33].
2.2. Dynamic Model of Return Multi-Flight-Phase
The attitude control and adjustment process are not considered in the motion of the center of mass. The three-degree-of-freedom (DOF) dynamic equation of the return multi-flight-phase for the vehicle can be derived from [6, 30], which can be expressed by the following differential equations:where these variables are illustrated in Figure 1, is the launch coordinate system, is the body coordinate system, and is the velocity coordinate system. are the position components in the launching coordinate system, is the magnitude of velocity, is the flight path angle, and is the yaw angle. The former three variables describe the position of the vehicle in the launching coordinate system and the latter three variables describe its velocity. is the current mass, is the Earth radius, is the radial distance from the center of the Earth, is the engine thrust, is the specific impulse of the engine, is the gravity acceleration at sea level, ( is the Earth gravity constant), is the attack angle, and is the sideslip angle. are the aerodynamic components in the velocity coordinate system, which are defined as follows:where are the axial force coefficients, normal force coefficients, and lateral force coefficients, respectively, and . is the vehicle reference surface area. is the current atmospheric density of the Earth , is the reference density, and is the reference height.
Let and be the state vector and control variable vector in equation (1), respectively. Then equation (1) can be written in the state space as follows:
2.3. Control and Path Constraints of Return Multi-Flight-Phase
The trajectory solution of equation (3), which represents the return state of the vehicle at time under the given initial state , needs to be feasible with the given control law.where represent attitude-adjusting phase, powered-descending phase, aero-braking phase, and vertical-landing phase, respectively. are discrete points of each flight phase, and the values of and can be obtained by linear interpolation between discrete points. is variable thrust coefficient of each flight phase. is the rated thrust of engine. are the lower and upper bounds of attack angle, sideslip angle, and thrust in each flight phase, respectively. are the start and end times of each flight phase.
The maximum overload the vehicle can bear during return flight phase should be considered:where are the thrust components in the velocity coordinate system and is the maximum overload of each flight phase the vehicle can bear.
The upper limit of dynamic pressure mainly depends on the strength of thermal protection material and the moment of pneumatic hinge:where is the upper limit of dynamic pressure of each flight phase.
In addition, the mass of the surplus propellant is limited because the vehicle has gone through several flight phases before landing. So it is necessary to ensure that the landing mass at the end of vertical-landing phase is greater than the structural mass:where is the structural mass of the first stage, that is, the lower limit of residual mass during return process.
What is more, each of the four phases in return trajectory is linked to the adjoining phases by a set of linkage conditions. These constraints force the state to be continuous and also account for the mass ejections, as
3. Establishment and Analysis of the Uncertainty Model for Return Trajectory Based on GPCE
The basic principle of GPCE is that, according to the distribution type of random variables, the stochastic state is expressed as the weighted sum of corresponding orthogonal polynomials, in which the weight function and the probability density function (PDF) of random variables have the same form. The stochastic state solution will be represented by basis function and coefficient of GPCE approximately. In this section, the uncertainty model of the first stage return mission is established based on the GPCE method. Then, the coefficient of GPCE is calculated based on SC method and the approximate solution of stochastic return dynamics is obtained. Finally, the statistical characteristics and Sobol index of trajectory state can be directly obtained through the coefficient of trajectory approximate solution.
3.1. Expansion of Generalized Polynomial Chaos
The return trajectory starts from the separation point, so the state of the separation point has an important influence on the trajectory. The key variable that affects the return trajectory is the flight path angle at separation time, which can reflect the changing trend of rocket flight, and is an important parameter in trajectory design [36]. Moreover, as mentioned in [24, 37, 38], the aerodynamic coefficient and the atmospheric density also have an important influence on the trajectory solution for the vehicle returning to the Earth’s atmosphere. On account of the first stage return being similar to ballistic return and the attack angle and sideslip angle being small for the return process of DRL, we can see from equation (2) that the influence of axial aerodynamic coefficient on the trajectory is far greater than that of normal aerodynamic coefficient and lateral aerodynamic coefficient. So the influence of aerodynamic coefficient on the trajectory solution can be considered as that of axial aerodynamic coefficient on the trajectory solution. Therefore, this paper mainly analyzes the uncertainty of the initial flight path angle at separation point, the axial aerodynamic coefficient , and the atmospheric density . For the convenience of expression, the symbol is used to replace in the following equation derivation. Therefore, the random variable dimension of the return model is three.
In practice, uncertainties in usually have their own lower and upper boundaries. The random variables with parametric uncertainty need to be parameterized in order to apply GPCE. These random variables , , and can be written aswhere , , and are the perturbations about the nominal values , , and respectively; , and are the corresponding independent random variable parameters. For the sake of generality, the interval of random variables is and all of them obey uniform distribution in this paper.
Let be the three-dimensional random variables vector. After introducing random variables in equation (3), the state equation of the determinate space will be transformed into the uncertain dynamic state equation of the stochastic space:
According to the theory of GPCE, the solution of each state variable in equation (10) is a stochastic process and can be approximately expressed by finite series as given bywhere is the basis function of orthogonal polynomials expanded for three-dimensional random variables, is the order of , is the coefficient of , and is determined by the dimension and maximum order of random variables. The number of terms of polynomial expansion can be derived by
It has been proved in [39] that, with approaching infinity, the approximate solution will converge to , so the accuracy of the approximate solution can be improved by increasing .
For the stochastic space composed of three-dimensional random variable vector , because the random variables are independent of each other, the polynomial basis function can be calculated by tensor product method through univariate polynomial basis functions , , and :where , , and are the orders of , , and , respectively, and satisfy
According to the relationship between the distribution type of random variables and Wiener-Askey hypergeometric orthogonal polynomials, Legendre orthogonal polynomials should be selected as the polynomial basis functions corresponding to uniformly distributed random variables to ensure the convergence speed of GPCE [39]. It should be noteworthy that if the random variable is subject to other types of probability distribution, such as Gaussian distribution, the interval of the random variable and the corresponding orthogonal polynomial basis function need to be replaced. The corresponding relationship between the random variable and the basis function can be found in [39].
3.2. Solution of Coefficients Based on SC
Considering that the basis function is determined by the distribution type of random variables, it can be seen from equation (11) that once the coefficient solution of is obtained, the approximate solution can be obtained. Therefore, the key of GPCE is to determine the coefficient of the polynomial basis function. In this paper, the SC method, which is one of the nonintrusive methods, is used to solve the problem of coefficient determination.
Any state coefficient can be obtained with Galerkin projection method [40]. According to equation (11), the stochastic mode in the expansion of is given bywhere represents the inner product.
By the virtue of orthogonality,
Replace with :
Thus, it can be seen thatwhere denotes the expectation operator.
According to the orthogonality of polynomial basis function, every one-dimensional orthogonal polynomial basis function satisfieswhere , and are the orders of and , respectively, is the PDF of orthogonal polynomials , is normalization factor, and is Kronecker delta function defined as
Therefore, in equation (15) can be expressed as
Let be substituted into equation (18); it can be seen thatwhere is the PDF of .
The SC method is to use the weighted integral operation to approximate the integral at the right of equation (22); that is,where is the collocation point generated for three-dimensional stochastic space, is the number of stochastic collocation points, and , , , and are the weights of , , , and , respectively.
For equation (23), the coefficient of GPCE can be obtained as long as the solution , polynomial basis function, and corresponding weight at each collocation point are determined. For a given collocation point and initial conditions, the solution is a deterministic problem, which are decoupled from each other. The collocation points of the SC method are chosen as the Gaussian quadrature points, and the references [40, 41] give points selection strategy details by quadrature rules. For the stochastic space composed of random variables vector , the number of collocation points corresponding to each one-dimensional random variable is , , and , respectively, and the corresponding stochastic collocation points are recorded as , , and , so all collocation points calculated by tensor product method can be written aswhere the number of total collocation points of tensor space is derived as . If the number of collocation points of each one-dimensional random variable is , the total number of collocation points in the three-dimensional stochastic space is .
As mentioned in Section 3.1, the type of random variables is uniform distribution, and the corresponding polynomial basis function is Legendre orthogonal polynomial, which satisfy the recurrence relation:with and . The term is the order of . So can be obtained by introducing the order and collocation point into equation (25).
According to the type of random variables with uniform distribution, the weights , , and can be obtained by Gauss Legendre equation, which is based on the following formula [41]:where is the number of collocation points corresponding to random variable . denotes the coordinate of the collocation point of , and is the corresponding weight. denotes the remainder of the quadrature. The weight is given bywhere is the first derivative of .
3.3. State Solutions and Statistical Properties Based on Coefficients
The approximate solution of equation (10) can be obtained after the coefficient of the approximate solution of the stochastic dynamic equation is obtained, so that the connection condition in stochastic space between the flight phases of equation (8) becomeswhere are the coefficients of basis functions of , respectively.
The start-up condition of the engine in powered-descending phase is time start-up, and shutdown condition is height shutdown. The shutdown condition under uncertainty can be written aswhere is the coefficient of at time in powered-descending phase.
The start-up condition of the engine in vertical-landing phase is height start-up, and the condition under uncertainty can be written aswhere is the coefficient of at time in vertical-landing phase.
The shutdown condition of the engine in vertical-landing phase is velocity shutdown, and the condition under uncertainty can be written aswhere is the coefficient of at time in vertical-landing phase.
The statistical properties including mean and variance of return trajectory solution can be directly obtained by polynomial coefficients [23, 42]. The mean value of can be obtained by
The variance of can be obtained by
3.4. Variance-Based Global Sensitivity Analysis with Sobol Index
Through uncertainty analysis, we can get the influence of uncertainty factors on the output parameters. But, in order to quantify the relative importance of uncertainty parameters inputted to the model, it is necessary to carry out parameter sensitivity analysis. Because of the expanded range of inputted uncertain parameters and the strong nonlinearity of the model in this paper, the local sensitivity would produce a large deviation. So, global sensitivity analysis with Sobol index based on variance is used to analyze the influence of inputted uncertainties on output parameters in this paper. Sobol index has been widely used in the fields of mathematics and engineering mechanics. The traditional Sobol index is obtained by a large number of Monte Carlo simulations, but it is difficult to apply to the calculation of complex model. The nonintrusive polynomial chaos expansion (NIPC) method has been used to calculate Sobol index [43, 44]. In this section, an efficient calculation method for Sobol index of return landing problem is given by conducting nonintrusive GPCE method. The sensitivity of landing height, velocity, longitude, and latitude to the three input uncertain parameters is analyzed comprehensively for VTVL RLV landing problem.
As described by Liu et al. [43] and Zhao et al. [44], the total variance can be decomposed into the sum of partial variances of different uncertain parameters:
Dividing both sides of equation (34) by , we can getwhere global sensitivity index is defined as
In this paper, the uncertainty dimension ; the items in equation (34) can be expressed as
The items in equation (35) can be expressed as
From the above analysis, considering the nonlinear correlation between input variables and output variables concerned, it can be seen that Sobol index can be used to provide the relative importance of each input uncertainty parameter to the overall change of output parameters. The aim of the current work is to use GPCE to calculate Sobol index and then evaluate the relative importance of each input uncertainty parameter.
3.5. Analysis Process
In this paper, polynomial chaos expansion is used as a high-order alternative model, which has good convergence, calculation accuracy, and calculation efficiency [23, 39]. The coefficient of expansion term can not only calculate the statistical characteristics such as mean value and variance but also directly calculate the sensitivity index. To sum up, the procedure of the uncertainty analysis based on GPCE is as follows. Step 1: Let the stochastic space contain three random variables , , and . For each random variable, according to its probability distribution type (uniform distribution in this paper), generate collocation points by quadrature rules. Step 2: Calculate three-dimensional collocation points ; calculate the polynomial basis functions according to corresponding collocation points and equation (25) and integral weights according to corresponding collocation points and equation (27) in one-dimensional stochastic space. Step 3: Input and initial condition into equation (11) for solution, and the values at each stochastic collocation can be outputted; calculate weights and the polynomial basis functions in three-dimensional stochastic space using tensor product method. Step 4: The output components , the corresponding weights , and the basis functions are substituted into equation (23), and the coefficient of GPEC can be calculated. Step 5: Analyze the uncertainty propagation of return trajectory by calculating the approximate solution of stochastic dynamic system according to equation (28), the mean value and variance characteristics according to equations (32) and (33), and the global sensitivity index according to equations (36)–(38).
As shown in Figure 2, the flow chart of uncertainty analysis is given.
4. Numerical Demonstration
4.1. Simulation Parameters
The computational assumptions are made as follows:(i)The range of attack angle and sideslip angle is in all flight phases except for the attitude-adjusting phase [6](ii)Only gravity is considered in the attitude-adjusting phase(iii)Considering that the attitude-adjusting phase is mainly in the exoatmosphere, wind disturbance is not considered in this phase(iv)The engine is installed symmetrically, and the thrust is along the axis of the vehicle body(v)The uncertainties of the three parameters are uniformly distributed
The simulation parameters are shown in Table 1.
4.2. Simulation Results and Analysis
In this section, the results of uncertainty propagation in a specific return scene are presented by using GPCE and MC method, mainly the influence of uncertain parameters on output landing parameters including height, velocity, longitude, and latitude. As shown in [31], specific landing conditions of height and velocity are prescribed, which have to be fulfilled:(i)Landing velocity: 0 m/s − max. 2.5 m/s(ii)Landing height:
Three kinds of uncertainty factors are considered, that is, initial flight path angle, axial aerodynamic coefficient, and atmospheric density. The deviations of the three factors are uniformly distributed, and, as referred by Huang et al. [24] and Long et al. [37], the upper and lower bounds of , , and are set as deg, of , and of , respectively. We will analyze the influence of multiple uncertainties coupling and single uncertainty on the return trajectory of the vehicle. In analyzing the propagation of the deviation of each uncertain factor, 6 collocation nodes are used to generate random samples. At the same time, in order to balance the calculation efficiency and cost, the 4-order Legendre orthogonal polynomials are used to obtain the GPCE approximate solution coefficients [24].
As shown in Figures 3–6, the variance evolution results calculated by GPCE are compared with those of MC, and it can be found that the results match well. A good convergence result is obtained by GPCE with 4-order Legendre orthogonal polynomials. The small calculation errors of GPCE prove that this method has a high approximation accuracy for the return trajectory dynamics. In addition, it should be noted that, under the simulation conditions in this section, the MC method needs 2500 groups of samples to achieve convergence. However, GPCE method only uses 6 samples under single uncertainty, with a time of about 8 s, and uses groups of samples under three composite uncertainties, with a time of about 230 s. But, with the same accuracy, MC method takes about 2800 s under 2500 samples. Compared with MC, GPCE improves computing efficiency by two orders of magnitude under a single uncertainty and one order of magnitude under three composite uncertainties. Therefore, GPCE method shows an obvious improvement in calculation efficiency under the premise of ensuring accuracy. In general, GPCE, which is based on the “black box” dynamic sampling method, can solve the statistical characteristics of the output response, such as variance, directly from the coefficients. Compared with MC method, GPCE has faster convergence speed, less calculation time, smaller sample, and higher approximation accuracy. The advantages of GPCE are obvious.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
4.2.1. Return Multi-Flight-Phase Trajectory Analysis under Composite Uncertainties
The stochastic deviation evolution under the coupling of initial flight path angle, axial aerodynamic coefficient, and atmospheric density is firstly considered. As shown in Figures 3 and 7, with the change of uncertain parameters, the region formed by the envelope of height, velocity, longitude, and latitude in different flight phases, variance evolution, and the probability distribution of these main output parameters at the end of each flight phase are changing accordingly. Therefore, the simulation results show that the uncertainties of three composite random parameters have an important influence on the evolution of height, velocity, longitude, and latitude.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
Since the attitude-adjusting phase has no thrust effect and does not consider the influence of aerodynamic force, the propagation of uncertainties of initial flight path angle, axial aerodynamic coefficient, and atmospheric density in dynamics is a linear process only under the action of gravity. Therefore, it can be seen that the variation of the height, longitude, and latitude variance of the attitude-adjusting phase with time is a parabola in Figure 3. Besides, the velocity variance is almost constant with time, approximately close to zero in Figure 3(b). It can be seen that the probability distribution of velocity and height in Figure 7(a) and longitude and latitude in Figure 7(e) at the end of attitude-adjusting phase is approximately linear. In the following analysis of trajectory propagation of initial flight path angle uncertainty, similar conclusions can be obtained, and no special explanation will be given. In powered-descending phase, it mainly depends on the engine thrust to decelerate, and the height variance increases first and then decreases. At the end of the powered-descending phase, the shutdown condition is height shutdown. Meanwhile, the impact of uncertainty will lead to a large deviation from the nominal value when reaching the same height, and the velocity, longitude, and latitude variance would suddenly increase. In the aero-braking phase and vertical-landing phase, under the joint control of aerodynamic force and engine thrust, the overall variance of velocity and height presents a downward trend, and the variance of longitude and latitude tends to be stable. Similar conclusions can be obtained in the subsequent separately analysis of the impact of initial flight path angle, axial aerodynamic coefficient, and atmospheric density uncertainties on the trajectory.
Figure 7 shows the probability distribution of height and velocity and longitude and latitude at the end of each flight phase under the condition of composite uncertainty. It can be seen from Figure 7 that there is velocity or height deviation at the end of each flight phase, especially at the landing point of the vertical-landing phase in Figure 7(d). The engine is shut down according to the specified velocity, but there is a great uncertainty of about 2800 m in height, which is seriously inconsistent with the landing requirements in [31]. There is uncertainty of 0.1 deg in longitude and 0.2 deg in latitude, which leads to an obvious violation of the landing constraints on the return trajectory. In other words, in order to ensure that the landing parameters still meet the predetermined nominal trajectory landing constraints under uncertainties, it is necessary to reformulate the control law. In addition, the design of return mission should be more conservative after considering the influence of uncertainty on return trajectory. In the subsequent analysis of the influence of single uncertain factors on the trajectory, the above similar conclusions can be obtained, so they will not be repeated.
4.2.2. Return Multi-Flight-Phase Trajectory Analysis under the Uncertainty of Initial Flight Path Angle
The results of the height, velocity, longitude, and latitude evolution, variance, and probability distribution at the end of the flight phase caused by the initial flight path angle uncertainty are shown in Figures 4 and 8. At the end of the attitude-adjusting phase, the uncertainty of deg in initial flight path angle would lead to the uncertainty of the height about 4000 m, the uncertainty of the velocity about 15 m/s, the uncertainty of the longitude about 0.03 deg, and the uncertainty of the latitude about 0.1 deg. Because of the longer trajectory integration time and larger control action range, the uncertainty of initial flight path angle increases gradually and propagates more significantly in dynamics than in other flight phases. In the front part of the powered-descending phase, the height variance still shows an upward trend with time. But, under the joint action of atmosphere and thrust, the height deviation in the second half of the powered-descending phase decreases gradually, and the variance shows a downward trend. The velocity at the end of the powered-descending phase has an uncertainty of about 20 m/s. At the end point of the aero-braking phase, there is an uncertainty of about 3 m/s in the velocity when the engine is started at a certain height, resulting in an uncertainty of about 200 m in the height of the landing point. Compared with the following analysis of the uncertainty of the axial aerodynamic coefficient and the atmospheric density, the height of the landing point is more robust to the uncertainty of the initial flight path angle. There is an uncertainty of 0.1 deg in longitude and 0.2 deg in latitude of landing point, which are greater than the uncertainty of axial aerodynamic coefficient and atmospheric density. So the uncertainty of initial flight path angle has greater influence on the evolution of longitude and latitude.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
4.2.3. Return Multi-Flight-Phase Trajectory Analysis under the Uncertainty of Axial Aerodynamic Coefficient
The results of the trajectory evolution caused by the uncertainty of in axial aerodynamic coefficient are shown in Figures 5 and 9. The influence of aerodynamic force is not considered in the attitude-adjusting phase, so the axial aerodynamic force coefficient does not cause deviation propagation in this phase, and the variances of height, velocity, longitude, and latitude are close to zero as shown in Figure 5, in which case the attitude-adjusting phase is the nominal trajectory. In the first half of the powered-descending phase, because the height is greater than 80 km, the influence of thin atmosphere is very small, and the uncertainty of aerodynamic coefficient has little influence on the trajectory. When the height is lower than 80 km, the atmospheric influence increases gradually, and the deviation of the trajectory caused by the uncertainty of axial aerodynamic coefficient increases gradually. The variance changes of height, velocity, longitude, and latitude are more obvious in aero-braking phase and vertical-landing phase. Similar conclusions can be obtained in the following trajectory propagation analysis of the uncertainty of atmospheric density. There is an uncertainty of about 1400 m in height and 0.04 deg in longitude and latitude at the landing point under the uncertainty of axial aerodynamic coefficient.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
4.2.4. Return Multi-Flight-Phase Trajectory Analysis under the Uncertainty of Atmospheric Density
Under the uncertainty of in atmospheric density, the evolution results of height, velocity, longitude, and latitude are as shown in Figures 6 and 10. Compared with the uncertainty of the axial aerodynamic coefficient, output parameters including height, velocity, longitude, and latitude have similar robustness in uncertain propagation process. The landing height is less robust to the uncertainties of axial aerodynamic coefficient and the atmospheric density than to the uncertainty of initial flight path angle, but longitude and latitude are less robust to the uncertainty of initial flight path angle. For the problem of low robustness, it can be improved by changing the thrust of engine and controlling the deflection angle of grid fins.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
4.2.5. Sensitivity Analysis
The simulation in this section is devoted to the quantitative analysis of the relative importance of input uncertain parameters on main landing parameters including height, velocity, longitude, and latitude. Table 2 shows the global sensitivity calculation results of landing point parameters under the single uncertainty and the multiple uncertainties coupling of initial flight path angle, axial aerodynamic coefficient, and atmospheric density. It can be seen from Table 2 that initial velocity, axial aerodynamic coefficient, and atmospheric density have different effects on these landing parameters. In terms of single uncertainty, landing altitude and velocity are more sensitive to the uncertainties of axial aerodynamic coefficients and atmospheric density. However, the longitude and latitude are more sensitive to the uncertainty of initial flight path angle, whose Sobol index is two orders of magnitude larger than the Sobol index of axial aerodynamic coefficient and atmospheric density. The sensitivity index of composite uncertainties with initial flight path angle is also relatively large. Although the uncertainties of the axial aerodynamic coefficient and the atmospheric density have little influence on the landing longitude and latitude, the influence cannot be ignored due to the strong nonlinearity of the return dynamics of VTVL RLV. It can also be seen from the table that the sensitivity index of the coupling interaction of two or three uncertainties is greater than that of single uncertainty. The higher the uncertainty dimension is, the greater the impact on the landing point parameters is. By analyzing the sensitivity of the input random parameters, we can know more about the nature of propagation itself. This kind of analysis is closely related to engineering, especially for the robustness optimization of the return trajectory of launch vehicle and the reduction of the influence of the deviation evolution caused by the uncertainties.
4.2.6. Application in the Return Mission Design
The results of the uncertainty analysis are very useful for analyzing the robustness of the return mission under the uncertain conditions, which is helpful for setting path constraints, selecting separation points, and landing points in the mission planning stage. The following shows the function of uncertainty analysis on mission design through a simple process, which is presented in Figure 11 for better illustration. The goal of this return mission is to select a separation point and ensure that at least 97% of the probability of reaching the range of scheduled landing point is achieved under the above-mentioned uncertain conditions.
For the return mission in this section, it is assumed that the VTVL RLV model is as above. Therefore, although there is uncertainty in the process of return, the nominal path constraints and should be satisfied with at least 97% reliability. In order to achieve this goal, it is feasible to use the results of uncertainty analysis to adjust the separation point, as shown in the following steps.(1)Given the VTVL RLV design, the nominal path constraints and are used as the initial design constraints of and , where the superscript represents the number of constraints updated, and the subscript represents the number of separation point states modified.(2)To ensure that the return environment is suitable and the landing point is reachable, the vehicle usually has certain orbit maneuverability before separation. Therefore, the separation point can be modified by orbit maneuver. It should be noted that the initial separation point , is used in the first iteration as shown in Table 1.(3)Considering the parameter uncertainties of , these uncertain variables can be modelled as shown in equation (9).(4)Given the design constraints and separation point states, the above method can be used to analyze the evolution of landing point state under uncertain conditions. If the design constraints do not meet the requirements, the next iteration is required. When the constraints meet the reliability standard and the probability of landing point within the specified range is less than 97%, the separation point needs to carry out iterative correction of the outer loop . It should be noted that the sensitivity of landing site state to different uncertainties is different. Therefore, the sensitivity analysis of individual uncertainties will help to update design constraints and modify separation point. At the same time, the uncertainty analysis is also helpful to the trajectory optimization under uncertainties, which is beyond the scope of this paper and will be discussed in future research.
5. Conclusion
This paper introduces the method of establishing and analyzing the uncertainty model for the return trajectory based on GPCE, and MC method is used as a validation. The uncertainties of initial flight path angle at separation point, axial aerodynamic coefficient, and atmospheric density are considered in the uncertainty analysis of the return trajectory. Among the output parameters of the return trajectory, the main parameters including height, velocity, longitude, and latitude are discussed in detail during return process and at key points at the end of flight phases due to their significant roles in recovery mission of VTVL RLV. The propagations of the return trajectory under composite uncertainties and single uncertainty are separately analyzed. Through the uncertainty analysis, the satisfaction and robustness of these main parameters under the uncertainty conditions can be quantified. Besides, the relative importance of uncertainty parameters is quantified by global sensitivity analysis with Sobol index based on GPCE. The safety, accuracy, and reliability of uncertain return mission can be improved by introducing uncertainty analysis. In addition, the proposed method not only provides a reliable tool for the uncertainty analysis of return trajectory but also extends the application of GPCE method.
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
The authors are grateful to Professor Li’s advice on how to best improve this paper. This research was supported by the National Natural Science Foundation of China (No. 11372345).