On Pre-Emptive In-Wheel Motor Control for Reducing the Longitudinal Acceleration Oscillations Caused by Road Irregularities

Road irregularities induce vertical and longitudinal vibrations of the sprung and unsprung masses, which affect vehicle comfort. While the vertical dynamics and related compensation techniques are extensively covered by the suspension control literature, the longitudinal dynamics on uneven road surfaces are less frequently addressed, and are significantly influenced by the tires and suspension systems. The relatively slow response of internal combustion engines does not allow any form of active compensation of the effect of road irregularities. However, in-wheel electric powertrains, in conjunction with pre-emptive control based on the information on the road profile ahead, have some potential for effective compensation, which, however, has not been explored yet. This paper presents a proof-of-concept nonlinear model predictive control (NMPC) implementation based on road preview, which is preliminarily assessed with a simulation model of an all-wheel drive electric vehicle with in-wheel motors, including a realistic tire model for ride comfort simulation. The major improvement brought by the proposed road preview controller is evaluated through objective performance indicators along multiple maneuvers, and is confirmed by the comparison with two benchmarking feedback controllers from the literature.


I. INTRODUCTION
LECTRIC vehicles (EVs) are becoming an important and increasing share of the automotive market. Electric powertrains usually provide enhanced wheel torque control performance, because of the accuracy and fast dynamics of electric motors, in comparison with conventional internal combustion engines and friction brakes. These benefits are particularly evident for in-wheel powertrain configurations [1], i.e., in which the electric motor is part of the unsprung mass, and is connected to the wheel hub either directly, which is the case of the direct drive in-wheel powertrains, or through a compact in-wheel single-speed transmission, which is the case of the so-called near-the-wheel configurations. In in-wheel powertrains, which are the object of this study, the absence of the torsional dynamics of the half-shafts: a) facilitates precise wheel torque modulation, and thus improves the performance of wheel slip and direct yaw moment controllers; and b) permits to omit the conventional anti-jerk control function [2]- [4], which, in on-board powertrains, attenuates the torsional halfshaft dynamics induced by the powertrain torque transients. The main source of torsional in-wheel drivetrain dynamics is the tire, and although in practical EV implementations no specific controller is strictly required to attenuate the associated vibrations, the study in [5] presents a dedicated fuzzy logic algorithm. The torque ripple of in-wheel machines, which is directly transmitted to the wheel, can also affect comfort on smooth road surfaces [6]- [7], and is attenuated through appropriate motor design or control [8]- [10]. Moreover, given the absence or the significantly reduced volume of the powertrain components installed on the EV sprung mass, inwheel powertrains allow enhanced freedom in the design of the vehicle interiors.
On the downside, in-wheel powertrains significantly increase the ratio between unsprung mass and sprung mass [11], which is detrimental to ride comfort, usually evaluated through acceleration-based performance indicators, and road holding, associated with the dynamic variation of the vertical tire-road contact force [12]- [14]. Extensive literature deals with the attenuation of the effect of road irregularities through active and semi-active suspension controllers [15], and a few studies present dedicated software and hardware solutions to address the unsprung mass increment associated with in-wheel powertrains. For example, in terms of suspension controllers, in [16] Katsuyama et al. propose the so-called 'unsprung negative skyhook' algorithm. In [17] Shao et al. discuss a fault-tolerant fuzzy H∞ suspension control approach for in-wheel motor driven EVs, providing robustness with respect to (w.r.t.) sprung mass variations, actuator failures, and control input constraints. Similarly, the H∞ implementation in [18] considers parameter uncertainties in terms of vehicle body mass, suspension spring stiffness, and suspension damping coefficient, while addressing control input saturation. The performance of the H∞ suspension controllers in [19] and [20] is facilitated by the inclusion of dynamic vibration absorbers (DVAs), which are a rather common technology in the studies on ride comfort with inwheel motors. In [21] a sliding mode set-up controls conventional active suspension actuators, as well as additional actuators located between the unsprung masses and DVA devices. In [22] and [23] the in-wheel motor is considered as an active DVA, which is isolated from the unsprung mass by a spring and a controllable damper. In [24] Ma et al. use a linear quadratic regulator to control an actuator between the in-wheel motor and the unsprung mass, in parallel with a passive spring and damper, with a cost function considering the vertical accelerations, relative displacements as well as the dynamic tire load. A similar philosophy is followed in [25], which uses skyhook control and two regenerative actuators on each corner. In [26] dedicated actuators compensate the effect of the vertical electromagnetic excitation forces of in-wheel machines, which are particularly significant in switched reluctance motors, see the analyses and remedial solutions in the form of suspension or motor controllers in [27]- [31].
Road irregularities also excite longitudinal tire force vibrations (on an irregular road the contact force inevitably has a longitudinal component in addition to the normal component), which are transmitted by the tire to the unsprung mass, and then to the sprung mass through the suspension arms and bushings, see the simulations in [32]- [34], where [33] and [34] include the experimental validation of the respective models.
Only a few authors have discussed this topic in detail, and in particular the option of using the powertrain torque to smoothen the related longitudinal acceleration oscillations. In [35] Bakirci et al. use a proportional integral controller; however, the benefits of their implementation are negligible as the application is an internal combustion engine driven vehicle. In the pioneering study in [36], Fukudome proposes a proportional controller based on the longitudinal speed difference between the unsprung and sprung masses of an EV with in-wheel motors. A control structure consisting of feedforward and feedback contributions is presented in [37], where the road input is an unknown disturbance. The controller is tuned with the longitudinal dynamics model validated in [34], and the feedback contribution includes a deadtime compensator observer to provide robustness against communications delays. Although the method could be beneficial to the suppression of longitudinal acceleration oscillations caused by any source, the results focus on the suppression of vibrations caused by the inwheel motor torque transients, rather than by the road irregularities. In [38] and [39] Walz et al. present a feedforward controller based on a two-point tire model to reduce the longitudinal acceleration oscillations over a known stepped road obstacle during very low speed (less than 10 km/h) maneuvering. In [39] a low level linear model predictive controller computes the engine and friction brake actuation to achieve the torque demand specified by the feedforward controller.
In conclusion, the previous analysis of the literature highlights a gap in model-based and predictive feedback control structures using the powertrains for the compensation of the longitudinal vehicle dynamics effect of road irregularities, starting from the information of the road profile ahead, which is already used in production vehicles for the operation of road preview suspension systems [40]- [41]. A fortiori, the literature misses simple yet accurate model formulations for model-based control design. This paper addresses the identified gap through the following contributions: • A novel proof-of-concept nonlinear model predictive controller (NMPC) formulation for a four-wheel-drive EV with in-wheel direct drive powertrains, using the road preview information to compensate for the longitudinal acceleration oscillations induced by road irregularities. • The preliminary performance comparison of the proposed controller with the formulations in [36] and [37], excluding road preview, by using an experimentally validated simulation model. The remainder of the paper is organized as follows: Section II describes the simulation environment; Section III deals with the novel NMPC formulation; Section IV discusses the selected performance indicators, the benchmarking controllers, and the controller tuning routine; Section V analyzes the simulation results; finally, Section VI summarizes the main conclusions.

