For the NARX model (Leontaritis and Billings, 1985), the current value of a time series is influenced not only by its past values but also by the current and past values of an exogenous series. Mathematically, this relationship can be expressed as: y t = F ϖ t + ε t (1a) ϖ t = y t − 1 , … , y t − n b , x t , … x t − n a T (1b) where y t is the output of the model at the tth step; ϖ t is the regression term; y t − 1 , … , y t − n b are the previous values of output before the tth step; x t , … x t − n a are external input for current and previous values; ε(t) is error term; and F ( ) represents nonlinear function. When employing the NARX model in Eq. 1b to formulate compensation method for delay, the predicted displacement d n p and command displacement d n c for the nth step are employed as its input and output, without the inclusion of an error term, denoted as x(t) and y(t) respectively. d n p = F d n − 1 p , ⋯ d n − i p , d n c , ⋯ d n − j c (2) where i and j are the number of predicted displacement and command displacement required. Derived from Eq. 2, the predicted displacement at the current time is computed based on predicted displacements at previous times and command displacements at the current and previous times. This resulting is referred to as NARX model-based compensation hereafter.

Eq. 2 presents a universal formulation of NARX model-based compensation, wherein the predicted displacements are interconnected, allowing for the application of any nonlinear function. However, under practical circumstances, this formulation of NARX model-based compensation may not be directly applicable. A simplified version is often necessary in such situations.

In contrast with calculated and measured displacements, predicted displacements serve as intermediate values and lack a direct correlation with structure responses in RTHS. When predicted displacements are interrelated, it poses a risk to the stability of RTHS, particularly if one predicted displacement falls outside a reasonable range. To address this issue, the predicted displacements from previous steps are ignored.

Among all possible types of functions F ( ), the efficacy of polynomial functions has been substantiated. In this context, Eq. 2 can be represented as: d n p = ∑ k = 0 j α k d n − k c + ∑ p = 0 j ∑ q = p j α l , m d n − p c d n − q c + ⋯ (3) where α k and α l , m are the corresponding coefficients of the NARX model for the first order and second order polynomial functions, respectively. As observed from Eq. 3, an infinite number of terms can be employed in the NARX model-based compensation method, representing infinite orders of polynomial functions. For the purpose of practical applications, Eq. 3 should be truncated to a certain order with finite terms. For example, Eq. 3 can be simplified as d n p = ∑ k = 0 j α k d n − k c (4a) when the order of the NARX model is two, and d n p = ∑ k = 0 j α k d n − k c + ∑ p = 0 j ∑ q = p j α l , m d n − p c d n − q c (4b)

When the order of the NARX model is one, Eq. 4a represents the simplest linear form among all possible NARX model-based compensations when compared with Eq. 4b. In this configuration, the NARX model degenerates into an ARX model. Eq. 4a is also referred to as the jth order NARX model based on the terms required. Table 1 outlines the coefficients for first, second, and third order NARX model-based compensations, denoted as F-, S-, and T-NARX, respectively (Xu et al., 2022).

In theory, the measured displacement d n m should be identical to the command displacement d n c when the delay compensation method works perfectly. In this scenario, the relationship between predicted displacements and measured displacements can be expressed as: d n p = F d n − 1 p , ⋯ d n − i p , d n m , ⋯ d n − j m (5a)

or d n − 1 p = F d n − 2 p , ⋯ d n − i − 1 p , d n − 1 m , ⋯ d n − j − 1 m (5b)

Given that all d p and d m values in Eq. 6b are accessible during RTHS, the NARX model can be derived through data regression between predicted and measured displacements. This model can then be applied to Eq. 2 to compensate time delay. For illustrative purposes, Eq. 4a is employed as an example to illustrate this procedure. By substituting d n c with d n m , Eq. 5a can also be written as d n p = ∑ k = 0 j α k d n − k m (6a)

or d n − 1 p = ∑ k = 0 j α k d n − k − 1 m (6b)

However, Eq. 6b comprises only one equation with j+1 undetermined parameters, leading to an under-determined system. Consequently, there exists an infinite number of possible solutions satisfying this equation. To address this challenge and arrive at an appropriate and unique solution, additional constraints must be introduced. In this study, the predicted displacements from previous steps are employed to estimate the undetermined parameters, i.e., d n − i p = ∑ k = 0 j α k d n − k − i m   i = 1 , 2 , … M (6c)

