Kalman Filter-Based Adaptive Delay Compensation for Benchmark Problem in Real-Time Hybrid Simulation

: Real-time hybrid simulation (RTHS) is a versatile, e ﬀ ective, and promising experimental method used to evaluate the structural performance under dynamic loads. In RTHS, the emulated structure is divided into a numerically simulated substructure (NS) and a physically tested substructure (PS), and a transfer system is used to ensure the force equilibrium and deformation compatibility between the substructures. Owing to the inherent dynamics of the PS and transfer system (referred to as a control plant in this study), there is a time-delay between the displacement command and measurement. This causes de-synchronization between the boundary of the PS and NS, and a ﬀ ects the stability and accuracy of the RTHS. In this study, a Kalman ﬁlter-based adaptive delay compensation (KF-ADC) method is proposed to address this issue. In this novel method, the control plant is represented by a discrete-time model, whose coe ﬃ cients are time-varying and are estimated online by the KF using the displacement commands and measurements. Based on this time-varying model, the delay compensator is constructed employing the desired displacements. The KF performance is investigated theoretically and numerically. To assess the performance of the proposed strategy, a series of virtual RTHSs are performed on the Benchmark problem in RTHS, which was based on an actual experimental system. Meanwhile, several promising delay-compensation strategies are employed for comparison. Results reveal that the proposed time-delay compensation method e ﬀ ectively enhances the accuracy, stability, and robustness of RTHS.


Introduction
Real-time hybrid simulation (RTHS) is a dynamic experimental method that deals with time-dependent specimens [1], and is an extension of hybrid simulation [2,3] originating from the pseudo-dynamic test [4]. In RTHS, the emulated structure, which is typically large-scale, is divided into numerical and physical portions. The former is numerically simulated by a computer program, and is usually called a numerical substructure (NS), while the latter is physically tested in a laboratory and is called a physical substructure (PS) [5,6]. Considering the cost and space requirement, RTHS has become an attractive alternative to the widely used quasi-static test and shaking table test to evaluate the performance of structures equipped with time-dependent members and/or devices [7][8][9].
A typical requirement of RTHS is that it must be executed in real time over a tiny time interval. However, owing to the computation, signal acquisition, and transmission processes, as well as the

Overview of the Benchmark Problem in RTHS
In RTHS, the emulated structure is divided into PS and NS, and they are related by a transfer system. The equation of motion for the NS is M N .. ..
where M, C, and K are the mass, damping, and stiffness matrices of the structure, respectively; the subscripts "N" and "P" represent the NS and PS, respectively; X denotes the displacement vector of the structure and Γ refers to the location vector; r is the restoring force vector and ..
x g is the ground acceleration.
Equation (1) is solved by employing a time-stepping integration algorithm, and the calculated displacement X P is forwarded to the PS to satisfy the compatibility condition using a displacement-control-based actuator. Then, the force at the boundary interface is measured directly and fed back into the NS to satisfy the equilibrium. The above procedure is repeated until the termination condition is reached.
However, owing to the dynamics of the control plant, the boundary conditions of NS and PS are often de-synchronized. Hence, in recent years, numerous approaches have been proposed and examined to treat this issue under certain conditions. To provide a unified performance evaluation standard for different methods and to advance the field of RTHS, a Benchmark problem in RTHS was established by Silva et al. [32]. In the Benchmark problem, the models of the servo-hydraulic transfer system, the sensing and control implementation hardware, control-structure interaction, and the PS are identified from the actual equipment [33] and modeled by transfer functions in MATLAB. The block diagram of the benchmark problem is depicted in Figure 1.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 3 of 19 converge to the true values by the KF, and there is a little fluctuation of the estimated parameters. Finally, the proposed method is compared with other delay compensation strategies in Section 6.

Overview of the Benchmark Problem in RTHS
In RTHS, the emulated structure is divided into PS and NS, and they are related by a transfer system. The equation of motion for the NS is      N N N N N N P P P P g + + + ( , , )=-x M X C X K X r X X X MΓ (1) where M, C, and K are the mass, damping, and stiffness matrices of the structure, respectively; the subscripts "N" and "P" represent the NS and PS, respectively; X denotes the displacement vector of the structure and Γ refers to the location vector; r is the restoring force vector and  g x is the ground acceleration.
Equation (1) is solved by employing a time-stepping integration algorithm, and the calculated displacement XP is forwarded to the PS to satisfy the compatibility condition using a displacement-control-based actuator. Then, the force at the boundary interface is measured directly and fed back into the NS to satisfy the equilibrium. The above procedure is repeated until the termination condition is reached.
However, owing to the dynamics of the control plant, the boundary conditions of NS and PS are often de-synchronized. Hence, in recent years, numerous approaches have been proposed and examined to treat this issue under certain conditions. To provide a unified performance evaluation standard for different methods and to advance the field of RTHS, a Benchmark problem in RTHS was established by Silva et al. [32]. In the Benchmark problem, the models of the servo-hydraulic transfer system, the sensing and control implementation hardware, control-structure interaction, and the PS are identified from the actual equipment [33] and modeled by transfer functions in MATLAB. The block diagram of the benchmark problem is depicted in Figure 1.

Reference Model
As illustrated in Figure 1a, the reference structure is a three-story, two-bay moment-resisting steel frame structure. The structure is made from A572 Grade 50 steel, with columns of commercial S3 × 5.7 sections and beams of built-up I sections that are welded by 50 × 6 webs and 38 × 6 flanges. With the assumptions of (1) negligible vertical translations of the nodes, (2) rigid diaphragm, and (3) a lumped mass at the middle of each span, and applying static condensation to the finite-element model, a three degree-of-freedom (DOF) model is reached to represent the reference structure, that is, only the horizontal DOFs are considered in the numerical model. The stiffness matrix of the reference structure is given by K =

