Fuzzy sliding mode controller design for semi-active seat suspension with neuro-inverse dynamics approximation for MR damper

2.1. Overview of quarter-car semi-active seat suspension model

Considering the quarter-car model to be with high accuracy in analyzing the suspension dynamics, it is employed to model the semi-active suspension in this paper. Fig. 1 presented a three DOF model of quarter-car suspension system containing the seat suspension with a MR damper. In this figure the car body, seat and human body are included as the sprung masses: $m_v$ and $m_s$, and the vertical dynamics of tire and axle is often considered by introducing an unsprung mass $m_t$ and spring $k_t$. MR damper is placed between the seat and the car body to form the seat suspension together with a spring and a damper.

Based on Newton second law, the dynamic equations of seat suspension system is:

where $m_t$, $m_v$ and $m_s$ are unsprung mass, quarter car body mass and seat (plus human body) mass respectively. $k_t$, $k_v$ and $k_s$ are the stiffness coefficients of the tire, quarter car suspension and seat suspension respectively. $c_v$ and $c_s$ are the damping coefficients of quarter car suspension and seat suspension respectively. $F_d$ is the semi-active damping force created by the MR damper. $z_0$, $z_t$, $z_v$ and $z_s$ are the road excitation, vertical displacements of car axle, body and seat, respectively.

Based on Eq. (1), the state equation of the system is:

where:

Fig. 1Model of quarter-car suspension

2.2. MR damper modeling and analysis

Bouc-wen hysteresis model paralleled with dashpot and spring is originally used to formulate the MR damper. It can describe the hysteretic nonlinearity of the MR damper, but it is unable to describe the nonlinear and saturation dependence of the magnetic field yielded by the direct drive current. A modified Bouc-wen hysteresis model, proposed by Spencer [32], effectively overcomes the above drawback and precisely describes the nonlinear saturated characteristic of the MR damper. The modified Bouc-wen hysteresis model is presented in Fig. 2 and it is finally obtained by adding a serial viscous damping with the Bouc-wen model and further paralleling a linear spring to the serialized structure.

Fig. 2Modified Bouc-Wen model of MR damper

In this paper, the modified Bouc-Wen model is used to describe the mechanical properties of the MR damper. The model introduces two internal variables, and it constructs a differential equation model with 14 parameters to be determined. According to Fig. 2, the mathematical equations of the modified Bouc-wen model are presented as follows:

where $\dot{u}=-\eta (u-v)$, $k_1$ is the stiffness of the damper accumulator. $c_0$ is the viscous damping observed when larger velocities is represented. $c_1$ is a dashpot, included in the model to produce the roll-off that was observed in the experimental data at low velocities. $k_0$ is presented to control the stiffness at large velocities, and $x_0$ is the initial displacement of spring $k_1$ is associated with the nominal damper force due to the accumulator. $u$ is given as the output of a first-order filter. $v$ is the commanded voltage sent to the current driver.

As for the RD-1005-3 damper produced by Lord Corporation, its parameters are chosen as follows. $\alpha =$ 963 N/cm, $c_0=$ 53N·S/cm, $k_0=$ 14 N/cm, $c_1=$ 930 N·s/cm, $k_1=$ 5.4 N/cm, $\gamma =$ 200 cm², $\beta =$ 200 cm^-2, $n=$ 2, $A=$ 207, and $x_0=$ 18.9 cm, the response of the proposed model at 2.5 Hz is obtained as shown in Fig. 3. It can be seen from Fig. 3 that MR damper can provide the damping effect in the plane of I and III quadrant of velocity-force plane, unlike the active actuator in all the four quadrants. Therefore, the output of MR damper has to track the desired damping force only when the expected force and velocity have same sign, otherwise it should output the least damping force, so the formula described is:

where $f\left(t\right)$ is the damping force of MR damper, $f_e\left(t\right)$ is the desired force which is obtained by using suspension controller, $f_{e\mathrm{m}\mathrm{i}\mathrm{n}}$ is the minimal damping force corresponding to the zero input current, however it is not a constant value changing with the instant velocity.

