# Acceleration Harmonic Estimation in a Hydraulic Shaking Table Using Water Cycle Algorithm

Authors: Jianjun Yao; Zhenshuai Wan; Yu Fu
Journal: Shock and Vibration
 (2018)

---

## Abstract

Hydraulic shaking table is mainly used to stimulate the desired vibration environment and evaluate the shock resistance of structure. However, due to the inherent nonlinearities of the hydraulic shaking table, the acceleration response displays amplitude attenuation and phase delay for sinusoidal excitation signal. The distorted response degrades the control performance and even leads to an increase in system instability. In this paper, the water cycle algorithm (WCA), a recently developed metaheuristic method, is developed to estimate the harmonic information such as amplitude and phase. The basic idea of the proposed algorithm is inspired from nature and based on the observation of water cycle process and how rivers and streams flow to the sea in the nature world. The estimation process based on WCA is sequentially updating the weight vector of the signal. The amplitude and phase of fundamental as well as each harmonic can be achieved when the objective function is minimized. Simulation and experimental results demonstrate that the proposed harmonic estimation algorithm has good real-time performance, fast convergence, and excellent accuracy.


---

## Body

## 1. Introduction

The shaking table is an essential experimental facility for the performance assessment of structures when they are subjected to ground motions [1, 2]. The drive types of the shaking table include mechanical drive, electrical drive, and hydraulic drive. Compared with other drive modes, hydraulic drive has such advantages as high power to weight ratio, load stiffness, response speed, and control precision [3]. What is more, the hydraulic shaking table is able to generate a large force at high speed within a relatively wide frequency range [4–6]. As a result, this shaking table has been widely applied in civil engineering, earthquake engineering, and industrial applications for vibration tests, where the desired acceleration is required to be reproduced exactly [7]. However, the acceleration response is hard to match the desired motion. One of the reasons is that hydraulic shaking table has many inherent nonlinearities, such as dead zone in servo valve, friction, and backlash in joints as well as the dynamic characteristic of hydraulic actuator [8]. This can cause harmonic distortion for sinusoidal vibration tests, resulting in lower control performance [9].To effectively cancel the harmonic distortion and improve the acceleration reproducibility of the hydraulic shaking table, it is necessary to know the harmonic information. The most used way to analyze harmonic is the Fast Fourier Transform (FFT) method, which transforms samples of time domain into the frequency domain, but the phenomena of aliasing, leakage, and picket fence effects may lead to inaccurate harmonic estimation [10]. Besides, the time delay of sampling and computation of Fourier coefficients make it difficult for real-time implementation [11]. To overcome the shortage of employing FFT method for analyzing harmonic, other advanced techniques have been developed. Yao et al. presented normalized least-mean-square (LMS) algorithm to identify acceleration harmonic of the shaking table. For the LMS-based harmonic identification method, the harmonic estimation accuracy, especially for the harmonic phase, descended along with the harmonic order [12]. The Kalman filter based method firstly modeled a state space model of the acceleration response which contains higher harmonics [13, 14]. By estimating the state of the model, the harmonic amplitudes and phases can be computed. The estimation accuracy is related to the state space model. Different acceleration models can lead to different estimation accuracy.In recent times, several artificial intelligence algorithms have been reported for harmonic estimation. Mori et al. used artificial neural network (ANN) based on back propagation learning technique to estimate voltage harmonic in power system [15]. Dash et al. presented an adaptive neural network method based adaptive perceptron to estimate the harmonic components of a power system. Bettayeb and Qidwai proposed genetic algorithm (GA) to estimate the amplitudes and phases of various frequency components present in the power system [16]. Ji et al. adopted an adaptive bacterial swarming algorithm to estimate the frequencies and phases of the harmonics of voltage or current waveform for power network [17]. Lu et al. utilized the particle swarm optimizer (PSO) with passive congregation to estimate the phases of the harmonic and the least square method that is used to estimate the amplitude [18]. Ray and Subudhi proposed bacterial foraging optimization (BFO) and recursive least square (RLS) hybrid technique to estimate the fundamental as well as harmonic component present in power system voltage and current waveform [19].From the above discussion, it can be seen that most of the current literature on harmonic estimation is focused on electric system, but it provides a way for harmonic estimation of hydraulic shaking table. The hydraulic shaking table needs to deal with mechanic signal, hydraulic signal, and electronic signal simultaneously, so the real time is to be highly required. However, the aforementioned harmonic estimation methods either converge slowly, like GA algorithm, or need prior knowledge to detect the amplitude and phase of each harmonic, like Kalman filter. In order to estimate the harmonic information fast and accurately, the WCA is proposed to estimate the amplitude and phase of fundamental as well as each harmonic. The WCA stimulates nature phenomenon: streams flow into rivers and rivers flow into the sea through several processes such as evaporation, transpiration, precipitation, condensation, and surface run-off [20]. The above phenomenon can make the proposed algorithm find a global optimal solution or a near-optimal solution via effective exploration and exploitation. What is more, the WCA just requires few insensitive user parameters, which means that it can address a wide range of optimization problems using the fixed user defined parameters.