Reference Model
As illustrated in Figure 1a, the reference structure is a three-story, two-bay moment-resisting steel frame structure.
The structure is made from A572 Grade 50 steel, with columns of commercial S3 × 5.7 sections and beams of built-up I sections that are welded by 50 × 6 webs and 38 × 6 flanges. With the assumptions of (1) negligible vertical translations of the nodes, (2) rigid diaphragm, and (3) a lumped mass at the middle of each span, and applying static condensation to the finite-element model, a three degree-of-freedom (DOF) model is reached to represent the reference structure, that is, only the horizontal DOFs are considered in the numerical model. The stiffness matrix of the reference structure is given by K = [2.6055 −2.3134 0.5937; −2.3134 3.2561 −1.4420; 0.5937 −1.4220 0.9267] × 10 7 N/m. The damping matrix of the reference structure is calculated using modal damping ratios.

Partitioning Schemes
It has been reported that the partitioning of the reference structure plays a critical role in the successful implementation of the RTHS [34]. The predictive stability indicator (PSI), which is defined as PSI = log 10 (τ cr ), where τ cr is the critical time delay, is used to quantify the sensitivity of different partitioning choices to de-synchronization at the interface between the substructures. A smaller value of PSI indicates that the RTHS is more sensitive to the de-synchronization. In the Benchmark problem, four different partitioning schemes were considered with PSI values ranging between 0.6 (moderately sensitive) and 1.05 (slightly sensitive), which is equivalent to τ cr values between 4 and 11 ms. The four partitioning cases are listed in Table 1. The PS is a portal frame that is built in the Intelligent Infrastructure Systems Lab (IISL), as described in Figure 1b. The mass, lateral stiffness, and damping coefficient of the PS are 29.12 kg, 1.19 × 10 6 N/m, and 114.6 Ns/m, respectively, which were identified using experimental data.

Control Plant and Control Problem Statement
Figure 1b depicts the model of the control plant, where the dynamics of the servo-valve, actuator, PS, and control-structure interaction (CSI, Dyke et al. [24]) are individually modeled using transfer functions. Specifically, the servo-valve is represented by a second-order transfer function with parameters of α 1 and β 0 to β 2 ; the dynamics of the hydraulic actuator are modeled by a first-order transfer function with parameters of α 3 ; the CSI that caused by the natural velocity feedback loop is modeled by the derivative with parameter α 2 ; the PS is represented by a transfer function with parameters of m P , c P , and k P . The parametric values of the model are identified using experimental data, which are presented in Table 2. Additionally, to validate the robust stability and performance of the designed controller or compensator, model uncertainties were accounted for using a random variable with a normal distribution of each parameter, and are also given in Table 2. Please note that the control plant assigned with the mean values of each parameter are denoted as a nominal plant, whereas that assigned using parameters that differ from the nominal ones are referred to as a perturbed model. As a result of series, parallel, and feedback connections of those components, the mathematical control plant model, from command displacement to the measured displacement, is given by where s is the Laplace transform variable, The magnitude, phase, and time delay of the nominal control plant (Equation (2)) are investigated and presented in Figure 2. In the figure, the time delay is evaluated using τ = ϕ/ω. From Figure 2a,b, it can be observed that when the frequency is low, e.g., the frequency <10 rad/s, both the magnitude and phase are minimal. As the frequency increases, the magnitude decreases while the phase increases, implying that the tracking performance is weakened. Hence, the compatibility and equilibrium cannot be satisfied for all interface boundary DOFs, and the compensation methods for amplitude and phase errors must be employed. Furthermore, the time delay is varying from 25 ms to 8 ms, as shown in Figure 2c. That is, the constant delay-based compensation method may no longer be suitable. Hence, to obtain an accurate RTHS result, the compensation strategy considering the variant time delay and amplitude error must be employed. where s is the Laplace transform variable, A0 = α1β0, B5 = mp, B4 = cp + mpα3 + mpβ1, B3 = (kp + cpα3 + α2) The magnitude, phase, and time delay of the nominal control plant (Equation (2)) are investigated and presented in Figure 2. In the figure, the time delay is evaluated using τ = φ/ω. From Figure 2a,b, it can be observed that when the frequency is low, e.g., the frequency <10 rad/s, both the magnitude and phase are minimal. As the frequency increases, the magnitude decreases while the phase increases, implying that the tracking performance is weakened. Hence, the compatibility and equilibrium cannot be satisfied for all interface boundary DOFs, and the compensation methods for amplitude and phase errors must be employed. Furthermore, the time delay is varying from 25 ms to 8 ms, as shown in Figure 2c. That is, the constant delay-based compensation method may no longer be suitable. Hence, to obtain an accurate RTHS result, the compensation strategy considering the variant time delay and amplitude error must be employed.

Implementation Constraints
As the control plant model is achieved from the physical devices available in the laboratory, virtual RTHS (vRTHS) must be subjected to the following implementation constraints: 1. The compensator has to be in the discrete form.

Implementation Constraints
As the control plant model is achieved from the physical devices available in the laboratory, virtual RTHS (vRTHS) must be subjected to the following implementation constraints:

1.
The compensator has to be in the discrete form.

2.
The sampling frequency of the vRTHS test is 4096 Hz.

