A bivariate Chebyshev polynomials method for nonlinear dynamic systems with interval uncertainties

A bivariate Chebyshev polynomials approach is proposed to estimate the dynamic response bounds of nonlinear systems with interval uncertainties. The existing collocation method directly searches the maximum and minimum values of the surrogate model in the entire interval space by the scanning method (SM). The presence of too many uncertain parameters will lead to expansive computational cost. To overcome this shortcoming, the dynamic response is decomposed by a bivariate function decomposition (BFD), established based on high-order Taylor expansion, into the sum of multiple univariate and bivariate response functions. The above univariate and bivariate functions are fitted using Chebyshev polynomials, and polynomial coefficients are obtained through one-dimensional (1D) and two-dimensional (2D) interpolation points. Thus, the solution of the nonlinear dynamic systems with uncertain parameters can be transformed into that of univariate and bivariate Chebyshev interval functions. The extremum values of the low-dimensional Chebyshev interval functions can be found by SM, and then the bounds of dynamic response are acquired by interval arithmetic. Since SM searches for extreme values only in 1D and 2D uncertain domains, the amount of calculation is reduced compared to searching the whole uncertain space. The efficiency, practicability and effectiveness of the proposed interval uncertainty analysis method are proved by three dynamic examples.


Introduction
Dynamic response forecasting of nonlinear systems is of great significance for structural design and analysis of most practical engineering. Traditional dynamic analysis is based on the deterministic model. However, any dynamic system has inevitable uncertainties, such as load fluctuation, manufacturing instability, and dispersion of material properties, which make structural parameters uncertain [1][2][3][4]. Hence, the dynamic response of nonlinear systems is also uncertain. Noticeably, small changes associated with parameter uncertainty may largely change the dynamic response [5]. As a result, the influence of uncertainty on dynamic systems must be considered.
In uncertainty analysis, probabilistic and non-probabilistic methods are often employed to quantify the effect of uncertainty on dynamic response. In general, the probability method is the most valuable solution under the premise that the probability density function (PDF) of uncertain parameters is clearly defined [6][7][8]. Unfortunately, the method has an obvious disadvantage, as it needs complete statistical data to accurately describe the PDF of random variables or fields, and these data are usually inaccessible and expensive in engineering problems. Under the circumstances, some non-probabilistic methods, including fuzzy set [9], evidence theory [10] and interval set [11], are developed to describe the uncertainty caused by limited data and knowledge.
In interval sets, uncertainty is represented by the upper and lower limits of uncertain parameters, rather than other stochastic data. Therefore, it is an effective assistance to the conventional probability method. In the past two decades, interval sets and interval analysis approaches have been widely used in structural dynamic response investigation. The interval arithmetic (IA), which was the first applied for dynamic T. Wei F. Li () G. Meng School of Mechanical and Aerospace Engineering, Jilin University, 5988 Renmin street Changchun 130025, China E-mail: Fengli@jlu.edu.cn response analysis, directly obtains the response interval based on interval algebraic operation rules [12][13][14][15][16]. Since the correlation between uncertain parameters is neglected, IA usually leads to overestimation of interval values, which is called the wrapping effect. This effect is more evident and serious in multidimensional and non-monotonic dynamic systems. There is extensive research on this topic [17,18]. Besides, the affine algorithm [19][20][21] introduced to remedy the wrapping effect can alleviate overestimation, but can cause serious errors in the case of highly nonlinearity. This is because nonlinear affine operation will induce more noise terms and overestimate response bounds [21]. As long as the uncertainty problem is monotone, the vertex method can get accurate interval results [22][23][24]. However, it can hardly guarantee or prove monotonicity, especially in the field of structural dynamics.
The Taylor expansion based method (TEBM) performs the first-or second-order Taylor series expansions at the midpoint of uncertain parameters and obtains the dynamic response interval [2,[25][26][27][28][29]. The advantage of TEBM is high computational efficiency, and the calculation consumption is raised linearly as the number of interval parameters rises. Meanwhile, it can reduce the wrapping effect due to the dependence among multiple interval parameters. However, TEBM has some disadvantages. First, it is accurate only when the range of uncertain parameters is narrow enough. The subinterval approach is a solution to this problem [30][31][32]. It divides the large interval into several small intervals, which improves the accuracy of interval analysis. But the amount of calculation will sharply increase as the number of subintervals increases. Second, TEBM requires the first-or even higher-order partial derivatives of the dynamic response concerning the interval variables. In some engineering conditions, the explicit functional relationship between input and output is too complicated, which may make it impossible to acquire first-and higherorder derivatives. Most importantly, TEBM can only moderately relieve overestimation, but may even lead to worse results in case of strong nonlinearity [26,33].
The dynamic response is always highly nonlinear relative to interval parameters. Compared with the above methods, the collocation method (CM) can solve large uncertainty and strongly nonlinear problems. Firstly, CM establishes a surrogate model according to the selected interpolation points. The common proxy models include Chebyshev polynomials [34][35][36], Legendre polynomials [37][38][39], and arbitrary polynomial chaos [40,41,11]. Then, the auxiliary algorithms are applied to find the extreme values of the surrogate model, and obtain the interval results of the dynamic response. The common auxiliary algorithms include IA [35,42,43], genetic algorithm (GA) [40,44,45] and SM [46][47][48]39]. IA has the fastest computing velocity, but often yields calculation errors due to the wrapping effect, especially when the dynamic response is largely uncertain and nonlinear. To this end, many improved IA-based methods have been studied, such as trigonometric function [49] and modified affine basis [50]. These methods only moderately improve the accuracy, but large errors still exist with the precise dynamic interval. Generally, GA achieves the predetermined accuracy through iteration, and its computation is expensive if there are many uncertain parameters. Moreover, when the nonlinearity of the surrogate model is strong, GA may fall into a local optimum. Compared with IA and GA, SM is the most precise, but its samples rise exponentially with the increased number of interval parameters.
This study concentrates on the following issues: When SM is employed to compute the maximum and minimum values of the surrogate model, the existing collocation methods need to search the entire uncertain domain. Thus, when there are many interval parameters, it easily falls into the curse of dimensionality. To increase efficiency and ensure accuracy, this paper proposes a bivariate Chebyshev polynomials method for interval analysis of dynamic response with many uncertain parameters. In this method, BFD decomposes the original dynamic function into the sum of univariate and bivariate response functions, which are then fitted by Chebyshev polynomials. In this way, the nonlinear dynamic systems described by ordinary differential equation (ODEs) can be transformed into 1D and 2D ODEs, and the whole interval space is divided into a linear combination of 1D and 2D spaces. Moreover, the ODEs can be solved by an appropriate numerical solver. Finally, the upper and lower bounds of the dynamic response at each time step are determined using SM and IA.
The statement of nonlinear dynamic system with uncertain parameters is given in Sect. 2. Section 3 illustrates the detailed derivation of Chebyshev interval method based on IA or SM. Then, a bivariate Chebyshev polynomials method is proposed in Sect. 4. In Sect. 5, the three dynamic cases are presented to verify the effectiveness and practicability of the proposed method.

