Proportional + integral + derivative control of nonlinear full-car electrohydraulic suspensions using global and evolutionary optimization techniques

Abstract

Resolving the trade-offs between suspension travel, ride comfort, road holding, vehicle handling and power consumption is the primary challenge in the design of active vehicle suspension system. Multi-loop proportional + integral + derivative controllers’ gains tuning with global and evolutionary optimization techniques is proposed to realize the best compromise between these conflicting criteria for a nonlinear full-car electrohydraulic active vehicle suspension system. Global and evolutionary optimization methods adopted include: controlled random search, differential evolution, particle swarm optimization, modified particle swarm optimization and modified controlled random search. The most improved performance was achieved with the differential evolution algorithm. The modified particle swarm optimization and modified controlled random search algorithms performed better than their predecessors, with modified controlled random search performing better than modified particle swarm optimization in all aspects of performance investigated both in time and frequency domain analyses.

Introduction

Good active vehicle suspension system (AVSS) is a control challenge in the design of vehicle suspension systems. The challenge includes the determination of the optimal trade-offs between conflicting suspension performance parameters like suspension travel, ride comfort, road holding and vehicle handling. For example, compromise must be reached when a hard suspension with limited suspension travel is required for good road holding, and a soft suspension is desired for a smooth and comfortable ride.

Another challenge faced in designing control system for AVSS is its nonlinear characteristics, especially due to complex dynamics of the electrohydraulic actuators.^1,2 Meanwhile the system is highly prone to operational and system uncertainties. The full-car is characterized by nonlinear coupling and multiple-variable control situation.^3–5

Whilst the superior performance of active suspension is not in doubt, it is unable to claim similar level of commercial acceptability like passive vehicle suspension systems (PVSSs), because it is more complex operationally and structurally.^6–8 Recent technological growth is creating opportunities for implementation of complex control methodologies in AVSS design. Linear optimal control schemes are well developed and already employed in the controller designs for AVSS and semi-active vehicle suspension systems.^8–17 Stability and robust characteristics are well established for these techniques, but these properties become highly limited when implemented in nonlinear AVSS. ¹⁸

Successful application usually depends on availability of good dynamic model of the system and availability of all the states. Implementation of nonlinear techniques like feedback linearization requires guarantee of stable zero dynamics in the system; backstepping requires repeated differentiation of the system’s nonlinear function and their implementation practically is usually characterized by chattering.^19–23 Combining nonlinear control schemes with computational intelligence techniques have largely been demonstrated to be effective, ²⁴ but it comes with the additional computational complexity that is associated with each scheme in the process. Demonstrating system stability can also be very challenging.^25,26

The case for application of proportional + integral + derivative (PID) control design to AVSS has been made and demonstrated. ⁵ Control design for full-car AVSS in this work requires a cascaded feedback loop arrangement with inner loops stabilizing the actuator dynamics with several PID feedback control loops to tune simultaneously in a system characterized by couplings amongst its parameters; employment of global and evolutionary algorithms for autotuning cannot be avoided. Its potentials have been documented already in the literature.^27–36 Benefit of this approach is that the optimization algorithms can use their objective functions to satisfy specified or required performance criterion.

Controlled random search (CRS) and differential evolution (DE) have been successfully used to compute optimal parameters for various nonlinear systems and engineering applications.^37–40 Classical particle swarm optimization (PSO)-tuned PID control has been employed in many control designs. Though it was initially designed for continuous space, it performs well when applied to discontinuous objective functions. However, it is prone to premature convergence in large-scale complex problems like most of the other stochastic algorithms.^35,36

The solution space in this research is expected to contain many peaks and troughs as the system is highly nonlinear with seven degree-of-freedom (DOF), thereby making the optimization process more challenging. Good AVSS designs require adaptive control properties in spite of the complex computation demands that arise through the use of sophisticated control methods. PID controllers tuned using global and evolutionary algorithms are able to take up this advantage, thereby controlling with some level of intelligence. ³⁵

The setbacks found in the performances of the algorithms are used in modifying them for better performances. These modifications are, however, uniquely suitable for particular problems. ³⁶ In this work, modified variants of CRS and PSO routines are introduced to address their shortcomings. For example, for PSO, for the ith particle x_i, instead of exploring the area around the vicinity of each particle best solution P_i, exploration will be conducted around any randomly selected particle R_i that has a better fitness value than the particle of interest, i.e. $f\left(R_i\right)<f\left(x_i\right)$ . This modification prevents particles from drifting to dead zones within the search space.

In the modified CRS, the underlining centre of gravity will be based on only a few randomly selected individuals instead of the whole population. For example, instead of creating a simplex using n + 1 solutions (where n solutions are chosen randomly), one being the overall best, a simplex using only three solutions is created. The worse of the three is reflected through the centre of gravity of the remaining two. This will add flexibility to the algorithm and improve the exploration capabilities of the algorithm. Hence, a comparative analysis will be performed to investigate the effects of these modifications.^41,42

The major contributions of the paper are as follows:

Selection of the cascaded PID controllers’ parameters for the seven DOF nonlinear electrohydraulic AVSS using global and evolutionary optimization techniques to address the trade-offs amongst the conflicting performance requirements; and

Frequency-domain analysis of the proposed DE-optimized, PID-based multi-loop controller within the whole-body vibration range of $0-80\text{H}z$ .

Employment of five different global and evolutionary optimization shows the level of complexity and computational difficulty used in achieving the design of the PID-based controller in this work with 24 optimized parameters.

The rest of the paper is organized as follows. The next section presents a brief description of the physical, mathematical and road disturbance input models for the AVSS. ‘Controller implementation’ section presents the controller design, the system performance specifications and evaluation criteria. ‘Evolutionary and global optimization algorithms’ section describes in detail the global and evolutionary optimization algorithms. Simulation results and their discussion are presented in ‘Simulation results and discussion’ section. Concluding remarks are given in the final section.

System overview and modelling

Mathematical modelling

A schematic of the full-car model is presented in Figure 1 and the governing equations are derived using Newton–Euler methods.

Figure 1.

Simplified full-car model for AVSS.

Figure 2.

Road disturbance input for the full-car AVSS.

Figure 3.

Control architecture of the proposed multi-loop PID controller. PID: proportional + integral + derivative.

Figure 4.

Convergence history plots for PSO, MPSO and DE. DE: differential evolution; MPSO: modified particle swarm optimization; PSO: particle swarm optimization.

Vertical forces acting on the chassis are given as

$$F_{\mathit{fr}}=F_{\mathit{kfr}}+F_{\mathit{bfr}}-F_{\mathit{afr}}$$

(1)

$$F_{\mathit{fl}}=F_{\mathit{kfl}}+F_{\mathit{bfl}}-F_{\mathit{afl}}$$

(2)

$$F_{\mathit{rr}}=F_{\mathit{krr}}+F_{\mathit{brr}}-F_{\mathit{arr}}\text{ and }$$

(3)

$$F_{\mathit{rl}}=F_{\mathit{krl}}+F_{\mathit{brl}}-F_{\mathit{arl}}$$

(4)

