While adaptive control schemes have been developed for SISO-RTHS (e.g., Strano and Terzo (2016), Huang et al. (2023)), they have yet to be applied in MIMO-RTHS. In this work, the proposed adaptive compensation framework is realized for coupled MIMO systems, such as those presented in Equations 1–4, accounting for system uncertainties ( Δ θ ) and capturing cross-coupling effects between actuators. A typical maRTHS scheme incorporating the proposed adaptive compensation framework is shown in Figure 1. The input signal to the system (i.e., external force), is denoted by F e x t ; the restoring force vector from the experimental substructure is denoted as f m ; the desired and measured signal vectors are x d and x m , respectively; u c is the command signal vector; and the control plant is represented in Laplace domain as H s θ , Δ θ [as shown in Section 2.1]. Hereafter, the control plant is simply denoted as H s .

The collection of blocks (I), (II), and (III) in Figure 1 generates the control law, which is realized by utilizing the inverse dynamics of the control plant derived from the transfer function in Equation 1. The control law incorporates adaptive parameters that are estimated by block (IV) at each time step. The parameter estimation process is described in Section 2.3. The remainder of this section explains the control generation process.

The transfer function matrices H s and G s relate the desired signal x d to the command u c and the measured signal x m as follows, X m = H s U c (5) U c = G s X d (6) where X m , X d , and U c denote the signals x m , x d , and u c in Laplace domain, respectively. From Equation 5 and Equation 6, it can be seen that the plant dynamics are cancelled when G s = H − 1 s , leading to the measured signal matching the desired signal X m = X d . However, a direct inverse of H s results in improper transfer functions for G s (Carrion et al., 2009; Phillips et al., 2014). To compute this inverse effectively, the proposed compensator makes use of the time derivatives of the desired signal. First, a common denominator polynomial of H s is found and denoted as ∆ s . The polynomial ∆ s incorporates poles and zeros from each actuator (2 actuators in this benchmark problem–see Section 3.1), capturing the coupling dynamics in the MIMO system. Then, by introducing X m = X d and G s = H − 1 s , Equation 5 can be rewritten as, X d = G − 1 s U c = 1 Δ s H ∼ s U c (7) where H ∼ s = Δ s H s , and thanks to Δ s , the i j -th entry of H ∼ s , H ∼ i j s , is a polynomial rather than a fraction and its inverse can be properly realized. It should be noted that the output in Equation 7 is now considered as X d , under the assumption that the desired and measured signals match ( X m = X d ). Further manipulation of Equation 7 yields the command signal vector U c to the plant, expressed as, U c = R s θ ^ k Δ s X d s = R s θ ^ k   Ψ d (8) where R s θ ^ k = H ∼ − 1 s is a matrix of transfer functions that can be properly realized. The term Δ s X d s is realized in the time domain by computing the function ψ d = g γ , θ ^ k , shown as the block (II) in Figure 1, where γ = x d , x ˙ d , x ¨ d , … denotes the collection of the time derivatives of x d [see block (I)]. The number of required time derivatives is determined by the order of the polynomial Δ s . θ ^ k denotes the estimated plant parameters at time step k , and Ψ d denotes the Laplace transform of the signal ψ d .

The presented approach addresses the issue of improper transfer functions in G s and facilitates the command generation. The following section describes the parameter estimation process for the updating of the plant poles and zeros ( θ ^ ) used in the control law of Equation 8.

The physical domain of RTHS is subject to disturbances such as imprecisions in experimental setup, unmodeled process and dynamics (e.g., friction), and process and measurement noises, which introduce uncertainties to the plant ( Δ θ ) during RTHS, as indicated in Section 2.1. To address the desynchronization issues caused by these uncertainties, this work employs online parameter estimation to update the plant parameters in real-time during the command generation. As observed from the two control plant formulations from Equations 1, 2 (for the TF model) and Equations 3, 4 (for the discrete-time SS model), the system formulations are nonlinear with respect to θ . Therefore, nonlinear estimators, EKF and UKF, are employed to estimate θ with the two introduced control plant models (see Section 2.1).

Both EKF and UKF are formulated for a general discrete-time system model expressed as, S t a t e   t r a n s i t i o n :   ϕ k = f ϕ k − 1 , μ k − 1 + w k − 1 (9) O b s e r v a t i o n   m o d e l :   y k = h ϕ k , μ k + v k (10) where f ϕ k − 1 , μ k − 1 is the process function; ϕ k is the state vector; μ k is the system input vector; w k is a white noise with zero mean and covariance matrix Q; y k is the system output vector; h ϕ k , μ k is the observation function; and v k is a white noise with zero mean and covariance matrix R. It should be noted from Equations 9, 10, that both process and measurement noises are considered as additive in this work.