The number of equations, denoted as M, in the parameter estimation process is thereafter referred to as the window length, which is an integer larger than j. It is then obvious from Eq. 6c that the calculated coefficients α k in NARX-based compensation are not instantaneous undetermined parameters rather parameters estimated through a fixed window length. When the window length is relatively small, the estimated parameters closely approach the current values. Defining α n = α 0 , α 1 , ⋯ , α j T , d n p = d n − 1 p , d n − 2 p , ⋯ , d n − L p T , d n m = d n − 1 m , d n − 2 m , ⋯ , d n − L m T and D= d n m , d n − 1 m , ⋯ , d n − j m , Eq. 7 can be written as d n p = D α n (7)

In Eq. 7, d p n and D can be obtained during RTHS, allowing for the estimation of parameters in NARX model-based compensation through these linear systems of equations. As the window length, denoted as L, is larger than j, Eq. 7 leads to an over-determined system, wherein the number of constraints exceeds the number of unknown variables. Consequently, a solution satisfying the linear system is generally nonexistent. Instead, an approximate solution is sought to minimize a given error. To obtain the estimated parameters, two data regression methods, namely, ordinary least squares and ridge regression, are introduced in this study.

Eqs 4a, 4b employ the linear and nonlinear formulation, respectively. The former has exhibited better applicability in real-time hybrid simulation compared to the latter, attributed to the inherent complexity and nonlinearity of hydraulic systems. The linear form can be adapted to various hydraulic systems, while the same cannot be said for the nonlinear form. In this section, the selection of nonlinear terms is very briefly discussed. The nonlinear terms are selected to serve as fine-tuning mechanisms to ensure the stability and accuracy of NARX model-based compensation composed of linear terms. The coefficients of the nonlinear terms are significantly smaller than those of the linear terms, thereby maintaining stability. It should also be noted that identifying suitable nonlinear terms for a specific hydraulic system is a challenging task, even with some preliminary knowledge. Nonetheless, a well-chosen nonlinear system can yield superior compensation performance, capitalizing on the inherent nonlinearity of the hydraulic system.

A novel maRTHS benchmark model has been proposed by Condori Uribe et al. (2023) focusing on a frame subjected to seismic loading at the base. As illustrated in Figure 1, the benchmark model features a steel moment-resisting frame with three bays and three stories. The mid-span at the bottom is designated as the experimental substructure, while the rest is considered as the numerical substructure. To simulate seismic loading, a scaled El Centro historic record serves as the input ground motion, with a scaling factor of 0.40 is applied to maintain linear elastic behavior in all structural components of the frame. For the experimental substructure, two actuators are strategically positioned to deliver equivalent translational and rotational motion. A coupler is employed to facilitate the coupling of the linear stroke of both actuators. A block diagram illustrating the key components, along with a properly designed and tuned control system, is presented in Figure 2A. The objective of this benchmark is to formulate a control system that ensures the output of the controlled plant accurately tracks the target displacement vector. In the benchmark model, an estimator is essential to filter out high-frequency noise before being sent to the control plant. However, this block is typically integrated into the data acquisition module for real RTHS. Consequently, the estimator is not employed for NARX model-based compensation in the benchmark model, as depicted in Figure 2B. Instead, a fifth order Butterworth filter with a cutoff frequency of 6 Hz is utilized in constructing the NARX model to balance the frequency components of the actuator response and noise mitigate performance. The time delay introduced by Butterworth filter is compensated simultaneously, since the displacements after Butterworth filter are sent to the compensator.

To evaluate the performance of NARX model-based compensation for maRTHS, various compensation methods, including traditional inverse compensation (IC), NARX compensation without nonlinear terms (NARX ^l ), NARX compensation with ridge regression (NARX ^r ), and NARX compensation with different nonlinear terms (NARX ⁿ ) are applied to the benchmark model.

