Adaptive disturbance rejection control of magnetorheological seat suspension systems based on adaptive neuro-fuzzy inference system inverse modeling and improved crow search optimization

Operators of agricultural vehicles are frequently exposed to prolonged vibrations, which pose serious health risks and degrade operational performance. Magnetorheological (MR) dampers have been increasingly adopted in seat suspension systems owing to their low power consumption and fast response. However, their inherent nonlinear hysteresis and time-varying dynamics significantly challenge the effectiveness of conventional control strategies. To overcome these limitations, this study proposes a hybrid ICSA-ANFIS-ADRC control framework that integrates an improved crow search algorithm (ICSA)-optimized adaptive neuro-fuzzy inference system (ANFIS) with an active disturbance rejection control (ADRC) strategy. The nonlinear behavior of the MR damper is first modeled using an improved Bouc–Wen approach. An inverse model is then constructed through ICSA-ANFIS training to predict the control current accurately. This inverse model is subsequently embedded within the ADRC framework to enable real-time multi-modal damping force regulation. Numerical simulations based on a 3-degrees-of-freedom seat suspension model demonstrate that the proposed method significantly outperforms conventional ANFIS-ADRC and CSA-ANFIS-ADRC controllers. Specifically, up to 32.9 % reduction in vertical vibration acceleration is achieved, while robust control performance is maintained under both random and shock road conditions. The inverse model attains a root mean square prediction error below 0.15 for control current, verifying its accuracy and adaptability. The proposed ICSA-ANFIS-ADRC control scheme thus provides a promising solution for intelligent seat suspension systems, effectively mitigating low-frequency resonance and enhancing ride comfort in agricultural vehicles.

1 Introduction

Operators of agricultural vehicles, such as tractors and specialized machinery, are frequently subjected to sustained vibrations caused by uneven terrain. Prolonged exposure to such vibrations has been shown to induce adverse health effects, including fatigue, reduced alertness, impaired driving performance, and long-term musculoskeletal disorders such as lower back pain and spinal discomfort. These health risks underscore the importance of effective seat suspension systems to improve both safety and ride comfort (Lecocq et al., 2022; Schneider et al., 2023; de la Hoz-Torres et al., 2021; Barač et al., 2025). In response to these risks, multiple seat suspension technologies have been proposed, encompassing passive (Samaroo et al., 2025; Stein et al., 2009), semi-active (Zhao et al., 2023; Liu et al., 2019), and active suspension systems (Maciejewski et al., 2023; Gheibollahi et al., 2024; Zhao et al., 2025). Despite their effectiveness, active suspensions are rarely adopted in agricultural vehicles owing to their complexity and high cost, which conflict with the sector's cost-sensitive design constraints. As a result, passive and semi-active seat suspension systems are more commonly employed (Zhao et al., 2025). Passive systems, while historically dominant due to their simplicity and reliability, lack adaptability to variable road conditions and complex excitations, limiting their effectiveness in improving ride comfort. Semi-active suspensions have thus emerged as a viable compromise, combining structural simplicity with real-time adaptability (Ferhath and Kasi, 2024; Lee and Oh, 2025).

As core elements in semi-active control systems, magnetorheological (MR) dampers exhibit strong potential for attenuating seat vibrations and mitigating whole-body vibrational exposure in vehicle occupants. With growing demand for enhanced seating comfort, the role of seat suspension systems in ensuring both ride comfort and safety has garnered substantial research interest (Hua, 2025; Gad et al., 2025). These systems utilize intelligent control strategies to dynamically modulate damping force in response to varying road and operational conditions, thereby enhancing ride comfort and handling stability without the high energy demands of active suspension systems. However, under shear-mode operation, MR dampers exhibit asymmetric hysteresis in their force–velocity characteristics, which introduces nonlinearities and causes time lag effects when conventional control methods are applied. Although methods such as active disturbance rejection control (ADRC) can partially address external disturbances and time delays, their reliance on fixed controller parameters limits adaptability under varying conditions (Wang et al., 2024; Yaghoubi and Ghanbarzadeh, 2024). Consequently, real-time optimization of key controller parameters has become a pivotal challenge in improving the damping performance of MR-based suspension systems.

Recent advancements have explored two primary directions for integrating intelligent optimization algorithms into control systems: (1) improving swarm intelligence optimization algorithms (e.g., particle swarm optimization, PSO) to prevent control performance degradation and (2) modifying conventional controllers to better accommodate the nonlinear and time-varying characteristics of MR dampers. For instance, Guo et al. (2024) developed an enhanced Jiles–Atherton hysteresis model optimized by an improved PSO algorithm, which significantly increased parameter identification accuracy (error was reduced by 40.2 %) and improved the prediction of damping force. Similarly, Zhang et al. (2024) applied NSGA-II to optimize the structural parameters of a conical-channel MR damper, demonstrating improved yaw stability under extreme maneuvers. Li et al. (2024) proposed a sliding-mode control strategy enhanced by an extended African vulture optimization algorithm (EAVOA) to suppress vibrations in semi-active suspension systems under stochastic road conditions. Population initialization was optimized via logistic chaotic mapping, and a nonlinear disturbance-based hunger value adjustment strategy was adopted to efficiently optimize control parameters. Road tests on B- and D-grade surfaces showed a maximum reduction of 42.5 % in vehicle body acceleration and a 25.9 % decrease in suspension dynamic travel, with smooth damping force response, providing a new approach to precise control of MR dampers. Guo et al. (2025) introduced a sensitivity analysis method based on PSO to identify parameters in MR damper dynamic models. The effects of particle number (N), iteration count (T), and learning factors ($c_{\mathrm{1}}/c_{\mathrm{2}}$) on the fitting accuracy of the hyperbolic tangent model were systematically quantified, and the optimal parameter set was determined through orthogonal experiments. As a result, the root mean square error (RMSE) of damping force prediction was reduced to 0.013 under the 0.8 A condition. Experimental results showed strong agreement between simulation and measured data, establishing a high-accuracy modeling foundation for semi-active suspension control. Zhao et al. (2024) proposed an adaptive neural-network-based backstepping control strategy to address uncertainties in sprung mass and time-varying input delays in MR semi-active air suspension systems. The nonlinear dynamics of the air spring were precisely modeled using a radial basis function neural network (RBFNN), and the sprung mass was continuously updated by a projection-constrained estimator. A Lyapunov–Krasovskii functional-based delay compensator was designed to suppress a 32.4 ms actuator delay. Under conditions involving step changes in sprung mass and variable delay, the vehicle body acceleration was reduced by 16.2 %, and the suspension dynamic deflection was reduced by 5.9 %, offering a viable solution for high-comfort suspension systems. Ji et al. (2024) proposed a variable universe fuzzy Proportional–Integral–Derivative (PID) control strategy based on a dynamic adjustment function (VUFP-DAF) to address the limitations of expert-dependent fuzzy rules and poor adaptability of fixed scaling factors in Continuous Damping Control (CDC) damper-based semi-active suspension systems. A dynamic scaling factor, using system error e(t) and its rate of change ec(t) as inputs, was designed for real-time universe adaptation, and system stability was rigorously proven. Verification on B- and C-grade random roads and bump surfaces showed a maximum reduction of 53.88 % in vertical acceleration and 45.65 % in suspension deflection, significantly outperforming conventional fuzzy PID and fixed-factor methods, offering a robust control solution under complex conditions. Zhan et al. (2025) addressed the adaptability challenges of MR semi-active suspension systems under strongly nonlinear and time-varying conditions by proposing a proximal policy optimization algorithm integrated with neural ordinary differential equations (PPO-NODEs). The improved Bouc–Wen model was used to accurately characterize damper dynamics, and an actor network enhanced by NODEs was designed to model continuous suspension state evolution. Simulations on C-grade random roads showed a 30.22 % reduction in vehicle acceleration compared to passive suspensions, while bench experiments confirmed an acceleration reduction of up to 77.51 %. The algorithm substantially improved control generalization under complex conditions, offering a new solution for intelligent vehicle suspension systems.

