A novel linear uncertainty propagation method for nonlinear dynamics with interval process

Interval process is a preferable model for time-varying uncertainty propagation of dynamic systems when only the range of uncertainties can be obtained. However, for nonlinear systems, except for Monte Carlo (MC) simulation, there are still few efficient uncertainty propagation methods under the interval process model. This paper develops a non-intrusive and semi-analytical uncertainty propagation method, named the “convex model linearization method (CMLM),” by constructing a linearization formulation of a nonlinear system in a non-probabilistic sense. First, the criterion to evaluate the difference between the original system and the linearization formulation is derived, represented by the discrepancy of middle point, radius and correlations of response. By minimizing these three parameters, the coefficients of linear equations will be optimized to obtain the linearization formulation of the original system. Then, analytical equations are built to calculate uncertainty response under the interval process, without time-consuming analysis of the original system. To further improve the efficiency of the linearization process, Chebyshev polynomial is introduced to approximate the nonlinear dynamic analysis. Two numerical examples of duffing oscillators and vehicle rides are set to test the proposed CMLM. Compared to the MC method, with comparable uncertainty response precision, the CMLM just needs 1–10% times of dynamic analyses of the nonlinear system. Furthermore, a practical launch vehicle ascent trajectory problem with black-box dynamics is solved by, respectively, the CMLM and MC method. The results verify the capacity of the CMLM to deal with black-box problems and show that the CMLM performs better in terms of accuracy, efficiency and robustness.


