An Uncertainty Analysis Method for Artillery Dynamics with Hybrid Stochastic and Interval Parameters

Liqun Wang,Zengtao Chen and Guolai Yang,＊

1School of Mechanical Engineering,Nanjing University of Science&Technology,Nanjing,210094,China

2Department of Mechanical Engineering,University of Alberta,Edmonton,T6G 1H9,Canada

ABSTRACT This paper proposes a non-intrusive uncertainty analysis method for artillery dynamics involving hybrid uncertainty using polynomial chaos expansion (PCE).The uncertainty parameters with sufficient information are regarded as stochastic variables,whereas the interval variables are used to treat the uncertainty parameters with limited stochastic knowledge.In this method,the PCE model is constructed through the Galerkin projection method,in which the sparse grid strategy is used to generate the integral points and the corresponding integral weights.Through the sampling in PCE,the original dynamic systems with hybrid stochastic and interval parameters can be transformed into deterministic dynamic systems,without changing their expressions.The yielded PCE model is utilized as a computationally efficient,surrogate model,and the supremum and infimum of the dynamic responses over all time iteration steps can be easily approximated through Monte Carlo simulation and percentile difference.A numerical example and an artillery exterior ballistic dynamics model are used to illustrate the feasibility and efficiency of this approach.The numerical results indicate that the dynamic response bounds obtained by the PCE approach almost match the results of the direct Monte Carlo simulation,but the computational efficiency of the PCE approach is much higher than direct Monte Carlo simulation.Moreover,the proposed method also exhibits fine precision even in high-dimensional uncertainty analysis problems.

KEYWORDS Uncertainty propagation and analysis;artillery dynamics;hybrid stochastic and interval uncertainty;polynomial chaos expansion;ordinary differential equations

1 Introduction

In practical dynamic systems,a variety of uncertainties associated with material properties,environmental factors,external loads,dimensional tolerances,and boundary conditions are ubiquitous.These uncertainties will inevitably affect the final system performances,and small variations associated with uncertainties might result in significant changes in the dynamic responses.The traditional way to investigate the dynamic response problems with uncertainty parameters is the probabilistic method,in which the uncertainty parameters are described as stochastic variables or processes with precise probability distributions.Among these probabilistic methods,Monte Carlo simulation (MCS) [1-3] is considered to be the most accurate and most versatile.Despite the simplicity,it suffers from the considerable computational burden,thus it is usually used to generate the reference solutions.To circumvent this drawback,the stochastic perturbation method [4-6]using first-order or high-order Taylor series expansion was incorporated into MCS to investigate the dynamic response problems.However,it is precise only for the small uncertainty level problems.Wan et al.[7] and Lu et al.[8] presented a general framework for analytical uncertainty quantification of model outputs using a Gaussian process metamodel,where the inputs are characterized as normal and/or uniform random variables,and characterized the uncertainty of modal frequencies of bridge systems.Lu et al.[9] employed modal decomposition together with a lowdimensional representation method to address the curse of dimensionality of frequency response functions output,subsequently provided a fast vector-output approximation of the random modal parameters through multi-output Gaussian process metamodel.Moreover,the spectral stochastic method [10-12] is also an efficient alternative for stochastic problems,of which the polynomial chaos expansion (PCE) method [13-15] has drawn much attention due to its strong power in functional representations of stochastic parameters.PCE has achieved significant successes in dealing with structural dynamic responses.Xiu et al.[13] presented the method for solving stochastic differential equations based on the Galerkin projections and extensions of Wiener’s polynomial chaos.Sandu et al.[16] Saha et al.[17] applied the generalized polynomial chaos to model the complex,nonlinear,multibody dynamic systems with stochastic uncertainty.Saha et al.[18]applied the generalized PCE to stochastic,dynamic response analysis of the base-isolated liquid storage tanks.Wan et al.[19] presented an efficient approach based on arbitrary polynomial chaos expansion for analytical,unified implementation of uncertainty quantification and global sensitivity analysis in structural dynamics.For the aforementioned probabilistic methods,determining precise probability distribution is required,which usually needs complete stochastic information or experimental data.Unfortunately,the objective information is usually insufficient and expensive to obtain in practical engineering problems.In this context,multiple types of non-probabilistic models,such as fuzzy model and interval model,have been used to quantify the uncertainty.

