Consider a reference structure subjected to a ground acceleration as shown in Figure 1. The reference structure is separated into a numerical substructure and a experimental substructure. The boundary conditions of the common degrees of freedom at the interface between substructures are imposed by actuators over the experimental substructure. In terms of equations, the substructuring is described as follow. The equation of motion (EOM) of the reference structure (subscript r) is expressed as:

where M_r, K_r, and C_r are the mass, stiffness, and damping matrices, respectively. x(t), ẋ(t), and ẍ(t) are the displacement, velocity, and accelerations vectors, respectively, all measured relative to the ground motion. ü_g(t) is the ground acceleration, and Γ is the seismic influence vector.

The reference structure is separated into a numerical (subscript n) and a experimental (subscript e) substructures, as shown in Equation (2):

Then, the EOM of the numerical substructure can be rewritten as shown in Equation (3), including the earthquake equivalent forces of the experimental substructure:

where x_n(t), ẋ_n(t), and ẍ_n(t) corresponds to the displacement, velocity, and acceleration of the numerical substructure, all relative to the ground motion. Meanwhile, F_e corresponds to the feedback forces from the experimental substructure described in Equation (4):

where x_e(t), ẋ_e(t), and ẍ_e(t) corresponds to the experimental displacement, velocity, and acceleration of the experimental substructure, and in the ideal case, the actuators imposes x_e = x_n. However, in an RTHS test, the experimental displacements are not necessarily equal to the numerical displacements due to synchronization errors introduced by the transfer system and interaction with the experimental substructure (Dyke et al., 1995).

Consider that the numerical substructure is an isolated system subjected to two external forces: (i) earthquake equivalent forces; and (ii) experimental feedback forces. Assuming that an appropriate numerical integration method is implemented, so that the calculated numerical substructure response is numerically stable (Chen et al., 2009; Bas and Moustafa, 2020), the equation of motion of the numerical substructure can be expressed as an energy balance, allowing to focus the stability analysis on the interaction between hybrid substructures due to experimental errors. If numerical substructure matrices M_n, C_n, K_n are symmetric;the energy balance can be obtained by taking the inner product with an infinitesimal numerical displacement trajectory, dx_n, and then integrating both sides of Equation (3) over the displacement trajectory:

Equation (3) can be later expressed as a scalar equation as shown in Equation (6):

where E_k, E_d, and E_s are the kinetic, dissipated and strain energy of numerical substructure, respectively, as described in Equation (7):

meanwhile, W_I is the work done by the earthquake forces (namely as input work), and W_F is the work done by the experimental forces (namely as feedback work):

It is worth mentioning that the energy balance from Equation (6) must be fulfilled during a test even if the experimental force is delayed, and the hybrid system is unstable. In this sense, the feedback force is just another external excitation supplied over the numerical substructure.

When a hybrid system becomes unstable, the mechanical energy of both substructures will grow exponentially even if the earthquake signal goes zero. This issue could happen if the experimental substructure introduces additional energy to the numerical substructure. Horiuchi et al. (1999) demonstrated that delay causes negative damping in the hybrid system, which is equivalent to adding external energy to the hybrid system. It is desirable to detect this problem before the mechanical energy grows up at levels that can damage the experimental equipment.

During an RTHS test, it is expected that the mechanical energy of the numerical substructure increases, because the earthquake forces introduces energy to the system, but the mechanical energy should not be greater than the input work. Reordering Equation (6) and grouping E_k and E_s in mechanical energy E_mec, Equation (9) is obtained:

If the feedback work, W_F, gets bigger than the dissipated energy E_d, the right side of Equation (9) becomes negative. Due to the energy balance, the left side becomes negative and indicates that mechanical energy gets bigger than the input energy. Thus, when the hybrid system's inherent dissipated energy is insufficient to counteract the added energy from the feedback force, the system will become unstable.

Therefore, the proposed analysis consists in comparing the feedback work with the dissipated energy through an indicator described in Equation (10):

