The maRTHS benchmark problem (Condori et al., 2023) is considered in this study for the development and validation of a proposed robust compensation algorithm that is capable of rejecting uncertainties mainly due to specimen-multiactuator interactions. In this benchmark problem, a three-story, three-bay, planar steel moment frame is chosen as the reference structure, as shown in Figure 1. Steel shapes and tributary gravity loads are provided to compute masses for dynamic transient analyses.

A multi-degree-of-freedom (MDOF) linear-time-invariant (LTI) numerical model is provided for the reference structure. The structural system has linear-elastic behavior and three degrees of freedom (DOFs) per node: two translational DOFs along the global x - and y -axes and one rotational DOF θ around the z -axis, perpendicular to the x y -plane. Therefore, the total number of DOFs in this model is 38, arranged in a displacement vector ψ according to the numeration in Figure 2, i.e., ψ = ψ 1 … ψ 38 T . Numeration begins with the horizontal displacements, followed by the vertical displacements, and finally, the rotations of each node. Moreover, only seismic loading in the horizontal direction is considered for this maRTHS study (i.e., no gravity loads are applied over the structural system).

The equation of motion (EOM) is given by Eq. 1: M ψ ¨ + C ψ ˙ + K ψ = − M Γ x ¨ g , (1) where M , C , and K are the mass, damping, and stiffness matrices of the reference structure, respectively. Γ is an influence vector that describes the inertial effects of the seismic excitation on the masses of the system, where a value of one is defined in the horizontal directions (i.e., the direction of the ground motion) and zero in the rest. x ¨ g is a scalar function associated with the ground acceleration in the x -direction. Finally, ψ ¨ , ψ ˙ , and ψ are the relative acceleration, velocity, and displacement vectors with respect to the ground for each DOF, respectively. Initial rest conditions are considered for this problem (i.e., ψ i 0 = ψ ˙ i 0 = 0 , where i = 1 , 38 ).

As a consequence, the first five natural frequencies are f 1 = 2.29 Hz, f 2 = 12.74 Hz, f 3 = 26.28 Hz, f 4 = 26.53 Hz, and f 5 = 29.91 Hz. Meanwhile, the damping matrix is calculated using the Rayleigh damping method with a damping ratio of 5% specified for the first and third modes, assumed to be proportional to the sum of the mass and stiffness matrices of the reference system.

In the context of RTHS, instead of solving the EOM for the entire reference domain, a process known as substructuring can be employed to subdivide the domain into smaller subdomains, such that the order of large and complex structural systems is reduced for efficient computations (Dermitzakis and Mahin, 1985). Each subdomain can be solved independently, provided that coupling between components is enforced through compatibility and equilibrium conditions at their interfaces.

In the benchmark case, the domain is decomposed into two non-overlapping subdomains, i.e., Ω = Ω N ∪ Ω E , where Ω N is called the numerical substructure and Ω E is the experimental substructure, as shown in Figure 3. The first-story central frame is selected as the experimental substructure, which is assumed to be the less-understood part of the entire structure.

If we assume the frame exhibits linear elastic behavior, the matrices for partitioned mass, damping, and stiffness can be expressed as the combination of numerical components (denoted with the superscript “ N ”) and experimental components (denoted with the superscript “ E ”): as shown in Eqs 2, 3: M N + M E ψ ¨ + C N + C E ψ ˙ + K N + K E ψ = − M Γ x ¨ g , (2) M = M N + M E , C = C N + C E , K = K N + K E . (3)

The terms associated with the experimental substructure can be arranged into the right-hand side of the EOM as described in Eqs 4, 5: M N ψ ¨ + C N ψ ˙ + K N ψ = − M Γ x ¨ g − f E , (4) f E = M E ψ ¨ + C E ψ ˙ + K E ψ . (5)

Here, f E is the feedback force vector produced by the experimental substructure during the dynamic testing, directly measured or estimated using measured data from the experimental domain.

Due to boundary conditions, both compatibility and equilibrium must be satisfied, as mentioned previously in this section. The first condition assures that the displacements at the interface of the subdomains are the same, allowing maRTHS to impose the numerical displacements (and their time derivatives) into the experimental substructure. In this case, they are given by the following DOFs: ψ N = ψ 4 N   ψ 7 N   ψ 16 N   ψ 19 N   ψ 28 N   ψ 31 N   T . The second condition establishes the relationship between the numerical and experimental forces as their sum at the interface must equal zero (i.e., force equilibrium). Then, in maRTHS, the generated physical restoring forces f E = F 4 E   F 7 E   F 16 E   F 19 E   M 28 E   M 31 E   T are measured and fed back to the numerical substructure. Both conditions are illustrated in Figure 3.