A. Simulation environment
The simulation environment of this proof-of-concept study consists of the following functional blocks, reported in the schematic in Fig. 1: • The driver model, which, in the context of the purely longitudinal dynamics analyses with symmetric road profiles of the considered maneuvers, defines the accelerator pedal position, . • The drivability layer, which computes the total electric motor torque at the EV level, , , through an appropriate map, function of and vehicle speed . • The front-to-total wheel torque distribution controller, which defines the individual reference motor torque values, , , , e.g., based on powertrain efficiency or vehicle dynamics criteria. In the previous notation and in the remainder, the subscript ' ' indicates the front or rear axle, while the subscript ' ' indicates the left or right vehicle side. In the preliminary implementation of this paper, focused on the longitudinal dynamics, the electric motor torque values are considered to be the same on the left and right EV sides.
• The longitudinal vibration controllers, which are the core of the study, and output the modified electric motor torque values, , , , , see the detailed discussion in Section III. Each controller is expected to be implemented within the inverter assembly of the respective in-wheel machine, which also facilitates the reduction of the pure time delays of the system, given the fact that the motor sensors are usually directly hard-wired to the inverter. • The nonlinear vehicle model of the plant, used for control system assessment. This is implemented in Matlab-Simulink, and considers the electromagnetic torque dynamics of the in-wheel machines through a transfer function formulation. The model includes consideration of the sprung and unsprung mass dynamics, the nonlinearity of suspension springs and dampers, and the longitudinal compliance characteristics related to the suspension bushings. For accurate tire-road contact modeling on irregular road surfaces, the vehicle model is interfaced with the well-known MF-Swift tire model [42]- [43], i.e., an extensively experimentally validated semi-empirical tire model, which is commercially available as a software package. MF-Swift uses the Pacejka magic formula for the computation of the tangential tire forces and moments, which is complemented by: a) a three-dimensional enveloping model, which implements geometric filtering to describe the capability of the tire to deform when rolling over road irregularities, by considering rows of elliptical cams vertically moving in accordance with the specific road profile; b) a rigid ring model of the tire structure, in which the belt (the rigid ring) is connected to the rim by means of springs and dampers, and its inertial, centrifugal, and gyroscopic effects are taken into account. This arrangement permits to accurately represent the dynamic behavior of the tire in the frequency range up to 100 Hz, in which the bending modes of the tire belt can still be neglected. The rigid ring model has been extensively validated by its developers through experimental measurements from rolling tires; and c) residual stiffness and damping models. These have been introduced between the contact patch and the rigid ring to ensure realistic modeling of the total quasistatic tire stiffness in vertical, longitudinal, lateral and yaw directions. In summary, the total tire model compliance is based on the carcass (ring suspension) compliance, the residual compliance (in reality a part of the total carcass compliance), and the tread compliance.