3.
The servo-hydraulic actuator responses should be within the maximum capacity of 8900 N, a stroke of 7 mm, and a maximum velocity of 25 mm/s.

4.
The A/D and D/A converters have an 18-bit precision and a span of 3.8 V, and are modeled as a saturation block combined with a quantizer in the Simulink software tool.

5.
Sensor conversion factors from voltage to displacement in mm and voltage to force in N are considered, and are 7.89 mm/V and 1096 N/V, respectively. 6.
The measurement noise is considered and modeled as a Gaussian process whose statistics are obtained from experimental observations. The noise level is approximately 0.05% of the full span of the A/D converter.

Evaluation Criteria
To quantify the performance of the vRTHS, nine evaluation criteria provided in Silva et al. [32] were used in this study. The evaluation criteria J 1 to J 9 are presented in Table 3. J 1 is the tracking time-delay, assessing the time delay between reference and output signals. J 2 is the root mean square error (RMSE), valuating the tracking error using the normalized root-mean square (RMS). J 3 is a normalized peak tracking error measure. J 4 to J 6 are the normalized RMS errors between the reference and RTHS floor displacements and J 7 to J 9 are the normalized peak tracking errors between the reference and RTHS floor displacements. J 4 to J 9 evaluate the performance of the vRTHS compared to the reference model. It should be noted that in this study, the reference solution is obtained by solving equation of motion employing the explicit fourth-order Runge-Kutta scheme with a fixed time step.
Normalized RMS relative displacement at the first floor Normalized RMS relative displacement at the second floor Normalized RMS relative displacement at the third floor Peak relative displacement error corresponding to the first floor Peak relative displacement error corresponding to the second floor Peak relative displacement error corresponding to the third floor It should be noted that d rj (j = 1,2,3) is the reference displacement, and X j (j = 2,3) is the floor displacement of RTHS.

Discrete-Time Adaptive Delay Compensator Based on Kalman Filter
To address the variant time delay in RTHS, numerous control strategies have been proposed and validated over the past 20 years. The feedforward-based approaches have received widespread attention for their fast compensation characteristics in terms of phase and amplitude. By taking the inverse of the control plant model as the feedforward controller, the time delay and amplitude discrepancy can be eliminated thoroughly. However, the mathematic model of the control plant is straightly proper, and the feedforward controller thus often consists of an additional low-pass filter and the inverse model of the control plant [35,36], which may cause design difficulties. Concerning adaptive delay compensation methods, their effectiveness largely depends in the performance of the parameter estimation method. This study proposes a KF-based discrete-time adaptive delay compensation (KF-ADC) method, which is an improvement of the general adaptive compensation method proposed by Wang et al. [13]. As shown in Figure 3, the KF-ADC method incorporates two modules, namely the discrete-time model online estimation module and a feedforward controller. Details of the proposed KF-ADC method are presented in this section.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 7 of 19 attention for their fast compensation characteristics in terms of phase and amplitude. By taking the inverse of the control plant model as the feedforward controller, the time delay and amplitude discrepancy can be eliminated thoroughly. However, the mathematic model of the control plant is straightly proper, and the feedforward controller thus often consists of an additional low-pass filter and the inverse model of the control plant [35,36], which may cause design difficulties. Concerning adaptive delay compensation methods, their effectiveness largely depends in the performance of the parameter estimation method. This study proposes a KF-based discrete-time adaptive delay compensation (KF-ADC) method, which is an improvement of the general adaptive compensation method proposed by Wang et al. [13]. As shown in Figure 3, the KF-ADC method incorporates two modules, namely the discrete-time model online estimation module and a feedforward controller. Details of the proposed KF-ADC method are presented in this section.

Discrete Model of Control Plant
Different from the widely investigated adaptive delay compensation methods that employ the inverse of the continuous model of the control plant [19,28,37], this proposed KF-ADC method employs the discrete-time model of the control plant because of its execution convenience on digital signal processing boards. Whether using the Auto-Regressive with Exogenous discrete-time model [38] or transforming the model from the Laplace domain to the time domain, one can always obtain a difference equation description of the control plant by where d c and d m represent the displacement command and measurement of the control plant, respectively, x is the system model parameter, p and q are two adjustable parameters determined by the orders of the control plant, and the subscript k denotes the k-th time step. Equation (3) can be further reorganized into a more compact form, namely, Please note that before establishing this discrete-time model, the orders of the control plant must be determined beforehand.

Discrete Model of Control Plant
Different from the widely investigated adaptive delay compensation methods that employ the inverse of the continuous model of the control plant [19,28,37], this proposed KF-ADC method employs the discrete-time model of the control plant because of its execution convenience on digital signal processing boards. Whether using the Auto-Regressive with Exogenous discrete-time model [38] or transforming the model from the Laplace domain to the time domain, one can always obtain a difference equation description of the control plant by where d c and d m represent the displacement command and measurement of the control plant, respectively, x is the system model parameter, p and q are two adjustable parameters determined by the orders of the control plant, and the subscript k denotes the k-th time step. Equation (3) can be further reorganized into a more compact form, namely, where Please note that before establishing this discrete-time model, the orders of the control plant must be determined beforehand.

