**3. Description of the Proposed Control Methodology**

The proposed control methodology comprises the following blocks:


Figure 5 shows a block diagram of the proposed control procedure, illustrating the interconnections between the previously enumerated modules. The proposed control methodology requires measuring the following variables: table (rod) displacement *xt*, rod acceleration .. *xt*, spool position *ysp*, pressure at actuator's chambers *P*<sup>1</sup> and *P*2, pressure and return pressures *PS* and *PR* at the servovalve's manifold and estimating the value of rod velocity, . *xt* [24,25].

**Figure 5.** Control architecture block diagram.

#### *3.1. Feedback Linearization*

Given a state-space representation of a multiple-input-multiple-output non-linear system: . *x* <sup>=</sup> *f*(*x*, *u*); *y* <sup>=</sup> *g*(*x*), the aim of the feedback linearization scheme is to find a state transformation, *z* <sup>=</sup> *z*(*x*) and an input transformation, *u* <sup>=</sup> *u*(*x*, *v*), such that the non-linear system is transformed into an equivalent linear system of the form . *z* <sup>=</sup> *Az* <sup>+</sup> *Bv* [26,27].

In the case under study, the non-linearities are present in servovalve flow-pressure (Equations (2) and (3)) and in chambers pressure time evolution expressions (Equations (4) and (5)). In this work, a direct approach has been employed to work out feedback linearization transformation.

By rearranging Equations (4) and (5), leaving the pressure derivatives on the left hand side, assuming β<sup>1</sup> = β<sup>2</sup> = β, and subtracting them, the time derivative of the pressure force acting on rod can be casted as:

$$\begin{aligned} \{\dot{P}\_1 - \dot{P}\_2\} A\_{\text{uv}} &= -\beta A\_{\text{uv}} \Big( \frac{Q\_1}{v\_{01} + A\_{\text{u}} x\_{\text{p}}} + \frac{Q\_2}{v\_{02} - A\_{\text{u}} x\_{\text{p}}} \Big) \\ &- \beta A\_{\text{uv}}^2 \dot{x}\_p \Big( \frac{1}{v\_{01} + A\_{\text{u}} x\_{\text{p}}} + \frac{1}{v\_{02} - A\_{\text{u}} x\_{\text{p}}} \Big) \\ &- \beta A\_{\text{u}} C\_{112} (P\_1 - P\_2) \Big( \frac{1}{v\_{01} + A\_{\text{u}} x\_{\text{p}}} + \frac{1}{v\_{02} - A\_{\text{u}} x\_{\text{p}}} \Big) \\ &- \beta A\_{\text{u}} \left( \frac{C\_{118} (P\_1 - P\_{B1})}{v\_{01} + A\_{\text{u}} x\_{\text{p}}} - \frac{C\_{228} (P\_2 - P\_{R2})}{v\_{02} - A\_{\text{u}} x\_{\text{p}}} \right) \end{aligned} \tag{10}$$

Now, for the case *ysp* ≥ 0, if Equations (2) and (3) are substituted in Equation (10), the following expression is obtained:

$$\left(\dot{P}\_1 - \dot{P}\_2\right)A\_{\text{\textquotedblleft}} = \ y\_{\text{sp}} \cdot F\_1^+ \left(P\_{1\prime}, P\_{2\prime}, P\_{\text{S}\prime}, P\_{\text{R}\prime}, \mathbf{x}\_p\right) + F\_2^+ \left(\mathbf{x}\_{p\prime} \cdot \dot{\mathbf{x}}\_p\right) + F\_3^+ \left(P\_{1\prime}, P\_{2\prime}, \mathbf{x}\_p\right),\tag{11}$$

where:

$${}\_{1}F\_{1}^{+}({}\_{1}P\_{2},P\_{S},P\_{R},\mathbf{x}\_{p}) = \beta \mathbf{A}\_{w} \mathbf{C}\_{d} \mathbf{K}\_{w} \left( \frac{\operatorname{sgn}({P\_{s}-P\_{1}})\sqrt{2|P\_{s}-P\_{1}|/\rho}}{v\_{01} + A\_{w}\mathbf{x}\_{p}} + \frac{\operatorname{sgn}({P\_{2}-P\_{R}})\sqrt{2|P\_{2}-P\_{R}|/\rho}}{v\_{02} - A\_{w}\mathbf{x}\_{p}} \right). \tag{12}$$

