Nonlinear model predictive control for yaw rate and body motion control through semi-active and active suspensions

Abstract

Active and semi-active suspensions for passenger cars traditionally enhance comfort through body control, and vehicle handling by reducing the tyre load variations induced by road irregularities. Active suspensions can also be designed to track a desired yaw rate profile through the control of the lateral load transfer distribution between the front and rear axles. This paper considers an integrated system including semi-active and active suspension actuation to control the yaw, roll, pitch and heave dynamics excited by the driving actions. To this purpose, two novel real-time-capable implicit nonlinear model predictive control (NMPC) formulations, excluding and including cost function weight adaptation, are proposed and compared with the passive vehicle, and the controlled vehicle with two combinations of skyhook and active roll control, the first based on a pseudoinverse decoupling transformation for obtaining the damping force contributions, and the second using an inverse formulation. The algorithms are assessed through an experimentally validated simulation model, along manoeuvres corresponding to sub-limit and limit handling operation, to analyse the trade-off between body motion reduction and cornering response enhancement. The results show that the adaptable NMPC configuration provides the best performance in all scenarios, also for significant variations of the main vehicle and tyre parameters.

KEYWORDS:

• Active suspension
• yaw rate control
• body motion control
• nonlinear model predictive control
• inverse
• pseudoinverse

1. Introduction

Semi-active and active suspensions are frequently adopted in modern cars to improve comfort and handling. Such systems could become more widely spread in the next generations of automated vehicles, to enable the users to comfortably carry out other activities during vehicle operation [1].

Controllable suspensions are based on actuators between the sprung and unsprung masses, and thus influence the motion of the sprung mass and the vertical tyre force distribution among the vehicle corners. Differently from active suspensions, semi-active actuators cannot introduce energy into the system, i.e. they are implemented in the form of shock absorbers with controllable damping characteristics, and therefore are effective only during transients. Semi-active implementations are more frequent than active suspensions because they represent a balanced compromise in terms of cost, component weight, power consumption, and performance [2,3]. However, active suspensions offer significantly enhanced functionality, including the capability of: a) shaping the level of vehicle understeer in steady-state conditions, and b) increasing yaw and sideslip damping during transients [4].

This paper proposes and compares controllers for a next-generation hydraulic suspension system set-up by Tenneco Automotive. The hardware architecture enables independent damping force control on each vehicle corner, and active anti-roll moment control at the axle level.

The research contributions are:

• A novel real-time-capable nonlinear model predictive control (NMPC) implementation for the case study integrated semi-active and active suspension system. The algorithm uses a five-degree-of-freedom internal model considering tyre nonlinearities, jacking forces, and anti-features of the suspension geometry, which enable meaningful model-based control of the pitch and heave dynamics induced by the driving actions. The NMPC benefits from the inclusion of a cost function weight adaptation mechanism to vary the control priorities based on the driving conditions. The NMPC architecture is modular, i.e. it can be interfaced with external control functions, in this specific case a skyhook controller targeting the compensation of the effect of road irregularities.

• A performance comparison of the proposed NMPCs, excluding and including weight adaptation, with benchmarking active roll moment and skyhook damping controllers, referred to as pseudo-inverse and inverse, depending on the involved decoupling transformation used to calculate the reference damping forces in the skyhook equations (as detailed in Section 4.5). The assessment is carried out through simulations with an experimentally validated vehicle model.

The remainder is organised as follows: Section 2 summarises the literature, and identifies the knowledge gap addressed by this study; Section 3 presents the simulation framework and control architecture; Section 4 describes the novel NMPC and benchmarking implementations; Section 5 discusses the simulated manoeuvres and considered performance indicators; Section 6 analyses the simulation results; finally, Section 7 summarises the main conclusions.

2. Literature review

Controllable suspensions are widely discussed in the literature, mostly to: i. mitigate the motion of the sprung mass caused by road irregularities; ii. reduce the vertical tyre load oscillations on irregular surfaces, and thus improve road holding; and iii. compensate the heave, pitch and roll motions of the sprung mass, excited by the driving actions, i.e. by traction/braking and cornering.

In many implementations targeting i.-iii., e.g. [5–11], a high-level algorithm calculates the total reference heave force at the centre of gravity, as well as the reference anti-pitch and anti-roll moments for the sprung mass. These are then converted into individual forces for the four actuators at the vehicle corners, through the so-called decoupling transformation, often based on a pseudo-inverse control allocation matrix. To achieve iii., the integrated chassis controller (ICC) in [12] uses a low-level control allocation algorithm that accounts for the anti-properties of the suspension system. Body motion control can also be implemented through algorithms directly generating the actuator control inputs, see the roll controllers in [13–15], the linear quadratic regulator (LQR) for heave, pitch and roll control in [16], and the LQRs and robust roll controllers in [17].

In [4] and [18–33], the controlled suspensions enhance the vehicle yaw rate response (objective iv. of controllable suspension systems) through variable front-to-total anti-roll moment distribution. In fact, for a given lateral acceleration, increased front anti-roll moment and decreased rear anti-roll moment increase the magnitude of the front slip angle and reduce the rear one, thus making the vehicle more understeering. On the one hand, because of the nonlinearity of the phenomenon, many front-to-total anti-roll moment distribution controllers are not model-based, see [18–26]. On the other hand, most of the available model-based anti-roll moment distribution designs use simplified vehicle models including lateral and yaw dynamics, where the lateral tyre force is linearly varying with the slip angle, and the cornering stiffness is parabolically dependent on the vertical tyre load, e.g. see [27–32]. Reference [33] combines the parabolic cornering stiffness variation with vertical load and a simplified version of the Pacejka magic formula. In [4] and [34] Ricco et al. highlight the significant limitations of the model in [27], propose a linearised model for front-to-total anti-roll moment distribution design in the frequency domain, and highlight the resulting performance benefits in limit handling conditions.

A recent trend is the implementation of model predictive controllers (MPCs) for suspension control, whose real-time operation is made possible by the increased computational capability of automotive microcontrollers, and the effectiveness of real-time optimisation solvers [35]. For example, the MPC in [36] addresses the heave, pitch and roll dynamics of the sprung mass induced by unknown road disturbances, while in [37] Shao et al. evaluate a distributed MPC for heave and pitch control on irregular roads. The road profile is considered a known input in the road preview controllers in [38–44], which target i., and, depending on the application, the improvement of other aspects, e.g. ii., while considering actuation effort and constraints. However, these MPCs neglect the vehicle cornering dynamics, and do not include tangential tyre force formulations in their prediction models. In [45] Zhu et al. propose a cascade control structure, with an upper level MPC that uses the simplified tyre model from [27] to reduce the roll motion and improve the yaw rate response. In [46] Adireddy et al. present an MPC-based ICC with an 8-degree-of-freedom prediction model, which neglects the heave and pitch dynamics, similarly to the real-time NMPC in [47]. Reference [48] shows the potential of linear time-varying MPC for ICC in the linear cornering response region, with maximum lateral accelerations of ∼4 m/s², in which controlled suspensions have very limited impact on vehicle dynamics. Similarly, the NMPC-based ICC in [49] is tested at low lateral accelerations, and the active suspension system enhances only ride comfort.

In conclusion, a real-time-capable MPC architecture that integrates the compensation of the sprung mass motions caused by the driving actions, ride comfort enhancement, and yaw rate control for operation at the limit of handling is still missing. Moreover, the available prediction model formulations for front-to-total anti-roll moment distribution and vehicle dynamics control neglect the anti-lift, -squat and -dive properties of the suspension system, as well as the jacking force effects. The identified gaps will be addressed in the remainder.

3. Case study vehicle and simulation framework

3.1. Reference vehicle and suspension actuation system

The case study vehicle is an Audi E-Tron prototype (see its main parameters in Table 1) equipped with: i) two on-board centralised electric motors, one per axle, each of them connected to the wheels through a single-speed transmission, mechanical differential, half-shafts, and constant velocity joints; and ii) the novel integrated CVSA2-Kinetic suspension system by Tenneco [50], including four hydraulically interconnected dampers, i.e. one per corner, according to the CVSA2 (continuously variable semi-active suspension solution) damping set-up, in which the hydraulic dampers are equipped with two external electro-hydraulic valves that independently control the damping characteristics for the rebound and compression motions, thus providing a semi-active contribution. The system is upgraded with a Kinetic roll control system, consisting of hydraulic lines between the dampers, replacing the mechanical anti-roll bars, with a pressure control unit including a pump and valves that are used to add an active contribution, to achieve the desired front and rear anti-roll moments. As a consequence, the force of each damper is given by the superposition of the independent damping from CVSA2 and the active contributions from the pump of the Kinetic system. The suspension controller therefore generates six inputs, i.e. four control currents for semi-active control, and two reference anti-roll moments for active roll control. The active anti-roll moments are expressed through active force contributions, which are constrained by the hydraulic arrangement to be equal and opposite for the hydraulic dampers on the same axle. The reference active forces are sent to the low-level controller of the actuation system, developed by Tenneco, which generates the expected actuation for the involved valves and pump, and whose detailed operation is beyond the scope of this study.

Table 1 Main vehicle parameters.

Description	Parameter	Value
Vehicle mass	$m$	2843 kg
Sprung mass	$m_s$	2593 kg
Yaw mass moment of inertia	$I_z$	5291 kg m^2
Roll mass moment of inertia of the sprung mass	$I_x$	550 kg m^2
Pitch mass moment of inertia of the sprung mass	$I_y$	2200 kg m^2
Front semi-wheelbase	$a_F$	1.47 m
Rear semi-wheelbase	$a_R$	1.46 m
Front and rear track widths	$t_F\approx t_R\approx t$	1.66 m
Centre of gravity height	$h$	0.63 m
Wheel radius	$R$	0.38 m

3.2. Simulation framework

The simulation framework, see Figure 1, consists of the following functional blocks:

• A virtual driver model, which generates the accelerator pedal position, $p_{\mathrm{acc}}$, brake pedal position, $p_{\mathrm{bk}}$, and steering wheel angle, ${\delta }_{\mathrm{sw}}$. For tracking the reference speed profiles, proportional integral (PI) controllers are used to output the desired pedal positions, while the ${\delta }_{\mathrm{sw}}$ profiles are directly generated with fixed look-up tables of ${\delta }_{\mathrm{sw}}$ as a function of time, without a feedback contribution.

• A drivability layer, which outputs the individual electric motor and braking torque levels, $T_{m,i}$ and $T_{b,\mathrm{ij}}$, where the subscript $i=F,R$ indicates the front or rear axles, and the subscript $j=L,R$ indicates the left or right sides.

• A reference generator, which defines the reference yaw, roll, pitch and heave rates (respectively ${\dot{\psi }}_{\mathrm{ref}}$, see [4] for its definition, ${\dot{\phi }}_{\mathrm{ref}}$, ${\dot{\theta }}_{\mathrm{ref}}$, and ${\dot{z}}_{\mathrm{hv},\mathrm{ref}}$), as well the reference anti-roll moment, anti-pitch moment, and anti-heave force ($M_{\phi }$, $M_{\theta }$, and $F_{\mathrm{hv}}$), starting from the driver steering input, ${\delta }_{\mathrm{sw}}$, and the measured/estimated signals, such as vehicle speed, $V$, and the variables in the state vector, $x$.