## 2. The Hydraulic Shaking Table

The hydraulic shaking table shown in Figure1 is a 3 m × 3 m six-degree-of-freedom shaking table. It is capable of simulating vibration events with acceleration of 5.6 g in vertical and 3.7 g in horizontal, respectively. Specimens up to 15 tons can be tested on the table at reduced acceleration. Four horizontal actuators rated at 70 KN each and 4 vertical actuators at 70 KN each provide the force input to the table. The supply pressure is 20.5 MPa and the system working frequency range is up to 100 Hz.Figure 1
The hydraulic shaking table.Three-variable controller (TVC) is used as basic controller because of hydraulic system poor damping ratio and lower natural frequency. The three variables in the TVC are displacement, velocity, and acceleration of the hydraulic shaking table. The reference (measured) displacement, velocity, and acceleration aredr,vr,ar(df,vf,af), respectively. The TVC consists of two parts: the feedforward part with three adjustable gains Kdr, Kvr, and Kar and the feedback forward part with three adjustable gains Kdf, Kvf, and Kaf. The feedforward part is used to improve the stability performance of the hydraulic shaking table and the feedback part is used to extend the frequency bandwidth and improve the tracking performance [21]. The TVC is designed based on the displacement control, but the hydraulic shaking table is under acceleration control [22]. Thus, a reference signal generator is designed to convert acceleration input to reference state. The transfer function of the reference signal generator can be expressed as(1)Gs=drR=Kus2+v0s+d0=Kud0s2/ω02+2ξ0/ω0s+1,where d0=ω02, v0=2ξ0ω0, dr is the reference displacement, and R is the reference acceleration. Ku is the acceleration gain, ξ0 is the damping ratio of the acceleration control, and ω0 is the original frequency of the acceleration.In general, it is difficult to simultaneously measure the displacement, velocity, and acceleration of the hydraulic shaking table. Thus, the velocity is synthesized from the measured displacement of the hydraulic actuator using a low-pass filter or the measured acceleration of the shaking table using a high-pass filter. The block diagram of TVC is shown in Figure2.Figure 2
Block diagram of the TVC.The transfer function of feedforward part can be expressed as(2)G1s=Kdr+Kvrs+Kars2.The transfer function of feedback part can be expressed as(3)G2s=Kdf+Kvfs+Kafs2.The control scheme of the hydraulic shaking table is shown in Figure3Figure 3
Block diagram of the control scheme.


## 3. Water Cycle Algorithm