Interval methods depict uncertainty through the upper and lower bounds of uncertain-butbounded parameters rather than other stochastic information,therefore,they are considered to be a powerful supplement to the classical probabilistic approaches and have been perfectly applied in dynamic structural response analysis.The interval arithmetic operation [20,21] was the first used for dynamic response analysis.We can obtain the approximate response bounds effectively through the specific interval operators.However,its inherent deficiency is the large overestimation of range enclosure namely wrapping effect [22,23],which becomes more obvious and severe in the multidimensional and non-monotonic problems.To this end,some other interval-based,dynamic response analysis methods,such as the Taylor model method [24,25],interval vertex method [26],and the heuristic optimization method [27],have been developed.Based on the first-order or second-order Taylor series expansion,Qiu et al.[28,29],Zhang et al.[30],and Han et al.[31] estimated the ranges of the nonlinear structural dynamic responses.However,it should be mentioned that the Taylor model methods are precise only if the uncertainty levels of interval parameters are small enough.To overcome this drawback,the subinterval method [32] was introduced as an important remedy to obtain more accurate results.Xia et al.[26] presented a method to produce the tighter dynamic response bounds,based on the vertex solution theorem for the first-order deviation of dynamic response.Liu et al.[27] employed a particle swarm optimization method to search for the dynamic response bounds of a vehicle-bridge interaction system.Wu et al.[33] proposed a Chebyshev polynomial series method to control the large overestimation.Wan et al.[34] proposed an efficient Bayesian optimization approach for estimating the dynamic response bounds of a train-bridge system.The biggest merit of the aforementioned interval methods is that it can produce a rigorous enclosure that contains all possible solutions by using limited information.However,the relevant research results [28,35] showed that the results of interval methods may be conservative.

As aforementioned,the probabilistic methods and interval methods both have their own merits,deficiencies,and limited scope of application.In practical structures,the stochastic and interval parameters may exist simultaneously due to the different availability of uncertainty information.In the early researches [36-38],only the external excitation is considered as a Gaussian random process,while the structure parameters are regarded as interval variables.These works basically consisted of combining random vibration theory with the first-order interval Taylor expansion of the mean-value and covariance vectors of the response,which can evaluate the lower and upper bounds of the second-order statistics of the response.Recently,the structural response analysis with respect to the hybrid stochastic and interval structure variables has attracted more attention,and some uncertainty techniques have been developed to treat such situations.In the relevant literatures [39-43],the statistics of the system response were derived based on the perturbation technique and random moment method,and the interval bounds of statistical moments of the system response were further calculated through the interval vertex method [39,40] and the interval perturbation method [41-43].Moreover,Xia et al.[44] proposed an inverse mapping hybrid perturbation method based on the Taylor series expansion,in which the functional representations of the response probability were obtained through the change-of-variable technique.However,it should be mentioned that these methods have some deficiencies.The stochastic and interval perturbation methods based on the Taylor series expansion are only suitable for problems of small uncertainty levels and systems of weakly nonlinearity.Besides,the applications of the vertex method are limited to monotonous cases,which is not suitable for the complex dynamics problems.In another type of methods,Wang et al.[45],and Xu et al.[46] proposed the polynomial chaos expansion methods for the prediction of frequency response of the structural-acoustic system with stochastic and interval parameters,whereby evaluating the interval bounds of statistical moments effectively.Wu et al.[47] integrated the PCE theory that accounts for the random uncertainty with the Chebyshev inclusion function theory that describes the interval uncertainty,and delivered a Polynomial-Chaos-Chebyshev-Interval (PCCI) method for the uncertainty analysis of vehicle dynamics.Wang et al.[48] proposed a non-intrusive uncertainty analysis method named polynomial chaos response surface (PCRS) for the frequency response analysis of acoustic field with hybrid uncertainty.In their work,the PCE model is employed to deal with the stochastic parameters,and the response surface method is used to handle the interval parameters.

These aforementioned methods are very heuristic,and have laid a solid methodological foundation for uncertainty analysis of artillery dynamics.However,it is frankly to say that few literatures considered uncertainties in artillery dynamics.Among them,Wang et al.[49] proposed an interval uncertain optimization method for structural dynamic responses of artillery system.From the overall perspective,research on the uncertainty analysis of artillery dynamics is in its preliminary stage,the potential and applicability of these methods in dealing with artillery dynamics system is required to explore.Developing a generally applicable method for estimation of the response bounds of artillery dynamics with hybrid uncertainty parameters is necessary.Polynomial chaos has been used extensively to model uncertainties in structural mechanics and in fluids,but to our knowledge they have yet to be applied to artillery dynamic simulations.The application of the polynomial chaos expansion method in this domain is promising,yet mostly unexplored.