• A vehicle stability controller (VSC) [51,52], based on the combination of heuristic rules and proportional integral derivative (PID) controllers. Emergency conditions are identified through the magnitude of the yaw rate error, i.e. the VSC is activated based on thresholds that are different depending on the understeering or oversteering condition of the vehicle, computed with respect to (w.r.t.) ${\dot{\psi }}_{\mathrm{ref}}$. The VSC algorithm specifies the: a) variation of the total longitudinal tyre force, e.g. the powertrain torque is brought down to zero independently from the driver input on the accelerator pedal, and the brake force demand can be increased or reduced w.r.t. the driver request; and/or b) generation of a direct yaw moment, through PID controllers with different tunings depending on the understeering or oversteering condition of the vehicle, according to the current industrial practice in VSC design. The PID outputs are saturated – including appropriate anti-windup – based on the estimated friction level on the individual corners, see [53]. a) and b) are achieved through the actuation of the friction brakes, where b) is generated mainly on the inner rear or outer front corners, up to their respective saturation levels. Once tyre force saturation is reached, also the other tyre located on the same vehicle side is actuated. The VSC outputs are the modified motor and braking torque levels, $T_{m,\mathrm{VSC},i}$ and $T_{b,\mathrm{VSC},\mathrm{ij}}$, the latter to be actuated by a brake-by-wire system.

• The proposed NMPC (see Sections 4.1–4.4), which outputs the reference value of the variation $\mathrm{\Delta }f$ of the front-to-total anti-roll moment distribution ratio for the active part of the suspension system w.r.t. its nominal value, $f_{\mathrm{nom}}$, as well as the current levels, ${\epsilon }_{\mathrm{NMPC},\mathrm{ij}}$, for the valves of the controllable dampers. The ${\epsilon }_{\mathrm{NMPC},\mathrm{ij}}$ values are computed from the reference currents, ${\epsilon }_{\mathrm{SH},\mathrm{ij}}$, generated by the benchmarking inverse skyhook algorithm, targeting ride comfort enhancement. Further damping control functions, external to the NMPC, can be integrated in the architecture, by summing the respective reference currents before they are sent to the NMPC. The ${\epsilon }_{\mathrm{NMPC},\mathrm{ij}}$ levels deviate from the externally generated values, to consider the NMPC objectives in terms of cornering and body motion control, while embedding actuation constraints.

• The front-to-total anti-roll moment distribution ratio function, which calculates the total active and semi-active anti-roll moment distribution factors, $\mathrm{\Delta }f$ and $\lambda$, for the benchmarking inverse controller.

• The benchmarking pseudo-inverse and inverse damping force controllers (see Section 4.5), which generate the reference damping forces, $F_{d,\mathrm{ij}}$, starting from $M_{\phi }$, $M_{\theta }$, $F_{\mathrm{hv}}$, and $\lambda$.

• The actuation algorithm for the benchmarking skyhook damping controllers, developed by the involved suspension component supplier, which calculates ${\epsilon }_{\mathrm{ij},\mathrm{SH}}$ starting from $F_{d,\mathrm{ij}}$ and the measured actuator speeds.

• The active suspension force distributor, which calculates the total active anti-roll moment starting from the lateral acceleration $a_y$, and, by combining it with the $\mathrm{\Delta }f$ value from the NMPC or the benchmarking inverse controller, generates the reference active suspension forces, $F_{\mathrm{act},i}$.

• The high-fidelity vehicle simulation model for control system assessment, implemented via the software package VSM by AVL [54], including consideration of the actuation dynamics of the semi-active and active suspension components, which are approximated through first order models with 20 and 35 ms time constants. The model was experimentally validated through vehicle tests executed in the framework of the European project EVC1000 [55]. Figure 2 reports examples of comparisons of experimental measurements with the model outputs in steady-state and transient cornering conditions (see Figure 3 for the definition of the main variables). The good match makes the model a valid tool for control performance assessment.

Figure 1 Simplified schematic of the simulation environment. The blocks with white background refer to the proposed suspension controllers; the blocks with grey background refer to the vehicle simulation model, and the powertrain and friction brake controllers.

Figure 2 Vehicle model validation results. Top: skidpad test results; bottom: transient steering manoeuvre at approximately constant torque demand, from an initial vehicle speed of ∼100 km/h. ‘Experimental’: vehicle measurements; ‘Simulation’: results from the model for control system assessment; ‘Internal’: results from the NMPC prediction model (see Section 4) running in open loop along the same test.

The block diagram in Figure 1 also includes the switches for activating/de-activating the NMPC, where in the latter case only the benchmarking algorithms are operational.

4. Controller formulations

4.1. NMPC internal model

Model concept

Figure 3 shows the top, rear and side views of the vehicle prediction model concept, also referred to as internal model, with indication of the main variables and sign conventions. The model formulation specifically targets suspension control applications, and considers the lateral and yaw vehicle dynamics, the heave, roll and pitch dynamics of the sprung mass, as well as the anti-features and jacking force effects of the front and rear suspensions. Since the involved industrial company is interested in suspension control only (typical suspension controllers operate independently from VSC systems), and expressed the requirement of limiting the computational load associated with the implicit NMPC, the longitudinal vehicle dynamics and rotational wheel dynamics have been neglected in the prediction model, although they are present in the high fidelity simulation model for control assessment. These dynamics could be included in the internal model of a model predictive controller for vehicle stability control based on the actuation of the friction brakes, e.g. see [56], which is not the topic of this manuscript. Similarly, as the NMPC addresses the body motions induced by the driving actions and yaw rate tracking, and works together with a skyhook algorithm (see Figure 1) that compensates the effect of road irregularities, the vertical motions of the unsprung masses have not been included in the internal model. The potential benefit of the specified additional degrees of freedom will be assessed in future research, focused on integrated direct yaw moment and active/semi-active suspension control.

Figure 3 Top, rear and left side views of the vehicle, with indication of the main variables and parameters used in the internal model formulation.

In the side view, the suspension is described through an equivalent swinging arm, connected to the sprung mass by a cylindrical joint, whose position, defined by the vertical and longitudinal distances $d_i$ and $e_i$ (see details in [57]), can be directly obtained from the geometry of the specific suspension arrangement, or the relevant anti-dive/-lift/-squat percentages. For the computation of the equivalent vertical forces transmitted by the suspension links to the chassis because of the friction braking torques, on the wheel side the joint of the equivalent suspension arm is considered to be located at the centre of the tyre contact patch, as in Figure 3. On the contrary, for the computation of the vertical suspension link force related to the powertrain torque contribution, the arm is considered connected to the wheel centre, to account for the effect of the reaction torque of the on-board powertrains, which is applied to the sprung mass. Hence, the longitudinal load transfer during traction and braking is split between a component through the equivalent suspension link/s, neglected in the available MPC implementations targeting cornering response enhancement, and a component through the suspension springs and actuators. This feature enables realistic prediction of the pitch and heave dynamics induced by the longitudinal tyre forces.

Similarly, the model accounts for the jacking forces, i.e. the vertical force components induced by the lateral tyre forces through the suspension links in case of non-zero roll centre height, see [58]. The resulting force, caused by the lateral tyre force difference among the two wheels of the same axle, is applied to the roll centre of each suspension. The model neglects the jacking force effects associated with roll axis migration, which will be the object of future analyses.

Model formulation

The prediction model dynamics are described by the following force and moment balance equations, under the assumption of small steering, sideslip and roll angles:

• Lateral force balance at the vehicle level (1) $\dot{\beta }=\frac{F_{y,\mathrm{FL}}+F_{y,\mathrm{FR}}+F_{y,\mathrm{RL}}+F_{y,\mathrm{RR}}}{\mathrm{mV}}-\dot{\psi }$(1)

• Yaw moment balance at the vehicle level (2) $\ddot{\psi }=\frac{a_F[F_{y,\mathrm{FL}}+F_{y,\mathrm{FR}}]-a_R[F_{y,\mathrm{RL}}+F_{y,\mathrm{RR}}]+M_{z,\mathrm{VDC}}}{I_z}$(2)

• Roll moment balance of the sprung mass $\begin{array}{rl}\ddot{\phi }& =\frac{m_s\{V[\dot{\beta }+\dot{\psi }]+\mathrm{g\phi }\}[h-h_{\mathrm{ra}}]}{I_x}+\end{array}$ $\begin{array}{rl}& +\frac{0.5t_F[F_{b,z,\mathrm{FL}}-F_{b,z,\mathrm{FR}}+F_{\mathrm{pwt},z,\mathrm{FR}}-F_{\mathrm{pwt},z,\mathrm{FL}}]}{I_x}+\end{array}$ $\begin{array}{rl}& +\frac{0.5t_R[F_{b,z,\mathrm{RR}}-F_{b,z,\mathrm{RL}}+F_{\mathrm{pwt},z,\mathrm{RL}}-F_{\mathrm{pwt},z,\mathrm{RR}}]}{I_x}+\end{array}$ (3) $\begin{array}{rl}& +\frac{0.5t_F[F_{\mathrm{tot},z,\mathrm{FL}}-F_{\mathrm{tot},z,\mathrm{FR}}]+0.5t_R[F_{\mathrm{tot},z,\mathrm{RL}}-F_{\mathrm{tot},z,\mathrm{RR}}]}{I_x}\end{array}$(3)

• Pitch moment balance of the sprung mass $\begin{array}{rl}\ddot{\theta }& =\frac{a_R[F_{\mathrm{tot},z,\mathrm{RL}}+F_{\mathrm{tot},z,\mathrm{RR}}]-a_F[F_{\mathrm{tot},z,\mathrm{FL}}+F_{\mathrm{tot},z,\mathrm{FR}}]}{I_y}+\end{array}$ $\begin{array}{rl}& +\frac{a_R[F_{y,z,\mathrm{RR}}-F_{y,z,\mathrm{RL}}]-a_F[F_{y,z,\mathrm{FR}}-F_{y,z,\mathrm{FL}}]}{I_y}+\end{array}$ $\begin{array}{rl}& +\frac{[F_{\mathrm{pwt},z,\mathrm{FL}}+F_{\mathrm{pwt},z,\mathrm{FR}}-F_{b,z,\mathrm{FL}}-F_{b,z,\mathrm{FR}}][a_F-e_F]}{I_y}+\end{array}$ $\begin{array}{rl}& +\frac{[F_{\mathrm{pwt},z,\mathrm{RL}}+F_{\mathrm{pwt},z,\mathrm{RR}}-F_{b,z,\mathrm{RL}}-F_{b,z,\mathrm{RR}}][a_R-e_R]}{I_y}+\end{array}$ $\begin{array}{rl}& +\frac{[F_{b,\mathrm{FL}}+F_{b,\mathrm{FR}}-F_{\mathrm{pwt},\mathrm{FL}}-F_{\mathrm{pwt},\mathrm{FR}}][h-d_F]}{I_y}+\end{array}$ (4) $\begin{array}{rl}& +\frac{[F_{b,\mathrm{RL}}+F_{b,\mathrm{RR}}-F_{\mathrm{pwt},\mathrm{RL}}-F_{\mathrm{pwt},\mathrm{RR}}][h-d_R]}{I_y}\end{array}$(4)