B. Case study vehicle and model validation
The case study EV is a sport utility vehicle, used as one of the demonstrators of the European Horizon 2020 EVC1000 project [44]. The main EV parameters are reported in Table I. The baseline version of the considered EV, i.e., prior to the installation of the in-wheel motors, was equipped with vertical and longitudinal accelerometers, and was experimentally tested on a typical ride comfort road. The time domain results were converted into the frequency domain, and are expressed in terms of power spectral densities (PSDs) in Fig. 2, which reports the vertical acceleration (̈, ) of the front left unsprung mass, and the vertical (̈, ) and longitudinal (̈, ) accelerations of the front left corner of the sprung mass.   implement the controller in a bespoke inverter unit, the pure time delays between the generation of the reference torque and the initial inverter reaction can be neglected. The good match between simulations and experiments makes the model a reliable tool for the preliminary assessment of control system performance.

C. Preliminary simulations: effect of road irregularities
To analyze the effect of road irregularities in terms of vertical and longitudinal acceleration oscillations of the sprung mass, the passive configuration of the considered case study EV is simulated on a typical ride comfort road profile. The results, reported in Fig. 4, show that the ratio between the peaks of the vertical and longitudinal accelerations ranges approximately from 1 to 2. This confirms that, although smaller than the vertical acceleration oscillations, the ̈ oscillations are not negligible, and thus have a perceivable impact on vehicle comfort. As the current ride comfort enhancement solutions are exclusively based on the compensation of the vertical components of the sprung mass acceleration, the potential benefit of this research, focused on the compensation of the longitudinal vehicle dynamics excited by road irregularities, is evident.

III. NONLINEAR MODEL PREDICTIVE CONTROL FORMULATION
A. Control structure Since model predictive control embeds, by definition, a prediction of the future behavior of the system, similarly to linear quadratic control, it is the ideal control method to implement road preview controllers, as highlighted in a recent review paper on suspension control based on road preview [41]. Moreover, the specific system of this study includes significant nonlinearities, the main ones being the longitudinal force / slip behavior of the tires, as well as the force / speed characteristics of the suspension shock absorbers. Also constraints are a very relevant limiting factor, as the control action should account for the motor torque limits. The concurrent requirements for achieving preview, accounting for system nonlinearities and formally considering constraints make NMPC the ideal control structure for the specific application of this study.
For each vehicle corner, the nonlinear model predictive control structure, see Fig. 1, includes: • A simplified two-dimensional tire enveloping model. Its input is the vertical road displacement vector along the prediction horizon, , = [ ,0 ,1 … , … , ℎ ], where is an individual vertical road displacement, the final number in the subscript of indicates the position of the element along the NMPC prediction horizon, and ℎ is the number of steps of the prediction horizon. , is generated from a spatial road profile map, output by the road preview sensor under the assumption of constant vehicle speed. In this respect, differently from the usual case of active or semi-active suspension control with road preview, the actuators, i.e., the specific direct drive in-wheel motors, are much more responsive, with a measured time constant of only ~6 ms (see Fig. 3), which is lowerby an order of magnitudethan the typical time constant of active or semiactive suspension actuators. This means that a very short prediction horizon, of the order of magnitude of a few tens of ms, against the hundreds of ms of the case of a controlled suspension system, is sufficient to provide the specific inwheel electric powertrains with preview information that can be used to effectively compensate the longitudinal vibrations of the vehicle body caused by the road irregularities. Such prediction horizon values make any consideration of the variation of vehicle speed along the horizon nearly irrelevant, in the on-line computation of the road profile ahead. Nevertheless, an assumption of constant acceleration could be adopted in future applications with less responsive powertrains. The enveloping model produces the vectors of the effective vertical road displacement and gradient along the prediction horizon, namely indicates an individual effective road gradient.
• The weight scheduling block, which adapts the NMPC cost function weights, , and , to the current operating condition of the EV.
• The implicit NMPC controllers, one for each EV corner, which calculate the online solution of the nonlinear optimal control problem, based on a prediction model of the EV corner dynamics and the defined set of external inputs, according to a distributed structure. The following subsections include the detailed description of the building blocks of the NMPC structure.