Overview of nonlinear dynamic system with interval parameters
In general, the nonlinear dynamic system is represented by a group of ODEs. In the case of l state vectors, the l-dimensional ODEs are demonstrated as: is a vector containing nonlinear functions,   T 12 , , , l y y y = y is a dynamic response state vector, and y0 is the initial condition vector; is the time range. Equation (1) is a deterministic model that can be computed by any numerical integration approach to obtain the change of y over time.
In actual engineering, the physical and geometric properties of the structure are uncertain to varying degrees. For uncertain input vector x, if the upper and lower bounds of the parameters can be estimated, then x can be expressed as: (2) where the superscript I is the interval; C k x and k x  represent the midpoint and radius of the kth uncertain parameter respectively; the transformation parameter I k  is a fixed interval defined as: where U k x and L k x represent the upper and lower bounds of k x , respectively. Since the parameter uncertainty is transitive, the original ODEs become a group of interval ODEs. Therefore, the solution of Eq.
(1) subjected to Eq. (2) is also an interval set: , , , Generally, the set  has a complex region. In interval mathematics [12], solving the nonlinear problem of Eq. (1) constrained by Eq. (2) is synonymous to find a multidimensional rectangle defined by Eq. (4). In other words, the upper and lower bounds of Eq. (4) will be calculated as follows: Obviously, the interval value of the dynamic response is essentially the supremum and infimum of Eqs. (6) and (7). The simplest and most direct method to solve Eq. (5) is SM [51], but it is computationally intensive. Additionally, IA and TEBM are simple, effective and commonly applied to evaluate interval functions. However, they may lead to overestimation as mentioned in the introduction. Moreover, TEBM is only suitable for small uncertainty problems. To solve the dynamic response problem with large uncertainty 12 [ , , , ] n X X X = X are transformed as follows: Let p T be the Chebyshev polynomial of degree p, then: For n-dimensional problems,  y X X (14) Given the difficulty to directly solve the integral formula of Eq. (14), the Gauss quadrature method [52] is introduced to calculate  y (15) where k j  is the interpolation point. The zeros of the Chebyshev polynomial of degree (p+1) are the interpolation points, which can be computed as: The number of interpolation points M for an n-dimensional problem is Obviously, M will be large for high-dimensional problems, leading to expensive calculation. To save costs, only a part of the interpolation points are utilized to build the Chebyshev polynomial. Therefore, the Chebyshev surrogate model in Eq. (8) (17) In Eq. (17) (19) After that, the bounds of ( ) I IC yx are computed. Two methods are commonly used: (1) improved Chebyshev interval model using interval arithmetic (IC+IA); (2) improved Chebyshev interval model using scanning method (IC+SM).  (20) 6