Introduction
Dynamic response evaluation of nonlinear systems plays an important role in various practical engineering. Conventionally, the evaluation is based on a deterministic condition; however, practically, any dynamic systems are in the presence of uncertainty. Under these uncertainties, the performances of the systems will be affected. Even a slight bias might cause a great change in the dynamic response. Therefore, in recent years, much research has recognized the importance of uncertainty problems and developed several uncertainty propagation methods for dynamic response evaluation [1][2][3].
Generally, the probabilistic model is widely applied in uncertainty problems, in which the uncertainties are quantified by stochastic variables with certain probability distributions. For uncertainty propagation of dynamic responses, some categories of probabilistic methods have been established, including Monte Carlo (MC) simulation [4], linear methods and nonlinear methods. MC simulation is capable of providing a result that approaches the true probability distribution but requires a high computational cost. Linear methods including local linearization [5,6] and statistical linearization [7] have advantages in simplicity and efficiency. Nonlinear methods, such as polynomial chaos expansion (PCE) [8][9][10], are developed for highly nonlinear systems. To apply all the aforementioned probabilistic methods, there is an essential premise that the precise probability distributions of uncertainties are available. However, in many practical engineering problems, due to limited experiment conditions or lack of knowledge, there are always restrictions on obtaining sufficient information of uncertainties. To remedy this limitation of probabilistic methods, many non-probabilistic methods have been developed [11][12][13][14].
Convex model theory [13] assumes that the imprecise parameters are enveloped into a convex set, and quantifies them by certain bounds rather than precise probability distributions. Thus, it is regarded as a powerful supplement to the traditional probabilistic method and, so far, has been widely investigated in a great number of uncertainty analysis areas. A growing number of researchers apply the convex model to evaluate the uncertain response of nonlinear systems, such as a wind turbine geared transmission system [15], rigid-flexible multibody systems [16] and an overhead crane [17]. Accordingly, various methods are developed. Wu et al. [18,19] introduced the Chebyshev inclusion functions into interval problems which can achieve high accuracy and efficiency. Li et al. [20] proposed a higher-efficiency sparse regression method to reduce the computational cost of the Chebyshev method. Fu et al. [21,22] proposed a polynomial surrogate method for rotor systems under interval uncertainties. Wang et al. [23,24] proposed a Legendre polynomial-based method and applied it to complicated multibody dynamic systems. Polynomial approximation approaches are shown to be effective to improve efficiency and guarantee precision for interval uncertainty propagation.
The underlying assumption of the aforementioned methods is that the uncertainties are fully independent and time-invariant. Practically, the uncertainties may be of strong or weak correlation or dynamic characteristics. To quantify the correlation of interval uncertainties, several improved convex models are proposed, including multidimensional parallelepiped models (MPMs) [25,26] and a multidimensional ellipsoidal model (MEM) [27]. These models describe the correlation of uncertainties through different sizes and shapes of convex models. Thus, they provide a potential approach to correlation analysis of dynamic responses under correlated interval uncertainties. Essentially, the response is a kind of time-varying uncertain variable, and its values at different times are correlated interval variables. In the probabilistic model, such similar time-varying uncertainty can be described as a stochastic process; however, there have been few similar effective mathematical tools for the convex model. To end this, Jiang [28] proposed a new non-stochastic quantification model for time-varying uncertainty, namely the ''interval process'' or ''nonprobabilistic convex model process.'' This model quantifies the time-varying uncertainties using lower and upper bound functions and describes the correlations of the values at different times using the MEM. Thus, it still presents the advantage that only the variation bounds rather than the precise probability distributions are required. Then the model is soon applied to dynamic response evaluation. The vibrations of single-or-multiple-degree-of-freedom systems and continuum structures are analyzed under the interval process [29]. Therefore, based on the interval process model, non-probabilistic methods for uncertainty propagation of dynamic systems under correlated uncertainties and time-varying uncertainties can be established. Up to now, most of the existing methods are only applicable to linear systems. For nonlinear systems, except for MC simulation [30], only a Karhunen-Loève (K-L) expansion method [31] has been proposed. In the method, an interval process is described by a superposition of infinite deterministic time-related functions with uncorrelated interval coefficients. For an interval process with weak time correlation and long duration, the K-L method may be computationally prohibitive due to dimension disaster.
In conclusion, the interval process model is a preferable non-probabilistic tool for uncertainty propagation of dynamic response. However, for nonlinear systems, there is a lack of an efficient method under the interval process model. In the probabilistic model, linear methods are known for their high efficiency. Among them, the statistical linearization method transforms a nonlinear system into an approximate linear formulation in some statistical sense. It presents advantages that the derivatives of dynamic systems are not required and the modification of the original systems can be avoided [32].
Inspired by the statistical linearization method, this work proposes an uncertainty propagation method, called the ''convex model linearization method (CMLM),'' for nonlinear dynamics under the interval process model. The following novelty is realized. 1) A novel linearization method is proposed for uncertain nonlinear systems under the convex model, which non-intrusively provides the optimal approximated linear formulations of nonlinear systems in a nonprobabilistic sense. 2) Based on the linearization method, a high-efficiency dynamic uncertainty propagation procedure is established for black-box nonlinear systems under the interval process. 3) The computational cost of the original nonlinear system in the procedure is further decreased by introducing the Chebyshev method.
The paper continues in Sect. 2 with a statement of the nonlinear dynamics with time-varying uncertainty under the interval process model. Section 3 introduces the proposed CMLM in detail. The methods are tested in two numerical examples and a launch vehicle (LV) ascent trajectory problem in Sect. 4. Finally, conclusions are given in Sect. 5.

Nonlinear dynamics with time-varying interval uncertainties
Consider nonlinear dynamics with time-varying interval uncertainties x I (t) and correlated initial interval uncertainties x 0 I . The expression can be represented as: where f(Á) are nonlinear functions; superscript I denotes interval; and B(t) is the corresponding input matrix. Under the interval uncertainties, the responses of the nonlinear system x I (t) are also time-varying interval variables. To quantify these interval variables, the following models are introduced to describe the correlated interval uncertainties and the time-varying interval uncertainties.

Correlated interval uncertainties description
Generally, interval variables are described by an interval model (IM) as shown in (a) and (c) of Fig. 1.
Since the IM is not effective to quantify the correlations between uncertainties, some other convex models are proposed. MEM is a well-performed model [33], as shown in (b) and (d) of Fig. 1. In the MEM, the correlation of the variables can be described by the shape of the ellipsoid. The ellipsoid of an n-dimensional MEM can be quantified as [27]: Fig. 1 Concepts of the interval model and the multidimensional ellipsoid model where C is the covariance matrix of the MEM: The covariance of two interval variables is defined as the geometric characteristics of the corresponding two-dimensional MEM [29], as shown in Fig. 2.
The value of covariance can be calculated as [29]: where h denotes a rotation angle of the ellipse from its normal state as shown in Fig. 2; meanwhile, r i and r j are the half lengths of the major axis and the minor axis of the ellipse. Naturally, the correlation coefficient [27] is defined as: The shapes of the ellipse that correspond to several different correlation coefficients are shown in Fig. 2. When | q |= 1, x I i and x I j are correlated and the ellipse degenerates into a line.

