A Global Sensitivity Analysis Framework for Hybrid Simulation with Stochastic Substructures

Hybrid simulation is an experimental method used to investigate the dynamic response of a reference prototype structure by decomposing it to physically-tested and numerically-simulated substructures. The latter substructures interact with each other in a real-time feedback loop and their coupling forms the hybrid model. In this study, we extend our previous work on metamodel-based sensitivity analysis of deterministic hybrid models to the practically more relevant case of stochastic hybrid models. The aim is to cover a more realistic situation where the physical substructure response is not deterministic, as nominally identical specimens are, in practice, never actually identical. A generalized lambda surrogate model recently developed by some of the authors is proposed to surrogate the hybrid model response, and Sobol’ sensitivity indices are computed for substructure quantity of interest response quantiles. Normally, several repetitions of every single sample of the inputs parameters would be required to replicate the response of a stochastic hybrid model. In this regard, a great advantage of the proposed framework is that the generalized lambda surrogate model does not require repeated evaluations of the same sample. The effectiveness of the proposed hybrid simulation global sensitivity analysis framework is demonstrated using an experiment.


Background and motivation
Hybrid simulation (HS) is used to investigate the experimental response of a structural component or subassembly subjected to a realistic loading scenario, which includes the interaction with a credible yet virtual structural system. Coupled physical and numerical substructures (PS and NS) form the so-called hybrid model. Specifically, the PS is tested using servo-controlled actuators provided with force transducers, while the NS is instantiated in a structural analysis software. A time-stepping analysis algorithm computes the hybrid model response on-the-fly. This ensures displacement compatibility and force balance between NS and PS throughout the entire experiment. When the PS response is rate-independent, it is good practice to perform HS with an extended time scale from 50−200 times slower than real-time in order to improve the accuracy of position/force control. HS is used to investigate the inner workings of a structural component beyond the linear regime without testing an entire structural assembly. As a result, the cost of experimentation is substantially reduced. A report by Schellenberg and co-workers provides a comprehensive review of HS methods and algorithms (Schellenberg et al., 2009).
In earthquake engineering, HS is the only viable solution for testing large structures such as bridges (e.g., (Bursi et al., 2017;Abbiati et al., 2019)). Similarly, HS has been recently proposed to test mooring lines of offshore structures for which hydrodynamic tests with sizeable scale are prohibitive (e.g., (Sauder et al., 2018;Vilsen et al., 2019)). HS is gaining popularity for component-level testing in fire engineering. The reason is that internal force redistribution occurring at the system level heavily influences the failure modes at the component level (e.g., (Whyte et al., 2016;Abbiati et al., 2020)). More recently, HS has been combined with centrifuge testing for investigating soil-structure interaction problems (Idinyang et al., 2019).
In all the cases reviewed above, NS and related excitation are conceived as deterministic. However, in the majority of cases encountered in structural engineering, loading is stochastic, while damping and boundary conditions are highly uncertain. An exhaustive exploration of all possible load cases is clearly not an option given the experimental cost associated with a single hybrid model evaluation. Accordingly, Abbiati and coworkers (Abbiati et al., 2021) proposed surrogate modeling to compute the variance-based global sensitivity analysis (GSA) of the response quantity of interest (QoI) of a given hybrid model with respect to a set of input parameters that characterize both substructures and loading excitations. In detail, polynomial chaos expansion (PCE) was used to construct a surrogate model (a.k.a. response surface) of the hybrid model response. Sobol' sensitivity indices of the QoIs were obtained as a by-product of polynomial coefficients as explained in (Sudret, 2008). The goal of such an approach was to reveal what influences what within the HS: This entails uncovering the inner workings of the PS, the unknown part of the hybrid model in both epistemic and aleatory sense. In that study, the PS was treated as deterministic, that is, aleatory uncertainty was neglected by assuming that nominally identical specimens have identical responses (plus some negligible measurement noise).
This assumption, however, is still far from a realistic scenario in structural testing. Nominally identical specimens are, in practice, never actually identical. Also, some sources of loading exerted through the PS are inherently stochastic (e.g., fire or hydrodynamic loading). As a result, uncertainty, both aleatory or epistemic, always affects the PS structure behavior.