If the MR damper controller employs the switching control method of Eq. (6), the input voltage of the MR damper switches between the minimum and maximum without being a continuous adjustable control signal, which limits the performance of the MR damper. In this paper neural network is used to simulate the inverse model of MR damper, and it is further be as the nonlinear controller of MR damper to create a continuous signal for MR damper, which will be discussed in next section.

Fig. 3Experimentally obtained response of the model

a) Force vs. displacement

b) Force vs. velocity

2.3. Neuro inverse model approximating of MR damper

The reverse model of MR damper is defined as solving the voltage corresponding to input MR damper after the desired damping force is obtained by FSMCEF control algorithm. The aim of the reverse model is to make MR damper track this desired force as possible. The inverse model can be described using the following nonlinear function:

where $\widehat{v}\left(t\right)$ is the input voltage, namely the output of MR damper reverse model. $\theta$ is the neural network weight vector, determined by training process. $\phi$ is the input vector as $\varphi \left(k\right)=\left[\widehat{v}\right(k-1),...,\widehat{v}(k-n_v),x(k),...,x(k-n_x),F(k),...,F(k-n_f\left)\right]$, in which $k$ is the $k$th(current) time step, $n_v$, $n_x$ and $n_f$ are the previous time step numbers of input voltage, damping force and displacement respectively.

The BP neural network can approximate an arbitrary nonlinear continuous function with arbitrary precision, and it is used for the inverse-dynamics approximation of MR damper is shown in Fig. 4.

Fig. 4Block diagram of neuro inverse dynamics model for MR damper

The typical BP network is divided into three layers, which are input layer, hidden layer and output layer. For the inverse dynamics model of Eq. (7), the network input layer is set with 9 nodes and the hidden layer with 20 nodes (see Fig. 5, $p=$ 20), the output layer has one node, standing for the input voltage of the MR damper.

The output of the hidden node is:

The output of the output node is:

Fig. 5Detailed neural network for the inverse dynamics approximation

This neural network training includes namely mode suitable transmission, error back propagation, memory training and learning convergence. The detailed training procedures are as follows.

(1) Initialization:

The weights of $\left\{w_{ij}\right\}$, $\left\{v_{jt}\right\}$ and threshold of $\left\{{\theta }_j\right\}$, $\left\{{\gamma }_t\right\}$ are all set as random values in (–1, 1).

(2) Randomly pick up a pair of samples for network training.

(3) Calculate the output of the hidden layer using Eq. (8).

(4) Calculate the output layer using Eq. (9).

(5) Calculate the average error of the output layer: $d_t^k=(y_t^k-C_t)C_t(1-C_t^{})$.

(6) Calculate the hidden layer general error: $e_j=(\sum _{t=1}^q{d}_t^k{v}_{jt})\cdot b_j(1-b_j^{})$.

(7) Modify the output layer weights and thresholds:

(8) Modify hidden layer weights and thresholds:

(9) Take out the next pair of sample and return Step (3) to repeat until completing all training samples.

(10) Determine whether a global error is less than the preset value, otherwise, to return to Step (2) to continue until meeting the requirements.

The displacement and the input voltage, as the network input, are generated by Gaussian white noise with frequency ranges of 0-3 Hz and 0-4 Hz respectively. A data set of 10000 points used for training and validation are created by 500 Hz sampling frequency in 20 s sample time. These data are used for training the network and another data set including 1200 points is further created to validate the training. The BP neural network training and validating process is shown in Fig. 6 and the control voltage between the predict output and desired output shown in Fig. 7. The BP network prediction error compared with the desired output is shown in Fig. 8.

Fig. 6BP neural network training and validating process

Fig. 7Control voltage between the predict output and the desired output

Fig. 8BP network prediction error

Table 1 shows the BP predict output data, the desired output data and their error. The sum of the absolute value of 1200 errors is 15.9525. Table 1 presented the detailed result of some samples, and the results show that the BP neural network predict output values substantially agree with the actual situation, and the errors are small, demonstrating that this model is highly with the project reality and the neural network model is effective.