Parameter Estimation
Although the mathematical model of the control plant can be established using the first principle modeling method, it is sometimes unable to suitably represent the actual relationship between the displacement command and measurement. Hence, it is necessary to use the actual testing data to revise the system model parameters.
It can be seen that Equation (4a) is a linear difference equation and Ning et al. [31] and Wang et al. [13] employed the least-square method to estimate the system model parameter. However, the convergence of the least-square method cannot be guaranteed, and it is sometimes time consuming. To overcome the previously mentioned problems, a promising estimator, the KF is used to estimate the system model parameter. In this study, the RTHS system is assumed to be a slow time-varying system, and the state and observed equation for this estimation problem can hence be given by where H k = ϕ k T , y k is the commanded displacement at the k-th time step, and v is the measurement where R is the measurement noise covariance. The KF estimates the system model parameter vector usingx wherex denotes the estimate of x, and G k is the Kalman gain, which is calculated by where P is the estimate error covariance, which is given by

Time-Delay Compensator
When forwarding the RTHS by one step, the model of control plant can be represented by Considering the slow time-varying system assumption in the last subsection, it can be assumed that there are minor changes in the system model parameter between two adjacent time steps. Hence, if the parameter estimation method is convergent, one can substitute x k+1 byx k , namely Thus, the commanded displacement can be predicted by substituting Equation (12) into Equation (11) based on known displacements. By driving the loading system using the predicted displacement, the discrepancy between the desired displacement and the measured displacement is minimized, and the measured displacement will be identical to the desired displacement in an ideal situation, i.e., where d is the desired displacement calculated by a stepwise integration method, such as the central difference method [39]. Therefore, the compensated command displacement can be evaluated by where However, owing to the noise in the measured displacement, the estimated parameters may fluctuate significantly, resulting in the compensated displacement command roughly, which is detrimental to RTHS. Moreover, before time step k+1, the measured displacements are assumed to be close to the desired displacements for the compensated RTHS system. Hence, to enhance the robustness of the compensated RTHS system, the desired displacements are used to predict the displacement command, resulting in the feedforward compensator eventually given by where

Performance of KF Used for Parameter Estimation in RTHS
As described in the last section, only the parameter estimation method is convergent, and the time delay can be well compensated by the proposed KF-ADC. Moreover, it can be seen from Equations (8)-(10) that the performance of the KF is affected by the measurement noise covariance R and the initial values of x 0 and P 0 . Hence, the performance of the KF in RTHS is investigated in this section.

Convergence of KF
It can be seen in the state equation and observed equation of the parameter estimation problem, namely Equations (5) and (6), that the process noise is absent. Hence, the stability of the KF cannot be proven based on the uniform complete controllability and observability of the line system [40]. Let e k denote the difference between the estimated parameter and actual value, namely Substituting Equations (6) and (18) into Equation (8), one can obtain the dynamics of the mathematical expectation of the estimation error, i.e., Here, defining the Lyapunov function for system Equation (19) as follows, then substituting Equations (9), (10), and (19) into Equation (20), the following equation can be obtained.
In RTHS, the elements in the observation matrix H k are measured from the test, and they are varying and linearly independent considering the measurement noise. If H k e k-1 = 0 holds for all values of k, one will obtain e k-1 = 0, namelyx k = x. If H k e k-1 0, then the Lyapunov function is decreased, and the system of Equation (19) is stable. Furthermore, the spectral radius of the matrix I-G k H k is not greater than 1, and is equals to 1 does not hold. Thus, e k-1 →0 as k→∞, namely the estimated parameters by KF can converge to the true value. Hence, the KF is convergent in parameter estimation problems in Section 3.2.

Tracking Performance with Prescribed Displacement Command
For convenience, the control plant is modeled by a second-order transfer function, which is given by where ω and ζ are the equivalent circular frequency and damping ratio, respectively, and the two parameters are set to be 50 rad/s and 0.8, respectively. By using the forward-difference method in Equation (22), a three-parameter discrete control plant model can be obtained, namely and the corresponding feedforward compensator is given by In this section, the sampling frequency is 4096 Hz, resulting in the true parameters for the three-parameter discrete model of the second-order control plant model in Equation (22) being x true = [0.6777 −1.3422 0.6646] T × 10 4 . The initial value of x, which is denoted by x 0 , is determined by employing the least-square method, where the displacement command and measurement originated from the control plant without using any compensation method. Thus, the initial value x 0 is [0.5817 −1.152 0.5704] T × 10 4 . The initial covariance matrix P 0 is assumed to be diagonal, which means that the initial values of the parameters are statistically uncorrelated. The diagonal entries of P 0 are computed as (q f ×x i ) 2 [41], where i = 1, 2, 3, and q f = 0.1. The measurement noise R is set to be 0.001. It should be noted that the measurement noise was modeled by the random number, and a KF was used to suppress the effect of noise.
The tracking performance of the second-order transfer system is shown in Figure 4, where the displacement command is a swept signal with frequency ranging from 0.1 Hz to 20 Hz, and has an amplitude of 100 mm. For comparison, the response without any compensation method is also given in the figure. It is seen that as the input frequency increases, the tracking performance without compensation decreases rapidly, which is detrimental to RTHS. On the contrary, it can be found in Figure 4a that the measurement displacement provided by the proposed KF-ADC method agrees well with the displacement command. The enlarged view is shown in Figure 4b, and it can be seen that the measurement displacement is almost identical to the displacement command. This reveals that by compensating using the proposed method, the control plant exhibits a very good tracking performance.