The dynamics systems described by ordinary differential equations (ODEs) which are quite common in artillery dynamics are considered in the present paper.To overcome the potential limitations for current research,the focus of this paper is to develop an accurate,efficient and generally applicable method for estimation of the response bounds of artillery dynamics with hybrid uncertainty parameters.Thus,a non-intrusive polynomial chaos expansion framework will be adopted due to its fine capacity and strong mathematical basis in quantifying the uncertainty of the dynamic systems.The PCE model is constructed through the Galerkin projection method,in which the sparse grid numerical integration method is used to generate the integral points and the corresponding integral weights.Through the sampling in PCE model,the original dynamic systems with hybrid stochastic and interval parameters can be transformed into ones with deterministic parameters,without changing their expressions.The yielded PCE model is utilized as a computationally efficient,surrogate model,and the supremum and infimum of the dynamic responses at each iteration time step can be easily approximated through MCS and percentile difference.

The remainder of this paper is organized as follows.A statement of the dynamic systems with hybrid uncertainty parameters is given in Section 2.Then in Sections 3 and 4,the proposed method and its solving process are introduced in detail,respectively.In Section 5,an artillery exterior ballistic dynamics model is given to show the effectiveness and feasibility of the proposed method in comparison with MCS.Finally,the conclusions are drawn in Section 6.

2 Problem Statement

The dynamics systems expressed as a set of ODEs are quite common in artillery dynamics.Suppose that there areqstate variables,theq-dimensional ODEs can be described as follows:

wherey=[y1,y2,...,yq]Tis aq-dimensional vector of dynamic responses,f=[f1,f2,...,fq]Tis aq-dimensional vector comprising nonlinear functions,y0is aq-dimensional initial condition vector,x=[x1,x2,...,xn]Tis ann-dimensional vector,t?[t0,te] is the time range.Eq.(1) represents a deterministic model.If no uncertainty exists,Eq.(1) can be solved by any numerical integration method,such as the Runge-Kutta method,to obtain the variation ofywith time.

Assuming that there are uncertainties in the structure,of which the uncertainty parameters with precise probability distributions are regarded as stochastic variables,given by

where superscriptRdenotes stochastic variables,D(xi) is the distribution function of thei-th stochastic parameter.Furthermore,the uncertainty parameters with limited stochastic information and knowledge are constrained to lie within an interval,given by

where superscriptIdenotes interval parameters,anddenote the upper and the lower bounds of thei-th interval parameter.

A dynamic system with hybrid stochastic and interval uncertainty parameters can be rewritten as

Considering the transitivity of uncertainty parameters,the solution of Eq.(4) subject to Eqs.(2) and (3) is also uncertain rather than deterministic values.The lower bound and upper bound ofycan be obtained by

In Eq.(4),the dynamic uncertainty propagation problem is to solve the ODEs with hybrid stochastic and interval parameters.The dynamic responses bounds are essentially the extrema at each integration time of the functional equations.Generally,the precise solution is difficult to compute directly through Eqs.(5) and (6).The most robust and accurate method may be MCS.However,it is well known that MCS is computationally intensive.To overcome this drawback,the PCE model is applied in the estimation of the dynamic responses bounds,which will be described in the following sections.

3 Polynomial Chaos Expansion for Uncertainty Functions with Hybrid Stochastic and Interval Parameters

Polynomial chaos employs a series of orthogonal polynomials with specific distribution in the uncertainty space to approximate the uncertainty processes.The original polynomial chaos proposed by Wiener [50] takes the Hermite polynomials in terms of Gaussian random variables as the trial basis to construct the PCE model.Cameron et al.[51] proved that such expansion converges to any second-order processes in theL2(Ω)sense,whereΩis the defined uncertainty space;this implies that any uncertain process can be approximated by the orthogonal polynomials of uncertainty variables inL2(Ω).The generalized polynomial chaos,or the Askey-chaos,was proposed in [13,52-54] and employed more types of orthogonal polynomials from the Askey scheme [54].

3.1 Basic Definition of Polynomial Chaos Expansions

The polynomial chaos expansion is essentially a representation of a functionf∈L2(Ω),whereΩis the properly defined probability space.Second-order random processes in [51] are processes with finite variance,and this can be applied to most uncertain event.Thus,a general second-order uncertain processy(x),viewed as a function ofx,i.e.,the uncertain event,can be represented in the form