where SW is called the Stability Warning Indicator, and H(·) is the Heaviside function used to compute SW only when W_F is positive. Notice that the dissipated energy is always strictly positive (i.e., E_d > 0); hence, E_d − W_F could be negative only if W_F is positive. Meanwhile, C_SW is a constant chosen to prevent large values of SW when the test starts (i.e., when E_d is close to zero). Notice that E_d = 0 at t = 0 (i.e., before the earthquake starts), so C_SW should be a constant selected that is bigger enough to make the denominator to be non-zero in the SW and to prevent large values of SW if W_F is dominated by noise in the force measurements. Moreover, C_SW should be small enough, so the SW denominator would be dominated by the E_d term during and after the earthquake. As a recommendation, both conditions could be satisfied if C_SW is selected as two magnitude orders below the maximum input work over the structural system. It is worth mentioning that this input work can not be calculated before the test since the complete hybrid system's response is unknown. However, an offline estimation of the input work magnitude order is enough to select a proper C_SW.

The explanation of the SW values is provided herein. When SW = 0%, due to negative values of W_F, there is no risk of instability because the experimental substructure is removing energy from the numerical substructure. Next, if 0% < SW < 100%, the experimental substructure is adding energy; but, the damping of the numerical substructure is enough to counter this effect. Finally, if SW ≥ 100%, then the numerical substructure's mechanical energy overcomes the input work, causing instability.

The main difference between EEI and SW indicators is that the former evaluates simulation accuracy when affected by synchronization errors, and the latter is focused on assessing the stability of the hybrid system. For example, if a hybrid system presents a significant delay in the transfer system, the experimental errors will cause the experimental substructure to add energy to the numerical substructure. Nevertheless, if the numerical substructure has enough damping to maintain the system stable, the SW will be <100%; meanwhile, EEI will show large values due to the synchronization errors produced by the delay.

The virtual RTHS benchmark problem from Silva et al. (2020) is considered as a practical example to prove the capacity of the proposed indicator to detect instability. This problem consists of a three-story moment frame with three lateral degrees of freedom as the reference structure (see Figure 2). This structure is separated into a 3DOF numerical substructure and an SDOF experimental substructure.

The reference structure properties are: mass per floor m = 1, 000 kg; natural frequencies f = [3.61, 16.00, 38.09] Hz; and modal damping ζ = 3% for each mode. These properties corresponds to Case IV from Silva et al. (2020). The reference structure is partitioned as described in Equation (11):

where M_e = diag(m_e, 0, 0); C_e = diag(c_e, 0, 0); and K_e = diag(k_e, 0, 0). The structural properties of the experimental substructure are m_e = 29.1 kg; c_e = 114.6 Nsec/m; and ke=1.19×106 N/m.

The reference structure and the numerical substructure are implemented in the state-space form using Simulink for direct integration. The solver utilized is a 4th Order Runge-Kutta with a fixed time step of Δt = 1/4, 096 s. The reference structure model is presented in Equations (12) and (13) with ground acceleration as input and reference displacements as output.

The numerical substructure model is presented in Equations (14) and (15) with ground acceleration and experimental force as inputs and numerical displacements, velocities, and accelerations as outputs.

where γ = [1;0;0]^T and f_e is experimental force calculated as shown in Equation (16):

and x_m, ẋ_m, and ẍ_m are the measured displacement, velocity and acceleration imposed on the experimental substructure.

To compute the Stability Warning Indicator (SW) it is necessary to evaluate different energy terms using numerical integration. To calculate the accumulated energy E[j] of a force P[i] over a displacement increment Δx[i], ∀i ∈ [1, j], a trapezoidal rule is formulated in Equation (17):

This integration method is implemented in Simulink, as shown in the block diagram of Figure 3, where Discrete Finite Impulse Response (FIR) Filters and Product blocks are employed to obtain the area of each trapezoid, and a Discrete-time Integrator block is required to obtain the accumulation over time. In the case of W_F the input force is P = f_e and the input displacement is x=xn(1) (scalar signals); meanwhile, to evaluate E_d, the input forces are P = C_nẋ_n and the displacements are x = x_n (vector signals).

Additionally, to compare the proposed stability indicator, the EEI (Ahmadizadeh and Mosqueda, 2009) is calculated for each simulation. The EEI can be obtained as follow:

where Ek*, Ed* are the kinetic and dissipated energy calculated with the first derivative of displacement instead of velocity of numerical substructure. EEmax is a nonzero elastic energy term utilized to normalize EEI when W_I is close to zero at the beginning of the test. It is important to declare that non-compensated delay implies added energy to the system, resulting in negative values of EEI. W_exp is the work done by experimental force over the measured displacement as shown in Equation (19):

In this study, the reference structure is subjected to the Kobe earthquake with a peak ground acceleration (PGA) scaled to 40%, and the maximum input work from the earthquake is around 100 Joules. This value can be computed from the integration of the reference structure. In a real RTHS experiment, the magnitude order of the input energy can be estimated from the ground motion and global properties of the structure, such as total mass and approximated damping and stiffness. Therefore, the denominator constants are selected considerably smaller, such as 1% of the maximum input work:

This value is bigger enough to prevent large values for each indicator at the beginning of the earthquake. Also, these constants will have a negligible effect on each indicator during and after the earthquake, since C_SW will be relatively small compared with E_d in the denominator of Equation (10), just like EEmax will be smaller than W_I in the denominator of Equation (18).

Before conducting the RTHS simulations with different delay values, an estimation of the critical delay is obtained with the method proposed by Gao and You (2019). This method consists in constructing the EOM of the hybrid system with some assumptions synchronization errors (delay and amplitude errors). The EOM of the hybrid system has the following equivalent mass, damping, and stiffness:

where Φ is the modal matrix of the reference system, ω is a diagonal matrix with the natural frequencies of the reference system, δ is a diagonal matrix with the delay for each mode, and Δ is a diagonal matrix with the amplitude error for each mode. With the equivalent properties of the hybrid system, the state matrix A_eq is constructed as shown in Equation (24):

The hybrid system is stable if the matrix A_eq has eigenvalues with negative real part; thus, the critical delay can be obtained by searching the minimum constant delay that produces eigenvalues with positive real part. For the partitioning of this problem, assuming Δ = I and δ = diag(τ_cr; τ_cr; τ_cr), the critical delay is:

This critical delay serves as a reference for the simulations with constant delay. It is expected that simulations with delays smaller than this critical value would remain stable while running the simulations. In contrast, the opposite will result in an unstable response.

The Simulink block diagram for the simulations with constant delay is presented in Figure 4, where the output of the numerical substructure, denoted as “target vector,” is defined as xt=[xn(1);x˙n(1);x¨n(1)]T (displacement, velocity and acceleration of first degree of freedom). Then, the Delay block generates a discrete delay of the input signal x_t by a specified number of samples d, so the time delay τ is a multiple of the simulation time step. Then, the delayed boundary conditions are:

Finally, the outputs of experimental substructure are described in Equations (27) and (28):

Four simulations are performed with different delay values from 0 to 9 ms. The structure is subjected to the Kobe earthquake with a peak ground acceleration (PGA) scaled to 40%. The measured displacement of experimental substructure x_m is compared with the first degree of freedom displacement x_r1 from the reference structure, as shown in Figure 5. Measured displacements are very close to reference displacements for lower values of delay (τ = 0 ms and τ = 0.98 ms). Then, for higher delay τ = 7.1 ms, a notorious synchronization error in displacements appear, but at least the measured displacement decays at the end of the simulation. This does not happen with τ = 9msec > τ_cr, where the measured displacement grows exponentially after the earthquake reaches zero.

In Figure 6, the error of displacement x_m respect to x_r1 for each delay is presented, including the normalized root mean square error (NRMSE) for each case. Clearly, for τ = 0 ms and τ = 0.98 ms, the system is stable; while, for τ = 9 ms is unstable (i.e., error grows exponentially). However, for τ = 7.1 ms, the error is considerable, with NRSME = 144%, but the system is still stable. It is easy to classify the case with τ = 9 ms as unstable after the test. Still, it would be beneficial to detect unstable behavior before the system reaches a large response.