The evaluation of compensation performance involves assessing the difference between measured and calculated displacements. However, tracking performance has limitations in evaluating maRTHS. To comprehensively gauge the effectiveness of various compensation methods, ten distinct evaluation criteria (J ₁∼J ₁₀) are employed. The initial six criteria focus on evaluating the tracking performance of the control system, while the remaining four concentrate on assessing the global performance of RTHS. Specifically, the ten evaluation criteria are defined as: J 1 , i = argmax r ∑ k = 1 N η n s , i k ⋅ η m , i k − r × 1000 / f s , i = 1 , 2 , (11a) J 2 , i = ∑ k = 1 N η m , i k − η n s , i k 2 ∑ k = 1 N η n s , i k 2 × 100 , i = 1 , 2 . (11b) J 3 , i = max η m , i k − η n s , i k max η n s , i k × 100 , i = 1 , 2 . (11c) J 4 , i = argmax r ∑ k = 1 N η n s , i k ⋅ η ^ m , i k − r × 1000 / f s , i = 1 , 2 , (11d) J 5 , i = ∑ k = 1 N ψ ^ m , i k − ψ n s , i k 2 ∑ k = 1 N ψ n s , i k 2 × 100 , i = 4 , 28 , (11e) J 6 , i = max ψ ^ m , i k − ψ n s , i k max ψ n s , i k × 100 , i = 4 , 28 . (11f) J 7 , i = ∑ k = 1 N ψ ^ m , i k − ψ i k 2 ∑ k = 1 N ψ i k 2 × 100 , i = 4 , 28 , (11g) J 8 , i = ∑ k = 1 N ψ n s , i k − ψ i k 2 ∑ k = 1 N ψ i k 2 × 100 ,   i = 2 , 26 , 3 , 27 . (11h) J 9 , i = max ψ ^ m , i k − ψ i k max ψ i k × 100 , i = 4 , 28 . (11i) J 10 , i = max ψ ^ n s , i k − ψ i k max ψ i k × 100 , i = 2 , 26 , 3 , 27 . (11j) where i represents the actuator number in Eqs 11a–11d, and represents the node freedom number in Eqs 11e–11i. In Eqs 11e–11i, ‘4/2/3′ and ‘28/26/27′ represent the x and θ direction responses for the first/second/third floor. ηns, ηm, and η ^ m represent the desired actuator displacements, measured actuator displacements, and estimated interface node displacements; ψ ^ m , Ψns and Ψ represent the frame target displacements, estimated interface node displacements, and numerical substructure response; f _s represents the sampling frequency and is set at 1024Hz. The ten criteria above can be divided into three categories: J ₁ to J ₃, J ₄ to J ₆, J ₇ to J ₁₀ represent tracking control, estimation results, and global performance, respectively (Uribe et al., 2023). Since the estimator is not adopted in NARX compensation, J ₁ is exactly the same as J ₄. Consequently, the results of J ₄ are not presented in this study.

The performance of traditional inverse compensation can be assessed through Case 1 and 2, with the reference and the measured displacements of the first floor depicted in Figure 3 for both cases. It is evident that the measured displacement for Case 2 is closer to the reference, indicating that Case 2 outperforms Case 1. The evaluation criteria are presented in Table 6 for inverse compensation. It can be observed that all criteria for Case 2 are superior to those for Case 1, which is consistent with observations in Figure 3. Given that a ₀ and a ₁ remain time-invariant throughout the simulation, the performance of inverse compensation is highly dependent on the initial estimation. From J _1,1 and J _1,2 for Case 2, inverse compensation proves effective in reducing time delay when the initial estimation is accurate. However, a notable time delay is observed when the initial estimation deviates from actuator dynamics, as observed in J _1,1 for Case 1.

An intriguing observation arises when comparing the tracking criteria of the second actuator, namely, J _1,2, J _2,2, and J _3,2, for the two cases. Despite employing the same initial estimations for both cases, the tracking criteria for Case 2 surpass those for Case 1. This can be attributed to the interaction between the two actuators in maRTHS. Consequently, full compensation of dynamics is imperative all actuators in maRTHS. This might stem from inadequate compensation altering the properties of the structural response, particularly the frequency components. Hence, compensation parameters that were initially suitable may become unsuitable. However, this issue is primarily associated with constant delay compensation methods, and it can be effectively addressed through adaptive compensation.

