High performance compensation using an adaptive strategy for real-time hybrid simulation

High performance compensation using an adaptive strategy for real-time hybrid simulation

Mechanical Systems and Signal Processing 133 (2019) 106262 Contents lists available at ScienceDirect Mechanical Systems and Signal Processing journa...

2MB Sizes 2 Downloads 15 Views

Mechanical Systems and Signal Processing 133 (2019) 106262

Contents lists available at ScienceDirect

Mechanical Systems and Signal Processing journal homepage: www.elsevier.com/locate/ymssp

High performance compensation using an adaptive strategy for real-time hybrid simulation Zhen Wang a,b,c,1, Xizhan Ning d,1, Guoshan Xu a,b,c,⇑, Huimeng Zhou a,e, Bin Wu f a

Key Lab of Structures Dynamic Behavior and Control, Ministry of Education, Harbin Institute of Technology, Harbin 150090, China Key Lab of Intelligent Disaster Mitigation, Ministry of Industry and Information Technology, Harbin 150090, China c School of Civil Engineering, Harbin Institute of Technology, Harbin 150090, China d College of Civil Engineering, Huaqiao University, Xiamen 361021, China e Key Laboratory of Earthquake Engineering and Engineering Vibration, Institute of Engineering Mechanics, China Earthquake Administration, Harbin 150080, China f School of Civil Engineering and Architecture, Wuhan University of Technology, Wuhan 430070, China b

a r t i c l e

i n f o

Article history: Received 17 February 2019 Received in revised form 14 June 2019 Accepted 20 July 2019 Available online 8 August 2019 Keywords: Real-time hybrid simulation Adaptive filter Adaptive delay compensation Robustness

a b s t r a c t Real-time hybrid simulation (RTHS) is an innovative and versatile technique for evaluating the dynamic responses of structural and mechanical systems. This technique separates the emulated system into numerical and physical substructures, which are analyzed by computers and loaded in laboratories, respectively. Ensuring the boundary conditions between the two substructures through a transfer system plays a significant role in obtaining reliable and accurate testing results. However, measurement noise and the delay between commands and responses due to the dynamic performance of the transfer system are inevitable in RTHS. To address these issues and to achieve outstanding tracking performance and excellent robustness, this paper proposes an adaptive Kalman-based noise filter and an adaptive two-stage delay compensation method. In particular, in the novel noise filter strategy, adaptive inverse compensation with parameters updated by the least squares method is adopted to accommodate the amplitude and phase errors induced by a traditional Kalman filter. In the proposed delay compensation method, classic polynomial extrapolation and an adaptive inverse strategy are employed for coarse and fine compensation, respectively. Virtual RTHS on a benchmark problem reveals the satisfactory tracking performance and robustness of the proposed methods. Comparisons with polynomial extrapolation and single-stage adaptive compensation indicate the superiority of the proposed two-stage delay compensation method. Ó 2019 Elsevier Ltd. All rights reserved.

1. Introduction Real-time hybrid simulation (RTHS) is an innovative and versatile technique for evaluating the dynamic responses of structural and mechanical systems subjected to earthquakes and other dynamic loads [1–5]. This technique, according to simulation complexities and research objectives, splits the emulated systems into a numerical substructures (NS) and a physical substructures (PS) [6]. The NS often denotes portions numerically modeled and simulated, while the PS represents portions physically fabricated and investigated in laboratories that are not suitable for conducting numerical analysis due to ⇑ Corresponding author at: School of Civil Engineering, Harbin Institute of Technology, Harbin 150090, China. 1

E-mail address: [email protected] (G. Xu). The first two authors, Zhen Wang and Xizhan Ning, contributed equally to this research work.

https://doi.org/10.1016/j.ymssp.2019.106262 0888-3270/Ó 2019 Elsevier Ltd. All rights reserved.

2

Z. Wang et al. / Mechanical Systems and Signal Processing 133 (2019) 106262