Depending on the plant model employed for the application of EKF or UKF, the corresponding estimation system models are explained below.

a) Transfer function model

For the control plant in transfer function (TF) form, the state vector and process function at time step k in Equations 9, 10 are defined as the plant parameters at that time step, resulting in the following state transition function, θ k ⏟ ϕ k = θ k − 1 ⏟ f ϕ k − 1 , μ k − 1 + w k − 1 TF ⏟ w k − 1 (11)

The observation function is set as the command signal u c to the control plant. Therefore, Equation 10 becomes, u c , k ⏟ y k = u ∼ c , k θ k , γ x m ⏟   h ϕ k , μ k + v k TF ⏟ v k (12)

In this formulation, the nonlinearity arises from Equation 12, where the computation of the system output vector u c , k ( u c at the k -th time step) involves the polynomial expansion of the plant poles and zeros, as described in Section 2.2 (Equation 8). Note that u ∼ c , k denotes the observation function for u c , k , and w k TF and v k TF are the additive noise vectors for the state vector (i.e., parameter θ k ) and observation vector (i.e., command signal   u c , k , respectively. Additionally, unlike in Section 2.2 where x d is used to calculate u c , the parameter estimation step utilizes the measured signal x m and its respective time derivatives γ x m [ i . e . , ; γ x m = x m , x ˙ m , x ¨ m , … ].

b) Discrete-time state-space model

When the discrete-time state-space (SS) model of the control plant is used for parameter estimation, the original states of the system from Equations 3, 4 are augmented with the parameters that need to be estimated. Thus, the state transition and observation functions for the augmented discrete-time SS model are formulated as, z a k ⏟ ϕ k = A a k − 1 θ k − 1 z a k − 1 + B a k − 1 θ k − 1 u c , k − 1 ⏟ f ϕ k − 1 , μ k − 1 + w k − 1 SS ⏟ w k − 1 (13) x m , k ⏟ y k = C a k θ k z a k + D a k θ k u c , k ⏟ h ϕ k , μ k + v k SS ⏟ v k (14) where z a k is the augmented state vector at the k -th time step, i.e., z a k = z k , θ k T , of dimension p + n , with p and n being the dimensions of the original state-space vector z k , and the parameter vector θ k , respectively. The input and output system vectors correspond to the command and measured signals u c , k and x m , k . The additive noise vectors for z a k and x m , k are denoted as w k SS and v k SS , respectively. A a k , B a k , C a k , and D a k are the augmented state-space matrices at the k -th time step, computed using the forward Euler method as, A a k θ = I p + n × p + n + Δ t ⋅ A θ p × p 0 p × n 0 n × p 0 n × n (15) B a k θ = Δ t ⋅ B θ p × h 0 n × h (16) C a k θ = D θ h × p 0 h × n (17) D a k θ = D θ h × h (18) where Δ t is the sampling time of the RTHS, and h is the dimension of u c , k or x m , k , which is assumed to be the same. This formulation introduces the nonlinearity through the dependence of the state-space matrices on θ k , as indicated in Equations 13–18.

Regardless of the selected system model (TF or SS), both EKF and UKF can be applied for parameter estimation. The estimation methods are depicted in Figure 2 utilizing the general system formulation of Equations 9, 10.

EKF and UKF are extensions of the standard Kalman filter, designed to handle nonlinear systems (Song, 2011). Both methods consist of two recursive steps: the prediction step and the update step, as shown in Figure 2. In the prediction step, the process function f is used to predict the prior state estimate ( ϕ ^ k − ) and its covariance ( P k − ). In the update step, this prediction is updated by incorporating the observation function h to obtain the posterior state estimate ( ϕ ^ k ) and its covariance ( P k ), while minimizing the mean square error of the estimates with the Kalman gain K k . The details for the implementing EKF and UFK algorithms are omitted here, but can be found in (Song, 2011; Song and Dyke, 2013; Song and Dyke, 2014).