3.1. Skyhook reference model

Herein the skyhook model is used as the reference to form the fuzzy sliding control algorithm, and it is presented in Fig. 10.

The dynamics equations are derived based on Newton second law as follows:

where $z_{rt}$, $z_{rv}$ and $z_{rs}$ are the road excitation, vertical displacements of unsprung mass, car body and seat respectively, which are the corresponding variables in the reference system compared with the plant system in Fig. 1. $c_{sh}$ is the damping coefficient of “sky-hook” damper.

Fig. 10Skyhook reference model

The state vector for the reference system is taken as $Z_r=[z_{rs},{\dot{z}}_{rs},z_{rv},{\dot{z}}_{rv},z_{rt},{\dot{z}}_{rt}{]}^T$, and the output vector is taken as $Y=[z_{rs},{\dot{z}}_{rs},z_{rv},{\dot{z}}_{rv}{]}^T$. According to Eq. (12) the state equations are established as:

where:

3.2. Error dynamics model used for SMC and FSMCEF

Both the sliding mode controller and the fuzzy sliding mode controller are designed to make the actual seat suspension motion to track the reference mode, and so they are based on the dynamic errors between the seat suspension and the skyhook reference model. Based on the above seat model and the reference model, the seat suspension displacement error, its integral and its differential (velocity error) are taken as the control variables, and they form the general tracking error vector$e$as $e=[\begin{array}{lll}{e}_1& e_2& e_3\end{array}{]}^T={\left[\begin{array}{lll}\int (z_s-z_{rs})& z_s-z_{rs}& {\dot{z}}_s-{\dot{z}}_{rs}\end{array}\right]}^T$ and its differential is $\dot{e}={\left[\begin{array}{lll}{z}_s-z_{rs}& {\dot{z}}_s-{\dot{z}}_{rs}& {\ddot{z}}_s-{\ddot{z}}_{rs}\end{array}\right]}^T$. so, the error dynamic equation is obtained as:

where:

3.3. Sliding mode control based on pole placement

The switching surface is taken as:

As for $s=c_1{e}_1+c_2{e}_2+e_3=0$, Eq. (14) can be written to partitioned matrix form as:

As for $\dot{s}=c_1{e}_2+c_2{e}_3+{\dot{e}}_3=0$, Eq. (16) can be written as:

The characteristic polynomial for Eq. (18) is $D\left(\lambda \right)={\lambda }^2+c_2\lambda +c_1$. To obtain the value of $c_1$ and $c_2$, its characteristic roots are equal to the given poles.

The chief problem of pole assignment is rationally to determine the desired closed-loop poles set. The standard form of the second-order system transfer function is:

The two closed-loop poles are $s_{\mathrm{1,2}}=-\zeta {\omega }_n\pm j{\omega }_n\sqrt{1-{\zeta }^2}$, and the system works in less damping state ($0<\zeta <1$) which makes these two poles to be conjugate complex roots located in the left half plane of $s$ domain and to have the appropriate oscillation and short transition process. For the three orders system are in Eqs. (16) and (17), the desired poles number are $n=$ 3. The conjugate pole pair of $s_1$ and $s_2$ are selected as the dominant poles, and the third is a non-dominant one. The poles placement of our closed-loop system is completed to ensure the two dynamic performance indices: peak time $t_p$ and overshoot$\sigma$ %. These two indices are set as $\sigma \le$ 15 %, $t_p\le$ 0.7, which determines the dominant poles be at –2.7326±4.4886i and the non-dominant pole at –20. The corresponding parameters are $\zeta =\text{0.52}$ and ${\omega }_n=\text{5.255}$. So, the switching function coefficient vector is $c=\left[\text{68}\begin{array}{l}\end{array}\text{4}\begin{array}{l}\end{array}\text{1}\right]$. The system performance is mainly determined by these two dominant poles and the non-dominant pole only produces minimal effect.

When system is in sliding mode motion, $s=0$, $ds/dt=0$ and:

The equivalent control for system into fuzzy sliding mode or sliding mode is $u^{*}$, and:

In order to improve the dynamic quality of the movement, the approaching mode employs the constant speed reaching law as:

where $\epsilon =$ 3.

The final sliding mode control law is taken as:

So, the desired real-time variable damping force is as:

where $F_{d_{sw}}=u_{sw}=\epsilon m_s\mathrm{s}\mathrm{g}\mathrm{n}\left(s\right)$, $\epsilon =$ 3.

3.4. Fuzzy sliding mode control

The fuzzy logic control is further added to overcome sliding mode controller “chattering” problem. The connection of the sliding mode controller and the fuzzy logic controller is shown in Fig. 11, which formed the final fuzzy sliding mode controller. The detailed content of the fuzzy control block in Fig. 11 is presented in Fig. 12. Its inputs are $s\left(e\right)$ and $\dot{s}\left(e\right)$, with one output $\epsilon$ sent to the sliding mode controller.

Fig. 11Block diagram of FSMCEF system

The fuzzy controller first change the range of $s\left(e\right)$ and $\dot{s}\left(e\right)$, namely from their original ranges of [–0.04, 0.03] and [–6×10^-3, +8×10^-3] both to the new range of [–6, +6] for further discretization and fuzzification. The corresponding conversion equation is:

Fig. 12The structure of fuzzy controller

The variables after conversion are $S$ and $SC$ (Fig. 12), and they are further discretized and fuzzified to form fuzzy sets $\underset{\_}{S}$, $\underset{\_}{S}C$. In this process $S$ and $SC$ are classified into seven grades, forming seven fuzzy subsets, including: NL (Negative Large), PL (Positive Large), NM (Negative Medium), PM (Positive Medium), NS (Negative Small), PS (Positive Small), NE (0). The $S$ and $SC$ of domain $X$ and $Y$ are belonging to 7 fuzzy subsets, respectively. Similarly, the output value $\underset{\_}{\epsilon }$ are also ranked into seven fuzzy subsets: NL, PL, NM, PM, NS, PS, NE.

For this double input and single output fuzzy controller, the control rules can be written as the following form: If $S={\underset{\_}{S}}_i$ and $SC=\underset{\_}{S}{C}_j$ then $U={\underset{\_}{U}}_{ij}$, ($i=$1, 2,…, 7, $j=i=$1, 2,.., 7), where ${\underset{\_}{S}}_i$, $\underset{\_}{S}{C}_j$ are input fuzzy sets, and ${\underset{\_}{U}}_{ij}$ is output fuzzy sets.

These fuzzy sets conditional statements can be summed up in a fuzzy relation $\underset{\_}{R}$, and $\underset{\_}{R}=\underset{ij}{\cup }({\underset{\_}{E}}_i\times \underset{\_}{E}{C}_j)\times {\underset{\_}{U}}_{ij}$. According to each inference rules, the corresponding fuzzy relations, ${\underset{\_}{R}}_1,{\underset{\_}{R}}_2,...,{\underset{\_}{R}}_n$ can be calculated. So, the total of the whole system corresponding fuzzy control rule $\underset{\_}{R}$ is:

The final fuzzy rules are shown in Table 2, and according to Table 2 the 3-D input-output relation diagram of the fuzzy controller is obtained as shown in Fig. 13.

Fig. 133-D diagram of Fuzzy control rules

When $\underset{\_}{R}$ is determined, according to $\underset{\_}{S}=\{-6,-5,\cdots +5,+6\}$, $\underset{\_}{S}C=\{-6,-5\cdots +5,+6\}$ and synthetic fuzzy reasoning rules, the corresponding fuzzy sets of controls is $\underset{\_}{U}=(\underset{\_}{E}\times \underset{\_}{E}C)\circ R$, and:

The final fuzzy sliding or sliding mode control law is taken as:

3.5. Fuzzy sliding mode control with expansion factor

