System identification of force transducers for dynamic measurements using particle swarm optimization

2.1. System set-up

The system set-up for dynamic calibrations of force transducers used in this investigation is shown in Fig. 1. This calibration facility is similar to the one developed by Kumme at Physikalisch-Technische Bundesanstalt (PTB). But different from Kumme’s method, the facility uses a sine mechanism instead of a shaker to generate sinusoidal motion. The sine mechanism has been successfully applied in the field of testing automobile shock absorber. When the crank disk is rotating, the frame and slider move orthogonally as illustrated in Fig. 1 [22]. With kinematic analysis, we can get the following functions:

where $s$ and $a$ are the displacement and acceleration of the frame respectively, $r$ is the rotation radius, $\omega$ is the rotational angular velocity, ${\phi }_0$ is the initial angle and $t$ is time. The force transducer to be calibrated is mounted on the frame and moves together with it according to Eq. (1).

Fig. 1System set-up for dynamic calibrations of force transducers

The displacement $s$ can be collected with a displacement transducer, and the acceleration $a$ is assumed to be the same for the frame and the base mass inside force transducer. The dynamic force could be obtained based on the determination of the inertia force as a known mass mounted on the force transducer [22].

2.1.1. Mathematical model

Force transducer can be described by the Voigt model, consisting of a spring mass system with stiffness $k_f$ and damping $b_f$ acting between the end mass $m_e$ and the base mass $m_b$ of the transducer [4]. When considering the influence of the coupling between the force transducer and the load mass $m_l$, another spring mass system with stiffness $k_c$ and damping $b_c$ could be employed. Therefore, in this investigation, the schematic diagram illustrated in Fig. 2 can be used to represent the force transducer which is mounted in the calibration set-up.

Fig. 2Schematic diagram of a force transducer

And the following differential equations could be established according to Fig. 2:

3a

$m_l{\ddot{x}}_l=-k_c\left(x_l-x_e\right)-b_c\left({\dot{x}}_l-{\dot{x}}_e\right)-m_{l}g,$

3b

$m_e{\ddot{x}}_e=k_c\left(x_l-x_e\right)+b_c\left({\dot{x}}_l-{\dot{x}}_e\right)-k_f\left(x_e-x_b\right)-b_f\left({\dot{x}}_e-{\dot{x}}_b\right)-m_{e}g,$

3c

$m_b{\ddot{x}}_b=k_f\left(x_e-x_b\right)+b_f\left({\dot{x}}_e-{\dot{x}}_b\right)+F_s-m_{b}g,$

where $x_l$, $x_e$, $x_b$ are the vertical movements of the load mass, the end mass and the base mass respectively, $F_s$ is the acting force. Since the measured signals in this investigation are collected from the force transducer to be calibrated and the displacement transducer connected to the frame, our goal is to obtain the relationship between $x_e-x_b$ and $x_b$. $x_e-x_b$ is proportional to the output signal $U_f$ of the force transducer to be calibrated, and $x_b$ can be collected with the displacement transducer.

To simplify system Eq. (3), take ${\ddot{X}}_l={\ddot{x}}_l+g$, ${\ddot{X}}_e={\ddot{x}}_e+g$, ${\ddot{X}}_b={\ddot{x}}_b+g$, ${Z_1=x}_l-x_e$ and ${Z_2=x}_e-x_b$. In addition, the problem can be further considered in the Laplace space by introducing the complex variable $s=j\omega$ to make the substitution $\dot{x}=sx$, and the abbreviations $K_f=k_f+s\bullet b_f$ and $K_c=k_c+s\bullet b_c$ will be employed, too [7]. Thus system Eq. (3) can be expressed as following:

4a

${\ddot{Z}}_1=-\left(\frac{K_c}{m_l}+\frac{K_c}{m_e}\right)\bullet Z_1+\frac{K_f}{m_e}\bullet Z_2,$

4b

${\ddot{Z}}_2=\frac{K_c}{m_e}\bullet Z_1-\left(\frac{K_f}{m_e}+\frac{K_f}{m_b}\right)\bullet Z_2-\frac{F_s}{m_b}.$