In Figure 7, the SW indicator is calculated for each model scenario. For τ = 0 ms, τ = 0.98 ms, and τ = 7.1 ms, the indicator stays under 100% during the entire test, showing the stable behavior of each simulation. In contrast, for the unstable case with τ = 9 ms, the indicator exceeds 100% near 10 s of simulation, much earlier than the instant when the system reaches the large displacements presented in Figure 5. Thus, the SW indicator could help stop the simulation in time, avoiding any dangerous behavior in the RTHS test.

For the sake of comparison, Figure 8 shows the absolute value of EEI for each case scenario. The EEI grows up with higher values of delay but is not evident in which case is stable and in which it is not. For τ = 0 ms, the EEI values are relatively small. Moreover, for τ = 0.98 ms, the EEI stays under 15% during the test; both cases are indeed stable. However, for τ = 7.1 ms, the behavior of EEI is similar to τ = 9 ms before 15 s of simulation, so it is not possible to detect instability until the earthquake signal dies out. If the EEI is utilized to stop unstable testing, it is quite challenging to define a threshold that relates EEI with instability. For example, for a threshold of EEI = 100%, the simulations with τ = 7.1 ms and τ = 9 ms should be stopped before 10 sec, even when the simulation with τ = 7.1 ms is stable and does not present any risk of large displacements. Nevertheless, it must be recognized that EEI is an excellent tool to assess the accuracy of the results from an RTHS test through a post-processing analysis.

In this subsection, the same partitioning as in the previous subsection is utilized, but in this case a transfer system is modeled as presented in Silva et al. (2020), and adaptive compensation with different design parameters is implemented to minimize the errors produced by the transfer system. The block diagram utilized to model this problem in Simulink is presented in Figure 9, where the output of the numerical substructure is the target displacement xt=xn(1) (displacement of first degree of freedom). Then, the target displacement is sent to the adaptive controller to determine the commanded displacement for the actuator. The plant (i.e., actuator connected to the experimental substructure) is modeled with the transfer function presented in Equation (29), which includes actuators parameters and interaction with the experimental substructure, and relates the commanded displacement x_c with the measured displacement x_m:

where s = iω, is the Laplace variable, i is the complex number and ω is the circular frequency. This transfer function has a frequency-dependent phase. The phase is another way to state that the measured and commanded displacements are delayed. Thus, the delay is frequency-dependent for stable and underdamped linear systems.

The interaction forces are obtained with Equation (30):

Besides, sensor noise has been added to the measured displacement x_m and the experimental force f_e using band-limited white noise (BLWN) blocks to model the physical sensors. Then, the signal x_m is sent to the controller for the adaptation process. Consequently, the signal f_e is fed back to the numerical substructure to be considered in the integration of the numerical EOM. Notice that force measurement noise could impact the experimental force power, but it has a small effect on the work done by this force due to the numerical integration process to calculate work. However, measurement noises can affect the hybrid loop, depending on the structural properties, transfer system, and dynamic compensation.

The adaptive compensation implemented in the controller consist in adaptive model-based compensation (Chen et al., 2015) with the modifications presented in Fermandois et al. (2020). A first-order adaptive feedforward determines the commanded displacement as shown in Equation (31):

where a₀ and a₁ are the adaptive control parameters, and ẋ_t is approximated using the backwards difference method. The adaptation laws of a₀ and a₁ are described in Equation (32):

where Γ = diag(Γ₀, Γ₁) is the adaptive gain matrix associated with the adaptation rate of parameters a₀ and a₁, respectively, and e is an estimation error of the adaptive parameters described in Equation (33):

Both signals x_c and x_m are filtered with a fourth-order Butterworth low-pass filter with a cutoff frequency of 20 Hz, and ẋ_m is obtained with the backward difference method.

This subsection aims to demonstrate the performance of the SW stability indicator for different compensation scenarios. Notice that the goal of the SW indicator is to detect instability during a test with a predefined compensation design, and not to tune the compensation/adaptation parameters. In these simulations, the initial adaptive parameters are selected as a₀(t = 0) = 1 and a₁(t = 0) = 10/1, 000 s. With these parameters, the controller cannot compensate for the control plant delay successfully; therefore, parameter adaptation is necessary.