Scope
This paper extends the GSA framework for HS proposed in (Abbiati et al., 2021) to the case of PSs with non-deterministic behavior. Similarly to the original framework, the idea is to surrogate the hybrid model response as a function of the input parameters that can be controlled by the experimenter and originate from substructures and loading (physical and numerical). However, latent variables that do not appear in the input parameter vector make the hybrid model response stochastic. To account for the latter, the generalized lambda model recently developed by Zhu and Sudret (Zhu and Sudret, 2020c,a) is used here to directly surrogate the probability density function (PDF) of the response QoI. This is achieved by means of the family of generalized lambda distributions, which are suitable to approximate a wide class of distributions commonly found in engineering contexts (Karian and Dudewicz, 2000). The parameters of the generalized lambda distributions are cast as functions of the input parameters of the hybrid model and approximated via PCE (Xiu and Karniadakis, 2002;Blatman and Sudret, 2011). Normally, several repetitions of every single sample of the inputs parameters would be required to replicate the response of a stochastic hybrid model. However, acquiring such repetitions would be impossible in a HS. Instead, by using the generalized lambda model presented in (Zhu and Sudret, 2020a) we can solve this problem, as the generalized lambda model can be computed in a non-intrusive manner and does not require repeated HSs for a single sample of input parameters. For these reasons, the generalized lambda model is well-suited to surrogate the response of a stochastic hybrid model. Variance-based GSA is uniquely defined for deterministic simulators (Saltelli, 2008). In the context of stochastic simulators, three alternative variants of Sobol' sensitivity indices are discussed in (Zhu and Sudret, 2020b), namely classical, quantile-based and trajectory-based Sobol' indices. In this work, quantile-based Sobol' indices are used for the GSA of the HS response QoI.
The effectiveness of the proposed framework is demonstrated for a 3-DoF hybrid model subjected to mechanical and thermal loading. Specifically, thermal loading is experimentally exerted on the PS so that temperature fluctuations (out of control of the experimenter) entail a stochastic response of the hybrid model . This paper is organized as follows. Section 2 introduces the generalized lambda model and describes the GSA framework for stochastic hybrid models. Section 3 presents the 3-DoF hybrid model used to test the framework. Section 4 discusses the results of the GSA of the stochastic hybrid model response. Finally, Section 5 presents the overall conclusions of this study.

Global sensitivity analysis framework
Let a random variable vector x = {x 1 , . . . , x M } ∈ D X ⊂ R M , where D X denotes the range of definition of x, represent the input parameters to a HS. Due to the random nature of the hybrid model response, for a given set of parameters x, the corresponding QoI Y (x) is a random variable rather than a deterministic value. This is because some latent variables z cannot be identified or measured in the process, which makes it impossible to include all the relevant variables in x. Therefore, a stochastic HS can be expressed as a mapping: (1) The latent variables z are grouped into a random vector Z. Note that with a fixed x and Z varying according to some probability distribution, the HS output Y (x) remains random.
The input parameters of vector x are treated as random, modeled by known probability distributions and grouped into a random vector X = {X 1 , . . . , X M }. X is characterized by its joint distribution with the PDF denoted by f X . Furthermore, we assume that X i 's are mutually independent, and thus the joint PDF is the product of marginal PDFs, i.e. f X (x) = M i=1 f Xi (x i ) with f Xi being the marginal PDF of the i-th variable. For a given set of input parameters x, the QoI Y (x) is a random variable characterized by an unknown conditional probability distribution. Therefore, representing the stochastic behavior of a HS consists in estimating the response distribution for any parameters x ∈ D X . However, one simulation for x does not provide the whole probability distribution but rather a single realization of Y (x). Hence, it is usually necessary to repeatedly conduct experiments for the same x (called replications) to have enough insight into the resulting hybrid model response probability distribution. This quickly becomes intractable when the number of x's to be investigated increases. To alleviate the burden, surrogate models can be constructed to emulate the stochastic behavior of a HS. Once a surrogate model is constructed, we can perform further analysis of the hybrid model response at a low cost, namely the GSA. The simplest surrogate model of a stochastic HS involves additive Gaussian noise: (2) To build such a surrogate, one needs to estimate the mean function h and the noise variance σ 2 . In this case, PCE (Xiu and Karniadakis, 2002;Berveiller et al., 2006) and Gaussian processes (Rasmussen and Williams, 2006) with a regression setup can be directly applied. However, Eq.
(2) can be rather restrictive. To cover a wider group of problems, we choose to use the recently developed statistical model called the generalized lambda model (Zhu and Sudret, 2020c,a).