B. Enveloping model
The tandem enveloping model with elliptical cams, proposed by Schmeitz in [45], is used in the controller implementation, as it represents a good trade-off between accuracy and computational effort. For a given wheel position, the model outputs the effective road profile, defined by and (the subscripts referring to a specific vehicle corner are omitted for simplicity of notation), which is calculated outside the NMPC and introduced in the controller as an external input. The tireroad contact model, see Fig. 5, consists of two identical rigid ellipses that move together with the wheel in the longitudinal direction without rotations, while they can independently translate vertically according to the road profile, and remain at a constant horizontal distance, . The effective road profile is calculated at 0 , which is aligned with the wheel center.  (1) and (2), starting from the vertical positions of the centers of the front and rear ellipses, , and , , expressed in the reference system in Fig. 5: where is the vertical semi-axis of the ellipses. By rearranging the equation of the front ellipse: where is the horizontal semi-axis of the ellipse, is the ellipse shape parameter, and , , is the local axis system of the front ellipse, the vertical distance between the center of the front ellipse and its bottom boundary, , can be expressed as a function of , ∈ [− , ]: Similar formulations are used for the computation of the corresponding distance, , for the rear ellipse. For a given longitudinal wheel coordinate , , and , are obtained as the maximum values of the sum of the road height, , and the distances and , across the possible range for , and , : where is a function of the longitudinal wheel position and the considered longitudinal position along the local axis of the respective ellipse, i.e., = ( , , / ).

C. Prediction model
A novel prediction model, also referred to as internal model, is embedded in the NMPC implementation of each corner, according to the schematic in Fig. 6. Each prediction model includes: a) the tire and wheel dynamics, the longitudinal and vertical unsprung mass dynamics, and the vertical dynamics of the sprung mass, all of them only for the considered corner; and b) the longitudinal dynamics of the sprung mass of the whole vehicle. As a consequence, from the mechanical viewpoint of the dynamics of the involved inertias, the system is characterized by five degrees of freedom, on top of which additional differential equations are considered, e.g., to describe the dynamics of the electric machines, and tire relaxation behavior. The prediction model neglects the longitudinal and vertical dynamics of the unsprung masses of the other vehicle corners. However, the electric motor torque values of the other three corners are sent to the internal model of the considered corner as external inputs, to ensure that the predicted average longitudinal acceleration profile of the sprung mass is realistic, while the impact of the high-frequency disturbance related to the road irregularity is considered only in terms of contribution generated by the specific corner. According to this approach, the apparent mass ( ) is used as the equivalent mass including the effect of the rotating parts of the vehicle that are not modeled in detail within the individual prediction model, i.e., the other wheels (including the brake discs, the rotating components of the hubs, and the rotors of the in-wheel machines). These components increase the overall vehicle inertia, but their dynamics for the other corners are not included in the prediction model for a specific corner.
In summary, in the prediction model concept, the road input affects both the longitudinal and vertical dynamics of the sprung mass. Hence, a possible future development is the consideration of suspension actuators to control the vertical dynamics of the sprung mass, concurrently with the control of its longitudinal dynamics through the innovative system proposed in this paper. The inclusion of the dynamics of a single corner within the prediction model significantly reduces the number of states, control inputs and parameters of the individual controller, w.r.t. a centralized architecture, thus providing computational efficiency, which facilitates the realtime implementation of the resulting implicit NMPC. The unsprung mass includes a central non-rotating part, corresponding to the suspension components, brake caliper and motor stator, and an outer rotating part, corresponding to the wheel rim, brake disc and motor rotor. The vertical stiffness and nonlinear damping elements between the sprung and unsprung masses emulate the suspension spring and shock absorber, while the longitudinal spring and damper, with stiffness , and damping coefficient , , reproduce the compliance effect of the suspension bushings of the considered corner.
The dynamics of the tire structure are modeled through two couples of spring-damper systems, the first one radial, with stiffness , and damping coefficient , , and a second one tangential, with stiffness , and damping coefficient , . The elements of the tire structure, which transmit their forces directly to the non-rotating part of the unsprung mass, can be supposed to be in contact with the road through a frictionless roller system. The longitudinal tire force caused by the longitudinal tire slip is considered to be independently applied to the rotating part of the sprung mass, and then to its nonrotating part through the bearing system. Although providing acceptable accuracy, this set-up significantly simplifies the formulations in comparison with the MF-Swift tire model, used for control system assessment.
The internal model dynamics are described by the following equations, which are reported for the front left corner (the same approach is used for the other vehicle corners): where , is the actual electro-magnetic motor torque, and is the respective time constant.
where is the total sprung mass, and , , and , , are the vertical suspension forces, associated with the spring and shock absorber. , , is given by: where , is the vertical suspension stiffness, and , is the vertical displacement of the unsprung mass. In the specific implementation, , is calculated with a continuous nonlinear function that approximates the passive damper behavior: where 1, , 1/2, , 1/2, , and 1/2, are constant coefficients.
where , is the relevant unsprung mass. , , and , , are the vertical components of the radial and tangential forces, , and , , associated with the tire structure, which are given by: where , / is the laden wheel radius. Although the longitudinal acceleration refers to the whole sprung mass, the subscript ' ' in ̈, indicates that (14) is expressed from the perspective of the specific corner, which, in its internal model formulation, is the only one transmitting its traction force to the sprung mass through the suspension bushings, while, for computational efficiency, the other corners apply their traction force directly to the sprung mass. , , and , , are the forces associated with the longitudinal compliance of the suspension system, and are given by: , , = , [̇, −̇, ] (15) and (16) refer to the specific corner targeted by the internal model, whichinsteadneglects the unsprung mass dynamics of the other three corners. In fact, these are considered by the respective internal models, according to the distributed control architecture in Fig. 1.
The aerodynamic drag force, , is calculated as: where is the air density, is the aerodynamic drag coefficient, and is the frontal area of the EV. The front and rear axle rolling resistance forces, , and , , are > REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) < 7 expressed as: where is the rolling resistance coefficient, given by: with 0 and 2 being constant coefficients. is the apparent sprung mass for the longitudinal force balance of the considered corner model, and includes the unsprung masses and mass moments of inertia ( , / ) of the rotating parts of the other corners: where , , and , , are the longitudinal components of where ̇ is the angular wheel acceleration, and , is the rolling resistance moment.
where , is the relaxation length, , is the delayed tire slip ratio, and ℎ, is the theoretical slip ratio, which is calculated as: The NMPC internal model, described by (7)- (25), is rearranged as a system of first order differential equations: where indicates a nonlinear function; is the state vector; = Δ , is the control action, i.e., the motor torque correction imposed by the longitudinal vibration controller, with , , , = , , + Δ , ; and is the parameter matrix.
is defined as: while includes the expected effective road profile, generated by the enveloping model, and the uncorrected torque demands: Further parameters could be added in future implementations, e.g., the estimated tire-road friction coefficient.