Four simulations are performed with different adaptive gains presented in Table 1, where Γ_a corresponds to a non-adaptive case, Γ_d to an optimally calibrated case (Fermandois et al., 2020), and Γ_b and Γ_c are intermediate cases of adaptation. The structure is subjected to the Kobe earthquake scaled to 40% PGA.

The measured displacement of experimental substructure x_m is compared with the first degree of freedom displacement x_r1 from the reference structure in Figure 10, where the case with Γ_a is unstable, while the case with Γ_b presents large measured displacements, but that decreases after the earthquake. Considering that in true RTHS tests, the reference displacement is not available, it would be difficult to distinguish if the measured displacements are correct or are associated with an unstable response. On the other hand, for Γ_c and Γ_d simulations are stable. In Figure 11, the error of displacement x_m respect to x_r1 for each adaptive gains is presented, including the normalized root mean square error (NRMSE) for each case. The cases with Γ_a present unacceptable errors, while for Γ_b, the errors are substantial, but at least decrease after the earthquake. For the cases with Γ_c and Γ_d, the errors are relatively small.

The difference between each case is explained by the adaptation process presented in Figure 12, where there is no adaptation for Γ_a, slow adaptation for Γ_b, and the same convergence for Γ_c and Γ_d but with different rates. The parameter a₁ is directly related to the delay compensation, so it is expected that for Γ_a and Γ_b, the non-compensated delay produces large errors in the simulations. However, the adaptation for Γ_b is apparently enough to stabilize the system during the test.

In Figure 13, the SW indicator for each case is presented. For Γ_c and Γ_d, the indicator is growing at the beginning of the earthquake. Then, the SW indicator decays due to the adaptation and good compensation, and stays under 100% during the entire simulation, confirming the stability of these simulations. For Γ_a and Γ_b, the indicator reaches 100% near 10 s of simulation, so much before each simulation presents large displacements. So, the stability warning indicator could help stop the simulation in time before any damage occurs in the test setup. In the particular case with Γ_b, the stability warning indicator exceeds 100% between 10 and 20 s, consistent with the large response of the system in that interval. After 20 s of simulation, when the controller has reached enough adaptation to compensate for the delay, the SW indicator decays below 100%. This result is consistent with the stable behavior of this case at the end of the simulation. However, this simulation could be stopped when the indicator reaches 100% to avoid the behavior shown between 10 and 20 s.

Moreover, notice that for Γ_b and Γ_c cases, the SW does not decay to zero at the end of the simulation. This situation occurs because of uncompensated delay during the earthquake, which causes a surplus of energy inserted into the numerical substructure. Later, when the delay is compensated, displacements and forces are too small, so the work done by the experimental forces remains constant. This effect does not happen for Γ_d, where the delay is compensated early when the earthquake strikes. The experimental substructure dissipates energy during the rest of the earthquake, leading W_F to zero and consequently, SW to zero as well.

Furthermore, the control plant in the benchmark problem is modeled as a linear time-invariant system, but in these simulations, the adaptive controller produces time-varying compensation. Therefore, the delay between target and measured displacements is a time-varying property. This scenario is a complicated situation to predict the stability of a system before the test. However, with the online analysis of the SW it is possible to detect which case results in stable behavior and which are not.

For comparison, Figure 14 shows the absolute value of EEI for each adaptive case scenario. For Γ_a and Γ_b, the EEI presents large values, while for Γ_d, the EEI keeps low during the test. So, the EEI could allow lab technicians to stop the simulations with Γ_a and Γ_b, and accept the simulation with Γ_d. However, for Γ_c case, the EEI reaches 100% before 10 s where the displacements are small, and the controller is adapting, and the test still has good chance to give acceptable results without the risk of instability. Thus, the EEI is an excellent tool to analyze the errors after the RTHS test, but it is not reliable for detecting instability online.

Additionally, in these simulations, displacement and force measurement noises do not have a notorious impact on the hybrid system response. The numerical substructure acts like a low-pass filter for the force measurement, and the adaptive compensator is not affected by the displacement noise.