Generalized lambda models
A generalized lambda model (GLAM) assumes that the probability distribution of hybrid model response QoI Y (x) can be approximated by a generalized lambda distribution (GLD). The latter is a highly flexible four-parameter distribution family, which is able to approximate many common distributions such as normal, lognormal, uniform and extreme value distributions (Freimer et al., 1988). A GLD is defined by its quantile function where u ∈ [0, 1] and λ = {λ l : l = 1, . . . , 4} are the four distribution parameters. More precisely, λ 1 is the location parameter, λ 2 is the scaling parameter, and λ 3 and λ 4 are the shape parameters. To have valid quantile functions, λ 2 should be positive. From Eq.
(3), we can derive the associated PDF: where 1 [0,1] is the indicator function. From the above equation, it is clear that evaluating the PDF for a particular y requires solving numerically the equation u = Q −1 (y; λ).
Under this setup, varying x will lead to Y (x) following a GLD with different distribution parameters λ. In other words, λ l 's are functions of x, which allows us to express the QoI as Recall that the input parameters x are modelled as independent random variables X with joint PDF Under appropriate assumptions (Ernst et al., 2012), each component of λ(x) admits a polynomial chaos (PC) representation where {ψ α : α ∈ N M } is a basis of multivariate polynomials that are mutually orthogonal with respect to the probability measure of X, α is a multi-index that identifies the polynomial degree in each of the input variables, A l ⊂ N M is a truncated set defining a finite set of basis functions for λ l (x) and c = {c 1,α , . . . , c l,α } denotes the associated coefficients (see (Zhu and Sudret, 2020c) for details). Note that the polynomial chaos expansion for λ 2 (x) is built on the logarithmic transform so as to ensure that λ PC 2 (x) is always positive. Combining Eq. (5) with Eq. (6), we define the generalized lambda surrogate model To build a GLAM, we need to determine the the associated coefficients c. To avoid the need for replications, which may result in a large number of experiments, we use the method developed by Zhu and Sudret (Zhu and Sudret, 2020a). We first generate a set of realizations of the input random vector X = x (1) , . . . , x (N ) , called the experimental design (ED). For each point of the ED, we conduct a HS and collect the corresponding QoI into Y = y (1) , . . . , y (N ) . Note that each HS may correspond to a different realization of the latent variable Z, which does not need to be explicitly known in the analysis. In a second step, we estimate c by maximizing the conditional likelihood, i.e. minimizing the negative log-likelihood, that iŝ where f GLD is the PDF of the GLD defined in Eq. (4). To solve the optimization problem, it is necessary to determine the support of c, which is equivalent to finding the truncation set A l for each λ PC l . To this end, we plug the hybrid-LAR algorithm (Blatman and Sudret, 2011) into the modified feasible generalized least-squares framework (Zhu and Sudret, 2020a). The latter fits the mean and the variance function in an alternative way. The basis functions selected for these two functions are then used to represent λ PC 1 (x) and λ PC 2 (x), respectively. As λ 3 and λ 4 mainly control the PDF shape of a GLD, which is expected not to change much when x is changed, we can pick polynomials with low degree, namely 0 or 1, for λ PC 3 (x) and λ PC 4 (x). After specifying the basis functions, we solve Eq. (4) to build the associated GLAM.