The primary difference between these estimation methods lies in handling the nonlinearities of the process and observation functions. The EKF approximates the nonlinear functions by linearizing them around the current state estimate using a first-order Taylor series expansion, which involves calculating their respective Jacobians ∇ f ϕ and ∇ h ϕ (see Figure 2). The calculation of the Jacobians is omitted here for succinctness. In contrast, UKF uses a deterministic sampling scheme known as the Unscented Transform (UT) to generate a set of “sigma points” that capture the mean (e.g., ϕ ^ k − ) and covariance of the state distribution (e.g., P k − ) without the Jacobian calculations (Julier, 2002). The sigma points, indicated in Figure 2 with the superscript “ S P ”, are carefully chosen to accurately reflect the state distribution. These sigma points are directly propagated through the nonlinear functions, resulting in transformed points ϕ ^ k * − and y ^ k * − for the state and output, respectively. The mean and covariance of the predicted state and observation are then computed as a weighted sum of their respective transformed sigma points (Wan and Van Der Merwe, 2000). Additionally, both EKF and UKF require an initialization of the state and covariance, indicated in Figure 2 as ϕ ^ 0 and P 0 , respectively, along with the selection of appropriate covariances Q and R , which involves a search process (Song et al., 2017; Song, 2018).

By replacing the general system formulation of Equations 9, 10 in Figure 2 with the selected control plant formulations from Equations 11, 12 (for the TF model) or Equations 13, 14 (for the discrete-time SS model), four (4) parameter estimation cases can be obtained: EKF-TF, UKF-TF, EKF-SS, and UKF-SS. These cases are illustrated in Figure 3, showing the architecture inside of block IV from Figure 1.

As shown in Figure 3, the same low-pass filter is applied to the signals u c and x m . The low-pass filter is indicated with dashed lines in block ( IV . b ) because this work also explores the option of not applying low-pass filtering for the SS model case (see Section 3.4). A difference between the TF and SS approaches is the definition of the system input and output vectors, μ k and y k . In the TF model (see Equations 11, 12), the input is γ x m [ i . e . , ; γ x m = x m , x ˙ m , x ¨ m , … ], which is related to the measured displacement of the plant x m and its time derivatives, and the output is the command signal u c . In contrast, the SS model (see Equations 13, 14) swaps the two with u c as the input and x m as the output, and therefore without the need for its time derivatives.