with F_kij, F_bij and F_aij denoting the spring, damping and actuator force components within the suspension. The subscripts i and j represent the position of the particular suspension system with i taking on either front f or rear r which pertains to the longitudinal position of the suspension on the vehicle body; j takes on right r or left position l on the lateral axis of the vehicle body. The vertical, pitch and roll equations of motion are given as ⁴³

$$M_s\ddot{z}=F_{\mathit{fr}}+F_{\mathit{fl}}+F_{\mathit{rr}}+F_{\mathit{rl}}$$

(5)

$$I_{\theta }\ddot{\theta }=-F_{\mathit{fr}}{l}_f-F_{\mathit{fl}}{l}_f+F_{\mathit{rr}}{l}_r+F_{\mathit{rl}}{l}_r$$

(6)

$$I_{\varphi }\ddot{\varphi }=\frac{a_f}{2}\left[F_{\mathit{fr}}-F_{\mathit{fl}}+F_{\mathit{rr}}-F_{\mathit{rl}}\right]$$

(7)

where M_s is the total mass of the vehicle; $I_{\theta }$ and $I_{\varphi }$ are the pitch and roll inertias of the vehicle, respectively; l_r and l_f are the longitudinal distances of the front and rear suspensions from the centre of gravity of the chassis, respectively; and a_f is the lateral distance of each suspension from the centre of gravity. The displacement induced at each corner of the chassis is as follows

$$z_{\mathit{lf}}=z-l_f\mathrm{\hspace{0.17em}sin}\theta -\frac{a_f}{2}\mathrm{\hspace{0.17em}sin}\varphi$$

(8)

$$z_{\mathit{rf}}=z-l_f\mathrm{\hspace{0.17em}sin}\theta +\frac{a_f}{2}\mathrm{\hspace{0.17em}sin}\varphi$$

(9)

$$z_{\mathit{lr}}=z+l_r\mathrm{\hspace{0.17em}sin}\theta -\frac{a_f}{2}\mathrm{\hspace{0.17em}sin}\varphi \text{ and }$$

(10)

$$z_{\mathit{rr}}=z+l_r\mathrm{\hspace{0.17em}sin}\theta +\frac{a_f}{2}\mathrm{\hspace{0.17em}sin}\varphi$$

(11)

The governing equation at each wheel is given as

$$m_{\mathit{uij}}{\ddot{z}}_{\mathit{tij}}=-F_{\mathit{kij}}-F_{\mathit{bij}}+F_{\mathit{aij}}+F_{\mathit{ktij}}+F_{\mathit{btij}}$$

(12)

where m_uij, z_tij, F_ktij and F_btij are the mass, vertical displacement, spring force and damping force of the $(i,j)\mathrm{th}$ wheel, respectively. The suspension subcomponents contain linear and nonlinear elements which are described as follows

$$F_{\mathit{kij}}=k_{\mathit{sij}}^l\left(z_{\mathit{tij}}-z_{\mathit{ij}}\right)+k_{\mathit{sij}}^{\mathit{nl}}{\left(z_{\mathit{tij}}-z_{\mathit{ij}}\right)}^3$$

(13)

$$F_{\mathit{bij}}=b_{\mathit{sij}}^l\left({\dot{z}}_{\mathit{tij}}-{\dot{z}}_{\mathit{ij}}\right)+b_{\mathit{sij}}^{\mathit{nl}}\sqrt{\left|{\dot{z}}_{\mathit{tij}}-{\dot{z}}_{\mathit{ij}}\right|}s\mathit{gn}\left({\dot{z}}_{\mathit{tij}}-{\dot{z}}_{\mathit{ij}}\right)-b_{\mathit{sij}}^{\mathit{sym}}\left|{\dot{z}}_{\mathit{tij}}-{\dot{z}}_{\mathit{ij}}\right|$$

(14)

with linear and nonlinear spring constants $k_{\mathit{sij}}^l$ and $k_{\mathit{sij}}^{\mathit{nl}}$ ; linear, nonlinear and symmetric damping constants $b_{\mathit{sij}}^l,b_{\mathit{sij}}^{\mathit{nl}}$ and $b_{\mathit{sij}}^{\mathit{sym}}$ . The wheels’ characteristics are assumed to be linear

$$F_{\mathit{ktij}}=k_{\mathit{tij}}\left(w_{\mathit{ij}}-z_{\mathit{tij}}\right)$$

(15)

$$F_{\mathit{btij}}=b_{\mathit{tij}}\left({\dot{w}}_{\mathit{ij}}-{\dot{z}}_{\mathit{tij}}\right)$$

(16)

where k_tij and b_tij are, respectively, the linear tyre stiffness and tyre damping constants.

The mathematical relations of the hydraulic actuator force are given as

$${\dot{F}}_{\mathit{aij}}=A_{\mathit{ij}}{\dot{P}}_{\mathit{Lij}}=A_{\mathit{ij}}\left[{\gamma }_{\mathit{ij}}{\Phi }_{\mathit{ij}}{x}_{\mathit{vij}}-{\beta }_{\mathit{ij}}{P}_{\mathit{Lij}}+{\Theta }_{\mathit{ij}}{A}_{\mathit{ij}}\left({\dot{z}}_{\mathit{ij}}-{\dot{z}}_{\mathit{tij}}\right)\right]$$

(17)

where A_ij are the hydraulic cylinders’ areas, and P_Lij are the pressure differences across the pistons. The pressures generated by the motions of the spool valves x_vij are ${\gamma }_{\mathit{ij}}{\Phi }_{\mathit{ij}}{x}_{\mathit{vij}},{\beta }_{\mathit{ij}}{P}_{\mathit{Lij}}$ are the pressure losses and ${\Theta }_{\mathit{ij}}{A}_{\mathit{ij}}\left({\dot{z}}_{\mathit{ij}}-{\dot{z}}_{\mathit{tij}}\right)$ are the pressure variations due to the damping nature of the actuators. Hydraulic parameters in these relations are as follows

$${\Theta }_{\mathit{ij}}=\frac{4{\beta }_{e_{\mathit{ij}}}}{V_{t_{\mathit{ij}}}},{\beta }_{\mathit{ij}}={\Theta }_{\mathit{ij}}{C}_{t{p}_{\mathit{ij}}},{\gamma }_{\mathit{ij}}=C_{d_{\mathit{ij}}}{S}_{\mathit{ij}}\frac{1}{{\rho }_{\mathit{ij}}}$$

where $V_{t_{\mathit{ij}}}$ are the volumes of the cylinders, ${\beta }_{e_{\mathit{ij}}}$ are the bulk moduli of the working fluids, $C_{t{p}_{\mathit{ij}}}$ are the leakage coefficients, $C_{d_{\mathit{ij}}}$ are the discharge coefficients between the supply lines and the hydraulic cylinders, S_ij are the spool valves’ area gradients and ρ_ij are the densities of fluids in each cylinder. The valves are manipulated through motors, which are modelled as first-order elements

$${\dot{x}}_{\mathit{vij}}=\frac{1}{{\tau }_{\mathit{ij}}}\left(K_{\mathit{vij}}{u}_{\mathit{ij}}-x_{\mathit{vij}}\right)$$

(18)

where u_ij are the control input voltages, K_vij are the valves’ gains and τ_ij are the time constants of the $(i,j)\mathrm{th}$ actuators.

The heave, pitch and roll motions of the vehicle are induced by travelling over a bump which varies in height along its lateral length as shown in Figure 2. The road profile at each wheel is given as

