Nonintrusive Uncertainty Quantification for automotive crash problems with VPS/Pamcrash
Marc Rocas, Alberto García-González, Sergio Zlotnik, Xabier Larráyoz, Pedro Díez
NNonintrusive Uncertainty Quantification forautomotive crash problems with VPS/Pamcrash
Marc Rocas(1,2), Alberto Garc´ıa-Gonz´alez(1), Sergio Zlotnik(1,3),Xabier Larr´ayoz(2), Pedro D´ıez(1,3)1- Laboratori de C`alcul Num`eric, E.T.S. de Ingenier´ıa de Caminos,Universitat Polit`ecnica de Catalunya – BarcelonaTech2- SEAT, Martorell, Barcelona3- The International Centre for NumericalMethods in Engineering, CIMNE, BarcelonaFebruary 16, 2021
Abstract
Uncertainty Quantification (UQ) is a key discipline for computationalmodeling of complex systems, enhancing reliability of engineering simula-tions. In crashworthiness, having an accurate assessment of the behaviorof the model uncertainty allows reducing the number of prototypes andassociated costs. Carrying out UQ in this framework is especially chal-lenging because it requires highly expensive simulations. In this context,surrogate models (metamodels) allow drastically reducing the computa-tional cost of Monte Carlo process. Different techniques to describe themetamodel are considered, Ordinary Kriging, Polynomial Response Sur-faces and a novel strategy (based on Proper Generalized Decomposition)denoted by Separated Response Surface (SRS). A large number of un-certain input parameters may jeopardize the efficiency of the metamod-els. Thus, previous to define a metamodel, kernel Principal ComponentAnalysis (kPCA) is found to be effective to simplify the model outcomedescription. A benchmark crash test is used to show the efficiency of com-bining metamodels with kPCA.Keywords: Uncertainty Quantification, Crashworthiness, SeparatedResponse Surface (SRS), kernel Principal Component Analysis (kPCA),Kriging, Surrogate modeling.
Uncertainty Quantification (UQ) in crashworthiness simulations is becomingan important asset to verify the design of vehicle structures. An important1 a r X i v : . [ s t a t . M E ] F e b hallenge for automotive industry is guaranteeing safety and reducing costs, andthis requires virtually testing a large number of tentative designs, and assessthe dispersion of the results of the virtual model due to the uncertain inputdata. UQ is particularly relevant for the crashworthiness analysis, where manyuncertainties have to be tackled together and then propagated in the simulationof this extremely complex, nearly chaotic, phenomenon.The growth and universal accessibility to computational resources and therobustness of codes, offer a good perspective on the possibilities of producingUQ for complex problems. However, in the case of crashworthiness simulations,a single simulation takes up to 18 CPU hours in a High Performance Comput-ing facility. Thus, the very large number of queries associated with a standardUQ process are practically unaffordable in this context. One viable alterna-tive is using a surrogate model (or metamodel) build upon a reduced numberof full-order simulations (denoted as training set ), see [20, 24, 17] for differ-ent approaches and comparative analyses. Still, the viability of metamodelsis limited by the number of input parameters: a large number of parametersresults in a highly multidimensional input space and therefore the engineer isafflicted by the so-called curse of dimensionality . If the parametric model isalready low-dimensional (3-4 design parameters), the actual threat is not the curse of dimensionality but dimensionality reduction is still necessary to com-putationally afford the simulation process for decision making . This is commonin crashworthiness and in the example included here for illustration.The idea is to determine a low number of relevant parameters (as combina-tions of the original ones) properly representing all the variability of the dataset.Principal Component Analysis (PCA) is the standard dimensionality reductiontechnique, to be used if the data structure is such that the low-dimensional sub-set where the data is contained (also referred as manifold) is linear. Other man-ifold learning techniques identify nonlinear low-dimensional structures. Amongthem, kernel Principal Component Analysis (kPCA) is considered here becauseit is one the simplest approaches, see [8] for a synthetic presentation. Thecombination of dimensionality reduction with surrogate modeling is a commonstrategy to carry out UQ in different disciplines and contexts [12, 13, 18].This paper analyzes the combination of different alternatives for dimen-sionality reduction techniques and surrogate models for UQ in crashworthinesssimulations with uncertain input parameters. Among the different techniquesexplores, the novel combination of kPCA and Separated Response Surface (SRS)demonstrates interesting properties. Other strategies are also used to define sur-rogate model fitting the training set like Polynomial Regression Surface (PRS)and Ordinary Kriging (OK). To do this, we used Monte Carlo sampling (asit is the simplest method for statistical analyses) in two steps of the proposedmethodology: 1) to obtain the training set for dimensionality reduction and sur-rogate model reconstrucction, and 2) once the low dimensional surrogate modelswere properly developed, standard Monte Carlo is performed (practically withno computational cost) to increase the probabilistic resolution in the descriptionof the QoI.This paper is structured as follows: section 2 presents the benchmark prob-2em for crashworthiness. The proposed UQ methodology is described in section3, which is divided in the three subsections. First, the main ideas of the kPCAtechnique are recalled. Next, the three different surrogates under considerationare briefly introduced (SRS, OK and PRS). Finally, the different surrogates arereadily used to quantify the uncertainty of the output. Section 4 illustrateshow the methodologies presented in the previous section perform for the crashproblem proposed in Section 2. Finally, section 5 includes some concludingremarks. The benchmark crash problem under consideration is illustrated in Fig. 1. Itcorresponds a reduced test model, ideally reproducing the main characteristicsof the simulation of a B-pillar, a well-known structural component of cars. Thisparticular benchmark test is used for different research studies in the Volkswa-gen Group. For the sake of saving computational cost and time, this model isoften used to test new materials, adhesives, welding spots or other conditionsbecause, due to the simplicity of the model, the numerical response requires acomputation of approximately 20 minutes. Even so, computing the solutions ofthe final training set with 2366 samples used here required around 788 hours(approximately 32 days) of computational time in one of the SEAT clusters. Be-sides, the computational time required for a standard full crash model is aroundone day per simulation. Thus, any effort in devising strategies to build a reliabletraining set with the minimum number of full-order solutions is worthwhile.The driving force in the model is provided by the impactor (green zone inFig. 1), that crashes at a speed of 50 mm/s against the vertical profile (redzone) during one second. This velocity might seem to be low, nevertheless thisis a standard benchmark test of interest for the Volkswagen Group accountingfor rate effects in the model. The reason is that the high deformations areconcentrated locally in the structure in specific areas and times, so that localrate effects are significant.The three structural parts are plates made of laminated steel sheet manufac-tured by cold folding. All the parts are joined with a structural adhesive bond,its material properties are characterized by Volkswagen.The structure is modelled using the Belytschko-Tsay shell element (which is astandard option in the VPS/Pamcrash package, very popular in the automotivesector, in particular for SEAT engineers) with one integration point in the plane.The impactor is considered to be a rigid body. The complete model has a totalof 13908 nodes (with 6 degrees of freedom).The Quantity of Interest (QoI) to be analyzed is the final plastic strain inthe 142 shell elements of the area depicted in black in Fig. 1.The numerical solver is implemented in VPS/Pamcrash [1], with the shellfinite element discretization mentioned above and an explicit time steppingscheme to solve the dynamical problem. The displacements of the points atthe ends of the horizontal profile are prescribed to zero (points marked with3igure 1 Crash benchmark. Thicknesses h , h and h are the three input ran-dom parameters corresponding to the vertical profile (red), horizontal profile(orange) and plate profile(blue). The impactor (green), and the area of ele-ments of interest (black) are also depicted.green arrows in Fig. 1). The contact between the different components of thestructure are treated with the surface-surface model defined in VPS/Pamcrash.Thicknesses h , h and h of the three parts of the structure are consid-ered to be stochastic parameters, that is random variables collected in vector h = [ h , h , h ] T . Their aleatory character is associated with the imperfec-tions produced during the manufacturing process. Random variables h , h and h are assumed to be normal and uncorrelated, that is h i ∼ N ( µ i , σ i ) andcov( h i , h j ) = 0, for i, j = 1 , ,
3. In each of the three parts, the correspondingthickness is considered to be constant. Besides, the three thicknesses h , h and h are modelled as having the same mean µ = µ = µ = 1 . σ = σ = σ = 0 .
12 mm.In order to build a training set, and as a first assessment of the stochas-tic behaviour of the system, a number of n s = 2366 Monte Carlo realizations(or samples) are performed. Thus, n s values of the input parameters h i , for i = 1 , , . . . , n s are generated with a random number generator and the corre-sponding VPS/Pamcrash solutions are obtained. These solutions (in particularthe vectors containing the plastic strain in the d = 142 elements of the zone ofinterest) are collected in a training set matrix X = [ x x · · · x n s ] ∈ IR d × n s , whereeach column x i = [ x i . . . x id ] T is the VPS/Pamcrash solution corresponding toinput h i . The actual QoI is the average plastic strain in the zone, is represented4y a form l ( · ), and for each h i and x i reads l ( x i ) = 1 d d (cid:88) j =1 x ij . Note that the Monte Carlo process with n s = 2366 samples is considered hereas a reference, and it is only obtained in the academic example under considera-tion. The number of full-scale computations affordable for a real problem in theautomotive industrial practice is much lower. A methodology to select the op-timal number of samples (the lowest providing some insights in the assessmentof uncertainty) is discussed in section 4.The points constituting the training are often taken as realizations of randomvariables (Monte Carlo) as indicated above. They can be also taken as pseudorandom (e.g. Latin hypercube) or deterministic (e.g. Hammersley points, Hal-ton sequences) [25]. The number of samples n s affordable in a real problem is generally not suffi-cient to produce a proper Monte Carlo assessment of the statistical propertiesof the output of the system. A review of nonintrusive UQ methodologies forcrashworthiness, see [21], demonstrates that the standard Monte Carlo samplingis extremely demanding and, in practice, beyond the possibilities of standardindustrial practitioners.As indicated in the previous section, the standard Monte Carlo approachconsists in generating random samples of the input, running the model andretrieving statistics of the output (or any QoI). This is what corresponds to theupper part (black arrows) in the scheme of Fig.2.However, the part of the standard model (also denoted as full-order, herecomputed with VPS/Pamcrash) is too computationally expensive to be per-formed for the number of samples providing statistical relevance. Thus, thealternative is to replace this full-order model by a surrogate, that is a simplefunctional transformation from h to x . The surrogate is created using a trainingset consisting in data generated by the full-order model.An additional difficulty is encountered due to the high-dimension of theoutcome of the model, x . It is complicated to create a high-dimensional func-tional approximation having a target space of d (here 142) dimensions. Thus,previous to undertake the determination of the surrogate, it is convenient toapply some dimensionality reduction technique. In the context of crashworthi-ness simulation, the data generated by the models are often adopting nonlineardata structures [8, 23]. Thus, it is expected to require nonlinear dimensionalityreduction (kPCA).The QoI is introduced as an essential indicator for decision making. The QoIsummarizes the information contained in x . Quantifying the uncertainty of the5oI is sufficient to take some decisions. For instance, to verify the crashwor-thiness response of the structural design. Uncertainty Quantification of high-dimensional objects like x is cumbersome and the outcome is difficult to useas a tool supporting decision making. In that sense, the stochastic assess-ment focuses in a low-dimensional (even purely scalar) QoI, rather than in ahigh-dimensional object like x . However, a deeper analysis of the phenomenonrequires understanding the underlying mechanisms associated with the overallmechanical response of the system. In that sense, all the information containedin x is pertinent. The fact that the model order reduction strategy is able torecover back the full-order object in as accurately as possible is therefore ex-tremely advantageous. In this aspect, kPCA behaves much better than PCA inmany cases: the simple QoI is fairly approximated by the PCA reduction, butkPCA improves the mapping back to the original variable x .All these aspects are covered in the methodology described in remainder ofthe paper, having the following steps: • Creation of training set.
Generate n s realizations of the input pa-rameters h i , for i = 1 , , . . . , n s , and compute the corresponding full-ordersolutions x i (that constitute the training set). • Dimensionality reduction.
Analyze the training set and find theprincipal components allowing to reduce the dimensionality of the familyof solutions. In practice, this boils down to apply PCA or kPCA anddetermine a mapping between the solutions x ∈ IR d and some new variable z (cid:63) ∈ IR k in a much lower-dimensional space ( k (cid:28) d ). The mappingbetween x and z (cid:63) is to be characterized forward and backward. ThekPCA backward mapping is found to be more accurate than with PCA.That is, kPCA recovers with much more accuracy a full-order x associatedwith a reduced-order coordinate z (cid:63) . Although this advantage is often notperceptible when assessing a low-dimensional (or scalar) QoI, a proper x recovery is crucial to deepen in the mechanical interpretation of the results.For instance, to identify the mechanisms associated with the differentmodes of the probability distribution. • Surrogate model.
The functional dependence z (cid:63) = F ( h ) is deter-mined from the data provided by the training set, and the dimensionalityreduction. • Complete Monte Carlo UQ (using surrogate).
Once the surrogate F ( · ) is available, for each input value h , the corresponding z (cid:63) is straight-forwardly computed as F ( h ). Then the backward mapping produces thecorresponding x , and l ( x ) is the associated QoI. The concatenation of thethree operations is computationally affordable. Therefore standard MonteCarlo can be performed with a sufficient number of realizations.The different aspects of the devised methodology are described in detail inthe following sections. It is important noting that, among the four conceptualsteps mentioned above, the computational cost is concentrated in the creation6igure 2 Schematic illustration of the methodology.of the training set. Obtaining this representative collection of solutions requiresa computational time in the range of weeks or months, depending on the type ofsimulation in crashworthiness. The other steps: dimensionality reduction, sur-rogate modeling and Monte Carlo UQ represent in practice a negligible amountof computational efforts (in the order of seconds). The training set matrix X = [ x x · · · x n s ] ∈ IR d × n s is seen as a set of n s pointsin a d -dimensional space. The idea of dimensionality reduction is to find asubspace of lower dimension k (cid:28) d where the set of points is contained.The PCA strategy consists in diagonalizing the square d × d matrix X T X (covariance matrix), that is finding U ∈ IR d × d such that X T X = UΛU T (1)where Λ is a diagonal matrix with eigenvalues λ ≥ λ ≥ · · · ≥ λ d . Thedimension is reduced from d to k if the last d − k eigenvalues are negligible withrespect to the k first. In this case, the new variable selected is z (cid:63) = U (cid:63) T x , (2)being U (cid:63) ∈ IR d × k the matrix with the first k columns of U . Eq. (2) describesthe forward mapping, that is how to map the high-dimensional vector x intothe element z (cid:63) reduced dimensional space, from dimension d to dimension k .The backward mapping goes in the opposite direction and reads x = U (cid:63) z (cid:63) . (3)Thus, PCA is a straightforward methodology relying in the fact that thetraining set is lying in a linear subspace. In many cases, the structure of the7ow-dimensional manifold where the solution ranges is nonlinear and more so-phisticated dimensionality reduction techniques are required. The kPCA is analternative reformulation based on the traditional PCA, in fact, it performsPCA in a new feature space where the data is transformed from the originalspace [22].The main characteristics of the kPCA are explained in detail in [8] and sum-marized here. It is assumed that some transformation Φ from IR d to a higher-dimensional space is able to flatten the training set. That is the transformedtraining set (cid:8) Φ ( x ) , Φ ( x ) , . . . , Φ ( x n s ) (cid:9) is such that PCA is able to discover alinear subspace of dimension k . In other words, the transformation Φ maps thenonlinear manifold (of dimension k ) where the training set ranges into a linearsubspace.The transformation Φ that produces this effect is a priori unknown. How-ever, it is worthy trying with some different alternatives and see what is thereduced dimension they propose: the best choice for Φ is the one producingthe lower value of k . Moreover, in practice Φ is indirectly characterized usingthe kernel trick . Thus, instead of describing directly Φ , an expression for thebivariate form κ ( · , · ) is provided, assuming that the following relation between κ ( · , · ) and Φ ( · ) holds κ ( x i , x j ) = Φ ( x i ) T Φ ( x j ) (4)for i, j = 1 , , . . . , n s .A classical choice for the kernel κ ( · , · ) is the so-called Gaussian kernel, definedas κ ( x i , x j ) = exp( − β (cid:13)(cid:13) x i − x j (cid:13)(cid:13) ) (5)where β is a parameter that, in the applications of this paper, is taken equal to0.1.Having the kernel at hand, one may compute a matrix equivalent to XX T for the samples transformed by Φ . This matrix is denoted by G ∈ IR n s × n s andhas generic coefficient [ G ] ij = κ ( x i , x j ) (6)It is worth noting that the eigenvalues of XX T are the same of those of X T X ,which are the ones extracted in (1). Actually, G is also readily diagonalized andthe following factorization is obtained G = V ˜ ΛV T (7)where ˜ Λ contains the same non-zero eigenvalues that would be obtained fromdiagonalizing the corresponding covariance matrix, which is not available.Thus, the eigenvalues λ ≥ λ ≥ · · · ≥ λ n s are computed, and the reduceddimension k is selected such that the last n s − k eigenvalues are negligible withrespect to the k first.Once k is obtained, the original variable x ∈ IR d is mapped into a variablein the reduced space, z (cid:63) ∈ IR k using the following the expression z (cid:63) = V (cid:63) T g ( x ) , (8)8here V (cid:63) ∈ IR n s × k is the matrix with the first k columns of V and g ( x )IR n s isa vector with generic component[ g ( x )] i = κ ( x i , x ) (9)for i = 1 , . . . , n s .As described in detail in [8], if the samples transformed by Φ are not centred,some corrections have to be done and both matrix G and vector g have to bemodified accordingly. These corrections are straightforward and are omittedhere for the sake of a simpler presentation.Equations (8) and (9) characterize the forward kPCA mapping, from x to z (cid:63) . The backward mapping for kPCA is not as simple as for the PCA versiondescribed in equation (3). A point z (cid:63) in the reduced space is mapped back toa point x which is recovered as a weighted average of the points of the trainingset, namely x = n s (cid:88) i =1 w i ( z (cid:63) ) x i , with weights such that n s (cid:88) i =1 w i ( z (cid:63) ) = 1 (10)The weights w i ( z (cid:63) ) are computed such that the forward mapping of x is as closeas possible to z (cid:63) . A popular strategy to compute these weights with a simpleapproach is to use a radial basis interpolation concept based on the distances of z (cid:63) to the images of the sample points, z (cid:63) i for i = 1 , . . . , n s . That is computing d i = (cid:107) z (cid:63) − z (cid:63) i (cid:107) and taking any value of w i ( z (cid:63) ) decreasing with d i , for example w i ( z (cid:63) ) ∝ d i . As already announced, the dimensionality reduction techniques presented aboveare a previous step to build a surrogate model. The training set is now usedto approximate the functional dependence associated with the full-order model.The final goal is to compute x as an easy-to-evaluate function of h , that isthe surrogate. With the dimensionality reduction, this is split in two steps: asurrogate from h to z (cid:63) plus the backward mapping from z (cid:63) to x , see Fig. 2.Here, the surrogates are presented as generic methodologies to establish afunctional dependency among some input h and some output function y ( h ) (weuse y to account for any output, that could be either x or z (cid:63) ). Obviously, doingthe surrogate with z (cid:63) has the advantage of dealing with a much lower dimension(number of components) of the model output. In practice, for the sake of asimpler presentation, a scalar output Y ( h ) is considered in the following, thatstand, for example, for any of the components of y ( h ). In the examples, Y coincides with the first component of the reduced space using kPCA, that is Y = [ z (cid:63) ] .In the following, the parameters describing the stochastic input space wherethe function takes values are collected in the vector h = [ h . . . h n d ] T ∈ IR n d .Where n d is the number of stochastic dimensions of the problem ( n d = 3 in thebenchmark under consideration). 9hus, the goal is to approximate the functional dependence Y ( h ) using theimages of the points of the training set y k = Y ( h k ), k = 1 , . . . , n s , where . Allthe sample points are collected in the vector y = [ y y ...y n s ] T . The idea of SRS is to find a separated approximation F ( h ) to Y ( h ). Theseparated character of F ( h ) means that it is a sum of rank-one terms, beingeach rank-one term the product of sectional modes (the adjective sectional isused to indicate that the mode depends only on one of the parameters). Thealgorithm employed to compute the SRS is based on the ideas of the least-squares PGD approximations described in detail in [3, 4, 15, 14].Thus, F ( h ) reads F ( h ) = n f (cid:88) j =1 σ j n d (cid:89) i =1 f ji ( h i ) (11)where each sectional mode f ji ( h i ) is represented in some discrete sectional space.The discrete sectional space is generated by a family of functions (cid:8) Ψ i ( h i ) Ψ i ( h i ) . . . Ψ i n i ( h i ) (cid:9) being n i the dimension of the sectional function space.Accordingly, sectional modes have the following expression f ji ( h i ) = n i (cid:88) m =1 a mi Ψ im ( h i ) (12)where the unknown coefficients a mi , for i = 1 , . . . , n d and m = 1 , . . . , n i , have tobe computed to determine the sectional mode f ji ( h i ), for j = 1 , , . . . Different alternatives are available as the approximation space defined by Ψ im ( h i ). Here we have considered a finite element discretization with n i nodes(in 1D domains, n i − F ( h ). Thus, F ( h ) ≈ Y ( h ) is taken such that it minimizes (cid:107) F ( h ) − Y ( h ) (cid:107) = (cid:104) F − Y, F − Y (cid:105) = n s (cid:88) k =1 w k ( F ( h k ) − y k ) (13)where the weights w k are introduced to assimilate the sum into an integral, thatis, to assume that (cid:90) Ω h F ( h ) d h ≈ n s (cid:88) k =1 w k F ( h k )note that the associated scalar product (cid:104)· , ·(cid:105) of two arbitrary functions F and G reads (cid:104) F, G (cid:105) = n s (cid:88) k =1 w k F ( h k ) G ( h k ) . (14)Note that weights w k , k = 1 , , . . . , n s must be selected corresponding to aquadrature having as integration points h k , where Y is known. Typically, the10istribution of points h k is provided by a stochastic sampling and cannot beenforced a priori by the user to obtain his/her preferred quadrature (e.g. aGauss-Legendre quadrature or a composite Simpson’s rule). Thus, the weightsof the quadrature are adapted to optimize the integration order in a (multidi-mensional) Newton-Cotes fashion.Thus, least-squares solution in a linear functional space V is readily com-puted as a projection, that is finding F ∈ V such that (cid:104) F, F (cid:63) (cid:105) = (cid:104) Y, F (cid:63) (cid:105) for all F (cid:63) ∈ V (15)Note that integral equation (15) has to be fulfilled for any weighting function(or test function) F (cid:63) in V , as in the standard weak form of a boundary valueproblem.The key aspect of any PGD algorithm is how to solve the rank-one approx-imation. That, is how to find an approximation to Y ( h ) with a function of theform F ( h ) = σ n d (cid:89) i =1 f i ( h i ) (16)which is a particular case of (11) with just one term.The standard PGD strategy consists in an alternate direction approach, thatis to compute the sectional mode f γ , the rest of the sectional modes f i for i (cid:54) = γ are assumed to be known. Thus, in practice, F ( h ) and F (cid:63) ( h ) are taken as F ( h ) = f γ ( x γ ) (cid:89) i (cid:54) = γ f i ( h i ) (17)and F (cid:63) ( h ) = f (cid:63)γ ( h γ ) (cid:89) i (cid:54) = γ f i ( h i ) . (18)This alternate directions strategy leads to a sectional problem, reduced to the γ coordinate. The family of sectional problem is to be solved sequentially for γ = 1 , , . . . , n s , and then iterated until convergence is reached.For the sake of simplifying the writing, the computable term depending onall the sectional modes but γ is denoted as T γ and T kγ when evaluated in h k ,namely T γ := (cid:89) i (cid:54) = γ f i ( h i ) and T kγ := (cid:89) i (cid:54) = γ f i ( h ki ) . (19)Thus, the sectional counterpart of (15) reads (cid:10) f γ T γ , f (cid:63)γ T γ (cid:11) = (cid:10) Y, f (cid:63)γ T γ (cid:11) for all f (cid:63)γ (20)that is n s (cid:88) k =1 w k ( T kγ ) f γ ( x kγ ) f (cid:63)γ ( h kγ ) = n s (cid:88) k =1 w k T kγ y k f (cid:63)γ ( x kγ ) . (21)11sing a particular case of (12), that is f γ ( h γ ) = n γ (cid:88) m =1 a mγ Ψ γm ( h i ) (22)in (21) and taking f (cid:63)γ ( h kγ ) = Ψ γ(cid:96) ( h kγ ) for all (cid:96) = 1 , . . . , n γ yields n s (cid:88) k =1 w k ( T kγ ) n γ (cid:88) m =1 a mγ Ψ γm ( h kγ ) Ψ γ(cid:96) ( h kγ ) = n s (cid:88) k =1 w k T kγ y k Ψ γ(cid:96) ( h kγ ) (23)or n γ (cid:88) m =1 (cid:34) n s (cid:88) k =1 w k ( T kγ ) Ψ γm ( h kγ ) Ψ γ(cid:96) ( h kγ ) (cid:35)(cid:124) (cid:123)(cid:122) (cid:125) M (cid:96) m a mγ = n s (cid:88) k =1 w k T kγ y k Ψ γ(cid:96) ( h kγ ) (cid:124) (cid:123)(cid:122) (cid:125) f (cid:96) (24)for all (cid:96) = 1 , . . . , n γ . That is, a linear system of n γ equation with n γ unknowns M a γ = f . (25)Once the sectional approximation is obtained solving (25), the loop in al-ternate directions iterations is continued until convergence and completion ofthe rank-one computation. As usual in PGD [9], once the rank-one solutionis obtained, the greedy approach aims at computing the next term (next j in(11)).As it is standard in this type of strategies, in order to compute an approx-imation having the separated form given in (11), there are three nested loops.First, the greedy approach (loop in j ) aims at computing rank-one terms hav-ing the form given in (16). Then an alternated direction iterative scheme isapplied consisting in two nested loops: the iterative loop to reach convergence(not described explicitly in this text with an iteration index) and an inner loopfor γ = 1 , , . . . n s , ranging all sectional dimensions. This is standard in thereferences describing any PGD scheme, see [4] for an algorithmic description.As mentioned above, functions Ψ im in (12) are chosen as classical C finiteelements shape functions. Contrary to other choices (e.g. high-order polyno-mials) these type of functions are more stable due to their local support butintroduce a lack of smoothness (jumps in the first derivatives, singularities inthe second derivatives). Consequently, when using a finite element approxima-tion for the sectional modes, it is important having the possibility of enforcingthe smoothness of the solution. This is equivalent to penalize in system (25),the non-smoothness of the sectional function described in (22).This requires, for instance, penalizing some postprocessed quantity of thesectional mode f γ ( h γ ), represented by the vector of nodal values, a γ . Thequantity to be penalized, the lack of smoothness, is represented by a matrix G mapping the nodal values of f γ ( h γ ) into the postprocessed quantity in somerepresentative points. In the following, G is taken as the standard gradient operator, computing the derivatives of f γ ( h γ ) in the integration points of the12lements of the mesh. Thus, G is a n G × n γ matrix, being n G the number ofintegration points in the mesh (assuming that the dimension of the sectionalspace is 1). Thus, the measure of the lack of smoothness that has to be reducedis given by a T γ G T Ga γ . Provided that system (25) is equivalent to minimize thefollowing functional 12 a T γ M a γ − f T a γ Enforcing the smoothness requires minimizing the perturbed functional12 a T γ M a γ − f T a γ + λ a T γ G T Ga γ for some value of the factor λ that states the importance of the smoothing. Thelarger is λ , the smoother is the recovered solution. This results in the followinglinear system, which is a modification of (25) accounting for the smoothing (cid:2) M + λ G T G (cid:3) a γ = f (26) Ordinary Kriging (OK) is an interpolation technique commonly used in engi-neering and originated for geostatistical problems [19]. The OK method deter-mines weights for a set of simulation points to calculate a prediction of a newsample. The weights are calculated with a variogram model that has the mainfeature to estimate variances for any distance. The kriging metamodel F ( h ) ofany point h is defined by: F ( h ) = n s (cid:88) i =1 w i y i , (27)where the unknowns w i are the weights and y i are the scalar values of thefunction to be interpolated. To determine the optimal values for the krigingweights, the variogram function plays an important role. There exist differentvariograms: Gaussian, exponential, linear among others [19]. The OK matrixsystem to obtain the weights reads, γ γ . . . γ n s γ n s γ n s . . . γ n s n s
11 1 . . . w ... w n s µ = γ ... γ n s . A specific condition for OK with respect to other Kriging methods is enforcingthe sum of weights equal to 1, (cid:80) n s i =1 w i = 1. This condition is achieved byintroducing the new unknown µ as Lagrange multiplier [16]. The entries of thematrix in the equation above depend on the variogram function γ evaluated foreach distance δ between a pair of samples, that is γ ij = γ ( δ ), being δ = (cid:107) h i − h j (cid:107) .The entries γ . . . γ n s are evaluations of the variogram γ between all the sample13oints with respect to the new (current) point. Here, we used the sphericalvariogram defined as: γ ( δ ) = (cid:40) C + C (cid:104) (cid:0) δa (cid:1) − (cid:16) δa (cid:17)(cid:105) , < δ ≤ aC + C , δ > a (28) C is the nugget constant representing the noise of the data, a is the rangeof the transition zone where the variogram levels off and the sill ( C + C ) isdefined as the total variance of the model. For the benchmark problem C = 0,in consequence C is the total variance of the model. In Fig.3 it is illustrated aspherical variogram function.Figure 3 Variogram with the three main parameters. The nugget C , the range a and the sill C + C . Polynomial Response Surface (PRS) has been applied in numerous studies tobuild metamodel for different engineering problems [7, 5, 10]. It consists in asimple multidimensional polynomial fitting. A second order polynomial model F ( h ) takes the form, F ( h ) = c + n d (cid:88) i =1 c i h i + n d (cid:88) i =1 c ii h i + n d − (cid:88) i =1 n d (cid:88) j = i +1 c ij h i h j , (29)where h i is the i -th stochastic input, the different coefficients c are the unknownsto be computed, collected in a vector c . If the approximation was able tointerpolate the data, the following linear system should be solved: Ac = y , (30)14here A is the matrix containing the values of the different interpolation func-tions in (29), that is, for n d = 3, (cid:8) , h , h , h , ( h ) , ( h ) , ( h ) , h h , h h , h h (cid:9) in the sample points h k , for k = 1 , . . . , n s . This results in A = h h h (cid:0) h (cid:1) (cid:0) h (cid:1) (cid:0) h (cid:1) h h h h h h h h h (cid:0) h (cid:1) (cid:0) h (cid:1) (cid:0) h (cid:1) h h h h h h ... ... ... ... ... ... ... ... ... ...1 h n s h n s h n s ( h n s ) ( h n s ) ( h n s ) h n s h n s h n s h n s h n s h n s , (31)The model is often non-interpolative, with more equations (points in the sam-ple) than unknowns (number of coefficients in c ). Therefore system (30) cannotbe solved exactly but using a least squares minimization criterion, that is min-imizing the Euclidean norm of the residual, namely (cid:107) y − A c (cid:107) . This results intaking the vector of unknown coefficients solution of( A T A ) c = A T y . (32)This method presents drawbacks for high-dimensional data and data with oscil-lations. Increasing the order of the Polynomials may improve accuracy. How-ever, for high-order approximations Runge’s phenomenon creates instabilitiesand wrong predictions [2]. Once the surrogate model is available, the Monte Carlo UQ assessment with alarge number of samples n MC is produced at an affordable computational cost.Thus, for each of the three metamodels introduced above (SRS, OK andPRS) n MC realizations h , h , . . . , h n MC are produced and the corresponding valueof of the mean, variance (and standard deviation) and Probability Density Func-tion (PDF, to be approximated as a histogram) is readily estimated:Mean = E [ F ( h )] = 1 n MC n MC (cid:88) k =1 F ( h k ) (33)Variance = σ = 1 n MC − n MC (cid:88) k =1 (cid:0) F ( h k ) − E [ F ( h )] (cid:1) , (34)The PDF corresponding to Y = F ( h ) is denoted by f Y ( y ) and it is approximatedby histogram p Y ( y ) computed on the basis of the n MC Monte Carlo samples y k = F ( h k ), for k = 1 , , . . . , n MC . Note that histogram p Y ( y ) is a piecewiseconstant function defined over a partition in uniform intervals of the Y domain,15 Y = (cid:83) n Y (cid:96) =1 I (cid:96) . Piecewise constant function p Y is such that for y ∈ I (cid:96) , p Y ( y ) isequal to the number of samples y k lying in I (cid:96) divided by n s .Each response surface is used to generate the images of the n MC = 50000samples of the input space h ∈ IR n d , that is h → Y ). The backward mappingtechnique ( Y → X → QoI ) described in Section 3.1 is used to obtain thestatistics of the QoI.Comparing the obtained values of mean and variance is straightforward be-cause they are scalar values. However, comparing PDFs is not as trivial. Here,the Kullback Leibler divergence technique is proposed as a criterion to comparethe PDF functions.
Kullback-Leibler (KL) divergence is used as a comparative criterion for the dif-ferent resulting Monte Carlo PDFs of each metamodel. This quantity measuresdifferences between two PDF functions [6]. Two random variables F and G havePDFs f and g . The KL divergence is introduced as a distance that quantifiesif the two random variables are similar enough. The random variable F and itsPDF f are taken as the reference and g is considered to be an approximationto f . Thus, KL divergence between the two continuous PDFs f and g reads: D KL ( F (cid:107) G ) = (cid:90) ∞−∞ f ( y ) log (cid:18) f ( y ) g ( y ) (cid:19) dy. (35)Note that equation (35) is associated with the notion of entropy and it is inter-preted as the relative entropy or the information gain from G to F .In the case the PDFs are replaced by their discrete counterparts, that ishistograms, instead of f and g , one has histograms p Y and q Y with the formatdescribed in the previous section. Thus, p Y and q Y are expressed as the valuesof the probability of being in each of the n Y bins, that is p Y ( y (cid:96) ) and q Y ( y (cid:96) ), for (cid:96) = 1 , , . . . , n Y and y (cid:96) ∈ I (cid:96) . The discrete counterpart of equation (35) reads D KL ( p Y (cid:107) q Y ) = n Y (cid:88) (cid:96) =1 p Y ( y (cid:96) ) log (cid:18) p Y ( y (cid:96) ) q Y ( y (cid:96) ) (cid:19) . (36)The values obtained using the discrete KL divergence introduced in equation(36) depend on the number of bins, n Y . In order to normalize these values, anormalizing constant is introduced providing a reference to understand wetherthe resulting discrete KL divergence is actually small enough. Note that theKL divergence is seen as a distance but it is not conceived as the norm of adifference. Thus, it is not possible to normalize dividing directly by the normof p Y (or f in (35)). In order to obtain a reference value, we propose taking the distance of p Y to the less informative distribution, that is the uniform histogram q U such that q U ( y (cid:96) ) = n Y , for (cid:96) = 1 , , . . . , n Y . The rationale behind this choiceis taking q U as the zero or absolute reference distribution. This value is denoted16s D KL and reads D KL = D KL ( p Y (cid:107) q U ) = n Y (cid:88) (cid:96) =1 p Y ( y (cid:96) ) log (cid:0) n Y p Y ( y (cid:96) ) (cid:1) . (37)Note that the quantity D KL defined in equation (37) is actually the entropyof p Y . Dividing the figures obtained with the KL divergence of (36) by D KL provides a relative value that allows elaborating a more informed criterion todecide if p Y and q Y are sufficiently close to each other. In this section, the combination methodology of kPCA + SRS is applied for thebenchmark crash problem. Also, OK and PRS are applied to compare the SRSresults.
Initially, a certain amount of VPS/Pamcrash FE simulations is required to re-construct the input matrix X for kPCA manifold analysis and dimensionalityreduction (the training set computed in an offline phase). These simulations arethe initial samples for the data analysis. A key issue is quantifying the numberof samples required to obtain enough and credible information to describe thelow-dimensional manifold containing the solution. As described above, we devisethe combined use of kPCA to reduce the dimensionality, different techniques tobuild a response surface, and the KL divergence as a measure to compare thedifferent probability distributions and stop enriching the sampling. This processis illustrated in the flowchart scheme shown in Fig.4, and it is detailed next:Figure 4 Flowchart scheme to select the number of simulation n s for the kPCAinput matrix X .The eigenvalues of the kernel matrix G = V ˜ ΛV T in kPCA (equation 7),show the quantity of information collected by each associated eigenvector, as17xplained in section 3.1 and more detailed in [8]. Summarizing, the largesteigenvalue measures the largest amount of information collected by the corre-sponding eigenvector. For instance, for the first eigenvalue of matrix ˜ Λ , theassociated eigenvector is the first column of matrix V . Adopting the first com-ponent as reduced model (keeping only one principal component, the first one),the d = 142 dimensions of the training set samples (each column of matrix X )are reduced to one scalar number and the n s samples are stored in the vector y = [ y y ...y n s ] T , see sections 3.1 and 3.2 for more details.In the case of our benchmark crash problem, reducing to one dimensioncollects more than 80% of information of the manifold where data belong. This80% figure allows considering as admissible in this context, and in agreementwith the resulting approximations, the very advantageous reduction to a singledimension. This figure is calculated along the sampling refinement process (fordifferent values of n s ) to check the behaviour of the quantity of informationretained in the one-dimensional reduction. This is shown in Fig.5, where it canbe noticed that even using only 100 samples, the first eigenvalue collects already80% of information.Figure 5 Quantity of information [%] stored by the first eigenvector by increas-ing the number of samples n s .At this point, the number of samples n s is required to guarantee some sta-tistical accuracy. The KL divergence is used to compare the subsequent distri-butions of probability of the QoI, obtained with the training sets correspondingto the different values of n s , see Fig. 4. That means, to validate the final valueof n s = 2366, to be used in the input matrix X for analysis (instead of thestarting value of n s = 100). That is, the histogram obtained by the values of y for a low number of samples n s is compared by the KL divergence criterionwith the histogram for a higher number of samples. This process is repeated byincreasing the number of samples until the value of the KL divergence becomes18maller than 10 − (considered low enough for our required accuracy). The re-sults obtained in this process are detailed in Fig.6, where it is clear how thehistograms become more stable increasing the number of samples, and conse-quently the KL-Div value decreases. For the final number of samples n s = 2366,the KL-Div value is below the prescribed tolerance. Recalling expression (37),the calculated value of D KL for the final histogram in Fig.6 is D KL = 0 . . n s = 2366.Additionally, the low dispersion of the results by using one reduced dimensionis confirmed (the first eigenvalue contains 82 .
61% of information, stored in y =[ y y ...y n s ] T ). Moreover, the consistency of the backward mapping from vector y x and then calculating the corresponding QoI is also confirmed by the resultsin Fig.8b. (a) (b) Figure 8 (a) Reference values of histogram, mean, variance and standard de-viation of the first principal component y of kPCA. This reference values areachieved with 2366 simulations in VPS/Pamcrash. (b) Scatter plot aroundthe identity function (red) of the QoI with respect to the approximated forthe first principal component y . The first principal component y is linked to its corresponding values of h =[ h , h , h ] T . Fig. 9 shows the scatter plot between the reduced space and theinputs h , h . Moreover, two clusters with different density are observed. Theinput h is discarded by the criterion of Spearman Correlation coefficient (SpC)[11]. The dependences between the first principal component with respect to20he inputs are: • SpC ( y , h ) = 0 . • SpC ( y , h ) = 0 . • SpC ( y , h ) = − . h shows a very small correlation, and therefore is discarded forsurrogate modeling.Figure 9 Scatter plot between the reduced space y and the inputs h , h . In Fig.10 it is shown the response surface of SRS, OK and PRS metamodelsbetween the first principal component of kPCA (storing 82 .
61% of information)and the inputs h , h . The metamodels F OK and F SRS show adaptive behaviourwith the sample points (in blue). However, the F PRS metamodel exhibits unstabletails in areas where there are few samples from the training set. This problem iscalled Runge’s phenomenon and is common due to lack of sampling in the tails ofthe distributions. To evaluate the behaviour and the robustness of the responsesurfaces, each surface is evaluated increasing new random samples until 50000points, aiming to compare histograms, means and standard deviations withrespect to the reference values plotted in Fig.8a.In Fig.11 the histogram, mean, variance and standard deviation results of50000 new random samples are illustrated for the evaluation of each metamodel.The three metamodels give approximations to the mean with an error around5% of the Standard Deviation. Since the current reference mean (Figure 8a) iszero, both positive and negative values can be expected. The histograms of the F SRS and F OK metamodels show similar bimodal distribution with respect to thereference histogram illustrated in Fig. 8a. However, the two modes of distri-bution are not captured with the F PRS metamodel. The Runge’s phenomenon of21 a)(b)(c)
Figure 10 In red, is shown the response surfaces of a) SRS, b) OK and c) PRS.In blue, the scattering samples.the response surface and a worse adaptation to the sampling points overlooksthe two distribution modes.In Fig. 12 the convergence of the three surrogate models are compared withthe reference values while new random points for each metamodel are increasedup. The results of this comparative study shows similar results in terms of meanand standard deviation for the three techniques. Meaning that these statisticalvalues are not sensitive for the criterion to select the best surface. However,the results for the KL divergence clearly shows a worse behaviour for the PRS22ethod caused for the tails of the response surface. In contrast, SRS and OKhave similar results for the KL divergence where a very good performance isobserved.
Once the surrogates F ( · ) are available, for each input value h , the correspond-ing z (cid:63) is straightforwardly computed as F ( h ). Then, the backward mappingproduces the corresponding input vector x of plastic deformation values in thearea of interest, bieng l ( x ) its the associated QoI. The concatenation of thethree operations is computationally negligible with respect to the cost of thetraining set of full order simulations. At this point, standard Monte Carlo isperformed with 50000 new random samples for h , h and h to evaluate eachmetamodel. In Fig. 13 it is presented the corresponding PDFs of the QoI forthe metamodels (SRS, OK and PRS). A bimodal function with approximately19% of probability for the small mode and 81% for the big mode can be appre-ciated for SRS and OK. Otherwise, PRS fails to capture such behaviour. Thestatistical variables of the QoI for each metamodel are presented in Table 1.Here the three variables present similar results, which means that any meta-model captures similar information in terms of mean, variance and standarddeviation.Recalling that kPCA improves the mapping back to the original variable x , a physical interpretation of the bimodal PDF can be performed. The cor-responding behaviour of the structure for each mode of the PDF is illustratedin Fig. 14. Clearly, two physical modes are observed. The first snapshot ofFig. 14 shows the higher values of plastic strain and a significant back-bend ofthe plate profile. Otherwise, the second snapshot shows lower values of plasticdeformation with a more rigid behaviour of the plate. The first physical casepresent 81% of probability and the second case 19% of probability occurrence.Table 1 Statistical variables for each metamodel. Mean Variance StDReference values
SRS OK PRS
A nonintrusive methodology to perform Uncertainty Quantification for crash-worthiness problems is presented. The basic idea is to combine some dimen-sionality reduction technique (here kPCA) with a surrogate model based on a23raining set of full-order solutions (ideally not too many, because of their com-putational cost). The dimensionality reduction eases the task of the surrogatemodel and enables the analyst to detect clusters and categorize the data. Thesurrogate model (or metamodel) substitutes at a negligible computational costthe original full-order model. It therefore permits producing multiple queriesto the model, corresponding the different parametric input values demanded bythe Monte Carlo strategies.In the benchmark problem considered, kPCA allows describing the full phe-nomenon with only one principal component, accounting for more than 82%of the total variance (that is, of the information). This problem is relevant inautomotive engineering (and often used as benchmark by SEAT engineers) and,despite the fact that only three input parameters are assumed to have stochasticnature (and 3 dimensions are not awakening the curse of dimensionality ), thedimensionality reduction is still pertinent to simplify the output of interest tobe analyzed. Actually, using kPCA, only one principal component is accountingfor more than 82% of the total information. It also detects two clusters cor-responding to two deformation modes and two different levels of the QoI. TheUQ methodology is also providing the probabilities of occurrence of these twomodes, which are 19% and 81%. This is reflected in a bi-modal PDF, one modehaving a probability four times larger than the other. Moreover, using kPCAas dimensionality reduction strategy the backward mapping from the reducedspace is more accurate and allows interpreting the mechanisms associated withthese two modes.In the presented methodology for a crash problem, the use of linear PCAis also a suitable option. On the one side, if a single scalar QoI is required fordecision making, PCA is simpler than kPCA. On the other side, if it is necessaryto find both an accurate QoI, as well as a more detailed approximation of the fulloriginal variable x map, corresponding to a complete deformation field, kPCAimproves the mapping back to the original input space, since it accounts on theintrinsic nonlinearities involved in the manifold of training set. Thus, both PCAor kPCA can be used for this methodology, depending on the main objectivesought. In this manuscript, for the reasons mentioned above, kPCA is used anddescribed in more detail. This allows dealing with problems representing morecomplex phenomena, where data lies in highly nonlinear manifold.Three formats of the surrogate models are taken into consideration, OrdinaryKriging (OK), Polynomial Response Surface (PRS) and a Separated ResponseSurface (SRS) approach, introduced here and based in the PGD methodology.The assessment of the mean and variance of the outcome (the QoI) is properlycomputed using the three alternative surrogates. However, when it comes toanalyze the PDFs (approximated by histograms), the SRS and OK surrogatesperform much better than the PRS. The PRS surrogate fails to capture thebi-modal character of the PDF.Being OK an interpolative methodology (the response surface passes throughthe data of the training test), it is pretty sensitive to the noise contained in thedata. In the current examples, this is not an important issue, because the datais not particularly noisy. However, it may be relevant in other cases. SRS24eing a least-squares fitting it is not suffering of this drawback. Moreover,SRS is proposing an explicit parametric solution, therefore it can be used tocompute derivatives or to integrate it analytically. This allows also to computethe statistical moments, probability density function and cumulative densityfunction with analytical methods, circumventing the Monte Carlo sampling.Another interesting feature of the SRS is its fair scalability with the number ofinput parameters (stochastic dimension).The combination of the kPCA manifold learning technique with the differentsurrogates offers an attractive framework to perform UQ in complex problems.Here, the application to parametric crashworthiness simulations opens new per-spectives. The available alternatives for the surrogates (in particular OK andSRS) and the dimensionality reduction techniques at hand, are a powerful tool-box allowing to attack challenging problems in science and engineering. Thecombination of dimensionality reduction and surrogate models produces ac-curate solutions at an affordable computational cost, accounting also for theuncertainty, that is assessing the credibility of the simulation. Particularly inthe context of crashworthiness UQ, the computational cost is a key issue and adriving force for the research developments in the field. Obviously, increasingaccuracy requires a higher computational effort. Finding a trade-off betweenthese two factors is a daily concern for research engineers. This paper intendsto provide tools to achieve accurate and credible crashworthiness industrial sim-ulations at an acceptable computational effort. Acknowledgments
This work is partially funded by Generalitat de Catalunya (Grant Number 1278SGR 2017-2019 and Pla de Doctorats Industrials 2017 DI 058) and Ministeriode Econom´ıa y Empresa and Ministerio de Ciencia, Innovaci´on y Universidades(Grant Number DPI2017-85139-C2-2-R).
Compliance with Ethical Standards
The authors declare that they have no conflict of interest.
References [1]
PAM-SCL - Theory Notes Manual . 2000.[2] John P Boyd and Fei Xu. Divergence (runge phenomenon) for least-squarespolynomial approximation on an equispaced grid and mock–chebyshev sub-set interpolation.
Applied Mathematics and Computation , 210(1):158–168,2009. 253] Pedro D´ıez, Sergio Zlotnik, Alberto Garc´ıa-Gonz´alez, and Antonio Huerta.Algebraic PGD for tensor separation and compression: an algorithmic ap-proach.
C. R. Mec. , 346(7):501–514, 2018.[4] Pedro D´ıez, Sergio Zlotnik, Alberto Garc´ıa-Gonz´alez, and Antonio Huerta.Encapsulated PGD Algebraic Toolbox Operating with High-DimensionalData.
Arch. Comput. Method Eng. , 2020. to appear.[5] Hongbing Fang, Masoud Rais-Rohani, Zheny Liu, and MF Horstemeyer. Acomparative study of metamodeling methods for multiobjective crashwor-thiness optimization.
Computers & structures , 83(25-26):2121–2136, 2005.[6] David J Galas, Gregory Dewey, James Kunert-Graf, and Nikita A Sakha-nenko. Expansion of the kullback-leibler divergence, and a new class ofinformation metrics.
Axioms , 6(2):8, 2017.[7] Shawn Gano, Harold Kim, and Don Brown. Comparison of three surrogatemodeling techniques: Datascape, kriging, and second order regression. In , page 7048, 2006.[8] Alberto Garc´ıa-Gonz´alez, Antonio Huerta, Sergio Zlotnik, and PedroD´ıez. A kernel principal component analysis (kpca) digest with a newbackward mapping (pre-image reconstruction) strategy. arXiv preprintarXiv:2001.01958 , 2020.[9] Hasini Garikapati, Sergio Zlotnik, Pedro D´ıez, Clemens V Verhoosel, andE Harald van Brummelen. A proper generalized decomposition (pgd) ap-proach to crack propagation in brittle materials: with application to ran-dom field material properties.
Computational Mechanics , 65(2):451–473,2020.[10] Anthony Giunta and Layne Watson. A comparison of approximationmodeling techniques-polynomial versus interpolating models. In , page 4758, 1998.[11] Jan Hauke and Tomasz Kossowski. Comparison of values of pearson’s andspearman’s correlation coefficients on the same sets of data.
Quaestionesgeographicae , 30(2):87–93, 2011.[12] Christos Lataniotis, Stefano Marelli, and Bruno Sudret. Extending clas-sical surrogate modelling to ultrahigh dimensional problems through su-pervised dimensionality reduction: a data-driven approach. arXiv preprintarXiv:1812.06309 , 2018.[13] Min Li, Ruo-Qian Wang, and Gaofeng Jia. Efficient dimension reductionand surrogate-based sensitivity analysis for expensive models with high-dimensional outputs.
Reliability Engineering & System Safety , 195:106725,2020. 2614] Y Lu, N Blal, and A Gravouil. Adaptive sparse grid based hopgd: Towarda nonintrusive strategy for constructing space-time welding computationalvademecum.
International Journal for Numerical Methods in Engineering ,114(13):1438–1461, 2018.[15] Ye Lu, Nawfal Blal, and Anthony Gravouil. Multi-parametric space-timecomputational vademecum for parametric studies: Application to real timewelding simulations.
Finite Elements in Analysis and Design , 139:62–72,2018.[16] Tomislav Malvi´c and Davorin Bali´c. Linearity and lagrange linear multipli-cator in the equations of ordinary kriging.
Nafta: exploration, production,processing, petrochemistry , 60(1):31–43, 2009.[17] Maliki Moustapha, Bruno Sudret, Jean-Marc Bourinet, and Benoˆıt Guil-laume. Metamodeling for crashworthiness design: comparative study ofkriging and support vector regression. In
Uncertainties 2014-proceedingsof the 2nd International Symposium on Uncertainty Quantification andStochastic Modeling, July 7th to July 11th, 2014, Rouen, France . ETH-Z¨urich, 2014.[18] Joseph B Nagel, J¨org Rieckermann, and Bruno Sudret. Uncertainty quan-tification in urban drainage simulation: fast surrogates for sensitivity anal-ysis and model calibration. arXiv preprint arXiv:1709.03283 , 2017.[19] MA Oliver and R Webster. A tutorial guide to geostatistics: Computingand modelling variograms and kriging.
Catena , 113:56–69, 2014.[20] Na Qiu, Yunkai Gao, Jianguang Fang, Guangyong Sun, Qing Li, andNam H Kim. Crashworthiness optimization with uncertainty from sur-rogate model and numerical error.
Thin-Walled Structures , 129:457–472,2018.[21] M Rocas, A Garc´ıa-Gonz´alez, X Larr´ayoz, and P D´ıez. Nonintrusivestochastic finite elements for crashworthiness with vps/pamcrash.
Archivesof Computational Methods in Engineering , pages 1–26, 2020.[22] Bernhard Sch¨olkopf, Alexander Smola, and Klaus-Robert M¨uller. Nonlin-ear component analysis as a kernel eigenvalue problem.
Neural computation ,10(5):1299–1319, 1998.[23] Laurens Van Der Maaten, Eric Postma, and Jaap Van den Herik. Dimen-sionality reduction: a comparative.
J Mach Learn Res , 10(66-71):13, 2009.[24] Tao Wang, Liangmo Wang, Chenzhi Wang, and Xiaojun Zou. Crashwor-thiness analysis and multi-objective optimization of a commercial vehicleframe: A mixed meta-modeling-based method.
Advances in mechanicalengineering , 10(5):1687814018778480, 2018.[25] Tien-Tsin Wong, Wai-Shing Luk, and Pheng-Ann Heng. Sampling withhammersley and halton points.
Journal of graphics tools , 2(2):9–24, 1997.27 a)(b)(c)
Figure 11 Histogram, mean and standard deviation results of the differentsurrogate models. a) SRS, b) OK and c) PRS. The results are obtained byevaluating each surrogate model with 50000 new random samples.28 a) (b)(c)
Figure 12 Convergence plots of SRS, OK and PRS evaluating KL divergence,mean and standard deviation with respect to the reference values(KL=0,mean=0, standard deviation=1.5428 plotted with the dashed line).29 a) (b)(c) (d)
Figure 13 a) shows the histogram of the 2366 reference samples. b), c) andd) shows the histograms of the QoI by evaluating 50000 random samples forSRS, OK and PRS metamodels. 30 a)(b)a)(b)