whereaiandφi(ξ)correspond to the polynomial coefficients and basis functions in Eq.(7),respectively.

Note that Eq.(8) is an infinite and convergent summation in the infinite dimensional space ofξ.To reduce computational effort,it is often truncated at thepth-order and rewritten in the following form

where (S+1) is the number of retained terms,which is given by

The most important aspect of the above expansions is that an uncertainty process has been decomposed into a set of deterministic functions,and the truncated expansion series,p,mainly depends on the complexity and the nonlinearity of the uncertainty process.A second-order approximation is recommended as a first attempt;this approximation can be refined further using higher order terms,depending on the accuracy needs and available computational resources.

The polynomialsφi(ξ)satisfy the following orthogonality relation:

whereδijdenotes Kronecker delta function,〈,〉denotes the inner product in the Hilbert space of variablesξwhich can be expressed as follows:

Here,W(ξ)is the weight function determined by the distribution of the corresponding uncertainty variablesξ.

3.2 Polynomial Basis Functions for Hybrid Uncertainty Parameters

In Eqs.(8) and (9),the multi-dimensional,hypergeometric polynomialsφi(ξ)can be expressed as the tensor product of the corresponding one-dimensional polynomials bases:

whereφik(ξk) denotes the one-dimensional polynomial of thek-dimensional uncertainty variableξk.

For the stochastic parameters satisfying normal distribution,the one-dimensional polynomialφ(ξ)is the Hermite polynomialHn(ξ),which can be expressed as follows:

The (n+1)-th order Hermite polynomial satisfies the following relation:

hence the first few Hermite polynomials can be expressed as

The inner product of two Hermite polynomials can be formulated as

whereδmnis Kronecker delta function,and

Moreover,the Legendre polynomialLn(ξ)is used to describe the interval parameters,which is given by

The (n+1)-th order Legendre polynomial satisfies the following relation:

hence the first few Legendre polynomials can be expressed as

The inner product of two Legendre polynomials can be formulated as

whereW(ξ)=1.More polynomial bases,such as Jacobi,Laguerre,and generalized Laguerre polynomials can be found in the relevant literatures [52,53].

3.3 Galerkin Projection Method for the Expansion Coefficients

The Galerkin projection method [13,54] is adopted to identify the coefficientaiin Eq.(8).Applying Galerkin projection,both sides of Eq.(8) are simultaneously projected onto the orthogonal polynomialφi(ξ),such that

Owing to the orthogonality,the expansion coefficientaiin Eq.(8) can be calculated as follows:

It is noted that the denominator term of Eq.(24) is only the inner product of the orthogonal polynomialφi(ξ),which is the tensor product of the corresponding one-dimensional polynomials bases.Hence the denominator term can be regarded as the product of the inner products of the one-dimensional,orthogonal polynomial bases:

The numerator term of Eq.(24) is determined by calculating the expectation ofy(ξ)φi(ξ),which can be expressed as

wheref (ξ)is the joint probability density function ofξ.It is noted that this is a numerical integration of multivariate functions,and will be extremely complicated for the complex,nonlinear functions.Here,we use the sparse grid,numerical integration method (SGNIM) [46,55] to compute Eq.(26),while the sparse grid collocation strategy is utilized to generate the integral collocation nodes.For thed-dimensional uncertainty variablesξ,the sparse grid with thek-th collocation level can be expressed as follows:

whereis the one-dimensional integral nodes for each uncertain variable;i1,...,ij,...,andidare the multi-indices determining the number of collocation nodes in each dimension (Nj),that is,Njis a function ofij;kis the collocation level,and larger values ofkcorrespond to higher accuracy.and the inequalityk+1≤|i|≤k+dlimits |i|to a certain interval,so that all integral collocation nodes satisfying this condition can be found.The detailed formulas are not provided herein;the interested readers are referred to [55].

SGNIM performs the tensor product only on the one-dimensional integration corresponding to the multi-indices that meetk+1≤|i|≤k+d,thus avoids the “dimension disaster” caused by the direct tensor product.Such constraint can effectively remove those integral points in the full factor grid [56] that do not contribute significantly to the improvement of accuracy,thereby greatly reducing the computational cost.

Assuming that we obtain a total ofN d-dimensional integral collocation nodes through Eq.(27),i.e.,{L1,...,Ll,...,LN},the corresponding integral weightwlof thel-th collocation nodeLl∈can be expressed as follows:

whereis the product of the corresponding weights in each dimension of thel-thd-dimensional integral collocation node.