Tracking Performance with Prescribed Displacement Command
For convenience, the control plant is modeled by a second-order transfer function, which is given by ( ) where ω and ζ are the equivalent circular frequency and damping ratio, respectively, and the two parameters are set to be 50 rad/s and 0.8, respectively. By using the forward-difference method in Equation (22), a three-parameter discrete control plant model can be obtained, namely and the corresponding feedforward compensator is given by In this section, the sampling frequency is 4096 Hz, resulting in the true parameters for the three-parameter discrete model of the second-order control plant model in Equation (22) being xtrue = [0.6777 −1.3422 0.6646] T × 10 4 . The initial value of x, which is denoted by x0, is determined by employing the least-square method, where the displacement command and measurement originated from the control plant without using any compensation method. Thus, the initial value x0 is [0.5817 −1.152 0.5704] T × 10 4 . The initial covariance matrix P0 is assumed to be diagonal, which means that the initial values of the parameters are statistically uncorrelated. The diagonal entries of P0 are computed as (qf×xi) 2 [41], where i = 1, 2, 3, and qf = 0.1. The measurement noise R is set to be 0.001. It should be noted that the measurement noise was modeled by the random number, and a KF was used to suppress the effect of noise.
The tracking performance of the second-order transfer system is shown in Figure 4, where the displacement command is a swept signal with frequency ranging from 0.1 Hz to 20 Hz, and has an amplitude of 100 mm. For comparison, the response without any compensation method is also given in the figure. It is seen that as the input frequency increases, the tracking performance without compensation decreases rapidly, which is detrimental to RTHS. On the contrary, it can be found in Figure 4a that the measurement displacement provided by the proposed KF-ADC method agrees well with the displacement command. The enlarged view is shown in Figure 4b, and it can be seen that the measurement displacement is almost identical to the displacement command. This reveals that by compensating using the proposed method, the control plant exhibits a very good tracking performance. The evolutions of the estimated system model parameter values are shown in Figure 5. It is seen in the figure that the estimated system model parameters vary slightly compared with the initial value after the first 0.1 s. Then, the estimated parameters move rapidly toward the true values The evolutions of the estimated system model parameter values are shown in Figure 5. It is seen in the figure that the estimated system model parameters vary slightly compared with the initial value after the first 0.1 s. Then, the estimated parameters move rapidly toward the true values between 0.1 s and 0.5 s. Subsequently, the estimated parameters remain relatively stable, and are almost identical to the true ones. This result shows that once the dynamics of the control plant are excited, the discrete-time model of the control plant can be quickly estimated using KFs with high-accuracy, and the time-delay can be compensated effectively. This indicates that the KF has a fast convergence speed and high estimation accuracy.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 11 of 19 between 0.1 s and 0.5 s. Subsequently, the estimated parameters remain relatively stable, and are almost identical to the true ones. This result shows that once the dynamics of the control plant are excited, the discrete-time model of the control plant can be quickly estimated using KFs with high-accuracy, and the time-delay can be compensated effectively. This indicates that the KF has a fast convergence speed and high estimation accuracy.

Parameter Study of KF
This subsection investigates the influence of the measurement noise covariance R and the initial values of x0 and P0 on the proposed time-delay compensation method. To perform this analysis, a series of different values of R, x0, and P0 are used. They are considered as Q = WfQ 0 , where Q represents R, x0, and P0, Q 0 represents the nominal values of R, x0, and P0 that were determined in the last subsection, and Wf is a factor varying from 0.7 to 1.3. The duration of the simulations is 10 s, and the stop frequencies of the swept signal are 10 Hz, 20 Hz, and 30 Hz, respectively.

Effect of Initial Guess x0
The time delay of the control plant compensated by the proposed adaptive delay scheme was calculated and presented in Table 4. The table shows that even though the initial guess x0 varied from −30% to 30%, the calculated time delay of the compensated system remained almost unchanged. Moreover, with the increase in the stop frequency of the swept input signal, the calculated time delay varied over a limited range. For all of the 39 simulations, almost the entire calculated time delay was 0.1 ms, of which the maximum value was 0.2 ms. This indicated that the proposed compensation method is effective and robust.

Parameter Study of KF
This subsection investigates the influence of the measurement noise covariance R and the initial values of x 0 and P 0 on the proposed time-delay compensation method. To perform this analysis, a series of different values of R, x 0 , and P 0 are used. They are considered as Q = W f Q 0 , where Q represents R, x 0 , and P 0 , Q 0 represents the nominal values of R, x 0 , and P 0 that were determined in the last subsection, and W f is a factor varying from 0.7 to 1.3. The duration of the simulations is 10 s, and the stop frequencies of the swept signal are 10 Hz, 20 Hz, and 30 Hz, respectively.