Despite these advances, most existing studies focus on refining controller structures or incorporating optimization algorithms for offline parameter tuning, whereas few works have addressed the real-time optimization of controller outputs in conjunction with adaptive control strategies. This gap limits the applicability of existing methods under complex and rapidly changing road conditions typically encountered by agricultural vehicles. To address these challenges, the present study makes the following contributions:

1. An improved crow search algorithm (ICSA) is developed by introducing a triangular probability distribution mechanism, which enhances population diversity and accelerates convergence to global optima.

2. A dynamic adaptive neuro-fuzzy inference system (ANFIS) structure with time-varying membership functions is constructed, enabling real-time adjustment of the damping control strategy and reducing the influence of the MR damper's time-varying properties.

3. A hybrid ICSA-ANFIS-ADRC control framework is proposed, where a Kalman filter is embedded in the ADRC observation layer to suppress noise, and the control signals are dynamically optimized by the ICSA-ANFIS inverse model. This integration achieves multi-modal damping control and robust vibration suppression across diverse operating conditions.

In summary, the proposed method not only improves modeling accuracy and control adaptability but also demonstrates superior robustness in vibration isolation compared with conventional ANFIS-ADRC and CSA-ANFIS-ADRC strategies.

Figure 1 The 3-DOF “human seat” model, with its coordinate system.

2 Dynamic model of the semi-active seat suspension

2.1 Tractor seat suspension model

Figure 1 illustrates a 3-degrees-of-freedom (3-DOF) lumped-parameter model designed to characterize the dynamic behavior of the seat suspension system equipped with an MR damper. In this model, m_b and m_s represent the masses of the human body and the seat, respectively; k_b and c_b correspond to the stiffness and damping of the human body; k_s and c_MR denote the stiffness and equivalent damping of the seat suspension and MR damper; and z_b, z_s, and z₀ represent the displacements of the human body, the seat, and the base excitation, respectively.

Based on Newtonian mechanics, the dynamic equations governing the 3-DOF seat suspension system are formulated as follows:

$$\begin{array}{}\text{(1)}& \left\{\begin{array}{l}{\mathit{m}}_{\mathrm{b}}{\ddot{z}}_{\mathrm{b}}+{\mathit{k}}_{\mathrm{b}}({\mathit{z}}_{\mathrm{1}}-{\mathit{z}}_{\mathrm{s}})+{\mathit{c}}_{\mathrm{b}}({\dot{z}}_{\mathrm{b}}-{\dot{z}}_{\mathrm{s}})=\mathrm{0}\\ {\mathit{m}}_{\mathrm{s}}{\ddot{z}}_{\mathrm{s}}-{\mathit{k}}_{\mathrm{b}}({\mathit{z}}_{\mathrm{1}}-{\mathit{z}}_{\mathrm{s}})-{\mathit{c}}_{\mathrm{b}}({\dot{z}}_{\mathrm{b}}-{\dot{z}}_{\mathrm{s}})\\ +{\mathit{k}}_{\mathrm{s}}({\mathit{z}}_{\mathrm{s}}-{\mathit{z}}_{\mathrm{c}})+{\mathbf{F}}_{\mathrm{MR}}=\mathrm{0}\\ {\mathit{m}}_{\mathrm{c}}{\ddot{z}}_{\mathrm{c}}+{\mathit{k}}_{\mathrm{c}}({\mathit{z}}_{\mathrm{c}}-{\mathit{z}}_{\mathrm{0}})+{\mathit{c}}_{\mathrm{c}}({\dot{z}}_{\mathrm{c}}-{\dot{z}}_{\mathrm{0}})\\ -{\mathit{k}}_{\mathrm{s}}({\mathit{z}}_{\mathrm{s}}-{\mathit{z}}_{\mathrm{c}})-{\mathbf{F}}_{\mathrm{MR}}=\mathrm{0}\end{array},\right.\end{array}$$

where F_MR represents the damping force of the MR damper of the seat.

To facilitate controller design and system simulation, the equations of motion are converted into a state space representation, as shown in Eq. (2):

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

where the state variable $\mathbf{X}=[{\mathit{z}}_{\mathrm{b}},{\mathit{z}}_{\mathrm{s}},{\mathit{z}}_{\mathrm{c}},{\dot{z}}_{\mathrm{1}},{\dot{z}}_{\mathrm{s}},{\dot{z}}_{\mathrm{c}}{]}^T$, the input variable $\mathbf{U}=[{\mathit{z}}_{\mathrm{0}},{\dot{z}}_{\mathrm{0}},{\mathbf{F}}_{\mathrm{MR}}{]}^{\mathbf{T}}$, and the output variable $\mathbf{Y}=[{\ddot{z}}_{\mathrm{b}},{\dot{z}}_{\mathrm{b}},{\mathit{z}}_{\mathrm{b}},{\ddot{z}}_{\mathrm{s}},{\dot{z}}_{\mathrm{s}},{\mathit{z}}_{\mathrm{s}}{]}^{\mathbf{T}}$. The equation of state of the system after substitution is given by Eq. (3).

$$\begin{array}{}\text{(3)}& \begin{array}{rl}& \dot{X}=\left[\begin{array}{ll}\begin{array}{lll}\mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}& \begin{array}{lll}\mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\\ \begin{array}{ccc}-\frac{k_{\mathrm{b}}}{m_{\mathrm{b}}}& \frac{k_{\mathrm{b}}}{m_{\mathrm{b}}}& \mathrm{0}\\ \frac{k_{\mathrm{b}}}{m_{\mathrm{s}}}& \frac{k_{\mathrm{s}}-k_{\mathrm{b}}}{m_{\mathrm{s}}}& \frac{k_{\mathrm{s}}}{m_{\mathrm{s}}}\\ \mathrm{0}& \frac{k_{\mathrm{s}}}{m_{\mathrm{c}}}& -\frac{k_{\mathrm{c}}+k_{\mathrm{s}}}{m_{\mathrm{c}}}\end{array}& \begin{array}{ccc}-\frac{c_{\mathrm{b}}}{m_{\mathrm{b}}}& \frac{c_{\mathrm{b}}}{m_{\mathrm{b}}}& \mathrm{0}\\ \frac{c_{\mathrm{b}}}{m_{\mathrm{s}}}& -\frac{c_{\mathrm{b}}}{m_{\mathrm{s}}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& -\frac{c_{\mathrm{c}}}{m_{\mathrm{c}}}\end{array}\end{array}\right]\\ & \left[\begin{array}{c}\begin{array}{c}{z}_{\mathrm{b}}\\ z_{\mathrm{s}}\\ z_{\mathrm{c}}\end{array}\\ \begin{array}{c}{\dot{z}}_{\mathrm{b}}\\ {\dot{z}}_{\mathrm{s}}\\ {\dot{z}}_{\mathrm{c}}\end{array}\end{array}\right]+\left[\begin{array}{c}\begin{array}{ccc}\mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\\ \begin{array}{ccc}\mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& -\frac{\mathrm{1}}{m_{\mathrm{s}}}\\ \frac{k_{\mathrm{c}}}{m_{\mathrm{c}}}& \frac{c_{\mathrm{c}}}{m_{\mathrm{c}}}& \frac{\mathrm{1}}{m_{\mathrm{c}}}\end{array}\end{array}\right]\left[\begin{array}{c}{z}_{\mathrm{0}}\\ {\dot{z}}_{\mathrm{0}}\\ F_{\mathrm{MR}}\end{array}\right]\\ & \mathbf{Y}=\left[\begin{array}{cc}\begin{array}{ccc}-\frac{{\mathit{k}}_{\mathrm{b}}}{{\mathit{m}}_{\mathrm{b}}}& \frac{{\mathit{k}}_{\mathrm{b}}}{{\mathit{m}}_{\mathrm{b}}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathbf{1}& \mathrm{0}& \mathrm{0}\end{array}& \begin{array}{ccc}-\frac{{\mathit{c}}_{\mathrm{b}}}{{\mathit{m}}_{\mathrm{b}}}& \frac{{\mathit{c}}_{\mathrm{b}}}{{\mathit{m}}_{\mathrm{b}}}& \mathrm{0}\\ \mathbf{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\\ \begin{array}{ccc}\frac{{\mathit{k}}_{\mathrm{b}}}{{\mathit{m}}_{\mathrm{s}}}& \frac{{\mathit{k}}_{\mathrm{b}}-{\mathit{k}}_{\mathrm{s}}}{{\mathit{m}}_{\mathrm{s}}}& \frac{{\mathit{k}}_{\mathrm{s}}}{{\mathit{m}}_{\mathrm{s}}}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}\end{array}& \begin{array}{ccc}\frac{{\mathit{c}}_{\mathrm{b}}}{{\mathit{m}}_{\mathrm{s}}}& -\frac{{\mathit{c}}_{\mathrm{b}}}{{\mathit{m}}_{\mathrm{s}}}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\end{array}\right]\\ & \left[\begin{array}{c}\begin{array}{c}{\mathit{z}}_{\mathrm{b}}\\ {\mathit{z}}_{\mathrm{s}}\\ {\mathit{z}}_{\mathrm{c}}\end{array}\\ \begin{array}{c}{\dot{z}}_{\mathrm{b}}\\ {\dot{z}}_{\mathrm{s}}\\ {\dot{z}}_{\mathrm{c}}\end{array}\end{array}\right]+\left[\begin{array}{c}\begin{array}{ccc}\mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\\ \begin{array}{ccc}\mathrm{0}& \mathrm{0}& -\frac{\mathrm{1}}{{\mathit{m}}_{\mathrm{s}}}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\end{array}\right]\left[\begin{array}{c}{\mathit{z}}_{\mathrm{0}}\\ {\dot{z}}_{\mathrm{0}}\\ {\mathbf{F}}_{\mathrm{MR}}\end{array}\right]\end{array}\end{array}$$

2.2 Improved Bouc–Wen model for MR dampers

Accurate mathematical modeling of MR dampers is essential for ensuring realistic system dynamics simulations. The Bouc–Wen model is widely adopted for its ability to characterize the nonlinear hysteresis behavior of MR dampers, particularly the relationship between relative displacement, velocity, and damping force. Thus, it has gained widespread application. While the conventional Bouc–Wen model comprises a spring, a viscous damper, and a hysteresis element arranged in parallel, an improved version introduces an additional spring and damper to enhance modeling accuracy and flexibility. In this study, the improved Bouc–Wen model is employed to represent the forward dynamic behavior of the MR damper, capturing its nonlinear and rate-dependent response characteristics, and its structure is shown in Fig. 2.

Figure 2 Schematic of the improved Bouc–Wen model structure.

To evaluate the force–velocity and force–displacement characteristics of the MR damper, experiments were conducted using an electro-hydraulic servo fatigue testing machine (WANGCE, China), as shown in Fig. 3. The test system consisted of a servo actuator, a load cell, displacement and velocity sensors, a programmable current driver, and a data acquisition module connected to a workstation. The actuator applied sinusoidal excitations with a displacement amplitude of 15 mm at a frequency of 2 Hz, while the control current was varied from 0 to 1.5 A in increments of 0.25 A. Displacement and velocity signals of the piston rod were recorded using high-precision sensors, while the input current signal was simultaneously monitored through the driver unit. The damping force was directly measured by a calibrated load cell attached in series with the damper. All signals were sampled at a frequency of 500 Hz, filtered to remove high-frequency noise, and subsequently processed in MATLAB/Simulink for parameter identification. This procedure provided a comprehensive dataset of piston displacement, velocity, current, and force responses, which was then used to calibrate the improved Bouc–Wen model through the Simulink Design Optimization toolbox.

Figure 3 Damping characteristic test of the MR damper.

The improved Bouc–Wen model can effectively fit the relationship between the output damping force of the MR damper and the relative displacement and velocity of the piston. The damping force is expressed as follows:

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

where F is the damping force of the MR damper, c₀ and c₁ are the viscous damping coefficients; k₀ and k₁ are the stiffness coefficient; α is the scaling parameter of the hysteresis force in the damping force; x₀ is the relative displacement of the equilibrium position offset; z is the hysteresis variable; x is the total displacement of the MR damper; y is the hysteresis displacement of c₁; n is the rounding coefficient; and γ, β, and A are the parameters used to adjust the smoothness of the model and the linear characteristics of the curves from the pre-yielding region to the post-yielding region. Parameters are used to adjust the smoothness of the model and the linear character of the curve from the pre-yield to the post-yield region.

Figure 4 Simulation model of the improved Bouc–Wen model.

The improved Bouc–Wen model is built in Simulink, and the constructed model can be obtained using Eq. (4), as shown in Fig. 4.

The improved Bouc–Wen model incorporates eight unknown parameters. While conventional identification methods such as genetic algorithms and least squares exist, they often entail computationally intensive procedures. Consequently, this study employs the Simulink Design Optimization toolbox for parameter identification. Measured piston relative displacement, velocity, and current serve as model inputs, with damping force designated as the output. Model predictions are subsequently compared against experimental results.

Following the implementation of the improved Bouc–Wen model in Simulink, parameter calibration was performed utilizing experimental datasets obtained under specified operating conditions: an excitation frequency of 2 Hz, displacement amplitude of 15 mm, and input currents ranging from 0 to 1.5 A. Sensitivity analysis indicated that parameters α and c₀ exhibit strong current dependency, while the remaining six parameters remain invariant across different operating states. Through iterative optimization, the current-dependent characteristics of α and c₀ were successfully approximated by quadratic polynomial functions, demonstrating close agreement with experimental observations. The finalized parameter identification results, encompassing both current-variant and current-invariant parameters, are comprehensively summarized in Table 1.

Figure 5 Comparison between model simulations and experimental data for MR damper characteristics under 2 Hz excitation.

Table 1 Improved Bouc–Wen model parameters.

The parameters of the improved Bouc–Wen model were identified using experimental datasets obtained under the specified excitation conditions. A sensitivity analysis was first conducted to determine the dependence of each parameter on the input current. The results indicated that the parameters α and c₀ exhibited strong current dependency as they directly influence the scaling of the hysteresis force and the viscous damping characteristics of the MR damper. Consequently, both parameters were expressed as quadratic polynomial functions of the input current, enabling accurate representation of their variation across different operating states. In contrast, the remaining parameters (β, γ, A, k₀, k₁, x₀, n) showed negligible dependence on current and were therefore treated as constants. This hybrid parameterization approach ensures that the improved Bouc–Wen model captures both the current-dependent nonlinear hysteresis effects and the invariant structural characteristics of the MR damper. As a result, the model achieves a high degree of consistency with the experimental measurements, as confirmed in Fig. 5.

Figure 5 illustrates the comparison between the experimental results (solid line) and the model predictions (dashed line). From the above figure, it can be intuitively seen that the curve obtained from the parameter fitting using the improved Bouc–Wen model and the applied Simulink parameter identification is basically consistent with the curve obtained from the damper characteristic test, and the fitting effect is good.

Figure 6 Flowchart of the crow search algorithm.

3 Inverse modeling of the MR damper based on ICSA-ANFIS

3.1 Standard crow search algorithm and its improvement

The crow search algorithm (CSA), introduced by Askarzadeh in 2016 (Askarzadeh, 2016), addresses global optimization problems by mimicking the foraging behavior of crows. It is a metaheuristic method characterized by few tunable parameters and strong convergence performance. It is considered to be a novel metaheuristic algorithm, and its process is illustrated in Fig. 6.

Specifically, the detailed procedure of the CSA algorithm is as follows:

• Step 1. A population of N crows is randomly initialized, with each crow's position and memory location being set to the same value. The initial positions of these crows (i.e., the initial solution set) are denoted as X₀= (X₁, X₂, …, X_N), and the memorized best hiding locations of food by the crows are denoted as M₀= (M₁, M₂, …, M_N), where M₀=X₀.

• Step 2. The position X_i of each crow is updated. The update principle assumes that crow i randomly selects crow j to follow. Two possible outcomes arise. (1) If crow j is unaware of being followed then crow i will move toward j's memorized best food-hiding location M_j. (2) If crow j is aware of being followed, it will intentionally mislead crow i to a random position. The position update strategy is mathematically expressed as

  $$\begin{array}{}\text{(5)}& X_i^{k+\mathrm{1}}\left(l\right)=\left\{\begin{array}{l}{X}_i^k\left(l\right)+\mathrm{rand}\cdot \mathrm{FL}\cdot \left[X_i^k\left(l\right)-X_j^k\left(l\right)\right],r_k\ge \mathrm{AP}\\ \text{random\hspace{0.17em} position},r_k<\mathrm{AP}\end{array}\right.,\end{array}$$

  where k denotes the current iteration number; AP and FL represent the awareness probability and flight length, respectively; and r_k and rand refer to uniformly distributed random numbers and random sequences within the range [0,1].

• Step 3. The feasibility of the solution X_i is evaluated. If the updated solution lies within the feasible domain, the update is accepted and replaces the original; otherwise, the original solution is retained.

• Step 4. The fitness value f of the solution X_i is calculated.

• Step 5. The food-hiding position M_i is updated. The specific update strategy is expressed as

  $$\begin{array}{}\text{(6)}& M_i^{k+\mathrm{1}}=\left\{\begin{array}{c}{X}_i^{k+\mathrm{1}},f\left(X_i^{k+\mathrm{1}}\right)>f\left(M_i^k\right)\\ M_i^k,f\left(X_i^{k+\mathrm{1}}\right)\le f\left(M_i^k\right)\end{array}.\right.\end{array}$$

• Step 6. Steps 2 to 5 are repeated to iteratively optimize X_i and M_i until the convergence criterion of the algorithm is satisfied. In the standard CSA described above, the crow j being followed in step 2 (i.e., the target crow) is randomly selected. This random selection often prevents many crows from improving their positions, which may cause the algorithm to become trapped in local optima. To address this issue, the standard CSA is improved by introducing a triangular probability distribution strategy to select the target crow (Lin et al., 2021). The triangular probability distribution is expressed as

  $$\begin{array}{}\text{(7)}& p_i={\displaystyle \frac{\mathrm{2}(N+\mathrm{1}-i)}{N(N+\mathrm{1})}},(i=\mathrm{1},\mathrm{2},\mathrm{\dots }N).\end{array}$$

  In each iteration, crows with a probability greater than p_t are selected, where p_t is calculated as

  $$\begin{array}{}\text{(8)}& p_t={\displaystyle \frac{\mathrm{2}+\mathrm{2}(N-\mathrm{1})\cdot {\mathrm{rand}}_{\mathrm{2}}}{N(N+\mathrm{1})}},\end{array}$$

  where rand₂ is a uniformly distributed random number in the range [0,1].

  The crows selected through this triangular probability method are provided with a greater competitive advantage. One target crow is then randomly selected from this group for tracking, ensuring population diversity within the crow swarm. Furthermore, to enhance the diversity of solutions, a mutation operation is applied to the updated solutions obtained from Eq. (7) (Lin et al., 2021). As an example of this operation, several consecutive elements in the encoding are randomly selected and then reversed to introduce variability.

  $$\begin{array}{}\text{(9)}& \begin{array}{rl}C& =C\mathit{\#}\left[\mathrm{spml}\right]=[\mathrm{1},\mathrm{3},\mathrm{4},\mathrm{7},\mathrm{9},\mathrm{11}]\mathit{\#}[\mathrm{2},\mathrm{3}]\\ & =[\mathrm{1},\mathrm{7},\mathrm{4},\mathrm{3},\mathrm{9},\mathrm{11}]\end{array}\end{array}$$

  In the above, C denotes the code segment to be mutated, # represents the inversion mutation operation, ml is the randomly selected mutation length, and sp is the randomly selected starting position of the mutation.

Figure 7 Structure of ANFIS.

3.2 ANFIS-based inverse modeling of the MR damper using ICSA

In constructing the ANFIS-based inverse model of the MR damper, the input dimensionality must be carefully selected to avoid the curse of dimensionality. Accordingly, a four-input, single-output architecture is employed (Fig. 7).

First, the damping force F(k) is calculated using the forward model described in Sect. 2.2 to calculate the damping force F(k) from three inputs: displacement x(k), velocity $\dot{x\left(k\right)}$, and control current I(k). This generates a training dataset comprising x(k), $\dot{x\left(k\right)}$, F(k); the previous current I(k−1); and the present current I(k), of which the first four are used as inputs for training the ANFIS network. The ANFIS structure is optimized using the ICSA (the specific optimization procedure is detailed later in this subsection), enabling the network to accurately output the predicted control current $\widehat{I}\left(k\right)$. Once the optimal inverse model of the MR damper has been trained, the forward model is applied again to compute the predicted damping force $\widehat{F}\left(k\right)$, thereby enabling further validation of the inverse model's effectiveness, as shown in Fig. 8.

Figure 8 Schematic diagram of the ICSA-ANFIS-based inverse modeling strategy.

For the ANFIS structure illustrated in Fig. 7, the system output y is given by

$$\begin{array}{}\text{(10)}& y={\sum }_{i=\mathrm{1}}^{\mathrm{4}}{}_i{f}_i=\mathbf{A}\mathit{\theta },\end{array}$$

where θ is the parameter vector, and A is the coefficient matrix constructed from the input data.

If the four inputs of the ANFIS model are denoted as x₁, x₂, x₃, and x₄ then the output of the fourth layer is given by

$$\begin{array}{}\text{(11)}& f_i=p_i{x}_{\mathrm{1}}+q_i{x}_{\mathrm{2}}+r_i{x}_{\mathrm{3}}+s_i{x}_{\mathrm{4}}+t_i.\end{array}$$

In Eq. (17), $\mathit{\theta }=\mathit{\{}{p}_i,q_i,r_i,s_i,t_i\mathit{\}}$ represents a set of consequent parameters, organized as a 20×1 dimensional vector. If N_t training samples are provided, the output vector y and the coefficient matrix A have dimensions of N_t×1 and N_t×20, respectively.

The standard hybrid learning algorithm of ANFIS consists of a combination of the gradient descent method and the least squares estimation (LSE) technique. During each iteration, the system input is first forward propagated using the least squares method to estimate the consequent parameters θ_e:

$$\begin{array}{}\text{(12)}& {\mathit{\theta }}_e=(A^{\mathrm{T}}A{)}^{-\mathrm{1}}{A}^{\mathrm{T}}y.\end{array}$$

Accordingly, the predicted output of the ANFIS, denoted as y_e, is given by

$$\begin{array}{}\text{(13)}& y_e=A{\mathit{\theta }}_e.\end{array}$$

Figure 9 Flowchart of the ICSA-ANFIS inverse modeling process.

Figure 10 Schematic diagram of ICSA-ANFIS-ADRC MR damper controller.

While keeping the estimated consequent parameters θ_e fixed, the antecedent parameters are optimized using a backpropagation algorithm based on gradient descent in order to derive the solution for the ICSA algorithm. Considering the tendency of the gradient descent method to fall into local optima, the ANFIS model is enhanced using the ICSA approach in this study. Specifically, ICSA is employed to optimize the antecedent parameters, while the least squares method is applied to identify the consequent parameters. The inverse modeling process is summarized in Fig. 9, which outlines the key stages: input selection, membership function definition, ICSA-based optimization, and model validation using forward simulation.

1. Model structure and data preparation.

   • Input selection. The inputs are selected to be displacement x(k), velocity $\dot{x\left(k\right)}$, damping force F(k), and the previous control current I(k−1). The predicted current is used as the output $\widehat{I}\left(k\right)$, forming a four-input, single-output ANFIS architecture.

   • Membership function selection. Gaussian membership functions are selected for each input variable. The general form of the Gaussian membership function is expressed as

     $$\begin{array}{}\text{(14)}& {\mathit{\mu }}_A\left(x\right)=\exp \left(-{\displaystyle \frac{(x-c_i{)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_i^{\mathrm{2}}}}\right),\end{array}$$

     where c_i denotes the center of the membership function, and σ_i represents its standard deviation. Three Gaussian membership functions are assigned to each input in order to cover the entire range of the input variables. The total number of antecedent parameters is 4 inputs ×3 functions ×2 parameters =24 antecedent parameters. The consequent parameters are defined as a linear combination of $(p_j{x}_{\mathrm{1}}+q_j{x}_{\mathrm{2}}+r_j{x}_{\mathrm{3}}+s_j{x}_{\mathrm{4}}+t_j)$, resulting in a total of 64 consequent parameters.

2. ICSA-ANFIS optimization.

   • Parameter encoding and initialization. The antecedent parameters (c_i, σ_i) and consequent parameters (p_i, q_i, r_i, s_i, t_i) are concatenated into a real-valued encoded vector. The ICSA parameters are initialized as population size N=50, maximum number of iterations K=500, awareness probability AP =0.1, and flight length FL =2.0. The boundaries for Gaussian membership function parameters are set as c_i∈ [−30,30], σ_i∈[0.1,10] according to the range of the input signals.

   • Fitness function design. The fitness function is defined by minimizing the root mean square error (RMSE) of the predicted control current:

     $$\begin{array}{}\text{(15)}& \begin{array}{rl}& f={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+\mathrm{RMSE}}},\\ & \mathrm{RMSE}=\sqrt{{\displaystyle \frac{\mathrm{1}}{N_{\mathrm{t}}}}{\sum }_{k=\mathrm{1}}^{N_{\mathrm{t}}}\left(\right(i\left(k\right)-\left(k\right){)}^{\mathrm{2}}},\end{array}\end{array}$$

     where N_t denotes the number of training samples.

   • Iterative optimization (per iteration loop).

     • Step 1. For each crow individual, the antecedent parameters are held constant, and the consequent parameters are estimated using the least squares method (as in Eq. 18). The fitness value f(X_i) is then computed.

     • Step 2: update the crow's position. A target crow is selected using the triangular probability distribution strategy (Eqs. 7 and 8) to avoid premature convergence caused by random selection. The position is updated according to Eq. (5): $X_i^{k+\mathrm{1}}=X_i^k+\mathrm{FL}\times \mathrm{rand}\times (M_j^k-X_i^k)$.

     • Step 3: Gaussian parameter diversity enhancement. A randomly selected segment of (c_i, σ_i) in the solution vector is subjected to inversion mutation (as in Eq. 9); i.e., the selected subsequence is reversed.

     • Step 4: boundary handling and memory update. Out-of-bound parameters are corrected to fall within the predefined limits. If the new solution is superior to the memorized one, the memory location is updated accordingly (as in Eq. 6).

     • Step 5. Record the global optimal solution (X_iθ_e), i.e., the antecedent and consequent parameters corresponding to the highest fitness value.

3. Advantageous Handling of Gaussian Membership Functions.

   • Targeted mutation. Inversion mutation is preferentially applied to the (c_i, σ_i) pairs of low-fitness individuals in order to prevent excessive overlap among Gaussian membership functions.

   • Convergence criterion. A Gaussian function similarity check is introduced as part of the convergence criterion: if the maximum variation in c_i and σ_i over 20 consecutive generations of the global best solution satisfies max|Δc_i|<01 and max|Δσ_i|<005, the optimization is terminated early.

Figure 11 Best convergence curves of training error for different inverse modeling algorithms.

Figure 12 Command current predicted by the ICSA-ANFIS inverse model of the MR damper for the training dataset.

Figure 13 Damping force predicted by the ICSA-ANFIS inverse model of the MR damper for the training dataset.

4 Control strategy and controller design

4.1 ADRC controller design

Based on the state space representation in Eq. (3), the system's output equation is derived as follows:

$$\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, where k₁ and c₁ are fixed, and Eq. (24) 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 we introduce the state variables x₁=y, $x_{\mathrm{2}}=\dot{y}$ and expand the state $x_{\mathrm{3}}=f\left(y\dot{y}\mathit{\omega }\right)$. Subsequently, Eq. (26) 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₃ refer to the system state variable, and $h=\dot{f}\left(y\dot{y}\mathit{\omega }\right)$ is the linearly expanding state observer established by Eq. (25).

$$\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 β₁β₂,β₃, the linear extended state observer (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\left(y\dot{y}\mathit{\omega }\right)$, 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}$$

The PD controller is then designed as follows:

$$\begin{array}{}\text{(22)}& u_{\mathrm{0}}=k_p\left(r-z_{\mathrm{1}}\right)-k_d{z}_{\mathrm{2}}.\end{array}$$

Through Eqs. (27) and (28), in which r is the reference signal and k_p and k_d denote the controller gain, the system closed-loop transfer function can be obtained.

$$\begin{array}{}\text{(23)}& G_{cl}=k_p/\left(s^{\mathrm{2}}+k_{d}s+k_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}$$

The ideal characteristic equation can then be selected:

$$\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}}$.

Table 2 Prediction error of training data for three types of inverse models.

Figure 14 Command current predicted by the ICSA-ANFIS inverse model of the MR damper for the validation dataset.

As stated in the literature (Jin et al., 2020; Zhong et al., 2022), the selection of ADRC parameters, particularly the observer bandwidth (ω₀) and controller bandwidth (ω_c), plays a critical role in determining system performance. In practice, ω₀ is chosen to be sufficiently larger than ω_c to ensure that the linear extended state observer (LESO) can accurately estimate disturbances and unmodeled dynamics in real time while avoiding excessive amplification of measurement noise. Typically, ω₀ is set to be 3–5 times greater than ω_c, which provides a balanced trade-off between estimation accuracy and noise sensitivity. The controller bandwidth ω_c is selected according to the desired closed-loop dynamics and can be tuned using the equivalent second-order model parameters ($k_p={\mathit{\omega }}_{\mathrm{c}}^{\mathrm{2}}$, k_d=2ξω_c), where ξ denotes the damping ratio. A sensitivity analysis revealed that overly small ω₀ values reduce disturbance rejection capability, while excessively large values amplify sensor noise. Similarly, a higher ω_c improves tracking performance but may lead to control signal saturation. Therefore, the final parameter set was determined by balancing these effects, ensuring robust vibration suppression under both random and shock excitations.

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

Figure 15 Damping force predicted by the ICSA-ANFIS inverse model of the MR damper for the validation dataset.

Table 4 Simulation comparison of seat suspension control performance.

4.2 Control system design

Given the inherent semi-active nature of MR dampers, this study employs a hierarchical control architecture comprising two cascaded modules:

• ADRC as the system controller, which computes the reference damping force F_d based on real-time suspension dynamics (displacement z_s, velocity ${\dot{z}}_{\mathrm{s}}$, and acceleration ${\ddot{z}}_{\mathrm{s}}$);

• ICSA-ANFIS inverse model as the damper controller, which maps F_d and suspension states to the optimal command current I_cmd, addressing the nonlinear hysteresis of the MR damper.

As illustrated in Fig. 10, the semi-active control framework operates through the following sequence:

• The ADRC generates an active control force F_c via its extended state observer (ESO), compensating for unmodeled disturbances (e.g., road excitation z₀).

  $$\begin{array}{}\text{(26)}& F_{\mathrm{c}}=k_p\left(r-z_{\mathrm{s}}\right)+k_d\left(\dot{r}-{\dot{z}}_{\mathrm{s}}\right)-{\displaystyle \frac{z_{\mathrm{3}}}{b_{\mathrm{0}}}}\end{array}$$

  In the above, z₃ estimates total disturbances, and k_p and k_d are PD gains.

• A force limiter constrains F_c to the MR damper's achievable force range F_d∈[F_minF_max], ensuring physical realizability.

• The improved ANFIS inverse model resolves the nonlinear mapping F_d→I_cmd using a hybrid learning algorithm, incorporating displacement z_s, velocity ${\dot{z}}_{\mathrm{s}}$, and historical force data F_d(k−1).

• The MR damper modulates its rheological properties based on I_cmd, generating the approximate damping force F_MR≈F_d.

This dual-loop structure synergizes ADRC's robust disturbance rejection (bandwidth ω_c=50 rad s⁻¹) with ICSA-ANFIS's nonlinear hysteresis compensation (rms current error <0.015 A), achieving 92 % vibration attenuation in tractor cab simulations under ISO 2631-1 vertical vibration criteria (International Organization for Standardization, 2025).

5 Simulation and analysis

5.1 Simulation analysis of MR damper inverse modeling

The displacement and current signals used in the training dataset (as shown in Fig. 5) were generated using band-limited white-noise signals with amplitude ranges of −20 to 20 mm and 0 to 3 A and frequency ranges of 0–7 and 0–6 Hz, respectively. The data were collected over a duration of 20 s at a sampling frequency of 500 Hz. The odd-indexed samples were used as the training dataset, while the even-indexed samples were used for validation. To evaluate the performance of the ICSA-ANFIS model, inverse models of the MR damper were also developed using standard ANFIS and CSA-ANFIS. Each of the three optimization algorithms was independently executed 20 times. The best convergence curves of the training errors for the three modeling algorithms are illustrated in Fig. 11. It can be observed that the ICSA-ANFIS algorithm effectively overcomes the tendency of standard ANFIS to become trapped in local optima and exhibits significantly superior convergence accuracy compared to the other two optimization methods.

Table 2 summarizes the minimum prediction errors of the control current and the corresponding damping force obtained using different inverse modeling approaches. As shown in Table 2, the ICSA-ANFIS model achieves the highest prediction accuracy for both current and damping force.

The ICSA-ANFIS model was validated using the test dataset. Figure 12 shows that the predicted current accurately tracks the target current, while Fig. 13 demonstrates that the predicted damping force, obtained from the ICSA-ANFIS inverse model, closely matches the target damping force.

Figure 16 Seat vertical acceleration response under ISO class-D road excitation at a vehicle speed of 1 m s⁻¹.

Figure 17 Seat vertical acceleration response under ISO class-D road excitation at a vehicle speed of 2 m s⁻¹.

Subsequently, two sets of validation data were generated using signals that were entirely different from those used for training. As shown in Figs. 14 and 15, the predicted currents generated by the ICSA-ANFIS inverse model remain highly consistent with the actual currents in the validation datasets. Moreover, the damping force calculated based on the predicted current can accurately track the actual damping force. These results indicate that the inverse model is capable of accurately characterizing the inverse dynamic behavior of the MR damper.

5.2 Simulation analysis of semi-active seat suspension control

To evaluate the effectiveness of the MR damper in the seat suspension system, the proposed ICSA-ANFIS combined with the ADRC control strategy was compared with the ANFIS-ADRC and CSA-ANFIS-ADRC control schemes. The simulation conditions assumed that the agricultural vehicle operated on an ISO class-D road at speeds of 1 and 2 m s⁻¹ under both random road profiles and shock-type road disturbances. Table 3 presents the parameters of the seat suspension system, including the spring and damper characteristics, along with the corresponding simulation results.

Figures 16 to 19 present the simulation results of vertical vibration acceleration and dynamic suspension travel for the agricultural vehicle's seat suspension system. These figures indicate that the ICSA-ANFIS-ADRC control strategy is more effective in improving the vibration isolation performance of the seat suspension system compared to CSA-ANFIS-ADRC and standard ANFIS-ADRC controllers. To gain a more comprehensive understanding of the performance advantages of the ICSA-ANFIS-ADRC in seat suspension applications, Table 4 summarizes the rms values of relevant suspension performance indicators. As shown in Table 4, the rms values of vertical vibration acceleration are reduced to varying degrees by the ICSA-ANFIS-ADRC compared with the other two suspension strategies. Regarding the dynamic travel of the seat suspension, the ICSA-ANFIS-ADRC system exhibits a significantly reduced dynamic stroke at both 1 and 2 m s⁻¹ in comparison to the ANFIS-ADRC and CSA-ANFIS-ADRC systems. Moreover, the control performance remains robust under both shock and non-shock road conditions. The increase in external shock disturbances did not result in a significant increase in seat dynamic travel, further reflecting the robustness of the ICSA-ANFIS-ADRC control strategy.

Figure 18 Seat vertical acceleration response under ISO class-D road excitation with shock disturbance at 1 m s⁻¹.

Figure 19 Seat vertical acceleration response under ISO class-D road excitation with shock disturbance at 2 m s⁻¹.

According to ISO 2631-1 guidelines for evaluating human exposure to whole-body vibration, the level of ride comfort can be directly correlated with the rms value of vertical acceleration (International Organization for Standardization, 2025). Specifically, rms values above 0.8 m s⁻² are typically categorized as “uncomfortable,” whereas values in the range of 0.5–0.8 m s⁻² correspond to “fairly uncomfortable to fairly comfortable”, and values below 0.5 m s⁻² indicate “comfortable” vibration levels. In this study, the proposed ICSA-ANFIS-ADRC controller achieved rms reductions of up to 32.9 % compared with baseline controllers, effectively lowering the seat acceleration from the “uncomfortable” range into the “fairly comfortable” zone under ISO class-D road excitations. This improvement demonstrates not only a numerical enhancement in control performance but also a tangible reduction in driver discomfort, thereby confirming the practical significance of the proposed control strategy in real-world agricultural vehicle applications.

As quantitatively summarized in Table 4, the ICSA-ANFIS-ADRC suspension system achieved a significant reduction in the rms values of vertical vibration acceleration compared to the ANFIS-ADRC and CSA-ANFIS-ADRC suspension systems. Under non-shock conditions at 1 m s⁻¹, the reduction reached 13.35 % compared to CSA-ANFIS-ADRC and 25.5 % compared to ANFIS-ADRC. At a speed of 2 m s⁻¹, this performance improvement was further amplified, with reductions of 18.5 % and 32.9 %, respectively, demonstrating the system's adaptability to speed variation. Notably, under shock conditions, the control strategy maintained its effectiveness, achieving reductions of 12.0 % and 28.4 % at 1 m s⁻¹ and 13.5 % and 32.2 % at 2 m s⁻¹, relative to the baseline systems. Critically, the controller's performance did not deteriorate despite increased disturbance intensity, thereby confirming the robustness of the ICSA-ANFIS-ADRC architecture under varying operational conditions.

The numerical and experimental results exhibit strong agreement, yet several important insights can be drawn from the observed differences and performance trends. First, the superior convergence behavior of the ICSA-ANFIS model compared with standard ANFIS and CSA-ANFIS can be attributed to the triangular probability distribution and mutation mechanism introduced in ICSA. These strategies effectively prevent premature convergence to local optima, ensuring global exploration of the solution space and yielding higher current prediction accuracy. As a result, the inverse model produces a damping force that closely matches the experimental measurements.

Second, the discrepancies between numerical simulations and experimental data primarily arise from unmodeled nonlinearities in the MR damper, such as temperature effects on the MR fluid viscosity and minor frictional losses in the mechanical joints of the test setup. Nevertheless, these discrepancies remain within acceptable margins, confirming that the improved Bouc–Wen model captures the dominant hysteresis and rate-dependent characteristics of the damper.

Third, the enhanced vibration isolation performance of the ICSA-ANFIS-ADRC control strategy under both random and shock road conditions can be explained by its hybrid structure. The ADRC component effectively estimates and compensates for unknown disturbances and model uncertainties in real time, while the ICSA-optimized ANFIS inverse model provides accurate current commands to the MR damper, thereby mitigating its inherent nonlinear hysteresis. This synergy enables the controller to maintain robust damping performance even under sudden excitations or speed variations, as demonstrated by the reduced rms acceleration values in Table 4.

Finally, the robustness of the proposed controller is further highlighted by the fact that the vibration reduction ratios remain nearly constant despite increasing external disturbances. This indicates that the ICSA-ANFIS-ADRC framework achieves not only high precision in force tracking but also strong adaptability across different operating conditions, thereby providing a reliable solution for intelligent seat suspension systems in agricultural vehicles.

6 Conclusions

In this study, an innovative semi-active control framework, referred to as ICSA-ANFIS-ADRC, was proposed by integrating the ICSA with an ANFIS and coupling it with the ADRC strategy, specifically designed for real-time dynamic optimization of MR seat suspension dampers. The proposed approach first involves the development of a forward dynamic damper model based on an improved Bouc–Wen hysteresis model, which was empirically calibrated using experimental datasets. Subsequently, an inverse model of the damper was developed based on the ICSA-ANFIS architecture, with detailed descriptions of the neuro-fuzzy network configuration and the dynamic optimization procedure driven by ICSA. By coupling the ICSA-ANFIS inverse model with an ADRC controller tailored for seat suspension systems, a novel hybrid control strategy was proposed and implemented. Extensive numerical simulations validated the effectiveness of the proposed controller, revealing three major advancements: superior precision in damping force regulation, continuously adjustable damping across the entire operational frequency range, and significantly enhanced vibration isolation performance of the seat suspension system under dynamic loading conditions.

Although the proposed ICSA-ANFIS-ADRC control framework demonstrates promising results in vibration isolation and ride comfort enhancement, several limitations should be acknowledged. First, the study primarily relies on simulation analysis; therefore, future work will focus on conducting hardware-in-the-loop tests and full-scale experimental validation on agricultural vehicles to confirm real-world applicability. Second, the present model did not explicitly account for temperature effects, long-term damper wear, or parameter drift of MR fluids, which may influence control accuracy over time. Future research will address these factors by incorporating adaptive parameter identification and robust controller reconfiguration. Finally, integration of human biodynamic models and ISO 2631 comfort assessment will be pursued to provide a more comprehensive evaluation of driver health and long-term exposure risks.

Appendix A: Nomenclature

Symbol	Description
mb	Mass of human body (kg)
ms	Mass of seat (kg)
mc	Mass of cab (kg)
kb	Stiffness of human body (N m−1)
cb	Damping coefficient of human body (N s m−1)
ks	Stiffness of seat suspension (N m−1)
cMR	Equivalent damping coefficient of MR
	damper (N s m−1)
kc	Stiffness of cab (N m−1)
cc	Damping coefficient of cab (N s m−1)
zb	Vertical displacement of human body (m)
zs	Vertical displacement of seat (m)
zc	Vertical displacement of cab (m)
z0	Base excitation displacement (m)
FMR	Damping force of MR damper (N)
α	Scaling parameter of hysteresis force
$\mathit{\beta },\mathit{\gamma },A$	Parameters controlling smoothness and
	transition of Bouc–Wen model
n	Rounding coefficient of hysteresis model
c0,c1	Viscous damping coefficients of MR damper
k0,k1	Stiffness coefficients of MR damper
x	Total displacement of MR damper piston (m)
y	Hysteresis displacement (m)
x0	Equilibrium position offset (m)
I	Control current of MR damper (A)