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

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.


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 modelbased 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.

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][6][7][8][9][10][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][14][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][19][20][21][22][23][24][25][26][27][28][29][30][31][32][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][19][20][21][22][23][24][25][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][28][29][30][31][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][39][40][41][42][43][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 MPCbased 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 2 , 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.

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.

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 acc , brake pedal position, p bk , and steering wheel angle, δ sw .For tracking the reference speed profiles, proportional integral (PI) controllers are used to output the desired pedal positions, while the δ sw profiles are directly generated with fixed look-up tables of δ 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,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 ψref , see [4] for its definition, φref , θref , and żhv,ref ), as well the reference anti-roll moment, anti-pitch moment, and anti-heave force (M ϕ , M θ , and F hv ), starting from the driver steering input, δ 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.) ψ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,VSC,i and T b,VSC,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 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 nom , as well as the current levels, ε NMPC,ij , for the valves of the controllable dampers.The ε NMPC,ij values are computed from the reference currents, ε SH,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 ε NMPC,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 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.

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.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 • Yaw moment balance at the vehicle level • Roll moment balance of the sprung mass • Pitch moment balance of the sprung mass • Heave force balance of the sprung mass where β, ψ, φ, θ, and zhv are the sideslip rate, yaw acceleration, roll acceleration, pitch acceleration, and heave acceleration; F y,ij is the lateral tyre force at the ij corner; F tot,z,ij is the total force through the deformable suspension components (i.e.springs and actuators) of the ij corner; F b,ij − F pwt,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,ij − F pwt,z,ij ; F y,z,ij is the corner's jacking force; M z,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,VSC,ij ); and h 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: where l = a F + a R is the wheelbase; and h rc,i is the roll centre height of the i axle.
The terms F tot,z,ij are given by: where F s,ij is the force of the ij passive spring; F d,ij is the virtual damping force of the ij suspension actuator; F 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 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,ij ± F act,i .The passive spring force is modelled as: where k s,i is the linear stiffness of the i passive spring, and s s,ij is the respective displacement, defined as: Within the prediction model, a continuous formulation, derived through curve fitting, approximates F d,ij as a function of the relative actuator speed, ṡact,ij , calculated with the time derivative of an expression equivalent to (9), and the current level of the relevant damping valve, ε ij : (10) provides a good approximation of the experimental behaviour of the available damping actuation hardware, see Figure 4.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 . 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 nom + f , which is used to control the cornering response.Hence, F act,F and F act,R are expressed as: where a ext y is the lateral acceleration measured by the inertial measurement unit (IMU) installed on the vehicle; and f nom is the nominal front-to-total active anti-roll moment ratio.f nom was set to 0.64, a value for which the cornering response is a desirable tradeoff 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 pwt,z,ij − F b,z,ij , are given by a moment balance about the centre of the tyre contact patch: where F pwt,ij − F b,ij is obtained from the electric motor and friction brake torques, T m,i and T b,ij , and the inertial force of the unsprung mass; a ext x is the longitudinal acceleration measured by the IMU; m us,i is the unsprung mass of a corner of the i axle; and τ i is the transmission gear ratio.In (12) the term F pwt,w,ij R, where F pwt,w,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 us,i a ext x h us,i is neglected in the implementation.Similarly, the suspension jacking forces are approximated as: where F ext y,ij is the externally estimated lateral tyre force, considered constant along the prediction horizon, for decreasing the computational load.
The lateral tyre forces, F y,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 α i , vertical tyre load F z,ij , tyre slip ratio σ ij,est , camber angle γ ij,est (its static value was used in the simulations), and tyre-road friction factor μ ij,est , where the subscript 'est' indicates estimated variables that are considered to remain constant along the prediction horizon, for computational efficiency: 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.α i is expressed through the following linearised formulations: where δ 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,ij , is calculated as the sum of the static load F 0 z,ij , the spring and actuator force F tot,z,ij , the vertical force components through the rigid links, F y,z,ij and F pwt,z,ij − F b,z,ij , and the load transfers, F z,us,ij , associated with the unsprung mass: where the static loads are: and the F z,us,ij terms are given by: where h 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): where the state vector, x, is: and the control input vector, u, which includes the control actions, is: 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, ψ, φ, and θ are directly measured by the IMU; ϕ, θ, z hv and ż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 β is estimated through sensor fusion, e.g.see the broad set of options in [60][61][62][63][64][65].

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: 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 = NT s , with T s being the discretization time; k is the discretization step; x and x are the lower and upper limits for x; u and ū are the lower and upper limits for u; x(k is the discretised version of the model in (19); N (x(N)) is the terminal cost; and (x(k), u(k))is the stage cost for each time step, defined as a least-squares function: where ψmax , φmax , θmax , żhv,max , ε max , and 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 ψ , W φ , W θ , W żhv , W ε , W f , and W ż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, ε SH,ij , are calculated through a pre-existing centralised skyhook formulation, see (29) in Section 4.5, which accounts for the front-to-total ratio, λ, 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 ε SH,ij = ε ij (F d,ij , ṡ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. ε min ≤ ε ij ≤ ε max , with ε min = 0.32 A and ε 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 nom ≤ f ≤ f max − f nom , with f min = 0.2 and f max = 0.8.• Limitation on the active force contributions, i.e.F act,min ≤ F act,i ≤ F act,max , with F act,min = − 7 kN and F act,max = 7 kN.

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 antiroll 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 tot,ξ , 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 tot,ξ , with ξ = 1,..,3, is the weighted linear combination of the root mean square ('rms' in the subscripts) values of the yaw rate error, and roll, pitch and heave rates: where: hv dt (25) with T in and T fin being the initial and final times of the relevant part of the manoeuvre; and e ψ,rms,passive , φrms,passive , θrms,passive and żhv,rms,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 ( φrms , θrms , żhv,rms ) indicators and handling indicator (e ψ,rms ).Based on the experience of the involved industrial partner, three sets of weights, defined by ξ , 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 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 tot,2 ; and iii) prioritise yaw rate tracking in limit handling operation, evaluated through J 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 2 , where such value approximately corresponds to the maximum |a y | in this manoeuvre, for which the objective was to minimise the cost function J tot,1 ; ii) test scenario 4 was used for optimising the weights in the 3.5-6 m/s 2 lateral acceleration range, typical of this test, with the objective to minimise J tot,2 ; and iii) test scenario 5 was used for the weights above 6 m/s 2 , as |a y | exceeds 10 m/s 2 in this manoeuvre, with the objective of minimising J 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 żhv , W f and W ε are not reported as they are constant).Whilst the increase in W e ψ 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 ψ value, an increase in W θ with |a y | tends to further reduce J tot,ξ , 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.

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.

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 hv , M θ and M ϕ , are computed as functions of the estimated or measured żhv , θ, and φ, 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.F hv , M θ and M ϕ can be expressed as functions of the damping forces at the corners: Hence, in the pseudoinverse case, F d,ij is computed through a right pseudoinverse transformation of ( 26) as: 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 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-tototal damping moment distribution λ: where λ is used to enhance vehicle dynamics in transient cornering conditions.If the inverse controller is implemented in isolation, λ 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, λ 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: 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 λ are progressively activated when |a y | exceeds 0.4 g, and are fully activated for |a y | > 0.6 g [4].

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 nom .• Inverse: the controlled vehicle with the damper forces computed through (29), and active anti-roll moment distribution ratio from a PI controller.• NMPC 1 : the NMPC formulation in (1)-( 23), with constant weights in the cost function in ( 22)-( 23).• NMPC 2 : 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 VSCbeing 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 2 , 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 μ ≈ 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, |δ 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.The comparison is based on the following indicators: • e ψ,rms , which accounts for the vehicle cornering performance.
• |α R,max |, i.e. the maximum rear axle sideslip angle magnitude, which is a cornering stability indicator.• F y z,rms , given by: which assesses the magnitude of the total lateral load transfer.• φrms , θrms , żhv,rms , as well the maximum roll angle magnitude, |ϕ max |, assessing the body control quality.• J tot,ξ , which assesses the overall performance in cornering conditions.In particular, test scenario 3 is evaluated through J tot,1 ; scenario 4 through J tot,2 ; and scenarios 5 and 6 through J tot,3 .
• zhv,rms and θrms , which assess ride comfort on irregular roads:

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 tot,ξ is the lowest for NMPC 2 , 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 zhv , and highlights the benefits of the pseudoinverse, inverse and NMPC 2 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 2 .Figure 8(b) and (c) report zhv,rms and θrms , and their percentage variations, zhv,rms and θrms , w.r.t. the passive configuration (the symbol ' ' will be used with the same meaning in the remainder).Overall, NMPC 2 brings a similar percentage reduction to the benchmarking controllers, and outperforms all other cases in terms of zhv,rms , with a > 25% reduction w.r.t. the passive configuration.
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, θ, and żhv in Figure 9(a).While the benchmarking controllers only marginally reduce the peak values of θ and żhv , the NMPCs bring a significant decrease.
Similar trends are visible in the histograms in Figure 9(b)-(c), reporting θrms , which is the most relevant indicator in this scenario, and żhv,rms , together with their percentage variations, θrms and żhv,rms , amounting to ∼ 19% for θ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 2 , 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 φrms reductions in excess of 45% (Figure 10(b)), and show substantially equivalent performance, see J tot,1 , with a very marginal penalty for NMPC 1 , since the benchmarking skyhook controllers have their damping terms specifically optimised for body control, similarly to the adaptive cost function weights of NMPC 2 , while NMPC 1 uses a trade-off calibration to meet the requirements for the whole operating range.
In test scenario 4, the car experiences a maximum value of |a y | of ∼ 6 m/s 2 , 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 |α 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 2 provides the highest reduction ( ∼ 43%) of the relevant cost function (J 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.  of roll angle, and ∼ 47% φ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 ψ,rms and |α 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 ψ,rms and |α R,max |, and a 3% reduction of ΔF y z,rms for NMPC 2 w.r.t. the state-of-the-art inverse controller.For completeness, Figure 12 shows the NMPC 2 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.As a robustness check, the test scenario 5 simulations were repeated by imposing μ = 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 2 brings the lowest J tot,3 value.
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 μ, 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 analysiswere 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 2 , Figure 14

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.

Figure 1 .
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 .
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.

Figure 3 .
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.

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).

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

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