The proposed methodology is applied to the maRTHS benchmark developed by Condori Uribe et al. (2023), which provides a virtual RTHS (vRTHS) platform that allows the implementation of customized compensators. In this benchmark problem a three-story, three-bay moment resisting frame is subjected to seismic base excitation corresponding to scaled ground acceleration records of El Centro 1940, Kobe 1995, and Morgan Hill 1984. The central frame of the structure (one-story, one-bay, simply supported) acts as the experimental substructure (see Figure 4), and the remaining structure is the numerical substructure. Displacement and rotational degrees of freedom (DOFs) are used as target signals at the NS-PS interface (control node). These DOFs are refer to as “frame coordinates”, and are represented in Figure 4 as ψ e s , 28 and ψ e s , 4 , where the subscripts 28 and 4 denote the corresponding rotational and translational DOFs defined by Condori Uribe et al. (2023). The NS-PS connection is accomplished through a transfer system that consists of two hydraulic actuators and a steel coupler (see Figure 4). The coupler transfers the linear displacement of the two actuators to produce the desired frame coordinates. The desired and measured actuator linear displacement vectors are referred to as “actuator coordinates”. Thus, a coordinate transformation is performed during the simulation to transition between these two coordinate systems. The virtual RTHS is run with a fixed sampling frequency of 1 , 024 Hz (i.e., the sampling period Δ t = 0.976  ms .

In this work, the compensation is performed in the actuator coordinates domain. Therefore, the desired signal is specified as x d = η n 1 η n 2 T , where η n 1 and η n 2 denote the desired linear displacements for actuators 1 and 2, respectively. The command signal is set as u c = u 1 u 2 T , where u 1 and u 2 correspond to the linear displacement commands to actuators 1 and 2, respectively. Likewise, the measured signal is indicated as x m = η m 1 η m 2 T , where η m 1 and η m 2 denote the measured linear displacements of actuators 1 and 2, respectively, as shown in Figure 4.

The control plant model for the maRTHS benchmark problem is given as Equation 19, H s = H 11 H 12 H 21 H 22 (19) where the subscripts represent the input-output pairs for actuator 1 and actuator 2, respectively. Specifically, u 1 and u 2 are the inputs to the system, while η m 1 are η m 2 are the outputs. The adaptive compensator is designed considering the nominal plant model identified by Condori Uribe et al. (2023), which has the following transfer functions, H 11 s = K 11 s + B 1   s + B 2 s + A 1   s + A 2   s + C 1 + C 2 j s + C 1 − C 2 j (20) H 12 s =   K 12   s + P 1 s + D 1   s + D 2   s + C 1 + C 2 j s + C 1 − C 2 j (21) H 21 s = K 21 s + B 1   s + B 2 s + A 1   s + A 2   s + C 1 + C 2 j s + C 1 − C 2 j (22) H 22 s =   K 12   s + P 2 s + D 1   s + D 2   s + C 1 + C 2 j s + C 1 − C 2 j (23)

From Equations 20–23, C 1 , and C 2 denote the poles of the numerical substructure; B 1 , B 2 , P 1 , and P 2 represent the zeros of the transfer functions of the transfer system (i.e., actuators and steel coupler); A 1 , A 2 , D 1 , and D 2 denote the poles of the transfer functions of the transfer system; and K 11 , K 12 , K 21 , and K 22 are the corresponding transfer function gains.

The parameters in Equations 20–23 are updated in real-time during RTHS, while the transfer function gains are considered as constants. This approach results in the definition of θ , as shown in Equation 24. θ = A 1 A 2 B 1 B 2 C 1 C 2 P 1 P 2 D 1 D 2 T (24)

Once the control plant model H s and the nominal parameter θ * are established, the four parameter estimation cases are formulated as described in Section 2.3.

As indicated in Section 2.3, a search process is performed to select the initial values of ϕ ^ 0 and P 0 , and the covariances Q and R . The details of the search are omitted here but readers can refer to (Song, 2011; Song et al., 2017; Song, 2018) for more information. This search process may not guarantee to yield the optimal selection, and other search methods can be also considered. The selected EKF and UKF parameters (covariances) are presented in Table 1.

In addition to the covariances presented in Table 1, scaling parameters used during the sigma points computation in UKF were adjusted to enhance the estimation performance. The values for these scaling parameters are set as follows: α = 0.4 , κ = 0 , and β = 2 . Details on how these parameters influence the sigma points distribution can be found in Wan and Van Der Merwe (2000) and Song (2011).

The proposed adaptive compensator is validated through the vRTHS platform in MATLAB/Simulink R2023b (Mathworks, 2024). The performance is evaluated with a set of 10 performance criteria defined in Condori Uribe et al. (2023). Performance indices (PI), J n , m , are utilized for each criterion, where n indicates the criteria number and m depends on the type of evaluation. For example, PIs from J 1 , m to J 4 , m evaluate the tracking performance in actuators coordinates ( m denotes the actuator number). PIs J 5 , m and J 6 , m indicate the tracking performance in frame coordinates ( m denotes the DOFs of the control node). PIs from J 7 , m and J 10 , m indicate the global RTHS performance, which seek to minimize the error between the reference structure response and the hybrid system response ( m denotes the DOFs of the control node and upper stories levels). Lower PIs values indicate better compensation performance, as they correspond to smaller errors in both tracking and global performance.

A total of 1,000 vRTHS simulations were executed to demonstrate the control robustness, each incorporating system uncertainties ( Δ θ ) through random variations in the plant parameters (i.e., poles and zeros), which were modeled as normal random variables. Details in generating the above random variations are described in the benchmark problem definition in Condori Uribe et al. (2023). The initial values of the estimated parameters are set to the values of the nominal plant provided in Condori Uribe et al. (2023). All simulations were performed on a personal computer equipped with an Intel Core i7-14700 processor and 64 GB of RAM, running Microsoft Windows 11.

The mean and standard deviation (SD) for each performance indicator (PI) across the four parameter estimation cases using 0.4 scaled El Centro Earthquake as input, are presented in Table 2. Additionally, the results from the linear quadratic regulator (LQR) controller presented in the benchmark (BM) by Condori Uribe et al. (2023) are included in the table for comparison. Results utilizing 0.4 scaled Kobe and 0.4 scaled Morgan earthquakes are presented in the Supplementary Material. Furthermore, Table 3 presents the average simulation time for each earthquake input to evaluate computational efficiency.

The results from Table 2 demonstrate the effectiveness and robustness of the parameter estimation cases for the proposed compensator using the 0.4 scaled El Centro Earthquake as input. The robustness is confirmed by the performance across the 1,000 simulations, which resulted in small PIs with small standard deviations. For time delay, all cases were effective, with the highest mean time delay of 0.47  ms , obtained in the UKF-TF case. The tracking performance was effective in both actuator and frame coordinates ( J 1 , m to J 6 , m ) with the largest PIs obtained for the UKF-TF case. For the global RTHS evaluation, all performance estimation cases outperformed the BM-LQG for the degrees of freedom at the control node, although higher errors were obtained at DOFs of the upper levels of the reference structure. Overall, the UKF-SS case exhibited the best tracking performance with PIs values of J 2 , 1 = 2.68 ± 0.02 and J 2 , 2 = 1.42 ± 0.01 (actuator coordinates), and with J 5 , 4 = 1.82 ± 0.01 and J 5 , 28 = 0.85 ± 0.03 (frame coordinates). The largest obtained PIs from the 1,000 simulations for the same UKF-SS case were J 2 , 1 = 2.75 and J 2 , 2 = 1.47 (actuator coordinates), and with J 5 , 4 = 1.86 and J 5 , 28 = 0.96 (frame coordinates). Figure 5 shows the desired and measured displacements in actuator coordinates of actuator 1 for the UKF-SS case with the largest J 2 , 1 , demonstrating the close match of the signals over time. Similar tracking performance is observed for actuator 2 and for signals in frame coordinates. Therefore, the time history plots for these signals are not included here. Figure 6 illustrates the parameters tracking history (all normalized to the initial nominal values) across the four parameter estimation cases, with the 0.4 scaled El Centro Earthquake as input. This figure demonstrates that meaningful parameter estimation begins around 5 s, when the earthquake input starts (see Figure 5). Parameters continue to update until reaching a stable value to minimize the estimation error. Although starting with the same initial values, the updating histories among the four cases are not consistent. For example, among all the parameters, P 1 shows the largest variation in the EKF cases, with changes of ∼ 12 % in EKF-SS and ∼ 30 % in EKF-TF. In contrast, the UKF cases exhibit the largest variations with D 2 ( ∼ 10 % ) in UKF-SS and with P 2 ( ∼ 2 % ) in UKF-TF. Figure 7 illustrates the tracking performance of the proposed compensation for all three input earthquakes in both actuator and frame coordinates. The height of each bar represents the mean value of the PI for the respective parameter estimation case, while the error bars indicate one standard deviation (+SD) of the corresponding PI value. Overall, the small PIs values in the bar graph demonstrates the effectiveness and robustness of the compensator across the three earthquake input scenarios. For each earthquake input, the UKF-SS yields the smallest PIs among all compensation cases, while UKF-TF shows the largest PIs and standard deviations, likely due to its less active updating in parameter tracking (see Figure 6D).

Although the UKF-SS shows superior tracking performance compared to the other cases studied, it is also more computationally expensive, as shown in Table 3. This increased computational cost is likely due to the requirement for the UKF algorithm to generate and propagate multiple sigma points, which involves additional matrix operations. Note that, the simulation time indicated in Table 3 is measured during the numerical simulations in MATLAB under Windows environment, not in a real-time computing environment. But by considering the same computational environment, Table 3 can still offer a relative comparison of the computational efficiency of each algorithm.

Results presented in Section 3.3 are obtained with the use of a 6th-order Butterworth low-pass filter with a cutoff frequency of 20 Hz, applied to u c and x m (see Figure 3). The low-pass filter helps attenuate high-frequency noise present in the measured signals. This noise is exacerbated during time differentiation when using the TF model, making the use of the low-pass filter crucial to ensure the stability of EKF-TF and UKF-TF cases. In contrast, the discrete-time SS model does not involve time derivatives of the signals and only deals with measurement and process noises. Therefore, the low-pass filter can potentially be removed when using the discrete-time SS model. Table 4 presents the tracking performance PIs ( J 1 , m to J 6 , m ) for EKF-SS and UKF-SS cases without low-pass filtering, using 0.4 scaled El Centro earthquake as input in 1,000 simulations. The EKF-SS results, compared to those in Table 2, showed significant decrease in performance, leading to large SDs in Table 4. However, for the UKF-SS case, the obtained PIs are consistent with the values in Table 2, with only a small increase in the PIs. Additionally, Figure 8 shows the bar charts for J 2 , 1 and J 2 , 2 (values from Tables 2, 4) for the mentioned EKF-SS and UKF-SS cases with and without low-pass filter. The bar height shows the mean PI value, while the error bars indicate one standard deviation (+SD) of the corresponding PI.

These results reveal that UKF can handle noises and nonlinearities in the system model more effectively than EKF when low-pass filter is not present. This capability is advantageous in RTHS because the use of low-pass filter can potentially introduce signal distortion leading to unwanted dynamics with potentially nonlinear systems, and therefore compromise the compensation performance and the accuracy of the simulated response.