D. Road preview
The NMPC internal model requires the effective road height and slope information from the enveloping model along the preview time = , which is the same as or shorter than the prediction horizon ℎ = ℎ , where is the number of preview points, and is the time step, see Fig. 7. Since in the specific application ℎ is short, vehicle speed can be considered constant along ℎ , and the wheel position is given by: where the subscript indicates the time step. If ℎ > , the effective road data estimated for ℎ is maintained constant for > , i.e.,

E. Nonlinear optimal control problem
The discrete form of the nonlinear optimal control problem formulation has the following structure, where the subscripts referring to the specific corner are omitted for simplicity of notation: is the vector of the control sequence along the prediction horizon; is the cost function, consisting of two terms, the first one aiming to minimize the response error at the end of the prediction horizon, and the second one targeting the response optmization along the prediction horizon; is the vector of the predicted system outputs, whose corresponding reference vector is ; and are positive diagonal matrices used for weighting; is the initial value of the state vector; is the discretized version of the model in (26); is the function expressing the system outputs; and and are the bounds of the control action . In the uncorrected torque demands are considered constant along the prediction horizon. In the specific implementation, includes the longitudinal acceleration of the vehicle body, which is the main variable related to longitudinal vehicle comfort, and the relative longitudinal speed between the unsprung mass of the considered corner and the sprung mass, which is at the origin of the longitudinal acceleration oscillations: while , is defined as: ̈, , is the reference acceleration, which is calculated with the following rigid vehicle model formulation: (33) is the result of the longitudinal force balance at the vehicle level, where the effect of road irregularities on the longitudinal acceleration profile is purposely neglected, i.e., the resulting acceleration is the longitudinal acceleration the vehicle would have in absence of road irregularities, which, given the scope of the proposed controller, is the required reference acceleration value. In particular, in (33) the motor torque values divided by the respective wheel radius represent the traction forces, from which the tire rolling resistance force and aerodynamic drag force are subtracted, before being divided by the total apparent mass of the vehicle. The lower and upper constraints in the control action are related to the motor torque limitations: where and are the lower and upper boundaries of the motor torque.
The controllers are set up through the ACADO toolkit [46], which uses Gauss-Newton iteration algorithms for fast NMPC with constraints, generating locally optimal solutions [47]. The selected solver parameters are: multiple shooting discretization method, fourth order Runge Kutta integrator, and qpOASES QP optimization algorithm. Only the first control input (or control move) is applied to the plant, according to the receding horizon approach.
The NMPC implementation requires a set of sensors and a state estimator for the real-time generation of . For example, the longitudinal and vertical velocities and displacements of the sprung and unsprung masses can be obtained through an appropriate state estimator, mainly based on the measured longitudinal and vertical accelerations of the bodies as well as the wheel speeds, as in [36], with similar methodology to that already in use for estimating the relevant vertical velocities for suspension control, starting from the vertical accelerations, e.g., see the benchmarking skyhook controller in [40]. The required state estimator for the operation of the NMPC, based on a model of the system similar to the prediction model, will be the object of a dedicated publication.