When comparing tracking criteria of the same actuator, namely, J ₁, J ₂, and J ₃ for the same case, the performance of the second actuator is significantly better. In Case 1, the substantial difference between the two actuators is attributed to the initial estimation. However, both J ₁ values are zero for both actuators in Case 2, while J ₂ and J ₃ are nearly three times higher for the first actuator compared to the second actuator. This discrepancy suggests that, in comparison with the second actuator, the actuator delay for the first actuator is more pronounced. In such cases, traditional constant delay compensation is unsuitable, and adaptive compensation becomes necessary. Moreover, the estimated results for x direction responses, i.e., J _5,4 and J _6,4, is reduced from about 11% for case 2%–3% for case 1. This implies the importance of initial estimation for inverse compensation. However, the estimated results for θ direction responses, i.e., J _5,28 and J _6,28, have little improvement. This however indicates the requirement of adaptive compensation.

The window length and the model order are two major factors that influence the performance of NARX model-based compensation. Cases 3 to 5 are employed to demonstrate the impact of window length on the performance of second order NARX model-based compensation without nonlinear terms. Figure 4 presents the comparison of measured first-floor displacements for these three cases. It can be observed in Figure 4A that it is evident that the measured displacement in cases 4 and 5 aligns quite well with the reference. In Case 3, the fluctuations in displacement over the last 15 s indicate instability. Additionally, in Figure 4B, none of the measured displacements exactly match the reference, suggesting that the delay compensation method can reduce but cannot eliminate time delay, even for NARX model-based compensation.

The evaluation criteria for NARX model-based compensation with different window lengths when considering tracking control and estimation results are presented in Table 7. Criteria are small and nearly identical for Case 4 and 5, while larger values criteria are observed for Case 3. Constructing the NARX model requires a certain amount of data, otherwise the NARX model may deviate from actuator dynamics. Thus, there is a limit on the minimum data length for NARX model-based compensation. When the window length is smaller than the limit, the performance of compensation deteriorates, as observed for Case 3. Once the minimum data length is reached, increasing the window length has little impact on performance, as observed from the minimal difference between Case 4 and 5. Although computational efficiency and response to time-varying delays might decrease with an increase of window length. When comparing Case 4 and 2, the NARX model-based compensation outperforms inverse compensation in all tracking evaluation criteria, i.e., J ₁ to J ₃, especially for J _2,1. Consequently the estimated results for θ, i.e., J _5,28 and J _6,28, direction responses is reduced from about 20% to 2% when NARX model-based compensation utilized.

The performance of NARX model-based compensation is determined by the parameters of the NARX model. Since the same compensation method is applied for both actuators, the parameters of the NARX model are nearly identical for Case 4 and 5. Therefore Figure 5 only presents the compensation parameters for the first actuator for Case 3 and 4. When nonlinear terms are not considered, the sum of a should be close to 1. From Figure 5, a ₀ and a ₂ are similar, while a ₁ is almost the opposite of the sum of a ₀ and a ₂ for both cases. In Figure 5A, the parameters undergo dramatic changes from 11 s to 15 s, 20 s to 25 s, and 28 s to the end of the simulation. Large value of a ₁ is observed more than 10,000 around 28 s, which corresponds to the time when measured displacement deviates from the reference and the simulation becomes unstable. On the contrary, the parameters in Figure 5B do not have a significant change, resulting in better compensation performance. This again indicates the importance of window length in NARX model-based compensation.

Case 4, 6, and seven all have window length of 1024 and are used to assess the impact of model order on NARX model-based compensation without nonlinear terms. Figure 7 displays the time history of the first-floor measured displacements for these cases compared with the reference. Instability is observed in Case 7 when a fourth order NARX model is employed, whereas Case 6, using a first order NARX model, closely resembles case 4 with a second order NARX model. Table 8 presents the criteria for Case 6 and 7, which is consistent with observations from Figure 6. This suggests that while higher order NARX models may become unstable for the same window length, stability can be maintained with lower order NARX models. Additionally, lower order models require fewer parameters, allowing for accurate parameter estimation with smaller window lengths. The compensation parameters are presented in Figure 7 for Case 6 and 7. In Figure 7A, a ₀ remains approximately 50 for the first order NARX compensation, except around 7 s and in the last 2 s. Furthermore, the sum of a ₀ and a ₁ is close to 1. In contrast, the parameters for the fourth order NARX model can reach up to 8×10⁷, resulting in significantly larger evaluation criteria in Table 8 for Case 7.