Effect of Initial Guess x 0
The time delay of the control plant compensated by the proposed adaptive delay scheme was calculated and presented in Table 4. The table shows that even though the initial guess x 0 varied from −30% to 30%, the calculated time delay of the compensated system remained almost unchanged. Moreover, with the increase in the stop frequency of the swept input signal, the calculated time delay varied over a limited range. For all of the 39 simulations, almost the entire calculated time delay was 0.1 ms, of which the maximum value was 0.2 ms. This indicated that the proposed compensation method is effective and robust. The RMSEs between displacement command and measurement are shown in Figure 6. The figure shows that as the stop frequency increases, the RMSE increases, and the RMSE of the highest stop frequency has the greatest value. Meanwhile, the RMSE decreases first and then increases as W f increases. Despite this, the maximum value of RMSE is less than 5.0%, indicating that the proposed method is still effective with a different initial guess x 0 .  The RMSEs between displacement command and measurement are shown in Figure 6. The figure shows that as the stop frequency increases, the RMSE increases, and the RMSE of the highest stop frequency has the greatest value. Meanwhile, the RMSE decreases first and then increases as Wf increases. Despite this, the maximum value of RMSE is less than 5.0%, indicating that the proposed method is still effective with a different initial guess x0. The convergence values for the discrete model are shown in Figure 7. For comparison convenience, the true values and an error threshold of 1% are also given in the figure. The figures show that the converged values remain almost unchanged with an increase in factor Wf. However, if focusing on one of the specifically prescribed displacements, e.g., a stop frequency of 10 Hz, one can find that the converged values increase with an increase in the factor Wf. Furthermore, the relative error between the converged value and the true value is within 1%, indicating that although there are notable errors in the initial value, the parameters of the control plant can also be significantly improved using KF.  The convergence values for the discrete model are shown in Figure 7. For comparison convenience, the true values and an error threshold of 1% are also given in the figure. The figures show that the converged values remain almost unchanged with an increase in factor W f . However, if focusing on one of the specifically prescribed displacements, e.g., a stop frequency of 10 Hz, one can find that the converged values increase with an increase in the factor W f . Furthermore, the relative error between the converged value and the true value is within 1%, indicating that although there are notable errors in the initial value, the parameters of the control plant can also be significantly improved using KF.  The RMSEs between displacement command and measurement are shown in Figure 6. The figure shows that as the stop frequency increases, the RMSE increases, and the RMSE of the highest stop frequency has the greatest value. Meanwhile, the RMSE decreases first and then increases as Wf increases. Despite this, the maximum value of RMSE is less than 5.0%, indicating that the proposed method is still effective with a different initial guess x0. The convergence values for the discrete model are shown in Figure 7. For comparison convenience, the true values and an error threshold of 1% are also given in the figure. The figures show that the converged values remain almost unchanged with an increase in factor Wf. However, if focusing on one of the specifically prescribed displacements, e.g., a stop frequency of 10 Hz, one can find that the converged values increase with an increase in the factor Wf. Furthermore, the relative error between the converged value and the true value is within 1%, indicating that although there are notable errors in the initial value, the parameters of the control plant can also be significantly improved using KF.

Effect of Initial Covariance Matrix P 0 and R
The relationships between the RMSE of the compensated system and the initial covariance and measurement noise covariance are shown in Figure 8. The figures show that for a specific displacement input, as the value of W f increases, RMSE remains almost the same. Moreover, with an increase in the stop frequency, the RMSE increases. However, the maximum value of the RMSE is less than 5%, indicating that the proposed method is still effective with varying initial guess values.

Effect of Initial Covariance Matrix P0 and R
The relationships between the RMSE of the compensated system and the initial covariance and measurement noise covariance are shown in Figure 8. The figures show that for a specific displacement input, as the value of Wf increases, RMSE remains almost the same. Moreover, with an increase in the stop frequency, the RMSE increases. However, the maximum value of the RMSE is less than 5%, indicating that the proposed method is still effective with varying initial guess values. The relationships between the converged value and the factor Wf are shown in Figures 9 and 10. Similar to that used in Figure 7, the true values and the error threshold of 1% are also given in the figure. It can be seen that for all of the considered cases, regardless of which of the three variables changes, i.e., R, P0, or Wf, the converged value remains at almost the same level. The maximum relative error between the convergency value and the true value is within 1%, indicating that the parameters of the control plant are improved notably, and the parameter estimation method exhibits strong robustness.

Numerical Validation
Virtual RTHSs were carried out with three different ground acceleration records, and the four partitioning cases shown in Table 1 were considered for each input. Moreover, one nominal vRTHS and 20 perturbed vRTHSs were included to consider the modeling uncertainties for the control plant in each partitioning configuration. A total of 252 vRTHSs (12 nominal vRTHSs and 240 perturbed vRTHSs) were carried out in this study. The relationships between the converged value and the factor W f are shown in Figures 9 and 10. Similar to that used in Figure 7, the true values and the error threshold of 1% are also given in the figure. It can be seen that for all of the considered cases, regardless of which of the three variables changes, i.e., R, P 0 , or W f , the converged value remains at almost the same level. The maximum relative error between the convergency value and the true value is within 1%, indicating that the parameters of the control plant are improved notably, and the parameter estimation method exhibits strong robustness.

Effect of Initial Covariance Matrix P0 and R
The relationships between the RMSE of the compensated system and the initial covariance and measurement noise covariance are shown in Figure 8. The figures show that for a specific displacement input, as the value of Wf increases, RMSE remains almost the same. Moreover, with an increase in the stop frequency, the RMSE increases. However, the maximum value of the RMSE is less than 5%, indicating that the proposed method is still effective with varying initial guess values. The relationships between the converged value and the factor Wf are shown in Figures 9 and 10. Similar to that used in Figure 7, the true values and the error threshold of 1% are also given in the figure. It can be seen that for all of the considered cases, regardless of which of the three variables changes, i.e., R, P0, or Wf, the converged value remains at almost the same level. The maximum relative error between the convergency value and the true value is within 1%, indicating that the parameters of the control plant are improved notably, and the parameter estimation method exhibits strong robustness.

Numerical Validation
Virtual RTHSs were carried out with three different ground acceleration records, and the four partitioning cases shown in Table 1 were considered for each input. Moreover, one nominal vRTHS and 20 perturbed vRTHSs were included to consider the modeling uncertainties for the control plant in each partitioning configuration. A total of 252 vRTHSs (12 nominal vRTHSs and 240 perturbed vRTHSs) were carried out in this study.