Take $u_1=m_l{m}_e/(m_l+m_e)\text{,}$$u_2=m_e{m}_b/(m_e+m_b)\text{,}$$F=-F_s/m_b\text{,}$ and transfer the algebraic equation system Eq. (4) into matrix formalism:

5

$\left[\begin{array}{cc}{s}^2+\frac{K_c}{u_1}& -\frac{K_f}{m_e}\\ -\frac{K_c}{m_e}& s^2+\frac{K_f}{u_2}\end{array}\right]\left[\begin{array}{c}{Z}_1\\ Z_2\end{array}\right]=\left[\begin{array}{c}0\\ F\end{array}\right].$

The solution of this system can be performed using Cramer’s rule [7]:

6

$Z_1=\frac{{∆}_1}{∆},Z_2=\frac{{∆}_2}{∆},$

where:

$∆=\left(s^2+\frac{K_c}{u_1}\right)\left(s^2+\frac{K_f}{u_2}\right)-\frac{K_f{K}_c}{m_e^2},{∆}_1=\frac{K_f}{m_e}\bullet F,{∆}_2=\left(s^2+\frac{K_c}{u_1}\right)\bullet F.$

To obtain the relationship between $x_e-x_b$ and $x_b$, take $Z_2/{\ddot{X}}_b$ as the target model:

7

$\frac{Z_2}{{\ddot{X}}_b}=1/\left(\frac{K_f}{m_b}-F\frac{∆}{{∆}_2}\right).$

When considering the influence of the coupling between the force transducer and the load mass, by substituting expressions of $K_f$, $F$, $∆$ and ${∆}_2$ into Eq. (7) we can obtain:

8

$\frac{Z_2\left(s\right)}{s^2{X}_b\left(s\right)}=-\frac{s^2+n_{1}s+n_2}{s^4+d_1{s}^3+d_2{s}^2+d_{3}s+d_4},$

where:

$d_1=\frac{b_c}{u_1}+\frac{b_f}{m_e},\mathrm{}\mathrm{}\mathrm{}\mathrm{}{d}_2=\frac{k_c}{u_1}+\frac{k_f}{m_e}+\frac{b_f{b}_c}{m_e{m}_l},\mathrm{}\mathrm{}\mathrm{}\mathrm{}{d}_3=\frac{k_f{b}_c+k_c{b}_f}{m_e{m}_l},\mathrm{}\mathrm{}\mathrm{}\mathrm{}$
$d_4=\frac{k_f{k}_c}{m_e{m}_l},\mathrm{}\mathrm{}\mathrm{}\mathrm{}{n}_1=\frac{b_c}{u_1},\mathrm{}\mathrm{}\mathrm{}\mathrm{}{n}_2=\frac{k_c}{u_1}.$

In the case where the coupling of the load mass to the force transducer is very stiff, the transfer function Eq. (7) can be simplified using the limits $k_c\to \infty$ and $b_c\to 0$ [7]:

9

$\frac{Z_2\left(s\right)}{s^2{X}_b\left(s\right)}\approx -\frac{1}{s^2+d_1^{\text{'}}s+d_2^{\text{'}}},$

where:

$d_1^{\text{'}}=\frac{b_f}{m_e+m_l},d_2^{\text{'}}=\frac{k_f}{m_e+m_l}.$

Since the system identification is finished in computer with application of advanced tools from Digital Signal Processing (DSP), we carry out $Z$ transform with bilinear method for $-\rho \bullet Z_2\left(s\right)/{\mathrm{s}}^2{X}_b\left(s\right)$, where the constant $\rho$ realizes the transformation of the displacement $Z_2$ into a force signal. And corresponding discrete transfer functions $H_1\left(z\right)$ and $H_2\left(z\right)$ for Eqs. (8) and (9) can be obtained as following:

10

$H_1\left(z\right)=\frac{b_0+b_1{z}^{-1}+b_2{z}^{-2}+b_3{z}^{-3}+b_4{z}^{-4}}{a_0+a_1{z}^{-1}+a_2{z}^{-2}+a_3{z}^{-3}+a_4{z}^{-4}},$

where:

$a_0=d_4{T}^4+2d_3{T}^3+4d_2{T}^2+8d_{1}T+16,a_1={4d}_4{T}^4+4d_3{T}^3-16d_{1}T-64,$
$a_2={6d}_4{T}^4-8d_2{T}^2+96,\mathrm{}\mathrm{}\mathrm{}\mathrm{}{a}_3={4d}_4{T}^4-4d_3{T}^3+16d_{1}T-64,$
$a_4=d_4{T}^4-2d_3{T}^3+4d_2{T}^2-8d_{1}T+16,b_0=n_2\rho T^4+2n_1\rho T^3+4\rho T^2,$
$b_1={4n}_2\rho T^4+4n_1\rho T^3,b_2={6n}_2\rho T^4-8\rho T^2,$
$b_3={4n}_2\rho T^4-4n_1\rho T^3,b_4=n_2\rho T^4-2n_1\rho T^3+4\rho T^2,$

and $T$ is the sampling time.

11

$H_2\left(z\right)=\frac{b_0^{\text{'}}+b_1^{\text{'}}{z}^{-1}+b_2^{\text{'}}{z}^{-2}}{a_0^{\text{'}}+a_1^{\text{'}}{z}^{-1}+a_2^{\text{'}}{z}^{-2}},$

where:

$a_0^{\text{'}}=d_2^{\text{'}}{T}^2+2d_1^{\text{'}}T+4,a_1^{\text{'}}=2d_2^{\text{'}}{T}^2-8,\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}{a}_2^{\text{'}}=d_2^{\text{'}}{T}^2-2d_1^{\text{'}}T+4$
$b_1^{\text{'}}=2{\rho T}^2,b_0^{\text{'}}=b_2^{\text{'}}=b_1^{\text{'}}/2$

3.1. Schematic for system identification