• Heave force balance of the sprung mass $\begin{array}{rl}{\ddot{z}}_{\mathrm{hv}}& =\frac{F_{\mathrm{tot},z,\mathrm{FL}}+F_{\mathrm{tot},z,\mathrm{FR}}+F_{\mathrm{tot},z,\mathrm{RL}}+F_{\mathrm{tot},z,\mathrm{RR}}}{m_s}+\end{array}$ $\begin{array}{rl}& +\frac{F_{y,z,\mathrm{FR}}-F_{y,z,\mathrm{FL}}+F_{y,z,\mathrm{RR}}-F_{y,z,\mathrm{RL}}+F_{b,z,\mathrm{FL}}+F_{b,z,\mathrm{FR}}-F_{b,z,\mathrm{RL}}-F_{b,z,\mathrm{RR}}}{m_s}\end{array}$ (5) $\begin{array}{rl}& +\frac{-F_{\mathrm{pwt},z,\mathrm{FL}}-F_{\mathrm{pwt},z,\mathrm{FR}}+F_{\mathrm{pwt},z,\mathrm{RL}}+F_{\mathrm{pwt},z,\mathrm{RR}}}{m_s}\end{array}$(5)

where $\dot{\beta }$, $\ddot{\psi }$, $\ddot{\phi }$, $\ddot{\theta }$, and ${\ddot{z}}_{\mathrm{hv}}$ are the sideslip rate, yaw acceleration, roll acceleration, pitch acceleration, and heave acceleration; $F_{y,\mathrm{ij}}$ is the lateral tyre force at the $\mathrm{ij}$ corner; $F_{\mathrm{tot},z,\mathrm{ij}}$ is the total force through the deformable suspension components (i.e. springs and actuators) of the $\mathrm{ij}$ corner; $F_{b,\mathrm{ij}}-F_{\mathrm{pwt},\mathrm{ij}}$ is the resulting longitudinal trailing arm force applied to the sprung mass, and caused by the friction brakes and powertrains, provoking the respective vertical component $F_{b,z,\mathrm{ij}}-F_{\mathrm{pwt},z,\mathrm{ij}}$; $F_{y,z,\mathrm{ij}}$ is the corner’s jacking force; $M_{z,\mathrm{VDC}}$ is the direct yaw moment contribution associated with the VSC intervention in emergency condition (which can be considered as an external input, obtained from $T_{b,\mathrm{VSC},\mathrm{ij}}$); and $h_{\mathrm{ra}}$ is the roll axis height at the longitudinal coordinate of the centre of gravity, computed through (6), neglecting the difference between the centre of gravity position of the vehicle and its sprung mass: (6) $h_{\mathrm{ra}}\approx \frac{h_{\mathrm{rc},F}{a}_R+h_{\mathrm{rc},R}{a}_F}{l}$(6) where $l=a_F+a_R$ is the wheelbase; and $h_{\mathrm{rc},i}$ is the roll centre height of the $i$ axle.

The terms $F_{\mathrm{tot},z,\mathrm{ij}}$ are given by: $\begin{array}{rl}{F}_{\mathrm{tot},z,\mathrm{FL}}& =i_{s,F}{F}_{s,\mathrm{FL}}+i_{\mathrm{act},F}[F_{d,\mathrm{FL}}-F_{\mathrm{act},F}]\end{array}$ $\begin{array}{rl}{F}_{\mathrm{tot},z,\mathrm{FR}}& =i_{s,F}{F}_{s,\mathrm{FR}}+i_{\mathrm{act},F}[F_{d,\mathrm{FR}}+F_{\mathrm{act},F}]\end{array}$ $\begin{array}{rl}{F}_{\mathrm{tot},z,\mathrm{RL}}& =i_{s,R}{F}_{s,\mathrm{RL}}+i_{\mathrm{act},R}[F_{d,\mathrm{RL}}-F_{\mathrm{act},R}]\end{array}$ (7) $\begin{array}{rl}{F}_{\mathrm{tot},z,\mathrm{RR}}& =i_{s,R}{F}_{s,\mathrm{RR}}+i_{\mathrm{act},R}[F_{d,\mathrm{RR}}+F_{\mathrm{act},R}]\end{array}$(7) where $F_{s,\mathrm{ij}}$ is the force of the $\mathrm{ij}$ passive spring; $F_{d,\mathrm{ij}}$ is the virtual damping force of the $\mathrm{ij}$ suspension actuator; $F_{\mathrm{act},i}$ is the virtual active force provided by the actuators on the $i$ axle, considered to be equal and opposite on the two vehicle sides, to match the characteristics of the considered electro-hydraulic hardware, and to generate the expected anti-roll moment; and $i_{s,i}$ and $i_{\mathrm{act},i}$ are the installation ratios – approximated by constant values – of the suspension springs and actuators, where, based on (7), the latter generate a total force given by $F_{d,\mathrm{ij}}\pm F_{\mathrm{act},i}$. The passive spring force is modelled as: (8) $F_{s,\mathrm{ij}}=k_{s,i}{s}_{s,\mathrm{ij}}$(8) where $k_{s,i}$ is the linear stiffness of the $i$ passive spring, and $s_{s,\mathrm{ij}}$ is the respective displacement, defined as: $\begin{array}{rl}{s}_{s,\mathrm{FL}}& =-i_{s,F}[z_{\mathrm{hv}}-\theta a_F+0.5t_F\phi ]\end{array}$ $\begin{array}{rl}{s}_{s,\mathrm{FR}}& =-i_{s,F}[z_{\mathrm{hv}}-\theta a_F-0.5t_F\phi ]\end{array}$ $\begin{array}{rl}{s}_{s,\mathrm{RL}}& =-i_{s,R}[z_{\mathrm{hv}}+\theta a_R+0.5t_R\phi ]\end{array}$ (9) $\begin{array}{rl}{s}_{s,\mathrm{RR}}& =-i_{s,R}[z_{\mathrm{hv}}+\theta a_R-0.5t_R\phi ]\end{array}$(9) Within the prediction model, a continuous formulation, derived through curve fitting, approximates $F_{d,\mathrm{ij}}$ as a function of the relative actuator speed, ${\dot{s}}_{\mathrm{act},\mathrm{ij}}$, calculated with the time derivative of an expression equivalent to (9), and the current level of the relevant damping valve, ${\epsilon }_{\mathrm{ij}}$: (10) $\begin{array}{rl}{F}_{d,\mathrm{ij}}& =p_{00}[p_{01}+\text{atan}(p_{02}{\dot{s}}_{\mathrm{act},\mathrm{ij}}{\epsilon }_{\mathrm{ij}}+p_{03})]+p_{11}{\dot{s}}_{\mathrm{act},\mathrm{ij}}{\epsilon }_{\mathrm{ij}}+p_{22}{\dot{s}}_{\mathrm{act},\mathrm{ij}}^2{\epsilon }_{\mathrm{ij}}\\ & +p_{10}{\dot{s}}_{\mathrm{act},\mathrm{ij}}+p_{20}{\dot{s}}_{\mathrm{act},\mathrm{ij}}^2+p_{100}{\epsilon }_{\mathrm{ij}}\end{array}$(10) (10) provides a good approximation of the experimental behaviour of the available damping actuation hardware, see Figure 4.

Figure 4 Static front and rear damping force characteristics of the case study actuators, as functions of actuator speed, for three current levels (minimum, medium and maximum). Comparison between experimental measurement points, indicated by the dots, on an actuator test rig, and the NMPC formulation in (10).

The total active force contribution is designed to ensure a desired level of roll angle compensation. The total roll moment associated with a lateral acceleration $a_y$ is approximately given by $\frac{m_s{a}_y[h-h_{\mathrm{ra}}]}{t}$. The active suspension force contribution targets a level of roll moment compensation defined by the gain $k_{a_y}$, while distributing the active anti-roll moments between the front and rear axles according to the front-to-total distribution ratio $f_{\mathrm{nom}}+\mathrm{\Delta }f$, which is used to control the cornering response. Hence, $F_{\mathrm{act},F}$ and $F_{\mathrm{act},R}$ are expressed as: $\begin{array}{rl}{F}_{\mathrm{act},F}& =\frac{k_{a_y}{m}_s{a}_y^{\mathrm{ext}}[h-h_{\mathrm{ra}}]}{t_F}[f_{\mathrm{nom}}+\mathrm{\Delta }f]\end{array}$ (11) $\begin{array}{rl}{F}_{\mathrm{act},R}& =\frac{k_{a_y}{m}_s{a}_y^{\mathrm{ext}}[h-h_{\mathrm{ra}}]}{t_R}[1-f_{\mathrm{nom}}-\mathrm{\Delta }f]\end{array}$(11) where $a_y^{\mathrm{ext}}$ is the lateral acceleration measured by the inertial measurement unit (IMU) installed on the vehicle; and $f_{\mathrm{nom}}$ is the nominal front-to-total active anti-roll moment ratio. $f_{\mathrm{nom}}$ was set to 0.64, a value for which the cornering response is a desirable trade-off between steady-state and transient manoeuvring requirements, and is aligned to that of the passive version of the same vehicle with mechanical anti-roll bars.

Under reasonable simplifications, the resulting vertical forces applied by the equivalent rigid links to the vehicle body, $F_{\mathrm{pwt},z,\mathrm{ij}}-F_{b,z,\mathrm{ij}}$, are given by a moment balance about the centre of the tyre contact patch: (12) $\begin{array}{rl}{F}_{\mathrm{pwt},z,\mathrm{ij}}-F_{b,z,\mathrm{ij}}& \approx \frac{[F_{\mathrm{pwt},\mathrm{ij}}-F_{b,\mathrm{ij}}]d_i-F_{\mathrm{pwt},w,\mathrm{ij}}R+m_{\mathrm{us},i}{a}_x^{\mathrm{ext}}{h}_{\mathrm{us},i}}{e_i}\\ & \approx \frac{\left[\frac{T_{m,i}{\tau }_i}{2R}-\frac{T_{b,\mathrm{ij}}}{R}-m_{\mathrm{us},i}{a}_x^{\mathrm{ext}}\right]d_i-\frac{T_{m,i}{\tau }_i}{2}}{e_i}\end{array}$(12) where $F_{\mathrm{pwt},\mathrm{ij}}-F_{b,\mathrm{ij}}$ is obtained from the electric motor and friction brake torques, $T_{m,i}$ and $T_{b,\mathrm{ij}}$, and the inertial force of the unsprung mass; $a_x^{\mathrm{ext}}$ is the longitudinal acceleration measured by the IMU; $m_{\mathrm{us},i}$ is the unsprung mass of a corner of the $i$ axle; and ${\tau }_i$ is the transmission gear ratio. In (12) the term $F_{\mathrm{pwt},w,\mathrm{ij}}R$, where $F_{\mathrm{pwt},w,\mathrm{ij}}$ is the longitudinal traction force caused by the powertrain, represents the effect of the reaction torque of the on-board powertrain on the sprung mass, while the moment $m_{\mathrm{us},i}{a}_x^{\mathrm{ext}}{h}_{\mathrm{us},i}$ is neglected in the implementation. Similarly, the suspension jacking forces are approximated as: (13) $F_{y,z,\mathrm{ij}}\approx \frac{2[F_{y,\mathrm{ij}}^{\mathrm{ext}}-m_{\mathrm{us},i}{a}_y^{\mathrm{ext}}]h_{\mathrm{rc},i}}{t_i}$(13) where $F_{y,\mathrm{ij}}^{\mathrm{ext}}$ is the externally estimated lateral tyre force, considered constant along the prediction horizon, for decreasing the computational load.