A. Considered road profiles
The results of this proof-of-concept study are based on the following road profiles: • Road profile 1, i.e., a 20 mm positive step.
• Road profile 2, i.e., a 55 mm high speed bump defined by [48]: where ℎ is the height of the bump, is its half width, and is the longitudinal coordinate of the road profile. • Road profile 3, which corresponds to a section of a ride comfort assessment road, normally used as a reference by the industrial partners of the project.

B. Key performance indicators
The following set of key performance indicators (KPIs) is used to evaluate the controllers:  • Δ̈, , i.e., the maximum longitudinal acceleration deviation from its reference value: • , i.e., the integral of the absolute value of the control action on the EV corners, averaged with time:

C. Benchmarking controllers
In the preliminary analysis of this study, the performance of the proposed road preview controller, referred to as NMPC (prev) in the remainder, has been compared with that of: • The passive configuration, excluding any form of longitudinal vibration control. • A version of the proposed NMPC excluding the road preview information, indicated as NMPC (w/o prev) . • The feedback controller in [36], abbreviated as FB, which is based on a torque correction, Δ , , proportional to the relative longitudinal speed between the sprung mass and unsprung mass of the considered corner, according to a virtual damping coefficient, ℎ, : • The combination of feedforward and feedback controllers (indicated as FF+FB) in [37], see Fig. 8, for which the road profile is an unknown disturbance. The input to the feedforward contribution is , , , which is converted into a local reference acceleration, ̈, , , through the gain . ̈, , is multiplied by: i) a low-pass filter, ; and ii) the transfer function of the inverse model of the system, i.e., ( )/̈( ), where is the equivalent traction force corresponding to the motor torque, and is the -transform operator. In the feedback loop, the estimated system states, ̂, are provided to the state feedback controller, , by a dedicated observer.

D. Controller tuning routine
The main tunable parameters of the NMPC (prev) algorithm are the cost function weights defined in the positive definite matrices and . The controller tuning was carried out through the genetic algorithm of the gamultiobj function of Matlab, which was coupled with the vehicle model for control system assessment, to minimize the cost function , considering the weighted average of ̈ and ̈. The NMPC weight optimization problem is defined as: where 1 and 2 are constant weights, and , , , , and define the possible range of the NMPC cost function weights. In the following analyses, NMPC (w/o prev) uses the same weights as NMPC (prev) . A similar optimization-based tuning routine was applied to the calibration of the benchmarking FB and FF+FB controllers.

V. SIMULATION RESULTS
The simulations of this section use the nonlinear vehicle model in Section II, including the MF-Swift tire model.

A. Results for optimized configuration for road profile 1
The target of this analysis is to investigate the potential performance of the proposed road preview NMPC when its calibration is optimized for a localized road event, i.e., the step of road profile 1. This approach is practically relevant, since in a smart control system implementation, it would be possible to modify the controller calibration depending on the measured (by the vision sensor) or expected (based on information from the cloud, in case of a connected vehicle) road profile ahead.  As a starting point, the sensitivity of the controller performance to the prediction horizon is assessed in conditions of zero torque demand on the front and rear powertrains. While ℎ is varied from 10 to 50 ms (with 1 ms selected as discretization time of the internal model), the cost function weights are maintained constant and equal to their value optimized for ℎ = 30 ms. Fig. 9 shows the percentage reduction of the selected KPIs w.r.t. their values for the passive EV. A progressive performance improvement is observed as a function of ℎ . For the specific implementation, ℎ = 30 ms represents a good trade-off between vibration attenuation and computational effort, with ~80% reduction of all comfort indicators. Hence, this prediction horizon value is used for all the following simulations. For completeness, to appreciate the impact of the extension of the prediction on the resulting computational effort, Fig. 10 shows the nondimensional computational time ̅ , which is normalized with respect to its value for the case of a single prediction step, as a function of the number of prediction steps, ℎ , with the controller running on a computer with an Intel Core i7-4700MQ 2.40 GHz processor and 16 GB RAM.
The optimization of the cost function weights of NMPC (prev) on road profile 1 was carried out for different individual wheel torque demands and initial vehicle speeds, according to the routine in Section IV, which provided weight scheduling maps. The longitudinal load transfer between the front and rear axles as well as the different suspension compliance characteristics cause a difference in the calibration maps of and for the front and rear controllers. Fig. 11 reports the time profiles of the main variables for , = 0 Nm and ̇, 0 = 40 km/h. In this figure and in the remainder, the plots refer to individual front and rear corners; given the symmetry of the selected road profile, the results are identical on the left and right vehicle sides. The road preview is particularly effective for localized events such as a step or a speed bump, as it allows the controller to prevent the first longitudinal acceleration peak when the car meets the bump, by requesting a control action in advance. On the contrary, given the localized nature of the road obstacle, NMPC (w/o prev) only brings a minor response improvement immediately after the step; the benefit becomes more evident after the first acceleration peak, when NMPC (w/o prev) facilitates the attenuation of the underdamped ̈ oscillations.   (prev) , with RMS values consistently lower than 11 m/s 3 , while the passive configuration always exceeds 30 m/s 3 . Similar comfort improvements are obtained by NMPC (prev) for , = 1200 and 2400 Nm. In summary, the important and novel conclusion is that the road preview controller based on a simplified corner model can achieve a very substantial comfort improvement at low-tomedium torque demands for a wide range of vehicle speeds, with nearly full compensation of the longitudinal acceleration peaks induced by the road disturbance.