Effect of Initial Covariance Matrix P0 and R
The relationships between the RMSE of the compensated system and the initial covariance and measurement noise covariance are shown in Figure 8. The figures show that for a specific displacement input, as the value of Wf increases, RMSE remains almost the same. Moreover, with an increase in the stop frequency, the RMSE increases. However, the maximum value of the RMSE is less than 5%, indicating that the proposed method is still effective with varying initial guess values. The relationships between the converged value and the factor Wf are shown in Figures 9 and 10. Similar to that used in Figure 7, the true values and the error threshold of 1% are also given in the figure. It can be seen that for all of the considered cases, regardless of which of the three variables changes, i.e., R, P0, or Wf, the converged value remains at almost the same level. The maximum relative error between the convergency value and the true value is within 1%, indicating that the parameters of the control plant are improved notably, and the parameter estimation method exhibits strong robustness.

Numerical Validation
Virtual RTHSs were carried out with three different ground acceleration records, and the four partitioning cases shown in Table 1 were considered for each input. Moreover, one nominal vRTHS and 20 perturbed vRTHSs were included to consider the modeling uncertainties for the control plant in each partitioning configuration. A total of 252 vRTHSs (12 nominal vRTHSs and 240 perturbed vRTHSs) were carried out in this study.

Numerical Validation
Virtual RTHSs were carried out with three different ground acceleration records, and the four partitioning cases shown in Table 1 were considered for each input. Moreover, one nominal vRTHS and 20 perturbed vRTHSs were included to consider the modeling uncertainties for the control plant in each partitioning configuration. A total of 252 vRTHSs (12 nominal vRTHSs and 240 perturbed vRTHSs) were carried out in this study.
In this section, a four-parameter discrete model of the control plant is adopted, namely and the corresponding feedforward controller is To determine the initial values of x 0 , P 0 , and R, a virtual RTHS (vRTHS) without any compensation strategy was carried out, where the El Centro ground acceleration is scaled to 0.2 g. By employing the least-square method, the initial value x 0 reaches [11.1688 −20.4226 13.6773 −3.4339] T . The measurement noise variance R is set to 1 × 10 -4 , which is determined by the characteristics of the noise level. Similar to that used in the last subsection, the initial covariance matrix is given by P 0 = diag([11.1688 2 20.4226 2 13.6773 2 3.4339 2 ]).
The displacement time history of the first floor from nominal virtual RTHS for Case 1 under the El Centro earthquake record is shown in Figure 11, with the reference response provided. From the figure, it can be seen that the displacement of the virtual RTHS is almost identical to the reference one, which indicates a good tracking performance of the proposed compensation strategy. The evolutions of the estimated system model parameters are given in Figure 12, correspondingly. It can be seen that once the dynamics of the control plant are excited and there are differences between the desired and measurement displacements, the KF responds rapidly to adjust the system model parameter, and then the parameters converge rapidly to a relatively steady value. It may be noted that the convergence trend of x 3 and x 4 is not very evident, but the compensated control plant exhibits a nearly perfect tracking performance, as shown in Figure 11b. This is because the parameters x 3 and x 4 are related to the dynamics of the control plant in the high-frequency range, which are not excited completely. Furthermore, it is seen that there are few large fluctuations during the loading period, implying the stable performance of the parameter estimation method.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 14 of 19 In this section, a four-parameter discrete model of the control plant is adopted, namely (25) and the corresponding feedforward controller is To determine the initial values of x0, P0, and R, a virtual RTHS (vRTHS) without any compensation strategy was carried out, where the El Centro ground acceleration is scaled to 0.2 g. By employing the least-square method, the initial value x0 reaches [11.1688 −20.4226 13.6773 −3.4339] T . The measurement noise variance R is set to 1 × 10 -4 , which is determined by the characteristics of the noise level. Similar to that used in the last subsection, the initial covariance matrix is given by P0 = diag([11.1688 2 20.4226 2 13.6773 2 3.4339 2 ]).
The displacement time history of the first floor from nominal virtual RTHS for Case 1 under the El Centro earthquake record is shown in Figure 11, with the reference response provided. From the figure, it can be seen that the displacement of the virtual RTHS is almost identical to the reference one, which indicates a good tracking performance of the proposed compensation strategy. The evolutions of the estimated system model parameters are given in Figure 12, correspondingly. It can be seen that once the dynamics of the control plant are excited and there are differences between the desired and measurement displacements, the KF responds rapidly to adjust the system model parameter, and then the parameters converge rapidly to a relatively steady value. It may be noted that the convergence trend of x3 and x4 is not very evident, but the compensated control plant exhibits a nearly perfect tracking performance, as shown in Figure 11b. This is because the parameters x3 and x4 are related to the dynamics of the control plant in the high-frequency range, which are not excited completely. Furthermore, it is seen that there are few large fluctuations during the loading period, implying the stable performance of the parameter estimation method.   The quantitative indicators that were introduced earlier for the nominal vRTHSs under El Centro, Kobe, and Morgan Hill ground motions are presented in Table 5. In the table, it is shown that the maximum value of J 1 , which is the calculated time delay, is 0.2 ms, which is much less than the critical time delay (approximately 4 ms) in all cases. This indicates that the proposed KF-ADC method is efficient in terms of providing a stable RTHS considering different partitioning schemes. Moreover, the maximum value of J 2 and J 3 is less than 2%, implying that the KF-ADC method also achieves good global agreement between the desired and measured displacements in all cases.    If focusing on the criteria between the results from vRTHS and the reference solution, namely the global errors (J 4 , J 5 , and J 6 ) and relative peak errors (J 7 , J 8 , and J 9 ), it can be seen that most of them are less than 5%. This demonstrates that the proposed KF-ADC method enables the RTHS to reliably assess the structural performance under different dynamic loads. Furthermore, it can also be observed that with the increase in the sensitivity of different partitioning schemes, both the global errors and relative peak errors generally tend to increase. However, the maximum value of the criteria J 4 to J 9 is still acceptable. This indicates that the proposed KF-ADC strategy enables RTHS to obtain reliable results when the partitioning configuration is in the moderately and slightly sensitive regions.
To further investigate the robustness of the KF-ADC strategy proposed in this study, the mean and standard deviation of the evaluation criteria for the 240 perturbed vRTHSs are shown in Tables 6 and 7, respectively. By comparing the mean of the nine criteria from the perturbed vRTHSs in Table 6 and those from the nominal vRTHSs in Table 5, similar trends can be observed. The maximum value of the criteria for J 4 to J 9 is less than 15%. Moreover, based on the standard deviation presented in Table 7, the less sensitive partitioning scheme is preferred in order to obtain more robust RTHS results. However, the standard deviation of the evaluation criteria of each partitioning configuration and load case is generally less than 2% and no more than 3%, which demonstrates that the proposed KF-ADC method exhibits good robustness when uncertainties are considered in the control plant.