The lateral tyre forces, $F_{y,\mathrm{ij}}$, which vary along the prediction horizon in (1) and (2), are computed through the version 5.2 of the magic formula as a function of the axle slip angle ${\alpha }_i$, vertical tyre load $F_{z,\mathrm{ij}}$, tyre slip ratio ${\sigma }_{\mathrm{ij},\mathrm{est}}$, camber angle ${\gamma }_{\mathrm{ij},\mathrm{est}}$ (its static value was used in the simulations), and tyre-road friction factor ${\mu }_{\mathrm{ij},\mathrm{est}}$, where the subscript ‘est’ indicates estimated variables that are considered to remain constant along the prediction horizon, for computational efficiency: (14) $F_{y,\mathrm{ij}}=F_{y,\mathrm{ij}}({\alpha }_{\mathrm{ij}},F_{z,\mathrm{ij}},{\sigma }_{\mathrm{ij},\mathrm{est}},{\gamma }_{\mathrm{ij},\mathrm{est}},{\mu }_{\mathrm{ij},\mathrm{est}})$(14) As the tyre model considers the effect of the interaction between the longitudinal and lateral slips [59], the formulation is suitable for modelling the vehicle response at and beyond the limit of handling. ${\alpha }_i$ is expressed through the following linearised formulations: $\begin{array}{rl}{\alpha }_F& \approx \beta +\frac{a_F\dot{\psi }}{V}-\delta \end{array}$ (15) $\begin{array}{rl}{\alpha }_R& \approx \beta -\frac{a_R\dot{\psi }}{V}\end{array}$(15) where $\delta$ is the average front steering angle. With the maximum slip angle magnitude seen in the broad range of test scenarios considered in Section 4 being only ∼10 deg, the approximation in (15) holds, which alleviates the computational load. The vertical tyre load, $F_{z,\mathrm{ij}}$, is calculated as the sum of the static load $F_{z,\mathrm{ij}}^0$, the spring and actuator force $F_{\mathrm{tot},z,\mathrm{ij}}$, the vertical force components through the rigid links, $F_{y,z,\mathrm{ij}}$ and $F_{\mathrm{pwt},z,\mathrm{ij}}-F_{b,z,\mathrm{ij}}$, and the load transfers, $F_{z,\mathrm{us},\mathrm{ij}}$, associated with the unsprung mass: $\begin{array}{rl}{F}_{z,\mathrm{FL}}& \approx F_{z,\mathrm{FL}}^0+F_{\mathrm{tot},z,\mathrm{FL}}-F_{y,z,\mathrm{FL}}+F_{b,z,\mathrm{FL}}-F_{\mathrm{pwt},z,\mathrm{FL}}+F_{z,\mathrm{us},\mathrm{FL}}\end{array}$ $\begin{array}{rl}{F}_{z,\mathrm{FR}}& \approx F_{z,\mathrm{FR}}^0+F_{\mathrm{tot},z,\mathrm{FR}}+F_{y,z,\mathrm{FR}}+F_{b,z,\mathrm{FR}}-F_{\mathrm{pwt},z,\mathrm{FR}}+F_{z,\mathrm{us},\mathrm{FR}}\end{array}$ $\begin{array}{rl}{F}_{z,\mathrm{RL}}& \approx F_{z,\mathrm{RL}}^0+F_{\mathrm{tot},z,\mathrm{RL}}-F_{y,z,\mathrm{RL}}-F_{b,z,\mathrm{RL}}+F_{\mathrm{pwt},z,\mathrm{RL}}+F_{z,\mathrm{us},\mathrm{RL}}\end{array}$ (16) $\begin{array}{rl}{F}_{z,\mathrm{RR}}& \approx F_{z,\mathrm{RR}}^0+F_{\mathrm{tot},z,\mathrm{RR}}+F_{y,z,\mathrm{RR}}-F_{b,z,\mathrm{RR}}+F_{\mathrm{pwt},z,\mathrm{RR}}+F_{z,\mathrm{us},\mathrm{RR}}\end{array}$(16) where the static loads are: $\begin{array}{rl}{F}_{z,\mathrm{FL}}^0& =F_{z,\mathrm{FR}}^0=\frac{\mathrm{mg}{a}_R}{2l}\end{array}$ (17) $\begin{array}{rl}{F}_{z,\mathrm{RL}}^0& =F_{z,\mathrm{RR}}^0=\frac{\mathrm{mg}{a}_F}{2l}\end{array}$(17) and the $F_{z,\mathrm{us},\mathrm{ij}}$ terms are given by: $\begin{array}{rl}{F}_{z,\mathrm{us},\mathrm{FL}}& \approx m_{\mathrm{us},F}{h}_{\mathrm{us},F}\left[-\frac{a_x^{\mathrm{ext}}}{l}-\frac{a_y^{\mathrm{ext}}}{t_F}\right]\end{array}$ $\begin{array}{rl}{F}_{z,\mathrm{us},\mathrm{FR}}& \approx m_{\mathrm{us},F}{h}_{\mathrm{us},F}\left[-\frac{a_x^{\mathrm{ext}}}{l}+\frac{a_y^{\mathrm{ext}}}{t_F}\right]\end{array}$ $\begin{array}{rl}{F}_{z,\mathrm{us},\mathrm{RL}}& \approx m_{\mathrm{us},R}{h}_{\mathrm{us},R}\left[\frac{a_x^{\mathrm{ext}}}{l}-\frac{a_y^{\mathrm{ext}}}{t_F}\right]\end{array}$ (18) $\begin{array}{rl}{F}_{z,\mathrm{us},\mathrm{RR}}& \approx m_{\mathrm{us},R}{h}_{\mathrm{us},R}\left[\frac{a_x^{\mathrm{ext}}}{l}+\frac{a_y^{\mathrm{ext}}}{t_R}\right]\end{array}$(18) where $h_{\mathrm{us},i}$ is the centre of gravity height of the respective unsprung mass.

The internal NMPC model is expressed through the following nonlinear continuous time formulation, based on the re-arrangement of (1)–(18): (19) $\dot{x}(t)=f(x(t),u(t))$(19) where the state vector, $x$, is: (20) $x=[\beta ,\dot{\psi },\phi ,\dot{\phi },\theta ,\dot{\theta },z_{\mathrm{hv}},{\dot{z}}_{\mathrm{hv}}{]}^T$(20) and the control input vector, $u$, which includes the control actions, is: (21) $u=[{\epsilon }_{\mathrm{NMPC},\mathrm{FL}},{\epsilon }_{\mathrm{NMPC},\mathrm{FR}},{\epsilon }_{\mathrm{NMPC},\mathrm{RL}},{\epsilon }_{\mathrm{NMPC},\mathrm{RR}},\mathrm{\Delta }f{]}^T$(21) The good match between the internal model running in open-loop and the experimental results in Figure 2 (Section 3.2) confirms the validity of the adopted assumptions.

In the following NMPC implementation, the current values of the states in (20) are used as initial conditions for the prediction. In particular, in a real vehicle application, $\dot{\psi }$, $\dot{\phi }$, and $\dot{\theta }$ are directly measured by the IMU; $\phi$, $\theta$, $z_{\mathrm{hv}}$ and ${\dot{z}}_{\mathrm{hv}}$ are inferred from the fusion of the displacement and speed measurements of the suspension actuators at the individual corners, according to the production algorithms already implemented by the involved industrial company; and $\beta$ is estimated through sensor fusion, e.g. see the broad set of options in [60–65].

4.2. Optimal control problem formulation

NMPC uses a model of the plant to predict and optimise the future behaviour of the system. The control action is obtained by solving, at each sampling time, a finite horizon optimal control problem, using the current value of the states. The optimisation yields an optimal control sequence, whose first element is applied to the plant, according to the receding horizon approach [66].

The proposed NMPC minimises the cost function $J_n$, subject to constraints, according to: $\underset{u}{min}{J}_n(x(0),u(\cdot )):=\sum _{k=0}^{N-1}\ell (x(k),u(k))+{\ell }_N(x(N))$ $\text{s}.\text{t}.$ $x(k+1)=f_d(x(k),u(k))$ $\underset{\_}{x}\le x(k)\le \overline{x}$ $\underset{\_}{u}\le u(k)\le \overline{u}$ (22) $u(\cdot ):[0,N-1]$(22) where $N$ is the number of steps of the prediction horizon $H_P$, in this implementation equal to the control horizon $H_c$, i.e. $H_c=H_p=N{T}_s$, with $T_s$ being the discretization time; $k$ is the discretization step; $\underset{\_}{x}$ and $\overline{x}$ are the lower and upper limits for $x$; $\underset{\_}{u}$ and $\overline{u}$ are the lower and upper limits for $u$; $x(k+1)=f_d(x(k),u(k))$ is the discretised version of the model in (19); ${\ell }_N(x(N))$ is the terminal cost; and $\ell (x(k),u(k))$is the stage cost for each time step, defined as a least-squares function: (23) $\begin{array}{l}\ell (x(k),u(k))=W_{e_{\dot{\psi }}}{\left[\frac{\dot{\psi }-{\dot{\psi }}_{\mathrm{ref}}}{\mathrm{\Delta }{\dot{\psi }}_{\max }}\right]}^2+W_{\dot{\phi }}{\left[\frac{\dot{\phi }}{{\dot{\phi }}_{\max }}\right]}^2+W_{\dot{\theta }}{\left[\frac{\dot{\theta }}{{\dot{\theta }}_{\max }}\right]}^2+W_{{\dot{z}}_{\mathrm{hv}}}{\left[\frac{{\dot{z}}_{\mathrm{hv}}}{{\dot{z}}_{\mathrm{hv},\max }}\right]}^2\\ +\frac{W_{\epsilon }}{\mathrm{\Delta }{{\epsilon }_{\max }}^2}\{{[{\epsilon }_{\mathrm{NMPC},\mathrm{FL}}-{\epsilon }_{\mathrm{SH},\mathrm{FL}}]}^2+{[{\epsilon }_{\mathrm{NMPC},\mathrm{FR}}-{\epsilon }_{\mathrm{SH},\mathrm{FR}}]}^2\\ +{[{\epsilon }_{\mathrm{NMPC},\mathrm{RL}}-{\epsilon }_{\mathrm{SH},\mathrm{RL}}]}^2+{[{\epsilon }_{\mathrm{NMPC},\mathrm{RR}}-{\epsilon }_{\mathrm{SH},\mathrm{RR}}]}^2\}+W_{\mathrm{\Delta }f}{\left[\frac{\mathrm{\Delta }f}{\mathrm{\Delta }{f}_{\max }}\right]}^2\end{array}$(23) where $\mathrm{\Delta }{\dot{\psi }}_{\max }$, ${\dot{\phi }}_{\max }$, ${\dot{\theta }}_{\max }$, ${\dot{z}}_{\mathrm{hv},\max }$, $\mathrm{\Delta }{\epsilon }_{\max }$, and $\mathrm{\Delta }{f}_{\max }$ are constant non-dimensionalisation factors, corresponding to the expected maximum value of the respective variable in the most critical conditions; and $W_{e_{\dot{\psi }}}$, $W_{\dot{\phi }}$, $W_{\dot{\theta }}$, $W_{{\dot{z}}_{\mathrm{hv}}}$, $W_{\epsilon }$, $W_{\mathrm{\Delta }f}$, and $W_{{\dot{z}}_{\mathrm{hv}}}$ are the cost function weights, respectively prioritising yaw rate error minimisation, roll rate reduction, pitch rate reduction, heave rate reduction, and control effort penalisation in terms of deviation of the damping valve currents and anti-roll moment distribution from their nominal values. The reference heave, pitch and roll rates for the sprung mass are set to zero. The reference currents, ${\epsilon }_{\mathrm{SH},\mathrm{ij}}$, are calculated through a pre-existing centralised skyhook formulation, see (29) in Section 4.5, which accounts for the front-to-total ratio, $\lambda$, computed from the damping force distribution output by the NMPC at the previous time step. More specifically, the skyhook forces from (29) are sent to inverse look-up tables corresponding to Figure 4, which generate ${\epsilon }_{\mathrm{SH},\mathrm{ij}}={\epsilon }_{\mathrm{ij}}(F_{d,\mathrm{ij}},{\dot{s}}_{\mathrm{ij}})$. The skyhook algorithm: a) enables the compensation of the effect of road irregularities, which is not a priority of the proposed NMPC formulation, and is already effectively carried out by the controller available at Tenneco Automotive. The NMPC-induced deviations from the reference skyhook currents address the body motions induced by the longitudinal and lateral accelerations, and the yaw rate and sideslip angle response; b) highlights the capability of the architecture to incorporate contributions from independent control functions, according to the current industrial practice in automotive suspension control; and c) steers the numerical optimisation of the implicit NMPC towards reasonable values, from which the controller can deviate.

