Acceleration Harmonic Estimation for Hydraulic Servo Shaking Table by Using Simulated Annealing Algorithm

Abstract

In order to study the anti-shock performance of equipment under different working conditions, a hydraulic servo shaking table is used to replicate the desired motion with high fidelity. However, the sinusoidal acceleration response waveform is not a pure sinusoid due to the inherent nonlinearities within the system. The generating higher harmonic deteriorates the control performance and leads to system instability. To suppress the harmonic distortion and accurately estimate harmonic information, the harmonic estimation scheme based on simulated annealing algorithm is proposed. The sum of error square between actual value and estimated value is defined as the objective function. The amplitude and phase of each harmonic can be directly extracted when the objective function is minimized. Simulation and experimental results indicate that the proposed algorithm has good convergence performance and high estimation precision.

1. Introduction

A shaking table is usually used to replicate various vibration conditions that occur in civil engineering, automation industry, and aseismatic performance tests [1,2,3,4,5]. According to the drive mode, a shaking table falls into categories of mechanical shaking tables, electrical shaking tables, and hydraulic servo shaking tables. Hydraulic servo shaking tables are widely used in the engineering field because of their advantages: higher force-to-weight ratio, higher response speed, higher load stiffness, and control precision [6]. However, they have many drawbacks, such as a dead zone in the servo valve, friction and backlash in joints, and the dynamic characteristic in hydraulic actuator [7]. Due to the above nonlinear factor, the sinusoid acceleration response shows amplitude attenuation and phase delay, which causes waveform distortion. The harmonic distortion deteriorates the control performance of system and leads to system instability [8]. In order to suppress harmonic distortion and precisely reproduce the desired signal, it is necessary to obtain the amplitude and phase of the each harmonic component. Furthermore, the estimated information can be used to compensate for the harmonic distortion [9].

In past decades, various estimation methods have been developed in order to accurately estimate harmonic information. The most widely used method is the fast Fourier transform (FFT). Nevertheless, this method has some problems, namely aliasing, leakage, and picket-fence effect, which could result in harmonic estimation inaccuracy [10]. Some recursive algorithms, like the least mean square (LMS) and recursive least square (RLS), are also proposed to solve harmonic estimation problem [11,12,13]. Ray et al. developed the ensemble Kalman filter to estimate the amplitude and phase of distorted power system signal [14,15]. Barros presented the wavelet-packet transform to estimate the harmonic groups in voltage and current waveforms [16]. Sahoo et al. studied a digital algorithm for harmonic identification based on complex H∞ filter [17]. Soliman et al. utilized a new algorithm based on fuzzy linear regression to estimate the frequency and harmonic component of the voltage signal in a power network [18]. With the development of intelligent optimization theory, it offers the possibility of estimating the harmonic component with modern optimal algorithm [19]. Almaita and Asumadu applied the radial basis function (RBF) neural network to duly estimate harmonic components of converter waveforms [20]. Wang used a Hopfield neural network to simultaneously determine the supply frequency variation, the amplitude variation and phase variation of fundamentals, as well as each harmonic [21]. Ray et al. proposed particle swarm optimization (PSO) algorithm to identify current harmonic associated with fundamental current from measurement [22]. Ji et al. adopted an adaptive bacterial swarming (ABS) algorithm to dynamically estimate the frequencies and phases of the voltage or current harmonics in power grids [23].

The difficulty in harmonic estimation is due to the fact that harmonic generating is dynamic and nonlinear in nature. Thus, fast and accurate harmonic estimation approaches are required. However, the aforementioned methods can only be implemented offline and have several limitations on convergence and accuracy. In this paper, the simulated annealing optimization algorithm is used for hydraulic servo shaking table acceleration harmonic estimation. Simulation and experimental results indicated that the proposed algorithm can accurately estimate the amplitude and phase of each harmonic. Furthermore, the estimation algorithm has better real time performance and convergence performance when compared with conventional algorithm, such as RLS.

2. Hydraulic Servo Shaking Table

Figure 1 shows the hydraulic servo shaking table developed by Harbin Engineering University, which includes three subsystems: mechanical system, hydraulic system, and control system. The mechanical system consist of a platform and a symmetric cylinder controlled by a two-stage symmetric servo valve. The hydraulic system is used to provide hydraulic energy and control the pressure and flows of servo actuator. The control system consists of a host computer and a target computer. The host computer serves as user interface and allows users to generate a command signal, and obtain current experiment data. The real time control algorithm and serial communication with the programmable logic controller are operated on the target computer. Some main parameters of the shaking table are listed in