Improved Chebyshev interval model using interval arithmetic
By expanding Eq. (20), the lower bounds of ( ) I IC_IA yx can be stated as: (21) The upper bounds of ( ) I IC_IA yx can be defined as: (22) Equations (21) and (22) can reduce the overestimation to a certain extent, but cannot completely eliminate the overestimation [35].

Improved Chebyshev interval model using scanning method
To reduce the error introduced in solving the bounds of II 11 cos cos nn ii  using IA, the interval value of ( ) First, all the sampling points enclosed in the interval I x are selected. In each interval be a sampling point, the set of all sampling points x is formulated as: (23) In Eq. (23), there are a total of e n sampling points. Second, each given sampling point is substituted into Eq. (17) Similarly, the upper bounds of ( ) x y x (26) Equations (25) and (26) are applied to scan the explicit IC constructed by Eq. (19), which improves the computational efficiency compared with the direct scanning of the original dynamic response function. However, considering that the sampling points of IC+SM is e n , the calculation cost is still unacceptable for the dynamic response problem with more sampling locations or more interval parameters.

Construction of Bivariate Chebyshev interval model
To reduce the calculating consumption of scanning the IC, we first decompose the n-dimensional dynamic response function into a sum of univariate and bivariate response functions by BFD, and then fit these functions based on Chebyshev polynomials.

Bivariate function decomposition
The original dynamic response function ( ) yx is processed by high-order Taylor expansion at the midpoint C x of the uncertain parameters, and then retained to the second-order term [53]: First, let Similarly, these univariate points Similarly, these bivariate points Second, adding each term in Eqs. (28) and (29), the sum of univariate function is expressed as: Similarly, the sum of the bivariate function can be obtained by adding Eqs. (30) and (31):  (33) By shifting the terms of Eqs. (32) and (33) separately, we can get: Finally, substituting Eqs. (34) and (35)

Bivariate Chebyshev interval model
The above ( )   (41) The interpolation points required by the Chebyshev model (CM, Eq. (8)), IC, and BC are given in Table  1. The number of interpolation points of CM increases exponentially with the rising number of uncertain parameters, but those of IC and BC are dramatically reduced, avoiding the huge calculation costs (Table 1). Furthermore, BC has fewer interpolation points for high-dimensional problems.