complexities and rate-dependent behavior. This ‘divide and conquer’ strategy takes advantage of numerical simulations and physical testing and endows versatile performance on the RTHS in terms of a shorter testing period, lower cost of specimens and setups, and accurate testing results. Since being proposed in 1992 [7], this technique has drawn increasing attention and has been extensively investigated worldwide. RTHS has been applied to various engineering projects regarding offshore platforms, wood buildings and steel frame structures [8–11]. The critical issue of RTHS is accurately replicating the boundary conditions between the PS and NS via a servo hydraulic loading system, and this issue is an indispensable prerequisite of reliable simulation results. However, the dynamic performance of the loading system inevitably introduces a time delay between the commanded and achieved displacements. Investigations [12–14] revealed that this delay has a negative impact on the accuracy of testing results and even induces instability problems for RTHS of one system with a spring specimen. Various methods have been proposed to address this issue; among others, see [12,15–25]. The most extensively applied method for delay compensation is the method based on polynomial extrapolation (PE) proposed by Horiuchi et al. [12] and improved by others [15–17]. This method assumes a constant delay in the loading system, and this delay is compensated for by imposing predicted structural responses after the delay time to the PS. To predict the response, a polynomial is fitted based on several past data points. However, the loading system delay often varies due to signal frequencies and specimen behavior; thus, compensation methods accommodating this variation, e.g., adaptive methods, are desirable. As an alternative, Ning et al. [18] mitigated this variation by H1 control and obtained excellent tracking and robustness performance using PE. Various adaptive strategies for compensating variable delay have been proposed. One category of these methods is based on online delay estimation, where estimated delay is applied to PE for delay compensation [15,19,20]. Another category adopts an inverse control strategy, where the inverse model for the loading system is identified online [21–23]. Additionally, adaptive compensation is conceived based on the synchronization error and nonlinear estimator in [14,24] and [26]. Research shows that adaptive compensation methods perform better than traditional PE. In addition to delay compensation, robustness control [18,27,28] and nonlinear control [29,30] have also drawn considerable attention in recent years. These methods aim to obtain reliable test results through sophisticated control schemes, even though the loading system is not accurately modeled and the system exhibits nonlinear behavior. For further investigation and development of RTHS, a benchmark problem [31] has recently been designed. Participants are encouraged to provide new types of solutions and present corresponding evaluation criteria relevant to tracking performance and robustness for this problem. Aiming to achieve high performance control in terms of outstanding tracking performance and robustness for this benchmark problem, this paper conceives a strategy with the assistance of a Kalman filter [32], least squares method and adaptive compensation. In particular, this paper proposes an adaptive filter to treat noise and adaptive two-stage compensation for delay. An investigation shows the comparatively satisfactory performance. The remainder of this paper is organized as follows. Section 2 presents a brief overview of this benchmark problem for RTHS. The proposed adaptive filter is described in Section 3, while Section 4 explains the proposed adaptive two-stage compensation. Virtual RTHS is performed and presented in Section 5. Section 6 provides a comparison between two-stage compensation with conventional PE and single-stage adaptive delay compensation (SADC). Section 7 concludes this research. 2. Overview of the benchmark problem for RTHS To provide a realistic problem of actuator tracking control in RTHS for exploring various solution methods and then further advancing this field, a flexible framework based on an actual RTHS experiment is designed and published [31]. This benchmark problem for RTHS contains a reference model, four partitioning cases, implementation constraints, a transfer system model with parametric uncertainties, and relevant MATLAB codes and Simulink models. Participants are expected to design a controller for this problem by considering not only accurate tracking performance but also strong robustness. Accordingly, this paper proposes a control strategy endowed with satisfactory tracking performance and robustness. In this benchmark problem, the emulated structure is an actual 3-story, 2-bay steel frame. By applying 4 assumptions [31], the dynamic model of this structure is simplified into 3 degrees of freedom in the horizontal direction. The leftbottom story is chosen as the PS, and the remainder is the NS, as shown in Fig. 1. To objectively evaluate the control schemes of participants, four partitioning cases are designed with various mass and modal damping of each floor of the reference structure, and these cases lead to distinct sensitivity of the testing system to desynchronization between the desired and actual displacement of the PS. Among the four cases, the fourth is the most sensitive to this loading error. The block diagram of this benchmark problem is depicted in Fig. 2. As shown in the top figure, the measured action force and earthquake record are fed to the NS, and then the desired displacement y at the interface of the two substructures is resolved and sent to the highlighted tracking controller, which is expected to design for ensuring the synchronization of the substructures. The commands yc from this controller are sent to the control plant, i.e., the transfer system and the PS. As shown in the bottom figure, the control plant is composed of a servo valve, an actuator, the PS and control-structure interactions. Moreover, to consider parametric uncertainty, actual variations of various parameters of the control plant are identified and published. Clearly, this benchmark problem provides a systematic and detailed model of an actual testing system that can serve as a platform for evaluating distinct control schemes. Additionally, this control problem is flexible to be modified and applied to more complicated RTHS. Consequently, this paper is devoted to designing a control strategy together with an adaptive filter for this problem.

Z. Wang et al. / Mechanical Systems and Signal Processing 133 (2019) 106262

3

Fig. 1. Reference model.

3. Adaptive filter for measurement noise rejection Measurement noise is always present in signals due to A/D converters and ambient electromagnetic interference. This noise affects the signal-to-noise ratio of the measured displacement. Moreover, delay compensation, e.g., PE, often amplifies signal components of higher frequency. In this sense, noise may induce risk to the loading system. Therefore, to ensure successful RTHS and reliable test results, noise must be effectively suppressed. Common filters, such as finite impulse response, cause phase lag to the signal and then increase the loading system delay and control complexity. As an alternative, a Kalman filter [32] has been employed in RTHS due to its promising performance [33]. However, as a model-based scheme, this filter’s performance depends on the model accuracy, and parameter variations of the loading system degrade the performance. For high performance control for RTHS, this paper proposes an adaptive filter consisting of a Kalman filter and adaptive error compensation, as shown in Fig. 3. In this figure, yn ; ym and yf represent the noise-contaminated displacement, recovered displacement, and Kalman filtered displacement, respectively. Additionally, ycl is the clean displacement, free of noise.

3.1. Noise rejection using a Kalman filter In RTHS, the control accuracy is compromised by considering the noise in the feedback measurement, especially when the signal-to-noise ratio is low. To obtain excellent tracking performance, the effect of noise must be suppressed. Hence, a state

Fig. 2. Block diagram of the benchmark problem.

4

Z. Wang et al. / Mechanical Systems and Signal Processing 133 (2019) 106262

Fig. 3. Feedback control based on the proposed adaptive filter for noise.

estimator, known as a Kalman filter [32], is adopted for RTHS in this study. The control plant in state-space form can be given by

x_ P ¼ AP xP þ BP ðuP þ eÞ yP ¼ C P xP þ DP uP þ v

ð1Þ

with

E½e ¼ 0; E½v ¼ 0;       E eeT ¼ Q; E eeT ¼ R; E evT ¼ 0 where the subscript p refers to the control plant, e and v are process and measurement noises, respectively, E represents expectation, and Q and R are the covariance of e and v, respectively. The steady-state Kalman gain Kf is given by

K f ¼ P e C p R1

ð2Þ

where the subscript f refers to the Kalman filter and Pe is the steady-state covariance. Hence, the steady-state Kalman filter can be recast into

^_ P ¼ AF x ^P þ BF uF y ^P þ DF uF ^P ¼ C F x x

ð3Þ

^P is the estimated state of the control plant; y ^P is the predicted measurement; uF is the input of the Kalman filter, where x consisting of the controller output and the measurement displacement, namely, uF = [uP yP]T; and the matrices AF, BF, CF and DF are the state-space realization of the Kalman filter given by

AF ¼ AP  K F C P ; BF ¼ ½ BP  K F DP

CF ¼ CP K F ; DF ¼ ½ DP

0

ð4Þ

In real-world applications, nominal parameters of loading systems are not often available; hence, this study regards the loading system as a black box. Consequently, a new second-order transfer function model is identified for this loading system using a swept sine signal with a frequency varying from 0.1 Hz to 10 Hz. In the simulation, the control plant runs with its nominal parameters. The resulting model is written as



9589 s2 þ 238:4s þ 9518

ð5Þ

Subsequently, the Kalman filter is designed using Q/R = 0.01 and R = 0.1. Investigations show that measurement noise is effectively suppressed by this Kalman filter; however, minor displacement errors are also introduced in terms of amplitude and phase discrepancies. These errors, which are acceptable for ordinary applications, have to be treated herein since this study aims to achieve high performance control. Additionally, these errors result in desynchronization for filtered displacement and specimen force because the Kalman filter is employed only for the displacement signal. To further improve the control performance, this paper devises an adaptive compensation for this Kalman filter error, as shown in Fig. 3 and described in the next subsection.

