Hybrid Simulation is a dynamic response simulation paradigm that merges physical experiments and com- putational models into a hybrid model. In earthquake engineering, it is used to investigate the response of structures to earthquake excitation. In the context of response to extreme loads, the structure, its boundary conditions, damping, and the ground motion excitation itself are all subjected to large parameter variability. However, in current seismic response testing practice, Hybrid Simulation campaigns rely on a few prototype structures with ﬁxed parameters subjected to one or two ground motions of diﬀerent intensity. While this approach eﬀectively reveals structural weaknesses, it does not reveal the sensitivity of structure’s response. This thus far missing information could support the planning of further experiments as well as drive modeling choices in subsequent analysis and evaluation phases of the structural design process. This paper describes a Global Sensitivity Analysis framework for Hybrid Simulation. This framework, based on Sobol’ sensitivity indices, is used to quantify the sensitivity of the response of a structure tested using the Hybrid Simulation approach due to the variability of the prototype structure and the excitation parameters. Polynomial Chaos Expansion is used to surrogate the hybrid model response. Thereafter, Sobol’ sensitivity indices are obtained as a by-product of polynomial coeﬃcients, entailing a reduced number of Hybrid Simulations compared to a crude Monte Carlo approach. An experimental veriﬁcation example highlights the excellent performance of Polynomial Chaos Expansion surrogates in terms of stable estimates of Sobol’ sensitivity indices in the presence of noise caused by random experimental errors. Keywords: Hybridsimulation; globalsensitivityanalysis; Sobol’indices; surrogatemodeling; polynomial chaos expansion. 1 Introduction Background and motivation Hybrid Simulation (HS) is a dynamic simulation paradigm that merges phys- ical experiments and computational models into a hybrid model. In mechanical and electrical engineering, such an approach is known as hardware-in-the-loop simulation. In earthquake engineering, hybrid simulation is used to investigate the behavior and eventual failure of a structure in response to extreme earthquake excitation. In detail, HS is conducted using a hybrid model, which combines physical and numerical substructures (NS and PS, respectively) that interact with each other in a real-time control loop, to simulate the time history response of a prototype structure subjected to a realistic excitation. The PS is tested in the laboratory using servo-controlled actuators equipped with force and displacement transducers. The NS is numerically simu- lated to avoid testing of those substructures that can be reliably replaced by numerical models (e.g., masses or components whose response remains in an easy-to-simulate range) or would exceed the testing capacities of the experimental facility. A real-time software ensures interface force balance and displacement compatibility between NS and PS by solving the equation of motion of the prototype structure online, that is, while struc- tural testing is performed. As a result, HS reproduces the time history response of the prototype structure subjected to a realistic excitation. When the PS response is rate-independent, it is good practice to perform HS with an extended time scale from 50 to 200 times slower than real-time in order to improve control accuracy. Time scaling and testing of geometrically scaled substructures also drastically reduce the costs of HS compared Corresponding Author: Dr. Giuseppe Abbiati, abbiati@cae.au.dk , ORCID: 0000-0002-5048-8505 . 1 DRAFTto full-scale experimentation. The seminal report by Schellenberg and Mahin [1] provides a comprehensive overview of state-of-art HS methods and algorithms. HS is suitable either to investigate the behavior of a prototype structure at the system level, where PSs replace components lacking validated numerical models, or to observe the component-level behavior of PSs ex- posed to a realistic excitation, which includes dynamic interaction with possibly multiple NSs. The overarching goal is to support the validation and calibration of computer models [2]. A system-level focus is pursued in earthquake engineering, where HS is preferred to shake table experiments for testing large structures such as bridges (e.g., [3]) and buildings (e.g., [4]). Similarly, testing of ﬂoating platforms with a sizeable scaling factor in hydrodynamic laboratories is challenging due to the extent of mooring line footprints, which can range from two to four times the water depth. For this reason, HS has been recently proposed to truncate mooring lines whose missing parts are simulated as NSs [5]. A component-level focus is pursued in ﬁre engineering, where hybrid ﬁre testing, the adaptation of HS to ﬁre loading, originated in response to limitations of standard ﬁre tests [6]. The latter consists of exposing a single structural component to a predeﬁned time-temperature loading curve inside a furnace and indeed does not account for the redistribution of internal forces caused, for example, by thermal expansion or partial collapse. For a similar motivation, Franza and Marshall [7] developed an HS platform to study soil-structure interaction problems where a numerical model of a building frame (i.e., the NS) is used to mimic a realistic structural loading exerted on the physical centrifuge model of a tunnel-soil-foundation system (i.e., the PS). The use of HS in other engineering and science disciplines is growing (e.g. [8], [9]). In earthquake engineering, which is the main ﬁeld of expertise of some of the authors, boundary conditions, damping, and ground motion are all subjected to large parameter variability [10]. However, in the current experimental earthquake engineering practice, HS campaigns rely on a few prototype structures with ﬁxed parameters subjected to one or two ground motions of diﬀerent intensity. This testing procedure eﬃciently identiﬁes the weaknesses of the tested prototype that could lead to the collapse of real structures in earth- quakes. However, such a truncated procedure is indeed not adequate to estimate the structural sensitivity of the prototype structure response, i.e. quantifying the fraction of the variability of a generic response quantity of interest (QoI) due to the variability of each of the input parameters. This thus far missing information could support the planning of further experiments as well as drive modeling choices in subsequent analysis and evaluation phases of the structural design process. Scope This paper describes a Global Sensitivity Analysis (GSA) framework for Hybrid Simulation that is intendedtoquantifythefractionofthevariabilityofagenericresponseQoIoftheprototypestructure, evaluated via HS, due to the variability of a set of input variables that provide a suitable parametrization of the prototype structure and the excitation. Among the many GSA options, Sobol’ sensitivity indices , a class of variance- decomposition-based methods [11], were selected. Methods of this class rely upon the decomposition of the variance of a QoI as a sum of contributions of each input variable taken singularly or in combination with others. Unlike regression-based methods, variance-based methods do not require any kind of linearity or monotonicity of the functional relationship between input variables and QoIs. An eﬃcient way to compute Sobol’ indices, originally introduced by Sudret [12, 13], is by post-processing the coeﬃcients of a Polynomial Chaos Expansion (PCE) of the QoI. PCE is a spectral decomposition method that projects the model response onto a functional basis of orthonormal multivariate polynomials in the input variables [14, 15]. Once the PCE representation is available, Sobol’ indices are readily available by simply post- processing polynomial coeﬃcients with no additional computational cost. The strategies for the calculation of the PCE can be broadly classiﬁed into two approaches: i) intrusive i.e., schemes that require the knowledge and adaptation of the governing equations of the model; and ii) non-intrusive i.e., schemes that (repeatedly) utilize the model to compute an outcome without modifying the model itself, which is considered as a black-box . HS and its PSs, and experimental testing in general, imply a non-intrusive PCE strategy. Hence, in this paper, the sparse adaptive least-square minimization based on Least Angle Regression (LAR) [16] was adopted. The proposed GSA framework for HS is illustrated using a case study based on a 3-degrees-of-freedom (3-DoFs) prototype structure excited by earthquake ground motions. HSs are performed using a 3-DoFs HS test rig installed at the Structural Laboratory of ETH Zurich. The GSA of the hybrid model response is performed for two response QoIs against four input variables, which parametrize both the NS and the excitation. While to the best knowledge of the authors this is a pioneering application of GSA to HS, it is important to remark that the goal of this study is not to quantify the generic sensitivity of the model response to uncertain inputs as in the context of uncertainty quantiﬁcation using Sobol’ sensitivity analysis. Instead, the proposed procedure aims to identify the primary drivers of the hybrid model response, that is, to better understand the inner workings of the hybrid model itself and thus inform the design of a HS campaign. For a discussion on 2 DRAFTsome recent applications of GSA outside structural engineering see e.g., [17, 18, 19]. The earthquake engineering validation case study was selected because this is the main applied ﬁeld of expertise of some of the authors and because of the available test setups. However, the proposed procedure treats the hybrid model as a black-box , and is thus applicable to any dynamic system response investigated via HS. Thepaperisorganizedasfollows. Section2describestheGSAframeworkbasedonPCE.Section3describes the case study, including a detailed description of the HS setup and the HSs. Section 4 presents the results of the GSA of the prototype structure response. Finally, conclusions and future outlook are given. 2 Global Sensitivity Analysis framework Let us denote with X=fX1;:::;XMganM-dimensional vector of input variables to parametrize the NSs, the PSs and the excitation (or a subset of these) of the hybrid model of the prototype structure, and with Y a generic response QoI of the prototype structure evaluated via HS of the response of the hybrid model. A mappingMbetween input variables and a single QoI reads, X2DXRM7!Y=M(X)R (1) It can be argued that the inherent variability of the response of diﬀerent samples corresponding to nominally identical PSs of the hybrid model, as well as measurement noise during HS, which both propagate through to the hybrid model response QoI, requires a stochastic mapping between input variables and QoIs. However, the scope of this paper is limited to the case of PSs with almost deterministic behavior. Moreover, relevant studies testify that measurement noise of standard structural testing equipment has a negligible eﬀect on the prototype structure response evaluated via HS [20, 21]. Therefore, Mis assumed to be deterministic. Since our target is to estimate the structural sensitivity of the prototype structure response, which means quantifying the fraction of the variability of a generic QoI due to the variability of the input variables, all input variables in Xare assumed to be statistically independent and uniformly distributed. The ﬁnite bounds of input variables compatible with force and displacement ranges allowed by the experimental facility are the only constraints of X. 2.1 Sobol’-Hoeﬀding decomposition and sensitivity indices Assume that the vector of input variables Xhas supportDXand follows an independent joint probability density function (PDF) fX(x) =QM ifXi(xi)wherefXiis the marginal PDF of the i-th input variable. Any square-integrable mapping Y=M(X)with respect to the probability measure associated with fX, can be written as a sum of functions of increasing dimension as [22]: M(X) =M0+MX i=1Mi(Xi) +X 1ijMMi;j(Xi;Xj) +:::+M1;2;:::;M (X); (2) or equivalently: M(X) =M0+X u6=;Mu(Xu); (3) whereM0is the mean value of Y,u=fi1;:::;isgf 1;:::;Mgare index sets, and Xudenotes a subvector of Xcontaining only the components indexed by u. The number of summands in the above equation is 2M1. The Sobol’ decomposition is unique under the condition: Z DXkMu(xu)fXk(xk)dxk= 0, ifk2u; (4) whereDXkandfXkdenote the support and marginal PDF of Xk, respectively, which leads to the orthogonality property: E[Mu(xu)Mv(xv)] = 0, ifu6=v: (5) The uniqueness and orthogonality properties allow for the following decomposition of the variance DofY: D=Var[M(X)] =X u6=;Du (6) 3 DRAFTwhereDudenotes the partial variance: Du=Var[Mu(xu)] =E M2 u(xu) : (7) The Sobol’ index Sucan be deﬁned as the fraction of the total variance Duthat corresponds to the set of input variables indexed by u: Su=Du D: (8) By construction,P u6=;Su= 1. First-order indices S(1) idescribe the inﬂuence of each parameters Xiconsidered separately. Second-order indices S(2) ijdescribe the inﬂuence from pairs of parameters fXi;Xjgnot already accounted for by XiorXjseparately. High-order indices describe combine inﬂuences from larger sets of parameters. The total sensitivity indices S(tot) irepresent the total eﬀect of an input variable Xiaccounting for its main eﬀect and all interactions with other input variables. It follows that S(tot) i= 1Si, whereSiis the sum of allSuwithunot including i. Sobol’ indices can be evaluated by Monte Carlo simulation [11], which requiresO(103)model evaluations for each index Su. Clearly, this is not aﬀordable when each evaluation of M(x)entails the execution of a HS. On the other hand, when the PCE of the QoI is available, Sobol’ indices can be obtained analytically at no additional cost. A concise description of PCE-based Sobol’ sensitivity indices estimation is given in the following section; for further details, the reader is referred to [12, 13]. 2.2 Sobol’ indices from Polynomial Chaos Expansion Sparse PCE [16] is a well-known surrogate modeling technique. Its close relation to Sobol’ variance decom- position [12] together with good denoising performance [23] make sparse PCE ideal to surrogate the proto- type structure response evaluated via HS. PCE relies on the decomposition of a random model response QoI Y=M(X)as a linear superposition of non-linear functions as follows: bY=MPCE(X) =X 2Ay (X); (9) wheref ; 2Agis a set of multivariate polynomials that are orthogonal with respect to the input vector with independent components XfX(x) =QM i=1fXi(xi), = ( 1;::: M)is a multi-index that identiﬁes the polynomial degree in each of the input variables, and y denotes the corresponding polynomial coeﬃcient (coordinate). Following the orthonormality condition, E MPCE(X) =y0. For practical purposes, the inﬁnite sum in (9) needs to be truncated to a ﬁnite series. This is commonly achieved by maximum-degree or hyperbolic norm truncation (for more details, see [16]). In this paper, a sparse least-square regression for the estimation of the set of coeﬃcients y=fy ; 2Agis adopted [16]. Aso-called experimental design (ED),consistingofasetofrealizationoftheinputvector X=fx(1);:::;x(N)g and the corresponding model evaluations Y=fy(1);:::;y(N)g, is generated ﬁrst using repeated HS. In a second step, the set of coeﬃcients y is estimated by minimizing the expected mean-square approximation error on the ED by solving: by = arg min y E YbY2 (10) In the present application, the ED was formed by sampling the input variable space with a Sobol’ low- discrepancy sequence. The minimization in (10) was solved using the hybrid LAR method originally proposed in [16]. A good measure of the accuracy of PCE is the mean-square residual error ErrG=E YbY2 , often simply referred to as generalization error . However, the evaluation of ErrGwould require an additional validation data set, which is typically unavailable. It is thus desirable to obtain an estimate of ErrGusing only the information included in the ED. One such error estimation method is the co-called leave-one-out cross-validation error (LOO). Denote by MPCEnithe PCE built from the subset of the ED obtained by removing the pair fx(i);y(i)g. The LOO error is deﬁned as: loo=1 NNX i=1h M x(i) MPCEni x(i)i2 : (11) 4 DRAFTIn the case of regression-based PCE, loocan be calculated directly from the coeﬃcients of the PCE computed with the entireED, without explicitly solving (10) multiple times [15]. This robust error estimate allows one to devise basis-adaptive PCE schemes. A truncation scheme can be gradually relaxed, e.g. by increasing the maximum polynomial degree, until some convergence criterion on loois achieved [16]. Due to the rapid convergence of the sparse PCE coeﬃcients, this approach provides a very eﬃcient means to propagate the variabilityintheinputvariablestoQoIswithouttheneedforseveralevaluationsofthemodelresponseotherwise required by standard MC techniques. It is straightforward to obtain the Sobol’ decomposition of Yin an analytical form by observing that the summandsMPCE u(Xu)in (3) can be written as: MPCE u(Xu) =X 2Auy (Xu); (12) whereAudenotes the set of multi-indices that depend only on u: Au=f 2A: k6= 0if and only if k2ug (13) Clearly,[Au=A. Consequently, due to the uniqueness of Sobol’-Hoeﬀding decomposition, there is an analyt- ical expression of MPCE u(12), which serves as a proxy of Mu. Due to the orthogonality of the basis, the total variance of a PCE is analytically given by [14]: D=Var MPCE(X) =X 2Ay2 : (14) Similarly, the partial variance Dureads: Du=Var MPCE u(Xu) =X 2Auy2 : (15) Accordingly, the Sobol’ indices of any order can be approximated by a simple combination of the squares of the PCE coeﬃcients by substituting (15) and (14) in (8). For instance, the ﬁrst-order Sobol’ indices, which describe the inﬂuence of each input variable Xiconsidered separately, read: S(1) i=P 2Aiy2 P 2Ay2 ,Ai=f 2A: i>0, i6=j= 0g (16) whereas the total Sobol’ indices, which represent the total eﬀect of an input variable Xiaccounting for its main eﬀect and all interaction with other input variables, are given by: S(tot) i=P 2Atot iy2 P 2Ay2 ,Atot i=f 2A: i>0g: (17) Sobol’ indices provide quantitative insight on the importance of an input variable. However, they do not include information about the direction in which an input variable aﬀects the QoI. So-called univariate eﬀects can answer this question [24, 18]. A univariate eﬀect is the expectation of a QoI conditioned to the value of a single input variable. Univariate eﬀects have a closed analytical form for PCE models, which is closely related to the ﬁrst-order Sobol’ decomposition: M(1) i(Xi) =X 2Aiy (Xi),Ai=f 2A: i>0; i6=j= 0g (18) In this paper, computations of both PCE and Sobol’ indices are performed using UQLab, which is a MATLAB toolbox for uncertainty quantiﬁcation developed by the Chair of Risk, Safety and Uncertainty Quantiﬁcation of ETH Zurich [25, 26, 27]. 3 Case study In order to provide experimental validation of the proposed GSA framework for HS, a case study was conducted using a 3-DoFs prototype structure subjected to earthquake ground motion excitation. 5 DRAFT3.1 Prototype structure and excitation The prototype structure consists of a simply-supported massive beam coupled at the ends to two linear-elastic rotation restraints, two lumped masses, and one dashpot. The beam is subjected to a ground motion excitation ag(t)acting in the direction orthogonal to its axis. Bending and axial deformation of the beam are coupled for large displacements, which means that the prototype structure is characterized by a geometrical nonlinearity. Figure 1 provides a schematic view of the prototype structure and its hybrid model substructured into one PS and two NS for the purpose of HS. (a) (b) Figure 1: Case study: (a) prototype structure; (b) hybrid model comprising two NS and one PS. In Figure 1 u1andu2indicate the two rotational DoFs of the simply-supported beam, while u3corresponds to its axial DoF. The PS is instantiated using a L= 470mm long iron plate (Young modulus E= 100GPa, density= 7850kg m3), with a 2002mm cross-section (area A= 400mm2and moment of inertia I= 133mm4). Since HSs are conducted with a testing time scale equal to 50, a numerical mass matrix accounts for the PS inertia. In detail, the PS mass matrix is estimated by condensing the six-by-six consistent mass matrix of a two-dimensional Bernoulli beam element to the three DoFs retained by the hybrid model. The out-of-diagonal entries of the reduced mass matrix couple the inertial acceleration of the ground motion to the rotation DoFs of the hybrid model. The NS of the hybrid model comprises two rotational masses J1=J2= 0:9kgm2, a translation mass M= 20;000kgall lumped at the supports. The linear dash-pot with a damping coeﬃcient C= 300;000Ns mis the only source of numerical damping and prevents the excitation of the axial eigenmode of the hybrid model. The two linear elastic rotation springs, whose stiﬀness K1andK2are deﬁned as input variables for GSA, represent the non-deterministic properties of the prototype structure supports. The velocity-pulse model proposed by [28], which provides an intuitive representation of near-fault earth- quake ground motions, deﬁnes the excitation ag(t). The velocity-pulse model consists of a cosine carrier wave with a period Tp, modulated by a truncated cosine function, namely the pulse function, with a period Tpand amplitudeVp. The latter corresponds to the peak velocity of the pulse occurring at time tmax;p. The expression for the velocity-pulse reads: vg(t) =Vp 2cos 2ttmax;p Tp + Dr Tp 1 + cos2 ttmax;p Tp ; (19) whereis the phase angle between the cosine wave and pulse function, while Dris the non-zero residual displacement at the end of the pulse, which is in turn given by: Dr=VpTpsin (+ )sin ( ) 4(1 2): (20) The transverse inertial acceleration ag(t)applied to the prototype structure corresponds to the time derivative of the velocity-pulse model. The velocity-pulse model is deﬁned for t2f0; Tpgand is zero outside this interval. 6 DRAFT (a) (b) Figure 2: Sample excitation used in the case study: (a) Velocity-pulse time history, and (b) the corresponding acceleration response spectrum Sacomputed using a 5% damping ratio. Table 1: Input variable intervals and assumed probability distributions. Parameter Probability distribution Lower bound Upper bound Unit K1 Uniform 20 120 [Nm rad] K2 Uniform 20 120 [Nm rad] Vp Uniform 3:50 4:50 [m s] Tp Uniform 0:50 1:50 [s] In order to reduce the number of free parameters, the velocity pulse peak is located at tmax;p = Tp 2and = 2 and= 0:1. Thus, pulse peak velocity VPand period TPare selected as input variables for GSA. Figure 2 depicts a realization of the velocity-pulse model described by (19) and (20) with VP= 2m s,TP= 1:25s, = 2:0 and= 0:1rad. The GSA of the hybrid model is evaluated using two response QoIs, u1;max, the maximum absolute rotation at the left support of the beam measured using the vertical actuator displacement, and r1;max, the maximum absolute bending moment at the left support of the beam measured using the vertical actuator load cells, as depicted in Figure 1. Following the notation introduced in Section 2, the vector of input variables and the vector of QoIs read: X=fVP;TP;K1;K2g Y =fu1;max;r1;maxg(21) Theinputvariablesuniformlydistributed. Inparticular, theintervalsreportedinTable1entailalinearresponse of the hybrid model, which is compatible with the force and displacement capacity of the HS test setup. IfK1=K2and both vary within the range proposed in Table 1, the period of the ﬁrst bending mode of the hybrid model ranges between 0:68and0:45s. This range overlaps with the period range of the carrier signal of the velocity-pulse excitation, which varies between 0:50and2:00s. Accordingly, the variable frequency tuning between the excitation and the ﬁrst bending eigenmode of the hybrid model is expected to aﬀect the dynamic ampliﬁcation of both response QoIs. There are two reasons for setting the domains of input variables to ensure a linear elastic PS behavior. Firstly, this choice allows for performing the entire experimental campaign using a single plate specimen as PS. Second, the case study focuses on the measurement noise as the sole source of uncertainty that aﬀects the intra-experiment variability of the PS response. Ultimately, the proposed GSA method for HS could be validated in this case study. 7 DRAFT Figure 3: Architecture of the 3-DoFs HS test rig. 3.2 Hybrid simulation setup The 3-DoFs HS test rig used to conduct the HSs is a stiﬀ loading frame equipped with four electro-mechanical actuators interfaced to an INDEL real-time computer [29]. The 3-DoFs HS test rig is designed to test plate specimens with an approximate footprint of 200500mm and thickness varying between 1and3mm. Figure 3 illustrates the architecture of the HS setup, including a close-up view of the plate specimen accommodation. The GINLink bus connects the actuator servo-driver INDEL SAC4 and all data acquisition modules INDEL COP-ADA to the real-time computer INDEL SAM4, which runs the HS software. The latter is developed in the MATLAB/SIMULINK environment and downloaded to the real-time computer INDEL SAM4 via Ethernet from the Host-PC. At each simulation time step, the HS software imposes displacements u1,u2andu3to the plate specimen, the PS, reads the corresponding restoring forces r1,r2andr3measured using force transducers, and solves the coupled equation of motion of the hybrid model. In addition, a laser sensor measures the out- of-plane deﬂection at the mid-span of the plate specimen (labeled UL in Figure 3). A detailed description of the time integration scheme adopted to solve the equation of motion is reported in the appendix of this paper. Figure 4 shows two axonometric views of the 3-DoFs HS test rig, including the main hardware components. In this ﬁgure, the gray-colored parts are ﬁxed to the reaction frame, which is not visible, while the moving parts of the 3-DoF HS test rig are colored in yellow. The two rack-pinion systems (10) convert motion along the vertical actuator axes y1andy2(1) to rotations u1andu2, respectively, which are imposed to the short edges of the plate specimen (6) through aluminum clamps (3). Horizontal actuators along axes x1andx2(2) control the position of the moving frame (4), which is mounted on two proﬁled rail guides using ball bearings (5). As a result, the position of the moving frame (4) equals the axial elongation of the plate specimen u3(6). 4 Case study results This section reports the results of the case study where the GSA procedure described in Section 2 was employed to evaluate and quantify the sensitivities of the response of the prototype structure deﬁned in Section 3 to earthquake ground motion excitation using the results of a number of HS conducted using the 3-DoFs test rig described in Section 3. 8 DRAFT Figure 4: Axonometric views of the 3-DoFs HS test rig with main components: (1) vertical actuators; (2) horizontal actuators; (3) installation clamps; (4) moving frame; (5) proﬁled rail guides; (6) plate specimen; (7) hinges; (8) vertical actuator load cells; (9) horizontal actuator load cells; (10) rack-pinion systems. The moving parts are colored in yellow while the grey parts are ﬁxed to the reaction frame, which is omitted in this ﬁgure for clarity. 4.1 Surrogate modeling using PCE Based on the input variable intervals and distributions speciﬁed in Table 1, a Sobol’ sequence was used to generate 200 samples of the input variable vector X. For each of these samples, a HS of the response of the hybrid model (Figure 1) of the prototype structure to a velocity-pulse ground motion excitation was conducted to determine the time history of the response QoIs in vector Y. An outcome of a sample HS of the response of the prototype structure with support rotation spring stiﬀness K1=K2= 70Nm/rad to the velocity pulse shown in Figure 2 in terms of time histories of the rotation u1(t)and moment r1(t)is reported in Figure 5, where the corresponding peak response QoIs u1;maxandr1;maxare also highlighted. According to the procedure outlined in Section 2.2, ED was used to train the PCEs of response QoI u1;maxandr1;max. The time history of support moment r1is much noisier than that of support rotation u1. In principle, when the stiﬀness of the specimen is much smaller than the stiﬀness of the actuation system, displacement control is fairly accurate, and restoring force measurement noise is the dominant component of experimental errors, which directly propagates to the acceleration response history of the hybrid model. Double integration from acceleration to displacement, done by the dynamic equation of motion of the hybrid model, acts as a low-pass ﬁlter that reduces the impact of noise on displacement measurements. As a consequence, measurement signal noise aﬀects r1;maxmore than it does u1;maxand, further negatively aﬀects the convergence of the corresponding PCE. Validation scatter plots of Figure 6 compare the measured QoIs to the corresponding PCE predictions and conﬁrm this hypothesis. In detail, for each QoIs, four PCEs were trained on four disjoint 150-sample partitions of the ED and evaluated on the remaining 50 samples. Figure 7 shows the convergence of LOO errors of both PCEs computed considering EDs of increasing size from 10 to 200 samples. It is worth noting that force measurement noise and displacement control errors are both sources of aleatory uncertainty that propagate through the hybrid model response, which is inherently stochastic. In this case, PCE acts as a data denoiser, and for this reason, LOO errors do not converge to zero. Figures 8 and 9 show the convergence plots of the mean, variance and CoV PCE-based estimates for both QoIs. Evidently, all PCE estimates converge in a stable manner, with reliable values attained for ED sizes larger than 50 samples. 9 DRAFT (a) (b) Figure 5: Outcomes of a sample HS of the response of the hybrid model of the prototype structure with K1=K2= 70Nm/rad to a velocity-pulse ground motion excitation with VP= 2m/s,TP= 1:25s in terms of left support rotation (a) and bending moment (b) response histories. QoIs are highlighted. (a) (b) Figure 6: Comparison between QoI measurements and corresponding PCE estimates for (a) u1;maxand (b) r1;max. 4.2 GSA based on Sobol’ sensitivity indices Based on PCEs of both QoIs, Sobol’ indices are obtained as explained in Subsection 2.2. Figure 10 reports ﬁrst- order and total Sobol’ indices for QoIs u1;maxandr1;max. The error bars of Figure 10 refer to 95% conﬁdence intervalsandarecomputedfrom50bootstrapreplicationsoftheunderlyingPCEbyusingthepercentilemethod introduced in [30]. It is important to remark that there is no theoretical guarantee that the percentile-based method provides unbiased quantile estimates. However, as discussed and demonstrated in [30], bias decays quickly with PCE, and, for a large enough ED (30-50 points), bootstrap provides an eﬀective estimator of the accuracy of the Sobol’ indices. According to Figure 10, the period Tpof the velocity pulse and the stiﬀness K1of the left torsional spring are the most sensitive input variables for both QoIs. This result ﬁnds an intuitive explanation. The range of velocity-pulse period Tpwas selected to cross the range of periods of the ﬁrst bending mode of the hybrid 10 DRAFT (a) (b) Figure 7: Convergence of LOO error estimates of (a) u1;maxand (b)r1;maxPCEs. (a) (b) (c) Figure 8: Convergence of PCE estimates of (a) mean, (b) variance and (c) CoV of QoI u1;max. (a) (b) (c) Figure 9: Convergence of PCE estimates of (a) mean, (b) variance and (c) CoV of QOI r1;max. model, given the range of the stiﬀness of the rotation springs. Further, rotation spring stiﬀness K1is a boundary condition at the left edge of the PS iron plate, corresponding to responses u1andr1. Therefore, QoIs are expected to be more sensitive to K1thanK2. It should be noted that the bending stiﬀness of the plate specimen has the same order of magnitude as the maximum stiﬀness value of the left-support rotation springK1(see Table 1), thus a ﬁxed boundary conﬁguration, which would heavily aﬀect the response of the plate, was not approached in this case study. It is, indeed, reasonable from the standpoint of resonance that the velocity-pulse period TPhas more inﬂuence on the prototype structure response than K1, and signiﬁcantly more inﬂuence than the velocity pulse amplitude VP. Finally, the small diﬀerence between pairs of ﬁrst-order and total Sobol’ indices indicates that the high-order interactions among input variables are negligible. In order to provide further insight into the prototype structure response, univariate eﬀects are reported in 11 DRAFT (a) (b) (c) (d) Figure 10: First-order and total Sobol’ indices for (a-b) u1;maxand (c-d)r1;maxPCEs. Figure 11 for input variables TPandK1, to which the response of the prototype structure is the most sensitive, as highlighted by the GSA. Related 95% conﬁdence intervals computed with 50 bootstrap replications of the underlying PCEs are also reported in the same ﬁgures. Evidently, univariate eﬀects decrease monotonically for PCEs of both QoIs over the ranges of the two considered input variables. 5 Conclusions A Global Sensitivity Analysis (GSA) framework for Hybrid Simulation (HS) was presented in this paper. The proposed framework, based on Sobol’ sensitivity indices, is intended to quantify the fraction of the variability of a generic response quantity of interest (QoI) of the prototype structure due to the variability of a set of input variables that provide a suitable parametrization of the hybrid model and the excitation. The response of the prototype structure is evaluated using HS, thus involving a hybrid model comprised of a few Physical Substructures (PSs) tested in laboratory test setups, and a potentially large number of Numerical Substructures (NSs), emulated using computer models. Such an approach provides for a more realistic evaluation of the response of the prototype structure but still precludes the conduct of a large number of tests due to cost and time considerations. To resolve this issue, Polynomial Chaos Expansion (PCE) was used to surrogate the hybrid model response. Thereafter, Sobol’ sensitivity indices were obtained as a by-product of polynomial coeﬃcients, entailing a reduced number of HSs compared to a crude Monte Carlo approach to achieve the same quality of response sensitivity estimates. A case study characterized by four uncertain input variables, two of which related to the NS while the other two related to the excitation, was used to validate the proposed framework. As a general observation, the noise ﬂoor of each response QoI of the prototype structure, which is determined by random experimental errors entailed in HS, prevents the generalization error of the corresponding PCE to converge to zero asymptotically. In the case of a soft specimen, such as the one used in the case study, 12 DRAFT (a) (b) (c) (d) Figure 11: Univariate eﬀects on PCEs for (a-b) u1;maxand (c-d)r1;max. Dashed lines indicate 95% conﬁdence intervals. displacement-based QoIs are less noisy than force-based QoIs. However, stable plateau values indicate that PCE acts as a powerful denoiser and, therefore, is very well-suited to surrogate the hybrid model response, making the proposed GSA framework for HS possible. Stable PCEs and, hence, stable Sobol’ sensitivity indices, were available with as few as 50 HSs in the conducted case study. Several such HSs can be easily aﬀorded when the PSs of the hybrid model can sustain repeated experiments without deteriorating appreciably (e.g. seismic isolation bearings). However, when the PSs accumulate damage during HS and need to be replaced experiment after experiment, several such HSs might not be aﬀordable. Therefore, future work should aim to reduce the number of HSs necessary to achieve stable convergence of PCEs by, for example, employing adaptive sampling strategies to build the ED. GSA adds a further dimension to HS to make it into a tool for exploring the dynamic response of the prototype structure in the domain of engineering relevance of the prototype structure response. Namely, a functional understanding of the hybrid model response to a dynamic excitation is sought (e.g. what inﬂuences what and how much?) within a domain deﬁned by reasonable variation of the input parameters. Furthermore, GSA enables quantiﬁcation of the impact of experimental uncertainties on surrogate models, a technique not documented in the literature to the best knowledge of the authors, opening a new ﬁeld of stochastic HS in structural earthquake engineering. Finally, the ﬁndings of the presented application example indicate that experimental campaigns designed according to the current earthquake engineering practice and, therefore, relying on a few HSs with degrading PSs, can only provide qualitative instead of quantitative information about the structural sensitivity of the prototype structure response. Yet, quantifying the response sensitivity is very valuable, both in terms of planning additional HSs for model validation and veriﬁcation, and in terms of modeling and analysis choices in the overall structural design process. This motivates our continuing work. 13 DRAFTAppendix A. Time integration algorithm for hybrid simulation A partitioned time integration algorithm tailored to state-space equations of motion computes the time history response of the hybrid model in a HS. The state-space formulation of the equation of motion of a generic mechanical system reads: M_Y+R(Y) =F(t) (A.1) with, Y=2 4u v s3 5;M=2 4I 0 0 0 m 0 0 0 I3 5;R=2 4v r(u;v;s) g(u;v;s)3 5;F(t) =2 40 f(t) 03 5 (A.2) where vectors u,vandscontain displacement, velocity and additional state variables, respectively. The former two vectors always appear as a pair in second-order mechanical systems, while the latter is used to model nonlinearities endowed with memory, e.g. hysteresis. In particular, r(u;v;s)is the nonlinear restoring force vector while g(u;v;s)is a nonlinear function that describes the evolution of additional state variables s;m is the mass matrix and f(t)is the external time varying load, while Iand0are identity and zero matrices, respectively. When (A.1)-(A.2) refer to a linear system, no additional state variable senters the state vector Yand the restoring force is deﬁned as: r(u;v) =ku+cv (A.3) where kandcare stiﬀness and damping matrices, respectively. For the sake of simplicity, time dependency is omitted hereafter while M,K,RandFare referred to as generalized mass, stiﬀness, restoring force and external loading. The partitioned time integration algorithm, which adopts the dual assembly strategy of the modiﬁed Prakash-Hjelmstad method [31], synchronizes two monolithic time integration processes characterized by diﬀer- ent time steps. In detail, a set of Lagrange multipliers enforces compatibility of PS and NS interface velocities at the coarser time step. The corresponding coupled system of equations of motion reads: 8 >>>< >>>:MN_YN k+1+RN YN k+1 =FN k+1+BNTk+1 MP_YP k+j m+RP YP k+j m =FP k+j m+BPT mj mk+j mk+1 GN_YN k+1+GP_YP k+1=0(A.4) where signed Boolean collocation matrices BandGlocalize interface forces and deﬁne compatibility equations, respectively, and are deﬁned as: B= 0 l 0 ;G= l 0 0 (A.5) Boolean collocation matrices llocalize interface DoF pairs across NS and PS. In principle, at each simulation step, the computation of the state vector (A.4) is split into freeandlinksolutions. The coupledsolution is the sum of freeandlinkcontributions. According to (A.4)-(A.5), two parameters deﬁne the setting of the algorithm, testing time scale , deﬁned as =tC tP=tS tN subcyclingm, deﬁned as m=tN tP=tS tC Time steps tNandtPrefer to simulation time, which is the virtual time axis deﬁned by the time integration process. As an example, the seismic ground motion history refers to the simulation time axis. In particular, tNis the time step adopted for integrating the NS response while tPrefers to the PS. On the other hand, tSandtCrefer to wall-clock time measured in the laboratory. The former, tS, deﬁnes the maximum solving time allocated to compute the NS response while the latter, tC, is the actuator controller time step, which typically ranges between 1 and 2msec. In order to minimize discontinuity of actuator trajectories, the PS equation should be evaluated within a single controller time step tC. The testing time scale deﬁnes the ratio between wall-clock and simulation time rates. In particular, = 1 corresponds to real-time HS. If  >1, the simulation time ﬂow is slower than the wall-clock time. This case corresponds to the pseudodynamic HS regime. When the response of the PS is not rate-dependent, usually rangesbetween 50and200. Thisapproachreducesdisplacementtrackingerrorsaswellasthedestabilizingeﬀect of actuator delay, which typically ranges from 10to20msec. The subcycling parameter mallows for adjusting tS, whichistheﬁxedtimewindowallocatedforsolvingtheNSresponse. Accordingly, optimaltuningoftesting time scaleand subcycling mis a trade-oﬀ between computational accuracy and experimental capacity. For a detailed description of the partitioned time integration algorithm, the reader should refer to [32]. 14 DRAFTAcknowledgments This work was supported by the Swiss Secretariat of Education Research and Innovation (SERI) - Swiss Space Oﬃce (SSO) [THERMICS MdP2016 Project (Thermo-Mechanical Virtualization of Hybrid Flax/Carbon Fiber Composite for Spacecraft Structures), grant number REF-1131-61001].