Moreover,the one-dimensional integral nodein Eq.(27) and integral weightin Eq.(28) (j=1,2,...,d) can be generally calculated through the moment-matching method [57].For the normal random parameters,andcan be easily obtained through Gaussian-Hermite quadrature formula,given by

whereandare the integral weights and nodes of Gaussian-Hermite integration,respectively.Similarly,theandof interval parameters can be computed by Gaussian-Legendre quadrature formulas.These weights and nodes of Gaussian quadrature can be easily found in the relevant literature [58],which is not discussed here.Taking the Gaussian-Hermite and Gaussian-Legendre as examples,the sparse grids in 2D space based on Gaussian abscissas are displayed in Fig.1.

Using thed-dimensional collocation nodes {L1,...,Ll,...,LN}and the corresponding weights{w1,...,wl,...,wN},Eq.(26) can be rewritten as follows:

Then substituting Eqs.(25) and (30) into Eq.(24),we can obtain the expansion coefficientai.

Figure 1:Sparse grids in 2D space based on Gaussian abscissas (a) Gaussian-Hermite,d=2,k=2.(b) Gaussian-Hermite,d=2,k=3.(c) Gaussian-Legendre,d=2,k=2.(d) Gaussian-Legendre,d=2,k=3

3.4 Evaluation of the Responses Bounds

After obtaining the PCE model,the statistical moments of the response can be directly derived through the expansion coefficients and the polynomial bases.The mean and standard deviation ofy(ξ)can be expressed as follows:

The lower and upper bounds of the response can be determined by the combination of MCS and PCE conveniently.In this process,the first step is to generate the samples of the hybrid uncertainty vectorξU.Then substitutingξUinto Eq.(9) to get the responses,we can obtain the probability density function (PDF) of the responses through the samples.In such circumstance,the percentile difference proposed by Du et al.[59] is usually used as an uncertainty measure,defined by

whereyα1andyα2are two values ofysatisfy the following probability relation:

α1andα2are,respectively,the cumulative distribution function (CDF) ofyat the left and right tail.For instance,the value ofα1can be taken as 0.05 or 0.01,while the value ofα2can be set as 0.95 or 0.99.yαiis the value ofythat corresponds to CDFαi.Fig.2 illustrates the concept of the percentile difference.More details can be found in [60].

Through the theory of percentile difference,the lower and upper bounds of the response can be regarded as the positions at the left and right tail of its distribution function,satisfying the given CDFαi(i=1,2).Hence the lower and upper bounds of the response can be obtained as follows:

where(·) is the inverse CDF ofy.

Figure 2:The concept of the percentile difference

4 Polynomial Chaos Expansion Approach for Dynamic Response Analysis of Interval ODEs

In this section,the direct MCS for dynamic response analysis will be reviewed briefly,and the polynomial chaos expansion approach for dynamic systems described by ODEs with hybrid uncertainty parameters will be introduced.A numerical example based on the two methods will then be compared to illustrate the advantages of the polynomial chaos expansion approach.

4.1 Direct MCS Method for Dynamic Response Analysis

When using numerical methods to solve the ODEs in Eq.(4),the range of timet?[t0,te] is discretized into a set of interpolation points (t0,t1,...,tj,...,te),and the solutions are calculated in turn at each time point.At any timetj,the direct MCS can be used to evaluate the response bounds of the dynamic systems with hybrid stochastic and interval uncertainty.The main solving process is summarized in Fig.3.

Figure 3:The solving process of the direct MCS method for dynamic response analysis

Step 1:Define the numerical solving method for ODEs,the number of samples for MCS,NM,the length of time period,te,the time step,Δt,the CDF of response at the left (right) tail,α1(α2),the distribution parameters of stochastic variables,xR,and the interval bounds of interval variables,xI.

Step 2:Initialize the cyclic counting indicesj=0 andi=1,wherejandiare introduced to count the iterations of time and the samples,respectively.

Step 3:Generate one group sample ofxRandxIfrom their distributions.

Step 4:Take the sample into the original ODEs,and use an appropriate numerical method to calculate the real response at timetj.

Step 5:Ifi=NM,continue.Else,assumei=i+1,and go to Step 3.

Step 6:Calculate the statistical moments,and obtain the PDF of the response at timetjthrough theNMsamples.

Step 7:Calculate the response bounds at timetjaccording to Eqs.(32)-(34).

