Design of a magnetorheological (MR) suspension damper for an agricultural tractor seat based on an adaptive neuro-fuzzy inference system (ANFIS) and active disturbance rejection control (ADRC)

To enhance the vibration isolation functionality of the seat suspension system, particularly in the context of magnetorheological (MR) seat dampers, a pioneering semi-active controller integrating an adaptive neuro-fuzzy inference system (ANFIS) and active disturbance rejection control (ADRC) was devised. Firstly, the Bouc–Wen model was employed to establish the dynamic model of the damper based on experimental data, thereby ensuring an accurate description of the actual physical behaviours. Subsequently, the ANFIS technology was employed to develop an inverse model, which detailed the structural design and training process of the inverse model. This resulted in the achievement of precise prediction and control of the damper behaviour. Ultimately, the ANFIS inverse model was integrated with the designed ADRC controller to create an innovative control scheme for the seat suspension system of a two-degrees-of-freedom dynamic model, and a simulation analysis was conducted. The simulation results demonstrate that the root-mean-square (rms) value of the vertical vibration acceleration of the ADRC-controlled suspension system decreased by 68.9 % and 34.4 % in comparison to proportional–integral–derivative (PID) control and passive control, respectively. The rms value of the dynamic disturbance of the ADRC-controlled suspension system decreased by 50.0 % and 28.6 % compared to PID control and passive control, respectively. This verifies the performance of the proposed controller, particularly in the precise control of damping force, demonstrating outstanding effectiveness.

1 Introduction

In recent years, there has been a notable shift in the engineering community's attention towards semi-active vibration control systems (Ahn et al., 2023; Bhowmik and Debnath, 2024; Luan et al., 2023). These systems are capable of providing vibration reduction effects that are comparable to those of active control devices while simultaneously maintaining low energy consumption and the reliability of passive devices. Among the numerous semi-active control systems, magnetorheological (MR) dampers have garnered considerable attention due to their simplicity, reliability, and robustness. Notably, semi-active control suspension systems utilizing MR dampers have been extensively investigated in the context of vehicle vibration comfort (Wang et al., 2024; Kumar and Bhushan, 2024; Xie and Hua, 2024).

The majority of research conducted on MR semi-active suspension systems is concentrated in three primary areas: the dynamic modelling of MR dampers, the modelling of suspension systems and their vibration transmission characteristics, and the simulation analysis and experimental validation of semi-active suspension control systems. Miao et al. (2024) employed the normalized Bouc–Wen model to ascertain the nonlinear characteristics of the MR rotary dampers. By means of numerical simulations, the effects of MR semi-active suspension and passive suspension on the dynamic performance of tracked vehicles under a variety of operating conditions were investigated. In comparison to passive suspension, the MR semi-active suspension exhibited enhanced vibration reduction capabilities across diverse road grades and vehicle speeds, thereby substantiating the superiority of this suspension system (Miao et al., 2024). The study by de Brett et al. (2023) proposed and validated a methodology for the analysis and modelling of the nonlinear vibration transfer of vehicle suspension dampers at low audio frequencies. Utilizing experimental data and model validation, they conducted nonlinear dynamic modelling and analysis of the front and rear suspension dampers of the test vehicle. The authors selected and fitted appropriate models based on experimental data and conducted a comprehensive analysis of various nonlinear effects, including the influence of friction and the mechanical characteristics of piston valves on the dynamic performance of the dampers (de Brett et al., 2023). Turgay Ergin and Meral Özarslan Yatak proposed and validated a proportional–integral–derivative (PID) controller based on particle swarm optimization (PSO) for a semi-active suspension system (SASS) of a quadricycle. They utilized artificial neural networks (ANNs) to simulate the model of the MR damper and conducted a series of experiments to assess the efficacy of the proposed algorithm, thereby confirming the superior performance of the controller in practical applications (Ergin and Yatak, 2023). The mechanical characteristics of MR dampers and the establishment of their forward and inverse models have constituted a research focus as the actuating elements in semi-active suspension systems. This is evidenced by the extensive literature on the subject (Gong et al., 2023; Tantray, 2023; Bahar et al., 2024). The forward model predicts the output damping force based on the input current and the relative motion state of the piston. The primary applications of this method are to elucidate the operational mechanism of MR dampers, to substitute for actual dampers in simulations, and to function as force sensors in actual control. The inverse model is employed to forecast the optimal control current in accordance with the desired control force and the piston's motion state, a process of considerable practical consequence in control applications (Abd Elwahed et al., 2024).

The forward dynamic models of MR dampers currently include the Bingham model (Zhao et al., 2024), the nonlinear hysteretic viscous model (Zhao et al., 2018), the phenomenological model (Jiang et al., 2023), the Dahl model (Aguirre et al., 2012), the LuGre model (Jiménez and Álvarez-Icaza, 2005), the sigmoid generalized hysteresis model, the hyperbolic tangent model, the polynomial model, the neural network model, and the neural-network–fuzzy model (Archakam and Muthuswamy, 2021; Lu et al., 2023; Wei et al., 2021). However, these models seldom take into account the impact of excitation characteristics on damper properties. Given the inherent variability of excitation characteristics in practical applications, models identified under specific excitation conditions may not fully align with the needs of real-world applications. In light of these considerations, this paper puts forth a dynamic model for MR dampers that is accurate, simple, highly adaptive, and reversible. Based on experimental data, a Bouc–Wen model was constructed to represent the forward mechanical performance of MR dampers, and an adaptive neuro-fuzzy inference system (ANFIS) was employed to develop an inverse dynamic model of MR dampers. The inverse model was subsequently integrated with the active disturbance rejection control (ADRC) for joint simulation and optimal control, and the resulting simulation results were subjected to analysis.

2 The Bouc–Wen model for MR dampers

In system dynamics simulation, it is essential to establish a parametric dynamic model of the magnetorheological damper to ensure that the simulation results align with the actual situation. The Bouc–Wen model is particularly effective in addressing the nonlinear hysteresis characteristics of the magnetorheological damper, and as a result, it has become a prevalent choice for simulating and analysing actual problems. The Bouc–Wen model comprises a hysteresis system, a spring, and a viscous damping element in parallel, as shown in Fig. 1.

Figure 1 Schematic diagram of the Bouc–Wen model structure.