•
Test scenario 4: a sinusoidal steering test with |δ sw,max | = 70 deg, starting from V = 80 km/h, for μ = 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 |δ sw,max | = 150 deg, and the fact that the simulations are run for both μ = 1 and μ = 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 |δ sw,max | = 150 deg, starting from a vehicle speed of 82 km/h in high friction conditions, with the same objective as scenario 5.

Figure 8 .
Figure 8. Summary plots for test scenario 1.(a) time histories of the heave acceleration zhv ; (b) values of θrms and zhv,rms ; and (c) values of θrms and zhv,rms , i.e. the percentage variations of the performance indicators w.r.t. the passive case.
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 ψ,rms , |α R,max |, and |ϕ 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

Figure 9 .
Figure 9. Summary plots for test scenario 2. (a) time histories of the vehicle speed V, pitch rate θ, and heave rate żhv ; (b) values of θrms and żhv,rms ; and (c) values of θrms and żhv,rms , i.e. the percentage variations of the performance indicators w.r.t. the passive case.
(b) reports the sensitivity of e ψ,rms and |α 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.The trends are confirmed by test scenario 6, where the e ψ,rms , |α R,max | and J tot,3 reductions w.r.t. the passive vehicle range from 20% for the pseudoinverse configuration to > 60% for NMPC 2 .

Figure 13 .
Figure 13.Summary plots for test scenario 5 with μ = 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.

Figure 14 .
Figure 14.Sensitivity analysis of the performance, expressed by e ψ,rms and |α R,max | along test scenario 5, of: (a) the proposed controllers to vehicle parameter variations; and (b) NMPC 2 to the number of prediction steps N.

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.