Step 8:Ifj＞te/Δt,output the statistical moments and bounds of dynamic response with time.Else,assumej=j+1 andi=1,then go to Step 3.

4.2 The Computational Procedure of the Polynomial Chaos Expansion Approach

At timetj,the dynamic response of the ODEs in Eq.(4) with respect to the hybrid uncertainty vector can be approximated by the truncated polynomial chaos as shown in Eq.(9).The polynomial chaos model is an uncertainty function since the ODEs are with hybrid stochastic and interval parameters.To generate a PCE model at timetj,we must preset the polynomial basis functions according to the distributions of the uncertainty parameters,and specify the expansion series of polynomial chaos as well as the collocation levelkin SGNIM.Afteraiis obtained through the Galerkin projection method,the PCE model is independent of the degree of freedom of the structure and therefore,a computationally efficient surrogate model.In this case,combined with MCS,the response bounds at timetjcan be efficiently evaluated with a low computational cost.The flowchart of the solving process is plotted in Fig.4,which can be described as follows:

Step 1:Define the numerical solving method for ODEs,the expansion series of polynomial chaos,p,the collocation level,k,the number of samples for MCS,NM,the length of time period,te,the time step,Δt,the CDF of response at the left (right) tail,α1(α2),the distribution parameters of stochastic variables,xR,and the interval bounds of interval variables,xI.

Step 2:Determine the polynomial basis functions according to the distributions of the uncertainty parameters.

Step 3:Utilize the sparse grid collocation strategy to generate the integral collocation nodes and the corresponding integral weights according to Eqs.(27) and (28).

Step 4:Initialize the cyclic counting indicesj=0.jis introduced to count the iterations of time.

Step 5:Take the collocation nodes into the original ODEs,and use an appropriate numerical method to calculate the real responses at timetj.

Step 6:Calculate the expansion coefficient of polynomial chaos through the Galerkin projection method according to Eqs.(23)-(26) and (30) to yield the polynomial chaos at timetj.

Step 7:Calculate the statistical moments of the response at timetjaccording to Eq.(31).

Step 8:GenerateNMgroup samples ofxRandxIin the corresponding standard uncertainty spaces.

Step 9:Take theNMgroup samples into the obtained PCE model,and obtain the approximate responses.

Step 10:Obtain the PDF of the responses at timetjthrough theNMsamples.

Step 11:Calculate the response bounds at timetjaccording to Eqs.(32)-(34).

Step 12:Ifj＞te/Δt,output the statistical moments and bounds of dynamic response with time.Else,assumej=j+1,and go to Step 5.

In a brief description,the proposed polynomial chaos expansion algorithm for solving the dynamic systems described by ODEs with hybrid stochastic and interval parameters is similar to a type of sampling method.It transforms the original ODEs with uncertainty parameters into ODEs with deterministic parameters through sampling,without changing their expressions.As a nonintrusive method,any numerical methods can be used to solve the differential equations.Thus,the polynomial chaos expansion approach can be regarded as a computationally convenient method.

Figure 4:The solving process of the polynomial chaos expansion approach for dynamic systems described by ODEs with hybrid uncertainty parameters

4.3 A Numerical Example and Discussions

A simple numerical example is used here to demonstrate the validity of the present method.The results obtained by the polynomial chaos method which is abbreviated as PCEM are compared with those of the direct MCS.The schematic of a double pendulum system is analyzed as depicted in Fig.5.The differential equations of this system can be expressed as follows:

wherem1andm2are the masses of the two pendulums,respectively;l1andl2are the lengths of the pendulum rods,respectively;θiandωi(i=1,2) are respectively the angles and the angular velocities of the pendulum rods,andgis the gravity acceleration.Here,we assume thatm1andm2satisfy Gaussian distribution,with the mean and standard deviation of 1.0 and 0.1 kg,respectively.Note that Gaussian distribution is unbounded distribution with infinite interval,it means that the parameter values will be negative to some extent.Therefore,the modeling of stochastic variables is an approximate situation.l1andl2are considered as interval parameters,and their interval ranges are [0.45,0.55] m,and [0.9,1,1] m,respectively.The initial conditions[θ1(t0),θ2(t0),ω1(t0),ω2(t0)] are [π/3,3π/5,0,0].

Figure 5:The schematic of a double pendulum system