To ensure the accurate and effective control of the output damping force of the magnetorheological damper, it is essential to guarantee the precision and reliability of the model. To assess the damping-force–velocity and damping-force–displacement characteristics of the magnetorheological damper, we utilize an electro-hydraulic servo fatigue machine, as illustrated in Fig. 2.

Figure 2 The damping characteristic test of the MR damper.

The Bouc–Wen model demonstrates superior efficacy in representing the relationship between the output damping force of the MR damper and the relative displacement and velocity of the piston. Its expression for the damping force is

$$\begin{array}{}\text{(1)}& \left\{\begin{array}{l}F=c_{\mathrm{0}}\dot{x}+k_{\mathrm{0}}\left(x-x_{\mathrm{0}}\right)+\mathit{\alpha }z,\\ \dot{z}=-\mathit{\gamma }\left|\dot{x}\right|z{\left|z\right|}^{n-\mathrm{1}}-\mathit{\beta }\dot{x}{\left|z\right|}^n+Ax,\end{array}\right.\end{array}$$

where F represents the damping of the MR damper, c₀ denotes the viscous damping coefficient, k₀ signifies the stiffness coefficient, α stands for the proportional adjustment parameter of the hysteresis force in the damping force, and x₀ is the relative equilibrium position offset displacement. Additionally, z represents the hysteresis variable. The variables x and $\dot{x}$ represent the relative displacement and relative velocity of the MR damper piston, respectively. The parameter n is referred to as the rounding coefficient, while γ is the width adjustment coefficient of the hysteresis model. The height adjustment coefficient of the hysteresis model is represented by β, and the scaling factor is designated as A. The adjustment coefficient of the hysteresis model, β, represents the height adjustment coefficient of the hysteresis model, whereas A denotes the scale factor.

The Bouc–Wen model comprises eight unknown parameters, and the code for identification methods such as the genetic algorithm and particle swarm algorithm is both cumbersome and complex. This paper therefore employs the nonlinear least squares method in the “Design Optimization” toolbox in Simulink to carry out the parameter identification. The iterative algorithms are automatically called using the experimental data to achieve a high degree of correlation between the experimental and simulation values, thus identifying all the parameters. The nonlinear least squares method identifies the optimal solution to the estimated objective function Q by locating the minimum value of the objective function Q of the error sum of squares. The identified model is expressed mathematically as follows:

$$\begin{array}{}\text{(2)}& y=f(x{}^{\prime }x{}^{\prime \prime },\mathrm{\cdots },\mathit{\theta }{}^{\prime }\mathit{\theta }{}^{\prime \prime },\mathrm{\cdots }),\end{array}$$

where y is the output of the system; $x{}^{\prime }x{}^{\prime \prime },\mathrm{\cdots }$ is the input series of the system; and $\mathit{\theta }{}^{\prime }\mathit{\theta }{}^{\prime \prime },\mathrm{\cdots }$ is the parameter series. Once the parameters have been estimated, the mathematical expression of the model, f, is known, and the data are obtained experimentally as ($x_{\mathrm{1}}{}^{\prime }{x}_{\mathrm{1}}{}^{\prime \prime }$, ⋯, y₁) ($x_{\mathrm{2}}{}^{\prime }{x}_{\mathrm{2}}{}^{\prime \prime }$, ⋯, y₂) ⋯ ($x_n{}^{\prime }{x}_n{}^{\prime \prime }$, ⋯, y_n). The objective function Q, which represents the sum of squared errors for the nonlinear model, is defined as follows:

$$\begin{array}{}\text{(3)}& Q=\sum _{n=\mathrm{1}}^N{\left[y_n-f(x_n{}^{\prime },x_n{}^{\prime \prime },\mathrm{\cdots },\mathit{\theta }{}^{\prime },\mathit{\theta }{}^{\prime \prime },\mathrm{\cdots })\right]}^{\mathrm{2}}.\end{array}$$

The eight unknown parameters in the Bouc–Wen model were identified through the use of measured piston relative displacements, velocities, and currents as inputs and the damping force as the output. The identification flowchart is shown in Fig. 3.

Figure 3 Bouc–Wen parameter identification flowchart.

Firstly, the constructed Bouc–Wen model was employed to import the experimental data, which were obtained through the “Parameter Estimation” toolbox. The input data consisted of displacement, velocity, and current values, while the output data represented the damping force. The initial values of the identification parameters were established, and the identification range was defined based on experimental data and prior experience, as illustrated in Table 1.

Table 1 Bouc–Wen parameter identification and initial values and ranges.

Subsequently, the optimization option was selected via the toolbox, and the nonlinear least squares method was selected for iterative computation. The cost function was then chosen, and the squared difference was employed for parameter estimation. The parameter estimation was terminated according to the proximity of the simulation curves to the experimental data, specifically when Q≤0.001.

Following the completion of the Simulink modelling of the Bouc–Wen model, the experimental data obtained from an input excitation frequency of 2 Hz; an amplitude of 10 mm; and input currents of 0, 0.5, 1.0, and 1.5 A were selected for the identification of the kinetic model parameters. The Bouc–Wen model comprises a total of eight unknown parameters, each of which exerts a distinct influence on the model's behaviour. The preliminary identification of model parameters revealed that the values of the two parameters α and c₀ undergo significant changes in response to variations in current. Consequently, α and c₀ were set to vary in accordance with changes in working conditions, while the remaining parameters were maintained at a constant value. The fitting of the Bouc–Wen model indicated that the relationship between α and c₀ and the current can be expressed by a quadratic polynomial. The synthesis of the other six fixed parameters yielded the Bouc–Wen parameter fitting results, which are presented in Table 2.

Table 2 Bouc–Wen model parameters.

Once the parameters had been identified, the resulting data were plotted and compared with the test results obtained at an excitation frequency of 2 Hz and an amplitude of 10 mm, as illustrated in Fig. 4.

Figure 4 Results of the MR damper model parameters and experimental results.

Figure 4 illustrates the experimental test results, represented by the solid line, and the model parameter fitting results, depicted by the dashed line. As can be seen from the figure, the curve obtained from parameter fitting using the Bouc–Wen model (applied using Simulink parameter identification) is in general agreement with the curve obtained from the characteristic test of the damper, and the fitting effect is satisfactory.

3 ANFIS inverse dynamics modelling

The damping force generated by the MR damper is mainly determined by the input current, the relative piston acceleration, and the relative piston displacement. Only when the input current can be directly controlled can the MR damper function properly. Therefore, it is important to obtain the command current according to the desired force in practical applications. Since the ANFIS has a general approximation ability for nonlinear systems (Karaboga and Kaya, 2019), in this section the ANFIS technique is applied to the inverse model building of the MR damper to build the inverse model of the control current and to calculate the current required to obtain the desired damping force.

3.1 Adaptive neuro-fuzzy inference system

The ANFIS is a combination of an adaptive network and a fuzzy inference system, which inherits the interpretability of a fuzzy inference system and the learning ability of an adaptive network; it utilizes the learning mechanism of a neural network to automatically extract the rules from the input and output sample data, constituting the adaptive neuro-fuzzy controller, and can autonomously adjust the parameters of the system through hybrid algorithms so as to make the output of the system closer to the actual output (Bilgundi et al., 2022). The ANFIS consists of five functional modules: input layer, affiliation layer, rule layer, decision layer, and output layer, as shown in Fig. 5.

Figure 5 Structure of ANFIS.

In Fig. 5, x and y are input variables that form two if–then rules. Rule 1: if x is A1 and y is B1, then continue; rule 2: if x is A2 and y is B2, then continue (where p₂, q_i, and r_i are variable coefficients).

The following is a description of the ANFIS layers+:

• Input layer. For adaptive nodes, fuzzy input features with an affiliation function to get the degree of affiliation in the interval and the output are denoted as

  $$\begin{array}{}\text{(4)}& {\mathit{\theta }}_i^{\mathrm{1}}\left(t\right)={\mathit{\mu }}_{Ai}\left(x\right),\end{array}$$

  where x is the input of node i, and μ_Ai(x) is the affiliation function, which usually has a bell curve, cosine function, or Gaussian function. In this paper, a Gaussian function is used, i.e.

  $$\begin{array}{}\text{(5)}& {\mathit{\mu }}_{Ai}\left(x\right)=e^{-{\left(\frac{x-c_i}{\mathrm{2}{\mathit{\sigma }}_i^{\mathrm{2}}}\right)}^{\mathrm{2}}},{\mathit{\mu }}_{Bi}\left(y\right)=e^{-{\left(\frac{y-c_i}{\mathrm{2}{\mathit{\sigma }}_i^{\mathrm{2}}}\right)}^{\mathrm{2}}},\end{array}$$

  where c_i and σ_i are parameters set by the Gaussian function.

• Affiliation layer. This is a fixed node, denoted by ∏. The trigger strength ω_i of each rule is obtained by multiplying the affiliation degree of each feature, and the output is

  $$\begin{array}{}\text{(6)}& {\mathit{\theta }}_i^{\mathrm{2}}={\mathit{\omega }}_i={\mathit{\mu }}_{Ai}\left(x\right)\times {\mathit{\mu }}_{Bi}\left(y\right).\end{array}$$

• Rule layer. The classification process is performed to normalize the trigger strength of each rule obtained from the previous layer (i.e. it is the corresponding probability of use of the rule), which is used to indicate the triggering ratio of the rule (Bilgundi et al., 2022).

  $$\begin{array}{}\text{(7)}& {\mathit{\theta }}_i^{\mathrm{3}}={\stackrel{\mathrm{‾}}{\mathit{\omega }}}_i={\displaystyle \frac{{\mathit{\omega }}_i}{{\sum }_{j=\mathrm{1}}^n{\mathit{\omega }}_j}},\end{array}$$

  where ${\stackrel{\mathrm{‾}}{\mathit{\omega }}}_i$ is the normalized trigger strength.

• Decision layer. For adaptive nodes, the output of each node is simply the product of the normalized trigger strength and a first-order polynomial, and the output is the following equation:

  $$\begin{array}{}\text{(8)}& {\mathit{\theta }}_i^{\mathrm{4}}={\stackrel{\mathrm{‾}}{\mathit{\omega }}}_i{f}_i={\stackrel{\mathrm{‾}}{\mathit{\omega }}}_i\left(p_{i}x+q_{i}y+r_i\right).\end{array}$$

• Output layer. It is a fixed node where the exact output is obtained by deburring, characterized by Σ. A weighted average of the results of each rule gives the final output of the model, as in the following equation:

  $$\begin{array}{}\text{(9)}& {\mathit{\theta }}_i^{\mathrm{5}}={\sum }_i{\stackrel{\mathrm{‾}}{\mathit{\omega }}}_i\cdot f_i.\end{array}$$

ANFIS overcomes the black-box characteristics of simple neural networks, as well as the incompleteness and roughness of the inference rules in the fuzzy inference process, and maximally simplifies the data processing with adaptive, self-organizing, and self-learning characteristics.

In this paper, a hybrid optimization method combining backpropagation and least squares is employed for training the model network. It is assumed that the parameters of the ANFIS system can be decomposed into two sets of outputs and inputs. In the event that the outputs are provided, the system will output a linear function. Consequently, the linear parameters can be identified through the least squares method. Once the inputs have been identified, the parameters in the output can be updated using the gradient method. The matrix form of the least squares method is presented in Eq. (10).

$$\begin{array}{}\text{(10)}& \begin{array}{l}\left[\begin{array}{llllll}{\stackrel{\mathrm{‾}}{\mathit{\omega }}}_{\mathrm{1}}^{\left(\mathrm{1}\right)}{x}^{\left(\mathrm{1}\right)}& {\stackrel{\mathrm{‾}}{\mathit{\omega }}}_{\mathrm{1}}^{\left(\mathrm{1}\right)}{y}^{\left(\mathrm{1}\right)}& {\stackrel{\mathrm{‾}}{\mathit{\omega }}}_{\mathrm{1}}^{\left(\mathrm{1}\right)}& {\stackrel{\mathrm{‾}}{\mathit{\omega }}}_{\mathrm{2}}^{\left(\mathrm{1}\right)}{x}^{\left(\mathrm{1}\right)}& {\stackrel{\mathrm{‾}}{\mathit{\omega }}}_{\mathrm{2}}^{\left(\mathrm{1}\right)}{y}^{\left(\mathrm{1}\right)}& {\stackrel{\mathrm{‾}}{\mathit{\omega }}}_{\mathrm{2}}^{\left(\mathrm{1}\right)}\\ {\stackrel{\mathrm{‾}}{\mathit{\omega }}}_{\mathrm{1}}^{\left(\mathrm{2}\right)}{x}^{\left(\mathrm{2}\right)}& {\stackrel{\mathrm{‾}}{\mathit{\omega }}}_{\mathrm{1}}^{\left(\mathrm{2}\right)}{y}^{\left(\mathrm{2}\right)}& {\stackrel{\mathrm{‾}}{\mathit{\omega }}}_{\mathrm{1}}^{\left(\mathrm{2}\right)}& {\stackrel{\mathrm{‾}}{\mathit{\omega }}}_{\mathrm{2}}^{\left(\mathrm{2}\right)}{x}^{\left(\mathrm{2}\right)}& {\stackrel{\mathrm{‾}}{\mathit{\omega }}}_{\mathrm{2}}^{\left(\mathrm{2}\right)}{y}^{\left(\mathrm{2}\right)}& {\stackrel{\mathrm{‾}}{\mathit{\omega }}}_{\mathrm{2}}^{\left(\mathrm{2}\right)}\\ \mathrm{\cdots }& \mathrm{\cdots }& \mathrm{\cdots }& \mathrm{\cdots }& \mathrm{\cdots }& \mathrm{\cdots }\\ {\stackrel{\mathrm{‾}}{\mathit{\omega }}}_{\mathrm{1}}^{\left(n\right)}{x}^{\left(n\right)}& {\stackrel{\mathrm{‾}}{\mathit{\omega }}}_{\mathrm{1}}^{\left(n\right)}{y}^{\left(n\right)}& {\stackrel{\mathrm{‾}}{\mathit{\omega }}}_{\mathrm{1}}^{\left(n\right)}& {\stackrel{\mathrm{‾}}{\mathit{\omega }}}_{\mathrm{2}}^{\left(n\right)}{x}^{\left(n\right)}& {\stackrel{\mathrm{‾}}{\mathit{\omega }}}_{\mathrm{2}}^{\left(n\right)}{y}^{\left(n\right)}& {\stackrel{\mathrm{‾}}{\mathit{\omega }}}_{\mathrm{2}}^{\left(n\right)}\end{array}\right]\\ \left[\begin{array}{l}{p}_{\mathrm{1}}\\ q_{\mathrm{1}}\\ r_{\mathrm{1}}\\ p_{\mathrm{2}}\\ q_{\mathrm{2}}\\ r_{\mathrm{2}}\end{array}\right]=\left[\begin{array}{l}{f}_{\mathrm{p}}^{\left(\mathrm{1}\right)}\\ f_{\mathrm{p}}^{\left(\mathrm{2}\right)}\\ \mathrm{\cdots }\\ f_{\mathrm{p}}^{\left(\mathrm{n}\right)}\end{array}\right],\end{array}\end{array}$$

where n is the total number of training data input–output pairs; $f_{\mathrm{p}}^{\left(\mathrm{n}\right)}$ is the predicted output of the network, where the backward parameters $[p_{\mathrm{1}},q_{\mathrm{1}},r_{\mathrm{1}},p_{\mathrm{2}},q_{\mathrm{2}},r_{\mathrm{2}}{]}^T$ are derived from Eq. (11); and $f_{\mathrm{d}}^{\left(n\right)}$ is the measured output.

$$\begin{array}{}\text{(11)}& \begin{array}{l}{\left[\begin{array}{llllll}{\stackrel{\mathrm{‾}}{\mathit{\omega }}}_{\mathrm{1}}^{\left(\mathrm{1}\right)}{x}^{\left(\mathrm{1}\right)}& {\stackrel{\mathrm{‾}}{\mathit{\omega }}}_{\mathrm{1}}^{\left(\mathrm{1}\right)}{y}^{\left(\mathrm{1}\right)}& {\stackrel{\mathrm{‾}}{\mathit{\omega }}}_{\mathrm{1}}^{\left(\mathrm{1}\right)}& {\stackrel{\mathrm{‾}}{\mathit{\omega }}}_{\mathrm{2}}^{\left(\mathrm{1}\right)}{x}^{\left(\mathrm{1}\right)}& {\stackrel{\mathrm{‾}}{\mathit{\omega }}}_{\mathrm{2}}^{\left(\mathrm{1}\right)}{y}^{\left(\mathrm{1}\right)}& {\stackrel{\mathrm{‾}}{\mathit{\omega }}}_{\mathrm{2}}^{\left(\mathrm{1}\right)}\\ {\stackrel{\mathrm{‾}}{\mathit{\omega }}}_{\mathrm{1}}^{\left(\mathrm{2}\right)}{x}^{\left(\mathrm{2}\right)}& {\stackrel{\mathrm{‾}}{\mathit{\omega }}}_{\mathrm{1}}^{\left(\mathrm{2}\right)}{y}^{\left(\mathrm{2}\right)}& {\stackrel{\mathrm{‾}}{\mathit{\omega }}}_{\mathrm{1}}^{\left(\mathrm{2}\right)}& {\stackrel{\mathrm{‾}}{\mathit{\omega }}}_{\mathrm{2}}^{\left(\mathrm{2}\right)}{x}^{\left(\mathrm{2}\right)}& {\stackrel{\mathrm{‾}}{\mathit{\omega }}}_{\mathrm{2}}^{\left(\mathrm{2}\right)}{y}^{\left(\mathrm{2}\right)}& {\stackrel{\mathrm{‾}}{\mathit{\omega }}}_{\mathrm{2}}^{\left(\mathrm{2}\right)}\\ \mathrm{\cdots }& \mathrm{\cdots }& \mathrm{\cdots }& \mathrm{\cdots }& \mathrm{\cdots }& \mathrm{\cdots }\\ {\stackrel{\mathrm{‾}}{\mathit{\omega }}}_{\mathrm{1}}^{\left(n\right)}{x}^{\left(n\right)}& {\stackrel{\mathrm{‾}}{\mathit{\omega }}}_{\mathrm{1}}^{\left(n\right)}{y}^{\left(n\right)}& {\stackrel{\mathrm{‾}}{\mathit{\omega }}}_{\mathrm{1}}^{\left(n\right)}& {\stackrel{\mathrm{‾}}{\mathit{\omega }}}_{\mathrm{2}}^{\left(n\right)}{x}^{\left(n\right)}& {\stackrel{\mathrm{‾}}{\mathit{\omega }}}_{\mathrm{2}}^{\left(n\right)}{y}^{\left(n\right)}& {\stackrel{\mathrm{‾}}{\mathit{\omega }}}_{\mathrm{2}}^{\left(n\right)}\end{array}\right]}^{-\mathrm{1}}\\ \left[\begin{array}{l}{f}_{\mathrm{d}}^{\left(\mathrm{1}\right)}\\ f_{\mathrm{d}}^{\mathrm{2}}\\ \mathrm{\cdots }\\ f_{\mathrm{d}}^{\left(\mathrm{n}\right)}\end{array}\right]=\left[\begin{array}{l}{p}_{\mathrm{1}}\\ q_{\mathrm{1}}\\ r_{\mathrm{1}}\\ p_{\mathrm{2}}\\ q_{\mathrm{2}}\\ r_{\mathrm{2}}\end{array}\right]\end{array}\end{array}$$

The discrepancy between the desired and predicted outputs is conveyed from the output layer to the input layer, where it is utilized to update the weights through the backpropagation algorithm, as illustrated in Eq. (12).

$$\begin{array}{}\text{(12)}& {\stackrel{\mathrm{‾}}{\mathit{\omega }}}_i^{\left(k\right)}(M+\mathrm{1})={\stackrel{\mathrm{‾}}{\mathit{\omega }}}_i^{\left(k\right)}\left(M\right)+(f_{\mathrm{d}}^{\left(k\right)}-f_{\mathrm{p}}^{\left(k\right)})\end{array}$$

According to Eq. (12), the weights of the input layer can be updated as shown in Eq. (13):

$$\begin{array}{}\text{(13)}& {\mathit{\omega }}_i^{\left(k\right)}(M+\mathrm{1})=\left\{\begin{array}{l}{\stackrel{\mathrm{‾}}{\mathit{\omega }}}_i^{\left(k\right)}\left(M\right)x+{\stackrel{\mathrm{‾}}{\mathit{\omega }}}_i^{\left(k\right)},\\ {\stackrel{\mathrm{‾}}{\mathit{\omega }}}_i^{\left(k\right)}\left(M\right)y+{\stackrel{\mathrm{‾}}{\mathit{\omega }}}_i^{\left(k\right)},\end{array}\right.\end{array}$$

where k is the input–output training pair, and M denotes each layer backward from the initial stage of output backpropagation.

Figure 6 The scheme of the ANFIS for modelling the inverse dynamics of the MR damper.

3.2 Training and validation of the ANFIS inverse model for MR dampers

In order to construct an accurate ANFIS inverse model, it is necessary to collect input–output data pairs that contain a sufficient amount of information. This information is then used to train the ANFIS model, after which the trained model is validated.

In the construction of an input–output dataset, the ANFIS develops a fuzzy inference system (FIS) whose membership function parameters are calibrated through the implementation of a hybrid algorithm. In general, the greater the number of inputs, the greater the accuracy of the inverse model. However, as the number of inputs increases, the inverse model will become increasingly complex, resulting in a significant increase in training time (Ghenai et al., 2022). In order to achieve an equilibrium between the accuracy and time consumption of the model, the inputs to the inverse model are selected to be the current displacement x(k), the current velocity $\dot{x}\left(k\right)$, the previous velocity $\dot{x}(k-\mathrm{1})$, the current desired damping force f(k), and the previous desired damping force f(k−1), while the output is the command current I(k). Figure 6 illustrates the application of the ANFIS approach to simulate the inverse dynamics of the MR damper. The displacement input is a Gaussian white noise signal with a frequency between 0 and 13 Hz and an amplitude of ±20 mm. The command current input was generated from a Gaussian white noise signal with a range of 0 to 1.5 A and a frequency of 0 to 2 Hz. The desired damping force was generated by the Bouc–Wen model, which was constructed in the previous section based on the displacement and command current inputs. The data were collected for a period of 10 s and subsequently sampled at a rate of 2000 Hz, resulting in the generation of 20 000 data points. The initial 10 000 data points were designated as training data, while the subsequent 10 000 data points were utilized as a verification set.

Figure 7 The validation flowchart of the ANFIS inverse model.

Figure 8 Validation of the ANFIS inverse model of the MR damper for training data: (a) the command current predicted by the ANFIS model and (b) the force predicted from the command current.

The construction of an accurate ANFIS inverse model hinges on the establishment of a comprehensive and precise training dataset, which serves as a cornerstone of the modelling process. In this study, data points are randomly selected from all test cases of the damper's dynamic mechanical characteristics experiment, including excitation amplitudes of 10, 20, 30, and 40 mm; excitation frequencies of 0.5 to 2 Hz; and input currents ranging from 0 to 1.5 A. This process resulted in a total dataset of 16 400 data points.

To enhance the ANFIS model's inference capability under varying current conditions, this paper, based on the test data obtained from experiments, establishes sweep frequency excitation with an excitation amplitude range of 10 to 40 mm and a frequency range of 0.1 to 2 Hz, according to the working range of the MR damper. By differentiating the displacement signal to obtain the velocity signal with a time step of 0.01 s, a total of 3000 datasets were generated through 30 s of simulation. The output current training data are generated by random signals with amplitudes ranging from 0 to 1.5 A. Based on the unchanged excitation data, five sets of random current arrays were generated, forming a total of 15 000 data points.

To validate the inverse dynamic neuro-fuzzy model, two validation datasets were analysed, namely training data and test data. The validation flowchart for these two cases is shown in Fig. 7. First, the target current, displacement, and velocity are input into Bouc-Wen_1 to generate the target force. Then, the target force, displacement, and velocity are input into the ANFIS inverse model to generate the predicted current, which is compared to the target current in the time domain. Finally, the predicted current, displacement, and velocity are input into Bouc-Wen_2 to generate the predicted force, which is compared to the target force in the time domain. The training data validation example is shown in Fig. 8. As seen in Fig. 8a, the predicted command current can reasonably track the target command current, and the damping force generated by the predicted command current is consistent with the damping force generated by the target command current, as shown in Fig. 8b.

The validation example for the test data is shown in Fig. 9. As shown in Fig. 9a, the accuracy of the test data is not as high as that of the training data. As seen in Fig. 9b, the damping force generated by the predicted command current can effectively track the damping force generated by the target command current. This meets the requirements of the MR damper inverse model, as the inverse model is primarily used to control the damping force of the MR damper.

Figure 9 Validation of the ANFIS inverse model of the MR damper for checking data: (a) the command current predicted by the ANFIS model and (b) the force predicted from the command current.

4 Semi-active seat suspension control design for agricultural tractors

4.1 Agricultural tractor seat suspension model

The two-degrees-of-freedom dynamic model of the MR damper seat suspension system is shown in Fig. 10. In the figure, m₁ and m_s represent the masses of the human body and the seat, respectively; k₁ and c₁ denote the stiffness and damping of the human body; k_s and c_MR indicate the stiffness of the seat suspension system and the equivalent damping of the MR damper, respectively; and z₁, z_s, and z₀ represent the displacements of the human body, seat, and base input excitation, respectively.

Figure 10 The two-degrees-of-freedom (2-DOF) human–seat model.

Figure 11 Schematic diagram of ANFIS-based MR damper and ADRC controller.

Table 3 Control parameter settings for PID and ADRC.

In accordance with Newton's second law, the seat system dynamics are represented by the model illustrated in Eq. (14).

$$\begin{array}{}\text{(14)}& \left\{\begin{array}{l}{m}_{\mathrm{1}}{\ddot{z}}_{\mathrm{1}}+k_{\mathrm{1}}\left(z_{\mathrm{1}}-z_{\mathrm{s}}\right)+c_{\mathrm{1}}\left({\dot{z}}_{\mathrm{1}}-{\dot{z}}_{\mathrm{s}}\right)=\mathrm{0},\\ m_{\mathrm{s}}{\ddot{z}}_{\mathrm{s}}-k_{\mathrm{1}}\left(z_{\mathrm{1}}-z_{\mathrm{s}}\right)-c_{\mathrm{1}}\left({\dot{z}}_{\mathrm{1}}-{\dot{z}}_{\mathrm{s}}\right)+k_{\mathrm{s}}\left(z_{\mathrm{s}}-z_{\mathrm{0}}\right)\\ +F_{\mathrm{MR}}=\mathrm{0},\end{array}\right.\end{array}$$

where F_MR represents the damping force of the seat MR damper. According to modern control theory, Eq. (14) is rewritten into the state equation form as shown in Eq. (15).

$$\begin{array}{}\text{(15)}& \left\{\begin{array}{c}\dot{X}=AX+BU,\\ Y=CX+DU,\end{array}\right.\end{array}$$

where state variable $X=[z_{\mathrm{2}},z_{\mathrm{1}},z_{\mathrm{s}},{\dot{z}}_{\mathrm{2}},{\dot{z}}_{\mathrm{1}},{\dot{z}}_{\mathrm{s}}{]}^T$, input variable $U=[z_{\mathrm{0}},F_{\mathrm{MR}}{]}^T$, and output variable $Y=[{\dot{z}}_{\mathrm{2}},{\dot{z}}_{\mathrm{1}},{\dot{z}}_{\mathrm{s}}{]}^T$.

$$\begin{array}{l}\dot{X}=\left[\begin{array}{cccc}\mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\\ -\frac{k_{\mathrm{1}}}{m_{\mathrm{1}}}& \frac{k_{\mathrm{1}}}{m_{\mathrm{1}}}& -\frac{c_{\mathrm{1}}}{m_{\mathrm{1}}}& \frac{c_{\mathrm{1}}}{m_{\mathrm{1}}}\\ \frac{k_{\mathrm{1}}}{m_{\mathrm{s}}}& -\frac{k_{\mathrm{1}}+k_{\mathrm{s}}}{m_{\mathrm{s}}}& \frac{c_{\mathrm{1}}}{m_{\mathrm{s}}}& -\frac{c_{\mathrm{1}}}{m_{\mathrm{s}}}\end{array}\right]\\ \left[\begin{array}{l}{z}_{\mathrm{1}}\\ z_{\mathrm{s}}\\ {\dot{z}}_{\mathrm{1}}\\ {\dot{z}}_{\mathrm{s}}\end{array}\right]+\left[\begin{array}{ll}\mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}\\ \frac{k_{\mathrm{s}}}{m_{\mathrm{s}}}& -\frac{\mathrm{1}}{m_{\mathrm{s}}}\end{array}\right]\left[\begin{array}{l}{z}_{\mathrm{0}}\\ F_{\mathrm{MR}}\end{array}\right],\end{array}$$

$$\begin{array}{l}Y=\left[\begin{array}{llll}-\frac{{\mathrm{k}}_{\mathrm{1}}}{m_{\mathrm{1}}}& \frac{{\mathrm{k}}_{\mathrm{1}}}{m_{\mathrm{1}}}& -\frac{{\mathrm{c}}_{\mathrm{1}}}{m_{\mathrm{1}}}& \frac{{\mathrm{c}}_{\mathrm{1}}}{m_{\mathrm{1}}}\\ \frac{{\mathrm{k}}_{\mathrm{1}}}{m_{\mathrm{s}}}& -\frac{{\mathrm{k}}_{\mathrm{1}}+k_{\mathrm{s}}}{m_{\mathrm{s}}}& \frac{{\mathrm{c}}_{\mathrm{1}}}{m_{\mathrm{s}}}& -\frac{{\mathrm{c}}_{\mathrm{1}}}{m_{\mathrm{s}}}\end{array}\right]\\ \left[\begin{array}{c}{z}_{\mathrm{1}}\\ z_{\mathrm{s}}\\ {\dot{z}}_{\mathrm{1}}\\ {\dot{z}}_{\mathrm{s}}\end{array}\right]+\left[\begin{array}{ll}\mathrm{0}& \mathrm{0}\\ \frac{k_{\mathrm{s}}}{m_{\mathrm{s}}}& -\frac{\mathrm{1}}{m_{\mathrm{s}}}\end{array}\right]\left[\begin{array}{l}{z}_{\mathrm{0}}\\ F_{\mathrm{MR}}\end{array}\right].\end{array}$$

Table 4 Main parameters of the human–seat model for agricultural tractors.

Figure 12 Road surface excitation.

Figure 13 Input current of the MR damper.

Figure 14 Comparison of actual damping force and desired damping force.

Figure 15 Seat dynamic disturbance comparison.

Figure 16 Vertical seat acceleration comparison.

Table 5 Seating suspension control rms value analysis.

4.2 ADRC controller design

From the two-degrees-of-freedom human–seat model, we have the following:

$$\begin{array}{}\text{(16)}& \left\{\begin{array}{l}{m}_{\mathrm{1}}{\ddot{z}}_{\mathrm{1}}+k_{\mathrm{1}}\left(z_{\mathrm{1}}-z_{\mathrm{s}}\right)+c_{\mathrm{1}}\left({\dot{z}}_{\mathrm{1}}-{\dot{z}}_{\mathrm{s}}\right)=\mathrm{0},\\ m_{\mathrm{s}}{\ddot{z}}_{\mathrm{s}}-k_{\mathrm{1}}\left(z_{\mathrm{1}}-z_{\mathrm{s}}\right)-c_{\mathrm{1}}\left({\dot{z}}_{\mathrm{1}}-{\dot{z}}_{\mathrm{s}}\right)+k_{\mathrm{s}}\left(z_{\mathrm{s}}-z_{\mathrm{0}}\right)\\ +F_{\mathrm{MR}}=\mathrm{0},\end{array}\right.\end{array}$$

where m₁ denotes the human body, and k₁ and c₁ are fixed; by forming f in Eq. (16), Eq. (16) is converted to a second-order model:

$$\begin{array}{}\text{(17)}& {\ddot{z}}_{\mathrm{s}}={\displaystyle \frac{k_{\mathrm{s}}{z}_{\mathrm{s}}}{m_{\mathrm{s}}}}-{\displaystyle \frac{\left(k_{\mathrm{s}}{z}_{\mathrm{0}}+F_b\right)}{m_{\mathrm{s}}}}-{\displaystyle \frac{F_{\mathrm{MR}}}{m_{\mathrm{s}}}}.\end{array}$$

Let $bu=-F_{\mathrm{MR}}/m_{\mathrm{s}}$, where b is the input gain whose value is unknown, b₀ is the nominal value, and u is the input signal of the MR damper; let y=z_s be the output of the MR damper, $a=k_{\mathrm{s}}/m_{\mathrm{s}}$ be the coefficients, and $\mathit{\omega }=-(k_{\mathrm{s}}{z}_{\mathrm{0}}+F_b)/m_{\mathrm{s}}$ be the external perturbation of the MR seat dynamic system. The final equation is simplified into the following equation:

$$\begin{array}{}\text{(18)}& \ddot{y}=ay+\mathit{\omega }+bu.\end{array}$$

The problem of second-order ADRC is to design the feedback controller such that y tracks the reference input signal r. Replacing the true value of b with the nominal value b₀ defines the total perturbation as $f=ay+\mathit{\omega }+(b-b_{\mathrm{0}})u$; then introducing the state variables x₁=y and $x_{\mathrm{2}}=\dot{y}$ and expanding the state $x_{\mathrm{3}}=f(y,\dot{y},\mathit{\omega })$, Eq. (18) can be rewritten as follows:

$$\begin{array}{}\text{(19)}& \left\{\begin{array}{l}{\dot{x}}_{\mathrm{1}}=x_{\mathrm{2}};\\ {\dot{x}}_{\mathrm{2}}=x_{\mathrm{3}}+b_{\mathrm{0}}u,y=x_{\mathrm{1}};\\ {\dot{x}}_{\mathrm{3}}=h;\end{array}\right.\end{array}$$

where x₁, x₂, and  x₃ are the system state variables, and $h=\dot{f}(y,\dot{y},\mathit{\omega })$.

A linear expansion state observer (LESO) is built for Eq. (19):

$$\begin{array}{}\text{(20)}& \begin{array}{l}{\dot{z}}_{\mathrm{1}}=z_{\mathrm{2}}-{\mathit{\beta }}_{\mathrm{1}}\left(z_{\mathrm{1}}-y\right),\\ {\dot{z}}_{\mathrm{2}}=z_{\mathrm{3}}-{\mathit{\beta }}_{\mathrm{2}}\left(z_{\mathrm{1}}-y\right)+b_{\mathrm{0}}u,\\ {\dot{z}}_{\mathrm{3}}=-{\mathit{\beta }}_{\mathrm{3}}\left(z_{\mathrm{1}}-y\right).\end{array}\end{array}$$

By choosing a suitable observer gain β₁,β₂, and β₃, the LESO is able to realize the real-time tracking of each state variable in the system. Taking $u=(-z_{\mathrm{3}}+u_{\mathrm{0}})/b_{\mathrm{0}}$ and neglecting the estimation error of z₃ on $f(y,\dot{y},\mathit{\omega })$, the system can be reduced to a double-integrated series structure:

$$\begin{array}{}\text{(21)}& \ddot{y}=\left[f\left(y,\dot{y},\mathit{\omega }\right)-z_{\mathrm{3}}\right]+u_{\mathrm{0}}\approx u_{\mathrm{0}}.\end{array}$$

We design the PD controller as follows:

$$\begin{array}{}\text{(22)}& u_{\mathrm{0}}=k_{\mathrm{p}}\left(r-z_{\mathrm{1}}\right)-k_{\mathrm{d}}{z}_{\mathrm{2}},\end{array}$$

where r is the reference signal, and k_p and k_d are the controller gains; according to Eqs. (21) and (22), the system's closed-loop transfer function can be obtained:

$$\begin{array}{}\text{(23)}& G_{\mathrm{cl}}=k_{\mathrm{p}}/\left(s^{\mathrm{2}}+k_{\mathrm{d}}s+k_{\mathrm{p}}\right).\end{array}$$

The characteristic equation of LESO is subsequently identified as follows:

$$\begin{array}{}\text{(24)}& \mathit{\lambda }\left(s\right)=s^{\mathrm{3}}+{\mathit{\beta }}_{\mathrm{1}}{s}^{\mathrm{2}}+{\mathit{\beta }}_{\mathrm{2}}s+{\mathit{\beta }}_{\mathrm{3}}.\end{array}$$

We select the ideal characteristic equation as follows:

$$\begin{array}{}\text{(25)}& \mathit{\lambda }\left(s\right)={\left(s+{\mathit{\omega }}_{\mathrm{0}}\right)}^{\mathrm{3}},\end{array}$$

where ω₀ is the observer bandwidth; then we have β₁=3ω₀, ${\mathit{\beta }}_{\mathrm{2}}=\mathrm{3}{\mathit{\omega }}_{\mathrm{0}}^{\mathrm{2}}$, and ${\mathit{\beta }}_{\mathrm{3}}={\mathit{\omega }}_{\mathrm{0}}^{\mathrm{3}}$.

As stated in the literature (Jin et al., 2020; Zhong et al., 2022), the parameters $k_{\mathrm{p}}={\mathit{\omega }}_{\mathrm{c}}^{\mathrm{2}}$, k_d=2ξω_c, and ω_c are selected as the controller bandwidth, and ξ is the damping ratio. This configuration simplifies the ADRC controller parameter selection process, as only the observer bandwidth ω₀ and the controller bandwidth ω_c need to be chosen. The configuration of reasonable b₀, ω₀, and ω_c in the ADRC controller can achieve the desired performance of the dynamic control system for seat vibration damping.

4.3 Control system design

Given that the MR damper is a semi-active device, it is most appropriate to utilize a semi-active control system. The semi-active control system comprises two principal components: a system controller and a damper controller. The system controller generates the desired damping force based on the dynamic response of the suspension, while the damper controller adjusts the input current in order to track the desired damping force. The semi-active controller described in this paper comprises an ADRC controller (system controller) and an ANFIS inverse model (damper controller). To simulate the vibration of the seat suspension in a tractor cab, a structure for a vehicle suspension semi-active controller with an MR damper is constructed, as illustrated in Fig. 11. The control principle is as follows: the ADRC controller calculates the active control force based on the measured output, generates the desired damping force through a force limiter based on the active control force, and then uses the ANFIS inverse model of the MR damper to adjust the command current according to the desired damping force and vehicle suspension response. Ultimately, the requisite input current is determined through the ANFIS inverse model, with the MR damper then approximating the necessary damping force.

5 Simulation and analysis

In order to ascertain the efficacy of the MR damper in seat suspension, the ADRC system was subjected to comparison with both passive control and PID control. The detailed parameter settings of the PID controller and ADRC controller are shown in Table 3.

The simulation conditions assumed a tractor travelling at a speed of 10 m s⁻¹ on an ISO C-class road profile, with road excitation taking into account the random characteristics of road conditions. In accordance with the prescribed expression for the power spectral density of the pavement outlined in GB/T 7031-2005, the mathematical model constructed in Simulink is illustrated in Eq. (26).

$$\begin{array}{}\text{(26)}& \dot{z}\left(t\right)=-\mathrm{2}\mathit{\pi }{f}_{\mathrm{0}}z\left(t\right)+\mathrm{2}\mathit{\pi }\sqrt{G_{\mathrm{0}}v}\mathit{\omega }\left(t\right),\end{array}$$

where ω(t) represents the Gaussian white noise signal, G₀ represents the pavement unevenness coefficient, and f₀ represents the cutoff frequency (0.01 Hz). For the ISO C-class road profile, take $G_{\mathrm{0}}=\mathrm{2.56}\times {\mathrm{10}}^{-\mathrm{4}}$; then the road profile simulation under v=10 m s⁻¹ obtains the input signal as shown in Fig. 12.

Table 4 enumerates the parameters of the seat suspension system, including the springs and dampers, along with the corresponding simulation results.

Figures 12–16 illustrate the random road profile; damping force; input current of the MR damper; and the system responses of the passive, PID control, and ADRC suspension systems under random excitation. Figure 12 illustrates the displacement of the road profile. Figure 13 illustrates the input current of the MR damper in the semi-active suspension system. Figure 14 illustrates the comparison between the actual and desired damping forces. Figures 15 and 16 illustrate that, in comparison to the passive suspension system and the PID control system, the ADRC suspension system exhibits a relatively lower vertical acceleration and seat dynamic disturbance. This evidence substantiates the efficacy of the MR damper with an ADRC controller in mitigating vibration within the suspension system. As illustrated in Fig. 14, the actual damping force generated by the MR damper exhibits a high degree of correlation with the desired damping force. This provides further confirmation of the efficacy of the ANFIS inverse model of the MR damper in controlling the damping force.

In order to gain a more comprehensive understanding of the performance advantages of ADRC in seat suspension, Table 5 lists the rms values of the relevant suspension performance indicators. As shown in Table 5, the rms value of the vertical vibration acceleration of the ADRC-controlled suspension system exhibited a notable decline of 68.9 % and 34.4 % in comparison to PID control and passive control, respectively. Additionally, the rms value of the dynamic disturbance of the ADRC-controlled suspension system demonstrated a considerable reduction of 50.0 % and 28.6 % in relation to PID control and passive control, respectively, which signifies a substantial vibration mitigation effect.

6 Conclusions

This paper proposes a novel semi-active controller that integrates ANFIS technology and ADRC, specifically designed for MR seat dampers. The article first establishes a forward dynamic model of the damper using the Bouc–Wen model based on experimental data. Then, it develops an inverse model using ANFIS technology and details the structure and training process of this inverse model. By combining the ANFIS inverse model with the designed ADRC controller, a new control scheme is proposed and applied to the seat suspension model. Numerical simulations verify the performance of the controller. The results show that this technology excels in precisely controlling the damping force, achieving continuous adjustability of the seat suspension damping, and effectively enhancing the vibration isolation effect of the seat suspension system. However, the computational complexity of ANFIS is relatively high, and the training time is lengthy. Furthermore, the combination of ANFIS and ADRC to determine the appropriate model structure parameters, such as the number and type of input affiliation functions and the number of fuzzy rules, necessitates specialized experience and experimentation. The different choices may result in a significant discrepancy in the model performance. To further enhance the performance of the control system, alternative optimization algorithms may be integrated with the ADRC algorithm.