$$w_{\mathit{fr}}\left(t\right)=\{\begin{array}{cc}{a}_1\frac{(1-\cos 2\pi \frac{V}{\lambda }t)}{2},& 0.45\le t\le 0.9,\\ 0,& \text{otherwise}\end{array}$$

(19)

$$w_{\mathit{fl}}\left(t\right)=\{\begin{array}{cc}{a}_2\frac{(1-\cos 2\pi \frac{V}{\lambda }t)}{2},& 0.45\le t\le 0.9,\\ 0,& \text{otherwise}\end{array}$$

(20)

$$w_{\mathit{rr}}\left(t\right)=\{\begin{array}{cc}{a}_1\frac{(1-\cos 2\pi \frac{V}{\lambda }t)}{2},& 0.45+\frac{\lambda }{V}\le t\le 0.9+\frac{\lambda }{V},\\ 0,& \text{otherwise}\end{array}$$

(21)

$$w_{\mathit{rl}}\left(t\right)=\{\begin{array}{cc}{a}_2\frac{\left(1-\cos 2\pi \frac{V}{\lambda }t\right)}{2},& 0.45+\frac{\lambda }{V}\le t\le 0.9+\frac{\lambda }{V},\\ 0,& \text{otherwise}\end{array}$$

(22)

where a₁, a₂ are the amplitudes of the bump, V is the horizontal speed and λ is the wavelength of the bump. Values of the system parameters are given in Table 1.

Table 1.

Parameters of the full-car model.

Parameters	Value
Sprung mass (Ms)	1060 kg
Unsprung masses ( $m_{u_{\mathit{ij}}}$ )	40, 40, 35, 35 kg
Pitch moment of inertia ( $I_{\theta }$ )	2200 kgm2
Roll moment of inertia ( $I_{\varphi }$ )	460 kgm2
Distance from vehicle front axle	1 m
to its centre of gravity ( $l_f$ )
Distance from vehicle rear axle	1.5 m
to its centre of gravity ( $l_r$ )
Vehicle width ( $a_f$ )	1.5 m
Linear suspension stiffness	2.35 × 104 N/m
at each wheel ( $k^l{\hspace{0.17em}}_{s_{\mathit{ij}}}$ )
Nonlinear suspension stiffness	2.35 × 106 N/m
at each wheel ( $k^{\mathit{nl}}{\hspace{0.17em}}_{s_{\mathit{ij}}}$ )
Tyre stiffness at each wheel ( $k_{t_{\mathit{ij}}}$ )	1.9 × 105 N/m
Linear suspension damping	700 N s/m
each wheel ( $b^l{\hspace{0.17em}}_{s_{\mathit{ij}}}$ )	700 N s/m
Nonlinear suspension damping	400 N s/m
each wheel ( $b^{\mathit{nl}}{\hspace{0.17em}}_{s_{\mathit{ij}}}$ )
Asymmetric suspension damping	400 N s/m
at each wheel ( $b^{\mathit{sym}}{\hspace{0.17em}}_{s_{\mathit{ij}}}$ )
Tyre damping at each wheel ( $b_{t_{\mathit{ij}}}$ )	80, 80, 70, 70 N s/m
Actuator parameter (αij)	4.515 × 1013
Actuator parameter (βij)	1
Actuator parameter (γij)	1.545 × 109
Piston area (Aij)	3.35 × 10−4 m2
Supply pressure ( $P_{s_{\mathit{ij}}}$ )	10,342,500 Pa
Time constant (τij)	3.33 × 10−2 s
Servo-valve gains ( $K_{v_{\mathit{ij}}}$ )	0.001 m/V
Bump amplitudes (a1, a2)	4, 6 cm
Bump wavelength (λ)	5 m
Vehicle speed (V)	40 km/h

Controller implementation

Figure 3 shows the architecture of the multi-loop PID control investigated in the present paper. It includes four outer loops to regulate the controlled variable and four inner force feedback loops are placed at each actuator to ensure stability of the hydraulic actuator. The control variable is the suspension travel as this output is the primary factor that captures the system dynamics. The outer control loops are modelled as follows

$$F_{d_{\mathit{ij}}}=K_{P_{\mathit{ij}}}{e}_{1_{\mathit{ij}}}+K_{D_{\mathit{ij}}}\frac{d{e}_{1_{\mathit{ij}}}}{\mathit{dt}}+K_{I_{\mathit{ij}}}{\displaystyle {\int }_0^T{e}_{1_{\mathit{ij}}}\text{d}t}$$

(23)

where $F_{d_{\mathit{ij}}},K_{P_{\mathit{ij}}},K_{D_{\mathit{ij}}},K_{I_{\mathit{ij}}},e_{1_{\mathit{ij}}}$ are the desired control forces computed by the outer-loop PID controllers; proportional, integral and derivative gains of the respective outer-loop PID controllers, and the control error signals, respectively. The control error signals $e_{1_{\mathit{ij}}}$ are given as

$$e_{1_{\mathit{ij}}}=R_{d_{\mathit{ij}}}\left(t\right)-y_{\mathit{ij}}$$

(24)

where $R_{d_{\mathit{ij}}}(t)$ is the suspension travel setpoint. This investigation addresses disturbance rejection problem, and hence the desired setpoints are set to zero.

In a similar manner, the inner control loop determines the control voltage sent to the hydraulic actuator. The governing equations of these inner-loop controllers are as follows

$$e_{2_{\mathit{ij}}}=F_{a_{\mathit{ij}}}-F_{d_{\mathit{ij}}}$$

(25)

$$u_{\mathit{ij}}=k_{p_{\mathit{ij}}}{e}_{2_{\mathit{ij}}}+k_{d_{\mathit{ij}}}\frac{d{e}_{2_{\mathit{ij}}}}{\mathit{dt}}+k_{i_{\mathit{ij}}}{\displaystyle {\int }_0^T{e}_{2_{\mathit{ij}}}\text{d}t}$$

(26)

where $e_{2_{\mathit{ij}}}$ are the deviations between the actual actuator forces generated and their desired setpoints; $k_{p_{\mathit{ij}}}$ , $k_{i_{\mathit{ij}}}$ and $k_{d_{\mathit{ij}}}$ are the proportional, integral and derivative gains of the respective inner-loop PID controllers; and u_ij are the control voltages supplied to the hydraulic actuators. The gains of each of the eight PID controllers are chosen with the following specifications in mind:

• 1. The controller should demonstrate a good low-frequency disturbance rejection;

• 2. Satisfactory transient response with minimal oscillations after the disturbance has disappeared; that is:

•  • the rise time should not be greater than 0.1 s,

•  • the maximum overshoot should be less than 5%

•  • and zero steady-state error.

• 3. Suspension travels y_ij are constrained to a maximum of ±0.1 m.

• 4. Control input voltages u_ij are constrained to ±10 V due to the limitations of the power supply.

• 5. The total actuation force must be less than the vehicle weight to ensure that the vehicle does not leave the ground, i.e. $\left|F_{a_{\mathit{ij}}}\right|<M_{s}g$ .

• 6. The body-heave acceleration ( $\ddot{z}$ ) should be less than 4.5 m/s² and the corresponding root-mean-square (RMS) value should be as low as possible, typically in the range $0-1$ m/s² for the ride comfort to lie in the ‘less discomfort’ region of the International Organization for Standardization (ISO) 2631. ⁴⁴

• 7. The performance index J which addresses each of the conflicting design criteria of an AVSS is to be minimized. This index has the following form

  $$J=J_1+J_2+J_3+J_4+J_5$$

  (27)
  $$J_1=\frac{1}{T}{\displaystyle {\int }_0^T\left[{\left(\frac{\ddot{z}}{{\ddot{z}}_{\mathit{max}}}\right)}^2+{\left(\frac{\ddot{\theta }}{{\ddot{\theta }}_{\mathit{max}}}\right)}^2+{\left(\frac{\ddot{\varphi }}{{\ddot{\varphi }}_{\mathit{max}}}\right)}^2\right]}\text{d}t$$

  (28)
  $$J_2={\displaystyle \sum _{j=\mathit{fr}}{\displaystyle \sum _{i=\mathit{rl}}\frac{1}{T}}}{\displaystyle {\int }_0^T\left[{\left(\frac{F_{k{t}_{\mathit{ij}}}+F_{b{t}_{\mathit{ij}}}}{{\left(F_{k{t}_{\mathit{ij}}}+F_{b{t}_{\mathit{ij}}}\right)}_{\mathit{max}}}\right)}^2\right]}\text{d}t$$

  (29)
  $$J_3={\displaystyle \sum _{i=1}^4\frac{1}{T}}{\displaystyle {\int }_0^T\left[{\left(\frac{y_i}{y_{i_{\mathit{max}}}}\right)}^2\right]}\text{d}t$$

  (30)
  $$J_4={\displaystyle \sum _{j=\mathit{fr}}{\displaystyle \sum _{i=\mathit{rl}}\frac{1}{T}}}{\displaystyle {\int }_0^T{\left(\frac{F_{a_{\mathit{ij}}}}{F_{a_{i{j}_{\mathit{max}}}}}\right)}^2}\text{d}t$$

  (31)
  $$J_5={\displaystyle \sum _{i=1}^4\frac{1}{T}}{\displaystyle {\int }_0^T\left[{\left(\frac{u_i}{u_{i_{\mathit{max}}}}\right)}^2\right]}\text{d}t$$

  (32)

where J is the performance index and J₁, J₂, J₃, J₄ and J₅ relate to the vehicle ride comfort and vehicle handling, road holding properties, suspension travel, actuation force and power consumption, respectively. ${\ddot{z}}_{\mathit{max}},{\ddot{\theta }}_{\mathit{max}},{\ddot{\varphi }}_{\mathit{max}},{\left(F_{k{t}_{\mathit{ij}}}+F_{b{t}_{\mathit{ij}}}\right)}_{\mathit{max}},y_{i_{\mathit{max}}},F_{a_{i{j}_{\mathit{max}}}},u_{i_{\mathit{max}}}$ are the maximum allowable body-heave, pitch and roll accelerations, tyre dynamic load, suspension travels, actuation forces and control input voltages, respectively. T is the period over which the simulation runs. The controller gains are determined manually and through the use of the selected global/evolutionary optimization algorithms.

Evolutionary and global optimization algorithms

Five optimal policies (CRS, DE, PSO, modified controlled random search (MCRS) and modified particle swarm optimization (MPSO)) will be used to select the controllers’ gains for the full-car AVSS. In this section the resulting AVSS performance and the performance of the optimal routines relative to one another will be studied. The objective function has the form presented in equations (27) to (32). The rest of this section will be devoted to describing the algorithms.

DE optimization algorithm

Evolutionary algorithms which include DE, GA and PSO are random search optimization methods where the optimal solution is produced through the evolution of a random population set $S=\left\{{\mathbf{x}}_0,{\mathbf{x}}_1,{\mathbf{x}}_2,\mathbf{\dots },{\mathbf{x}}_{\mathbf{N}}\right\}$ with each individual denoted as ${\mathbf{x}}_i$ . These individuals are vectors with each of its elements pertaining to a specific controller gain. Each of the proposed algorithms only differs in the manner in which the population changes through each generation as well as the conditions that must be satisfied in order for the respective individuals to change. The search space may be predefined to operate within a feasible region in order to improve computation time and convergence characteristics. Such a search space is chosen through intuitive reasoning and experience gained through manual tuning as in the case of PID controller tuning. ⁴⁵

In each generation step of DE, an associated trial individual ${\mathbf{y}}_i$ is generated for each targeted individual ${\mathbf{x}}_i$ in S. This trial individual ${\mathbf{y}}_i$ is a function of ${\mathbf{x}}_i$ and a mutated individual ${\widehat{\mathbf{x}}}_i$ . ${\widehat{\mathbf{x}}}_i$ is generated from three distinct randomly chosen individuals ${\mathbf{x}}_{\alpha }$ , ${\mathbf{x}}_{\beta }$ and ${\mathbf{x}}_{\gamma }$ from the population set S in the order of ascending fitness value, respectively. Each problem variable k (i.e., each PID gain) of the mutated individual ${\widehat{\mathbf{x}}}_i$ is determined using the following relation

$${\widehat{\mathbf{x}}}_{i_k}={\mathbf{x}}_{{\alpha }_k}+F\left({\mathbf{x}}_{{\beta }_k}-{\mathbf{x}}_{{\gamma }_k}\right)$$

(33)

where, α, β and γ are representative of the three distinctive randomly selected individuals from the population S, k signifies the kth parameter of these individuals, i is the ith individual that is currently being mutated and F is the primary DE scaling parameter.

Crossover is thereafter performed between the targeted individual ${\mathbf{x}}_i$ and its mutated counterpart ${\widehat{\mathbf{x}}}_i$ to produce a corresponding new trial point ${\mathbf{y}}_i$ . In crossover, each parameter k within ${\mathbf{y}}_i$ either equals the associating kth parameter in ${\mathbf{x}}_i$ or ${\widehat{\mathbf{x}}}_i$ . This decision is made by producing a random number R^k and if R^k exceeds a predefined threshold number C_R the kth parameter of ${\mathbf{y}}_i$ correlates to the kth parameter of ${\mathbf{x}}_i$ . Conversely, if this condition is not met the kth parameter of ${\mathbf{y}}_i$ is given the same value as the kth parameter of ${\widehat{\mathbf{x}}}_i$ . These relationships are stated as follows

$${\mathbf{y}}_{i_k}=\{\begin{array}{l}{\widehat{\mathbf{x}}}_{i_k}\text{ if } R^k\le C_R\hfill \\ {\mathbf{x}}_{i_k}\text{otherwise}\hfill \end{array}$$

(34)

The trial vector ${\mathbf{y}}_i$ of the targeted individual ${\mathbf{x}}_i$ from S is only accepted and replaces its predecessor or previous values ${\mathbf{x}}_i\left(t-1\right)$ if and only if it possesses a superior fitness to that of ${\mathbf{x}}_i\left(t-1\right)$ which may be summarized as

$${\mathbf{x}}_i=\{\begin{array}{l}{\mathbf{y}}_i\text{ if } f\left({\mathbf{y}}_i\right)<f\left({\mathbf{x}}_i\right)\hfill \\ {\mathbf{x}}_i\text{otherwise}\hfill \end{array}$$

(35)

where $f({\mathbf{y}}_i)$ is the fitness value of the trial vector, and $f({\mathbf{x}}_i)$ is the fitness of the targeted individual. After the whole population S has evolved the process is repeated until the stopping criterion is met. Thereafter, the optimal solution is chosen to be the individual in S with the best fitness value. In the case of this investigation, the algorithm was set to stop after a predefined number of iterations. This routine is summarized in the following steps:

• Step 1 Create a random population $S({\mathbf{x}}_0,{\mathbf{x}}_1,\dots ,{\mathbf{x}}_{\mathbf{N}})$ and let ${\mathbf{x}}_0$ be the initial condition which in the case of controller tuning matches the set of controller gains attained through manual tuning.

• Step 2 If stopping criterion has been met, select the fittest individual in S as the optimal solution; otherwise continue to the next step.

• Step 3 In each iteration, evolve the ith individual ${\mathbf{x}}_i$ in population S as follows: Randomly choose three mutually independent individuals from S and denote them as ${\mathbf{x}}_{\alpha },{\mathbf{x}}_{\beta }$ and ${\mathbf{x}}_{\gamma }$ .

• Compute the trial candidate ${\mathbf{y}}_i$ by applying crossover with the individual ${\mathbf{x}}_i$ and the mutated individual ${\widehat{\mathbf{x}}}_i$ according to equation (34):

• Replace the ith individual ${\mathbf{x}}_i$ in S with the candidate solution ${\mathbf{y}}_i$ if and only if ${\mathbf{y}}_i$ contains a better fitness value than its predecessor.

• Repeat the process for each individual within S and return to step 2.

We have used $C_R=0.5$ , F = 2, k_max = 150 and N = 100.

PSO algorithm

In the case of PSO, the search space is reflected as an n-dimensional world (n represents the number of parameters being optimized) where swarms of animals or particles of a random population set $S=\{{\mathbf{x}}_0,{\mathbf{x}}_1,{\mathbf{x}}_2,\dots ,{\mathbf{x}}_{\mathbf{N}}\}$ flock to search for food or in this case the optimal solution. After the initial population or swarm has been defined, the fittest individual or particle is registered and every particle in the swarm is programmed to travel to a new position described as

$$\mathbf{x}\left(t+1\right)=\mathbf{x}\left(t\right)+\mathbf{V}\left(t+1\right)$$

(36)

where x is the population position matrix, with each row containing the set of problem variables. V is the matrix of particle velocity and holds the velocities at which each particle travels in each direction in the n-dimensional world. In essence, this refers to the amount at which the particles increment by in each direction. t represents the information relating to the previous iteration and t + 1 denotes the information concerning the next iteration.

The velocity at which each particle travels in each direction depends on the particle’s previous velocity in that specified direction, the corresponding position of the fittest particle in that direction and on the corresponding position of the particle’s personal best position. To record the personal best results of each particle, the matrix ${\mathbf{P}}_{\mathit{best}}$ which contains the personal best parameter values for each particle is created and updated for all iterations.

In each iteration, every particle converges to some extent towards both the fittest individual in the population ${\mathbf{G}}_{\mathit{best}}$ , and in the direction of its personal best solution ${\mathbf{P}}_{\mathit{best}}$ . The convergence towards the global best particle is known as the global search and that towards its corresponding personal best particle is known as local search. The rate of convergence in these neighbourhoods is primarily factors of how much weighting is placed on local and global search. The PSO weighting parameters may be adjusted to attain optimal convergence. The velocity matrix is constructed on local and global search vectors and their associated weighting. It is described as follows

$$\mathbf{V}\left(t+1\right)=w_1\mathbf{V}\left(t\right)+w_{2}r\mathit{and}1\left({\mathbf{P}}_{\mathit{best}}-\mathbf{x}\left(t\right)\right)+w_{3}r\mathit{and}2\left({\mathbf{G}}_{\mathit{best}}-\mathbf{x}\left(t\right)\right)$$

(37)

where ${\mathbf{G}}_{\mathit{best}}$ is the position of global best solution, ${\mathbf{P}}_{\mathit{best}}$ is a matrix containing positions of the personal best solution for each particle in the swarm, $\left({\mathbf{P}}_{\mathit{best}}-\mathbf{x}\left(t\right)\right)$ represents the local search and $\left({\mathbf{G}}_{\mathit{best}}-\mathbf{x}\left(t\right)\right)$ denotes the global search contributions. w₁ is the weighting factor of the previous velocity and determines how much the particle’s previous velocity affects the particle’s new velocity, w₂ is the social weighting and regulates how much the local search affects the particle’s new velocity, w₃ is the global weighting and controls how much the global search affects the particle’s new velocity, and rand1 and rand2 are vectors of random numbers with a size of that of the swarm and these vectors ensure the local and global searches for each particle occurs at different rates.

By changing the position of each particle according to the aforementioned equations, each particle has the potential to improve upon its personal best location and has the ability to become the global best particle. Hence after each iteration, the fitness value of each particle is analysed and if it improves from its personal best solution, its personal best location is replaced with its current location. This is further explained as follows

$${\mathbf{P}}_{\mathit{bes}{t}_i}=\{\begin{array}{l}{\mathbf{x}}_i\left(t+1\right)\text{ if } f\left({\mathbf{x}}_i\left(t+1\right)\right)<f\left({\mathbf{P}}_{\mathit{bes}{t}_i}\right)\hfill \\ {\mathbf{P}}_{\mathit{bes}{t}_i}\text{otherwise}\hfill \end{array}$$

(38)

where ${\mathbf{P}}_{\mathit{bes}{t}_i}$ is the personal best location of the ith particle in the swarm, and ${\mathbf{x}}_i\left(t+1\right)$ is the corresponding position of the ith particle that has been computed in the most recent iteration. Moreover, the global best particle may be defined as the fittest particle from personal best matrix as follows

$${\mathbf{G}}_{\mathit{bes}{t}_i}\left(t+1\right)=\arg \underset{{\mathbf{P}}_{\mathit{bes}{t}_i}}{\min }f\left({\mathbf{P}}_{\mathit{bes}{t}_i}\left(t+1\right)\right)$$

(39)

where $f(\dots )$ denotes the respective performance index J of the various particles. The process is repeated in each iteration and terminates once the stopping criterion has been met and the optimal solution is chosen as the global best particle at the end of that iteration. The tasks involved in this algorithm are further clarified in the following steps:

1. Produce a random swarm of particles $S\left({\mathbf{x}}_0,{\mathbf{x}}_1,\dots ,{\mathbf{x}}_{\mathbf{N}}\right)$ .

2. Define the global best particle as the fittest particle in the swarm and let the personal best particles be the same as the initial population.

3. If stopping criterion is met, i.e. K_max iterations are performed, advance to step 7, or else carry on to next step.

4. Calculate the new set of positions $\mathbf{x}(t+1)$ for the various particles using equations (36) and (37).

5. For each particle perform the following actions: If the fitness of the newly computed particle is better than its personal best location, then replace the personal best particle’s location with those of the newly computed particle as described by equation (38).

6. Register the best particle in the personal best matrix as the global best particle using equation (39).

7. Use the global best particle as the optimal solution.

We have used $w_1=0.5,w_2=2,w_3=2$ , k_max=100 and N = 100.

CRS optimization algorithm

First, an initial population $S\left\{{\mathbf{x}}_0,{\mathbf{x}}_1,\dots {\mathbf{x}}_{10\mathbf{n}}\right\}$ is generated as in the case of DE and PSO algorithms. However, the size of this population is set to be exactly 10 times larger than the number of variables present in the problem. Each iteration consists of several steps, the first of which involves randomly selecting n distinct individuals $\{{\mathbf{v}}_1,{\mathbf{v}}_2,\dots {\mathbf{v}}_{\mathbf{n}+1}\}$ from the population set, ${\mathbf{v}}_1$ being the best individual in S. Next, the centre of gravity G of the first n selected individuals $\{{\mathbf{v}}_1,{\mathbf{v}}_2,\dots {\mathbf{v}}_{\mathbf{n}}\}$ is computed. Thereafter, a candidate individual y is calculated based on G and the individual ${\mathbf{v}}_{n+1}$ that was earlier chosen from the population set and not incorporated in calculating G. The governing equation for the candidate individual is as follows

$$\mathbf{y}=2\mathbf{G}-{\mathbf{v}}_{n+1}$$

(40)

This candidate solution y is accepted and replaces the weakest individual in the population if and only if its fitness is superior to that of the weakest individual ${\mathbf{x}}_w$ in the population. This reasoning can be modelled as follows

$${\mathbf{x}}_w=\mathit{min}(\mathbf{S})$$

(41)

$${\mathbf{x}}_w=\{\begin{array}{l}\mathbf{y}\text{ if } f\left(\mathbf{y}\right)\le f\left({\mathbf{x}}_w\right)\hfill \\ {\mathbf{x}}_w\text{otherwise}\hfill \end{array}$$

(42)

After the above steps are completed, the procedure is continued until the stopping criterion is met. ⁴⁶ The procedure for CRS global optimization algorithm is summarized in the following steps:

Step 1 Generate a randomly distributed population set that uniformly spans the search space.

Step 2 If the stopping criterion is met, then proceed to step 7; otherwise continue to the next step.

Step 3 Randomly select n distinct individuals from the population set S, ${\mathbf{v}}_1$ being the current best individual.

Step 4 Compute the median for each problem variable from the first n individuals that were chosen in the previous step.

Step 5 Formulate a candidate individual y according to equation (40).

Step 6 Replace the weakest individual in the population ${\mathbf{x}}_w$ with the candidate individual y if the y incurs a lowest cost than ${\mathbf{x}}_w$ , and thereafter return to step 2.

Step 7 Select the individual with the best fitness value as the optimal solution.

We have used n = 60, k_max = 6000 and N = 240.

MPSO

To explain the modification and its resulting impact in comparison to the PSO, the analogy of a swarm of particles in 2D search space with two variables is presented. Further study of both the PSO governing equations presented in equations (36) and (37), several inferences can be made regarding convergence basis of the local search.

According to the PSO equations, each particle is programmed to search for the optimal solution in the vicinity of its personal best solution, in the region around the particle in the solution space which has the best solution, whilst subjected to the momentum of its velocity from the previous iteration. The resultant position vector is denoted as ${\mathbf{x}}_i(t+1)$ .

The resultant vector and hence the resulting particle position may be considerably swayed by the local search vector based on the ratio $\frac{\mathit{rand}1}{\mathit{rand}2}$ . Moreover, if the personal best position of the particle happens to be poor, the resulting position of the particle may fall even further away from the optimal solution. Such a scenario is highly possible as early iterations have shown poor convergence. ⁴²

In accordance with this reasoning, it would be appropriate to alter equation (37) by replacing the personal best position ${\mathbf{x}}_i(t)$ of the particle of interest with personal position ${\mathbf{P}}_i(t)$ of any particle in the solution space which has a superior personal best result than that of the particle of interest. In light of these alterations, equation (37) now becomes

$${\mathbf{V}}_i\left(t+1\right)=w_1{\mathbf{V}}_i\left(t\right)+w_{2}r\mathit{and}1\left({\mathbf{P}}_i\left(t\right)-{\mathbf{x}}_i\left(t\right)\right)+w_{3}r\mathit{and}2\left({\mathbf{G}}_{\mathit{best}}-{\mathbf{x}}_i\left(t\right)\right)$$

(43)

where ${\mathbf{V}}_i$ and ${\mathbf{x}}_i$ are the velocity and position of the ith particle in the solution space; ${\mathbf{P}}_i(t)$ is the personal best position of any particle in the solution space whose personal best position is fitter than the ith particle, $f\left({\mathbf{P}}_i\left(t\right))<f({\mathbf{x}}_i\left(t\right)\right)$ . Such a modification would tend to bend the resultant vector closer towards the optimal solution and hence improve convergence. The algorithm is summarized in the following steps:

We have used $w_1=0.5,w_2=2,w_3=2$ , k_max = 100 and N = 100.

MCRS

In the fundamental CRS equations presented in equations (40) to (42) the centre of gravity G has a major impact on the candidate solution y. As G is primarily dependent on the n selected individuals $({\mathbf{v}}_1,{\mathbf{v}}_2,\dots {\mathbf{v}}_{\mathbf{n}})$ , its value may become fixed and hence the routine will lack flexibility if the individuals in the solution space become relatively cluttered. Consequently, this will lead to early convergence and limit the success rate (number of times the weakest individual ${\mathbf{x}}_w$ is replaced) of the algorithm.

To overcome this shortfall, three random individuals will be selected from the solution space S, as opposed to the n + 1 individuals that were previously chosen, these individuals will be in ascending order according to their fitness values with ${\mathbf{x}}_1$ being the fittest individual followed by ${\mathbf{x}}_2$ and ${\mathbf{x}}_3$ , respectively. G will correspond to the mean of ${\mathbf{x}}_1$ and ${\mathbf{x}}_2$ , and the candidate solution y will be computed as follows

$$\mathbf{y}=2\mathbf{G}-{\mathbf{x}}_3$$

(44)

By doing so, the flexibility of the G and hence the flexibility in the candidate solution will improve. Furthermore, the weakest of the randomly selected individual ${\mathbf{x}}_3$ shows greater potential in improving. The candidate solution y is driven closer to the fitter solutions of ${\mathbf{x}}_1$ and ${\mathbf{x}}_2$ and may randomly fall closer to the optimal solution O. Hence, the success rate of the algorithm will improve. This MCRS algorithm is summarized in the following steps:

Step 1 Generate a randomly distributed population set S that uniformly spans the search space.

Step 2 Check if the stopping criterion is met and if not continue to the next step, otherwise proceed to step 7.

Step 3 Randomly select three distinct individuals from the population set S.

Step 4 Order these individuals in ascending order according to their fitness values with ${\mathbf{x}}_1$ being the fittest individual followed by ${\mathbf{x}}_2$ and ${\mathbf{x}}_3$ , respectively.

Step 5 Set the centre of gravity G to be the mean of the two fittest individuals ${\mathbf{x}}_1$ and ${\mathbf{x}}_2$ .

Step 6 Determine the trial individual using equation (44).

Step 7 Replace the weakest individual in the population ${\mathbf{x}}_w$ with the candidate individual y; if the candidate solution y incurs a lower cost than ${\mathbf{x}}_w$ , then return to step 2.

Step 8 Select the individual with the best fitness value as the optimal solution.

We have used n = 3, k_max = 6000 and N = 240.

Analysis of the convergence of the various routines

The convergence histories of the fitness value using the proposed optimal routines are plotted in Figures 4 and 5. The MCRS and MPSO algorithms outperformed their predecessors with improved fitness value and quicker convergence. Hence, it may be concluded that the suggested modifications made improved the respective policies, with the CRS becoming more flexible and the PSO showing better convergence of weaker particles.

Performance of the CRS routines may also be evaluated in terms of success rate, which is measured as the ratio of how often the weakest individual in the population ${\mathbf{x}}_w$ is replaced

$$\text{success} \text{rate}=\frac{\text{number} \text{of} \text{times} {\mathbf{x}}_w \text{replaced}}{\text{number} \text{of} \text{iterations}}$$

(45)

The resulting success rate for the CRS was 0.075 and that of the MCRS was 0.15. This infers that the MCRS produces more improved solutions than the CRS.

Figures 4 and 5 show that the suggested modifications in CRS and PSO routines yielded superior results. Altering the local search characteristics of each particle in the PSO method improved the effectiveness and efficiency of the algorithm. It may also be concluded that increasing the flexibility of CRS using the proposed method outlined in ‘CRS optimization algorithm’ section added value to the algorithm with a better success rate and prevented early convergence and produced a better resulting fitness value. The DE algorithm produced the best results followed by MPSO and MCRS, respectively. The controller gains computed by these superior algorithms in addition to the manually tuned controller are listed in Tables 2 and 3. The inner-loop derivative controllers’ gains k_d are of a small order of magnitude and may be ignored for practical purposes.

Figure 5.

Convergence history plots for CRS and MCRS. CRS: controlled random search; MCRS: modified controlled random search.

Figure 6.

Bar graphs depicting the variation in suspension travel for each control law. DE: differential evolution; MCRS: modified controlled random search; MPSO: modified particle swarm optimization; RMS: root-mean-square.

Table 2.

Gains computed using the various optimization algorithms for front suspensions of the PID-controlled AVSS.

	Outer PID loop gains			Inner PID loop gains
Technique	KP	KI	KD	kp	ki	kd
Front right suspension system
Manual	1100	360	140	0.002	0.001	0
DE	1692	267	166	0.0038	0.0010	3 × 10−9
MCRS	17,251	−145	−92	0.0029	0.0700	5 × 10−9
MPSO	7270	320	−769	0.0061	−0.00063	3 × 10−9
Front left suspension system
Manual	1050	170	220	0.002	0.001	0
DE	1692	267	166	0.0038	0.0010	3 × 10−9
MCRS	16,003	184	158	0.0028	0.0665	4 × 10−9
MPSO	5025	−1250	−1350	0.0080	0.0024	3 × 10−9

AVSS: active vehicle suspension system; DE: differential evolution; MCRS: modified controlled random search; MPSO: modified particle swarm optimization; PID: proportional + integral + derivative.

Table 3.

Gains computed using the various optimization algorithms for rear suspensions of the PID-controlled AVSS.

	Outer PID loop gains			Inner PID loop gains
Technique	KP	KI	KD	kp	ki	kd
Rear right suspension system
Manual	1200	340	150	0.002	0.001	0
DE	1692	267	166	0.0038	0.0010	3 × 10−9
MCRS	6723	−87	19	0.0020	0.0182	3 × 10−9
MPSO	8566	−1274	−5	0.0033	0.0046	3 × 10−9
Rear left suspension system
Manual	1000	200	200	0.002	0.001	0
DE	1692	267	166	0.0038	0.0010	3 × 10−9
MCRS	4579	302	178	0.0030	0.0233	2 × 10−9
MPSO	–788	−372	950	0.0055	−0.0002	6 × 10−9

AVSS: active vehicle suspension system; DE: differential evolution; MCRS: modified controlled random search; MPSO: modified particle swarm optimization; PID: proportional + integral + derivative.

Simulation results and discussion

The evolutions of the performance index through the use of the proposed tuning algorithms presented in Figures 4 and 5 show that the DE routine gave the best performance index followed by the MPSO and MCRS algorithms, respectively. However, these plots cannot provide information on how well the suspension trade-offs have been resolved. Hence plots for each suspension performance criterion will be plotted for the non-optimized, DE, MCRS and MPSO cases, respectively.

Simulations were performed in the Matlab/Simulink environment and the simulation time was set to 5 s. Suspension travel performance, road holding (tyre dynamic load) and power consumption (control input voltage) are plotted in Figures 10 to 15, respectively. These plots focus on the rear left suspension corner only as the worst performance behaviour was observed at this location. Vehicle handling (roll and pitch accelerations) is presented in Figures 12 and 13, respectively. Ride comfort (body-heave acceleration) is shown in Figure 14 and the cumulative hydraulic force applied to the vehicle body is plotted in Figure 16. The peak and RMS values pertaining to these quantities are compared in Figures 6 to 9.

The suspension travel response obtained when implementing each of the optimal routines displayed reduced peak and RMS values in comparison with PVSSs and manually tuned PID cases, respectively. They showed better transient behaviour by damping out with no further peaks immediately after the road disturbance was removed. These results were anticipated as the performance index did indeed address suspension travel with a fair and considerable weighting factor.

The road holding capabilities acquired for DE and MCRS cases were similar and superior to those of the PVSSs and manually tuned cases with lower peaks, better RMS values, quicker settling times and improved transient response that had fewer peaks and reduced oscillations. On the other hand, the MPSO case did manage to improve the RMS value, transient behaviour and settling time, but still produced the largest peak values. However, even with the inter-relationship between RMS and peak values, it is not guaranteed that reducing the RMS value will always reduce the peak values and hence such a shortfall is possible.

Regarding vehicle handling and ride comfort, the optimal policies produced lower peak and RMS values with quicker settling times than those of the manually tuned AVSS and PVSSs. However, each policy contained a greater degree of chattering which would tend to deteriorate system components. In terms of ride comfort and handling performance, the MPSO case was the worst from the optimal methods followed by the MCRS and DE cases, respectively. Its weak performance is due to the fact that this routine produced a lower performance index value as compared to the DE case. It is worth noting that although MCRS had a similar performance index, it did perform adequately in vehicle handling, but lacked quality in suspension travel. In conclusion, both algorithms produced a weaker performance index than DE, but they did exhibit desired responses in certain performance aspects whilst performing weaker in other aspects. This implies that these algorithms were not as good as the DE in resolving the trade-offs between the various performance criteria. Furthermore, DE produced the best performance index which was approximately 40% better than its counterparts.

The comparative plot pertaining to control input voltage and cumulative hydraulic force showed that all the optimal tuning policies produced lower peak and RMS values than the manually tuned case whilst at the same time they were able to enhance the RMS values and response of ride comfort, road holding and vehicle handling criteria. From a computational standpoint, this is projected as control input voltage and supplied hydraulic forces are substantial factors of the performance index. However, from an engineering point of view, such data are rather contradictory to both typical quarter-car models and linear control techniques as a larger force or voltage is demanded to improve the various performance benchmarks. These results infer that the coupling and nonlinearities of the full-car nonlinear system is a major factor that supplements the outcomes of the system and it is thus imperative that they are thoroughly investigated.

Frequency-domain analysis

The ISO 2631 ⁴⁴ states that human exposure to frequencies ranging from 0.5 to 80 Hz significantly affects human comfort. Thus, it is imperative that relevant AVSS performance criteria be analysed in this range. Bode plots pertaining to ride comfort, vehicle handling and road holding for the PVSSs, manually tuned PID and DE-optimized PID cases are illustrated in Figures 17 to 20. These plots were generated using the power spectral density estimates based on Welch algorithm in the Matlab/Simulink signal processing toolbox. The following parameters are used in computing the Welch’s periodograms: the windowing function–Hanning window function; the number of points used in forming each fast Fourier transform, NFFT $=1024$ ; length of the window, NWind $=256$ ; and the sampling frequency of the windows was set at 80 Hz to accommodate the whole body vibration range.

Figure 7.

Bar graphs illustrating the difference in control input voltage for the proposed controllers. DE: differential evolution; MCRS: modified controlled random search; MPSO: modified particle swarm optimization; RMS: root-mean-square.

Figure 8.

Bar graphs summarizing the road holding aspect for each control case. DE: differential evolution; MCRS: modified controlled random search; MPSO: modified particle swarm optimization; RMS: root-mean-square.

Figure 9.

Quantitative information pertaining to vehicle handling, ride comfort and actuator force supplied. DE: differential evolution; MCRS: modified controlled random search; MPSO: modified particle swarm optimization; RMS: root-mean-square.

Figure 10.

Illustration of the variation in suspension travel response using the proposed tuning methods. DE: differential evolution; MCRS: modified controlled random search; MPSO: modified particle swarm optimization; PID: proportional + integral + derivative.

Figure 11.

Variation in the tyre dynamic load experienced at the rear left suspension system for the various tuning policies. DE: differential evolution; MCRS: modified controlled random search; MPSO: modified particle swarm optimization; PID: proportional + integral + derivative.

Figure 12.

Vehicle body pitch acceleration for each of the tuning routines. DE: differential evolution; MCRS: modified controlled random search; MPSO: modified particle swarm optimization; PID: proportional + integral + derivative.

Figure 13.

Vehicle body roll acceleration for the various tuning policies. DE: differential evolution; MCRS: modified controlled random search; MPSO: modified particle swarm optimization; PID: proportional + integral + derivative.

Figure 14.

Ride comfort experienced for each of the recommended tuning methods. DE: differential evolution; MCRS: modified controlled random search; MPSO: modified particle swarm optimization; PID: proportional + integral + derivative.

Figure 15.

Difference in control input voltage produced using the suggested tuning approaches. DE: differential evolution; MCRS: modified controlled random search; MPSO: modified particle swarm optimization; PID: proportional + integral + derivative.

Figure 16.

Effective hydraulic force applied to the vehicle chassis for the various tuning algorithms. DE: differential evolution; MCRS: modified controlled random search; MPSO: modified particle swarm optimization; PID: proportional + integral + derivative.

Figure 17.

Frequency response comparative plot for the body-heave acceleration. DE: differential evolution; PID: proportional + integral + derivative.

Figure 18.

Comparative plot of frequency response pertaining to body pitch acceleration. DE: differential evolution; PID: proportional + integral + derivative.

Figure 19.

Roll acceleration frequency response comparative plot for relevant cases. DE: differential evolution; PID: proportional + integral + derivative.

Figure 20.

Frequency response for road holding. DE: differential evolution; PID: proportional + integral + derivative.

In terms of ride comfort, the PVSSs effectively behaved as a high-pass filter, whose high frequency signals were successfully attenuated. It possessed the worst attenuation at the onset where it had a magnitude of approximately 1 in the range of 0.01–1 Hz. On the contrary, the AVSS methods produced a significant improvement in this range as they were able to successfully attenuate signals therein, with the optimal controller performing best. Thereafter, the PVSS case possessed superior attenuation properties up until 2 Hz. Afterwards, the AVSS cases performed similarly and produced better attenuation results up to 10 Hz, where a resonant peak developed.

In relation to pitch acceleration, the PVSSs performed the worst at the onset which is by the standard regarded by the European Commission ⁴⁷ as the most sensitive frequencies experienced by humans. The AVSS cases showed an improvement in comparison with the PVSSs in this regard, with the optimal case performing the best. In this range the power ratio did increase with increasing frequency but was still able to perform better than the PVSS case. At the onset of 1 Hz the power ratios of each case dropped off sharply inferring successful signal attenuation. The AVSS cases did however experience a resonance peak later at around 12 Hz, but the power ratios were comparatively low at these high frequencies, which implies successful attenuation.

With regards to roll acceleration, there was only a marginal variation between the PID-controlled and PVSS cases at the sensitive low frequencies. The DE-optimized case performed significantly better than its predecessors in this range but it produced an increasingly higher power ratio with increasing frequency. Thereafter, at around 1 Hz it began to perform worse than the other cases. From 11 Hz onwards the power ratio dropped off substantially indicating that the signals were henceforth attenuated considerably. A resonant peak did however develop for the AVSS cases at 10.5 Hz but the power ratio was low enough which saw the continued attenuation of input signals.

Road holding frequency plots resemble that of a high-pass filter where the sensitive low frequencies produced the largest peaks and the less sensitive high frequency signals were successfully attenuated. The PID-controlled and PVSS cases produced similar results for the whole range of frequencies where the optimal DE case performed significantly better. The continued superior performance of the optimal PID controller in relation to its counterparts for all suspension criteria in the sensitive low frequency ranges suggests two things. First, it was more successful than both its predecessor in resolving the conflicting trade-offs in ride comfort, road holding and vehicle handling. Second, it was successful in managing the sensitive low frequency signals for all suspension criteria. These results imply that controller tuning through global optimization methods does play a significant role in both resolving AVSS trade-offs and improving robustness to parameter variations.

Conclusion and future work

This paper presented applications of global and evolutionary optimization algorithms (CRS, DE, PSO, MCRS and MPSO) to tune the gains of multi-loop PID controllers for full-car nonlinear electrohydraulic AVSS to resolve the conflicting performance criteria. Multiple control loops were formulated to both regulate the controlled variable and maintain actuator integrity (stability).

The DE algorithm proved to be the most consistent as it was the only algorithm that maintained a desirable compromise and satisfactory performance in ride comfort, vehicle handling, suspension travel, road holding and power consumption. System performance through optimization achieved comparatively superior results for more severe disturbance than those that were reported in previous studies. Even though MPSO and MCRS algorithms significantly outperformed their predecessors, they did not however achieve the effectiveness of the DE algorithm. In the frequency domain, the DE still managed to maintain an edge in ride comfort and road holding over its manually tuned and PVSS counterparts. The successful performance of DE algorithm was due to the structure of its algorithm, which permits efficient exploration of the search space, and through certain acceptance conditions only replaces weaker solutions if and only if better solutions are found.

Further works must be extended to fault tolerant controller design for active suspension systems and experimental validation. Real world models contain more complexities that are ignored in the course of numerical simulations.