Comparison with Existing Methods
To comprehensively investigate the KF-ADC strategy proposed in this study, several other promising delay compensation methods were employed to provide the reference for evaluating the proposed compensator strategy. These are the widely used proportional integral (PI) control with a phase lead compensator (PI, Silva et al. [32]), the model-based controller (MBC, Phillips and Spencer Jr [19], the adaptive model reference controller (aMRC, Najafi and Spencer [42]), and the tracking error-based adaptive delay compensation method (TE-AC, Tao and Mercan [43]). The simulation results provided in Najafi and Spencer [42] and Tao and Mercan [43] are cited here.
The evaluation criteria for vRTHS employing nominal parameters with different partitioning configurations under the El Centro earthquake are shown in Figure 13. It can be seen in the figure that the nine criteria provided by PI are generally the largest, indicating that the traditional control method may not be very effective to deal with variant time delay in RTHS. With the exception of the tracking criteria of J 2 and J 3 , the simulation results generated by aMRC are better than that by MBC, implying that the time delay in RTHS is varying, and constant delay compensation is not sufficient. The aMRC exhibits better performance compared to TE-AC, which is mainly caused by the large number of parameters that should be adjusted in TE-AC.
Furthermore, it can be seen in Figure 13 that the proposed delay compensation scheme is competitive compared with the promising aMRC strategy. With the exception of the calculated time delay of the proposed method is a little larger than that of the aMFC, the other eight evaluation criteria provided by the proposed method are generally smaller, especially for J 2 and J 3 . The reason is that the low-pass filter and the reference model in the aMRC strategy should be determined empirically, which are extremely important. In summary, the proposed method is effective and practical for use in RTHS. results provided in Najafi and Spencer [42] and Tao and Mercan [43] are cited here.
The evaluation criteria for vRTHS employing nominal parameters with different partitioning configurations under the El Centro earthquake are shown in Figure 13. It can be seen in the figure that the nine criteria provided by PI are generally the largest, indicating that the traditional control method may not be very effective to deal with variant time delay in RTHS. With the exception of the tracking criteria of J2 and J3, the simulation results generated by aMRC are better than that by MBC, implying that the time delay in RTHS is varying, and constant delay compensation is not sufficient. The aMRC exhibits better performance compared to TE-AC, which is mainly caused by the large number of parameters that should be adjusted in TE-AC.

Conclusions
To deal with the variant time delay caused by the dynamics of control plant in RTHS, a KF-based discrete-time adaptive delay compensation (KF-ADC) method was proposed. Numerical simulations on the Benchmark problem in RTHS were performed to examine the performance of the proposed method. The main conclusions are as follows.

1.
The principle of the proposed delay compensation method for RTHS was provided. The control plant is represented by a discrete model with time-varying coefficients, which are estimated and updated by KF using the displacement commands and measurements. Based on this model, the compensator is constructed by employing the desired displacements to deal with the variant time delay. 2.
The performance of the KF used for parameter estimation in RTHS was investigated. By employing the Lyapunov stability theory, it was confirmed that the KF can converge to the true value. Simulations revealed that the KF converged rapidly. When the initial settings of the KF, namely R, x 0 , and P 0 , varied with 30% of the nominal value, the relative errors of the converged values of the system parameters were less than 1%. It means errors in the initial settings of the KF had a limited effect on the convergent value, and the KF was robust to initial settings.

3.
The effectiveness of the proposed method was validated through the Benchmark problem in RTHS. The results of 252 vRTHSs revealed a satisfactory tracking performance and robustness of the proposed strategy. Specifically, the maximum value of mean of the estimated time delay was less than 0.5 ms, the maximum value of tracking errors (J 2 and J 3 ) was less than 2%, and the global errors (J 4 , J 5 , and J 6 ) and relative peak errors (J 7 , J 8 , and J 9 ) were less than 5%. Moreover, the proposed method was compared with the other four promising schemes. Results revealed that adaptive delay compensations were better than the constant delay-based compensation strategy. The proposed method was competitive compared with the adaptive model reference controller. In general, the proposed method exhibited a superior performance.
Despite the performance of the proposed strategy can be well illustrated by the virtual RTHSs, where the measurement noise, uncertainties of loading system and specimen, real earthquake records and multiple degrees-of-freedom structures are considered, physical tests are also needed to further validate the actual performance. These actual tests are already planned.