5

Z. Wang et al. / Mechanical Systems and Signal Processing 133 (2019) 106262

3.2. Adaptive compensation for Kalman filter error This subsection proposes filter compensation based on the least squares method to effectively mitigate the aforementioned filter errors. Similar schemes have been widely adopted for delay estimation and compensation in RTHS [20,22]. The relation between Kalman filtered displacement and noisy measurement displacement is simplified to a pure delay together with an amplitude factor, namely,

yn ðtÞ ¼ kf yf ðt þ sf Þ

ð6Þ

where y, k, t and s denote the displacement, amplitude factor, time and delay, respectively, and the subscripts n and f represent noisy and filtered, respectively. The filtered displacement is replaced with its first-order approximation obtained using Taylor expansion, yielding

yn ðtÞ ¼ kf yf ðtÞ þ kf y_ f ðtÞsf

ð7Þ

Although it is possible to directly obtain the filtered velocity from the Kalman filter, this velocity is evaluated by means of finite difference, namely,

y_ f ðtÞ ¼

yf ðtÞ  yf ðt - DtÞ Dt

ð8Þ

where Dt is the time interval of data points. Substituting Eq. (8) into Eq. (7) yields

yn ðtÞ ¼ kf ð1 þ

sf Dt

Þyf ðtÞ  kf

sf Dt

yf ðt - DtÞ

ð9Þ

This result indicates a simplified relation between the noisy and filtered displacements. For simplicity, this relation can be recast into a discrete-time equation as

yn;i ¼ hn1 yf;i þ hn2 yf;i1

ð10Þ

in which hn1 ; hn2 are two coefficients and i denotes the i-th sampling. With a group of filtered and measured displacements at different time instants, the parameters in Eq. (10) can be identified as a linear regression problem. With the assistance of the estimated parameters, the uncontaminated signal can be recovered by means of

ym;i ¼ ^hn1 yf;i þ ^hn2 yf;i1

ð11Þ

in which ^ hn1 ; ^ hn2 are estimated parameters. Therefore, ym is less contaminated by noise since yf is almost free of noise, and the model parameters ^ hn1 ; ^ hn2 are less influenced by noise in y due to the upcoming identification method. Therefore, y is a n

m

reasonable representation of ycl . The main procedure of this adaptive filter is depicted in Fig. 3. Clearly, the identification method plays a significant role because of its influence on the quality of the recovered displacement. In view of its characteristics, the least squares algorithm is adopted to estimate these parameters. Given possible parameter variation, the recursive least squares algorithm with a forgetting factor [34] is a more favorable candidate, and this algorithm is performed with weighted data and requires an initial guess to start up the process. This method is suitable for online estimation of time-varying parameters because of its small storage size and low calculation efforts as well as its weighted data. The recursive formulas of the method are _

_

h i ¼ h i1 þ

P i1 Wi WTi P i1 Wi þ

"

h

q

P i1 Wi WTi P i1 P i1  T Pi ¼ q Wi P i1 Wi þ q 1

 Wi ¼ yf;i

yf;i1

T

_

ym;i ðt i Þ  WTi h i1

i

ð12Þ

# ð13Þ

ð14Þ

where q is the forgetting factor, 0.9 < q  1. The greater the forgetting factor q is, the greater the effect that the previous data have on the current estimated delay. When q = 1, the algorithm degenerates to the recursive least squares algorithm. Consideration regarding the time interval Dt and initial parameters is discussed here. By comparing Eqs. (10) and (9), one can state that a smaller Dt leads to larger fluctuations in the estimated parameters under the same other conditions. Moreover, a smaller Dt renders the estimation sensitive to measurement noise and A/D conversion error. To obtain more stable estimation results, the time interval has to be set larger. This paper chooses Dt to be 40/4096 s. On the other hand, theoret^ can be derived ically, the amplitude gain and delay are close to 1 and 0, respectively; accordingly, the initial parameter for h

as [1 0]T. In this research, P is initialized by an identity matrix multiplied by 1000. To demonstrate the performance of the adaptive filter, Fig. 4 depicts the results of the virtual RTHS of partitioning Case 1 subjected to the El Centro earthquake compensated by the proposed method in the next section. Nominal plant parameters are adopted here. In Fig. 4a), the dotted, dashed and thick solid lines mean the relative errors between the clean and noisy

6

Z. Wang et al. / Mechanical Systems and Signal Processing 133 (2019) 106262

Fig. 4. Performance of the adaptive filter for suppressing measurement noise.

displacements, between the clean and Kalman filtered displacements, and between the clean and filtered compensated displacements, respectively. The three lines represent the measurement noise, the performance of the Kalman filter and the proposed adaptive filter, respectively. In the first 5 s, the Kalman filter effectively suppresses the noise to a very low level, in that the commands of the loading system are zero. After that, the Kalman filter performs well except for some exceptions. Conversely, the adaptive filter exhibits consistent and wonderful performance. Fig. 4b) shows that the estimated parameters vary stably. 4. Adaptive two-stage delay compensation for RTHS PE is suitable for constant delay; in most scenarios, however, system delay varies for different reasons. Adaptive delay compensation can accommodate delay change; however, its performance always depends on parameter identification accuracy. Especially for a case with large estimated delay variation, the identification requires more time before convergence, which actually affects compensation efficacy. Additionally, a delay often means a strong nonlinear relation between the command and response of loading systems, and a larger delay implies a larger truncation error when the delay is approximated. Therefore, for high performance control, this study combines these two schemes to form a new strategy, i.e., adaptive twostage delay compensation, as shown in Fig. 5. Specifically, the first stage means compensation for the main parts of the system delay by means of PE, referred to as course compensation; the second stage represents compensation with adaptive inverse compensation for the remaining delay, called fine compensation. This new scheme is presented in this section, followed by validation in the subsequent sections. 4.1. Stage one: delay compensation by polynomial extrapolation PE is extensively adopted for delay compensation in RTHS due to its simplicity and effectiveness, which encourages us to implement this method for this benchmark problem. Three parameters have to be determined for this method before its application, namely, the polynomial order, the time interval of displacement sampling and the system delay. Generally, the higher the order is, the worse the stability of the testing is considering the time integration. For this benchmark problem, the time interval of the integration is so small that there is no concern regarding instability due to the integration algorithm [33]. Additionally, higher order extrapolation often means less prediction error. Therefore, third-order extrapolation is employed herein, namely,

    yc;i ¼ 1 þ 11 g þ g2 þ 16 g3 yac;i  3g þ 32 g2 þ 12 g3 yac;i1 3 6    þ 2 g þ 2g2 þ 12 g3 yac;i2  13 g þ 12 g2 þ 16 g3 yac;i3