The WCA is a novel metaheuristic optimization method proposed by Eskandar et al. in 2012 [23]. The inspiration of the WCA is based on the observation of water cycle process and how rivers and streams flow downhill towards the sea in the nature world. Figure 4 is the simplified schematic diagram of the water cycle. The water in lakes and rivers is evaporated into the atmosphere and leads to the formation of clouds that condense in the colder air above. During photosynthesis, the same process happens in the plant of transpiration. Then, the water is back to the earth in the form of rain or precipitation [24, 25].Figure 4
The process of water cycle.As with other metaheuristic algorithms, the proposed method begins with an initial population, that is, the streams [26]. First, suppose that we have rain or precipitation. The best individual (best stream) is chosen as a sea. Then, some good streams are chosen as rivers and the rest of steams flow into the rivers and sea. Each river absorbs water from the streams according to their magnitude of flow. In fact, the amount of water in a stream entering a river and/or sea varies from other streams. Moreover, rivers directly flow to the sea at the most downhill location [27].To solve an optimization problem using the proposed algorithm, it is necessary that the values of problem variables must be formed as an array. In aDvar dimensional optimization problem, a stream is an array of 1×Dvar. This array is defined as follows:(4)Stream=x1,x2,…,xD.To start the optimization algorithm, a candidate representing a matrix of streams of sizeNpop×Dvar is generated. The matrix generated randomly is given as(5)Totalpopulation=SeaRiver1River2⋮StreamNsr+1StreamNsr+2⋮StreamNpop=x11x21⋯xDvar1x12x22⋯xDvar2⋮⋮⋮⋮x1Npopx2Npop⋯xDvarNpop,where Npop and Dvar are represented by the number of streams (initial population) and design variables, respectively.The stream cost could be determined by calculating the function of cost (C) as(6)Ci=Costi=∫x1i,x2i,…,xNvarii=1,2,…,Npop.First, Npop streams are created. A number of Nsr are chosen as the sea and rivers from the best individuals (minimum values). The stream with the least value among the rest is taken as a sea. Actually, Nsr represents the total number of rivers and a single sea as shown in (7). The remainder of the population is determined by using (8).(7)Nsr=NumberofRivers+1︸Sea,(8)Nstream=Npop-Nsr.The following equation is used to assign streams into the sea or the rivers regarding the strength of the flow:(9)NSn=roundCostn∑i=1NsrCosti×NStreams,n=1,2,…,Nsr,where NSn represents the number of streams that flow into certain sea or rivers.The streams are linked with each other to generate new rivers, where some of the streams may directly flow into the sea. All the streams and rivers ultimately end in the sea. This idea can also be applied on rivers that flow into the sea, so the new position for the rivers and streams can be given as [28](10)XStreami+1=Xstreami+rand×2×XSeai-XStreami,XStreami+1=Xstreami+rand×2×XRiveri-XStreami,XRiveri+1=XRiveri+rand×2×XSeai-XRiveri,where i is iteration times and rand is a uniformly distributed random number between 0 and 1. If the solution given by a stream is better than its connecting river, the positions of river and stream are exchanged (stream becomes river and vice versa). Such exchange can similarly happen for rivers and sea.If the distance between the sea and the river is less thandmax (closer to zero), it signifies that the river arrived at or linked with the sea. The evaporation process is taken into consideration in this situation and as can be observed in nature, after ample evaporation has taken place, it will begin to rain or precipitate. A large dmax value will lower the search but a small value will encourage an intensification of the search close to the sea. As such, the intensity of the search close to the sea is controlled by dmax. The value of dmax adapts accordingly and decreases as(11)dmaxi+1=dmaxi-dmaximaxiteration.After satisfying the evaporation process, the raining process is applied. In the raining process, the new raindrops form streams in the different locations. For specifying the new locations of the newly formed streams, the following equation is used:(12)XStreamnew=LB+rand×UB-LB,where LB and UB are the lower and upper bounds, respectively.Again, the best newly formed stream is considered as a river flowing to the sea. The rest of new streams flow to the rivers or may directly flow to the sea. In order to enhance the convergence rate and computational performance of the algorithm for engineering problems, (13) is used only for the streams which directly flow to the sea.(13)XStreamnew=XSea+μ×randn1,Nvar,where μ is a coefficient which shows the range of searching region near the sea and randn is the normally distributed random number. The larger value for μ increases the possibility of exiting from feasible region. On the other hand, the smaller value for μ leads the algorithm to search in smaller region near the sea. A suitable value for μ is set to 0.1.The procedure for the proposed WCA is shown in Figure5.Figure 5
The flowchart of WCA.


## 4. Harmonic Estimation Problem

A distorted signal can be modeled as the sum of higher order harmonic of unknown amplitude and phase and can be expressed as(14)at=∑n=1NAnsin⁡ωnt+φn,where n=1,2,…,N represents the order of the harmonic; An, φn, and ωn are the amplitude, phase, and angular frequency of the nth harmonic, respectively. The angular frequency can be expressed as ωn=2πnf1, where f1 is the fundamental frequency.The discrete time version of (14) can be represented as(15)ak=∑n=1NAnsin⁡ωnkTs+φn,where Ts is sampling period.For the amplitude and phase estimation, (15) can be rewritten as(16)ak=∑n=1NAncos⁡φnsin⁡ωnkTs+Ansin⁡φncos⁡ωnkTs.Signal in parametric form becomes(17)ak=Hkθk.According to the structure of the WCA estimation scheme, the input vector and weight vector are defined as follows:(18)Hk=sin⁡ω1kTs,cos⁡ω1kTs,sin⁡ω2kTs,cos⁡ω2kTs,…,sin⁡ωNkTs,cos⁡ωNkTsTθk=A1cos⁡φ1,A1sin⁡φ1,A2cos⁡φ2,A2sin⁡φ2,…,ANcos⁡φN,ANsin⁡φNT.Letxn1=Ancos⁡φn, xn2=Ansin⁡φn; weight vector can be rewritten as θ(k)=x11,x12,x21,x22,…,xN1,xN2T.The sum of error square between the actual acceleration and the estimated acceleration is defined as the objective function, which is used for optimizing the weight vector and can be written as(19)J=∑k=1Nse2k=∑k=1Nsak-a^k,where a^(k) represents the estimated output of the harmonic signal, which is optimized using WCA algorithm.The WCA based harmonic estimation scheme can be illustrated with Figure6. After updating the weight vector using the WCA, the amplitudes and phases of fundamental and each harmonic can be computed using (20) as follows:(20)An=xn12+xn22,φn=tan-1⁡xn2xn1.Figure 6
Structure of WCA estimation scheme.Moreover, fundamental as well as each harmonic can be directly extracted from the final updated weight:(21)a^nk=xn1sin⁡ωnkTs+xn2cos⁡ωnkTs.