$$F\_2^+ \left( \mathbf{x}\_{p\prime} \cdot \dot{\mathbf{x}}\_p \right) = -\beta A\_w^2 \dot{\mathbf{x}}\_p \Big( \frac{1}{\upsilon\_{01} + A\_w \mathbf{x}\_p} + \frac{1}{\upsilon\_{02} - A\_w \mathbf{x}\_p} \Big) \tag{13}$$

and

$$F\_{\mathcal{B}}^{+}\left(P\_{1},P\_{2},\chi\_{\mathcal{P}}\right) = -\beta A\_{w} \left[ C\_{\mathrm{II2}} \left( P\_{1} - P\_{2} \right) \left( \frac{1}{\mathrm{u}\chi + A\_{\mathrm{B}}\chi\_{\mathcal{P}}} + \frac{1}{\mathrm{v}\chi\_{\mathrm{2}} - A\_{\mathrm{B}}\chi\_{\mathcal{P}}} \right) + \left( \frac{C\_{\mathrm{III}} \left( P\_{1} - P\_{\mathrm{II}} \right)}{\mathrm{u}\chi + A\_{\mathrm{B}}\chi\_{\mathcal{P}}} - \frac{C\_{\mathrm{III}} \left( P\_{2} - P\_{\mathrm{II}} \right)}{\mathrm{u}\chi\_{\mathcal{2}} - A\_{\mathrm{B}}\chi\_{\mathcal{P}}} \right) \right].\tag{14}$$

When *ysp* < 0, Equation (3) is substituted in Equation (10) and the time derivative of pressure force results in:

$$\dot{\mathbf{F}}\_{1}\left(\dot{\mathbf{P}}\_{1} - \dot{\mathbf{P}}\_{2}\right)\mathbf{A}\_{\mathbf{w}} = \mathbf{y}\_{\text{sp}} \cdot \mathbf{F}\_{1}^{-} \begin{Bmatrix} P\_{1}, P\_{2}, P\_{S}, P\_{R}, \ \mathbf{x}\_{p} \end{Bmatrix} + F\_{2}^{-} \begin{Bmatrix} \mathbf{x}\_{p\prime} \ \dot{\mathbf{x}}\_{p} \end{Bmatrix} + F\_{3}^{-} \begin{Bmatrix} P\_{1}, P\_{2}, \mathbf{x}\_{p} \end{Bmatrix},\tag{15}$$

where:

$$F\_1^\*(\mathcal{P}\_1, \mathcal{P}\_2, \mathcal{P}\_{\mathcal{S}}, \mathcal{P}\_{\mathcal{R}}, \mathbf{x}\_{\mathcal{P}}) = \ \beta A\_{\mathcal{W}} \mathcal{C}\_4 K\_{\mathcal{W}} \left( \frac{s\_{\mathcal{S}} n (\mathcal{P}\_1 - \mathcal{P}\_{\mathcal{R}}) \sqrt{2 \mathcal{P}\_1 - \mathcal{P}\_{\mathcal{R}} \mathcal{I} / \rho}}{\upsilon\_{\mathcal{W}} + A\_{\mathcal{W}} \mathbf{x}\_{\mathcal{P}}} + \frac{s\_{\mathcal{S}} n (\mathcal{P}\_{\mathcal{S}} - \mathcal{P}\_2) \sqrt{2 |\mathcal{P}\_{\mathcal{S}} - \mathcal{P}\_2 \mathcal{I} / \rho}}{\upsilon\_{\mathcal{W}} + A\_{\mathcal{W}} \mathbf{x}\_{\mathcal{P}}} \right), \tag{16}$$

$$F\_2^- \left( \mathbf{x}\_{\mathcal{P}}, \dot{\mathbf{x}}\_{\mathcal{P}} \right) = \; F\_2^+ \left( \mathbf{x}\_{\mathcal{P}}, \dot{\mathbf{x}}\_{\mathcal{P}} \right), \tag{17}$$

and

$$F\_3^-(P\_1, P\_2, \mathbf{x}\_p) = F\_3^+(P\_1, P\_2, \mathbf{x}\_p). \tag{18}$$

Let us now define a desired change in pressure force time derivative, *Aw*<sup>Δ</sup> . *Pdes*. Then, if the spool position *ysp* could be forced to instantaneously take the values defined by:

$$\begin{array}{lcl} \; y\_{sp} &=& \begin{cases} & \frac{A\_{\text{w}} \Delta \dot{P}\_{\text{dav}} - F\_{2}^{+} \left( \mathbf{x}\_{p}, \dot{\mathbf{x}}\_{p} \right) - F\_{3}^{+} \left( P\_{1}, P\_{2}, \mathbf{x}\_{p} \right)}{F\_{1}^{+} \left( P\_{1}, P\_{2}, P\_{S}, P\_{R}, \mathbf{x}\_{p} \right)} & \; ; \; y\_{sp, prv} \ge 0 \\ & \frac{A\_{\text{w}} \Delta \dot{P}\_{\text{dav}} - F\_{2}^{-} \left( \mathbf{x}\_{p}, \dot{\mathbf{x}}\_{p} \right) - F\_{3}^{-} \left( P\_{1}, P\_{2}, \mathbf{x}\_{p} \right)}{F\_{1}^{-} \left( P\_{1}, P\_{2}, P\_{S}, P\_{R}, \mathbf{x}\_{p} \right)} & \; ; \; y\_{sp, prv} < 0 \end{cases} \end{array} \tag{19}$$

( . *<sup>P</sup>*<sup>1</sup> <sup>−</sup> . *<sup>P</sup>*2)*Aw* <sup>=</sup> *Aw*<sup>Δ</sup> . *Pdes* would hold, and an arbitrary time variation of pressure force could be imposed on the servoactuator. Above, *ysp*,*prev* denotes the value of the spool position at the previous iteration of the real-time control system, in which feedback linearization scheme is implemented.

This is the core idea to the control procedure presented in this work: finding the required time derivative of pressure force on the actuator's rod, so that the acceleration reference profile is fulfilled. Due to the feedback linearization transformation found (Equation (19)), this value of *Aw*<sup>Δ</sup> . *Pdes* will be effectively imposed on servoactuator's rod.

In order to force the spool position to accurately track the value determined by Equation (19), a servovalve spool dynamics inversion algorithm must be implemented. Assuming that spool motion is governed by the first order system in Equation (1), the dynamics inversion can be expressed as:

$$
\mathcal{L}\_{sp}\mu\_{\text{sv},\text{FL}} = \tau\_{sp}\dot{y}\_{sp,\text{des}} + y\_{sp,\text{des}}.\tag{20}
$$

in which *usv*,*FL* is the voltage output of the spool dynamics inversion algorithm and *ysp*,*des* is the desired spool position obtained from Equation (19). The calculation of the required voltage input to the servovalve therefore implies calculating the time derivative of *ysp*,*des*. For that purpose, a fourth order backward differentiation scheme has been utilized here [28]:

$$\dot{y}\_n = \frac{3y\_{n-4} - 16y\_{n-3} + 36y\_{n-2} - 48y\_{n-1} + 25y\_n}{12\Delta t} \tag{21}$$

where *n* denotes the actual time step, *yn* the function evaluated at that time step and Δ*t* the time step which equals <sup>1</sup> <sup>×</sup> <sup>10</sup>−<sup>4</sup> s (see Table 1). Figure <sup>6</sup> illustrates the block diagram of the described feedback linearization scheme.

## *3.2. System Identification*

As mentioned before, the aim of the system identification module is twofold: (i) to work out an estimate of system's IF and (ii) to identify hydraulic system parameters which are required to carry out feedback linearization. These procedures are dealt with in next two subsubsections.

**Figure 6.** Feedback linearization scheme.

#### 3.2.1. Impedance Function Identification Procedure

The first step taken in finding a suitable approximation of the IF of the system has been to work out an estimate of its AF to later derive an appropriate IF. Two scenarios for IF identification are considered in what follows: (i) a noise-free environment and (ii) a more realistic situation in which noise contaminates servovalve voltage input and force and acceleration measurements. The former is presented for theoretical validation purposes, while the latter constitutes a robustness check of the IF identification procedure necessary to correctly assess the potential of the proposed methodology.