ð15Þ

where g ¼ Dst; yac is the output displacement of Stage 2 compensation, shown in Fig. 5; s is the time delay of the closedloop system; and Dt is the time interval. Regarding the interval, the case with too small intervals implies that data of a very short time period are applied to fit a polynomial, which is then used to predict the displacement command after a relatively large time period, namely, the delay. Inevitably, these data have been influenced by force measurement noise and other uncertainties. Due to the smaller time interval, the actual displacement increments are small, and the uncertainties can result in considerable fluctuations around

7

Z. Wang et al. / Mechanical Systems and Signal Processing 133 (2019) 106262

Fig. 5. Feedback control with the proposed adaptive two-stage delay compensation.

the displacement peaks. Conversely, an excessively large time interval indicates the strong impact of very old data on the prediction. Therefore, a relatively reasonable method is to set the interval as the delay level. In this study, Dt is set to the sampling period multiplied by 40. 4.2. Stage 2: adaptive fine compensation for loading errors Extrapolation for delay compensation performs well together with a sophisticated controller for the loading system; among others, see [18]. However, the delay of a system with a classic controller, e.g., PI in this benchmark example, varies with respect to the input frequency. In addition, system parameter uncertainty has to be taken into account and also induces delay variation. In these scenarios, the performance of polynomial delay compensation is somewhat unsatisfactory, as presented in Section 6. To compensate for this error, this paper proposes to attach an adaptive strategy to the classic PE. This scheme outperforms SADC, as shown in forthcoming simulations. To save space, we directly suppose a difference equation with three parameters as the model, written as

yac;i ¼ ha1 ym;i þ ha2 ym;i1 þ ha3 ym;i2

ð16Þ

where y represents the displacement and subscripts m and a denote the measurement and adaptive compensator, respectively. One can derive this equation by supposing that the system is a second-order system and substituting the differential terms with a difference in displacement with respect to time. Since our objective is to compensate for the dynamics of this model, i.e., to render the system response close to the desired displacement, inverse control can help achieve this objective; that is, the command is evaluated by means of

yac;i ¼ ^ha1 ya;i þ ^ha2 ya;i

- 1

þ ^ha3 ya;i

ð17Þ

- 2

where ya represents the desired displacement of the adaptive compensator. If ym;i1 ¼ ya;i1 ; ym;i1 ¼ ya;i1 , then ym;i ¼ ya;i . This result reveals the reason why this scheme is effective for ideal cases; for more general cases, the reason is explained herein. Subtracting Eq. (17) from Eq. (16) gives

ym;i  ya;i ¼ 

h2 h3 ½y  ya;i1   ½ym;i2 ¼ ya;i2  h1 m;i1 h1

ð18Þ

This equation can be recast as

ei ¼ a1 ei

- 1

þ a2 ei

ð19Þ

- 2

where ei ¼ ym;i - ya;i ; a1 ¼ 



ei ei

- 1



 ¼

a1 a2 1

0



h2 h1

; a2 ¼ 

ei

- 1

ei

- 2



h3 . h1

 ei ¼A ei

Eq. (19) is a recursive equation that can be rewritten in compact form, that is, - 1

 ð20Þ

- 2

in which A is called the amplification matrix. Similar to stability analysis of numerical integration methods, if and only if the absolute values of the eigenvalues of this amplification matrix are less than one, then the error vanishes as time approaches infinity. Fortunately, these eigenvalues can be solved analytically, namely,

k1;2

h2 1 ¼   2h1 2

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi h22 4h3  h1 h21

ð21Þ

8

Z. Wang et al. / Mechanical Systems and Signal Processing 133 (2019) 106262

Hence, to ensure convergence, this condition k1;2 < 1 has to be checked. This analysis shows that under this condition, the achieved displacement approaches the desired displacement, indicating complete compensation for the system dynamics. To obtain the model parameters, the least squares method presented in Section 3 is also adopted. In addition, for ideal cases where the delay is completely compensated by PE, the model parameters in Eq. (17) are [1 0 0]T, which is used for parameter initialization.

5. Virtual RTHS To validate the performance of the proposed controller and filter, numerous virtual tests are performed on the benchmark problem, of which results are presented in this and the next sections. Three earthquakes and 4 substructure partitioning cases provided by the benchmark model are taken into account. Each partitioning case is performed for 20 runs with perturbed system parameters to consider the robustness, along with one run using nominal system parameters. 5.1. System model The Simulink model of the proposed method is shown in Fig. 6. Discrete PI control with Kp = 1.65 and Ki = 95 is performed as the plant controller. As mentioned in Section 3, one second-order transfer function for the loading system is identified, and then the Kalman filter is designed to suppress noise. To accurately recover the real specimen displacement, adaptive compensation to this filter output is added using the least squares method with the forgetting factor q = 1.0. To compensate for the loading system delay and dynamics, PE for a delay of 10.3 ms is performed as course compensation, while adaptive compensation based on the least squares method with the forgetting factor q = 1.0 is performed as fine compensation. Another extrapolation for the computation delay of sampling time is also performed, whose necessity is elaborated in the next section. 5.2. Results Table 1 lists the criteria results of virtual RTHS using the proposed adaptive two-stage compensation together with the adaptive filter. J1 to J9 are evaluation criteria defined in [31]. In the table, Nominal denotes the case where nominal loading and specimen parameters are deployed. Mean and STD represent the mean and standard deviation values of each evaluation criterion in the 20 runs of the partitioning case using perturbed parameters for robustness evaluation, respectively. Average and Maximum stand for the average and maximum values of the criterion of all cases, respectively. For example, Average Nominal J2 is the average of 12 Nominal values of 4 cases in that column, i.e., 0.76, 0.78, 0.61, 0.70, 0.89 and so on. From the definitions of these criteria, one can see that J1, J2 and J3 reveal the actuator control performance, while all other criteria indicate test errors relative to reference responses. Since the delay originates in not only actuator dynamics but also computation models, as revealed in Section 6, the first three criteria cannot clearly reflect the test performance. From Maximum Nominal results, the largest value of nine criteria in simulations using nominal parameters is J5 = 1.86%. Analogously, the largest mean and STD values of nine criteria in simulations using perturbed parameters are, J4 = 3.32%, and J7 = 1.86%, respectively. From the Average results, the largest Nominal, Mean and STD values corresponding to nine criteria are J4 = 1.52%, J4 = 2.31%, and J6 = J7 = 1.16%, respectively. Especially for J6 to J9, the Average results in simulations using nominal parameters are J6 = J7 = 1.18%, J8 = 0.85%, J9 = 0.86%, which are on the same level as the results of the ideal cases described in Section 6. These results show the outstanding tracking performance and excellent robustness of this method.