Due to smaller window length than the minimum for the fourth order NARX model, From Figure 6 and Table 8, it is evident that the measured displacement for Case 7 becomes unstable. This results in inability for high order NARX model-based compensation. It is apparent that traditional ordinary least square has limitations, especially when the number of undetermined parameters is more than three in the NARX model. Rigid regression, instead of ordinary least square, is employed in Case 8–11 to calculate parameters in the fourth order NARX model with different values of parameter γ. The time histories of compensation parameters are presented in Figure 8 and the corresponding values of evaluation criteria are summarized in Table 9.

In Figure 8A, it can be observed that Case 11 remains unstable due to γ value as small as 1 × 10⁻⁹, indicating that the rigid regression does not take effect. For very small value of γ, rigid regression is same as ordinary least square regression. With the increase of γ, cases 8–10 become stable such as case 10 when γ equals 1 × 10⁻⁷. Compared with case 4 in Table 7, the evaluation criteria in Table 9 for case 10 are better for almost all criteria. This suggests the advantages of the higher order NARX model-based compensation with ridge regression. When γ increases to 0.01, the tracking error visibly increases for the first actuator, while a slight increase is observed for the second actuator as shown in Table 9. Moreover, a significant increase can be observed for all evaluation criteria when γ increases to 1. Consequently, ridge regression would reduce the performance of compensation when γ reaches a threshold. However, the simulations remain stable for larger γ. The parameters for Case 11 are like those of Case 7, while Case 9 falls between Case 8 and 10. Thus, only the parameters of the NARX model for the first actuator for Case 8 and Case 10 are presented in Figure 9. It can be observed that the absolute value of each compensation parameter decreases with the increase of γ. Influence of nonlinear term.

Nonlinear terms are explored in this study for the NARX model-based compensation. In Figure 10A, clear instability is observed for Case 15, while the rest cases remain stable. In the close-up view of Figures 10A, B slight difference between the measured displacement and reference is observed for Case 16, while the measured displacements in the remaining cases agree well with the reference. This implies that different nonlinear term leads to distinct compensation performance. The evaluation criteria are summarized in Table 10 for the five cases with nonlinear terms. It is worth noting that Case 12 and 13, which involve x _n-2 ² with little difference from x _n-1×x _n-2, have almost identical evaluation criteria. In comparison with corresponding case 4 with only linear terms, these criteria are slightly larger. Conversely, Case 14 demonstrates improved evaluation criteria particularly with smaller tracking control when compared with corresponding linear case. This observation suggests that the inclusion of a carefully selected nonlinear term can further improve the performance of the NARX model-based compensation. However, large evaluation criteria are noticeable for Case 15 and 16, consistent with observations in Figure 10.

In the context of utilizing a second order NARX model for compensation, three key parameters (a ₀, a ₁, a ₂) are derived during the simulation. While the behavior of a ₀ and a ₁ aligns with the first order NARX model compensation as shown in Figure 7A, the focus in Figure 11 is exclusively on the time history of the third parameter a ₂. From Figure 11A, it can be observed that a ₂ exhibits slight variation ranging between −0.015 and 0.01. This magnitude is notably smaller when compared to a ₀ and a ₁ from Figure 7A. This observation underscores the nuanced and fine-tuning role played by the nonlinear term in shaping the compensation performance. In particular, the trajectories of a ₂ for cases 12, 13, and 16 remain relatively small throughout the simulation with occasional jumps. Conversely, for case 15, the absolute values of a ₂ surpass 200, indicative of significant fluctuations, ultimately contributing to worse compensation performance.

As maRTHS incorporates a minimum of two actuators, a comprehensive global evaluation becomes imperative. Given that cases 3, 7, 8, 11, and 15 have been established as unstable through simulation, the remaining 11 cases are used to evaluate the performance of NARX model-based compensation. Figure 12 presents the global performance evaluation criteria (J ₇ to J ₁₀) for these cases.

It can be observed that case 2 exhibits the best global performance among all stable cases, while case 1 shows the worst global performance evaluation criteria, with the exception of the θ direction responses for the first floor. This observation confirms previous finding that the performance of constant delay compensation significantly rely on initial estimation. Notably, the values for the second and third floors (J ₈ and J ₁₀) in case 2 are surprisingly lower than those in NARX-model-based compensation. This might be attributed to the reduction in testing errors for the x and θ directions when calculating the response of the second and third floors. Consequently, it becomes apparent that improving the tracking performance of individual actuators does not necessarily translate to an enhancement in the global performance of maRTHS, which is distinct from saRTHS. Figure 13.