B. Controller evaluation on different road profiles
While Section V.A focused on a controller calibration optimized for road profile 1, under the assumption of an implementation capable of adapting its tuning based on the expected road profile ahead, this section presents the results for a controller calibration that is a trade-off across the considered relevant spectrum of road profiles. Therefore, in the following analyses the calibration weights, functions of motor torque demand and speed, remain the same for road scenarios 1-3. Also, with respect to Section V.A, the tire structure parameters in the prediction model, namely , , , , , and , , are modified to achieve a generally good match between the NMPC prediction model and the MF-Swift tire model across a wider range of conditions, at the price of an increased mismatch for road profile 1.
Figs. 12-14 are the resulting time histories of the main relevant variables, for the three considered road profiles, in conditions of zero motor torque demand, from 40 km/h. The figures include the passive configuration, NMPC (prev) , and the benchmarking FB and FF+FB implementations. For road profile 1, the performance decrease of NMPC (prev) in comparison with its calibration discussed in Section V.A is evident, even if the controller still manages to more than halve the longitudinal acceleration peaks that follow the application of the road step input to each axle. In the same conditions, FB and FF+FB provide only a marginal improvement w.r.t. the passive vehicle, and mostly only after the first and largest acceleration peak. Similar trends are visible in Figs. 13 and 14.
The torque profiles of Figs. 12 and 13 highlight the preemptive nature of NMPC (prev) , which starts modifying the torque demand before the other controllers, and for which the torque corrections are localized in the corner involved in the specific road excitation, i.e., the torque correction is first applied to the front corner and then to the corresponding rear corner. This corner-related behavior, typical of a distributed control structure, is visible also for the FB controller, even if the dynamics and delays between the road disturbance and the generation of the longitudinal suspension deflection rate (̇, − , used by the controller for computing the control input), as well as the in-wheel powertrain torque dynamics, compromise the effectiveness of FB. Increased performance of the FB controller was observed in further simulations neglecting the electric motor dynamics. Differently from the other controllers, as FF+FB has a centralized structure based on ̈, its interventions are not directly associated with a road excitation on a specific axle, and therefore its torque corrections are evenly distributed among the front and rear powertrains.
The KPIs for all considered configurations, including NMPC (w/o prev) , are reported in Table III. For FB and FF+FB, the performance improvement in some indicators is often associated with increased values of the other indicators. For example, with respect to the passive case, FB effectively reduces the values of the first two indicators for road profiles 1 and 3; however, it marginally worsens them on the speed bump, and the longitudinal jerk is increased in all considered scenarios. In terms of longitudinal acceleration indicators, FF+FB brings a performance improvement for the localized events of road profiles 1 and 2, while on the ride comfort road the longitudinal acceleration indicator values are higher than for the passive set-up. However, FF+FB is very effective as an antijerk controller, since it decreases the RMS value of the longitudinal jerk by ~25% in all tests.    The conclusion is that the for the specific application and maneuvers, the benchmarking controllers from the literature can only bring limited performance improvement, as they are reactive rather than pre-emptive, and, although being modelbased, they only indirectly rely on model formulations of the system dynamics. Both NMPC versions significantly outperform the passive configuration and the benchmarking controllers across the three maneuvers. In particular, the top subplot of Fig. 15 displays the KPI percentage reduction brought by NMPC (w/o prev) w.r.t. the passive configuration, while the bottom subplot shows the KPI reduction caused by NMPC (prev) w.r.t. NMPC (w/o prev) . The NMPC without preview reduces the acceleration-related indicators from ~15% to >40%, and is especially effective on the ride comfort road. The reduction of the longitudinal jerk is very limited in the step case, because of the nearly impulsive nature of the event and the reactive nature of NMPC (w/o prev) ; however, for road profile 3, the jerk reduction exceeds 45%, which is a major result. NMPC (prev) further improves all indicators by >20% to >50%, depending on the maneuver and the selected KPI. The effectiveness of the proposed NMPC formulations is further confirmed by the values, which are aligned with those of the benchmarking implementations.