Fig. 6. Simulink model of the proposed method.

9

Z. Wang et al. / Mechanical Systems and Signal Processing 133 (2019) 106262 Table 1 Criteria of virtual RTHS using the proposed method. Earthquake

Case No.

Type

J1 (ms)

J2 (%)

J3 (%)

J4 (%)

J5 (%)

J6 (%)

J7 (%)

J8 (%)

J9 (%)

El Centro

Case 1

Nominal Mean STD Nominal Mean STD Nominal Mean STD Nominal Mean STD

0.00 0.09 0.10 0.20 0.07 0.10 0.20 0.20 0.00 0.20 0.19 0.04

0.76 0.79 0.04 0.78 0.77 0.02 0.61 0.62 0.02 0.70 0.70 0.02

1.03 1.03 0.05 0.96 1.04 0.06 0.90 0.88 0.05 0.89 0.91 0.06

1.43 2.24 0.73 1.30 1.95 1.22 1.41 1.92 0.88 1.72 2.47 1.42

1.72 2.04 0.58 1.18 1.67 0.82 1.50 1.85 0.59 1.86 2.29 0.93

1.10 1.96 0.76 0.94 1.63 1.25 1.09 1.64 0.90 1.51 2.24 1.47

1.10 1.96 0.76 0.95 1.64 1.25 1.09 1.64 0.90 1.51 2.25 1.47

0.85 1.58 0.60 0.88 1.29 0.82 0.77 1.28 0.69 1.13 1.73 1.01

0.85 1.60 0.61 0.89 1.30 0.82 0.75 1.31 0.69 1.16 1.77 1.04

Nominal Mean STD Nominal Mean STD Nominal Mean STD Nominal Mean STD

0.00 0.07 0.10 0.20 0.12 0.10 0.20 0.18 0.06 0.20 0.12 0.10

0.89 0.91 0.05 0.78 0.80 0.03 0.69 0.71 0.03 0.81 0.85 0.05

0.57 0.63 0.05 0.58 0.62 0.05 0.48 0.53 0.04 0.58 0.65 0.06

1.33 1.65 0.52 1.48 1.92 0.60 1.22 1.87 0.76 1.67 2.07 0.79

0.85 1.11 0.39 1.22 1.44 0.41 0.80 1.38 0.64 1.06 1.34 0.46

0.87 1.27 0.57 1.09 1.59 0.64 0.90 1.67 0.80 1.29 1.74 0.84

0.87 1.27 0.57 1.09 1.59 0.64 0.89 1.66 0.79 1.29 1.74 0.84

0.57 0.89 0.42 0.63 1.02 0.45 0.48 1.19 0.66 0.79 1.08 0.48

0.57 0.88 0.41 0.63 1.01 0.44 0.50 1.19 0.66 0.79 1.08 0.48

Nominal Mean STD Nominal Mean STD Nominal Mean STD Nominal Mean STD

0.00 0.05 0.09 0.00 0.04 0.08 0.00 0.09 0.10 0.20 0.16 0.08

1.05 1.02 0.06 1.03 1.02 0.06 0.98 0.99 0.05 0.88 0.92 0.11

0.97 0.98 0.09 0.79 0.79 0.05 0.75 0.76 0.06 0.85 0.87 0.09

1.58 2.95 1.57 1.67 2.53 1.38 1.68 2.83 1.60 1.76 3.32 1.77

1.34 2.56 1.39 1.36 1.99 1.07 1.25 1.92 0.96 1.48 2.96 1.59

1.21 2.66 1.66 1.35 2.24 1.49 1.33 2.54 1.70 1.52 3.06 1.85

1.21 2.67 1.66 1.35 2.25 1.50 1.33 2.55 1.70 1.52 3.06 1.86

0.99 2.26 1.36 0.94 1.66 1.14 1.02 1.66 1.02 1.19 2.65 1.62

1.03 2.25 1.35 0.92 1.67 1.13 1.01 1.66 1.03 1.20 2.65 1.62

Average

Nominal Mean STD

0.12 0.12 0.08

0.83 0.84 0.05

0.78 0.81 0.06

1.52 2.31 1.10

1.30 1.88 0.82

1.18 2.02 1.16

1.18 2.02 1.16

0.85 1.53 0.86

0.86 1.53 0.86

Maximum

Nominal Mean STD

0.00 0.04 0.10

1.05 1.02 0.11

1.03 1.04 0.09

1.76 3.32 1.77

1.86 2.96 1.59

1.52 3.06 1.85

1.52 3.06 1.86

1.19 2.65 1.62

1.20 2.65 1.62

Case 2

Case 3

Case 4

Kobe

Case 1

Case 2

Case 3

Case 4

Morgan

Case 1

Case 2

Case 3

Case 4

The displacement responses of the first floor of the virtual RTHS are depicted in Fig. 7 in comparison with the reference. Moreover, Fig. 7i) shows a close-up view of Fig. 7a). In the test, nominal plant parameters and the El Centro earthquake are adopted together with the proposed adaptive filter and adaptive two-stage delay compensation. In the figure, the solid lines labeled Discrepancy in the legend indicate the difference between the reference and the test displacement. Fig. 7a)–d) show that for all cases, the test responses match the reference very well. As shown in Fig. 7e)–h), the largest discrepancy is approximately 0.09 mm, indicating the remarkable agreement and excellent tracking performance of the proposed schemes. The close-up view in Fig. 7i) shows good agreement as well. Virtual RTHS using chirp excitation has also been implemented but is not presented herein to save space. These simulations show that the proposed method performs relatively well for an excitation with a frequency range of [0.1 40] Hz.

6. Comparisons with conventional methods This section describes the virtual RTHS results of various scenarios, including an ideal case and cases compensated by classic PE and conventional adaptive delay compensation. These simulations and results provide a reference for evaluating the proposed compensation scheme.