Bivariate Chebyshev interval model using auxiliary methods
This section first uses the SM to calculate the limits of ( ) be a 2D sampling point, then the set of all 2D sampling points is:  (43) Therefore, each bivariate function requires e 2 sampling points. The set of all sampling points BĈ x can be stated as: (44) Obviously, scanning BC requires n(n-1)e 2 /2+en sample points in total.
Second, ˆk x is substituted into Eq. (38) to obtain the value of ( )  (46) y y x y y x (47) Similarly, the upper bounds of ( ) y y x y y x (48) Finally, substituting Eqs. (47) and (48)  Compared with IC+SM in Eqs. (25) and (26), which requires e n samples, BCM only needs n(n-1)e 2 /2+en. The sample distribution comparison between IC+SM and BCM in a three-dimensional uncertain space is shown in Figure 1. Clearly, BCM has a relatively low computational cost (Fig. 1).

BCM to solve nonlinear dynamic systems with interval parameters
BCM is introduced for nonlinear dynamic systems expressed by ODE with interval uncertainty. When employing numerical methods to solve differential equations, the time range 0 , f t t t    is discretized into a In Eq. (51), the original dynamic response is approximated as the sum of multiple univariate and bivariate functions. Then, these functions are explicitly fitted by Chebyshev polynomials, and the expansion coefficients are obtained through solving 1D and 2D ODEs. When the surrogate model is built, the interval of y at time b t can be determined simply by auxiliary methods. Moreover, any numerical method can be used to solve 1D and 2D ODEs. The flowchart of the solution is depicted in Fig. 2: Step (1). Define the degree p of the Chebyshev polynomial, the interval of x, the sampling position e of each interval parameter, the length of the time period tf, and the time step Δt.
Step (2). Decompose the original dynamic response into the sum of multiple univariate and bivariate functions by BDE (Eq. (36)); Step (3). Generate 1D and 2D interpolation points by Eq. (16) for fitting the above univariate and bivariate functions. Initialize the loop count index 0 b = , which is used to calculate the iterations of time; Step (4). Substitute the 1D and 2D interpolation points separately into the original differential equation, and use the numerical method to compute the actual 1D and 2D dynamic responses at the time Step (5). Calculate the 1D and 2D expansion coefficients, Step (6). Insert the 1D and 2D sampling points into the univariate and bivariate Chebyshev interval functions, and obtain the sets Γ k and 12 Γ kk ; Step (7). Scan to find the maximum and minimum values of Γ k and 12 Γ kk , and get the bounds of univariate and bivariate Chebyshev interval functions; Step (8). Use IA to compute the interval value of y at the time

Numerical cases
Three numerical cases are used to verify the effectiveness and accuracy of the proposed method: (1) a double pendulum problem with three interval parameters, (2) a damped spring-mass system with five interval parameters, and (3) a two-degree-of-freedom (DOF) vehicle model with seven interval parameters. SM (regarded as the exact solution), IC+IA, IC+SM, and BCM are used for comparison. All the algorithms are programmed using MATLAB R2019a operated on the same computer. The numerical calculation adopts the Runge-Kutta method.

Double pendulum problem
In this double pendulum system (Fig. 3), 1  and 2  are the angles of the pendulum, 1 m and 2 m are the masses of the two rods, 1 L and 2 L are the lengths of the two rods, and g is the local acceleration of gravity. The ODEs of the highly nonlinear double pendulum system are expressed as [23]: where 1 w and 2 w represent the angular velocity of the two pendulum rods respectively. In this example, 1 L ,  Tables 2 and 3. IC+SM and BCM provide more accurate results than IC+IA, and IC+SM and BCM almost consistent with SM (Figures 4 and 5, Tables 2 and 3). Moreover, the results of IC+IA diverge slowly with time. In contrast, IC+SM and BCM can effectively reduce the wrapping effect throughout the numerical process.
Compared with SM, IC+SM and BCM are both searching for the maximum and minimum values of the surrogate models, and thus take less computational consumption. IC+SM requires 20 3 =8000 sampling points at each time point for the double pendulum problem with 3 interval parameters. However, BCM only needs 1261 sample points, which is more efficient than IC+SM.