Interval process model
Based on the MEM, the interval process model [28,34] is proposed to describe a time-varying interval variable x I (t). At every time t, the value of x(t) is an interval variable. Thus, the variation of x I (t) is enveloped within a pair of upper bound and lower bound functions, x U (t) and x L (t), as shown in Fig. 3.
The middle-point function x M (t) and radius function x R (t) [29] are defined as: For any times t 1 and t 2 , the auto-covariance function of x I (t) is defined as the covariance [28] of x I (t 1 ) and x I (t 2 ) based on MEM, as shown in Fig. 3.
Similarly, the correlation coefficient function [28] is defined as: Moreover, for a set of interval processes, x I (t) = {x I 1 (t), x I 2 (t),…, x I n (t)}, the cross-covariance function matrix [29] is defined as:  where the elements c ij is the covariance of x I i (t 1 ) and x I j (t 2 ) based on MEM as: The radius function x R (t) can be obtained using the cross-covariance function matrix as: Finally, under the interval process model, timevarying uncertainties x I (t) in problem (1) can be described by interval processes, and correlated interval initial uncertainties x 0 I can be expressed as Under these uncertainties, the responses of the nonlinear system (1) are also interval processes represented by the upper bound functions x U (t) and lower bound functions x L (t), and the goal of uncertainty propagation is calculating x U (t) and x L (t). To achieve uncertainty propagation through the CMLM, the criterion to evaluate the degree of approximation of a linear formulation to the original system is derived at first. Then the basic theory of time-varying uncertainty propagation for linear dynamics is established. Finally, the complete uncertainty propagation process can be established.