The following constraints have been implemented in (22), to consider the actuation hardware capabilities:

• Limitation on the individual current levels, i.e. ${\epsilon }_{\min }\le {\epsilon }_{\mathrm{ij}}\le {\epsilon }_{\max }$, with ${\epsilon }_{\min }=$ 0.32 A and ${\epsilon }_{\max }=$ 1.6 A.

• Limitation on the variation of the front-to-total anti-roll moment distribution ratio w.r.t. the nominal value, i.e. $f_{\min }-f_{\mathrm{nom}}\le \mathrm{\Delta }f\le f_{\max }-f_{\mathrm{nom}}$, with $f_{\min }=$ 0.2 and $f_{\max }=$ 0.8.

• Limitation on the active force contributions, i.e. $F_{\mathrm{act},\min }\le F_{\mathrm{act},i}\le F_{\mathrm{act},\max }$, with $F_{\mathrm{act},\min }=-$ 7 kN and $F_{\mathrm{act},\max }=$ 7 kN.

4.3. Weight adaptation

Through (23), the NMPC incorporates different objectives in the same cost function, e.g. body control and cornering response enhancement. To maximise effectiveness, the controller can adapt its priorities to the driving scenarios by modifying its weights, e.g. to focus on body control during normal driving, and on yaw rate tracking in limit handling. Following extensive testing, the lateral acceleration magnitude, $|a_y|$, was selected as the weight scheduling variable, which is continuously modified. In fact, the effect of the anti-roll moment distribution on the cornering response is directly related to the lateral load transfers, which are proportional to $|a_y|$.

A simulation-based brute force algorithm was used to heuristically calibrate the NMPC weights for high tyre-road friction conditions, while considering a set of manoeuvres, each of them designed to target only low, medium or high $|a_y|$. A calibration cost function, $J_{\mathrm{tot},\xi }$, was minimised for each manoeuvre (see their list in the following Section 5), while producing a vehicle response that is free of any unwanted oscillations. $J_{\mathrm{tot},\xi }$, with $\xi =$ 1,..,3, is the weighted linear combination of the root mean square (‘$\mathrm{rms}$’ in the subscripts) values of the yaw rate error, and roll, pitch and heave rates: (24) $J_{\mathrm{tot},\xi }=W_{1,\xi }\frac{e_{\dot{\psi },\mathrm{rms}}}{e_{\dot{\psi },\mathrm{rms},\mathrm{passive}}}+W_{2,\xi }\frac{{\dot{\phi }}_{\mathrm{rms}}}{{\dot{\phi }}_{\mathrm{rms},\mathrm{passive}}}+W_{3,\xi }\frac{{\dot{\theta }}_{\mathrm{rms}}}{{\dot{\theta }}_{\mathrm{rms},\mathrm{passive}}}+W_{4,\xi }\frac{{\dot{z}}_{\mathrm{hv},\mathrm{rms}}}{{\dot{z}}_{\mathrm{hv},\mathrm{rms},\mathrm{passive}}}$(24) where: $e_{\dot{\psi },\mathrm{rms}}=\sqrt{\frac{1}{T_{\mathrm{fin}}-T_{\mathrm{in}}}\underset{T_{\mathrm{in}}}{\overset{T_{\mathrm{fin}}}{\int }}{[\dot{\psi }-{\dot{\psi }}_{\mathrm{ref}}]}^2\mathrm{dt}}$ ${\dot{\phi }}_{\mathrm{rms}}=\sqrt{\frac{1}{T_{\mathrm{fin}}-T_{\mathrm{in}}}\underset{T_{\mathrm{in}}}{\overset{T_{\mathrm{fin}}}{\int }}{\dot{\phi }}^2\mathrm{dt}}$ ${\dot{\theta }}_{\mathrm{rms}}=\sqrt{\frac{1}{T_{\mathrm{fin}}-T_{\mathrm{in}}}\underset{T_{\mathrm{in}}}{\overset{T_{\mathrm{fin}}}{\int }}{\dot{\theta }}^2\mathrm{dt}}$ (25) ${\dot{z}}_{\mathrm{hv},\mathrm{rms}}=\sqrt{\frac{1}{T_{\mathrm{fin}}-T_{\mathrm{in}}}\underset{T_{\mathrm{in}}}{\overset{T_{\mathrm{fin}}}{\int }}{{\dot{z}}_{\mathrm{hv}}}^2\mathrm{dt}}$(25) with $T_{\mathrm{in}}$ and $T_{\mathrm{fin}}$ being the initial and final times of the relevant part of the manoeuvre; and $e_{\dot{\psi },\mathrm{rms},\mathrm{passive}}$, ${\dot{\phi }}_{\mathrm{rms},\mathrm{passive}}$, ${\dot{\theta }}_{\mathrm{rms},\mathrm{passive}}$ and ${\dot{z}}_{\mathrm{hv},\mathrm{rms},\mathrm{passive}}$ being normalisation factors, defined as the values of the indicators for the passive vehicle, obtained along test scenario 5 in Section 5. In (24) the weights define the relative significance of the body control $({\dot{\phi }}_{\mathrm{rms}},{\dot{\theta }}_{\mathrm{rms}},{\dot{z}}_{\mathrm{hv},\mathrm{rms}})$ indicators and handling indicator $(e_{\dot{\psi },\mathrm{rms}})$. Based on the experience of the involved industrial partner, three sets of weights, defined by $\xi$, were used, depending on the nature of the considered test scenario, to: i) prioritise body control in non-critical vehicle dynamics conditions, which corresponds to the cost function $J_{\mathrm{tot},1}$; ii) achieve a balance between body control and vehicle dynamics in manoeuvres with significant lateral accelerations, although well within the cornering limit, which corresponds to $J_{\mathrm{tot},2}$; and iii) prioritise yaw rate tracking in limit handling operation, evaluated through $J_{\mathrm{tot},3}$.

Therefore, the three sets of weights were chosen through the aforementioned calibration method; more specifically, referring to the test scenarios defined in Section 5: i) test scenario 3 was used for selecting the weights below a lateral acceleration magnitude of 3.5 m/s², where such value approximately corresponds to the maximum $|a_y|$ in this manoeuvre, for which the objective was to minimise the cost function $J_{\mathrm{tot},1}$; ii) test scenario 4 was used for optimising the weights in the 3.5–6 m/s² lateral acceleration range, typical of this test, with the objective to minimise $J_{\mathrm{tot},2}$; and iii) test scenario 5 was used for the weights above 6 m/s², as $|a_y|$ exceeds 10 m/s² in this manoeuvre, with the objective of minimising $J_{\mathrm{tot},3}.$ In the online implementation, the weights are organised into maps, and vary linearly between the three sets of calibrated data to simplify the scheme and minimise the computational requirements, see Figure 5 ($W_{{\dot{z}}_{\mathrm{hv}}}$, $W_{\mathrm{\Delta }f}$ and $W_{\epsilon }$ are not reported as they are constant). Whilst the increase in $W_{e_{\dot{\psi }}}$ with $|a_y|$ is self-explanatory due to the prioritisation of yaw rate tracking during significant cornering, it was also found that for a given $W_{e_{\dot{\psi }}}$ value, an increase in $W_{\dot{\theta }}$ with $|a_y|$ tends to further reduce $J_{\mathrm{tot},\xi }$, because of the impact of the pitch rate on the longitudinal load transfer dynamics. The resulting maps were tested with excellent results across the test scenarios in Section 5. Future developments will evaluate more complex weight scheduling schemes.

Figure 5 Example of normalised NMPC weight adaptation map.

4.4. Real-time implementation of the NMPC algorithm

The proposed NMPC was implemented in implicit form via the ACADO toolkit [35], which offers a powerful interface for NMPC development, with the following settings: Gauss Newton Hessian approximation, multiple shooting discretisation, fourth order implicit Runge Kutta integrator, and qpOASES solver. The qpOASES solver, see its details in [67], is responsible for solving the optimal control problem of the NMPC as described in Section 4.2. It efficiently finds the optimal values of the control actions, as per (21), at each time step, to minimise the cost function in (22), while meeting the constraints outlined at the end of Section 4.2.

The controller sampling time, $T_s$, and number of optimisation steps, $N$, were set to 11 ms and 2, corresponding to a 22 ms prediction horizon. The number of iterations of the optimiser was set to 1, to ensure computational efficiency. Diminishing performance improvements were found with increasing $N$, and given that the execution time increases with increasing $N$, a low value of $N$ was used to ensure the real-time implementability of the controller, see Section 6 for the sensitivity analysis on $N$. The discretisation time of the internal model is 1 ms, which ensures its numerical stability without significantly affecting the computational time.

The indicated parametrisations enable the controller to run in real-time with good margin on the dSPACE MicroAutoBox II system (900 MHz, 16 Mb flash memory) in Figure 6. For example, during a sinusoidal steering manoeuvre at the limit of the handling, a maximum execution time of 6.8 ms was achieved.

Figure 6 Real-time implementation of the proposed NMPC on a dSPACE MicroAutoBox II unit.

4.5. Skyhook controllers

The benchmarking skyhook controllers use pseudoinverse and inverse formulations for computing the damping forces at the individual corners. The equivalent reference heave damping force at the centre of gravity of the sprung mass and the reference pitch and roll damping moments, respectively $F_{\mathrm{hv}}$, $M_{\theta }$ and $M_{\phi }$, are computed as functions of the estimated or measured ${\dot{z}}_{\mathrm{hv}}$, $\dot{\theta }$, and $\dot{\phi }$, through nonlinear look-up tables, see Figure 7, according to a skyhook approach based on the ride comfort and body motion control criteria deriving from the experience of Tenneco Automotive in production suspension controllers.