C. Robustness to the road preview error
The target of this analysis is to verify through simulations the robustness of the controller w.r.t. estimation inaccuracies of the road profile ahead, and compare the resulting performance with those of NMPC (w/o prev) and the passive configuration, without modifying the calibration of the cost function weights of NMPC (prev) . The first kind of considered road preview information error is modeled as a zero-mean Gaussian white noise process, with its magnitude being expressed in terms of percentage of the average absolute value of the actual road profile deviation from its mean, similarly to the approach in [49]. The effect of the relative magnitude of sensor noise on the performance of the preview controller is reported in Fig. 16 for road profile 3, which shows that NMPC (prev) brings a benefit in all indicators over the corresponding version without preview up to ~25% noise intensity, and over the passive configuration up to ~45% noise intensity. The second typology of considered road preview information error is in the form a time shift of the road profile signal, which can result either into an advance or delay with respect to the correct road profile, with a maximum considered value of 20 ms, see Fig. 17. The results highlight that NMPC (prev) performs better than NMPC (w/o prev) for a magnitude of the time shift up to nearly 10 ms.

D. Real-time implementation
As a proof-of-concept, it was verified that the proposed controller with road preview runs in real-time on a dSPACE MicroAutoBox II system (900 MHz, 16 Mb flash memory), see Fig. 18, with 2 control steps along a 30 ms prediction horizon. The discretization time of the internal model is fixed to 1 ms, which ensures numerical stability without significantly affecting the computational time. For example, with these parametrizations, despite the lower number of prediction steps than in the simulations of the previous sections, NMPC (prev) still reduces the considered acceleration-and jerk-related indicators from >20% to ~45% and from ~33% to ~49% w.r.t. the passive configuration, along road profile 3, respectively from initial speeds equal to 20 and 40 km/h. As much more powerful realtime control hardware is currently available (see the most recent rapid control prototyping units and the computing platforms for automated driving applications), there is significant further margin for improvement of the real-time setting and its performance.

VI. CONCLUSIONS
This study presented a novel proof-of-concept nonlinear model predictive control formulation for electric vehicles with in-wheel motors, to attenuate the longitudinal acceleration oscillations caused by irregular road profiles, through electric powertrain control based on the road preview information. The preliminary results, obtained through simulations with an experimentally validated model, bring the following conclusions: • The novel prediction model concept adopted in the NMPC formulations is sufficiently accurate to provide effective predictive control for the considered range of road inputs and operating conditions. • If optimized for the positive step input, the NMPC formulation with road preview can achieve nearly full compensation of the longitudinal acceleration peaks in conditions of low-to-medium torque demand, with a reduction of the longitudinal acceleration performance indicators by an order of magnitude, with respect to the passive case. Also, the root mean square value of the longitudinal jerk is reduced by a factor >3. • The trade-off calibration of the NMPC formulations across the selected range of road profiles confirms the potential of the proposed controllers. In particular, the controller version without road preview reduces the acceleration-related indicators from ~15% to >40%, while the longitudinal jerk reduction ranges from 2.5% to >45%. The inclusion of road preview further improves all indicators by >20% to >50%. • The comparative analysis with two recent controllers from the literature for the compensation of the longitudinal acceleration oscillations in electric vehicles with in-wheel motors highlights the very significant benefit of the proposed model-based formulations. In fact, because of their reactive nature, in the specific vehicle the benchmarking controllers are unable to simultaneously improve the longitudinal acceleration and jerk performance indicators, and in any case their benefit is rather limited with respect to the passive configuration. • A setting of the proposed algorithm with road preview was preliminarily run in real-time on a rapid control prototyping unit, and is able to reduce the acceleration-and jerk-related indicators from >20% to ~45% and from ~33% to ~49% w.r.t. the passive vehicle, along the considered ride comfort road, respectively from 20 and 40 km/h. The extent of the potential performance improvement associated with the proposed novel controllers, highlighted by the simulation results of this study, could originate a research line targeting the further development of the technology, including: i) the analysis of further and more complex test cases, e.g., with asymmetric road inputs on the two vehicle sides; ii) the evaluation of the required state estimators; and iii) the experimental demonstration on a vehicle with in-wheel powertrains.