Linearization of nonlinear dynamic system
Consider the nonlinear part of the system (1), f(x, u, t), under the interval process inputs, it becomes an interval process and its value at time t is an interval variable f I as: Now try to replace the nonlinear part with an approximate linear equation, which is selected in a way to minimize the difference between it and the original equation. The linear equation, under interval uncertainties, is also an interval variable at time t as: where x I = r I ? m; thus, r I is a vector of interval parts with zero middle-point values, and N T is the approximate linear state matrix of interval parts; meanwhile, m is a vector of middle-point parts with zero radius values, and M T is the approximate linear state matrix of middle-point parts.
If the two interval responses of nonlinear and linear equations to the same interval inputs have the minimum difference, the linear equation can be defined as the optimal approximation to the nonlinear equation in some non-probabilistic sense or called convex model sense. In the MEM, the differences between two interval variables can be evaluated by: where e 1 evaluates the difference in middle-point values, e 2 evaluates the difference in radius values and e 3 evaluates the difference in correlation. Thus, the optimal approximate linear equation represented by N * and M * can be obtained by minimizing e 1 , e 2 and e 3 simultaneously as: However, it is not convenient to calculate the correlation coefficient in the e 3 by the definition (5), because the geometrical properties of the MEM cannot be obtained analytically. The literature [33] has proposed the sample correlation coefficient (SCC) which is capable of calculating the correlation coefficient between any two interval variables without the constraint of the form of convex models: where x I i and x I j are two interval variables; and x i (s) and x j (s) are sampling points of x I i and x I j . q s = q, when the sampling points are uniformly distributed within a MEM.
Then, to obtain N * and M * , the muti-objective optimization problem (15) can be transformed into a single-objective optimization problem by weighted sum (see Appendix 1 for the determination of weighting coefficients) as: After simplification (see Appendix 1 for details), the expression can be obtained as: Minimization of e can be obtained when: Thus, N * and M * can be calculated as: The optimal approximate linear function can be expressed as: The approximate linear equation can be obtained by calculating the middle point f M and the radius f R of the nonlinear function, as well as the correlations q s (f I , x I ) between the interval response and interval inputs.
However, since the calculation of these values is based on the sampling points, a large amount of computational cost of nonlinear functions is required. To reduce the cost, the Chebyshev polynomial approximation of the original function is applied. A where Then, a k-dimensional function f(x) can be approximated by Chebyshev polynomials with order no more than n as: where the coefficients c can usually be determined by two methods. The first one is the Chebyshev tensor product method (CTPM) [35], and in the CTPM, the number of sampling points is no less than (n ? 1) k . The CTPM can achieve the highest precision of the approximation, but for high-dimensional problems, computational efficiency will be extremely low. Another alternative construction form for Chebyshev polynomial approximation is the Chebyshev collocation method (CCM) [35] which requires only a portion of samples from the CTPM. The lowest number of samples m equals: In practice, the number of samples is usually selected as at least 2 m. Therefore, for low-dimensional problems (especially when (n ? 1) k \ 2 m) CTPM is a better choice to construct a Chebyshev polynomial approximation; meanwhile, for high-dimensional conditions, CCM is a satisfactory alternative method.
Afterward, for a MEM domain, the Chebyshev polynomial approximation should be constructed in a corresponding IM domain at first. Due to that the MEM domain is a subset of the IM domain for the same interval variables, the Chebyshev polynomial approximation is also effective in the MEM domain. Next through a coordinate transformation method [28] (see Appendix 2 for details), take samples in the MEM domain and scan them to find the radius and middle point of the response of the nonlinear function. Meanwhile, using the sampling points, the SCCs between the response and inputs can also be obtained by (16). Thus, the approximate linear function can be calculated through (21), and the whole process of linearization is concluded as follows: Step 1 Define the order of Chebyshev polynomial n. According to the dimension of interval variables k, determine the number of interpolation points m for the construction of Chebyshev polynomial approximation.
Step 2 Take m sets of interpolation points in the IM domain of interval variables without consideration of correlation for the time being.
Step 3 Calculate values of the nonlinear function at interpolation points. Construct a Chebyshev polynomial approximation of the nonlinear function in the IM domain.
Step 4 Take samples in the ME domain through the coordinate transformation method (see Appendix 2 for details). Scan to calculate radius f R and middle point f M using the Chebyshev polynomial approximation, and calculate the SCCs q s (f I , x I ) by (16).
Step 5 Calculate the approximate linear equation by (20).

Time-varying uncertainty propagation of linear dynamics
The general expression of a linear time-varying system [36] can be expressed as: is an input vector and B(t) is the corresponding input matrix. The solution to the problem with initial states x(t 0 ) can be expressed as [36]: where U(t,t 0 ) is defined as a state transition matrix, which satisfies the following property [36]: When the inputs to the system are uncertain and described as interval processes u I (t), the responses are also transformed into interval processes [37] as: Based on the theories of differential and integral of the interval process [29,37], the middle-point functions x M (t) can be calculated as: The cross-covariance function matrix can be calculated as: where C u I u I Á ð Þ is the cross-covariance function matrix of u I (t). Finally, the radius functions can be obtained when t 1 = t 2 = t as:

Uncertainty propagation process of nonlinear dynamics
Through the CMLM, the nonlinear dynamics with time-varying interval uncertainties can be transformed into an approximate linear problem as: The middle-point function and radius function of the linear problem can be obtained through (30) and (32). However, the known middle point and covariance are required previously before the linearization through CMLM. Therefore, the problem can only be solved discretely and iteratively.
To illustrate the iterative process, it is assumed that the linearization and calculation of previous i-1 discrete time points have been finished, and N 1 * , N 2 * , …, N i-1 * are obtained as: Thus, the system can be regarded as a piecewise linear problem, and the state transition matrix, from t 0 to t i-1 , can be expressed as: Then, calculate the response of the system at the time t i , and accomplish the linearization. First, scan to calculate the middle-point function at time t i , x M (t i ). Take samples in the MEM domain of x I (t i-1 ), and solve the ordinary differential equations (ODEs) from t i-1 to t i at these sampling points, by numerical ODE methods as: where f ODE (Á) denotes an iterative step of an arbitrary numerical ODE method. x M (t i ) can be obtained by scanning method as: It should be noted that the Chebyshev polynomial approximation can also be applied here to reduce the computational cost of the scanning method. Consequently, calculate the covariance matrix function. Initialize N i , and let: The state transition matrix, from t 0 to t i , can be updated, by the property (28), as: Accordingly, the covariance matrix function can be calculated, by (31), as: It is not necessary to calculate the integration from initial time t 0 , and the results can be derived from the obtained covariance matrix at the previous time point t i-1 as: Thus, only the calculations related to the integration from t i-1 to t i are required in (41). Afterward, the linearization can be updated, by CMLM, as: The state transition matrix can also be updated by the latest N, repeat the calculation (39)-(42), until N satisfies the following criterion: Finally, the optimal approximate linear equation and the corresponding covariance matrix, at time t i , are obtained.
By applying the above iterative process, from the initial time t 0 to the final time T, the uncertainty propagation problem can be solved completely. The flowchart of the whole process is presented in Fig. 4; meanwhile, the main process is illustrated as follows: Step 1 Provide the initial conditions: x M (t 0 ), x R (t 0 ), C(t 0 , t 0 ). Define the time span and other options of the ODE numerical method. Define an initial value of the state matrix of approximate linear system N 0 * . Define the order of Chebyshev polynomial n, and according to the dimension of the dynamic problem, determine the number of interpolation points m for the construction of Chebyshev polynomial approximation. Let i = 0; Step 2 Refer to the process established in Sect. 3.1. Take m sets of interpolation points of x I (t 0 ) in the IM domain, according to the initial conditions, for the construction of the Chebyshev polynomial. Then calculate the response to these points at the next time t 1 through the numerical ODE method; Step 3 Let i = i ? 1. Construct Chebyshev polynomial approximations of the nonlinear function between x(t i-1 ) and x(t i ); then, scan to find x M (t i ). Meanwhile, according to the sampling points, evaluate an approximated radius R app (x I (t i )) as: where a is an amplification coefficient, which is in the range [2,5], to ensure the domain to construct Chebyshev polynomial approximation, for the following step, totally envelope the exact domain of x(t i ); Step 4 Refer to the process established in Sect. 3.1. Take m sets of interpolation points of x I (t i ) in the IM domain, according to the approximated domain obtained in Step 3. Then calculate the responses to the sampling points at next time t i?1 through the numerical ODE method. Meanwhile, due to the property of ODE numerical methods, the value of nonlinear functions f(Á) at x (s) (t i ) can also be obtained simultaneously; Step 5 Construct Chebyshev polynomial approximation of the nonlinear function between x(t i ) and f(x(t i )). Initialize the approximate linear system from t i-1 to t i , and let N i (0) = N i-1 * . Let j = 1; Step 6 Carry out the calculation (39)-(42), however, all the calculations of nonlinear function are replaced by the Chebyshev polynomial approximation constructed in step 5. If N i (j) satisfies the criterion (43), go to the next step, else let j = j ? 1 and repeat Step 6; Step 7 If t i \ T, go to step 3, else output x L (t) and x U (t) at every discrete time point.
The time step of the ODE numerical method should be determined according to the simulation precision of the middle-point function. In other words, the time step equals one of the simulations of the nominal trajectory.
For Step 3, there exists an alternative approach to calculating the middle-point function. For systems with small uncertainties, the middle-point function can be approximated by the nominal trajectory that can be calculated previously. In the iteration, the calculation