Table 1. The main parameters of shaking table.

Component	Parameter
Piston diameter	40 mm
Rod diameter	35 mm
Stroke	25 mm
Supply pressure	8 Mpa
Frequency range	0~50 Hz
Maximum velocity	1.4 m/s
Maximum acceleration	10 m/s2

.

2.1. Dynamic Model

The power mechanism of the hydraulic servo shaking table includes servo valve, hydraulic cylinder, and load, which can be simplified as symmetric hydraulic cylinder controlled by four-way valve. The schematic diagram of the power mechanism is shown in Figure 2, where $p_s$ and $p_b$ are the supply/return pressure, $q_i$ and $q_o$ are the input/output oil flow of the cylinder, $p_i$ and $p_o$ are the input/output pressure, $V_i$ and $V_o$ are the input/output oil volume, m is the mass, $x_v$ is the valve position, B is the piston’s viscous damping coefficient, y is the load position, A is the effective area of the cylinder. The physical features of shaking table’s hydraulic system can be described by the following three equations.

The linearized flow equation of the servo valve can be expressed as

$$q_L=K_q{x}_v-K_c{p}_L$$

(1)

where $q_L=(q_i+q_o)/2$ is load flow, $K_q=\partial q_L/\partial x_v$ is flow gain, $p_L=p_i-p_o$ is pressure burden, $K_c=\partial q_L/\partial p_L$ is flow-pressure coefficient.

The flow continuity equation of the cylinder is

$$q_L=A\frac{dy}{dt}+C_{tp}{p}_L+\frac{V_t}{4{\beta }_e}·\frac{d{p}_L}{dt}$$

(2)

where $C_{tp}$ is leakage coefficient, $V_t=V_i+V_o$ is the total volume of the actuator’s two chambers, ${\beta }_e$ is the effective bulk modulus.

The force balance equation is

$$F_g=A\left(p_i-p_o\right)=A{p}_L=m\frac{d^{2}y}{d{t}^2}+B\frac{dy}{dt}+F_L$$

(3)

Performing Laplace transform from Equation (1) to Equation (3), the transfer function $G_p(s)$ from the valve displacement $x_v$ to the hydraulic cylinder displacement y can be expressed as

$$G_p(s)=\frac{Y}{X_v}=\frac{K_q/A}{s\left(\frac{s^2}{{\omega }_h^2}+\frac{2{\zeta }_h}{{\omega }_h}s+1\right)}$$

(4)

where ${\omega }_h=\sqrt{\frac{4{\beta }_e{A}^2}{m{V}_t}}$ is load inherent resonant frequency, ${\zeta }_h=\frac{K_c}{A}\sqrt{\frac{{\beta }_{e}m}{V_t}}+\frac{B_c}{4A}\sqrt{\frac{V_t}{{\beta }_{e}m}}$ is load damping ratio.

2.2. Control Principle

Figure 3 shows the schematic diagram of the hydraulic servo shaking table. The three variable controller (TVC) is a commonly used controller for shaking table [24]. The three variable of TVC represent displacement, velocity, and acceleration of the shaking table. The TVC feedforward gain ($K_{dr}$, $K_{vr}$, $K_{ar}$) is used to improve system stability and the feedback gain ($K_{df}$, $K_{vf}$, $K_{af}$) is used to extend frequency bandwidth and reduce tracking error. The velocity feedback signal is usually synthesized from the measured displacement and acceleration, as it is difficult to measure the displacement, velocity, and acceleration simultaneously. The velocity in the high/low frequency range is obtained from the displacement/acceleration using a high pass filter (HPF) or low pass filter (LPF), respectively. The measured displacement, velocity, and acceleration are acquired by the ACL-8316 A/D port, transforming analogue signal into digital signal. The input signal is obtained by input filter. The TVC compares the feedback signal ($d_f$, $v_f$, $a_f$) with the given signal ($d_r$, $v_r$, $a_r$) to generate the control signal e, which is sent to hydraulic actuator by means of ACL-8316 D/A port that transforms digital signal into analogue signal. Hydraulic power is supplied by the hydraulic pump to drive the hydraulic cylinder, causing the desired movement of hydraulic servo shaking table.

3. Simulated Annealing Algorithm