10

Z. Wang et al. / Mechanical Systems and Signal Processing 133 (2019) 106262

Fig. 7. Displacement comparisons of the first floor subjected to the El Centro earthquake.

6.1. Ideal results An investigation of the ideal model, where no actuator and measurement noise are considered, is significant since it can provide model examination and limitations of any control schemes, namely, possible excellent results. Hence, before PE has been implemented for delay compensation, simulations are carried out on a spring specimen where the transfer system of

Z. Wang et al. / Mechanical Systems and Signal Processing 133 (2019) 106262

11

Fig. 7 (continued)

the plant is assumed to be 1. Criteria corresponding to this case are presented in Table 2, indicating considerable errors between the testing and reference responses. To reveal the error source, a sampling time delay, i.e., 1/4096 s, has been compensated for by means of the PE method, with the results also shown in Table 2. Clearly, the criteria of J4 to J9 have been substantially reduced by this compensation, while the criteria of J1 to J3 show discrepancies between the desired and measured displacements. One can conclude from these two simulations that (i) this Simulink model contains a delay of a sampling time, called computational delay in this paper; (ii) the criteria of J1 to J3 are designed to evaluate loading performance, and J1 = J2 = J3 = 0 does not mean perfect testing results with this computational delay in mind; and (iii) the criteria of J4 to J9 are more straightforward for evaluating test results. It is worth mentioning that the code for evaluating J1 in the m file entitled ‘‘F4_evaluation.m” was redesigned to accommodate both negative and positive delays. Finally, simulations of an ideal case with a dynamic specimen and PE compensation for a delay of 1/4096 s are performed, as shown in the last row in Table 2. These results differ little from those of the above case, and this result can be attributed to the small mass of the dynamic specimen. Since this case, where the transfer function of the loading system is assumed to be 1, is the limitation of the loading system control, the results presented here can be regarded as an ideal reference for virtual RTHS. Compared with the results in Table 1 associated with the El Centro earthquake and the nominal parameters of Case 1, the proposed scheme performs very well. This analysis also reveals the necessity of compensating for a sampling interval.

6.2. Comparison with polynomial extrapolation Prior to implementation of PE compensation, the system delay has to be measured accurately. In addition, the system delay varies due to system parameter uncertainty. In this benchmark problem, parameters, including actuator parameters

12

Z. Wang et al. / Mechanical Systems and Signal Processing 133 (2019) 106262

Table 2 Simulation results of virtual RTHS of ideal cases. Case characteristics

J1 (ms)

J2 (%)

J3 (%)

J4 (%)

J5 (%)

J6 (%)

J7 (%)

J8 (%)

J9 (%)

Spring specimen Spring & compensation Spring, mass & compensation

0.0 0.2 0.2

0.00 0.54 0.54

0.00 0.54 0.57

3.05 1.03 1.00

1.53 0.57 0.56

3.06 1.00 0.96

3.07 1.00 0.96

1.53 0.50 0.48

1.53 0.50 0.48

and specimen parameters, are taken with variations. To measure the nominal system delay, loading simulations for the four cases considering nominal parameters are carried out on the benchmark problem together with an adaptive filter for loading the emulated structural responses. J1 is evaluated for each run and averaged, resulting in a measured delay of 10.3 ms. This delay is compensated by the PE method in this subsection and in Section 5. Virtual RTHS is compensated by the PE method for both the loading delay of 10.3 ms and the computational delay of 1/4096 s mentioned in Subsection 6.1 and is carried out, with the statistics of the criteria results presented in Table 3. From the Maximum results, the largest values of all criteria corresponding to Nominal, Mean and STD are J4 = 6.74%, J4 = 7.48%, and J4 = 3.62%, respectively. From the Average results, the largest values of all criteria corresponding to Nominal, Mean and STD are J4 = 4.38%, J4 = 4.66%, and J4 = 1.75%, respectively. Compared with these results, the proposed delay compensation whose simulation results are detailed in Table 1 performs more satisfactorily.

6.3. Comparison with single-stage adaptive delay compensation Single-stage adaptive delay compensation (SADC) is implemented in this subsection. For this case, only the adaptive delay compensation described in Section 4.2 is implemented to compensate for the loading system delay, and this approach is different from the proposed method, which contains not only adaptive delay compensation but also PE compensation. For this adaptive compensation, the initial estimated parameters play a significant role in the accurate testing results. For the proposed method, the system delay is estimated offline and utilized for PE. To be fair, this estimated system delay has to be used for this SADC. Following a similar procedure to determine the initial parameters for the adaptive filter, one can obtain the formula of the initial SADC parameters, namely,

^ha1 ¼ 1 þ

s Dt

þ

s2 2 Dt

2

^ha2 ¼ 

s Dt



s2 Dt

2

^ha3 ¼

s2

ð22Þ

2 Dt 2

Eq. (22) provides applicable initial values of the parameters; however, these values are not accurate in that the real system often differs from a pure delay system, and truncation and algorithmic errors influence the parameter values. Virtual RTHS compensated by SADC was performed with the initial parameters evaluated by Eq. (22) for a delay of 10.3 ms. Note that the PE compensation for a sampling time described in Subsection 6.1 has also been carried out. The results are presented in Table 4. Additionally, to save space, only the statistics are listed in this table. Generally, the nominal and mean values in this table are greater than the counterpart values in Table 1, indicating the improvement due to the two-stage strategy. Additionally, from the Maximum results, the largest values of all criteria corresponding to Nominal, Mean and STD are J5 = 3.80%, J5 = 4.07%, and J6 = J7 = 1.04%, respectively. From the Average results, the largest values of all criteria corresponding to Nominal, Mean and STD are J5 = 2.77%, J4 = 3.08%, and J6 = J7 = 0.54%, respectively. The reason why the proposed adaptive two-stage compensation outperforms the SADC is elaborated here. For pure delay, the Taylor expansion truncation error depends on the number of expansion terms and the delay value. This effect means that for the same expansion, the truncation error is a function of the delay values. The larger the delay is, the greater the truncation error is. In other words, the proposed method is endowed with smaller model error due to the smaller delay needed to be compensated by the adaptive scheme. Additionally, the adaptive compensation for the proposed scheme needs to treat delay variation and compensation errors of PE. SADC has to update parameters for delay variation and compensation error for initial parameters. These initial parameters are often not as effective as the PE method. In summary, the proposed delay compensation simplifies the parameter initialization of the adaptive process, reduces the model error, and hence is superior to the conventional SADC method.