For this problem,the fifth-order polynomial chaos approach is used to solve the differential equations in the time period of 0-10 s.The collocation levelkin SGNIM is set as 4,thus the total number of samples in PCEM are 999,and the number of samples for MCS is 10000.The direct MCS is used to generate the reference solutions for comparison,in which the sample size is 10000.Moreover,α1andα2in the two algorithms are set as 0.001 and 0.999,respectively.The results of the angles and the angular velocities of the two pendulum rods are plotted in Figs.6 and 7,respectively.To show it in a more intuitive way,detailed information of the angular velocities of the first and second pendulum rods are listed in Tabs.1 and 2 at 10 time instants from 1 to 10 s.All calculations are completed on a desktop personal computer with 16 GB RAM and a CPU core of Intel(R) Core(TM) i5-4590 CPU clocked at 3.30 GHz.

As shown in Figs.6 and 7 and Tabs.1 and 2,the response bounds of the double pendulum system yielded by PCEM almost match the results of the direct MCS in most time period.Therefore,we can summarize that the PCEM can be considered as an alternative algorithm with fine precision when calculating the response bounds of dynamic systems with hybrid stochastic and interval uncertainty.Furthermore,Figs.6 and 7 show that the ranges of the upper and lower bounds ofθ1,θ2,ω1,andω2gradually expand with time.The dynamic responses of the double pendulum system are very sensitive to the uncertainty parameters.Therefore,the uncertainty of the dynamic response deserves more attention in structure design and analysis.

Moreover,the computations of PCEM and direct MCS take 128.67 s and 1347.42 s,respectively.Clearly,the computational efficiency of PCEM is much higher than that of the direct MCS.Therefore,the proposed method can be considered as an effective and efficient uncertainty analysis method for the dynamic systems with hybrid stochastic and interval uncertainty parameters.

Figure 6:The response bounds of the angles in a time period 10 s (a) the angle of the first pendulum rod (b) the angle of the second pendulum rod

Figure 7:The bounds of the angular velocities in a time period 10 s (a) the angular velocity of the first pendulum rod (b) the angular velocity of the second pendulum rod

Table 2:The detailed information of ω2 (rad/s) in a time period 10 s

5 The Uncertainty Analysis of Artillery Exterior Ballistic Dynamics

The schematic of the artillery exterior ballistic dynamics model is analyzed as depicted in Fig.8.The artillery exterior ballistic is the most typical uncertainty process in the artillery launching process,which is affected by various uncertainties such as propellant parameters,projectile characteristic parameters,artillery muzzle vibration,meteorological environment,operation error,and so on.These uncertainties will inevitably affect the final dynamics performance,thus the flight trajectory and motion attitude of each projectile is different.These differences eventually converge to the projectile impact point,and affect a very important tactical indicator of artillery system,namely artillery firing dispersion.The existence of uncertainties is the primary cause projectile dispersion.The uncertainty analysis of artillery exterior ballistic dynamics model can explore the impact of the uncertainty parameters on the projectile flight,which is of great significance for achieving effective control of the artillery firing dispersion.

Figure 8:The schematic of the artillery exterior ballistic dynamics model

A general six-degree of freedom exterior ballistic model can be formulated as follows:

wherevis the projectile velocity,θis the path angle,ψ1andψ2are,respectively,the vertical and the horizontal deflection angle;(X,Y,Z) indicate the position of centroid of the projectile in the ground coordinate corresponding to the firing distance,height,and offset,respectively;δ1andδ2are,respectively,the vertical and horizontal attack angle;γis the rotation angle,andβis the rotation angular velocity;Ω1andΩ2are the vertical and the horizontal angle of oscillation,respectively,while ˙Ω1and ˙Ω2are the corresponding angular velocities;ICis the polar moment of inertia;IAis the equatorial moment of inertia;bx,by,bz,kxz,kzz,ky,andkzare the aerodynamic parameters of the projectile,which are related to the shape and size of the projectile,such as the projectile massmp,the maximum diameter of projectiledmax,and the distance from the projectile centroid to its toplcg.Detailed auxiliary equations of the aerodynamic parameters can be found in [61].wx2,wy2,andwz2are the projections of wind speed on the velocity coordinate,which are given by

wherewxandwzare the longitudinal wind and the horizontal wind,respectively.

Through the above equations,we can obtain the flight trajectory and motion attitude of the projectile.Assuming that there are uncertain parameters in this dynamics model,the detailed uncertainty types and values are shown in Tab.3.This process has a total of 14 uncertainty parameters,which is a typical,high-dimensional,uncertainty analysis problem.

Table 3:Uncertainty parameters for an exterior ballistic model