## 5. Simulation Results

In order to validate the proposed harmonic estimation algorithm, the simulation is firstly carried out. The simulation signal is(22)at=2sin⁡12πt-0.1+0.8sin⁡24πt+0.2+0.6sin⁡36πt-0.3+0.5sin⁡48πt+0.4+0.4sin⁡60πt-0.5+0.3sin⁡72πt+0.6+0.2sin⁡84πt-0.8+0.1sin⁡96πt+0.8.The parameters of the algorithm are selected asNsr=8, dmax=0.001, and max⁡iteration=24 in each trial. The sampling time is 1 ms, and the initial values are all set to zero. Figures 7 and 8 are the estimated amplitude and phase, respectively. Large fluctuation can be seen in the initial estimation, but they converge to stabilize after 1 s. This is because estimation starts from randomly generated solutions. It is perceived that the estimated value is close to its original value, except some slight difference. This indicates that harmonic estimation scheme based on WCA is able to estimate the harmonic information efficiently.Figure 7
Estimated amplitude of harmonic.Figure 8
Estimated phase of harmonic.


## 6. Experimental Results

In order to investigate the WCA harmonic estimation algorithm further, the shaking test is performed on the hydraulic shaking table in Figure1. When the acceleration excitation is 3sin⁡(12πt)m/s2, the frequency spectrum is shown in Figure 9. The harmonics present in a distorted waveform depends on the inherent nonlinearities of the system. So it is very difficult to predict the order of harmonics present in the waveform only by observing the signal waveform. For all practical purposes, we need to only think about specified number of harmonics that presented in the signal. It can be seen from the figure that the frequency of harmonics (up to eight harmonic) is integer multiples of the fundamental frequency. That is, the second harmonic is at 12 Hz, the third harmonic is at 18 Hz, the fourth harmonic at 24 Hz, and so on. The amplitude of each harmonic is presented in Table 1. They are computed offline by FFT from a period of sampling data.Table 1
Each harmonic’s amplitude.
A
1 (6 Hz)
A
2 (12 Hz)
A
3 (18 Hz)
A
4 (24 Hz)
A
5 (30 Hz)
A
6 (36 Hz)
A
7 (42 Hz)
A
8 (48 Hz)
2.426
0.067
0.144
0.017
0.249
0.055
0.182
0.119Figure 9
The frequency spectrum of response.To evaluate the effects of harmonic on response, total harmonic distortion (THD) is used to measure the level of distortion and it is defined as(23)THD=A22+A32+⋯+AN2A1.The value of THD based on FFT is 15.3%. Table2 gives the estimated amplitude and phase of each harmonic when they are in steady sate. The THD based on WCA is 5.4%, which is smaller than FFT.Table 2
The estimated amplitude and phase.
Harmonic order
Amplitude (m/s2)
Phase (rad)
Fundamental
2.68
-
0.67
Second harmonic
0.07
-
0.86
Third harmonic
0.15
0.68
Fourth harmonic
0.02
0.95
Fifth harmonic
0.25
1.37
Sixth harmonic
0.06
-
0.09
Seventh harmonic
0.19
-
0.47
Eighth harmonic
0.12
-
0.78Figures10 and 11 represent the estimated amplitude and phase of each harmonic. The results show that the estimation efficiency of amplitude is better than phase, but they converge to its steady state after 1 s.Figure 10
Estimated amplitude of harmonic.Figure 11
Estimated phase of harmonic.The waveform of each harmonic reconstructed from the estimated weight vector is shown in Figure12. As can be seen from the diagram, each harmonic is well estimated after a large fluctuation. This is helpful to compensate harmonic distortion.Figure 12
The estimated harmonic.Figure13 shows the actual signal and signal reconstructed from the harmonic weight vector estimated by WCA. It can be seen from the figure that the estimated signal generated based on estimated parameters matches the actual signal very well. This demonstrates that the proposed method can precisely estimate the harmonic waveform.Figure 13
The estimated acceleration.


## 7. Conclusions

This paper proposed the WCA based harmonic estimation algorithm. The algorithm can accurately estimate the amplitude and phase of each harmonic by minimizing the objective function. Unlike other harmonic estimation algorithms, the proposed scheme does not need system’s prior knowledge and this makes the algorithm suitable for online applications. Simulation and experimental results have shown that WCA is able to estimate the parameters of the harmonics in dynamic environment. What is more, the convergence performance of WCA is high enough for the algorithm to be designed for real-time estimation.



---

*Source: 7278589-2018-01-15.xml*