Equidistant domain partitioning method is used in fuzzy control generally. When the error is large, the system has sufficient error resolution, and it is shown as the “big error” dotted line in Fig. 14. When the error is small, the system response only changes around “ZO” corresponding to the original fuzzy partition, and other fuzzy subsets obviously do not work. Ideally, when the error is reduced, the domain of the fuzzy controller should be able to make self-adaptation adjustment. The accuracy of the fuzzy controller is related to the number of output variables and fuzzy rules. Supposing the input is $n$-dimensional, the fuzzy control rule of the universe is divided into $m$; the total rule number is $m^n$, if the fuzzy subset of the domain is divided into smaller, the number of fuzzy control rules will increase exponentially, increasing the difficulty of making rules. Therefore, under the premise of not affecting the control effect, we should try to use less fuzzy subset to reduce the number of fuzzy control rules. With the same form of rules, the universe shrinks as the error becomes smaller, and the universe expands as the error increases. The contraction of the domain is equivalent to increasing the fuzzy control rules to improve the control accuracy. The function of the scaling factor $\alpha \left(x\right)$ transforms the domain into $[-\alpha (x)E,\alpha (x\left)E\right]$, where $\alpha \left(x\right)$ is a continuous function of the error variable $x$. The appropriate domain expansion factor $\alpha \left(x\right)$ is chosen so that the range of the universe changes with the error, which can realize the adaptive implementation of the expansion of the domain, without the need for other auxiliary algorithms and increase the control rules.

Fig. 14Adjustment of the domain

Let $X_i=[-E,E]$ be the universe of input variable $x_i(i=\mathrm{1,2},\dots ,n)$, $Y=[-U,U]$ is the universe of output variable $y$; ${\psi }_i=\left\{A_{ij}\right\}$ is the fuzzy partition on $X_i$, ${\varphi }_j=\left\{B_j\right\}$ is the fuzzy partition on $Y$, $1\le j\le m$. As ${\psi }_i$, ${\varphi }_j$ are the linguistic variables, fuzzy inference rule $R$ can be formed:

where $x_i$ is $A_{ij}$ peak, and $y_j$ is $B_j$ peak $(i=\mathrm{1,2},\dots ,n)$, $(j=\mathrm{1,2},\dots ,m)$, The fuzzy control system of Eq. (31) can be expressed as an $n$-piece piecewise interpolation function:

The variable universe is that domain $X_i$ and $Y$ can be adjusted independently with the change of input variable $x_i$ and output $y$, respectively:

where ${\alpha }_i\left(x_i\right)$ and $\beta \left(y\right)$ is the domain expansion factor. In contrast to the variable universe, the original universe $X_i$ and $Y$ is called the initial universe. Eq. (32) can be expressed as $n$-piece dynamic interpolation function:

where $x\left(t\right)\triangleq [x_1\left(t\right),$$x_2\left(t\right),\dots ,x_n\left(t\right){]}^T$, and ${\alpha }_i\left(x_i\right)$ is chosen as:

or:

In this paper $\alpha \left(x\right)=1-1/\sqrt{1+k{x}^2}$, $k=$ 10⁴. $\beta \left(y\right)$ is chosen as:

where $K_i$ is the proportionality constant, and $\beta \left(0\right)$ is according to the actual situation, usually try $\beta \left(0\right)=$ 1. The control law of the variable universe fuzzy controller is:

$X=\left[-E,E\right]$, $Y=\left[-D,D\right]$ is two-dimensional input domain, respectively, and $Z=\left[-U,U\right]$ is the output domain. When $X$ and $Y$ are relatively independent, we can get the expansion factor $\alpha \left(x\right)$, $\beta \left(y\right)$ and $Z$ expansion factor $\gamma \left(z\right)$. But in most cases, $Y$ and $X$ are related. If $X$ is the error domain and $Y$ is often the domain of error variation, $Y=\left(-\dot{E},\dot{E}\right)$ and $\beta \left(y\right)$ should be defined on $X\times Y$ and $\beta \left(y\right)=\beta (x,y)$, then the input and output of the domain expansion factor are:

Since the error variation depends on the error, $\beta$ can be simply taken as $\beta \left(y\right)$, Eq. (41) can be rewritten as:

3.6. Stability analysis