Figure 7 Nonlinear look-up tables for the computation of the reference heave force, and pitch and roll moments.

$F_{\mathrm{hv}}$, $M_{\theta }$ and $M_{\phi }$ can be expressed as functions of the damping forces at the corners: $F_{\mathrm{hv}}=-i_{d,F}[F_{d,\mathrm{FL}}+F_{d,\mathrm{FR}}]-i_{d,R}[F_{d,\mathrm{RL}}+F_{d,\mathrm{RR}}]$ $M_{\theta }=i_{d,F}{a}_F[F_{d,\mathrm{FL}}+F_{d,\mathrm{FR}}]-i_{d,R}{a}_R[F_{d,\mathrm{RL}}+F_{d,\mathrm{RR}}]$ (26) $M_{\phi }=0.5\{i_{d,F}{t}_F[F_{d,\mathrm{FR}}-F_{d,\mathrm{FL}}]+i_{d,R}{t}_R[F_{d,\mathrm{RR}}-F_{d,\mathrm{RL}}]\}$(26)

Hence, in the pseudoinverse case, $F_{d,\mathrm{ij}}$ is computed through a right pseudoinverse transformation of (26) as: (27) $\left[\begin{array}{c}{F}_{d,\mathrm{FL}}\\ F_{d,\mathrm{FR}}\\ F_{d,\mathrm{RL}}\\ F_{d,\mathrm{RR}}\end{array}\right]=\left[\begin{array}{ccc}{\displaystyle \frac{-0.5a_R}{i_{d,F}l}}& {\displaystyle \frac{0.5}{i_{d,F}l}}& {\displaystyle \frac{-i_{d,F}{t}_F}{i_{d,F}^2{t}_F^2+i_{d,R}^2{t}_R^2}}\\ {\displaystyle \frac{-0.5a_R}{i_{d,F}l}}& {\displaystyle \frac{0.5}{i_{d,F}l}}& {\displaystyle \frac{i_{d,F}{t}_F}{i_{d,F}^2{t}_F^2+i_{d,R}^2{t}_R^2}}\\ {\displaystyle \frac{-0.5a_F}{i_{d,F}l}}& {\displaystyle \frac{-0.5}{i_{d,R}l}}& {\displaystyle \frac{-i_{D,R}{t}_R}{i_{d,F}^2{t}_F^2+i_{d,R}^2{t}_R^2}}\\ {\displaystyle \frac{-0.5a_F}{i_{d,F}l}}& {\displaystyle \frac{-0.5}{i_{d,R}l}}& {\displaystyle \frac{i_{D,R}{t}_R}{i_{d,F}^2{t}_F^2+i_{d,R}^2{t}_R^2}}\end{array}\right]\left[\begin{array}{c}{F}_{\mathrm{hv}}\\ M_{\theta }\\ M_{\phi }\end{array}\right]$(27) with appropriate damping force saturations, functions of the individual actuator speeds, being imposed on the outputs of (27). In the pseudoinverse controller, which corresponds to the conventional suspension control implementation commercially proposed by the involved industrial partner, the active forces are obtained through (11) by imposing $\mathrm{\Delta }f=0$, i.e. the active anti-roll moment contribution is only used for roll angle compensation.

In the second benchmarking controller, i.e. the so-called inverse controller, an additional condition, on top of those in (26), is imposed, to account for the desirable front-to-total damping moment distribution $\lambda$: (28) $\lambda M_{\phi }=0.5i_{d,F}{t}_F[F_{d,\mathrm{FR}}-F_{d,\mathrm{FL}}]$(28) where $\lambda$ is used to enhance vehicle dynamics in transient cornering conditions. If the inverse controller is implemented in isolation, $\lambda$ is obtained as the sum of a nominal value and a contribution proportional to the yaw rate error. If the inverse controller operates in conjunction with the NMPC, $\lambda$ is equal to the resulting value from the NMPC at the previous time step. By adding (28) to (26), the system of equations includes four conditions to be achieved through four actuators, and therefore the damper forces are calculated through an inverse formulation: (29) $\left[\begin{array}{c}{F}_{d,\mathrm{FL}}\\ F_{d,\mathrm{FR}}\\ F_{d,\mathrm{RL}}\\ F_{d,\mathrm{RR}}\end{array}\right]=\left[\begin{array}{cccc}{\displaystyle \frac{-0.5a_R}{i_{d,F}l}}& {\displaystyle \frac{0.5}{i_{d,F}l}}& 0& -{\displaystyle \frac{1}{i_{d,F}{t}_F}}\\ {\displaystyle \frac{-0.5a_R}{i_{d,F}l}}& {\displaystyle \frac{0.5}{i_{d,F}l}}& 0& {\displaystyle \frac{1}{i_{d,F}{t}_F}}\\ {\displaystyle \frac{-0.5a_F}{i_{d,R}l}}& {\displaystyle \frac{-0.5}{i_{d,R}l}}& -{\displaystyle \frac{1}{i_{d,R}{t}_R}}& {\displaystyle \frac{1}{i_{d,R}{t}_R}}\\ {\displaystyle \frac{-0.5a_F}{i_{d,R}l}}& {\displaystyle \frac{-0.5}{i_{d,R}l}}& {\displaystyle \frac{1}{i_{d,R}{t}_R}}& -{\displaystyle \frac{1}{i_{d,R}{t}_R}}\end{array}\right]\left[\begin{array}{c}{F}_{\mathrm{hv}}\\ M_{\theta }\\ M_{\phi }\\ \lambda M_{\phi }\end{array}\right]$(29) The inverse formulation of the skyhook damping force contribution is coupled with an active contribution using an anti-roll moment distribution calculated through a PI controller of the yaw rate error, see the state-of-the-art design methodology in [4] and [34]. To prevent undesired oscillations during normal driving, the yaw rate based variations of $f$ and $\lambda$ are progressively activated when $|a_y|$ exceeds 0.4 g, and are fully activated for $|a_y|>$ 0.6 g [4].

5. Considered configurations, manoeuvres and performance indicators

The simulation analysis in Section 6, using the model for control system assessment in Section 3, compares the following vehicle and controller configurations:

• Passive: the case study vehicle without controllable dampers and active roll control, but including conventional passive anti-roll bars.

• Pseudoinverse: the controlled vehicle, i.e. without anti-roll bars, with the damper forces computed through (27), and using a constant active anti-roll moment distribution ratio, equal to $f_{\mathrm{nom}}$.

• Inverse: the controlled vehicle with the damper forces computed through (29), and active anti-roll moment distribution ratio from a PI controller.

• NMPC₁: the NMPC formulation in (1)-(23), with constant weights in the cost function in (22)-(23).

• NMPC₂: the proposed NMPC including the heuristic weight adaptation in Section 4.3.

To highlight the performance of the controlled suspension system, the simulations are carried out both with the de-activated and activated VSC. When it is switched on, the VSC – being a controller for emergency conditions – intervenes only in the two most extreme test scenarios (scenarios 5 and 6) of this section.

The following transient manoeuvres are simulated:

• Test scenario 1: a straight line manoeuvre at the constant speed $V=$ 50 km/h, on an experimentally available ride comfort road profile used by the involved company. This scenario is introduced to show that the additional functions of the NMPC architecture – which focuses on vehicle dynamics – do not compromise ride comfort w.r.t. the benchmarking skyhook controllers.

• Test scenario 2: a straight line hard braking test with an average deceleration of 9.3 m/s², on the verge of the anti-lock braking system (ABS) intervention, from a vehicle speed of 110 km/h, in high tyre-road friction conditions, i.e. with a friction factor $\mu \approx$ 1. The manoeuvre allows the assessment of the reduction – brought by the controlled suspension configurations – of the pitch and heave rates during the initial brake application and final brake release. As the test does not involve the ABS intervention, and the vehicle speed profile is an external input to be tracked by the driver model, the longitudinal vehicle dynamics are the same for all configurations, and the slip ratio behaviour can be excluded from the comparison.

• Test scenario 3: a sinusoidal steering test with a maximum steering wheel angle amplitude, $|{\delta }_{\mathrm{sw},\max }|$, of 30 deg, starting from a vehicle speed of 80 km/h, in high tyre-road friction conditions. The torque demand is constant and equal to the value required to keep the vehicle at the initial speed in straight line conditions, i.e. $V$ decreases during the test. This also applies to test scenarios 4, 5, and 6. The focus is the assessment of the reduction of the roll, pitch and heave rates through the semi-active system, and the roll angle peaks through the active anti-roll moment.

• Test scenario 4: a sinusoidal steering test with $|{\delta }_{\mathrm{sw},\max }|=$ 70 deg, starting from $V=$ 80 km/h, for $\mu =$ 1. The target is to achieve a good trade-off between handling and body control.

• Test scenario 5: the same as scenario 3, apart from $|{\delta }_{\mathrm{sw},\max }|=$ 150 deg, and the fact that the simulations are run for both $\mu =$ 1 and $\mu =$ 0.75, on a flat road profile, and on the same road profile as in test scenario 1. The main objective is vehicle stabilisation through yaw rate tracking.

• Test scenario 6: a sine-with-dwell test with $|{\delta }_{\mathrm{sw},\max }|=$ 150 deg, starting from a vehicle speed of 82 km/h in high friction conditions, with the same objective as scenario 5.

The comparison is based on the following indicators:

• $e_{\dot{\psi },\mathrm{rms}}$, which accounts for the vehicle cornering performance.

• $|{\alpha }_{R,\max }|$, i.e. the maximum rear axle sideslip angle magnitude, which is a cornering stability indicator.

• $\mathrm{\Delta }{F}_{z,\mathrm{rms}}^y$, given by: (30) $\mathrm{\Delta }{F}_{z,\mathrm{rms}}^y=\sqrt{\frac{1}{T_{\mathrm{fin}}-T_{\mathrm{in}}}\underset{T_{\mathrm{in}}}{\overset{T_{\mathrm{fin}}}{\int }}{[F_{z,\mathrm{FL}}+F_{z,\mathrm{RL}}-F_{z,\mathrm{FR}}-F_{z,\mathrm{RR}}]}^2\mathrm{dt}}$(30) which assesses the magnitude of the total lateral load transfer.

• ${\dot{\phi }}_{\mathrm{rms}}$, ${\dot{\theta }}_{\mathrm{rms}}$, ${\dot{z}}_{\mathrm{hv},\mathrm{rms}}$, as well the maximum roll angle magnitude, $|{\phi }_{\max }|$, assessing the body control quality.

• $J_{\mathrm{tot},\xi }$, which assesses the overall performance in cornering conditions. In particular, test scenario 3 is evaluated through $J_{\mathrm{tot},1}$; scenario 4 through $J_{\mathrm{tot},2}$; and scenarios 5 and 6 through $J_{\mathrm{tot},3}$.

• ${\ddot{z}}_{\mathrm{hv},\mathrm{rms}}$ and ${\ddot{\theta }}_{\mathrm{rms}}$, which assess ride comfort on irregular roads: ${\ddot{z}}_{\mathrm{hv},\mathrm{rms}}=\sqrt{\frac{1}{T_{\mathrm{fin}}-T_{\mathrm{in}}}\underset{T_{\mathrm{in}}}{\overset{T_{\mathrm{fin}}}{\int }}{\ddot{z}}_{\mathrm{hv}}^2\mathrm{dt}}$ (31) ${\ddot{\theta }}_{\mathrm{rms}}=\sqrt{\frac{1}{T_{\mathrm{fin}}-T_{\mathrm{in}}}\underset{T_{\mathrm{in}}}{\overset{T_{\mathrm{fin}}}{\int }}{\ddot{\theta }}^2\mathrm{dt}}$(31)