Sobol' sensitivity indices
Variance-based sensitivity analysis has been extensively studied and successfully developed in the context of deterministic models (Abbiati et al., 2021;Saltelli, 2008). For a random vector X with independent components, any deterministic mapping Y = M s (X) with Var [Y ] < +∞ can be decomposed as (Sobol', 1993) where M 0 is a constant, and M u (x u ) is a function of the subset of input variables x u specified by u = ∅, u ⊂ {1, . . . , p}. This decomposition is unique, and the elementary functions are defined by conditional expectations where |u| gives the cardinality of u. Following this definition, Moreover, the various terms M u are orthogonal with each other w.r.t. the inner product induced by the input PDF. Thus we can decompose the variance of Y as: Additionally, the definition in Eq. (10) allows for calculating the variance of M u (X u ) by The Sobol' index S u is defined as the ratio of V u to the total variance Var [Y ] For |u| = 1, we obtain the first-order Sobol' indices {S i : i = 1, . . . , p}, which represent the main effect of each input variable. Higher-order indices quantify the interactive effect within a given group of input variables. The total Sobol' index S Ti account for all the effect related to X i Due to the random nature of stochastic simulators, decomposition similar to Eq. (9) is generally impossible. Therefore, it is necessary to represent a stochastic model by a deterministic function to obtain the associated Sobol' indices. Based on the choice of the deterministic representation, we can have different extensions of the Sobol' indices (Zhu and Sudret, 2020b).
The most straightforward way is to include the latent variables within the input variables (Iooss and Ribatet, 2009). This leads to the underlying (yet unknown in practice) deterministic model M s defined in Eq. (1). Decomposing this function, we have As the definition of M u is the same as Eq. (10), the equation (12) for V u still holds. This implies that S u can be determined by the statistical dependence between Y and X u . The definition of the total index S Ti requires including the interactive effect between X i and the latent variables Z. However, interactions with Z cannot be determined by the response distribution but rather depend on the precise data generation process (i.e. how Z is present in the function M s ). Since the latent variables Z are generally impossible to characterize and control in a real experiment, the total Sobol' indices cannot be assessed.
Alternatively to Eq. (15), certain summary quantities of the response random variable Y (x) can be employed as a deterministic representation of Y (x). This is particularly helpful when the selected summary quantity itself is of interest. Typical quantities are mean m( (Iooss and Ribatet, 2009), and α-quantiles q α (x) (Browne, 2017). As these functions are well-defined as deterministic functions of x (since the effect of the latent variables Z has been marginalized), the associated Sobol' indices follow directly from Eq. (9).
A generalized lambda surrogate model emulates the response distribution of a stochastic model, which fully captures the statistical dependence between the input variables X and the QoI Y . Therefore, such a surrogate allows for evaluating both types of Sobol' indices mentioned above. More precisely, we can apply either Monte Carlo simulations or polynomial chaos expansions to the easy-to-evaluate emulator (see (Zhu and Sudret, 2020b) for more details).

Stochastic hybrid model
The proposed GSA framework is illustrated using a stochastic 3-DoFs hybrid model subjected to both thermal and mechanical loading, illustrated in Figure 1. As can be appreciated from Figure 1a, the hybrid model consists of a simply-supported beam provided with rotational elastic restraints. In this regard, u 1 and u 2 indicate the two rotational DoFs, while u 3 is the axial DoF. Figure 1b describes the substructuring of the hybrid model into PS and NS. The NS comprises the two elastic rotational restraints, the axial spring and the linear dashpot whereas the PS coincides with the beam element. Specifically, the axial spring is characterized by a constant stiffness of K 3 = 8100 × 10 3 N/m. Two lumped rotational masses are defined by J 1 = J 2 = 10 kgm while the lumped translational mass is defined by M 3 = 5000 kg. The linear dash-pot is characterized by C 3 = 1129 × 10 3 Ns/m. The PS consists of an aluminium plate of 0.2 × 0.002 m rectangular cross-section and length L = 0.47 m. Accordingly, the cross-section of the plate is characterized by an area A = 4 × 10 −4 m 2 and a moment of inertia I = 66.67 × 10 −12 m 4 . The Young's modulus, density and thermal expansion coefficient of aluminum are E = 69.5 GPa, ρ = 2700 kg/m 3 , and α = 23 × 10 −6 • C -1 , respectively. Since HS are conducted in pseudodynamic mode using a testing time scale equal to 50, the PS mass (rotational and translational) does not contribute to the hybrid model inertia. For the sake of clarity, all equations and plots in the following refer to simulation time, which is virtual and 50 times slower than the wall-clock time.
Mechanical loading is supplied as a bending moment history F (t) applied to the right rotational DoF u 2 while thermal loading, applied by a heating lamp, is defined by a ramp & hold temperature history θ(t). The expressions of both read: where F max and T are peak value and period of the half-side bending moment pulse applied to DoF u 2 with a time shift t 0 = 0.1 s;θ and θ max are temperature rate and plateau characterizing the temperature history imposed to the PS. For the sake of this example, Figure 2a depicts the bending moment history computed for F max = 45 Nm, T = 1 s. Similarly, Figure 2b depicts the temperature history computed forθ = 21.5 • C/s and θ max = 120 • C.
Consistent with the motivations underlying the development of the GSA framework, the stiffness of the two elastic rotational springs, which play the role of boundary conditions to the PS, as well as the loading parameters, are selected as input parameters for the surrogate modelling phase. In line with the procedure described in Section 2, the input parameters are described by independent uniform distributions, whose bounds are summarized in Table 1.
It is important to remark that, in order to reduce the experimental effort required to validate the proposed framework, the hybrid model and loading excitation were designed such that the PS always remained in the linear response regime. As a result, HS were conducted using a single aluminum plate.
The response QoI selected for the GSA corresponds to the maximum absolute out-of-plane deflection of the tested aluminum plate, which is denoted as u L,max .

Hybrid simulation setup
The 3-DoFs HS test rig used to conduct HS is a stiff loading frame equipped with four electro-mechanical actuators and an infrared (IR) lamp module interfaced to an INDEL real-time computer (Abbiati et al., 2018). The 3-DoFs HS test rig is designed to test plate specimens with an approximate footprint of 200 × 500 mm and thickness varying between 1 and 3 mm. Figure 3 illustrates the architecture of the HS setup, including a close-up view of the plate specimen accommodation. Two axonometric views of the 3-DoF test rig, consisting of the main hardware components are shown in Figure 4.
The moving parts of the test rig are colored in yellow, the plate specimen in brown and the fixed parts in gray. The latter are fixed to a reaction frame. In order to impose the u 1 and u 2 rotations, two rack-pinion systems (10) are installed along the vertical actuator y 1 and y 2 (1). The rack-pinion systems aim at transforming the commanded displacements from the actuators to rotational DoFs, applied to the short edges of the plate specimen (6) through aluminium clamps (3). The two horizontal actuators x 1 and x 2 control the position of the moving frame mounted on profiled rail guides (4) and corresponding to the axial DoF u 3 of the plate specimen (6). A linear variable differential transformer (LVDT) measures the out-of-plane deflection at the mid-span of the plate specimen (labeled u L in Figure 3). A Type K thermocouple installed at the center of the plate provides the feedback signal for the control of the IR lamp, which imposes the temperature history The GINLink bus connects actuator servo-drivers INDEL SAC4, IR lamp and all DAQ modules to the real-time computer INDEL SAM4, which executes the HS software. The latter is developed in MAT-LAB/SIMULINK, compiled, and downloaded to the INDEL SAM4 from the Host-PC. At each simulation time step, the HS software imposes displacements u 1 , u 2 and u 3 to the plate specimen, the PS, reads the corresponding restoring forces r 1 , r 2 and r 3 measured using force transducers, and solves the coupled equation of motion of the hybrid model. Also, the HS software generates the temperature command for the IR lamp. A comprehensive description of the time integration scheme used for HS is reported in (Abbiati et al., 2019).
Forces were manually set to zero before starting HS. An electric fan cools down the PS at the end of each test. Room temperature was quite stable and equal to 30 • C for the entire testing campaign.  (1) vertical actuators, (2) horizontal actuators, (3) installation clamps, (4) moving frame, (5) profiled rail guides, (6) plate specimen, (7) hinges, (8) vertical actuator load cells, (9) horizontal actuator load cells and (10) rack-pinion systems. The moving parts are colored in yellow, the plate specimen is brown, while those parts fixed to the reaction frame are gray. The latter is omitted in this figure.

Results and discussion
The response of the hybrid model described in Section 3.1 was evaluated using the HS setup described in Section 3.2 on 200 samples of the input parameter vector generated using latin hypercube sampling (Mckay et al., 2000). The resulting ED X, Y was used for computing surrogate models. In this regard, Figure 5 reports out-of-plane and axial displacement histories obtained via HS for a single sample of input parameters, where u L,max , u 3,max , u L,0 and u 3,0 scalar quantities are highlighted. Specifically, u L,max corresponds to the QoI (see Eq. (18)) and u 3,max to the absolute maximum axial displacement while u L,0 and u 3,0 indicate the initial position of the out-of-plane displacement and axial axes respectively, relative to the value measured during the first HS. Figure 6 describes the evolution of u L,0 and u 3,0 scalar quantities over the entire experimental campaign.

Drift observed in measurement data
From Figure 6 it is clear that both u L,0 and u 3,0 quantities have a constant drift, which results in a total accumulated out-of-plane and axial displacements of 3.0 mm and 0.3 mm, respectively. Such a drift can be reasonably ascribed to cumulative slippage of plate fixtures produced by heating/cooling cycles. This drift occurs regardless of the type of analysis that the HS setup was used for, and hence it should be removed from the acquired raw data before any further post-processing. Accordingly, prior to the calculation of surrogate models, the effect of drift on u L,max was eliminated via linear detrending with respect to u L,0 . The detrended QoI is referred to asû L,max and compared to original values in Figure 7. Notably, u L,0 is independent from the parameters of the hybrid model since the initial position of the PS was set by zeroing actuator forces.
Consistent with the notation introduced in Section 2, surrogate modeling was performed considering the following input parameter vector and response QoI,

Global sensitivity analysis framework results
A GLAM of the hybrid model response QoI was computed as explained in Section 2. In this study, we set the candidate degrees up to 5 for λ PC 1 and 3 for λ PC 2 . In order to validate the GLAM, 10 repeated HS were performed for two validation ED points, namely 58 and 157, associated with different regions of the input parameter space and characterized by appreciably different QoI values. For each validation ED point, Figure 8 compares the GLAM prediction to the empirical distribution of the 10 related repetitions. It is observed that the GLAM correctly captures the empirical distribution of the QoI for both points. It is interesting to note that the computed GLAM model converged to zero-order polynomials for λ PC 2 , and that the two PDFs of Figure 8 are similar in shape, thus suggesting an homoscedastic stochasticity of the hybrid model, which was further verified using a PCE surrogate model. Specifically, a residual analysis was performed on the difference between the measured QoI and its PCE.
As highlighted by Torre and co-workers (Torre et al., 2019), in presence of noisy data, PCE is a powerful denoiser, which naturally provides a surrogate of the average model response E Z [Y |x]. In this regard, the Tukey-Anscombe plot (Anscombe and Tukey, 1963) of Figure 9a compares the PCE output to the corresponding residual for each sample of the ED. The Q-Q plot of Figure 9b compares the empirical quantiles of the residuals normalized to unit standard deviation to the theoretical values of a standard normal distribution N (0, 1). The zero-average uniform scattering of residuals highlighted by the Tukey-Anscombe plot and the fairly good agreement between empirical and theoretical quantiles highlighted by the Q-Q plot confirms that the hybrid model response was affected by a Gaussian homoscedastic additive noise. As reported in Section 2.2, only first-and higher-order Sobol' indices but not total Sobol' indices can be obtained from the GLAM of the QoI. Instead, first and total Sobol' indices can be computed for the QoI quantiles. Accordingly, Figure 10 provides first and total Sobol' indices of the 5, 50 and 95 % quantiles of u L,max . The results of the GSA indicate that the temperature plateau value θ max is the most sensitive input parameter for the selected QoI. Additionally, the equal Sobol' indices values for each quantile unveils the homoscedastic response of the stochastic surrogate.
The development and implementation of the surrogate modelling, as well as the GSA, was performed using the UQLab software framework developed by the Chair of Risk, Safety and Uncertainty Quantification in ETH Zurich (Marelli and Sudret, 2014).

Conclusions
This paper described a framework for global sensitivity analysis of stochastic hybrid models. A generalized lambda surrogate modeling technique recently developed by some of the authors is used to compute Sobol' sensitivity indices for the quantiles of a response quantity of interest. The idea of using surrogate modeling to enable global sensitivity studies with few expensive-to-evaluate hybrid simulations was already presented in a previous work of the authors. The novelty of this work lies in the extension of global sensitivity analysis to the case of stochastic hybrid models, to cover the more realistic situation where the hybrid model response for two nominally identical physical substructures is not repeatable.
The effectiveness of the proposed framework is demonstrated in an experimental application consisting of a hybrid model with five parameters and subjected to mechanical and thermal loading. A great advantage of the proposed framework is that the generalized lambda surrogate model does not require repeated evaluations of the same sample. The results of the demonstration study highlight that the stochasticity of the particular hybrid model under consideration is homoscedastic with respect to the hybrid model parameters. Accordingly, both first-order and total Sobol' sensitivity indices of 5, 50, and 95 % quantiles are almost identical. Moreover, the temperature plateau value of the thermal loading is the most sensitive parameter for the selected response quantity of interest. The outcome of the experiment demonstrates the effectiveness of the proposed GSA framework in revealing the inner workings of the hybrid model.
Future research will address the issue of adaptive sampling of the parameter space of the hybrid model to minimize the experimental cost necessary to compute an accurate surrogate model for global sensitivity analysis.

Data Availability Statement
All data and code that support the findings of this study are available from the corresponding author upon reasonable request.