The simulated annealing is a population based stochastic algorithm, which was first introduced by Metropolis in 1953 [25]. Since then, it has been widely used to solve combination optimization problem. The algorithm stimulates the process of heating up a solid to a high temperature followed by slow cooling [26,27,28,29,30,31]. As the temperature goes up, the particles of solid change into disordered state, causing an increase in internal energy. However, the internal energy of solid decreases with a drop in temperature, and the particles gradually back in order. If the cooling rate is cooled slowly enough for the solid, the particles can reach thermal equilibrium at each temperature [32]. The minimum energy is obtained by reaching the ground state of room temperature. According to the Metropolis guidelines, the probability of particles tend to become equilibrium state $\exp (-\Delta f/kT)$ at a given temperature T, where k denotes the Boltzmann constant. Assuming a solid in current state m₀ with energy level f(m) and the next state m′ with energy level f(m′), if the difference between the two energy levels is less than or equal to zero, the new state m is accepted. Otherwise, if the difference $\Delta f=f(m^{\prime })-f(m)$ is greater than zero, the new state is accepted with Metropolis criteria. The acceptance probability is compared to a number $P_r\in \left[0,1\right]$ generated randomly and m is accepted whenever $P>P_r$. The initial temperature T₀ is supposed to be high enough to allow acceptance of any new state in the first step. In each step, the procedure generates a fixed number of neighborhood solutions and evaluates them using the current temperature value $T_s=\alpha T_0$, where $\alpha (0<\alpha <1)$ is an attenuation factor [33,34]. The proposed simulated annealing algorithm used in this paper can be summarized as Figure 4.

4. Harmonic Estimation Scheme

The distorted signal can be considered as the sum of fundamental response and harmonic component. Each harmonic has its own amplitude and phase, and a frequency that is integral multiple of the fundamental response. Therefore, the sinusoidal acceleration response of the shaking table at time t can be expressed by

$$a(t)={\displaystyle \sum _{n=1}^N{A}_n}\sin (n\omega t+{\varphi }_n)$$

(5)

where N represents the order of the harmonic; ω is the fundamental frequency; A_n and φ_n are the amplitude and phase of nth harmonic, respectively.

The discrete time version of Equation (5) can be represented as

$$a(k)={\displaystyle \sum _{n=1}^N{A}_n}\sin (n\omega k+{\varphi }_n)$$

(6)

For estimating amplitude and phase, Equation (6) can be rewritten as

$$a(k)={\displaystyle \sum _{n=1}^N(A_n}\cos {\varphi }_n\sin n\omega k+A_n\sin {\varphi }_n\cos n\omega k)$$

(7)

Define $x_{n1}=A_n\cos {\varphi }_n$, $x_{n2}=A_n\sin {\varphi }_n$. Then, Equation (7) can be written as

$$a(k)={\displaystyle \sum _{n=1}^N(x_{n1}}\sin n\omega k+x_{n2}\cos n\omega k)$$

(8)

Sinusoidal signal in parametric form becomes

$$a(k)=\mathbf{h}(k)\mathbf{x}(k)$$

(9)

where $\mathbf{h}(k)={[\sin \omega k,\cos \omega k,\sin 2\omega k,\cos 2\omega k,\cdots ,\sin N\omega k,\cos N\omega k]}^{\mathrm{T}}$ is the input vector and $\mathbf{x}(k)={[x_{11},x_{12},x_{21},x_{22},\cdots x_{N1},x_{N2}]}^{\mathrm{T}}$ is the weight vector.

The estimation model for the algorithm is

$$\widehat{a}(k)=\mathbf{h}(k)\widehat{\mathbf{x}}(k)$$

(10)

where $\widehat{\mathbf{x}}(k)$ is the estimation weight vector.

The estimate value $\widehat{\mathbf{x}}(k)$ for the requested parameters $\widehat{\mathbf{x}}(k)$ can be obtained by minimizing the objective function

$$f(\mathbf{x})={\displaystyle \sum _{n=1}^N{[e(k)]}^2}={\displaystyle \sum _{n=1}^N{\left[a(k)-\widehat{a}(k)\right]}^2}$$

(11)

The harmonic estimation scheme based on simulated annealing algorithm can be illustrated with Figure 5. The sinusoidal acceleration response is used as desired signal, and the reference harmonic generates the input vector. The error between the actual value and the estimated value is optimized using simulated annealing algorithm. The problem is formulated as an optimization problem, where the goal is to minimize the objective function. After updating the weight vector, the amplitude and phase of the fundamental and nth harmonic are derived as follows