The second-order polynomial chaos approach is used to solve the differential equations in the time period of 0-107 s.The collocation levelkin SGNIM is set as 2,thus the sample size for PCEM is 632,and the number of samples for MCS is 10000.The direct MCS is used to generate the reference solutions for comparison,in which the sample size is 10000.α1andα2in the two algorithms are set as 0.001 and 0.999,respectively.Partial results including the velocityv,the rotation angular velocityβ,and the horizontal deflection angleψ2of the projectile are shown in Fig.9.Detailed information ofvandψ2is listed in Tabs.4 and 5,respectively.Moreover,Fig.10 gives the spatial ranges of the projectile flight trajectory.

Figs.9 and 10 show that the response bounds of the artillery exterior ballistic dynamics model yielded by PCEM perfectly match the results of the direct MCS perfectly,which is further indicated intuitively by Tabs.4 and 5.Moreover,it is noted that the number of hybrid uncertainty parameters is up to 14.The results show that the proposed method has fine precision even in the high-dimensional,uncertainty problem.Furthermore,the computational time of direct MCS is 124879.612 s,while the computing time of PCEM is only 4521.29 s.This example demonstrates that high precision and low computational cost can be achieved by the proposed method when calculating the response bounds of artillery dynamics with hybrid stochastic and interval uncertainty.

Figure 9:The response bounds of an artillery exterior ballistic dynamics model in a time period 107 s (a) the velocity of the projectile (b) the rotation angular velocity of the projectile (c) the horizontal deflection angle of the projectile

Table 4:The detailed information of v (m/s) in a time period 107 s

Table 5:The detailed information of ψ2 (rad) in a time period 107 s

Furthermore,Fig.9 shows that the dynamic behaviors obtained by the proposed method are ranges rather than deterministic ones.It fully reflects the fluctuation of these dynamic behaviors caused by uncertainty.Hence the upper and lower bounds of these dynamic behaviors can provide a good reference for the artillery exterior ballistic design.More importantly,the uncertainty analysis of artillery external ballistic dynamics can predict the projectile flight trajectory under uncertain parameters.Fig.10a gives the spatial ranges of the projectile flight trajectory derived by the two methods,and the shadow area in Fig.10b is the possible range of the projectile impact point,so that the projectile dispersion caused by uncertainty parameters can be predicted.

Figure 10:The spatial ranges of the projectile flight trajectory (a) flight trajectory bounds in ground coordinate (b) projection of projectile flight in XZ plane

6 Conclusions

In this paper,a hybrid stochastic and interval uncertainty analysis method with polynomial chaos expansion is proposed to evaluate the response bounds of artillery dynamics.In this method,the Hermite polynomial in stochastic space and the Legendre polynomial in interval space are employed,respectively,as the trial basis to expand the Gaussian and interval processes.The polynomial coefficients are calculated through the Galerkin projection method,in which the sparse grid,numerical integration method is used to generate the integral points and the corresponding integral weights.Through the sampling in polynomial chaos expansion,the original artillery dynamics systems with hybrid stochastic and interval parameters can be transformed into ones with deterministic parameters,without changing their expressions.The yielded polynomial chaos expansion model is utilized as a computationally efficient surrogate model,and the supremum and infimum of the dynamic responses at each iteration time step can be easily approximated through Monte Carlo simulation and percentile difference.A numerical example and an artillery exterior ballistic dynamics model were used to illustrate the feasibility and efficiency of this approach.The numerical results indicate that the dynamic response bounds obtained by the polynomial chaos expansion approach almost match the results of the direct Monte Carlo simulation,but the computational efficiency of the polynomial chaos expansion approach is much higher than direct Monte Carlo simulation.Moreover,the proposed method also exhibits fine precision even in highdimensional uncertainty analysis problems.Another advantage of the polynomial chaos expansion approach is that it is non-intrusive.In a brief description,the proposed algorithm is similar to a type of sampling method.It has no special restrictions on numerical methods for solving the differential equations.Therefore,this method can be potentially applied to other artillery dynamics systems with hybrid uncertainty parameters,such as the artillery multi-body dynamics system described by the differential-algebraic equations (DAEs).However,further research is required to explore the potential of the proposed method in dealing with other artillery dynamics system.

Funding Statement:This research was financially supported by the National Natural Science Foundation of China [Grant Nos.301070603 and 11572158].Besides,the authors wish to express their many thanks to the reviewers for their useful and constructive comments.

Conflicts of Interest:The authors declare that they have no conflict of intrest to report regarding the present study.