Duffing oscillator analysis
AT first, a single-degree-of-freedom duffing oscillator system is considered as: Under an uncertain excitation u I (t) described as an interval process, the problem can be expressed as: where v and x are, respectively, the velocity and displacement of the oscillator system; the initial condition is given as [v(t 0 ), x(t 0 )] T = [0, 0] T ; m, c and k equal 1, 0.5p and 4p 2 , respectively; and e equals 1. Moreover, the middle-point function, radius function and correlation coefficient function of u I (t) are given as: The problem is solved in the period of 0-5, and the time step of RKM is 0.025 s. The lower and upper bounds obtained by CMLM and 10,000-time-MC samples are shown in Fig. 5. Obviously, the samples are tightly enveloped by the lower and upper bounds.
For a further analysis of precision, 20 discrete time points, at 0.25 s intervals, are selected, and the value of radius function at these points yielded by the CMLM, 100-time MC, 500-time MC and 1000-time MC are compared to the ones of 10,000-time MC. The sum of squares of error can be calculated as: where N, the number of discrete time points, equals 20 in the case; and x R std represents a standard solution, which is the result of 10,000-time MC. Then the CMLM, 100-time MC, 500-time MC and 1000-time MC are, respectively, repeated 50 times. Subsequently, values of e 2 of x obtained by the above methods at all run-times are represented as a boxplot in Fig. 6. The error of the CMLM seems to be comparable to 100-time MC and be larger than the other analyses. However, when t [ 1.25 s, the error of CMLM drops significantly, as shown in (b) of Fig. 6. The error of the CMLM at the time points, which are larger than 1.25 s, is much smaller than that of 1000-time MC. Meanwhile, the CMLM presents greater robustness due to the analytical solution of covariance. In brief, the CMLM performs well in precision and robustness except for a slight error in a very short initial stage.
Finally, the precision and efficiency are concluded in Table 1. Since the CTPM is applied to construct Chebyshev polynomial approximation, for the twodimensional problem, 9 interpolation points are required. Accordingly, the solving times of the nonlinear ODEs, N ODE , is also 9. CMLM can achieve higher accuracy with much less computational cost than no less than 1000 times MC analyses.

Vehicle ride analysis
IN the second case, a two-degree-of-freedom quartercar model [18,20,23] is analyzed, as shown in Fig. 7.   Fig. 7 Schematic of a quarter-car model with two degrees of freedom and the roughness of the road where the initial condition is given as [x s (t 0 ), x u (t 0 ), v s (t 0 ), v u (t 0 )] T = [0, 0, 0, 0] T ; the sprung mass m s and unsprung mass m u equal 400 kg and 600 kg, respectively; the suspension damping rate c s equals 1000; the linear stiffness characteristics of the suspension and tire, k s and k t , equal 1.5 9 10 4 and 2 9 10 5 , respectively; and the cubic stiffness characteristics of the suspension and tire, K s and K u , equal 1.5 9 10 6 and 2 9 10 7 , respectively. It is supposed that the vehicle drives through a standard triangular road block with a speed of v = 10 m/s, x r can then be given by 0 t\0:02 0:24 À 6t 0:02 t\0:04 0 t ! 0:04 Then the standard triangular roadblock is supposed to be rough. Under the interval process model, due to an additional roughness x I (t), x r (t) is transformed into an interval process x r I (t) as where x I (t) is defined as the following interval process that Afterward, by neglecting the effect of x I (t) on the three-order term, the nonlinear dynamics with timevarying interval uncertainties problem can be defined as where f s and f u denote the nonlinear functions of v s and v u in (49).
The problem is solved in the period of 0-1 s, and the time step of RKM is 2 9 10 -3 s. The lower and upper bounds obtained by CMLM and 10,000-time MC samples are shown in Fig. 8. Similarly, 20 discrete time points, at 0.05 s intervals, are also selected. And by comparing the results of 10,000-time MC, the sum of squares of error of x s R , obtained by the CMLM, 100-time MC, 500-time MC and 1000-time MC, are, respectively, calculated. The results of a 50-time repeat are represented as a boxplot in Fig. 9. In comparison with 1000-time MC, the CMLM has a comparable but slightly larger error from 0 to 1 s as shown in (a) of Fig. 9. However, when 0.35 B t B 1, the CMLM achieves higher precision than other approaches as shown in (b) of Fig. 9.
At last, the details of precision and efficiency are listed in Table 2. The CTPM is selected to construct the Chebyshev polynomial approximation, and 81 times numerical solution of ODE is required. From the comparison, the CMLM achieves comparable precision to 1000-time with much higher efficiency.

Application in uncertainty propagation of LV ascent trajectory
The concept of the LV flight dynamic model is shown in Fig. 10, and corresponding three-degree-of-freedom dynamic equations can be expressed in vector form as (54), which describes the motion of the center of mass of an LV.
where In practical engineering, the values of the above forces cannot be obtained analytically and are always provided in the form of complex discrete tables. Therefore, the dynamic model is usually too complex to be modified arbitrarily and the calculation of LV trajectory is generally regarded as a black-box problem.

Uncertainty propagation problems
Consider a three-stage LV, and the basic parameters of the LV are provided in Table 3. Meanwhile, the parameters of the flight conditions are listed in Table 4. The flight program angle, the control program formulated to make the LV fly according to a certain flight trajectory, is given as (55). Based on the above parameters, the baseline trajectory is simulated as shown in Fig. 11.
The actual flight of the LV is usually affected by various uncertainties [38][39][40]. The most common one is the time-varying uncertainties of the atmospheric environment. At an earlier phase of LV design, precise atmosphere information is hard to be obtained. Thus, the interval process model is proper to describe the uncertainties.
where f(Á) represents the nonlinear function vector in (54) and t is the flight time. Then the uncertainty propagation problems can be established as follows.
Problem 1: During the first-stage fight, the wind will cause additional aerodynamic acceleration due to variations in dense atmosphere environments according to ?tic=?>Fig. 12. The problem can be formulated as: Fig. 11 Baseline trajectory of the LV Fig. 12 Comparison of the variation of flight height and atmospheric density Problem 2: After that, the uncertainty propagation of the following two-stage flight is considered. Since the LV has left the area with a dense atmosphere from ?tic=?>Fig. 12, there are few external time-varying uncertainties. The flight is mainly impacted by uncertain initial conditions left by the first stage. Thus, the uncertainty propagation of the last two stages can be expressed as: In conclusion, the uncertainty propagation of LV trajectory is defined, which is divided into two subproblems: One is a problem with time-varying uncertainties of the first stage formulated as (57), and the other is a problem with correlated initial uncertainties of the second and third stages formulated as (59).

Results and discussions
Firstly, Problem 1 is solved through the CMLM, in the period of 0-60 s. The RKM is applied as the solver, and the time step is 0.5 s. The order of Chebyshev polynomials for the approximation is defined as 2. The lower and upper bounds obtained by CMLM in comparison with 10,000-time-MC samples are shown in Fig. 13.
Then 12 discrete points, at 5-s intervals, are selected from 0 to 60 s. At these points, the e 2 of v z and z obtained by the CMLM,   is shown that the CMLM performs better in robustness, meanwhile providing a result with comparable precision to 1000-time MC. Finally, the precision and efficiency are concluded in Table 6. For a six-dimensional problem, the CCM is chosen to construct the Chebyshev polynomial approximation. Therefore, at least 56 interpolation points are required, while the minimum times of numerical solution of ODE are 56. It is shown that CMLM achieves higher precision than 1000-time MC in v y , v z , y and z with a considerably smaller computational cost. Thus, it can be concluded that the CMLM performs well in precision, efficiency and robustness for a practical LV trajectory problem under time-varying uncertainties.
After the uncertainty propagation of the first-stage flight, a correlated uncertain domain of the final state is caused. The values of the middle point and radius are listed in Table 6; meanwhile, the correlation coefficient matrix is expressed as follows: Under the above initial value uncertainties, Problem 2 can also be solved by CMLM, and the results are shown in Fig. 15. The problem is solved in the period of 60 s-185 s, and the RKM is still employed as the solver. The order of Chebyshev polynomials for the approximation is also defined as 2. The lower and upper bounds of v z and z obtained by CMLM in comparison with 10,000-time-MC samples are shown in Fig. 15. Then 25 discrete points, at 5-s intervals, are selected from 60 to 185 s. At these points, the e 2 of v z and z to 10,000-time MC obtained by the CMLM,  Fig. 16, and the precision and efficiency are concluded in Table 7. It is shown that the CMLM performance a better in precision and robustness than 1000-time MC with only 56-time calculation of the original functions.
Eventually, through solving Problem 1 and Problem 2, the uncertainty propagation of the whole trajectory is obtained, as shown in Fig. 17. The results demonstrate the effectiveness of the proposed CMLM in solving practical engineering nonlinear dynamics with time-varying uncertainties and correlated initial uncertainties under the interval process model.

Conclusion
This work applies the interval process model to quantify the uncertain excitations and responses of nonlinear dynamics. The interval process works well in describing the bounds and correlations of timevarying uncertainties without precise probability distributions required. To efficiently propagate uncertainties of nonlinear dynamics under the interval process, a novel method, named CMLM, is proposed in this work.
The CMLM realizes uncertainty propagation by linearizing nonlinear systems non-intrusively and utilizing analytical equations to calculate the bounds of system responses. During the linearization, representing the differences between two intervals by middle point, radius and correlation takes effect in constructing optimal linear approximations of nonlinear dynamics under interval uncertainties. Moreover, the Chebyshev polynomial effectively reduces the computational cost of nonlinear system analyses while providing the CMLM with satisfactory precision. The performances of the CMLM are verified by test cases involving two numerical cases and a practical engineering application of LV ascent trajectory. The results show that the CMLM achieves comparable precision and greater robustness than 1000-time MC simulation with only no more than 10% of the computational cost for the MC simulation required. It is also demonstrated that the CMLM is effective in solving complicated black-box problems in practical engineering.
Funding The present work was supported by the major advanced research project of Civil Aerospace from State Administration of Science, Technology and Industry of China.

Derivation of CMLM criterion
By the weighted sum of e 1 , e 2 and e 3 , a comprehensive index e can be established. The weighting coefficients are first determined. The index e 1 is independent because it is only related to M and not to N. Thus, its weighting coefficient can be defined as 1, and the value does not affect the optimal value of M and N. Then, to achieve the weighted sum of e 2 and e 3 in an equivalent order of magnitude, they are compared in a similar form represented by covariances: From the comparison, it is obvious that e 3 should be weighted by Then e 2 and e 3 can be weighted equivalently, and e is finally expressed as: where the SCC q s is defined as: where a and b are the angle of attack and sideslip angle, respectively. C x is the drag coefficient, and C y a is the derivative of the lift coefficient with respect to a; S R is the reference surface area; and q is the dynamic pressure, which is calculated as: where q is the atmospheric density and v is the resultant velocity of the LV flight as Then, G B and G V are coordinate transform matrixes as: cos ucosw À sin u cos u sin w sin u cos w cos u sin u sin w À sin w 0 cos w where u and w are, respectively, the pitch angle and yaw angle, which describe the flight attitude of the LV.
Meanwhile, h and r are, respectively, flight path angle and flight path azimuth angle, which describe the flight direction of the LV. These angles can be derived as: Moreover, A and B, in(75), are the matrixes to describe inertial force caused by the rotation of the earth as: x ex x ey x ex x ez x ex x ey x 2 ey À x 2 e x ey x ez x ex x ez x ey x ez where x e is the earth rotation rate and [x ex, x ey, x ez ] T are the components of the vector x e . Afterward, [R 0x , R 0y , R 0z ] T , in (75), are the components of the vector R 0 to describe the position of the launch point. Next, g r and g xe are the components of gravitational acceleration and can be calculated as: where l and J are the constant characteristics of gravity. a e is the length of the semimajor axis of the earth under an ellipsoid model; besides, the semiminor axis is symbolled by b e . And, r, the geocentric distance of the LV, is calculated as: Meanwhile, /, the geocentric latitudinal, can be derived from: In addition, the flight height of the LV can also be obtained by r and / as: h ¼ r À a e b e ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a 2 e sin 2 / þ b 2 e cos 2 / q ð85Þ Finally, these equations can be solved according to a given flight program angle, and generally, they are provided in the form: Accordingly, to achieve the flight program, the corresponding a and b are expressed as a ¼ A u ½ðu PR À x ez t À hÞ b ¼ A w ½ðu ex sin u À x ey cos uÞt À r ( ð87Þ where A u and A w are both constant coefficients.