$$\{\begin{array}{l}{A}_n=\sqrt{x_{n1}^2+x_{n2}^2}\\ {\varphi }_n=\tan (x_{n2}/x_{n1})\end{array}$$

(12)

Moreover, fundamental as well as each harmonic can be directly extracted from the final updated weights

$${\widehat{a}}_n(k)=x_{n1}\sin n\omega k+x_{n2}\cos n\omega k$$

(13)

5. Simulation and Results

In order to verify the effectiveness of the developed acceleration harmonic estimation algorithm, the simulation is firstly implemented. The simulation signal contains six harmonic. It is

$$\begin{array}{cc}\hfill a(t)=& 10\sin (10\pi t)+8\sin (20\pi \mathrm{t}-1.2)+6\sin (30\pi \mathrm{t}+0.55)+4\sin (40\pi t-0.8)\hfill \\ & +2\sin (50\pi t+1.4)+1\sin (60\pi t+1)\hfill \end{array}$$

Figure 6 and Figure 7 are the estimated amplitudes and phases, respectively. It is seen that the estimated values converge to their corresponding original values after about 0.4 s. Figure 8 shows the error between the given signal $a(t)$ and the estimated signal $\widehat{a}(t)$ along with time. The estimation error is asymptotically converged to zero within 0.4 s, which means that the given signal is estimated precisely. The amplitude and phase of each estimated harmonic when they in stable states are listed in

Table 2. The estimated value compared with given value.

Harmonic Order	Given Value		Estimated Value
	Amplitude (m/s2)	Phase (rad)	Amplitude (m/s2)	Phase (rad)
Fundamental response	10	0	9.999816	−0.000004
Second harmonic	8	−1.2	8.000417	−1.19999
Third harmonic	6	0.55	5.999853	0.550035
Fourth harmonic	4	−0.8	4.000689	−0.79981
Fifth harmonic	2	1.4	1.999284	1.399997
Sixth harmonic	1	1	0.999776	0.999728

with its given value. At the same time, it is clear that the estimated amplitude and phases are very close to its given value for each harmonic. Thus the developed acceleration harmonic estimation scheme can estimate harmonic information efficiently.

6. Experiment and Results

When the acceleration excitation is $4\sin (2\pi \times 5t){\mathrm{m}/\mathrm{s}}^2$, the excitation signal and acceleration response in the time domain are shown in Figure 9. It can be seen that the acceleration response is seriously distorted, due to the effects of nonlinearities in the hydraulic servo shaking table. Since higher harmonics are less dominant in all harmonics, the harmonic estimation is usually up to sixth harmonic in practical applications. The amplitude frequency diagram computed by FFT in Figure 10 indicates that the fundamental response is at 5 Hz, the second harmonic is at 10 Hz, the third harmonic is at 15 Hz, and so on. The total harmonic distortion (THD) is a very important criterion to evaluate its waveform performance and it is given by

$$\mathrm{THD}=\frac{\sqrt{A_2^2+A_3^2+\cdots A_N^2}}{A_1}\times 100\%$$

(14)

where A₁ is the amplitude of the fundamental, A₂ the amplitude of the second harmonic, A₃ the amplitude of the third harmonic, and so on.

The result of THD based on FFT is calculated off line as shown in

Table 3. Total harmonic distortion analysis result.

THD	Harmonic Amplitude (m/s2)
22.22%	A1	A2	A3	A4	A5	A6
	3.9830	0.4162	0.3510	0.1971	0.5824	0.3306

. It can be seen that the fifth harmonic plays a vital role in all harmonics, and the fourth harmonic is the least harmonic. The THD is 22.22%, which is a large value.

Figure 11 and Figure 12 represent the amplitude and phase estimation results respectively using simulated annealing and RLS algorithm of the given acceleration response. As can be seen from the figures, RLS based estimation method has large estimation error and oscillatory in the initial stage of estimation process, though it converges to the stable value. From the perspective of harmonic estimation, the estimated algorithm by simulated annealing is better than RLS in real time performance and convergence performance. In addition, each harmonic can be directly extracted from the distortion signal. Its results are shown in Figure 13, from which it can be seen that there are large fluctuations at the beginning of the estimation, but converge quickly to its steady states within 0.5 s.

7. Conclusions

Due to the inherent nonlinearities within the hydraulic servo shaking table, the acceleration response produced is a time varying signal. Thus, fast and accurate harmonic estimation algorithm is required. In this paper, the simulated annealing optimization algorithm is performed on the acceleration sinusoidal response replication. Simulation and experimental results show that the harmonic estimation scheme can estimate the amplitude and phase of each harmonic and reconstruct individual harmonics at the same time. As an intelligent optimization technique, the proposed algorithm does not require a priori knowledge of the system.

Compared with conventional RLS algorithms, experimental results show that the proposed simulated annealing algorithm is better in convergence performance and real time performance. This is important for harmonic estimation and harmonic suppression. In addition, the proposed technique can be extended to many areas: online estimation, optimization, and control problems.