The main task of system identification is to search iteratively for the parameters of the modeled system such that the input-output relationship matches closely to that of the actual system [11]. The basic block diagram for system identification is shown in Fig. 3. System input $x\left(k\right)$ is given to both the unknown system to be identified and the modeled system. The output $y^{\text{'}}\left(k\right)$ mixed with a noise signal gives the final output $y\left(k\right)$ to the actual system$\mathrm{}$. On the other hand, the modeled system has an output of $\widehat{y}\left(k\right)$ for the same input. The difference $e\left(k\right)$ between the two output signals is used by the identifier to adjust the parameters [22]. Since traditional algorithms usually have difficulties in optimizing complex nonlinear systems, the identification of transfer Eqs. (10) and (11) will be finished with particle swarm optimization algorithm in this work.

Fig. 3Block diagram for PSO based system identification

3.2. PSO algorithm

Particle swarm optimization algorithm is initialized with a population of possible solutions called particles. Particles start flying from the initial positions through the search space with velocities exploring almost all optimal solutions. The velocity of each particle is dynamically adjusted according to its own flying experience and its companions’ flying experience, and the performance of each particle position is evaluated by a fitness function. During the flights, the best previous experience for each particle is stored in its memory and called the personal best (Pbest), and the best previous position among all particles is called the global best (Gbest) [10, 12, 22, 23].

Let $X_i=\left(x_{i1},x_{i2},\cdots ,x_{iN}\right)$ be the position of the $i$th particle, for Eq. (10) $N=10$, and for Eq. (11) $N=4$. Similarly, the velocity $V_i$, the personal best $X_i^{Pbest}$ and the global best $X^{Gbest}$ are represented as $V_i=\left(v_{i1},v_{i2},\cdots ,v_{iN}\right)\text{,}$$X_i^{Pbest}=\left(x_{i1}^{Pbest},x_{i2}^{Pbest},\cdots ,x_{iN}^{Pbest}\right)\text{,}$ and $X^{Gbest}=\left(x_1^{Gbest},x_2^{Gbest},\cdots ,x_N^{Gbest}\right)$ respectively. The particles are manipulated according to the following equations [10, 24]:

where acceleration coefficients $c_1\text{,}$$c_2$ are positive constant parameters with the constraint $c_1+c_2\le 4$ that control the maximum step size. Typically, $c_1$, $c_2$ are both suggested to be 2.0. In [22] the authors also have studied the influence of acceleration coefficients $c_1$, $c_2$ on PSO performance, and the results indicate that parameters $c_1=\text{2.0}$ and $c_2=\text{2.0}$ can get the tradeoff between identification accuracy and convergence speed. Therefore, $c_1=\text{2.0}$ and $c_2=\text{2.0}$ will be taken in this investigation. $r_1$, $r_2$ are uniformly distributed random variables in the range [0, 1]. $w$ is the inertia weight and controls the impact of the previous velocity of the particle on its current one. In general, a bigger $w$ can guarantee the globally search ability, and avoid sinking into a local optimal solution. On the contrary, a smaller $w$ is good for local search, and can guarantee the convergence of the arithmetic. Thus, a self-adaptive adjusted strategy for inertia weight $w$ as follows can be taken [13]:

where $j$ is the current iteration generation, ${iter}_{\mathrm{m}\mathrm{a}\mathrm{x}}$ is the maximum iteration times, $w_{\mathrm{m}\mathrm{i}\mathrm{n}}$ and $w_{\mathrm{m}\mathrm{a}\mathrm{x}}$ are the minimum inertia weight and maximum inertia weight, respectively. The value of inertia weight $w$ in PSO is linearly decreased according to Eq. (14), and $w_{\mathrm{m}\mathrm{a}\mathrm{x}}$, $w_{\mathrm{m}\mathrm{i}\mathrm{n}}$ are initial and final values for the inertia weight [23]. In this way, PSO tends to have more global search ability at the beginning of the run while having more local search ability near the end of the run [22]. $w_{\mathrm{m}\mathrm{a}\mathrm{x}}=\text{0.9}$ and $w_{\mathrm{m}\mathrm{i}\mathrm{n}}=\text{0.4}$ are generally taken in previous studies [14, 22].

In the system identification problem, PSO algorithm tries to minimize the error $e\left(k\right)$ by adjusting the parameters of the modeled system. And mean square error (MSE) of time samples between the outputs of the actual system and the designed system as given by Eq. (15) is always considered as the fitness function [9, 11]. In mathematical statistics, mean square error is the mathematical expectation of differences between the estimated values and the truth values, and it is a convenient method to measure the average error:

where $M$ is the number of data points, $y\left(k\right)$ is the output of the actual system, $\widehat{y}\left(k\right)=\frac{\widehat{B}\left(z^{-1}\right)}{\widehat{A}\left(z^{-1}\right)}x\left(k\right)$ is the output of the modeled system with $\widehat{A}\left(z^{-1}\right)={\widehat{a}}_0+{\widehat{a}}_1{z}^{-1}+\cdots +{\widehat{a}}_n{z}^{-n}\text{,}$$\widehat{B}\left(z^{-1}\right)={\widehat{b}}_0+{\widehat{b}}_1{z}^{-1}+\cdots +{\widehat{b}}_m{z}^{-m}$. The smaller the MSE value is, the higher accuracy of the modeled system is.

In general, the main problems of the PSO algorithm are slow convergence and premature convergence [25]. To deal with the problems, exponential transformation is introduced to the fitness function for helping quickly search the global optimal solution:

And with Eq. (16), the fitness function can be converted to a positive indicator in the range [0, 1].

3.3. Algorithm procedure

The implementation of PSO algorithm with exponential transformation (PSO+ algorithm) applied to system identification in this work is summarized as follows:

Step 1. Generate $M$ input-output data points to form the system to be identified.

Step 2. Set acceleration coefficients $c_1$, $c_2$, the maximum iteration times ${iter}_{\mathrm{m}\mathrm{a}\mathrm{x}}$, and the scope of the inertial factors $w_{\mathrm{m}\mathrm{i}\mathrm{n}}$, $w_{\mathrm{m}\mathrm{a}\mathrm{x}}$ that appear in Eqs. (12) and (14).

Step 3. Initialize the position and velocity for each particle. In this investigation, particle position corresponds to the parameters of transfer functions (10) and (11).

Step 4. Remember current particle position as Pbest of each particle. Find the particle with the best value of fitness function and remember its position as Gbest.

Step 5. Generate inertia weight $w$ and random numbers $r_1$, $r_2$in Eq. (12), and update the velocities and positions according to Eqs. (12) and (13).

Step 6. Evaluate fitness values and update Pbest and Gbest.

Step 7. Return to step 5 until the total number of iterations is achieved.

Step 8. Output the identified parameters namely Gbest and end the program.

4.1. Validation of PSO+ algorithm

For testing the performance of PSO+ algorithm, two benchmark functions listed in Table 1 which are similar to transfer functions Eqs. (10) and (11) are considered for the case studies. The inputs for the two studies are both $5\mathrm{s}\mathrm{i}\mathrm{n}(2\pi \times 4t)$. And the two cases are simulated with PSO and PSO+ algorithms in MATLAB with the following parameters: $c_1=c_2=\text{2}$, $w_{\mathrm{m}\mathrm{i}\mathrm{n}}=\text{0.4}$, $w_{\mathrm{m}\mathrm{a}\mathrm{x}}=\text{0.9}$, population size, maximum position and maximum velocity are 150, 2 and 0.1 respectively. For case I, $N=\text{4}$ and the maximum generation is set to 1000. For case II, $N=\text{10}$ and the maximum generation is 2000. The transfer functions and their identification results are shown in Table 1. The outputs of actual systems and modeled systems for cases I and II are illustrated in Figs. 4 and 7 respectively. Figs. 5 and 6 are the curves of fitness functions for PSO and PSO+ in case I while Figs. 8 and 9 are those in case II.

From Figs. 4-9 and Table 1 we can see that, identification results of the two cases are in accordance with each other. Identification accuracies of PSO in the two cases are slightly higher than those of PSO+, but identification errors of PSO+ are very small too, since the MSE values listed in Table 1 are in the 10^-17 order or less. From this point we know that, PSO+ algorithm is effective for system identification.

When comparing the convergence processes shown in Figs. 5-6 and Figs. 8-9, we can find that the convergence speed of PSO+ algorithm is obviously faster than that of PSO. For case I, PSO takes almost 880 iterations to converge to the global best while PSO+ takes less than 100 iterations. And for case II, PSO takes 1416 iterations to converge to the global best while PSO+ takes only 850 iterations. Consequently, we reach a conclusion that PSO+ algorithm can be used for system identification and it converges faster than conventional PSO.

Fig. 4Outputs for case I

Fig. 5Fitness value of global best particle for PSO in case I

Fig. 6Fitness value of global best particle for PSO+ in case I

Fig. 7Outputs for case II

Fig. 8Fitness value of global best particle for PSO in case II

Fig. 9Fitness value of global best particle for PSO+ in case II

4.2. Performance study of PSO+ in the presence of the additive noise

For studying the ability of PSO+ algorithm to operate in the presence of the additive noise, the benchmark function of case I listed in Table 1 is employed to finish another simulation study. The input for this study is still $5\mathrm{s}\mathrm{i}\mathrm{n}(2\pi \times 4t)$, and white noises of different signal to noise ratios (SNRs) are added to the output. The parameters of PSO+ algorithm is set as those in case I.

Fig. 10Output in the case of 50 dB SNR

Fig. 11Output in the case of 40 dB SNR

Fig. 12Output in the case of 30 dB SNR

Fig. 13Output in the case of 20 dB SNR

The outputs of actual systems and modeled systems are illustrated in Figs. 10-14 respectively. The MSE values are shown in Fig. 15.

From Figs. 10-14 we can see that, the higher the SNR is, the smoother the output curve is. And in this situation, the identification accuracy of PSO+ algorithm is higher, too. The relationship between identification accuracy and SNR is shown in Fig. 15. The MSE values in the cases of 10 dB-50 dB SNRs are 80.1880, 7.9021, 0.9229, 0.0921 and 0.0092 respectively. In all these situations, the outputs of the modeled systems can follow the changes of the actual outputs. But in order to obtain the higher identification accuracy, filtering process for the actual outputs is suggested to be finished before identifying with PSO+ algorithm.

Fig. 14Output in the case of 10 dB SNR

Fig. 15MSEs in the cases of different SNRs

5.1. Experiment I

Firstly, the experiment is carried out with a force transducer of 5 kN. The sine mechanism is operating at a frequency of $4\text{Hz}$ and the transducer signals are sampled at 20 kHz. After filtering and resampling, PSO+ algorithm is employed to identify the simplified function Eq. (11) (Case III) and function Eq. (10) (Case IV). 200 input-output data points are taken for system identification and another 400 data points are taken for testing. The parameters of PSO+ algorithm are set as follows: $c_1=c_2=\text{2}$, $w_{\mathrm{m}\mathrm{i}\mathrm{n}}=\text{0.4}$, $w_{\mathrm{m}\mathrm{a}\mathrm{x}}=\text{0.9}$, population size and maximum position are 200 and 100 respectively. When identifying the simplified function Eq. (11), Recursive Least Squares (RLS) method is also used to do the identification besides the PSO+ algorithm. Alternative form of function Eq. (11) for RLS identification is shown below:

where:

Fig. 17 shows the identification errors of RLS and PSO+ for cases III and IV, and the convergence processes of fitness values for PSO+ algorithm are illustrated in Fig. 18. Testing results of the identified models for cases III and IV are shown in Fig. 19, and Fig. 20 shows the testing errors of the identified models with PSO+ algorithm. And for the empirical analysis of the testing results, MSE is employed as the accuracy evaluation standard which is shown in Table 2.

Fig. 17Identification errors for cases III and IV

Fig. 18Fitness values of PSO+ for cases III and IV

Fig. 19Testing results for cases III and IV

Fig. 20Testing errors of PSO+ for cases III and IV

The results shown in Figs. 17 and 19 indicate that, RLS is not suitable for the identification of function Eq. (11) since big errors exist in the processes of identifying and testing. The testing MSE values shown in Table 2 also tell this, as MSE of RLS is more than 1000 times bigger than MSE of PSO+. So, when compared with RLS, the identification performance of PSO+ for case III is encouraging.

When comparing the results of PSO+ for cases III and IV shown in Figs. 19-20 and Table 2, we can find that simplification between functions Eqs. (10) and (11) almost has no influence on identification accuracy. The testing MSE values of PSO+ for case III and IV are 0.0191 and 0.0186 respectively, identification accuracy in case IV is slightly better than that in case III. Nevertheless, in Fig. 18, we can see that the convergence speed of PSO+ for case III is faster than that for case IV. Therefore, by comprehensively considering the identification accuracy and convergence speed, we reach a conclusion that the simplification between functions Eq. (10) and (11) is feasible. And with the identified parameters obtained in case III we can calculate out the stiffness and damping of the force transducer conveniently. In this experiment, the stiffness and damping of the 5 kN force transducer are 0.291×10⁸ Nm^-1 and 17136 Nsm^-1 respectively.

5.2. Experiment II

For investigating the generality of this identification method, another force transducer of 12.5 kN is taken to carry out a same experiment. The parameters of PSO+ algorithm is set as those in experiment I. The identifications of 12.5 kN transducer for function Eq. (11) and function Eq. (10) are denoted as case V and case VI respectively. The identification results are shown in Figs. 21-24 and Table 3.

Fig. 21Identification errors for cases V and VI

Fig. 22Fitness values of PSO+ for cases V and VI

Fig. 23Testing results for cases V and VI

Fig. 24Testing errors of PSO+ for cases V and VI

The results of experiment II are basically in accordance with those of experiment I. Identification error of RLS in Fig. 21 is still very big, and this leads to that the testing performance of RLS illustrated in Fig. 23 is not acceptable. Its MSE value is up to 0.9304 while MSE of PSO+ is only 0.0020. Then compare the identification performances of PSO+ for cases V and VI. As shown in Fig. 22, the convergence speed of PSO+ for case V is also faster than that for case VI obviously. The testing accuracies of PSO+ algorithm for cases V and VI are almost the same since MSE values of cases V and VI are 0.0020 and 0.0017 respectively. And with the identified parameters obtained in case V, the stiffness and damping of 12.5 kN force transducer are calculated as 0.246×10⁸ Nm^-1 and 3032 Nsm^-1 respectively.