In maRTHS, the previously mentioned process is repeated until the simulation reaches the final simulation time t f , with real-time constraints. In that regard, Figure 3 can be considered the ideal maRTHS since there are no additional components such as hydraulic actuators, inertial effects (e.g., the mass of loading connectors), or measured signals contaminated with high-frequency noise. In an actual RTHS implementation, the experimental substructure is physically built and connected to the numerical substructure through a transfer system (i.e., loading assembly) to ensure dynamic synchronized motion at the substructure interface nodes during the real-time experiment. The insertion of the transfer system in the experimental domain changes the dynamic properties of the experimental substructure, and a control approach must be included to compensate for these added dynamic and other negative effects. The virtual simulation of this realistic computational implementation is called virtual RTHS (vRTHS), which is considered the closest realization of an actual RTHS for training and development purposes.

In this benchmark problem, the EOM is integrated using the fourth-order Runge–Kutta integration scheme incorporated in Simulink (ode4), and the sampling frequency F s is fixed at 1024 Hz.

In this section, we explain the transfer system, which consists of a set of servo-hydraulic actuators. First, the vertical DOFs along the global coordinate y are ignored due to their insignificant values compared to the horizontal displacements. The main reason for this assumption is the negligible axial deformations in columns since gravity loads are not incorporated in the maRTHS simulation and columns show high axial stiffness. On the other hand, due to experimental setup limitations, only the boundary conditions for node 4 of the experimental substructure are prescribed and implemented in the maRTHS benchmark problem. Therefore, only two actuators provide translational and rotational motion to node 4 of the experimental substructure with the aid of a coupler system. The experimental setup is shown in Figure 4, assembled in the Intelligent Infrastructure Systems Laboratory (IISL) at Purdue University (Condori et al., 2023). It should be noted that DOF ψ 16 , which represents the vertical displacement of node 4, has been excluded from this study, as elucidated earlier, and consequently is not depicted in Figure 4. In addition, node 7 DOFs should be measured in practice. However, in this virtual implementation, they are estimated from the ideal properties of the experimental substructure and the estimated node 4 DOFs from the displacements of the actuators.

The servo-hydraulic actuators are Shore Western 910D series, with a nominal force capacity of 9.34 kN and a stroke of ±63 mm, with a built-in linear variable differential transformer (LVDT) transducer that collects measurements of linear displacements, and two load cells (Interface, 1000 series) with a nominal force capacity of 11.2 kN, providing instantaneous force measurements. These hydraulic actuators operate with a hydraulic power supply (MTS pump) with a capacity of up to 680 L/min at 206 Bar. Meanwhile, the coupler, made from SAE 1018 low-carbon steel plates, has a mass of 17.9 kg. Furthermore, it is worth mentioning that the benchmark problem statement does not explicitly report the moment of inertia with respect to the center of mass of the coupler since it is assumed that these inertial effects were indeed incorporated when they conducted system identification of the control plant (i.e., transfer system connected to the test specimen).

In this study, the recursive least squares adaptive compensation (RLS-AC) initially proposed by Galmez and Fermandois (2022) for maRTHS systems is considered for application in the context of the maRTHS benchmark problem to provide a detailed comparison in terms of its performance and robustness with the benchmark results. The compensation is employed in a decentralized manner for reference tracking of multiple actuators in maRTHS. In this way, each actuator is treated as an independent control plant, and the goal of each compensator is to minimize the tracking error between a target signal x t and the measured signal x m for a single actuator, i.e., e = x m − x t ≈ 0 .

The RLS-AC technique consists of the design of two fundamental components: (i) a feedforward compensator, responsible for reference-tracking performance; and (ii) a recursive parameter estimator, which aims to provide sufficient robustness to the closed-loop system under unwanted disturbances and measurement noises. Both the feedforward and parameter estimators are designed independently, based on the nominal control plant model defined in Section 2.3. A schematic representation of the RLS-AC architecture employed in RTHS is presented in Figure 7.

In this study, three discrete-time signals are considered in the RLS-AC architecture: (i) target signal x t i coming from the numerical substructure; (ii) command signal x c i that is generated from a feedforward controller and applied to the control plant; and (iii) measured signal x m i from the control plant, including the specimen-actuator interaction and sensor noise. Note that the relationship between continuous and discrete time is given by t i = i T s , where i is the discrete-time index and T s = 1 / F s is the sampling time. The adaptation process is made through four adaptive parameters represented by the vector F F i = f f 0 i   f f 1 i   f f 2 i   f f 3 i   . Thus, the RLS-AC cost function is built considering only four tapped delays of the present time signals. Indeed, sampling time and size in RLS-AC may have an influence on the performance of the controller. However, this study does not consider the effects of different sampling times and sizes on the adaptive parameters’ evolution.