Damped spring-mass system
In this two-DOF damped spring-mass model (Fig. 6), 1 M and 2 M are masses, C is the damping coefficient, and 1 K and 2 K are stiffness of the springs. Besides, the external forces F=1 N acts on the mass 2 M . The ODEs of the system can be expressed as [54]:  (Table 4).  The y1 obtained by IC+IA, IC+SM, and BCM has a smaller error in the linear region (i.e. in the red dashed box) (Fig. 7), which can also be illustrated by y2 in the red dashed box of Fig. 8. However, in the highly nonlinear region (the peaks and valleys of the curve), the error between the dynamic response intervals obtained by IC+IA and SM is considerably large. The displacement intervals of IC+SM and BCM are always close to the accurate solution without major deviation (Figs. 7 and 8, Tables 3 and 4). Therefore, IC+SM and BCM can acquire correct and dependable dynamic response interval throughout the numerical period.
For the damped spring-mass system with five interval parameters, IC+SM requires 10 5 sample points at each time point, while BCM only needs 1051 sample points, so the computational efficiency is much higher than IC+SM.

Two-DOF vehicle model
To prove that the proposed method is applicable to practical engineering, the third example adopts a two-DOF vehicle model ignoring the influence of steering and suspension system (Fig. 9). In Fig. 9, O is the center of mass; v and u are the lateral and longitudinal velocities of the center of mass respectively;  is the sideslip angle of the center of mass; r w is the yaw velocity; a and b are the distances from the center of mass to the front and rear axles respectively;  is the front wheel steering angle.
The ODEs of the two-DOF vehicle can be written as [55]: where k1 and k2 are the front and rear wheels cornering stiffness respectively; m is the vehicle mass, z I is the moment of inertia.
 can be denoted as: This example considers the dynamic response of the steering angle step input: the step input is defined as the front tire rotating to 2° in 0.2 s, and then keeping the angle unchanged, which can be expressed as: Assuming that all parameters of the vehicle model are uncertain, the information of uncertain parameters is listed in Table 7. The midpoint values of these parameters are all from Ref. [55]. Moreover, the initial condition is Clearly, this dynamic model is relatively complex. SM, IC+IA, IC+SM and BCM are all applied to solve this dynamic systems in a time period of 0-3.5 s. The time step is 0.001 s t = . Given the degree p of the Chebyshev polynomial is 2, each interval parameter uses 10 symmetrical sampling positions. The interval bounds of r w are shown in Fig. 10. Besides, ten time points are chosen ( Table 8). The interval bounds of r w generated by IC+SM and BCM are basically consistent with that of SM, while IC+IA has a large error (Fig. 10). Table 8 and the red dotted box in Fig. 10 demonstrate this result more intuitively. Therefore, BCM has a good accuracy even for complex dynamic systems. However, due to the use of IA in BCM, some overestimation cannot be completely avoided.
Given that IC+SM needs 10 7 sample points at each time point for the vehicle model with 7 interval parameters and that BCM only requires 2171 samples, BCM has a very high computational efficiency. Furthermore, the example clearly shows that BCM can achieve high accuracy and low calculation amount when computing the response limits of a nonlinear dynamic system with interval uncertainty.

Conclusions
A bivariate Chebyshev polynomials method is put forward to estimate the upper and lower bounds of nonlinear dynamic systems with interval parameters. This method can effectively alleviate the curse of dimensionality caused by scanning the entire uncertain domain. Based on BFD and Chebyshev polynomials, the solving of high-dimensional ODEs with interval parameters is converted into that of 1D and 2D ODEs, and the entire uncertain space is decomposed into a linear combination of 1D and 2D spaces. The upper and lower bounds of the dynamic response at each time step can be easily approximated by SM and IA. Three numerical examples are used to illustrate the feasibility and effectiveness of the method, and the following conclusions can be further summarized: BCM provide more precise and reliable dynamic response limits than IC-IA. Meantime, the interval value of BCM is almost consistent with those of IC-SM and SM in the whole numerical period, but its computational efficiency is much higher. Hence, BCM is more accurate with less computational cost in solving nonlinear dynamic systems with interval parameters. BCM is also non-intrusive and only requires 1D and 2D interpolation points to obtain a surrogate model of the original response. Thus, BCM is more easy-to-use in practical engineering, particularly for implicit dynamic response with multiple uncertain parameters.

Interval arithmetic
Let LU , aa   and LU , bb   be two intervals, then interval arithmetic can be written as: a b a b a b a b  a b a b a b