6. Simulation results

Table 2 in the Appendix reports the performance indicator values for test scenarios 3–6 (scenarios 1 and 2 are trivial given the absence of cornering). For all manoeuvres, $J_{\mathrm{tot},\xi }$ is the lowest for NMPC₂, which by itself provides a clear indication on the best configuration. A more detailed analysis is included in the following subsections.

Test scenario 1: straight ride comfort road at constant speed. In this test, the lateral acceleration is negligible, and therefore the performance is only determined by the semi-active contributions. Figure 8(a) shows a sample of the time history of ${\ddot{z}}_{\mathrm{hv}}$, and highlights the benefits of the pseudoinverse, inverse and NMPC₂ configurations w.r.t. the passive vehicle, particularly at the peak at ∼1 s, when the heave acceleration is reduced from almost 6 to 2.5 m/s². Figure 8(b) and (c) report ${\ddot{z}}_{\mathrm{hv},\mathrm{rms}}$ and ${\ddot{\theta }}_{\mathrm{rms}}$, and their percentage variations, $\mathrm{\Delta }{\ddot{z}}_{\mathrm{hv},\mathrm{rms}}$ and $\mathrm{\Delta }{\ddot{\theta }}_{\mathrm{rms}}$, w.r.t. the passive configuration (the symbol ‘$\mathrm{\Delta }$’ will be used with the same meaning in the remainder). Overall, NMPC₂ brings a similar percentage reduction to the benchmarking controllers, and outperforms all other cases in terms of ${\ddot{z}}_{\mathrm{hv},\mathrm{rms}}$, with a > 25% reduction w.r.t. the passive configuration.

Figure 8 Summary plots for test scenario 1. (a) time histories of the heave acceleration ${\ddot{z}}_{\mathrm{hv}}$; (b) values of ${\ddot{\theta }}_{\mathrm{rms}}$ and ${\ddot{z}}_{\mathrm{hv},\mathrm{rms}}$; and (c) values of $\mathrm{\Delta }{\ddot{\theta }}_{\mathrm{rms}}$ and $\mathrm{\Delta }{\ddot{z}}_{\mathrm{hv},\mathrm{rms}}$, i.e. the percentage variations of the performance indicators w.r.t. the passive case.

Test scenario 2: hard braking test. Given that also in this scenario only the damping force contributions are influential, the benefits are visible only during the transients occurring after 0.5 and 3.5 s, corresponding to the initial brake input application and final brake pressure release, see the profiles of $V$, $\dot{\theta }$, and ${\dot{z}}_{\mathrm{hv}}$ in Figure 9(a). While the benchmarking controllers only marginally reduce the peak values of $\dot{\theta }$ and ${\dot{z}}_{\mathrm{hv}}$, the NMPCs bring a significant decrease.

Figure 9 Summary plots for test scenario 2. (a) time histories of the vehicle speed $V$, pitch rate $\dot{\theta }$, and heave rate ${\dot{z}}_{\mathrm{hv}}$; (b) values of ${\dot{\theta }}_{\mathrm{rms}}$ and ${\dot{z}}_{\mathrm{hv},\mathrm{rms}}$; and (c) values of $\mathrm{\Delta }{\dot{\theta }}_{\mathrm{rms}}$ and $\mathrm{\Delta }{\dot{z}}_{\mathrm{hv},\mathrm{rms}}$, i.e. the percentage variations of the performance indicators w.r.t. the passive case.

Similar trends are visible in the histograms in Figure 9(b)-(c), reporting ${\dot{\theta }}_{\mathrm{rms}}$, which is the most relevant indicator in this scenario, and ${\dot{z}}_{\mathrm{hv},\mathrm{rms}}$, together with their percentage variations, $\mathrm{\Delta }{\dot{\theta }}_{\mathrm{rms}}$ and $\mathrm{\Delta }{\dot{z}}_{\mathrm{hv},\mathrm{rms}}$, amounting to ∼19% for ${\dot{\theta }}_{\mathrm{rms}}$ for the pseudoinverse and inverse cases, and exceeding 28% for the NMPCs.

Test scenarios 3 and 4: sinusoidal steering tests with low and medium steering amplitudes. Test scenario 3 implies maximum values of $|a_y|$ of ∼3 m/s², for which only body control – and especially roll control – is relevant. The active roll control contribution, which has the same total anti-roll moment control law in all configurations, nearly halves the roll angle peaks in Figure 10(a), according to the compensation level defined by $k_{a_y}$. Moreover, all controllers bring ${\dot{\phi }}_{\mathrm{rms}}$ reductions in excess of 45% (Figure 10(b)), and show substantially equivalent performance, see $J_{\mathrm{tot},1}$, with a very marginal penalty for NMPC₁, since the benchmarking skyhook controllers have their damping terms specifically optimised for body control, similarly to the adaptive cost function weights of NMPC₂, while NMPC₁ uses a trade-off calibration to meet the requirements for the whole operating range.

Figure 10 Summary plots for test scenario 3 ((a) time histories of roll angle $\phi$, roll rate $\dot{\phi }$, pitch angle $\theta$, and pitch rate $\dot{\theta }$; and (b) values of ${\dot{\phi }}_{\mathrm{rms}}$, ${\dot{\theta }}_{\mathrm{rms}}$, ${\dot{z}}_{\mathrm{hv},\mathrm{rms}}$, and $J_{\mathrm{tot},1}$); and for test scenario 4 ((c) values of $\mathrm{\Delta }|{\alpha }_{R,\max }|$, $\mathrm{\Delta }{\dot{\phi }}_{\mathrm{rms}}$, $\mathrm{\Delta }{\dot{\theta }}_{\mathrm{rms}}$, and $\mathrm{\Delta }{J}_{\mathrm{tot},2}$).

In test scenario 4, the car experiences a maximum value of $|a_y|$ of ∼6 m/s², which is beyond the range of typical driving conditions, but still well within the cornering limit. Therefore, the results, see Table 2 and Figure 10(c), are similar to those in scenario 3. Interestingly, the improved roll dynamics and consequent marginally reduced peaks of lateral load transfer bring an enhancement of the cornering response of all controlled configurations w.r.t. the passive vehicle. This is associated with a reduction of $|{\alpha }_{R,\max }|$, symptomatic of vehicle stability improvement, which is slightly more moderate (∼20%) for the pseudoinverse configuration, and the highest for the NMPCs. Like in all scenarios, NMPC₂ provides the highest reduction (∼43%) of the relevant cost function ($J_{\mathrm{tot},2}$), and in the specific manoeuvre is followed by the inverse controller.

Test scenarios 5 and 6: sinusoidal steering and sine-with-dwell tests at the limit of handling. Figure 11 reports the profiles of the main variables during test scenario 5, which is critical in terms of vehicle stability, since the passive configuration has $e_{\dot{\psi },\mathrm{rms}}$, $|{\alpha }_{R,\max }|$, and $|{\phi }_{\max }|$ values respectively exceeding 15 deg/s, 10 deg, and 6 deg. The pseudoinverse controller guarantees good body control performance, with a 3.6 deg peak value of roll angle, and ∼47% ${\dot{\phi }}_{\mathrm{rms}}$ reduction; however, the yaw rate errors and maximum rear sideslip angles, respectively ∼12 deg/s and ∼9 deg, are still well beyond acceptable thresholds. The configurations with controllable front-to-total anti-roll moment distribution bring a significant improvement of the cornering response, namely reduced $e_{\dot{\psi },\mathrm{rms}}$ and $|{\alpha }_{R,\max }|$ in comparison with the pseudoinverse and passive cases, and enable safe behaviour. In particular, these extreme conditions highlight the superior performance of the proposed NMPCs, which, for example, results in a substantial 13% decrease of $e_{\dot{\psi },\mathrm{rms}}$ and $|{\alpha }_{R,\max }|$, and a 3% reduction of ${\mathrm{\Delta F}}_{z,\mathrm{rms}}^y$ for NMPC₂ w.r.t. the state-of-the-art inverse controller. For completeness, Figure 12 shows the NMPC₂ weight profiles along the manoeuvre, which prioritise yaw rate tracking at high lateral accelerations, according to the adaption mechanism of Section 4.3, tailoring performance to the driving scenario. The trends remain the same when performing the manoeuvre on the irregular road profile, see Table 2.

Figure 11 Summary plots for test scenario 5 (VSC de-activated): (a) time histories of lateral acceleration $a_y$, yaw rate $\dot{\psi }$, rear axle sideslip angle ${\alpha }_R$, roll angle $\phi$, roll rate $\dot{\phi }$, pitch rate $\dot{\theta }$, corner damping forces $F_{d,\mathrm{ij}}$, and active axle forces $F_{\mathrm{act},i}$; (b) values of $e_{\dot{\psi },\mathrm{rms}}$, $|{\alpha }_{R,\max }|$, and ${\dot{\phi }}_{\mathrm{rms}}$; and (c) values of $\mathrm{\Delta }{e}_{\dot{\psi },\mathrm{rms}}$, $\mathrm{\Delta }|{\alpha }_{R,\max }|$, $\mathrm{\Delta }{\dot{\phi }}_{\mathrm{rms}}$, and $\mathrm{\Delta }{J}_{\mathrm{tot},3}$.

Figure 12 Weight adaptation profiles for NMPC₂ during test scenario 5.

As a robustness check, the test scenario 5 simulations were repeated by imposing $\mu =$ 0.75 in the model for control system assessment, but without modifying any controller parameters, i.e. by leaving the tyre-road friction factor set to 1 in the internal models of the NMPCs, and for computing the reference yaw rate. Despite this, the NMPC performance remains safe, and clearly exceeds the one of the passive configuration as well as those of the benchmarking skyhook controllers, e.g. see the reduction of the peaks of yaw rate and rear sideslip angle in Figure 13(a), and the performance indicators in Figure 13(b)-(c). When also the VSC system is active, see Figure 13(d)–(f), the direct yaw moment actuation brings a major reduction of the magnitude of the oscillations and overshoots, especially in terms of sideslip angle, which enhances vehicle stability. Nevertheless, the ranking of the configurations w.r.t. the suspension controllers remains the same as for the cases without VSC, i.e. NMPC₂ brings the lowest $J_{\mathrm{tot},3}$ value.

Figure 13 Summary plots for test scenario 5 with $\mu =$ 0.75, and incorrect tyre-road friction information sent to the NMPCs: in (a), (b) and (c) only the suspension controllers are active, while in (d), (e) and (f) both the suspension controllers and the VSC are active.

Interestingly, based on Table 2, in high tyre-road friction conditions (see also the experimental results in [4]), the controlled suspension system on its own is more effective that the VSC in isolation in controlling the yaw rate and sideslip response. This is caused by: i) the high nonlinearity of the specific tyres at high lateral load transfers; and ii) the continuous operation of the proposed active suspension system, whereas the VSC intervenes only in case of significant yaw rate errors, and – in the specific tuning of this paper – not very intrusively, since the brake interventions are clearly perceived by the driver and passengers. In case of reduced $\mu$, because of the decreased lateral load transfer, the suspension on its own is not as effective as in high tyre-road friction conditions, nevertheless it supports the VSC.

