Applications of the RST Algorithm to Nonlinear Systems in Real-Time Hybrid Simulation

Real-time substructure testing (RST) algorithm is a newly developed integration method for real-time hybrid simulation (RTHS) which has structure-dependent and explicit formulations for both displacement and velocity. (e most favourable characteristics of the RSTalgorithm is unconditionally stable for linear and no iterations are needed. In order to fully evaluate the performance of the RST method in solving dynamic problems for nonlinear systems, stability, numerical dispersion, energy dissipation, and overshooting properties are discussed. Stability analysis shows that the RSTmethod is only conditionally stable when applied to nonlinear systems. (e upper stability limit increases for stiffness-softening systems with an increasing value of the instantaneous degree of nonlinearity while decreases for stiffness-hardening systems when the instantaneous degree of nonlinearity becomes larger. Meanwhile, the initial damping ratio of the system has a negative impact on the upper stability limit especially for instantaneous stiffness softening systems, and a larger value of the damping ratio will significantly decrease the upper stability limit of the RST method. It is shown in the accuracy analysis that the RST method has relatively smaller period errors and numerical damping ratios for nonlinear systems when compared with other two well-developed algorithms. (ree simplified engineering cases are presented to investigate the dynamic performance of the RSTmethod, and the numerical results indicate that this method has a more desirable accuracy than other methods in solving dynamic problems for both linear and nonliner systems.


Introduction
Real-time hybrid simulation (RTHS) is an advanced technique for assessing the seismic behaviours of structures, especially large-scale and ultralimit complex systems [1][2][3][4]. Similar to other hybrid simulation techniques, the original system in a RTHS is divided into two components, including numerical substructures and experimental substructures. e displacement-coordination and force-equilibrium are two principles for the coupling between two components at the interfaces [5]. Both the restoring force calculated from numerical substructures and the resisting force measured from experimental substructures at each time step are fed back to an integration algorithm. en, the subsequent timestep target displacement applied to the experimental substructure is calculated by solving the step-by-step equation of motion under external excitations and the reacting forces measured from experimental substructures [6]. "Real-time" is a significant characteristic that makes RTHS different from other hybrid simulations, and the corresponding applications include "real-time calculating of numerical substructures", "real-time loadings on experiment substructures," and "realtime transferring the interactions between these substructures." So far, a few works have focused on the latter two applications. However, with an increasing size and complexity of structural systems, "real-time calculating" becomes more important and the requirement of numerical integration algorithm is higher.
In RTHS, integration algorithms are effective methods to solve discrete equations of motion for structural dynamics [7]. Integration algorithms are usually classified as explicit and implicit based on the expressions of the displacement and the velocity. In a pseudodynamic testing that quasistatically imposes the target displacement on the experimental substructures, an integration algorithm is explicit when the displacement for the next time step can be obtained from the response quantities (displacement, velocity, and acceleration) at the current and previous time steps, and otherwise it is implicit [8]. However, the definitions of explicit and implicit algorithms have changed in RTHS, as the loading rate becomes faster or real-time and the velocity of experimental substructures should be feedback for integration. erefore, an integration algorithm is explicit when both the displacement and velocity for the next time step can be determined as a function of the response quantities at the current and previous time steps. Although various integration algorithms have been used to conduct RTHS, explicit integration algorithms [9][10][11][12][13][14][15][16][17][18] are generally preferred over implicit ones [19][20][21][22][23] as they do not involve nonlinear iterations and show a higher calculation efficiency. However, conditional stability limits the application of conventional explicit algorithms in solving structural dynamic problems because it may lead to an extremely small step size for largenumber multidegree-of-freedom (MDOF) systems which have high frequency modes.
More and more researchers focus on developing new types of integration algorithms, which can achieve both explicit formulation and unconditional stability at the same time. Chang [24] developed an explicit integration algorithm which is referred to as Chang explicit method (CEM) for pseudodynamic testing, where the expression of displacement is an explicit form. However, when used in RTHS, CEM becomes semiexplicit since the determination of the velocity for the next time step requires knowledge of the unknown acceleration at the next time step. e major difference between CEM and other algorithms lies in the coefficients of difference equations introduced in displacement increment which are functions of initial structural properties and time step rather than constants, making it a structure-dependent algorithm. Chen and Ricles [12] proposed another type of structure-dependent integration algorithm named CR method (CRM) based on discrete control theory. Different from CEM, the expressions for displacement and velocity are both in explicit forms and the coefficients for both displacement increment and velocity increment are functions of initial structural properties and time step. It is noted that CEM and CRM are unconditionally stable for linear and instantaneous stiffness softening systems and possess second order accuracy [25,26]. Kolay and Ricles [14] introduced a family of explicit direct integration algorithms using discrete control theory, which is named as explicit KR-α algorithm. is method is unconditionally stable for linear and instantaneous stiffness softening systems, and it can control the amount of numerical damping by a single parameter. However, all these favourable attributes are achieved at the expense of degraded accuracy for general structural dynamics problems [27]. Tang and Lou [15] presented a new type of structure-dependent explicit integration algorithm for RTHS using a pole-mapping rule from the discrete domain and named it real-time substructure testing (RST) algorithm. For the RST method, the structure-dependent coefficients are introduced only in displacement increment while those for velocity increment are constant. Kim and Lee [28] also proposed an explicit integration algorithm with controllable numerical dissipation in the high frequency range, and an extended stability limit and improved spectral characteristics are achieved as well. is method has second-order accuracy but is conditionally stable.
ough the performance of the RST method in solving dynamic problems for linear systems has been well studied, its application to nonlinear systems need to be further discussed and evaluated. In this work, the numerical properties of the RST method in RTHS are discussed for nonlinear structural dynamic problems and compared with other two structure-dependent integration algorithms (i.e., CEM and CRM) since they have quite similar formulations and numerical properties. Firstly, the stability and the upper stability limit of the RST method for nonlinear systems with softening and hardening stiffness are discussed using the spectra radius method. en, the numerical dispersion and energy dissipation characteristics caused by the numerical damping are evaluated. In addition, the overshooting behaviour [25,29] of the RST method in displacement is addressed as well.
ree different simplified engineering structural systems are presented as numerical examples to demonstrate the computational stability and accuracy of the RST method when compared with other algorithms.

Formulations of Structure-Dependent
Integration Algorithms 2.1. Single-Degree-of-Freedom (SDOF) Systems. e differential equation of motion for a SDOF system under external excitations can be expressed in time domain as [30,31] where m, c, and k are the mass, viscous damping, and stiffness, respectively; € u, _ u, and u are the acceleration, velocity, and displacement, respectively; r N is the restoring force of the numerical substructure, where r N � k N u(t) is for linear systems; R E (u, _ u) is the reacting force measured from the experimental substructure; subscript N and E are the numerical substructure and the experimental substructure; and f (t) is the external exciting force. When discretized by time steps, equation (1) can be written in its discrete form: where € u i , _ u i , and u i are the acceleration, velocity, and displacement at ith time step, respectively; R E (u i , _ u i ) is the reacting force of the experimental substructure at ith time step; and fi is the external exciting force at ith time step. e general expressions of the displacement and velocity increments for the RST, CEM, and CRM methods is found to be 2 Mathematical Problems in Engineering where ∆t is the chosen integration time step. e coefficients α 1 , α 2 , α 3 , and α 4 for the three methods are and ω n is the natural frequency of the system determined from the initial system.
Obviously, the coefficients in each method are not all constant, and some of them depend on the initial structural properties (ξ and ω n ) and time step ∆t. Substituting equation (6) into equation (5), it is seen that CEM is implicit as there is € u i+1 in the velocity formulation, which is unknown at (i+1)th time step. Substituting equations (4)-(6) into equation (3) gives where c and k are the viscous damping coefficient and stiffness of the global structural system, which are equal to (c N + c E ) and (k N + k E ), respectively. e velocity of the CEM method at (i + 1)th time step can be rewritten in an explicit form as follows:

Multidegree-of-Freedom (MDOF) Systems.
For a MDOF, the coupled differential equation of motion should be written in a matrix form as follows: where M, C, and K are mass matrix, viscous damping coefficients matrix, and stiffness matrix, respectively; € U i , _ U i , and U i are the acceleration vector, velocity vector, and displacement vector at ith time step, respectively; R N is the restoring force matrix of the numerical substructure, where R N � K N U i for linear systems; and R E and F i are the reacting force matrix of the experimental substructure and the external exciting force vector, respectively.
For simplicity, C is supposed to be proportional and the orthogonality of modes is available in the following equations. For a MDOF linear system, rewrite equations (4), (5), and (11) in the modal coordinate system as follows: in which, Φ � ϕ 1 ϕ 2 · · · ϕ n is the modal matrix derived from solving the eigenvalue problem; n is the number of modes; M * � Φ T MΦ, C * � Φ T CΦ, and K * � Φ T KΦ are modal mass matrix, modal diagonal damping matrix, and modal diagonal stiffness matrix, respectively; € Y, _ Y, and Y are the acceleration vector, velocity vector, and displacement vector in modal coordinates, respectively; α * j � Φ −1 α j Φ (j � 1, 2, 3, 4) are the diagonal integration coefficient matrices, and α j is expressed as follows: , Mathematical Problems in Engineering 3 where For a MDOF nonlinear system, equations (12)- (14) are not valid, and the general formulation and coefficients for the three methods should be is the damping matrix generally determined from the initial structural properties; and K 0 is the initial stiffness matrix.

Stability Analysis of the RST Method in Nonlinear Systems.
e unconditional stability of the RST method for linear elastic systems has been verified by authors [15,24,32]. In this section, stability analysis of the RST method for nonlinear systems is mainly carried out and compared with that of the CEM and CRM methods.
In order to realistically reflect the change of stiffness during the test, a parameter, defined as the ratios of the stiffness at the end of each time step k i over the initial stiffness k 0 , is introduced to monitor this change. at is, where δ i + 1 is the instantaneous degree of nonlinearity. When δ i + 1 � 1, it means that the system is linear, and no stiffness change occurs during the test. When δ i + 1 >1, it represents that the instantaneous stiffness becomes harder at the end of the (i + 1)th time step, while 0 < δ i + 1 < 1 denotes instantaneous stiffness soften. From equations (3)-(5), the computing sequence at (i + 1)th time step can be written in a recursive matrix form as follows: ] T is state vector values at the (i + 1)th time step; load factor L � [0, 0, Δt 2 /m] T ; and amplification matrix A for each algorithm is found to be Based on the algorithm stability analysis theory, a stable computation of an integration algorithm can be obtained when the spectral radius ρ(A) � max(|λ 1,2,3 |) ≤ 1, where the eigenvalues λ 1,2,3 can be determined by solving the following eigenvalue problem [31,32]: where I is the identity matrix and A 1 , A 2 , and A 3 are three coefficients indicating half of the trace, sum of principal minors, and determinant of A, respectively. After substituting equations (21) Figure 1. It illustrates that the RST method is conditionally stable for nonlinear systems except for instantaneous stiffness softening systems when ξ � 0.0. e upper stability limit [Ω] of the RST method increases with an increasing value of δ when 0 <δ < 1 while decreases with a growth of δ value when δ > 1. e initial damping ratio has a great impact on [Ω] as well especially for instantaneous stiffness softening systems, and a larger value of ξ will significantly decrease the upper stability limit of the RST method. e variations of [Ω] with δ for the CEM and CRM methods are also plotted in Figure 1, and both these methods are unconditionally stable when 0 <δ ≤ 1 while conditionally stable when δ > 1. In particular, the curves for the CRM method are overlapped together for different ξ while that for the RST and CEM methods seem to move upward as the value of ξ increases.

Accuracy Analysis
Two principal eigenvalues at the (i + 1)th time step for the three structure-dependent integration algorithms are complex conjugate and can be expressed as [33] in which ξ i+1 and ω i+1 are the numerical damping ratio and frequency related to integration algorithms. In general, numerical dispersion and energy dissipation characteristics are two indexes evaluating the accuracy of an integration algorithm. e former characteristic is usually expressed by the relative period error PE � ((T− T n )/T n ) � (ω n /ω − 1), where T � 2π/ω and T n � 2π/ω n referred to the numerical and exact periods of the systems. e later one can be expressed by the numerical damping ratio ξ, which is related to amplitude decay in one cycle of freevibration. e variations of PE are shown in Figure 2 (ξ � 0) and Figure 3 (ξ � 0.1) for different δ. It is worth noting that values of PE in all the figures are positive and increase with increasing values of Ω for a given δ, which results in period elongation for the three methods. It is also seen that the curves are totally overlapped for the RST, CEM, and CRM methods for different values of δ as ξ � 0 as well as for δ � 1 as ξ � 0.1. Meanwhile, for a specific value of Ω, the RST method possesses the minimum period elongation while the CRM method produces the maximum period elongation for δ � 0.5, 2 as ξ � 0.1. e variations of ξ with Ω for the RST, CRM, and CEM methods can be found in equation (26) with A 1 and A 2 derived from equations (21)- (24). It is noteworthy that for ξ � 0, A 2 equals to 1 for all the three methods, which leads to ρ � 1 and ξ � 0 in equation (26). is implies that no energy dissipation occurred for the three methods if a zero viscous damping ratio adopted. On the contrary, there will be no amplitude decay for an undamped free vibration. When viscous damping ratio is considered, ρ ∞ � 1 and the variations of ξ with Ω for different δ as ξ � 0.1 are shown in Figure 4. For a linear elastic system with viscous damping ratio (δ �1), ξ � ξ for Ω � 0 and then reduces below ξ with increasing Ω, indicating that the viscous damping ratio is not effective in reducing spurious participation of higher modes. Also, the curves for the three methods are totally overlapped in this case, as shown in Figure 4(b). Meanwhile, drastically different phenomena are found for the case of δ � 0.5 and 2.0. In Figure 4(a), the variations of ξ with Ω are as similar as the case of δ � 1, but ξ > ξ for a small value of produces undesired dissipation in lower modes. In Figure 4(c), ξ reduces far below ξ with increasing omega and then increases again in a specific value of Ω. It implies that favourable energy dissipation can be obtained to filter out the spurious participation of higher modes.

Overshooting
In general, two primary factors are taken into account for choosing an appropriate time step Δt in the step-by-step integration algorithm in order to perform well in solving structural dynamic problems. e first factor is the highest frequency f h that is noteworthy in the dynamic analysis, and its value depends both on the higher modal frequency of the system and the effective high frequency of the dynamic load [15]. e second factor is the stability requirement of an integration algorithm. Since the value of a conditionally stable algorithm is strictly limited by upper stability limit, the value of Δt can be easily determined. If an algorithm is unconditionally stable, the value of Δt depends largely on the first factor.
According to the basic theory of structural dynamics, reliable solutions can be achieved when Δt/T ≤1/10 (T �1/ f h ) [31].
at is, a large time step involving in dynamic analysis, properly speaking, a large value of the ratio of Δt to T will lead to calculation distortion or so called "overshooting." Actually, the phenomenon that an integration algorithm might overshoot the exact solution in the initial stage of the transient response or in the steadystate response when Ω ⟶ ∞ has been found and discussed in some papers [25,29]. Tang and Lou [34] provided a detailed explanation for the issue of overshooting by solving free and forced vibration problems of three Mathematical Problems in Engineering 5 single-degree-of-freedom systems with different initial stiffness (k 0 � 10 2 , 10 4 , 10 6 ) and different values of Δt/T. It is reasonable to discuss the overshooting issue in zero damping systems as undesirable contaminated responses will be eliminated due to vibration damping. It is noted that overshooting in the steady-state response depends mainly on the rationality of integration time step selection. In a real-time hybrid testing, overshooting will not appear in the entire response for any unconditionally stable integration algorithm as the value of time step Δt is strictly limited by the "real-time" requirement. erefore, the overshooting property of the three methods will not be discussed in detail herein, and the conclusions made by Tang and Lou [34] will be demonstrated by the following numerical examples.

Numerical Simulation of RTHS
In this section, three numerical examples are applied to examine the numerical properties of the RST method when compared with the CRM and CEM methods for both linear elastic and nonlinear systems.

Example 1. An undamped 4-story building
A linear 4-story shear building with no structural damping is shown in Figure 5. e initial structural properties are m 1 � 10 8 kg, m 2 � 10 5 kg, m 3 � 10 5 kg, m 4 � 10 3 kg, and k � 360 × 10 7 N/m. e natural frequencies of the system are 5.9940, 116.95, 306.59, and 1906.9 rad/s for the four modes, respectively; the normalized modal matrix is given by Φ � are calculated. e time step of Δt � 0.01s is employed for RST, CRM, and CEM, and the ratio of Δt/T 4 � 3.036. e displacement response obtained from the undamped displacement equation based on the structural dynamics theory [31] is considered as the exact solution. Defining the relative error e (%) as in which A is the peak value of responses derived from each integration algorithm; subscript "METHOD" refers to the RST, CRM, and CEM methods, while "EXACT" refers to the exact solutions, which are calculated by the Newmark method with a time step of 0.005 s for the overall structural system. It is shown in Figure 6 that all the three methods provide reliable solutions, and the relative errors of the displacement peaks are 1.46% (RST), 2.64% (CRM), and 0.37% (CEM), respectively. e three methods result in a slight period elongation and zero amplification decay, as shown in Figure 6(a), which are consistent with the previous conclusions.
e results also indicate that overshooting behaviours will not appear in a zero damping system if the time step is determined from the perspective of the highest frequency. Figure 7. In this case, the frame is taken as the numerical substructure while the vibration isolator is supposed to be the experimental substructure. e initial structural properties are m � 160 × 10 3 kg, k N � 360 × 10 7 N/ m, k E � 9 × e seismic performance of the system is studied by exciting the ground acceleration record of EL Centro (1940 NS) with peak ground acceleration (PGA) scaled to 0.8 g.

Example 2. A 3-story building with a vibration isolator A linear 3-story shear building with a vibration isolator is shown in
e system is subjected to zero initial velocity and initial displacement. e time step of Δt � 0.02 s is employed for the Example 3. A 5-story building with softening stiffness A 5-story shear building (as shown in Figure 9) with softening stiffness is applied to investigate the performance of the RST, CRM, and CEM methods in a nonlinear system. In this case, the top three stories are taken as the numerical substructure while the first two stories are supposed to be the experimental substructure. e initial structural properties are m � 10 5 kg, k 0 � 10 8 N/m, and ξ � 0.05. e softening stiffness of the system after deformation is intentionally designated as k i � k i,0 (1 − ���� |Δu i | ), where subscript i refers to the ith story and Δu i is the interstory drift for the ith story. e natural initial frequencies of the whole system are 9.001, 26 [31], the damping matrix of the system and the proportionality factors are given by e seismic responses of the system with zero initial conditions is studied by exciting the ground acceleration record of EL Centro (1940 NS) with peak ground acceleration (PGA) scaled to 0.9 g. e results obtained from N4 algorithm with a time step Δt � 0.005 s is considered as the exact solution for comparisons. e time step of Δt � 0.02 s is employed for the RST, CRM, and CEM  methods, and the ratio of Δt/T 5 � 0.1933. Seismic responses for the 5th story and the shear force-displacement curves for the 1st story are plotted in Figures 10 and 11. It is seen that all the methods can result in acceptable solutions for the nonlinear variation of the stiffness and overshooting behaviours do not appear as well. e    relative errors of the displacement peaks for the RST, CRM, and CEM methods are 3.98%, 4.05%, and 3.97%, respectively.

Conclusions
Applications of the RST method to nonlinear systems in RTHS have been evaluated and compared with other two structure-dependent explicit integration algorithms, i.e., the CRM and CEM methods, since they have quite similar formulations and numerical properties. Stability analysis indicates that the RST method is only conditionally stable when applied to nonlinear systems, and the instantaneous degree of nonlinearity will exert great influence on the upper stability limit of the method. For stiffness-softening systems, the upper stability limit increases with an increasing value of the instantaneous degree of nonlinearity, while for stiffness-hardening systems, it decreases when the instantaneous degree of nonlinearity becomes larger. Meanwhile, the initial damping ratio of the system has a negative impact on the upper stability limit especially for instantaneous stiffness softening systems, and a larger value of the damping ratio will significantly decrease the upper stability limit of the RST method. ree simplified engineering structural systems are presented as numerical examples to demonstrate the computational stability, accuracy, and overshooting behaviours of the three methods. Numerical results demonstrate that overshooting behaviours will not appear if the time step is reasonably chosen for integration. It is also illustrated that the RST and CEM method have a better accuracy than the CRM method in solving dynamic problems for both linear and nonlinear systems. However, as the velocity of the CEM method is in an implicit form, the RST method provides more benefits in computational efficiency as both the displacement and velocity are in an explicit form.

Data Availability
e data used to support the findings of this study are included within the article. e results of this paper can be verified according to the data and methods given in the paper.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.