Hybrid simulation is a method used to investigate the dynamic response of a system subjected to a realistic loading scenario. The system under consideration is divided into multiple individual loading-rate-sensitive substructures, out of which one or more are tested physically, whereas the remaining are simulated numerically. The coupling of all substructures forms the so-called hybrid model. Although hybrid simulation has beenextensivelyusedacrossvariousengineeringdisciplines,itisoftenthecasethatthe hybridmodelandrelatedexcitationisconceivedasdeterministic. However,associated uncertainties are present whilst simulation deviation due to their presence could be signiﬁcant. To this regard, global sensitivity analysis based on Sobol’ indices can be used to determine the sensitivity of the hybrid model response due to the presence of the associated uncertainties. Nonetheless, estimation of the Sobol’ sensitivity indices requires unaﬀordable amount of hybrid simulation evaluations. Therefore, surrogate modelingtechniquesareusedtoalleviatethisburden. Inthispaper,threediﬀerentsur- rogate modeling methods are examined, namely polynomial chaos expansion, Kriging and polynomial chaos Kriging. A case study encompassing a virtual hybrid model is employed and hybrid model response quantities of interest are selected. Their respec- tive surrogates are developed using all three aforementioned techniques. The Sobol’ indices obtained utilizing each examined surrogate are compared with each other and results highlight potential deviations. Keywords: real-timehybridsimulation,surrogatemodeling,polynomialchaosexpan- sion, Kriging, polynomial chaos Kriging, global sensitivity analysis, Sobol’ indices, uncertainty quantiﬁcation. 1 Introduction 1 Hybrid simulation (HS), also known as Hardware-in-the-Loop (HiL) or model-based 2 system testing, is a method used to investigate the dynamic response of a system 3 1 DRAFTsubjectedtoarealisticloadingscenario. Thesystemunderconsiderationisdividedinto 4 multiple individual loading-rate-sensitive substructures, out of which one or more are 5 tested physically, and thus correspond to the physical substructures (PS), whereas the 6 remainingsubstructuresaresimulatednumerically,namelythenumericalsubstructures 7 (NS).ThecouplingofPSandNSformstheso-calledhybridmodel. Theresponseofthe 8 latterisobtainedfromastep-by-stepnumericalsolutionoftheequationsgoverningthe 9 motion of the underlying hybrid model, combined with measurements acquired from 10 the PS. Therefore, the calculation of the next HS time step is computed online. An 11 arrangement of transfer systems, e.g. servo-controlled motors or actuators, is used to 12 deliver the motion between NS and PS. The employed transfer system accounts for the 13 substructure coupling and hence for synchronizing their dynamic boundary conditions 14 in every time step of the HS. From the coupling of substructures, several challenges 15 arise,e.g. timedelaysduetoinherenttransfersystemdynamicsorduetocomputational 16 power needed to compute the NS response. Advance control techniques [1, 2] and 17 model order reduction methods [3, 4] have been used to tackle such issues. Albeit 18 the challenges, the HS approach is beneﬁcial since it can be used to experimentally 19 study the inner workings of speciﬁc substructures over their linear regime and hence 20 acquirerealisticresultsbutwithoutconstructingtheentireconsideredsystemnorrisking 21 damagingit. AdetailedoverviewofHScanbefoundinthereportbySchellenbergand 22 Mahin [5]. 23 HS has been widely used across diﬀerent engineering disciplines [6, 7, 8, 9, 10]. 24 However,itisoftenthecasethatthehybridmodelandrelatedexcitationareconceived 25 as deterministic. Nonetheless, this assumption does not reﬂect to real-world struc- 26 tures and systems as loading is stochastic whereas hybrid model parameters may be 27 highly uncertain. A thorough investigation of all possible parameter and excitation 28 variations is not feasible considering the cost related to an individual HS evaluation. 29 Still, uncertainties are present whilst simulation deviation due to their presence could 30 be signiﬁcant. 31 In this regard, surrogate modeling has been proposed to perform global sensitivity 32 analysis (GSA) of a quantity of interest (QoI) hybrid model response with respect to 33 a set of input parameters originating from both substructures and excitation [11, 12]. 34 GSA aims to quantitatively determine the degree each input parameter aﬀects the se- 35 lected QoI of the hybrid model response. Sobol’ indices is a popular methodology to 36 performvariance-basedGSA,inwhichthevarianceoftheresponseQoIisdecomposed 37 asasumofcontributionsrelatedtoeachinputvariableorcombinationwitheachother 38 [13]. Monte Carlo simulations can be used to estimate the Sobol’ indices. Nonethe- 39 less, for each index estimation, O(103)model evaluations would be required [11]. To 40 alleviate this burden, surrogate models of the QoI response can be developed. Surro- 41 gates (a.k.a. metamodels) are used to substitute an expensive-to-compute model with 42 an inexpensive-to-compute surrogate. Therefore, Monte Carlo simulations using the 43 developedsurrogatescanbeperformedataﬀordablecostandSobol’indicesestimation 44 becomeseasiertocompute. Especially,forthecasethatapolynomialchaosexpansion 45 (PCE)hasbeenconstructed,thentheSobol’indicescanbecomputedanalyticallyatno 46 additional cost [14, 15]. 47 Inthisstudy,threediﬀerentsurrogatemodelingmethods,namelyPCE,Krigingand 48 polynomial chaos Kriging (PCK), are used to perform GSA via Sobol’ indices in HS. 49 2 DRAFTAcasestudyencompassingavirtualhybridmodelisemployedandspeciﬁcQoIofthe 50 hybrid model response are selected. Surrogates for the selected QoI are constructed 51 andtheirperformanceisevaluated. Sobol’indicesarecomputedusingeachdeveloped 52 surrogate and comparison with each other is made to highlight potential deviations. 53 Results demonstrate the eﬀectiveness of each examined surrogate modeling method. 54 This paper is organized as follows. Sections 2 and 3 present the background of 55 the examined surrogate modeling methods as well as the GSA with Sobol’ indices. 56 Section4presentsthecasestudywiththevirtualhybridmodelandSection5discusses 57 the respective results. Finally, Section 6 presents the overall conclusions of this study. 58 2 Surrogate modeling methods 59 Considering an input vector X∈ DX⊂RNand a computational model Y= 60 M(X)withY∈RN, machine-learning data-driven regression algorithms formu- 61 late a mapMs:X/mapsto→Ybased on an obtained sample set of input points X= 62/braceleftbig x(1),...,x(N)/bracerightbigTand of the respective output values, i.e. model evaluations or ex- 63 perimental measurements, Y=/braceleftbig y(1),...,y(N)/bracerightbigT. The set ofX,Yrealizations 64 corresponds to the so-called experimental design (ED). 65 2.1 Polynomial chaos expansion 66 PCE is a well-known uncertainty quantiﬁcation spectral method used to substitute the 67 dynamicsofanexpensive-to-computenumericalmodelwithaninexpensive-to-compute 68 surrogate (a.k.a. metamodel), representing the outputs of the model by a polynomial 69 functionofitsinputs[16,17]. Itisproventobeapowerfulsurrogatetechniqueusedin 70 a wide variety of engineering contexts, to replicate the dynamic response of complex 71 high-dimensional models [11, 18, 19]. 72 In more detail, given a random input vector Xwith independent components 73 expressedbythejointPDF fXandaﬁnitevariancecomputationalmodel Y=M(X), 74 such that E[Y2] =/integraltext DXM2(x)fX(x)dx<∞, the PCE ofM(X)follows: 75 Y=M(X)≈MPCE(X) =/summationdisplay α∈NNyαΨα(X). (1) The PCE function is built on the Ψα(X)multivariate orthonormal polynomial basis 76 with respect to the input vector fX. The degree of the Ψαpolynomials components 77 is identiﬁed by the α= (α1,...,αN),α∈NNmulti-index for each of the input 78 variables, while yαcorresponds to the polynomial coeﬃcients. 79 2.1.1 Polynomial basis 80 The multivariate polynomials are constructed as a tensor product of their univariate 81 orthonormal polynomials φ(i) k(x(i)), such that: 82 3 DRAFTΨα(X) =N/productdisplay i=1φ(i) αi(x(i)). (2) The latter meets the orthonormality criteria: 83 /angbracketleftBig φ(i) j(x(i)),φ(i) k(x(i))/angbracketrightBig =/integraldisplay DXiφ(i) j(x(i))φ(i) k(x(i))fXi(x(i))dxi=δjk,(3) whereicorresponds to the input variable, jandkto the polynomial degree, fXi(x(i)) 84 to theith-input marginal PDF and δjkto the Kronecker symbol. The selection of the 85 univariateorthonormalpolynomialfamiliesdependsonthemarginalPDFofeachinput 86 variable to which they are orthogonal, e.g. if an input variable follows the uniform/ 87 GaussiandistributionthentheLegendre/Hermiteorthogonalpolynomialfamilyisused 88 respectively for this speciﬁc input variable [17, 20]. 89 2.1.2 Truncation schemes 90 Once the univariate polynomials families are selected for each input variable, the next 91 step is the construction of the PCE following Eq. (1). However, because the sum of 92 Eq. (1) consists of inﬁnite terms, it is often truncated to a ﬁnite number of terms for 93 practical reasons. Hence, the truncated basis is deﬁned as A⊂NNand the PCE of 94 M(X)admits: 95 MPCE(X) =/summationdisplay α∈AyαΨα(X). (4) The performance of PCE is closely connected with the truncation scheme used. Un- 96 derﬁtting or overﬁtting is possible to occur when too many terms are discarded or 97 introduced respectively [19]. In the standard basis truncation scheme [17] the maxi- 98 mum degree p∈N+is deﬁned, such that the degree of each polynomial is capped by 99 thisvalue. Therefore,thestandardbasistruncationschemeconsistsof/parenleftbigN+p p/parenrightbig =(N+p)! N!p!100 elements and follows: 101 AN,p=/braceleftbig α∈NN:|α|≤p/bracerightbig , (5) whereNcorresponds to the number of input variables of M(X)and|α|=/summationtextN i=1αi102 tothetotaldegreeofallpolynomialsin Ψ. Tofurtherreducethepolynomialbasissize, 103 additional truncation schemes have been developed, namely the maximum interaction 104 and hyperbolic truncation scheme. 105 Inthemaximuminteractiontruncationscheme[20],thebasisofEq.(5)isreduced 106 such that theαindices to include at most rnon-zero components and thus the rank of 107 αis decreased [20]. Accordingly, Eq. (5) is written as: 108 AN,p,r=/braceleftbig α∈AN,p:||α||0≤r/bracerightbig , (6) where||α||0=/summationtextN i=11{αi>0}and 1{αi>0}is the indicator function. 109 4 DRAFTThehyperbolic(a.k.a. q-norm)truncationscheme[21]reformsthebasisofEq.(5) 110 such that: 111 AN,p,q=/braceleftbig α∈AN,p:||α||q≤p/bracerightbig , (7) where||α||q=/parenleftBig/summationtextN i=1aq i/parenrightBig1 qandq∈(0,1]. 112 2.1.3 PCE coeﬃcient calculation 113 Oncethetruncationschemeisdetermined,thenextstepisthecalculationofcoeﬃcients 114 y={yα,α∈A}of Eq. (4). Various methods have been developed to do this [20]. 115 Oneofthesetechniquesistheleast-squaresminimization[22]. Thisisannon-intrusive 116 method,meaningthatthecoeﬃcientsareobtainedafterpost-processingtheEDpoints. 117 More speciﬁcally, Eq. (4) is reformed as: 118 MPCE(X) =P−1/summationdisplay j=0yjΨj(X) +εP=yTΨ(X) +εP, (8) wherey= (y0,...,yP−1)Tdenote the PCE coeﬃcients, εPthe truncation error, 119 Ψ(X) ={Ψ0(X),..., ΨP−1(X)}Tthe multivariate orthonormal polynomials, and 120 P=/parenleftbigN+p p/parenrightbig . Then, the PCE coeﬃcients are obtained by solving: 121 ˆy= arg min y∈R|A|E/bracketleftBig/parenleftbig yTΨ(X)−M (X)/parenrightbig2/bracketrightBig . (9) The solution of Eq. (9) is attained by Ordinary Least-Squares (OLS) and follows: 122 ˆy=/parenleftbig ATA/parenrightbig−1ATY, (10) where the experimental matrix is Aij= Ψj(x(i))withi= 1,...,Nandj= 123 0,...,P−1. 124 However, for N < P−1theATAmatrix is non-invertible and thus the OLS 125 problem cannot be solved [19, 23]. Nonetheless, through sparse regression accurate 126 surrogates can be constructed with fewer coeﬃcients. A way to achieve sparse regres- 127 sioninhigh-dimensionsisbyamodiﬁcationoftheleast-squaresminimizationmethod. 128 In more detail, a penalty term is introduced in Eq. (9) and the latter is modiﬁed as: 129 ˆy= arg min y∈R|A|E/bracketleftBig/parenleftbig yTΨ(X)−M (X)/parenrightbig2/bracketrightBig +λ||y||1, (11) where||y||1=/summationtext α∈A|yα|andλa parameter of the penalty term. This penalty 130 term causes the minimization to facilitate sparse solutions. One of the widely utilized 131 techniques, used to solve the minimization problem of Eq. (11) is the Least Angle 132 Regression (LAR) algorithm [21]. The latter is also employed in the case study in 133 Section 4 for the respective PCE surrogate. 134 5 DRAFT2.2 Kriging 135 Kriging (a.k.a. Gaussian process) is a widely used surrogate modeling method. It 136 considers the output of a model Y=M(X)as a realization of a Gaussian process, 137 withx∈DX⊂RNand follows [24, 25]: 138 MKrig(x) =βTf(x)+σ2Z(x,ω), (12) whereβ= (β1,...,βP)andf(x) = (f1(x),...,fP(x))are coeﬃcients and regres- 139 sion functions, respectively. The product of the latter two corresponds to the mean 140 valuesoftheGaussianprocess,i.e. theKrigingtrend,while σ2toitsvariance. Z(x,ω) 141 is a stationary Gaussian process with zero mean and unit variance, deﬁned by an au- 142 tocorrelation function R/parenleftbig x,x/prime;θ/parenrightbig and its hyperparameters θ. The coeﬃcients βare 143 computed from the generalized least-square solution following: 144 β=/parenleftbig FTR−1F/parenrightbig−1FTR−1Y, (13) whereFij=fj(x(i))andRik=R/parenleftbig x(i),x(k);θ/parenrightbig the correlation matrix with i,k= 145 1,...,Nandj= 1,...,P. 146 In Kriging, the model Yand surrogate/hatwideY(x) =MKrig(x)response is assumed 147 to have a joint Gaussian distribution and therefore the mean value and variance of the 148 prediction/hatwideY(x)read [24, 26]: 149 µ/hatwideY(x) =f(x)Tβ+r(x)TR−1/parenleftbig Y−Fβ/parenrightbig , (14) σ2 /hatwideY(x) =σ2/parenleftBig 1−r(x)TR−1r(x) +u(x)T/parenleftbig FTR−1F/parenrightbig−1u(x)/parenrightBig ,(15) whereri(x) =R/parenleftbig x,x(i);θ/parenrightbig consists of the cross-correlations between predictions x 150 and sample points x(i)withi= 1,...,Nandu(x) =FTR−1r(x)−f(x). 151 2.2.1 Trend families 152 The ﬁrst step in constructing a Kriging surrogate is to determine its trend family. 153 Three common family types can be found in the literature [26, 27], namely the simple 154 Kriging expressed by βTf(x)=/summationtextP j=1fj(x), the ordinary with βTf(x)=β0and 155 the universal which follows βTf(x)=/summationtextP j=1βjfj(x). The latter is a generic type 156 admitting a variety of formations, e.g. multivariate polynomials (see Section 2.3). 157 2.2.2 Autocorrelation functions 158 The autocorrelation functions used in Kriging metamodeling represent the relative 159 position of two input sample points xandx/prime. Some common autocorrelation function 160 families [26, 27] are the linear: 161 R/parenleftbig x,x/prime;θ/parenrightbig =N/productdisplay i=1max/parenleftbigg 0,1−|x(i)−x/prime(i)| θi/parenrightbigg , (16) 6 DRAFTthe exponential: 162 R/parenleftbig x,x/prime;θ/parenrightbig = exp/parenleftbigg −N/summationdisplay i=1|x(i)−x/prime(i)| θi/parenrightbigg , (17) the Gaussian: 163 R/parenleftbig x,x/prime;θ/parenrightbig = exp/parenleftBigg −N/summationdisplay i=1/parenleftbigg|x(i)−x/prime(i)| θi/parenrightbigg2/parenrightBigg (18) and the Matérn: 164 R/parenleftbig x,x/prime;θ,v/parenrightbig =N/productdisplay i=11 2v−1Γ(v)/parenleftbigg 2√v|x(i)−x/prime(i)| θi/parenrightbiggv Kv/parenleftbigg 2√v|x(i)−x/prime(i)| θi/parenrightbigg , (19) wherev≥1/2correspondstotheshapeparameter, ΓistheEulerGammafunctionand 165 Kvis the modiﬁed Bessel function of the second kind. 166 2.2.3 Hyperparameters estimation 167 Followingtheselectionoftrendandautocorrelationfamily,thenextstepinbuildingthe 168 Kriging surrogate consists of estimating the unknown hyperparameters θ. One of the 169 existingmethodologiesisthemaximum-likelihoodestimation. Thegoalofthelatteris 170 to maximize the likelihood of model evaluations Yby determining the appropriate set 171 of the Kriging parameters β,σ2andθ. The likelihood function of Yadmits: 172 L(β,σ2,θ;Y) =/parenleftbig detR/parenrightbig−1/2 /parenleftbig 2πσ2/parenrightbigN/2exp/parenleftBig −1 2σ2/parenleftbig Y−Fβ/parenrightbigTR−1/parenleftbig Y−Fβ/parenrightbig/parenrightBig .(20) The hyperparameters θare computed by solving the optimization problem: 173 θ= arg min θ∈Dθ/parenleftbig −logL(θ;Y)/parenrightbig = arg min θ∈Dθ1 2/parenleftbig log(detR) +Nlog(2πσ2) +N/parenrightbig . (21) For a more detailed description of the derivation of Eq. (21) along with additional 174 methods for hyperparameter estimation, e.g. cross-validation estimation, the reader 175 should refer to [24, 27]. 176 2.2.4 Optimization methods 177 The solution of Eq. (21) is obtained by employing an optimization algorithm. Two 178 main categories of such algorithms are the local, e.g. gradient-based, and global, 179 7 DRAFTe.g. evolutionary algorithms, methods. Local methods need less objective functions 180 evaluations and can converge faster. However, their performance may be comprised 181 due to multiple local minima. Integrating characteristics from both local and global 182 methods,resultsinathirdcategory,thehybridmethods. Suchamethodisforexample 183 the hybrid genetic algorithm [26, 28]. 184 TheKrigingsurrogateusedlaterinthecasestudiesofthischapterisoftheordinary 185 trendtypeusingtheMatérnautocorrelationfunctionfamily,whilethehyperparameters 186 are estimated with the maximum-likelihood method and optimized using the hybrid 187 genetic algorithm method. 188 2.3 Polynomial chaos Kriging 189 PCK is another surrogate modeling method that combines features from PCE and 190 Kriging into one. In more detail, PCK is a Kriging model with a universal trend type. 191 Itstrendthoughisformulatedfromsparseorthogonalpolynomials. TheoutputofPCK 192 combines Eq. (4) and Eq. (12) and follows [29, 30, 31]: 193 MPCK(x) =/summationdisplay α∈AyαΨα(X) +σ2Z(x,ω). (22) The ﬁrst term,/summationtext α∈AyαΨα(X), is the trend of the Kriging model as described in 194 Section 2.2 and is computed from sparse PCE as described in Section 2.1. The second 195 term of Eq. (22) corresponds to the variance of the Gaussian process σ2andZ(x,ω) 196 is a stationary Gaussian process with zero mean and unit variance as of Eq. (12). 197 Therefore,buildingaPCKinvolvesastepofconstructingthePCEandanadditional 198 one to determine the Kriging model, i.e. to estimate its parameters. Two ways of 199 combining the above two steps are the so-called sequential and optimal approaches 200 [29, 30, 32]. 201 The PCK surrogate used later in the case studies of this chapter is developed using 202 the optimal approach and the Matérn autocorrelation function family is used along 203 withthemaximum-likelihoodestimationandthehybridgeneticalgorithmoptimization 204 method. 205 2.4 Leave-one-out cross-validation error 206 To assess the performance of the developed surrogates the leave-one-out (LOO) cross- 207 validation error is used, deﬁned by: 208 εLOO =1 NN/summationdisplay i=1/parenleftbigg M/parenleftbig x(i)/parenrightbig −Ms\i/parenleftbig x(i)/parenrightbig/parenrightbigg2 , (23) whereMs\i/parenleftbig x(i)/parenrightbig correspondstothesurrogatesdevelopedfromthesubsetoftheED, 209 acquired by removing the x(i)point. 210 8 DRAFT3 Global sensitivity analysis with Sobol’ indices 211 Broadly speaking, GSA aims to quantitatively identify the level each input variable 212 inX=/braceleftbig x(1),...,x(N)/bracerightbigT,X∈DX⊂RNaﬀects the response of Y=M(X), 213 Y∈R. As mentioned previously, in this dissertation, the Sobol’ indices are used for 214 GSA. It should be noted that this method is only valid for independent input variables. 215 3.1 Sobol’-Hoeﬀding decomposition 216 ForX∈DX⊂RNwithindependentcomponents, Y=M(X)withVar [Y]<+∞ 217 can be rewritten as [33, 11]: 218 Ms(X) =M0+N/summationdisplay i=1Mi(x(i)) +/summationdisplay 1≤i≤j≤NMij(x(i),x(j)) +···+M1,2,...,N(X) =M0+/summationdisplay u⊂{1,...,N} u/negationslash=∅Mu/parenleftbig Xu/parenrightbig , (24) whereM0is constant and denotes the mean value of Y,u/negationslash=∅,u⊂{1,...,N}and 219 Xuis a sub-vector of X, including the u-indexed elements. The summands in the 220 latter equation corresponds to 2N−1terms. The elementary functions are deﬁned by 221 conditional expectations: 222 Mu(Xu) =/summationdisplay v⊂u(−1)|u|−|v|E[Y|Xv=xv], (25) where|u|gives the cardinality of u. When the following holds: 223 /integraldisplay DXmMu/parenleftbig Xu/parenrightbig fXm(xm)dxm= 0,ifm∈u, (26) theSobol’-HoeﬀdingdecompositionofEq.(24)isunique. Duetothisuniqueness,the 224 orthogonality property: 225 E/bracketleftbig Mu(Xu)·Mv(Xv)/bracketrightbig = 0,ifu/negationslash=v, (27) also holds. Because of Eq. (26) and (27) the total variance of M(X)decomposes as: 226 D= Var [M(X)] =/summationdisplay u⊂{1,...,N} u/negationslash=∅Du, (28) where 227 Du= Var [Mu(Xu)] =E/bracketleftbig M2 u(Xu)/bracketrightbig , (29) corresponds to the partial variance. 228 9 DRAFT3.2 Sobol’ indices 229 The Sobol’ index Suis determined as the ratio between the partial and total variance 230 that represents the u-indexed set of input variables and follows: 231 Su=Du D. (30) For|u|= 1,theseparateimpactofeachinputvariable x(i)totheresponse Y=M(X) 232 is described by the indices respective to this speciﬁc variable, namely the ﬁrst-order 233 indicesS(1) i=Di/D. The impact originating from the interaction of pairs of the 234 input variables (x(i),x(j))and not already accounted in the Si,Sjis described by 235 the second-order indices S(2) ij=Dij/D. The impact from interaction of a larger 236 set of input variables is described from the higher-order indices. Total indices STi237 describe the overall impact of each input variable, considering its ﬁrst-order index and 238 all interactions with the other input variables and admit: 239 STi= 1−S−i, (31) whereS−idenotes the sum of all Suofubuti. It follows that/summationtext u/negationslash=∅Su= 1. 240 3.3 Sobol’ indices estimation 241 3.3.1 Monte Carlo-based estimation 242 The variances of Eq (28) and (29) can be obtained by estimators from Monte Carlo 243 simulations. These read [34]: 244 ˆM0=1 NN/summationdisplay i=1Ms/parenleftbig x(i)/parenrightbig ˆD=1 NN/summationdisplay i=1Ms/parenleftbig x(i)/parenrightbig2−ˆM2 0 ˆDu=1 NN/summationdisplay i=1Ms/parenleftbig x(i) n,x(i) ∼n/parenrightbig Ms/parenleftbig x(i) n,x/prime(i) ∼n/parenrightbig −ˆM2 0,(32) wherex(i) ∼ndenotes thei-th realization of xthat does not involve the ninput variable 245 andx/primecorresponds to a realization of Xwhich is independent of x=/parenleftbig x(i) n,x(i) ∼n/parenrightbigT. 246 Once the total and partial variance estimators are computed, the Sobol’ indices can be 247 derived from Eq. (30) and (31). 248 3.3.2 Sobol’ indices from polynomial chaos expansion 249 Recall from Eq. (4), that the PCE of M(X)admits: 250 MPCE(X) =/summationdisplay α∈AyαΨα(X). (33) 10 DRAFTDuetotheorthonormalityofthePCEbasis,themeanandvarianceof MPCE(X)can 251 be computed analytically from the ycoeﬃcients at no extra cost and follow [15]: 252 ˆµ=E/bracketleftbig MPCE(X)/bracketrightbig =E/bracketleftBigg/summationdisplay α∈AyαΨα(X)/bracketrightBigg =y0 (34) ˆσ2= Var/bracketleftbig MPCE(X)/bracketrightbig =E/bracketleftBigg/parenleftbigg MPCE(X)−y0/parenrightbigg2/bracketrightBigg =/summationdisplay α∈A α/negationslash=0y2 α.(35) The above equations hold since Ψ0≡1. Therefore, for the case that a PCE of Yis 253 already constructed, the Sobol’-Hoeﬀding decomposition of Eq. (24) can be rewritten 254 as: 255 MPCE(X) =y0+/summationdisplay u⊂{1,...,N} u/negationslash=∅MPCE u(Xu), (36) whereAu={α∈A:αm/negationslash= 0if and only if m∈u}corresponds to the set of multi- 256 indices including only uandMPCE u(Xu) =/summationtext α∈AuyαΨα(Xu), as of Eq. (33). 257 Thus,becausetheSobol’-Hoeﬀdingdecompositionisunique,thereexistsananalytical 258 representation ofMPCE u[14]. Accordingly to Eqs. (34), (35), the total and partial 259 variances ofMPCE(X)admit: 260 D= Var/bracketleftbig MPCE(X)/bracketrightbig =/summationdisplay α∈Ay2 α Du= Var/bracketleftbig MPCE u(X)/bracketrightbig =/summationdisplay α∈Auy2 α.(37) Hence, the ﬁrst-order and total Sobol’ indices are obtained from: 261 S(1) i=/summationtext α∈Aiy2 α/summationtext α∈Ay2α,Ai=/braceleftbig α∈A:αi>0, αi/negationslash=j= 0/bracerightbig STi=/summationtext α∈AT iy2 α/summationtext α∈Ay2α,AT i=/braceleftbig α∈A:αi>0/bracerightbig ,(38) respectively. Therefore, the Sobol’ indices can be obtained as a by-product of the PCE 262 coeﬃcients at no extra cost [14]. 263 4 Case study 264 This case study employs a virtual hybrid model to conduct the respective RTHS. 265 Therefore, a virtual PS (vPS) is utilized instead of physical specimen and the overall 266 RTHS is performed numerically, resulting thus in a virtual RTHS (vRTHS). 267 11 DRAFTThe hybrid model represents a prototype motorcycle. It consists of one vPS, the 268 electronically-controlledcontinuouslyvariabletransmission(eCVT)ofthemotorcycle, 269 andfourNS,namelyi)theengine,ii)themotorcyclebodydynamics,iii)therearwheel 270 braking system and iv) the front wheel braking system. Figure 1 depicts how the NS 271 and vPS are interconnected, forming the hybrid model. 272 hybrid_model_GSA.jpg Figure 1. Motorcycle hybrid model block diagram. The eCVT vPS is a multiple-input-multiple-output (MIMO) model including two 273 sets of one input and one output each. The ﬁrst set is interconnected to the motorcycle 274 engine NS, which simulates the dynamics of the combustion engine, while the second 275 settothemotorcyclebodyNS.Thelatterconnectionrepresentsthemotorcycle’strans- 276 missionoutputshaft. ThemotorcyclebodyNSincludesthemotorcycle’sbody/chassis 277 dynamicsalongwiththedynamicsofthewheels,tiresandsuspensions. Theproﬁleof 278 the road as well as the environmental driving conditions are also included in this NS. 279 TheengineNSisexpressedbyamultiple-input-single-output(MISO)model. Itsinputs 280 12 DRAFTare the throttle percentage thrand the angular velocity of the engine ωen, whereas 281 its output is the torque of the engine τen. The motorcycle body NS is expressed by 282 a MIMO model, which includes three sets of one input and one output each. The 283 ﬁrst set of input/output is the torque τvdand angular velocity ωvdof the transmission 284 output shaft, respectively. The last two sets are connected to the rear and front wheel 285 braking systems NS, respectively. The latter NS are both MISO models. In detail, 286 the inputs of the rear wheel braking system NS are the angular velocity of the rear 287 wheelωrwand the applied force on the brake pedal Fbrrw, whereas its output is the 288 brakingtorqueoftherearwheel τrw. Accordingly,theinputsofthefrontwheelbraking 289 system NS are the angular velocity of the front wheel ωfwand the applied force on 290 the brake lever Fbrfw, whereas its output is the braking torque of the front wheel τfw. 291 The aforementioned substructures are developed using the multi-physics simulation 292 software Simcenter Amesim . The report of [35] provides a detailed description of the 293 development of these substructures, explaining the equations that govern their motion. 294 Figure 2 depicts the real eCVT PS, located at the testing facilities of Siemens Industry 295 Software, that would be utilized in a non-virtual RTHS. 296 To examine the performance of the motorcycle prototype, testing of its virtual 297 hybrid model under predeﬁned driving, road and wind scenarios takes place. To co- 298 simulate the employed substructures and to coordinate the overall RTHS algorithm, 299 theSimcenter real-time platform was utilized. For the time integration scheme of the 300 RTHS, the fourth-order Runge–Kutta (RK4) method is used with a ﬁxed time-step of 301 0.1 msec. 302 Thesimulationdurationofthecasestudyis45seclong,onagivendrivingscenario 303 andwind/roadconditions. Inparticular,thedrivingscenarioisdeﬁnedbythevariations 304 of the throttle thrand of the forces applied to the brakes. It is assumed that the force 305 appliedinthebrakeleverisequaltotheforceappliedinthebrakepedalofthemotorcycle 306 ineachtime-stepoftheRTHS,namely Fbrrw=Fbrfw. Eqs.(39)and(40)describethe 307 values that the throttle and braking forces obtain, in each time-step of the simulation, 308 whileFigure3illustratestheirvariationsintime. Itshouldbenotedthatthemaximum 309 appliedthrottleis 0.5(50%),whichcorrespondstoahalf-openthrottle,andthebraking 310 forcesareexpressedinNewton. Eq.(41)describestheroadproﬁle,namelytheheightof 311 theground,thatisassumedinthiscasestudy. Itisexpressedbythesinusoidalfunction 312 h(x), wherexdenotes the current position of the motorcycle in meters. Additionally, 313 the ambient wind velocity is considered to be zero. 314 thr(t) =  0.125t,5≤t<9 0.5 ,9≤t≤13 −0.125t,13<t≤17 0 , elsewhere(39) 13 DRAFTFbrrw(t) =  10t,20≤t<25 50,25≤t≤32 −10t,32<t≤37 0, elsewhere(40) h(x) =/braceleftBigg 0 ,0≤x≤2 0.02/parenleftBig cos/parenleftbig2πx 3/parenrightbig −1/parenrightBig , elsewhere(41) Figure 2. The eCVT test bench at the testing facilities of Siemens Industry Software. Figure 3. Case study driving scenario. For the GSA, two hybrid model response QoI are selected, namely the maximum 315 vmaxandmeanvmeanvaluesofthemotorcyclevelocity. BothQoIareexpressedinkm/h. 316 Scalar quantities from the dynamic response of the motorcycle’s velocity were chosen 317 in particular, since the latter is a global response parameter that can characterize the 318 overalldynamicbehaviorofthehybridmodeloftheprototypemotorcycle. Inaddition, 319 the input variables to GSA are the 11 parameters which are listed in Table 1. The 320 probabilitydistributionsthatareassignedtoeachparameteralongwiththeirrespective 321 moments are also presented in Table 1. The bounds of each parameter were identiﬁed 322 from[36,37,38]toreﬂectarangeofpossibleparametervariationsofthecorresponding 323 motorcycle components. According to the notation presented in Sections 2 and 3, the 324 experimental design (ED) vectors for this case study follow: 325 X={Krt,Zrt,Kft,Zft,Krs,Zrs,Kfs,Zfs,M,J,µ} Y={vmax,vmean},(42) 14 DRAFTFromtheinputparametersofTable1,50pointsweregeneratedusingtheLatinhy- 326 percubesampling(LHS)method[39]andaccordingly50vRTHSwereconducted,col- 327 lectingtheQoIvaluesinthe YvectorofEq.42. WiththeobtaineddatainED {X,Y}, 328 a PCE, a Kriging and a PCK surrogate was developed for each QoI, as described in 329 Sections2.1,2.2and2.3,respectively. Figure4illustratesasampleRTHStimehistory 330 responseofthehybridmodel,indicatingtheQoIofthecasestudy,whicharehighlighted 331 in red color. Additionally, Figure 5 depicts some indicative dynamic responses of the 332 variables that couple the substructures in the virtual hybrid model. The responses in 333 the latter two ﬁgures are obtained using the mean values of Table 1, while the hybrid 334 modelisexcitedwiththedrivingscenarioandtheroadproﬁleandwindconditionsthat 335 are described by Eqs. (39)-(41). 336 Param.Prob. Distrib.Mean ValueStand. Dev.CV (%)Parameter DescriptionUnits KrtUnif. 58570 11714 20Vertical stiﬀness rear tireN m ZrtUnif. 11650 3495 30Vertical damping rear tireNs m KftUnif. 25000 5000 20Vertical stiﬀness front tireN m ZftUnif. 2134 640.230Vertical damping front tireNs m KrsUnif. 125000 25000 20Stiﬀness rear suspensionN m ZrsUnif. 10000 3000 30Damping rear suspensionNs m KfsUnif. 19000 3800 20Stiﬀness front suspensionN m ZfsUnif. 1250 375 30Damping front suspensionNs m M Unif. 300 15 5Motorcycle massKg J Unif. 0.0115 0.0023 20Engine moment of inertiaKgm2 µ Unif. 0.0012 0.00018 15Engine coeﬃcient of viscous frictionNm rev/min Table1.InputvariablesofthehybridmodelforGSAandtheirprobabilitydistributions. In addition, the time history responses of Figures 4 and 5 can be intuitively inter- 337 preted. Figure5bshowstheresponseofthefrontwheelbrakingtorque τfw,generated 338 by the utilized driving scenario. Recall from Figure 3 and Eq. (40), that the brakes of 339 themotorcycleareenabledinthetimeintervalbetweenthe 20thand37thsecondofthe 340 simulation and as a result the braking torques are non-zero only then. According to 341 Figure 3 and Eq. (39), the throttle applied in the motorcycle’s engine in the ﬁve ﬁrst 342 seconds of the simulation is zero and therefore its velocity is also zero (Figure 4). In 343 15 DRAFT Figure 4. Sample RTHS time history response of the motorcycle prototype hybrid model, using the mean values of Table 1. The QoI are highlighted in red color. thetimeintervalbetweenthe 5thand17thsecondofthesimulation,thethrottlestartsto 344 increases and hence the motorcycle starts to accelerate. From Figure 4 it can be seen 345 that the velocity is also increasing. After the 20thsecond of the simulation the braking 346 is starting (Eq. (40)) and as a result, the velocity of the motorcycle begins to decrease, 347 while it is brought to almost a full stop in the 45thsecond of the simulation. 348 5 Results and discussions 349 Figure 6 displays the convergence of the LOO errors of the PCE, PCK and Kriging 350 surrogates for both QoI. The LOO errors are computed using EDs of increasing size 351 from2to50samples. FromFigure6,itcanbeappreciatedthattheLOOerrorsofPCE 352 andPCKarequiteclose,whiletheerrorofKrigingisslightlylarger. Yet,alltheerrors 353 arenegligiblysmallastheyconvergetovaluesinthemagnitudeof 10−3orsmaller. The 354 satisfying performance of the surrogates can be also conﬁrmed from Figure 7, which 355 illustrates validation scatter plots by comparing the measurements of the two QoI with 356 the respective predictions from the corresponding PCE, PCK and Kriging surrogates. 357 Inmoredetail,foreachQoI,thesurrogatesaretrainedon40samplesoftheEDandthe 358 remaining 10 samples are used to validate them. 359 Inaddition,Figures9and10depicttheconvergenceplotsofthePCE-basedmoments 360 estimates, i.e. mean, standard deviation and CV values, for both QoI. Evidently, each 361 PCE-based moment estimate is converging in a stable manner for EDs larger than 20 362 samples. RecallthatthePCE-basedmomentsestimatescanbecomputedfromthePCE 363 coeﬃcientsaspost-processatnoadditionalcomputationalcost[14]. However,thelatter 364 is not a feature of Kriging or PCK. Thus, Figures 9 and 10 do not include moments 365 estimates from these two surrogates. However, to further evaluate the performance 366 of the developed surrogates, Figure 8 compares the prediction PDF of each surrogate 367 16 DRAFTagainstthehistogramsofthemeasurementsforeachQoI.Thesurrogatepredictionsare 368 evaluated based on 50 random samples, generated from the distributions described in 369 Table 1 using the LHS method. From Figure 8 it can be acknowledged that all three 370 developed surrogates can capture quite well both the mean value and the variance of 371 the corresponding measurements. 372 From Figures 6, 7 and 8 it can be appreciated that the performance and accuracy 373 prediction of PCE, PCK and Kriging are similar. However, it should be noted that 374 the computational time associated to train each surrogate diﬀers. For the case of 375 PCE and Kriging, the training computational time was approximately 1.3and1.6sec, 376 respectively. Nevertheless, for the case of PCK, it was approximately 12.7sec. The 377 latter ﬁnding can be intuitively interpreted as in the PCK case, the training phase 378 consists of two steps; one constructing the PCE and a second one to determine the 379 Kriging model. Note that these timings depend heavily on the ED sample size as well 380 as on the settings each surrogate is developed. 381 OncethesurrogatesofeachQoIaredeveloped,theSobol’indicescanbecomputed. 382 ForthecaseofPCE,theSobol’indicesarecomputedanalyticallyatnoextracost[14], 383 as described in Section 3.3.2. For the case of PCK and Kriging though, this is not 384 possibleandthustheSobol’indicesarecomputedusingMonteCarloevaluationsofthe 385 surrogateresponseQoI,asdescribedinSection3.3.1. Itshouldbenotedthatmorethan 386 104Monte Carlo evaluations of the Kriging and PCK surrogates are needed to ensure 387 convergence of the Sobol’ indices. In this case study, 105evaluations of the response 388 QoI of the Kriging and PCK surrogates were performed to compute the respective 389 Sobol’indices. However,itisworthmentioningthatthecomputationalcostofrunning 390 these evaluations was negligibly small ( ≈1sec for each surrogate). 391 In Figure 11 the ﬁrst-order and total Sobol’ indices for both QoI vmaxandvmean392 are reported. Clearly, both ﬁrst-order and total Sobol’ indices computed utilizing the 393 diﬀerent examined surrogates are admittedly close with each other. This conﬁrms 394 that the indices obtained from the PCK and Kriging converge to the analytical values 395 computed from the PCE coeﬃcients. From Figure 11 it can be stated that the input 396 variables related to the motorcycle mass Mas well as to the engine coeﬃcient of 397 viscous friction µ, are the most sensitive variables for both QoI. In particular, for the 398 QoI related to the maximum motorcycle velocity, namely vmax, the motorcycle mass 399 hasagreatereﬀectthantheenginefrictioncoeﬃcient. Whereas,fortheQoIrelatedto 400 the mean motorcycle velocity, namely vmean, the engine friction coeﬃcient contributes 401 morethanthemotorcyclemass. Onthecontrary,theremaininginputvariables,namely 402 thestiﬀnessanddampingofthefrontandreartireandsuspensionaswellastheengine 403 momentofinertia,donotaﬀecttheselectedQoI.Furthermore,sincetheﬁrst-orderand 404 total Sobol’ indices for both QoI are quite close, it can be noticed that the eﬀect from 405 high-order interactions between the input variables is negligible. 406 Finally, as mentioned before, in this case study surrogates with negligible small 407 validationerrorweredevelopedwithEDsizelargerthan20samples. However,itshould 408 benotedthattheinvestigatedcasestudywasnumericalandhencenomeasurementnoise 409 waspresentintheresults. Asreportedin[11],forthecasethatlargemeasurementnoise 410 is present in the obtained data, it is possible that ED of larger size would be needed to 411 train surrogates with an accepted validation error. 412 The development and implementation of the surrogate modelling, as well as the 413 17 DRAFTGSA, was performed with the UQLab software framework developed by the Chair of 414 Risk, Safety and Uncertainty Quantiﬁcation in ETH Zurich [40]. 415 6 Conclusions 416 Thispaperservesasasurveyofinvestigatinghowdiﬀerentsurrogatemodelingmethods 417 can be used to perform global sensitivity analysis via Sobol’ indices in hybrid simu- 418 lation. Global sensitivity analysis using Sobol’ indices is established as a powerful 419 tool to unveil the sensitivity of a structural system subjected to parameter and loading 420 variations. However,estimationoftheSobol’indicesthroughMonteCarlosimulations 421 oftheoriginalhybridmodelisrarelyaﬀordable,sinceseveralhybridmodelevaluations 422 are essential, while each evaluation relies on a single hybrid simulation. To alleviate 423 this burden surrogate models of selected hybrid model response quantities of interest 424 can be constructed and used to estimate the Sobol’ indices at much lower cost. A case 425 study encompassing a virtual hybrid model is employed to investigate this approach 426 by using a variety of diﬀerent surrogate modeling techniques. In this regard, polyno- 427 mialchaosexpansion,Kriging(a.k.aGaussianprocess)andpolynomialchaosKriging 428 are used to develop surrogates of the selected quantities of interest. For the case of 429 polynomial chaos expansion, the Sobol’ indices are computed as a by-product of the 430 polynomial coeﬃcients and hence no additional computational cost is required. How- 431 ever, the latter is not a feature of Kriging and polynomial chaos Kriging. Therefore, 432 Monte Carlo simulations using these surrogates are performed to estimate the Sobol’ 433 indices. Performance comparisons between each developed surrogate indicated that 434 the respective validation errors were similar. Additionally, negligibly small deviations 435 were observed in the Sobol’ indices that were obtained using the diﬀerent surrogate 436 modeling techniques. Future work aims at adaptive sampling of the input parameter 437 space of the hybrid model to minimize the experimentation cost necessary to compute 438 an accurate surrogate model for global sensitivity analysis. 439 7 Data Availability Statement 440 All data and code that support the ﬁndings of this study are available from the corre- 441 sponding author upon reasonable request. 442 Acknowledgments 443 This project has received funding from the European Union’s Horizon 2020 research 444 and innovation programme under the Marie Skłodowska-Curie grant agreement No. 445 764547. Thesoleresponsibilityofthispublicationlieswiththeauthor(s). TheEuropean 446 Union is not responsible for any use that may be made of the information contained 447 herein. 448 18 DRAFT