Table 3 Statistics of the evaluation criteria of virtual RTHS compensated by PE. Item

Type

J1 (ms)

J2 (%)

J3 (%)

J4 (%)

J5 (%)

J6 (%)

J7 (%)

J8 (%)

J9 (%)

Average

Nominal Mean STD

0.20 0.16 0.14

1.15 1.19 0.19

1.22 1.27 0.22

4.38 4.66 1.75

2.98 3.19 1.20

3.91 4.20 1.54

3.92 4.21 1.55

2.49 2.71 1.03

2.48 2.71 1.03

Maximum

Nominal Mean STD

0.20 0.23 0.20

1.36 1.42 0.38

1.48 1.56 0.37

6.74 7.48 3.62

5.80 6.34 3.06

6.12 6.87 3.19

6.13 6.89 3.21

5.14 5.66 2.64

5.14 5.65 2.65

13

Z. Wang et al. / Mechanical Systems and Signal Processing 133 (2019) 106262 Table 4 Statistics of the evaluation criteria of virtual RTHS compensated by SADC. Item

Type

J1 (ms)

J2 (%)

J3 (%)

J4 (%)

J5 (%)

J6 (%)

J7 (%)

J8 (%)

J9 (%)

Average

Nominal Mean STD

0.20 0.23 0.08

1.55 1.55 0.07

1.76 1.77 0.08

2.71 3.08 0.53

2.77 3.05 0.41

1.86 2.32 0.54

1.87 2.33 0.54

1.52 1.90 0.44

1.53 1.91 0.44

Maximum

Nominal Mean STD

0.20 0.13 0.15

1.97 1.98 0.14

2.01 2.01 0.15

3.65 3.76 0.89

3.80 4.07 0.88

2.82 3.03 1.04

2.83 3.04 1.04

2.31 2.44 0.91

2.32 2.48 0.89

7. Conclusions Aiming to achieve high performance control of RTHS, this paper conceives a noise filter and delay compensation based on an adaptive inverse strategy. In particular, to address measurement noise, an adaptive filter consisting of a Kalman filter and its compensation using an adaptive inverse strategy is proposed and examined. Moreover, adaptive two-stage delay compensation is presented, and this approach applies polynomial extrapolation (PE) for the main delay part and adaptive delay compensation for the remaining delay and errors of the PE. The strategy behind this method is to coarsely compensate delay with the PE method and finely compensate with an adaptive strategy. Virtual RTHS considering 4 partitioning cases, 3 earthquakes and nominal and perturbed parameters is carried out using the proposed delay compensation, conventional PE and single-stage adaptive compensation. The results show that the proposed method outperforms the other two schemes and achieves nearly ideal results for nominal cases. Specifically, for the proposed delay compensation, the largest values of all criteria for all cases corresponding to Nominal, Mean and STD are J5 = 1.86%, J4 = 3.32%, and J7 = 1.86%, respectively. For the Average, the largest values of all criteria corresponding to Nominal, Mean and STD are J4 = 1.52%, J4 = 2.31%, and J6 = J7 = 1.16%, respectively. Acknowledgments The authors gratefully acknowledge the financial support of the National Key Research and Development Program of China, China (Grant No. 2016YFC0701106) and the National Natural Science Foundation of China, China (Grant Nos. 51878525, 51778190, and 51408157). This study is also supported by the Foundation of Key Laboratory of Structures Dynamic Behavior and Control (Ministry of Education) at the Harbin Institute of Technology, China. References [1] J.-T. Wang, Y. Gui, F. Zhu, F. Jin, M.-X. Zhou, Real-time hybrid simulation of multi-story structures installed with tuned liquid damper, Struct. Control Heal. Monit. 23 (2016) 1015–1031, https://doi.org/10.1002/stc.182. [2] T. Sauder, V. Chabaud, M. Thys, E.E. Bachynski, L.O. Sæther, Real-time hybrid model testing of a braceless semi-submersible wind turbine: part I — the hybrid approach p. V006T09A039, in: Vol. 6 Ocean Sp. Util. Ocean Renew. Energy, ASME, 2016, https://doi.org/10.1115/OMAE2016-54435. [3] S.-H. Eem, J.-H. Koo, H.-J. Jung, Feasibility study of an adaptive mount system based on magnetorheological elastomer using real-time hybrid simulation, J. Intell. Mater. Syst. Struct. 1045389X1875434 (2018), https://doi.org/10.1177/1045389X18754347. [4] W. Su, W. Song, V. Hill, Real-time hybrid simulation and experiment for aeroelastic testing of flexible wings, in: AIAA Scitech 2019 Forum, American Institute of Aeronautics and Astronautics, Reston, Virginia, 2019, https://doi.org/10.2514/6.2019-2032. [5] Z. Zhang, B. Basu, S.R.K. Nielsen, Real-time hybrid aeroelastic simulation of wind turbines with various types of full-scale tuned liquid dampers, Wind Energy. 22 (2019) 239–256, https://doi.org/10.1002/we.2281. [6] S.N. Dermitzakis, S.A. Mahin. Development of substructuring techniques for on-line computer controlled seismic performance testing. UCB/EERC-85/ 04. [7] M. Nakashima, H. Kato, E. Takaoka, Development of real-time pseudo dynamic testing, Earthq. Eng. Struct. Dyn. 21 (1992) 79–92, https://doi.org/ 10.1002/eqe.4290210106. [8] C.Y. Chae, M. Park, C.-Y. Kim, Y.S. Park, Experimental study on the rate-dependency of reinforced concrete structures using slow and real-time hybrid simulations, Eng. Struct. 132 (2017) 648–658, https://doi.org/10.1016/J.ENGSTRUCT.2016.11.065. [9] B. Wu, P. Shi, Q. Wang, X. Guan, J. Qu, Performance of an offshore platform with MR dampers subjected to ice and earthquake, Struct. Control Heal. Monit. 18 (2011) 682–697, https://doi.org/10.1002/stc.398. [10] X. Shao, J. van de Lindt, P. Bahmani, W. Pang, E. Ziaei, M. Symans, J. Tian, T. Dao, Real-time hybrid simulation of a multi-story wood shear wall with first-story experimental substructure incorporating a rate-dependent seismic energy dissipation device, Smart Struct. Syst. 14 (2014) 1031–1054, https://doi.org/10.12989/sss.2014.14.6.1031. [11] Y. Chae, J.M. Ricles, R. Sause, Large-scale real-time hybrid simulation of a three-story steel frame building with magneto-rheological dampers, Earthq. Eng. Struct. Dyn. 43 (2014) 1915–1933, https://doi.org/10.1002/eqe.2429. [12] T. Horiuchi, M. Inoue, T. Konno, Y. Namita, Real-time hybrid experimental system with actuator delay compensation and its application to a piping system with energy absorber, Earthq. Eng. Struct. Dyn. 28 (1999) 1121–1141, https://doi.org/10.1002/(SICI)1096-9845(199910)28:10<1121::AIDEQE858>3.0.CO;2-O. [13] L. Huang, C. Chen, T. Guo, M. Chen, Stability analysis of real-time hybrid simulation for time-varying actuator delay using the lyapunov-krasovskii functional approach, J. Eng. Mech. 145 (2019) 04018124, https://doi.org/10.1061/(ASCE)EM.1943-7889.0001550. [14] M.I. Wallace, J. Sieber, S.A. Neild, D.J. Wagg, B. Krauskopf, Stability analysis of real-time dynamic substructuring using delay differential equation models, Earthq. Eng. Struct. Dyn. 34 (2005) 1817–1832, https://doi.org/10.1002/eqe.513. [15] A.P. Darby, M.S. Williams, A. Blakeborough, Stability and delay compensation for real-time substructure testing, J. Eng. Mech. 128 (2002) 1276–1284, https://doi.org/10.1061/(ASCE)0733-9399(2002) 128:12(1276). [16] M. Nakashima, N. Masaoka, Real-time on-line test for MDOF systems, Earthq. Eng. Struct. Dyn. 28 (1999) 393–420, https://doi.org/10.1002/(SICI)10969845(199904)28:4<393::AID-EQE823>3.0.CO;2-C.