For a fair assessment of the results, it must be considered that the active suspension effect on the cornering response is very dependent on the nonlinearity level of the lateral tyre force with respect to the vertical tyre load, i.e. the benefit is caused by the fact that tyre performance decreases more on the inner corner than it increases on the outer corner of the same axle, because of the lateral load transfers. Therefore, the design and tuning of the active suspension contribution to enhance the cornering response is very subtle (see [4]) and reliant of the properties of the individual tyre makes, which can change during the vehicle life span, e.g. when the tyres are replaced with different ones. On the contrary, the reference direct yaw moment output by a VSC can be predictably generated, rather independently from the installed tyres, up to the tyre saturation level.

Given the intrinsic reliance of the NMPC implementations on their internal model, controller robustness was also evaluated by varying, within the model for control system assessment: i) the tyre cornering stiffness $C_y$ (by up to ±20%), achieved through the calibration of the respective magic formula scaling factor; and ii) the vehicle mass and yaw mass moment of inertia, which were increased by up to 20%. In the meantime, all controller parameters – including those subject to variations in the sensitivity analysis – were kept unchanged. The results are in Figure 14(a). For both i) and ii) the cornering dynamics variation is more significant for the pseudoinverse and inverse configurations than for the NMPCs, which provide consistently better results than the benchmarking implementations. Moreover, for NMPC₂, Figure 14(b) reports the sensitivity of $\mathrm{\Delta }{e}_{\dot{\psi },\mathrm{rms}}$ and $\mathrm{\Delta }|{\alpha }_{R,\max }|$ to the number of prediction steps $N$. A moderate increase of $N$ brings a marginal performance improvement, at the price of increased computational effort, which compromises the real-time implementability on the available control hardware, see Section 4.4.

Figure 14 Sensitivity analysis of the performance, expressed by $\mathrm{\Delta }{e}_{\dot{\psi },\mathrm{rms}}$ and $\mathrm{\Delta }|{\alpha }_{R,\max }|$ along test scenario 5, of: (a) the proposed controllers to vehicle parameter variations; and (b) NMPC₂ to the number of prediction steps N.

The trends are confirmed by test scenario 6, where the $e_{\dot{\psi },\mathrm{rms}}$, $|{\alpha }_{R,\max }|$ and $J_{\mathrm{tot},3}$ reductions w.r.t. the passive vehicle range from 20% for the pseudoinverse configuration to >60% for NMPC₂.

7. Conclusions

This study presented two novel real-time-capable NMPC implementations for combined semi-active and active suspension control, tailored to a new hydraulic actuation system developed by Tenneco Automotive. The algorithms target the compensation of the body motions induced by the driving actions, and cornering response enhancement. They are based on nonlinear 5-degree-of-freedom vehicle models, including consideration of tyre slip behaviour as well as suspension anti-properties and jacking force effects. The NMPCs can be integrated with other algorithms, based on different control technologies, targeting complementary aspects, e.g. ride comfort enhancement on irregular roads. The novel controllers were compared with: a) a benchmarking skyhook formulation using a pseudoinverse transformation, combined with an active roll control contribution with fixed front-to-rear moment distribution, corresponding to the current level of implementation on production vehicles; and b) a second skyhook controller, i.e. the so-called inverse formulation, coupled with a front-to-total anti-roll moment distribution algorithm for yaw rate control, the latter designed through a recently published state-of-the-art linear model-based methodology.

The simulation analysis, using an experimentally validated model of a case study electric vehicle, brings the following conclusions:

• The NMPCs enhance the level of body control during acceleration and braking transients, covered by test scenario 2, with >28% reduction of the root mean square value of the pitch rate (${\dot{\theta }}_{\mathrm{rms}}$) w.r.t. the passive configuration, against the ∼19% decrease brought by the benchmarking controllers.

• In transient cornering at low-to-medium lateral acceleration levels, corresponding to test scenarios 3 and 4, the NMPCs and skyhook controllers tend to provide rather similar body control performance, with a marginal benefit for the NMPC formulation with weight adaptation (NMPC₂), and a reduction (e.g. ≥ 20% in scenario 4) of the peak values of the rear axle sideslip angle ($|{\alpha }_{R,\max }|$).

• During limit handling operation, see test scenarios 5 and 6, all controllers with front-to-total anti-roll moment distribution capability bring a major enhancement in the cornering response and active safety, e.g. with >65% and >50% reductions of the root mean square value of the yaw rate error, $e_{\dot{\psi },\mathrm{rms}}$, and $|{\alpha }_{R,\max }|$. Importantly, the NMPCs enable better performance than the inverse skyhook controller both in terms of vehicle dynamics and body control, e.g. see the 13% decrease of $e_{\dot{\psi },\mathrm{rms}}$ and $|{\alpha }_{R,\max }|$, and the 14% and 8% decrements of the peak values of roll angle magnitude ($|{\phi }_{\max }|$) and root mean square roll rate (${\dot{\phi }}_{\mathrm{rms}}$) for NMPC₂ and NMPC₁ in test scenario 5. The trends are confirmed also in case of VSC activation.

• Despite embedding a detailed vehicle model in their online algorithms, both NMPCs provide higher robustness than the benchmarking controllers w.r.t. significant parameter variations, namely of vehicle mass and yaw mass moment of inertia, tyre cornering stiffness, as well as estimated tyre-road friction factor.

• The weight adaptation mechanism of NMPC₂ allows this set-up to be the best among the considered suspension configurations in all scenarios, see the values of the cost functions $J_{\mathrm{tot},\xi }$, without the compromises, typical of the calibration of the other algorithms, related to the different operating conditions and respective control objectives. Also, on irregular road profiles, see test scenario 1, NMPC₂ provides comparable ride comfort results to those of the benchmarking skyhook controllers.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Appendix

Table 2. Performance indicator values for test scenarios 3–6 (different indicators are used for test scenarios 1 and 2). The bold fonts highlight the controller configuration with the lowest cost function value in the respective scenario.

	$e_{\dot{\psi },\mathrm{rms}}$ [deg/s]	$|{\alpha }_{R,\max }|$ [deg]	$\mathrm{\Delta }{F}_{z,\mathrm{rms}}^y$ [kN]	$|{\phi }_{\max }|$ [deg]	${\dot{\phi }}_{\mathrm{rms}}$ [deg/s]	${\dot{\theta }}_{\mathrm{rms}}$ [deg/s]	${\dot{z}}_{\mathrm{hv},\mathrm{rms}}$ [m/s]	$J_{\mathrm{tot},\xi }$ [-]
Test scenario 3
Passive	0.29	0.82	3.74	1.58	3.69	0.10	0.0025	0.128
Pseudoinverse	0.43	0.70	3.50	0.88	2.00	0.04	0.0012	0.071
Inverse	0.43	0.70	3.50	0.88	2.00	0.05	0.0012	0.072
NMPC1	0.43	0.71	3.50	0.88	1.98	0.05	0.0012	0.072
NMPC2	0.43	0.70	3.50	0.86	1.96	0.04	0.0012	0.070
Test scenario 4
Passive	0.68	1.77	6.93	2.87	7.06	0.13	0.0066	0.200
Pseudoinverse	0.84	1.41	6.50	1.64	3.74	0.05	0.0031	0.117
Inverse	0.78	1.39	6.44	1.59	3.60	0.09	0.0031	0.114
NMPC1	0.78	1.37	6.41	1.62	3.69	0.12	0.0032	0.119
NMPC2	0.74	1.37	6.36	1.58	3.57	0.10	0.0028	0.113
Test scenario 5
Passive	16.07	10.05	14.38	6.29	16.13	3.13	0.0281	1.000
Passive + VSC	6.12	5.87	12.74	5.09	13.66	2.06	0.0173	0.488
Pseudoinverse	11.97	8.18	14.27	3.60	8.53	0.50	0.0224	0.672
Pseudoinv.+VSC	5.50	5.13	12.70	3.35	7.96	1.48	0.0169	0.379
Inverse	5.10	4.91	13.03	4.08	9.36	1.25	0.0292	0.374
Inverse + VSC	3.92	4.00	12.54	3.13	8.10	1.25	0.0202	0.304
NMPC1	4.96	4.87	13.01	3.50	8.62	0.63	0.0272	0.349
NMPC1 + VSC	3.78	4.02	12.53	3.03	7.99	0.93	0.0200	0.290
NMPC2	4.49	4.63	12.91	3.89	9.20	0.69	0.0251	0.335
NMPC2 + VSC	3.68	3.99	12.59	3.20	8.42	0.84	0.0198	0.289
Test scenario 5 – on irregular road profile
Passive	16.04	9.93	14.52	6.30	16.30	2.52	0.0283	0.991
Pseudoinverse	12.44	8.58	14.35	3.56	8.61	1.08	0.0251	0.705
Inverse	5.52	4.98	13.03	4.00	9.20	1.34	0.0300	0.393
NMPC1	5.37	4.95	13.00	3.54	8.66	1.22	0.0293	0.378
NMPC2	4.65	4.76	13.00	4.09	9.40	1.16	0.0298	0.353
Test scenario 5 – $\mu =$ 0.75 only in the plant model
Passive	16.72	9.46	12.20	4.85	11.85	1.32	0.0082	0.948
Passive + VSC	6.30	5.10	10.27	3.90	10.15	1.86	0.0064	0.450
Pseudoinverse	15.66	8.85	12.05	2.75	6.15	0.38	0.0111	0.813
Pseudoinv.+VSC	5.95	4.52	10.21	2.36	5.74	1.43	0.0118	0.372
Inverse	10.80	7.17	11.64	2.73	6.50	0.50	0.0152	0.593
Inverse + VSC	5.03	3.94	10.25	2.35	5.85	1.58	0.0092	0.333
NMPC1	9.99	6.82	11.56	2.76	6.53	0.41	0.0182	0.554
NMPC1 + VSC	4.84	3.88	10.22	2.46	6.10	1.28	0.0101	0.322
NMPC2	9.64	6.54	11.53	2.98	6.93	0.48	0.0196	0.544
NMPC2 + VSC	4.80	3.88	10.21	2.38	5.91	1.32	0.0091	0.318
Test scenario 6
Passive	7.72	10.38	14.33	5.70	10.42	1.21	0.0100	0.509
Passive + VSC	4.56	6.34	12.89	4.99	9.32	1.37	0.0131	0.351
Pseudoinverse	6.10	7.68	14.00	3.98	5.82	0.75	0.0308	0.369
Pseudoinv.+VSC	4.06	5.48	12.88	3.08	5.18	1.22	0.0189	0.273
Inverse	3.03	4.26	12.50	3.47	5.39	0.57	0.0231	0.217
Inverse + VSC	2.64	3.85	12.19	3.02	5.25	0.88	0.0190	0.201
NMPC1	2.85	4.07	12.19	3.26	5.32	0.84	0.0238	0.212
NMPC1 + VSC	2.47	3.72	11.99	3.26	5.28	0.88	0.0199	0.195
NMPC2	2.59	3.91	12.22	3.32	5.39	0.66	0.0233	0.198
NMPC2 + VSC	2.29	3.64	12.03	3.0	5.20	0.70	0.0160	0.183