The followed procedure has been essentially the same for both cases and consists in feeding the servovalve with several (voltage) blocks of random stimuli and recording simultaneously the force on the table and the acceleration at the control point. These constitute, respectively, the input and the output of the shake table-SuT system (see Figure 5). Later on, the AF has been estimated making use of classical FRF estimation algorithms [29]. In the noise-free environment the *H*<sup>1</sup> estimator has been employed. With this approach the AF is expressed as *AF*(ω) = *GAF*(ω)/*GFF*(ω), where *GAF*(ω) is the cross-spectrum between force and acceleration and *GFF*(ω) is force autospectrum. For the noise-contaminated scenario, the *Hv* estimator, which minimizes the effect of the noise at both the input and the output of the system, has been employed.

The input to the servovalve has been selected to be a gaussian random waveform with a duration of 20 s, a flat frequency content between 0.1 Hz and 100 Hz and a maximum amplitude of 20 mV. This signal has been windowed with a unit square signal with a duty cycle of 50%. In this way, excitation is only effective during the first half of the block, allowing for system response (acceleration) to attenuate towards the end of the block, therefore minimizing leakage errors. Figure 7 shows signals obtained in one iteration of the identification stage in time and frequency domains.

**Figure 7.** Identification signals (one block): (**a**) Force on table and table acceleration in time domain; (**b**) Servovalve voltage, Force on table and table acceleration in frequency domain.

The AF approximation has been computed by linearly averaging the results obtained with sixteen input-output blocks. Figure 8 shows the achieved estimate for AF and its theoretical shape calculated analytically by transforming into frequency domain Equation (9) and performing the appropriate manipulations.

**Figure 8.** Identified and analytical AF (noise-free scenario): (**a**) Magnitude; (**b**) Phase.

As it was explained in Section 3.1, the feedback linearization scheme, in theory, allows for the imposition of an arbitrary time derivative of pressure force exerted on actuator's rod. Consequently, the IF sought must relate table acceleration to the time derivative of pressure force. This frequency function can be easily obtained by differentiating in frequency domain, without resorting to perform numerical derivatives on the desired pressure force obtained in time domain. Figures 9 and 10 show the FRF (acceleration over pressure force time derivative) and the IF (pressure force time derivative over acceleration) finally used by the drive calculation module, along with their theoretical value. The quality of the identified inverse model is quite acceptable within the complete frequency range of interest, except for the very low frequencies and in the neighborhood of the first modal frequency of the system under control, where small differences can be observed.

**Figure 9.** Identified and analytical FRF (noise-free scenario): (**a**) Magnitude; (**b**) Phase.

**Figure 10.** Identified and analytical IF (noise-free scenario): (**a**) Magnitude; (**b**) Phase.

In order to correctly assess the potential of the proposed methodology in a more realistic scenario, in what follows, the outcomes of an IF identification procedure in which servovalve voltage, force and acceleration signals have been influenced by electrical noise is presented.

The noise in measurements has been simulated by adding gaussian noise to servovalve input voltage and force and acceleration signals. The peak magnitude of the noise has been set to *un*,*peak* =0.01 V (see Section 2), which is an attainable value, when good industrial practices for low distance voltage signals wiring and shielding are observed. Later, the value of the noise affecting the physical quantities has been calculated by multiplying the electrical noise by each sensor's gain as explained in Section 2. Figures 11 and 12 show, in the presence of noise, the same information as Figures 9 and 10. Clearly, the quality of the estimates of FRF and IF decreases; nevertheless, identification error remains within reasonable limits and the obtained estimates are sufficiently good in the whole frequency range of interest. So as to better compare the FRF estimates, coherence functions associated to FRF identification, both in the noise-free and noise-contaminated cases are shown in Figure 13. The coherence function is defined as γ<sup>2</sup> = *GFA*(ω) 2 /[*GFF*(ω)*GAA*(ω)] and measures the degree of linear relationship between two signals. Even though coherence in the noise-contaminated case is evidently worse than its noise-free counterpart, its values remain quite close to unity in the frequency range of interest, consequently confirming the validity of the FRF and IF estimates in the presence of noise of a reasonable magnitude.

**Figure 11.** Identified and analytical FRF (noise-contaminated scenario): (**a**) Magnitude; (**b**) Phase.

**Figure 12.** Identified and analytical IF (noise-contaminated scenario): (**a**) Magnitude; (**b**) Phase.

**Figure 13.** Coherence functions comparison.

#### 3.2.2. Hydraulic Parameters Identification Procedure

The implementation of the feedback linearization scheme implies knowing accurate estimates of hydraulic parameters. In order to identify the sought values, and taking advance of the collection of data available from IF identification stage, a linear state-space model of the servoactuator has been identified. The inputs to this state-space model are, on the one hand, the voltage input to the servovalve, *usv*, and the force exerted on piston rod, *Fp*, by the shake table, on the other. The latter can be calculated by means of:

$$F\_p = \; m\_p \ddot{\mathbf{x}}\_t - (P\_1 - P\_2)A\_w + F\_{fr, \nu^\*} \tag{22}$$

which would correspond to the reading of a load cell installed between rod tip and shake table.

The state variables of the model have been selected as the velocity of the piston rod, which, as mentioned before, is identical to table velocity, . *xt*, the difference of pressures across chambers denoted by Δ*P* and the servovalve's main spool position, *ysp*. It has been assumed that Bulk moduli of each compartment are identical and are denoted by β. As mentioned in Section 2, the leakage flows between chambers and from chambers to bearings have been neglected. Accounting for these assumptions and linearizing Equations (1)–(5) around the mid-stroke operating point of the hydraulic cylinder, leads to the following analytical form of the state equations of the servoactuator system:

$$
\begin{Bmatrix} \dot{\bar{x}}\_{l} \\ \dot{\Delta P} \\ \dot{y}\_{sp} \end{Bmatrix} = \begin{bmatrix} -\mathbb{C}\_{\mathbb{P}}/m\_{\mathbb{P}} & A\_{\text{if}}/m\_{\mathbb{P}} & 0 \\ -2A\_{\text{if}}\mathbb{R}/\mathbb{r}\_{0} & 0 & 2\mathbb{C}\_{d}K\_{\mathbb{S}0}\mathbb{R}\sqrt{\mathbb{P}\_{\mathbb{S}}/\rho}/\mathbb{r}\_{0} \\ 0 & 0 & -1/\mathbb{r}\_{\mathbb{S}0} \end{bmatrix} \begin{Bmatrix} \dot{\bar{x}}\_{l} \\ \Delta P \\ \dot{y}\_{sp} \end{Bmatrix} + \begin{Bmatrix} 0 & 1/m\_{\mathbb{P}} \\ 0 & 0 \\ \mathbb{C}\_{\mathbb{S}p}/\mathbb{r}\_{sp} & 0 \end{Bmatrix} \begin{Bmatrix} \boldsymbol{u}\_{\mathbb{S}0} \\ \boldsymbol{f}\_{p} \end{Bmatrix} \tag{23}
$$

By means of a least squares procedure, the components of the matrices in Equation (23) have been identified. This process implies approximating the values of the rod velocity . *xt* and time derivatives of state variables Δ*P* and *ysp* (see Equation (21)).

Once the estimates of matrices components are available, it is possible to estimate directly the values of *Aw*, β/*v*0, *CdKsv* 2/ρ, τ*sp* and *Csp* used in the feedback linearization scheme and also the values of *mp* and *Cp*. In this work, it has been considered that the initial volumes of actuators chambers are known, and therefore, β can be estimated from β/*v*<sup>0</sup> value.

A check of the robustness against noise of the hydraulic parameters estimation process has been performed in the same way as for the IF estimation case. Table 2 shows the nominal and identified values of the hydraulic parameters and the identification relative error both for the noise-free and noise-contaminated identification cases. Despite the fact that parameters estimation quality decreases when a noisy environment is considered, the obtained values still represent with reasonable accuracy system actual parameters.


**Table 2.** Identified hydraulic parameters.

#### *3.3. Drive Calculation Module*

The drive calculation module computes the desired pressure force time derivative to be injected to the feedback linearization module. Firstly, the reference acceleration profile is transformed into frequency domain by means of an FFT process. Then, the drive signal is calculated multiplying the identified IF by the transformed acceleration profile. Finally, the result is transformed back into time domain by means of an IFFT process. Figure 14 schematizes the described process.

**Figure 14.** Drive calculation process.

#### *3.4. Three Variable Controller*

A feedback controller has been implemented, in parallel with the previously described architecture, with the aim to cope with the unavoidable errors occurring within AF and hydraulic parameters identification processes. TVC philosophy has been adopted so as to provide real-time simultaneous corrections to rod displacement, velocity and acceleration errors. The control law of the TVC is defined as follows:

$$u\_{TVC} = K\_d(d\_{ref} - \mathbf{x}\_p) + K\_v(v\_{ref} - \dot{\mathbf{x}}\_p) + K\_d(a\_{ref} - \ddot{\mathbf{x}}\_p)\_\prime \tag{24}$$

where *uTVC* is the control voltage output by TVC, *dre f* , *vre f* , *are f* are respectively the displacement, velocity and acceleration reference waveforms and *Kd*, *Kv* and *Ka* are the control gains for displacement, velocity and acceleration errors. As it can be noticed, the implementation of this controller

requires calculating, by integration, the reference displacements and velocities from the given reference acceleration.

#### **4. Numerical Simulations Results and Control Methods Comparison**

In this section, the numerical results obtained with the model described in Section 2 are presented. Section 4.1. shows and discusses simulation results for the new proposed control method in scenarios with and without noise present in measurements from sensors. Section 4.2 compares the performance of the suggested control procedure to that of the classical iterative control approach illustrating the main differences.

#### *4.1. Numerical Simulations Results*

The chosen acceleration reference in all the presented cases is a random gaussian waveform with a duration of 20 s and a flat frequency content between 1 Hz and 80 Hz (see Figure 17). A Hanning window has been applied to the reference profile to ensure null values at the block ends. The approximate peak displacement, velocity and acceleration values are 30 mm, 0.4 m/s and 25 m/s2. A fixed step solver and a time step of 1 <sup>×</sup> <sup>10</sup>−<sup>4</sup> s has been used for all the simulations in this section.

A first set of simulations has been carried out with the TVC feature disabled in a noise-free environment. As explained in Section 3.1, the feedback linearization module calculates an instantaneous spool position, which is attempted to be imposed on servovalve's main stage by means of a spool dynamics inversion algorithm (see Equation (20)). Figure 15 shows the reference and achieved servovalve spool position, both in time and frequency domain. It can be concluded that the tracking achieved by the dynamic inversion algorithm is accurate within the whole frequency range of interest.

**Figure 15.** Spool position tracking (noise-free scenario). TVC disabled: (**a**) Time domain; (**b**) Frequency domain.

The time derivative of pressure force on actuator's rod, synthesized by the drive calculation module explained in Section 3.3, and the actual one effectively imposed on servoactuator owing to the feedback linearization process, are shown in Figure 16. Figure 16a has been zoomed around the area where maximum error takes place. Tracking is acceptable in the entire frequency range of interest with larger errors at low frequencies and around the first modal frequency of the table-SuT system, due to the poorer IF estimate obtained at those frequency values (see Figure 10).

**Figure 16.** Pressure force time derivative tracking (noise-free scenario). TVC disabled: (**a**) Time domain; (**b**) Frequency domain.

Figure 17 shows the target acceleration together with the one attained in numerical simulations. A tracking error of approximately 0.55 m/s2 rms has been achieved. As in the previous case, and due to the same reasons stated there, reference acceleration tracking is reasonably satisfactory, except for the low frequencies and at the neighborhood of the first modal frequency of the system. Acceleration tracking error is shown in more detail in Figure 18. Despite the fact that tracking can be deemed acceptable, accumulation of errors within the low frequency region lead to increased velocity errors and displacement drifts which may hinder successful test execution due to limited servovalve flow rate capacity and actuator stroke. Therefore, it seems mandatory to enhance the control architecture with a parallel controller able to keep, simultaneously, all kinematic variables tracking errors within reasonable limits.

**Figure 17.** Shake table acceleration tracking (noise-free scenario). TVC disabled: (**a**) Time domain; (**b**) Frequency domain.

A second round of simulations has been carried out with the TVC feature enabled in a noise-free environment. Figures 19–21 show the same information as that offered in Figures 16–18. The employed values of TVC control gains have been *Kd* = 1, *Kv* = 0.5 and *Ka* = 0.25. Figures below show a drastic decrease in acceleration tracking error of from 0.55 m/s<sup>2</sup> rms in the previous case to 0.087 m/s2. This reduction especially marked in the low frequency range, down to 0.6 Hz, confirming the effectiveness of the TVC feedback controller.

**Figure 18.** Shake table acceleration error (noise-free scenario). TVC disabled: (**a**) Time domain; (**b**) Frequency domain.

**Figure 19.** Pressure force time derivative tracking (noise-free scenario). TVC Enabled: (**a**) Time domain; (**b**) Frequency domain.

**Figure 20.** Shake table acceleration tracking (noise-free scenario). TVC Enabled: (**a**) Time domain; (**b**) Frequency domain.

**Figure 21.** Shake table acceleration error (noise-free scenario). TVC Enabled: (**a**) Time domain; (**b**) Frequency domain.

A third set of simulations has been carried out to assess the suggested control method in a scenario where noise is present in sensors measurements. The electrical noise has been modelled as a gaussian waveform with a peak value of 10 mV according to considerations made in Section 3.2.1. Sources of noise of this magnitude have been added to all the transducers present in the model (displacement, acceleration, chamber pressures and spool position) and to servovalve input voltage. Noise in voltage has been multiplied by the appropriate gains to translate it into the physical quantities affecting the model (see Section 2). The drive estimate has been synthesized making use of the IF obtained in a noisy environment (see Section 3.2.1) and the hydraulic parameters utilized by feedback linearization scheme have been those estimated in presence of noise. The TVC feature has been enabled and the control gains used have been the same as in the previous case, that is, *Kd* = 1, *Kv* = 0.5 and *Ka* = 0.25.

Figures 22–24 show the same information as that in Figures 19–21. An overall tracking error of 0.403 m/s2 rms has been achieved. Tracking error has increased in the whole frequency range and specially around the higher frequency limit of the reference profile. The effect of noise is obviously more accused at lower target acceleration values due to the reduced signal to noise ratio at those sections. Nevertheless, despite the fact that electrical noise clearly affects negatively tracking quality, performance is still reasonably good and the stability of the system is maintained, therefore confirming the robustness of the proposed method when electrical noise of a reasonable magnitude contaminates sensors measurements.

**Figure 22.** Pressure force time derivative tracking (noise-contaminated scenario). TVC Enabled: (**a**) Time domain; (**b**) Frequency domain.

**Figure 23.** Shake table acceleration tracking (noise-contaminated scenario). TVC Enabled: (**a**) Time domain; (**b**) Frequency domain.

**Figure 24.** Shake table acceleration error (noise-contaminated scenario). TVC Enabled: (**a**) Time domain; (**b**) Frequency domain.

Finally, a fourth round of simulations has been conducted to explore the trend of system behavior in a noisy environment, when the values of the control gains of the TVC are increased while keeping the rest of the parameters unaltered. The values of the control gains employed have been *Kd* = 3, *Kv* = 2 and *Ka* = 1. Figures 25 and 26 show, respectively, acceleration and acceleration error both in time and frequency domain. With the employed control parameters, the influence of the noise in the system is remarkably reduced yielding similar results as in the second set of simulations. In particular, a tracking error of 0.1154 m/s2 rms has been attained.

**Figure 25.** Shake table acceleration tracking (noise-contaminated scenario). TVC Enabled with improved parameters: (**a**) Time domain; (**b**) Frequency domain.

**Figure 26.** Shake table acceleration error (noise-contaminated scenario). TVC Enabled with improved parameters: (**a**) Time domain; (**b**) Frequency domain.

#### *4.2. Comparison between Control Methods*

In this subsection, a comparison between the classical iterative control approach, which constitutes the industry standard for shake table testing, and the new proposed method is presented. A generic iterative scheme (see Figure 2), in a noise-free scenario, has been employed to obtain qualitative results representative of the classical method performance. First off, the FRF of the system has been identified making use of the *H*<sup>1</sup> estimator, according to traditional approach (acceleration over voltage). Then, it has been inverted to obtain the IF (voltage over acceleration). An estimate of the initial drive to be fed to the system has been calculated by means of: *d*0(ω) = *IF*(ω)*are f*(ω), where *d*0(ω) is the initial drive in frequency domain. The result has been transformed into time domain by means of an IFFT and has been injected into the servovalve. After the initial iteration, the drive signal would be updated in successive runs by an iterative scheme of the form: *dn*+1(ω) = *dn*(ω) + *kIF*(ω) !.. *xt <sup>n</sup>*(ω) − *are f*(ω) " , formulated in frequency domain, where *n* denotes iteration number and *k* is the correction gain. However, only the acceleration results obtained at the first iteration have been considered, with the aim of evaluating the control methods in comparable operation conditions. In these simulations, the TVC feature in the new control procedure has been enabled and the values of control gains used have been: *Kd* = 1, *Kv* = 0.5 and *Ka* = 0.25.

Figures 27 and 28 show acceleration response and tracking error for both methods. The first iteration of classical approach reaches a tracking error of 1.484 m/s<sup>2</sup> rms as opposed to the 0.087 m/s<sup>2</sup> rms featured by the new implementation. The proposed method shows much better behavior than the first iteration of the classical approach over the complete frequency range. Nevertheless, this difference in performances is likely to decrease if a certain number of control iterations were carried out. According to Figure 28a, the error of the classical approach increases with the magnitude of the target acceleration. This tracking error rise is caused by the fact that this method relies on a linearization of a non-linear system around an operating point, which may no longer be valid when the target acceleration profile implies reaching large values of forces and displacements. In opposition, the new suggested procedure performs well even at high accelerations, in part, due to the fact that the implemented feedback linearization scheme excludes the non-linearities associated to hydraulic system from the control loop.

**Figure 27.** Comparison between classical and new control methods. Acceleration tracking: (**a**) Time domain; (**b**) Frequency domain.

**Figure 28.** Comparison between classical and new control methods. Acceleration error: (**a**) Time domain; (**b**) Frequency domain.

#### **5. Conclusions**

This paper presents a novel mixed time-frequency acceleration control method for shake table systems. The suggested procedure includes identifying an FRF-IF pair which relates table acceleration to the time derivative of the pressure force acting on servoactuator's piston rod prior to the test. This approach allows for the calculation of a time variation of pressure force waveform which can

be directly imposed by means of a feedback linearization scheme, which approximately cancels out non-linearities associated hydraulic actuation system leaving only those related to SuT present in the control loop. Modeling errors, feedback linearization imperfections and external perturbations are dealt with by a TVC implemented in parallel with the previously outlined architecture. Consequently, this method includes features belonging to time and frequency domain methods usually employed in shake table control systems.

The potential effectiveness of the methodology was assessed by means of numerical simulations carried over a model of the shake table loaded with a two stories shear building. Four groups of numerical simulations were performed:


Results corresponding to the first group show quite acceptable acceleration tracking; however, errors at low frequencies may lead to undesired table drifts. When the TVC feature is enabled and electrical noise is not considered (group 2), tracking errors are drastically reduced, leading to almost-perfect acceleration tracking with a low control burden placed on the TVC controller. Noise affects negatively the performance in the whole frequency range, as it was demonstrated by the third group of simulations; however, tracking error remains within acceptable limits and the stability of the system is preserved, thus confirming the robustness of the proposed control procedure when electrical noise of a reasonable magnitude contaminates sensors measurements. Finally, the fourth group of simulations demonstrates that, by a proper tuning of the TVC control parameters, almost-perfect tracking is possible even in the presence of noise. The performance of the new proposed method is better than that of the classical iterative approaches when these operate on a single iteration basis.

The presented method thus appears quite promising for its implementation in real systems and features the following advantages over traditional iterative methods: (i) no iterations are required in test execution stage, (ii) non-linearities associated to hydraulic actuation are excluded from the control loop, improving tracking characteristics and (iii) the method is less sensitive to uncertainties in IF identification than the traditional control approaches due to the parallel TVC feature.

The proposed methodology requires, however, measuring or estimating rod displacement, velocity and acceleration, pressures at both actuator's chambers, pressures at supply and return ports of actuator manifolds and position of servovalve's main stage spool position. Therefore, its implementation implies increased instrumentation needs with respect to that used in traditional control methods.

Current works are focused on the implementation of the proposed architecture in an actual shake table system, paying special attention to the following points:


**Author Contributions:** Conceptualization, J.R.S., J.H.G.-P. and I.M.D.; methodology, J.R.S.; software, J.R.S.; validation, J.R.S.; formal analysis, J.R.S.; investigation, J.R.S. and I.M.D.; resources, J.R.S.; data curation, J.R.S.; writing—original draft preparation, J.R.S.; writing—review and editing, J.R.S., J.H.G.-P. and I.M.D.; visualization, J.R.S.; supervision, J.H.G.-P. and I.M.D.; project administration, J.H.G.-P. and I.M.D.; funding acquisition, J.H.G.-P. and I.M.D. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Spanish Ministry of Science, Innovation and Universities through the project SEED-SD (RTI2018-099639-B-I00).

**Acknowledgments:** The authors would like to express their gratitude to Vzero Engineering Solutions, S.L. for providing several pictures.

**Conflicts of Interest:** The authors declare no conflict of interest.