14

Z. Wang et al. / Mechanical Systems and Signal Processing 133 (2019) 106262

[17] P.A. Bonnet, C.N. Lim, M.S. Williams, A. Blakeborough, S.A. Neild, D.P. Stoten, C.A. Taylor, Real-time hybrid experiments with Newmark integration, MCSmd outer-loop control and multi-tasking strategies, Earthq. Eng. Struct. Dyn. 36 (2007) 119–141, https://doi.org/10.1002/eqe.628. [18] X. Ning, Z. Wang, H. Zhou, B. Wu, Y. Ding, X. Bin, Robust actuator dynamics compensation method for real-time hybrid simulation, Mech. Syst. Signal Process 131 (2019) 49–70. [19] M. Ahmadizadeh, G. Mosqueda, A.M. Reinhorn, Compensation of actuator delay and dynamics for real-time hybrid structural simulation, Earthq. Eng. Struct. Dyn. 37 (2008) 21–42, https://doi.org/10.1002/eqe.743. [20] Z. Wang, B. Wu, O.S. Bursi, G. Xu, Y. Ding, An effective online delay estimation method based on a simplified physical system model for real-time hybrid simulation, Smart Struct. Syst. 14 (2014) 1247–1267, https://doi.org/10.12989/sss.2014.14.6.1247. [21] T.V. Nguyen, E.U. Dorka, Phase lag compensation in real-time substructure testing based on online system identification, Proc. 14th World Conf. Earthq. Eng., 2008, pp. 12-01. [22] Y. Chae, K. Kazemibidokhti, J.M. Ricles, Adaptive time series compensator for delay compensation of servo-hydraulic actuator systems for real-time hybrid simulation, Earthq. Eng. Struct. Dyn. 42 (2013) 1697–1715, https://doi.org/10.1002/eqe.2294. [23] C. Chen, J.M. Ricles, T. Guo, Improved adaptive inverse compensation technique for real-time hybrid simulation, J. Eng. Mech. 138 (2012) 1432–1446, https://doi.org/10.1061/(ASCE)EM.1943-7889.0000450. [24] H. Zhou, D.J. Wagg, M. Li, Equivalent force control combined with adaptive polynomial-based forward prediction for real-time hybrid simulation, Struct. Control Heal. Monit. 24 (2017) e2018, https://doi.org/10.1002/stc.2018. [25] B. Wu, Z. Wang, O.S. Bursi, Actuator dynamics compensation based on upper bound delay for real-time hybrid simulation, Earthq. Eng. Struct. Dyn. 42 (2013) 1749–1765, https://doi.org/10.1002/eqe.2296. [26] S. Strano, M. Terzo, Actuator dynamics compensation for real-time hybrid simulation: an adaptive approach by means of a nonlinear estimator, Nonlinear Dyn. 85 (2016) 2353–2368, https://doi.org/10.1007/s11071-016-2831-0. [27] X. Gao, N. Castaneda, S.J. Dyke, Real time hybrid simulation: from dynamic system, motion control to experimental error, Earthq. Eng. Struct. Dyn. 42 (2013) 815–832, https://doi.org/10.1002/eqe.2246. [28] G. Ou, A.I. Ozdagli, S.J. Dyke, B. Wu, Robust integrated actuator control: experimental verification and real-time hybrid-simulation implementation, Earthq. Eng. Struct. Dyn. 44 (2015) 441–460, https://doi.org/10.1002/eqe.2479. [29] N. Nakata, Effective force testing using a robust loop shaping controller, Earthq. Eng. Struct. Dyn. 42 (2013) 261–275, https://doi.org/10.1002/eqe.2207. [30] B. Wu, H. Zhou, Sliding mode for equivalent force control in real-time substructure testing, Struct. Control Heal. Monit. 21 (2014) 1284–1303, https:// doi.org/10.1002/stc.1648. [31] C.E. Silva, D. Gomez, A. Maghareh, S.J. Dyke, B.F. Spencer, Benchmark control problem for real-time hybrid simulation, Mech. Syst. Signal Process. Under revi (2019). [32] C.K. Chui, G. Chen, Kalman Filtering with Real-time Applications (5th edn), Springer, 2017. [33] P. Shi, B. Wu, B.F. Spencer, B.M. Phillips, C.-M. Chang, Real-time hybrid testing with equivalent force control method incorporating Kalman filter, Struct. Control Heal. Monit. 23 (2016) 735–748, https://doi.org/10.1002/stc.1808. [34] T. SÖDERSTRÖM, P. STOICA, System identification, Prentice Hall International, Page 324, 1989.