Here, the command signal to the actuator can be obtained from the Eq. 11: x c i = F F i X t t a p i = f f 0 i   f f 1 i   f f 2 i   f f 3 i x t i x t i − 1 x t i − 2 x t i − 3 , (11) where x c i is the command signal at the time step i , F F i are the adaptive parameters updated during the test with a recursive least squares (RLS) method (Söderström and Stoica, 1988), and X t t a p i = x t i   x t i − 1   x t i − 2   x t i − 3   T are the tapped target displacements. The parameters F F i are updated using the adaptation law provided in Eqs 12, 13: F F i = F F i − 1 + P i − 1 X m t a p i ρ + X m t a p i T P i − 1 X m t a p i e i , (12) e i = x c i − F F i X m t a p i , (13) where X m t a p i is the tapped measured displacement, e i is the compensator error, ρ is the forgetting factor (see Section 4.2), and P i is the 4 × 4 covariance matrix, which is updated according to the difference equation in Eq. 14: P i = 1 ρ I − P i − 1 X m t a p i ρ + X m t a p i T P i − 1 X m t a p i X m t a p i T P i − 1 , (14) where I is the 4 × 4 identity matrix.

The initial adaptive parameter F F 0 and covariance matrix P 0 must be defined to initialize the RLS method. A standard least-squares method with offline test data is employed to evaluate the initial parameters according to Söderström and Stoica (1988) and Wang et al. (2020), as described in Eqs 15–18: P 0 = Φ T Φ − 1 , (15) F F 0 = Φ T Φ − 1 Φ T Y , (16) Φ = X m t a p 1 , X m t a p 2 , … , X m t a p L T , (17) Y = x t 1 , x t 2 , … , x t L T , (18) where L indicates the length of the data sequence employed for the initial calibration. To obtain the offline test data, one can evaluate and impose the response of the numerical substructure on the control plant, with or without a specimen. In addition, a fourth-order Butterworth filter with a cut-off frequency of 20 Hz is chosen for noise rejection. Note that the associated time delay only affects the adaptation process and does not affect compensation directly. The previously mentioned process is schematized in Figure 8.

The algorithm mentioned before is applied to each actuator in a decentralized manner; this means that for each actuator, there is a different compensator that works independently from each other, with unique initial parameters and adaptation processes, as shown in the Simulink diagram in Figure 9A. The independence of this compensation algorithm in maRTHS, even assuming identical actuators, is proved by Galmez and Fermandois (2022) as the parameter adaptation process is different for each one of the actuators. Figure 9B shows a block diagram of the implemented RoDeAC for virtual maRTHS.

Offline tests are required to provide the necessary data for initializing the RLS-AC method and assessing the method’s performance for prescribed desired displacement loading. Then, a Simulink implementation is used, as shown in Figure 10, where the structural response is measured in an open-loop system without the feedback of the experimental forces back to the numerical substructure. In this case, we use the multi-actuator models described in Section 2.3.2, which provides information about the control system with the numerical substructure subjected to seismic excitation. Notably, this process can be carried out even without any interaction with the test specimen. Then, the structural response is imposed on the specimen ( x t ), and the displacement response is measured ( x m ). Finally, the values of the parameters P 0 and F F 0 are obtained using Eqs 15 and 16.

Different scenarios are used to analyze the initial parameters’ dispersion for the offline calibration. Therefore, three different earthquakes (El Centro, 1940; Kobe, 1995; and Morgan, 1984) with different seismic intensity measures (20%, 40%, and 60% PGA scales) are used to calculate different P 0 and F F 0 parameters. Furthermore, we considered uncertainty in the control plant as described in Section 2.2.2, with 10 simulations for each combination of earthquake ground motion and intensity measure. The uncertainty propagation in P 0 and F F 0 initial parameters is shown in boxplots in Figures 11, 12.

Note that, under the assumption of different conditions, the values for the initial parameter values change. Moreover, in this case, different values for parameters P 0 and F F 0 are obtained for slightly different actuators modeled as described in Section 2.2.2.

Finally, the mean values are considered in this study for the initialization of each parameter and are presented in Eqs 19–22.

• Actuator 1:

The forgetting factor is used to weaken the influence of older data on the estimated parameter values. This means that as simulation time increases, the old data are discarded exponentially by the parameter estimator of the adaptive compensator (Ioannou and Sun, 2012). A small forgetting factor decreases the impact of older data on the identified parameters and is thus more suitable for time-varying systems. However, small forgetting factors may cause serious fluctuations in the identified parameters (Wang et al., 2020). To observe how this parameter affects our virtual maRTHS, graphs illustrating the J 2 and J 5 error indicators (defined in Appendix A of the Supplementary Material) are presented in Figure 13. As reflected in these graphs, the lowest error is achieved with a fixed value of the forgetting factor ρ = 1 since lower values seem to cause instabilities and fluctuations in this approach, indicating that the identified parameters depend heavily on the older data. Therefore, a forgetting factor ρ = 1 is adopted in this study.

To study the tracking performance of the RoDeAC in the time domain, the seismic excitation chosen was the 1940 EL Centro earthquake, with 40% scaled PGA intensity. Figure 14 shows the target ( x t ) and measured ( x m ) displacements. It can be observed that both target and measured displacements have an excellent match at all simulation times, especially for significant time reversals. In the same way, Figure 15 shows the estimation RTHS accuracy of the interface node. Small differences between calculated and actual signals can be appreciated, which could be attributed to the non-linear nature of the equations for the coupler coordinate transformations or a deficiency in the proposed algorithm. More information on the accuracy of the RoDeAC (measured compared to reference signals) can be found in Appendix B of the Supplementary Material.

Furthermore, Figure 16 presents the evolution of parameter F F for each actuator. These results provide evidence of different evolutions in decentralized adaptive compensation and fast adaptation at the beginning of the earthquake (i.e., after 5 [sec]), allowing us to estimate the plant dynamics and achieve reasonable compensation.

To evaluate the accuracy and efficacy of the tracking control, 10 performance indices were defined in the benchmark problem (Condori et al., 2023), which are divided into two groups: (1) the tracking control performance indices and (2) the global RTHS experiment performance indices. The first three criteria ( J 1 , J 2 , and J 3 ) are measures of the reference tracking properties of the chosen controller, comparing the measured and target displacements. On the other hand, criteria J 4 to J 6 compare the target and commanded displacements. Meanwhile, criteria J 7 to J 10 establish the error between the measured displacement responses of the hybrid system, compared to the reference structure, at all floor locations of the three-story moment frame. The performance indices J 1 and J 4 have units of time (ms), while the rest are in percentages. Lower scores on all performance indices are desired. Due to the MDOF characteristic of the problem, each performance index is a vector and could be related to the number of the actuator or the node at which it is evaluated. All the equations for the performance indices are included in Appendix A of the Supplementary Material.

Table 2 summarizes the mean calculated indices ( J ¯ ) of 100 consecutive experiments of the proposed method compared with the values declared in the benchmark problem. It is assumed that 100 realizations are a sufficient sample space for a high return rate and to ensure higher robustness on the problem. Moreover, it should be noticed that some differences were encountered between the uncertainty declared on the benchmark problem and the actual MATLAB codes. Therefore, high instabilities are obtained using the benchmark controller with the stated uncertainty.

The median value of the performance indices of 100 consecutive experiments using RoDeAC, along with the interquartile range, for the El Centro earthquake, scaled to 40% of the PGA intensity, and compared with 100 consecutive experiments with the benchmark controller, are plotted in Figure 17A. Tracking control and estimation indices show that RoDeAC has lower errors than the benchmark controller. For the global RTHS accuracy, the normalized and peak errors are lower at the controlled levels but slightly higher at the upper levels. Along with its high level of robustness, RoDeAC demonstrates excellent overall performance. More information on the performance comparison of the RoDeAC with the benchmark controller for other ground motions can be found in Appendix B of the Supplementary Material. Finally, Figure 17B summarizes the performance indices for the proposed RoDeAC, showing the controller’s stable and consistent behavior for the three chosen earthquakes and with 100 different perturbations of the control plant. In particular, the worst performance indices were obtained from the Morgan earthquake with 40% PGA.

On the other hand, the synchronization subspace plots (SSPs) of displacements comparing the benchmark and RoDeAC controllers are shown in Figure 18 for all three earthquakes as inputs of the systems scaled to 40% of the PGA intensity. This visualization allows us to directly compare the system’s synchronization, where a straight line with a 1:1 slope represents a perfect synchronization. The plots show that the best synchronization is achieved for the RoDeAC controller in all cases with thinner SSPs. Nevertheless, the Cartesian SSPs show slightly worse tracking performance compared with the actuator SSPs, especially for the rotation degree of freedom.