Analysis of the Influence of Suspension Actuator Limitations on Ride Comfort in Passenger Cars Using Model Predictive Control

Abstract

Active suspension systems help to deliver superior ride comfort and can be used to resolve the objective conflict between ride comfort and road-holding. Currently, there exists no method for analyzing the influence of actuator limitations, such as maximum force and maximum rate of change, on the achievable ride comfort. This research paper presents a method that is capable of doing this. It uses model predictive control to eliminate the influence of feedback controller performance and to integrate both actuator limitations and necessary constraints on dynamic wheel-load variation and suspension travel. Various scenarios are simulated, such as driving over a speed bump and inner city driving, as well as driving on a country road and motorway driving, using a state-of-the-art quarter-car model, parameterized for a luxury class vehicle. It is analyzed how comfort, or in one scenario road-holding, can be improved with consideration for the actuator limitations. The results indicate that actuator rate limitation has a strong influence on vertical vehicle dynamics control system performance, and that relatively small maximum forces of around 1000 to $2000$ N are sufficient to successfully reject disturbances from road irregularities, provided the actuator is capable of supplying the forces at a sufficiently high rate of change.

1. Introduction

Together, automated and autonomous driving, and new propulsion technologies, make up one of the biggest research fields in the automotive industry. The switch in roles from driver to passenger, offers both opportunities and pitfalls for the near future, at least in some scenarios. Many characteristics of today’s road vehicles are designed to please the driver, as he or she is usually the customer. A typical development methodology is the so called V-Model [1,2], which is defined for mechatronic systems in VDI 2206 [3]. It prescribes a framework for development according to requirements defined for the product, subsystem and component levels. Vehicles thus developed must fulfill all requirements prescribed by law. In automotive software development, ISO 26262 [4] is applied to secure functional safety. Other requirements are more difficult to define, as they are dependent on many internal and external factors, such as customer demand, market strategy, self-set quality standards and price targets. Having set these vehicle-specific requirements on product level, they need to be broken down to subsystem and component level. In order to do this, it is necessary to understand the interdependency between subsystem performance and product performance. Objective measures for quantifying subsystem and product performance are also required. In the case of vehicle dynamics and suspension development for passenger cars, these matters are well-understood thanks to the decades of experience and research carried out by the manufactures and universities [5,6]. However, this all changes when it comes to automated and autonomous vehicles, as there is a lack of accumulated experience, and it is not clear how the customer preferences might change. Previous studies of the authors indicate a demand for increased comfort during an automated drive [7,8,9]. This is in conflict with traditional driving-dynamics requirements where requirements placed on longitudinal and lateral dynamics generally have priority over comfort, except in the case of luxury vehicles. Furthermore, an ideal vehicle offering automated driving functionality must be able to fulfill both automated and non-automated requirements in the best possible manner. Active suspension systems help to resolve this conflict, because they can be tuned to suit the different objectives for each driving mode. However, the precise relationship between active suspension system actuator performance and comfort requirements at vehicle level are currently neither known nor understood, either in academia, or in the automotive industry. This research paper aims to close this research gap.

2. State of the Art

This section gives an overview of the current state of the art. Section 2.1 outlines vertical dynamics components their properties. Section 2.2 summarizes the objectification methods, while Section 2.3 and Section 2.4 provide an overview of passive, semi-active and active suspension systems and the current research into model predictive control in the context of vertical vehicle dynamics control. Finally, the gap in research is outlined in Section 2.5.

2.1. Suspension Components

The most important suspension components responsible for vertical vehicle dynamics are the springs and dampers, because it is these that have a major effect on the vehicle’s response to road irregularities and the reaction of the vehicle body to longitudinal and lateral accelerations. In combination with the mass of the vehicle body, the unsprung mass and the vertical stiffness of the tires, they define the frequency response of the vehicle in the range of $0\mathrm{Hz}$ to $20\mathrm{Hz}$, which includes the body-mode eigenfrequency and the wheel-mode eigenfrequency. A more detailed description of this is given in Section 2.3. The general characteristics of passive and (semi-)active suspension springs and dampers will be explained in the following. Further information can be found in [10] (pp. 426–474).

The majority of main suspension springs fitted to road cars are either coil-, leaf- or air-type. Most active springs are air-type that can either switch between different airchambers to change suspension stiffness [11], or are controlled by air-pressure, mainly for ride-height adjustments in a low-bandwidth system [5]. An example of a high-bandwidth active spring system for controlling vehicle response is presented in [12]. Anti-roll bars (ARBs) are another type of springs fitted to vehicles, which connect the wheels (unsprung masses) of an axle to reduce body roll under lateral acceleration. The wheels are coupled, and any road-excitation from one side is transferred to the other. Active ARBs were developed as a means of providing roll stiffness under lateral acceleration, varying the roll-stiffness distribution, and minimizing coupling during straight-line driving. Examples of such systems are described in [13,14]. Besides springs, dampers are another integral part of a vehicle suspension. The design and working principles of suspension dampers are described in detail by Dixon [15]. As will be shown in Section 2.3, there is a conflict between ride comfort and road-holding. Electronically controlled, semi-active continuously variable dampers (CVD) were invented in an attempt to resolve this conflict. CVDs are usually magnetorheological [16] or electromechanical [17] systems. An implementation of a control system for CVDs is described in [18]. Active dampers are able to generate a desired force independently of the relative velocity of the damper, as introduced by Daimler [19]. An overview of passive, semi-active and active suspension dampers is given in [17]. There are also other active elements that cannot be directly classified as springs, dampers or ARBs, for example the fully active system developed by Lotus [20], or the system introduced by Ovalo and Audi, which can be installed in place of an ARB [21].

2.2. Objective Performance Measures

To quantify the performance of a vehicle suspension by simulation, it is necessary to define objective performance indicators that reflect human perception while driving a vehicle. In the context of this research, objective measures for vertical vehicle dynamics are needed. In terms of comfort, the main task of the suspension system is to insulate the vehicle occupant from road vibration induced by road irregularities. In order to quantify suspension performance, it is therefore necessary to quantify the vibration of the car body for a given road excitation, as it is this vibration that is transferred to the vehicle occupant through the seat, and it is also the output quantity that is solely dependent on the performance of the suspension system.

Human sense vibration through their visual, vestibular, somatic and auditory systems [22]. For measuring ride comfort in the context of vehicle suspension, the frequency range of interest is between $0\mathrm{Hz}$ to $20\mathrm{Hz}$ as stated in Section 2.1. As this is under the auditory threshold, the auditory system can be neglected in terms of suspension ride comfort, contrary to the field of noise, vibration and harshness. The mechanism of perception is still complex and there is a genuine difference between comfort and discomfort. According to Looze et al. [23] comfort has not yet been clearly defined, and the exact relationship between comfort and discomfort is still subject to debate. One theory is that comfort and discomfort are independent factors, and so a reduction in discomfort will not automatically yield an increase in comfort [24]. Furthermore, discomfort is a result of several factors, of which vibration is only one [25] (p. 44). The automotive term ’ride comfort’ is therefore imprecise, because it is usually suspension vibration discomfort that is the subject of analysis. Other influences, such as noise, can still have an unwanted effect on the subjective rating of suspension vibration discomfort in real vehicles [26].

There are three standards by which vibration discomfort is generally assessed, these being ISO 2631 [27], BS 6841 [28] and VDI 2057 [29]. They all apply the same method of using acceleration measurements to calculate a weighted root mean square (RMS) acceleration comfort value. Weighting is performed in the frequency domain with factors representing the frequency sensitivity. RMS acceleration values for each direction are calculated from the weighted signals and combined to a discomfort value using directional weighting factors representing the directional sensitivity of human perception. VDI 2057 [29] extends the method of ISO 2631 [27] for vibration containing shocks, while BS 6841 [28] uses slightly different frequency weighting factors. The authors performed a comparison of objectification methods for ride comfort. The results indicate, that the standardized methods are more robust and reliable than methods developed by other researchers [9]. All of these methods have in common that they objectify the ride comfort (vibration discomfort) as perceived by an untrained customer. Development engineers specializing in vertical vehicle dynamics have a different perception because of their background knowledge. They try to identify such different vibration phenomena as wheel hop, engine shake and copying and then correlate them to different parts and systems of the car. An overview of subjective development criteria is given in [30] (pp. 530–531).

2.3. Passive Suspension Systems

The quarter-vehicle model is used to describe the vertical dynamics of a passenger car. Here, the full vehicle is reduced to one corner and modeled as a two-mass spring-damper system, as illustrated in Figure 1. The mass $m_{\mathrm{Bo}}$ represents the body, $m_{\mathrm{Wh}}$ the unsprung mass. The two are connected by a spring and a damper representing the suspension stiffness $c_{\mathrm{Bo}}$ and suspension damping $d_{\mathrm{Bo}}$. The wheel and the road are also connected by a spring and a damper, where $c_{\mathrm{Wh}}$ is the tire stiffness and $d_{\mathrm{Bo}}$ the tire damping. The displacements of the two masses and the road-profile are defined as $z_{\mathrm{Bo}}$, $z_{\mathrm{Wh}}$ and $z_{\mathrm{Ro}}$. For a longitudinally symmetric vehicle, the quarter body mass is half of the corresponding axle weight, which itself depends on the front-to-rear weight distribution. The unsprung mass represents the mass of the wheel assembly. Tire stiffness is usually one magnitude higher than suspension stiffness, and tire damping is virtually negligible. The physical derivation of this model is described in [30] (pp. 295–299). This model has only two degrees of freedom and cannot model the pitch and roll behavior of vehicles. However, it is sufficient for describing the essential vertical dynamic behavior of a vehicle, which is characterized by the objective conflict between ride comfort, suspension travel and road-holding [31] (p. 112). The advantage of this model is its simplicity as a linear lumped parameter model, which makes it well suited to be used in control systems and for simulation purposes.

The conflict between ride comfort and road-holding is best described in the frequency domain. Figure 2 shows the frequency response of a quarter-vehicle model in terms of dynamic wheel-load and body acceleration to road excitation for different suspension damping coefficients. The first peak of around $1\mathrm{Hz}$ corresponds to the body eigenmode, the second peak around $12\mathrm{Hz}$ to the wheel eigenmode. Low damping results in a high response amplitude of the body acceleration at the body eigenfrequency and low response amplitude between the invariant points [32]. For high damping, the opposite applies. The response of the dynamic wheel-load with low damping is high at the body and the wheel eigenfrequencies. Because humans are sensitive to frequencies between $4\mathrm{Hz}$ and $10\mathrm{Hz}$, it is desirable to have low suspension damping, but this results in high dynamic wheel-load variations, which in turn results in bad road-holding, because the side force potential of the tire decreases with increasing wheel-load variations. Higher damping would be beneficial for good road-holding, but this entails higher body accelerations and less ride comfort. Active systems can resolve this conflict but still have to find a balance between these objectives [33] (p. 267).

2.4. Semi-Active and Active Suspension Control

To resolve the conflict between ride comfort and road-holding, a suitable control algorithm needs to be implemented that outputs the required forces, which in turn need to be supplied by the actuators. Suspension systems can be classified as, passive, slowly-variable/adaptive, semi-active and active systems [34]. In the following, research on control strategies for semi-active and active systems will be summarized, including the application of model predictive control.

The first concept for the vertical dynamics control of active vehicle suspensions was developed by Karnopp [35]. This approach was later extended to semi-active force generators [36]. Ten years after the introduction of the skyhook principle, a modal control concept for fully active suspensions was developed by Wright and Williams [20]. A comparison of SHC and LMC, as well as more information on each controller can be found in [37]. Linear optimal control theory, developed by Anderson [38], was also applied to vertical vehicle dynamics control [39,40,41]. A thorough literature review on the application of LQ-optimal control can be found in [42]. Another important research topic is the use of preview information for suspension control, which is investigated in [43,44,45,46]. Furthermore, model and parameter uncertainties can be incorporated by using ${\mathrm{H}}_{\infty }$-control [47,48]. An example of adaptive non-linear control is given in [49]. The skyhook control logic is extended to triple-skyhook control in [50], while an alternative to the skyhook principle is acceleration-driven damping, which was developed by Savaresi et al. [51,52]. General summaries of control methods for semi-active suspensions can be found in [53,54] and for active suspensions in [55]. An analytic derivation of performance constraints and limitations in vertical vehicle dynamics control is given in [56].

The latest research in the field focuses on new signal differentiation methods [57], improved vehicle stability control with active suspension systems [58], optimal control application [59,60], preview strategies [61,62,63], incorporating physiological and psychological effects [64], adaptive controllers [65,66,67,68] and fault tolerant control [69]. A multi-objective optimization method for semi-active suspensions is shown in [70]. A non-linear controller considering actuator and state constraints is developed in [71]. Actuator time delays are included in [72,73], while actuator dynamics are introduced in [74,75,76]. Clipped optimal control for vehicle suspensions is presented in [77]. Another way to incorporate system and actuator constrains is model predictive control.

Model predictive control (MPC) can be regarded as an advancement of optimal control. An overview of the history and industrial application of MPC is given in [78]. The first application of MPC to vertical vehicle dynamics control of which the authors are aware is presented in [79,80]. Robust MPC [81] and fast MPC [82,83] were applied later. A comparison of MPC to clipped optimal control is given in [84]. One of the main problems of real-time MPC in vertical vehicle dynamics control is the requirement of optimization convergence within sampling times of 1–10 ms and the required quality of the road-profile preview. A recent implementation of MPC in low-bandwidth suspension systems is presented by Göhrle [85,86]. He also compares the performance of MPC to other control strategies [61,87] and shows different methods of road-profile estimation [62]. For real-time capable implementations explicit MPC is used, as presented in [88,89,90]. Other recent publications are focusing especially on MPC for semi-active suspension [91,92]. Fault estimation for electro-rheological semi-active dampers is present in [93].

2.5. Research Gap

Recapitulating the previous sections, it is clear that a lot of effort has been invested in developing various kinds of control strategies for active and semi-active suspension systems. A number of different suspension actuators exist, and research has been carried out to understand actuator system properties and performance. However, the current state-of-the-art lacks information on the effect of actuator system properties, such as maximum force and maximum rate of change, on achievable ride comfort. This is not trivial to investigate, because the control system performance depends on actuator performance in conjunction with the performance of the control logic. To fill this research gap, the authors present a new method in this article, which makes it possible to quantify the effect of actuator limitations on achievable ride comfort, while assuring best possible control logic performance through the application of MPC. This approach allows to deduct the actuator characteristics required to achieve desired ride comfort targets and displays the sensitivity of ride comfort with respect to these actuator attributes. It shows the pure influence of the actuator properties for the best possible control input. To the best knowledge of the authors, no method and results comparable to this have previously been published in the literature.

3. Method

This section describes the method used to investigate the influence of actuator limitations on achievable ride comfort. Section 3.1 introduces the simulation model, while Section 3.2 describes the different test scenarios and Section 3.3 presents the data processing and analysis procedures. The simulations were performed in MATLAB© using the Model Predictive Control Toolbox. The underlying theory of MPC can be found in [94,95,96]. The simulation can also be built with non-proprietary software such as Octave or Python using the equations given in Section 3.1 and a suitable solver for convex cost functions.

In optimal control, a feedback law based on the system properties and disturbances is derived for an infinite time horizon. In model predictive control, this time horizon is finite and set to a certain preview time referred to as prediction horizon (PH). The control horizon defines the number of preview steps for which an optimal control input is calculated. Whereas in optimal control, a constant feedback matrix is derived for the controller, in MPC the optimization problem is solved in each time step through the PH. At the end of a time step, the optimal control for the next time step is outputted as control demand to the actuator. In the next time step, the system states and PH are updated and the optimization is solved again. This is referred to as receding horizon optimization. As it involves the repeated calculation of an updated feed-forward control problem, MPC is a form of feedback control, because it always starts from measured or estimated states that are independent of model inaccuracies and dependent on non-modelled disturbances. Convergence of the optimization solution within one sampling interval is a prerequisite for real-time implementation of MPC, which is hard to achieve and a reason for research work into explicit MPC. The performance of MPC in a real vehicle is also highly dependent on quality of the preview information, which has to be estimated from sensor measurements of the vehicle (camera, lidar, etc.), or needs to be supplied from a back-end (cloud) for the current position of the vehicle. In this research work here, MPC is used as a tool to obtain the theoretically best solution, considering given actuator limitations. There is no need for real-time feasibility, and preview information as well as system states are known exactly in all the simulations. Detailed information about the structure of the MPC in this research is given in Section 3.1.3.

3.1. System Modelling

This subsection introduces the different parts of the simulation model. First, the vehicle model and the derivation of the discrete-time state-space model is explained in Section 3.1.1, the generation of the stochastic road-profile and the single-obstacle are presented in Section 3.1.2. The synthesis of the MPC controller including the derivation of the cost function, is described in Section 3.1.3. Finally, the implementation of the actuator limitations (input constraints) and boundary conditions (output constraints) is presented in Section 3.1.4, along with the weighting of the optimization objectives.

3.1.1. Vehicle Model

The equations of motion for the quarter-vehicle model in Figure 1 are shown in Equation (1). The model has two degrees of freedom, which results in two equations of motion.

$$\begin{array}{cc}\hfill m_{\mathrm{Bo}}{\ddot{z}}_{\mathrm{Bo}}& =-c_{\mathrm{Bo}}\left(z_{\mathrm{Bo}}-z_{\mathrm{Wh}}\right)-d_{\mathrm{Bo}}\left({\dot{z}}_{\mathrm{Bo}}-{\dot{z}}_{\mathrm{Wh}}\right)+F_{\mathrm{Ac}}\hfill \end{array}$$

(1a)

$$\begin{array}{cc}\hfill m_{\mathrm{Wh}}{\ddot{z}}_{\mathrm{Wh}}& =c_{\mathrm{Bo}}\left(z_{\mathrm{Bo}}-z_{\mathrm{Wh}}\right)+d_{\mathrm{Bo}}\left({\dot{z}}_{\mathrm{Bo}}-{\dot{z}}_{\mathrm{Wh}}\right)-c_{\mathrm{Wh}}\left(z_{\mathrm{Wh}}-z_{\mathrm{Ro}}\right)-d_{\mathrm{Wh}}\left({\dot{z}}_{\mathrm{Wh}}-{\dot{z}}_{\mathrm{Ro}}\right)\hfill \end{array}$$

(1b)

These can be formulated as a continuous-time state-space model (Equation (2)) with system matrix ${\mathit{A}}_{\mathrm{c}}$, input matrix ${\mathit{B}}_{\mathrm{c}}$, transition matrix ${\mathit{C}}_{\mathrm{c}}$, feedthrough matrix ${\mathit{D}}_{\mathrm{c}}$, state vector $\mathit{x}\left(t\right)$, input vector $\mathit{u}\left(t\right)$ and output vector $\mathit{y}\left(t\right)$.

$$\begin{array}{cc}\hfill \dot{\mathit{x}}\left(t\right)& ={\mathit{A}}_{\mathrm{c}}\mathit{x}\left(t\right)+{\mathit{B}}_{\mathrm{c}}\mathit{u}\left(t\right)\hfill \end{array}$$

(2a)

$$\begin{array}{cc}\hfill \mathit{y}\left(t\right)& ={\mathit{C}}_{\mathrm{c}}\mathit{x}\left(t\right)+{\mathit{D}}_{\mathrm{c}}\mathit{u}\left(t\right)\hfill \end{array}$$

(2b)

Equation (2a) is the state equation and Equation (2b) is the output equation. The matrices and vectors of Equation (2) can be found in Appendix A. This form of the state-space representation with the chosen states was also used in [65]. The model is transformed into a discrete-time state-space model with sample time $t_{\mathrm{s}}$, according to Equation (3).

$$\begin{array}{cc}\hfill {\mathit{A}}_{\mathrm{d}}& =\exp \left({\mathit{A}}_{\mathrm{c}}{t}_{\mathrm{s}}\right)\hfill \end{array}$$

(3a)

$$\begin{array}{cc}\hfill {\mathit{B}}_{\mathrm{d}}& =\left[{\mathit{A}}_{\mathrm{d}}-{\mathrm{I}}_4\right]{\mathit{A}}_{\mathrm{c}}^{-1}{\mathit{B}}_{\mathrm{c}}\hfill \end{array}$$

(3b)

$$\begin{array}{cc}\hfill {\mathit{C}}_{\mathrm{d}}& ={\mathit{C}}_{\mathrm{c}}\hfill \end{array}$$

(3c)

$$\begin{array}{cc}\hfill {\mathit{D}}_{\mathrm{d}}& ={\mathit{D}}_{\mathrm{c}}\hfill \end{array}$$

(3d)

Another option would be to use the forward Euler method to transform the model. The delay of the actuator force is absorbed into the state-space model by defining a new state vector containing the delayed actuator force, and a new input vector containing the actuator force which is calculated in the current step.

$${\mathit{x}}_k^{*}=\left[\begin{array}{c}{z}_{\mathrm{Bo},k}-z_{\mathrm{Wh},k}\\ {\dot{z}}_{\mathrm{Bo},k}\\ z_{\mathrm{Wh},k}-z_{\mathrm{Ro},k}\\ {\dot{z}}_{\mathrm{Ro},k}\\ F_{\mathrm{Ac},k-1}\end{array}\right];{\mathit{u}}_k^{*}=\left[\begin{array}{c}{F}_{\mathrm{Ac},k}\\ {\dot{z}}_{\mathrm{Ro},k}\end{array}\right]$$

The state-space matrices then have to be transformed accordingly. The state matrix is changed so that the current force becomes the delayed force in the next time step, and the input matrix is changed so that only the delayed input is used to calculate the new states.

$$\mathit{A}=\left[\begin{array}{cc}{\mathit{A}}_{ij}^{\left(\mathrm{d}\right)6pt}& {\mathit{B}}_{i1}^{\left(\mathrm{d}\right)6pt}\\ {\mathbf{0}}_{1,4}& 0\end{array}\right];\mathit{B}=\left[\begin{array}{cc}{\mathbf{0}}_{4,1}& {\mathit{B}}_{i2}^{\left(\mathrm{d}\right)6pt}\\ 1& 0\end{array}\right];\mathit{C}=\left[\begin{array}{cc}{\mathit{C}}_{ij}^{\left(\mathrm{d}\right)6pt}& {\mathit{D}}_{i1}^{\left(\mathrm{d}\right)6pt}\end{array}\right];\mathit{D}=\left[\begin{array}{cc}{\mathbf{0}}_{3,1}& {\mathit{D}}_{i2}^{\left(\mathrm{d}\right)6pt}\end{array}\right]$$

The discrete state-space model is then given as

$$\begin{array}{cc}\hfill {\mathit{x}}_{k+1}& =\mathit{A}{\mathit{x}}_k+\mathit{B}{\mathit{u}}_k\hfill \end{array}$$

(4a)

$$\begin{array}{cc}\hfill {\mathit{y}}_k& =\mathit{C}{\mathit{x}}_k+\mathit{D}{\mathit{u}}_k\hfill \end{array}$$

(4b)

with the state solution

$${\mathit{x}}_k={\mathit{A}}^k{\mathit{x}}_0+\sum _{i=0}^{k-1}{\mathit{A}}^{k-i-1}\mathit{B}{\mathit{u}}_k$$

(5)

and the output solution

$${\mathit{y}}_k=\mathit{C}\left[{\mathit{A}}^k{\mathit{x}}_0+\sum _{i=0}^{k-1}{\mathit{A}}^{k-i-1}\mathit{B}{\mathit{u}}_k\right]+\mathit{D}{\mathit{u}}_k$$

(6)

3.1.2. Road Model

The inverse Discrete Fourier Transform (DFT) method of Cebon and Newland [97] is used to generate a road surface profile, which is subsequently used as disturbance input. It is implemented according to Agostinacchio et al. [98], where the height profile is calculated as:

$$z_{\mathrm{Ro}}\left(x\right)=\sum _{i=0}^{\mathrm{N}}\sqrt{\Delta n}{2}^{{\mathrm{k}}_{\mathrm{r}}}·{10}^{-3}\left(\frac{n_0}{i\Delta n}\right)cos\left(2\pi i\Delta nx+{\phi }_i\right)$$

(7)

$$n_{\max }=\frac{1}{B}\Delta n=\frac{1}{L}\mathrm{N}=\frac{n_{\max }}{\Delta n}$$

In Equation (7), the total length of the road-profile is L (in $\mathrm{m}$), with a sampling interval of length B (in $\mathrm{m}$). This defines the maximum spatial sampling frequency resolution ($n_{\max }$), which produces $\mathrm{N}$ equally spaced spacial frequency domain points in the interval from 0 to $\Delta n$. Using the reference spatial frequency $n_0=0.1/\mathrm{m}$ and a set of random and uniformly distributed phase angles ${\phi }_i$ in the interval of $\left[0;2\pi \right]$, a road-profile can be generated in dependency of the parameter ${\mathrm{k}}_{\mathrm{r}}$, which corresponds to an ISO 8608 [99] road roughness classification according to

Table 1. Corresponding ISO 8608 [99] class ranges for the different values of ${\mathrm{k}}_{\mathrm{r}}$ [98].

${\mathbf{k}}_{\mathbf{r}}$	3	4	5	6	7	8	9
ISO class	A-B	B-C	C-D	D-E	E-F	F-G	G-H

.

Following generation by inverse DFT, the profile is filtered with a moving average filter of length $l_{\mathrm{cp}}$, which is the assumed contact patch length of the tire. Finally, the profile is converted to a height profile as a function of time, using the vehicle speed as in Equation (8):

$$z_{\mathrm{Ro}}\left(t\right)=z_{\mathrm{Ro}}(x/v_{x,\mathrm{veh}})$$

(8)

The single-obstacle (bump) is modelled according to Equation (9) as a half-sine wave [30] (pp. 345–346) with height $2\widehat{h}=0.1\mathrm{m}$ and length $L=3.8\mathrm{m}$.

$$f\left(x\right):=\left\{\begin{array}{cc}0\hfill & \mathrm{for}\mathrm{x}\le 0\mathrm{and}\mathrm{x}\ge \mathrm{L}\hfill \\ \widehat{h}\left(1-cos\left(\frac{2\pi }{L}x\right)\right)\hfill & \mathrm{for}\mathrm{x}\ge 0\mathrm{and}\mathrm{x}\le \mathrm{L}\hfill \end{array}\right.$$

(9)

3.1.3. Controller Design

In order to investigate the influence of actuator limitations on achievable ride comfort, it is necessary to exclude to the influence of controller performance. In general, an optimal control problem is formulated, which results in the best possible control input for the actuator. In the case of vertical vehicle dynamics control, it is necessary to incorporate constraints on the dynamic wheel-load and suspension travel (see Section 3.1.4). In optimal control, these constraints can only be considered indirectly via the cost function, and the optimization problem needs to be solved over the whole length of the simulation. MPC allows explicit consideration of these constraints and divides the problem into sub-problems of the size of the prediction horizon, which results in faster computation time, and was therefore used in this investigation.

Figure 3 shows a schematic overview of the simulation model with the MPC controller. The vehicle plant model from Section 3.1.1 is also used in the MPC controller. At time step k, the road disturbance $v_k$ is inputted to the plant, and the vector of road disturbances over the preview horizon ${\tilde{v}}_k$ is inputted to the MPC controller. The current state vector $x_k$ of the plant is also inputted to the MPC controller. This controller computes the optimal sequence of control inputs over the prediction horizon ${\tilde{u}}_k$ with respect to the cost function (Equation (17)), which is depending on the vector of outputvectors ${\tilde{y}}_k$ over the prediction horizon (PH). The first input of the computed optimal control sequence is then inputted to the plant for the next time step. To derive the cost function, it is necessary to split the input into actuator input u and road disturbance v, as shown in Equation (10).

$$\begin{array}{cc}\hfill {\mathit{x}}_{k+1}& =\mathit{A}{\mathit{x}}_k+{\mathit{b}}^{\left(\mathrm{u}\right)}{u}_k+{\mathit{b}}^{\left(\mathrm{v}\right)}{v}_k\hfill \end{array}$$

(10a)

$$\begin{array}{cc}\hfill {\mathit{y}}_k& =\mathit{C}{\mathit{x}}_k+\mathit{d}{v}_k\hfill \end{array}$$

(10b)

$${\mathit{b}}^{\left(\mathrm{u}\right)}={\mathit{B}}_{i1};{\mathit{b}}^{\left(\mathrm{v}\right)}={\mathit{B}}_{i2};\mathit{d}={\mathit{D}}_{i2};u_k=F_{\mathrm{Ac},k};v_k={\dot{z}}_{\mathrm{Ro},k}$$

The state solution becomes

$${\mathit{x}}_k={\mathit{A}}^k{\mathit{x}}_0+\sum _{i=0}^{k-1}{\mathit{A}}^{k-i-1}\left({\mathit{b}}^{\left(\mathrm{u}\right)}{u}_k+{\mathit{b}}^{\left(\mathrm{v}\right)}{v}_k\right)$$

(11)

and the output solution

$${\mathit{y}}_k=\mathit{C}\left[{\mathit{A}}^k{\mathit{x}}_0+\sum _{i=0}^{k-1}{\mathit{A}}^{k-i-1}\left({\mathit{b}}^{\left(\mathrm{u}\right)}{u}_k+{\mathit{b}}^{\left(\mathrm{v}\right)}{v}_k\right)\right]+\mathit{d}{v}_k$$

(12)

The solutions for the next p steps over the PH can be written in matrix form as

$$\begin{array}{c}\left[\begin{array}{c}{\mathit{y}}_{k+1}\\ {\mathit{y}}_{k+2}\\ \dots \\ {\mathit{y}}_{k+\mathrm{p}}\end{array}\right]=\left[\begin{array}{c}\mathit{CA}\\ \mathit{C}{\mathit{A}}^2\\ \dots \\ \mathit{C}{\mathit{A}}^{\mathrm{p}}\end{array}\right]{\mathit{x}}_k+\left[\begin{array}{cccc}\mathit{C}{\mathit{b}}^{\left(\mathrm{u}\right)}& 0& \dots & 0\\ \mathit{CA}{\mathit{b}}^{\left(\mathrm{u}\right)}& \mathit{C}{\mathit{b}}^{\left(\mathrm{u}\right)}& \dots & 0\\ \dots & \dots & \dots & \dots \\ \mathit{C}{\mathit{A}}^{\mathrm{p}-1}{\mathit{b}}^{\left(\mathrm{u}\right)}& \mathit{C}{\mathit{A}}^{\mathrm{p}-2}{\mathit{b}}^{\left(\mathrm{u}\right)}& \dots & \mathit{C}{\mathit{b}}^{\left(\mathrm{u}\right)}\end{array}\right]\left[\begin{array}{c}{u}_k\\ u_{k+1}\\ \dots \\ u_{k+\mathrm{p}-1}\end{array}\right]+\dots \hfill \\ \hfill \dots +\left[\begin{array}{ccccc}\mathit{C}{\mathit{b}}^{\left(\mathrm{v}\right)}& \mathit{d}& 0& \dots & 0\\ \mathit{CA}{\mathit{b}}^{\left(\mathrm{v}\right)}& \mathit{C}{\mathit{b}}^{\left(\mathrm{v}\right)}& \mathit{d}& \dots & 0\\ \dots & \dots & \dots & \dots & \dots \\ \mathit{C}{\mathit{A}}^{p-1}{\mathit{b}}^{\left(\mathrm{v}\right)}& \mathit{C}{\mathit{A}}^{p-2}{\mathit{b}}^{\left(\mathrm{v}\right)}& \mathit{C}{\mathit{A}}^{p-3}{\mathit{b}}^{\left(\mathrm{v}\right)}& \dots & \mathit{d}\end{array}\right]\left[\begin{array}{c}{v}_k\\ v_{k+1}\\ \dots \\ v_{k+\mathrm{p}}\end{array}\right]\end{array}$$

(13)

and by defining

$$\tilde{\mathit{y}}=\left[\begin{array}{c}{\mathit{y}}_{k+1}\\ {\mathit{y}}_{k+2}\\ \dots \\ {\mathit{y}}_{k+\mathrm{p}}\end{array}\right];\tilde{\mathit{u}}=\left[\begin{array}{c}{u}_k\\ u_{k+1}\\ \dots \\ u_{k+\mathrm{p}-1}\end{array}\right];\tilde{\mathit{v}}=\left[\begin{array}{c}{v}_k\\ v_{k+1}\\ \dots \\ v_{k+\mathrm{p}}\end{array}\right];\mathbf{\Omega }=\left[\begin{array}{c}\mathit{CA}\\ \mathit{C}{\mathit{A}}^2\\ \dots \\ \mathit{C}{\mathit{A}}^{\mathrm{p}}\end{array}\right];$$

$${\mathbf{\Gamma }}_{\mathrm{u}}=\left[\begin{array}{cccc}\mathit{C}{\mathit{b}}^{\left(\mathrm{u}\right)}& 0& \dots & 0\\ \mathit{CA}{\mathit{b}}^{\left(\mathrm{u}\right)}& \mathit{C}{\mathit{b}}^{\left(\mathrm{u}\right)}& \dots & 0\\ \dots & \dots & \dots & \dots \\ \mathit{C}{\mathit{A}}^{\mathrm{p}-1}{\mathit{b}}^{\left(\mathrm{u}\right)}& \mathit{C}{\mathit{A}}^{\mathrm{p}-2}{\mathit{b}}^{\left(\mathrm{u}\right)}& \dots & \mathit{C}{\mathit{b}}^{\left(\mathrm{u}\right)}\end{array}\right];{\mathbf{\Gamma }}_{\mathrm{v}}=\left[\begin{array}{ccccc}\mathit{C}{\mathit{b}}^{\left(\mathrm{v}\right)}& \mathit{d}& 0& \dots & 0\\ \mathit{C}\mathit{A}{\mathit{b}}^{\left(\mathrm{v}\right)}& \mathit{C}{\mathit{b}}^{\left(\mathrm{v}\right)}& \mathit{d}& \dots & 0\\ \dots & \dots & \dots & \dots & \dots \\ \mathit{C}{\mathit{A}}^{p-1}{\mathit{b}}^{\left(\mathrm{v}\right)}& \mathit{C}{\mathit{A}}^{p-2}{\mathit{b}}^{\left(\mathrm{v}\right)}& \mathit{C}{\mathit{A}}^{p-3}{\mathit{b}}^{\left(\mathrm{v}\right)}& \dots & \mathit{d}\end{array}\right];$$

Equation (13) becomes

$${\tilde{\mathit{y}}}_k=\mathbf{\Omega }{\mathit{x}}_k+{\mathbf{\Gamma }}_{\mathrm{u}}{\tilde{\mathit{u}}}_k+{\mathbf{\Gamma }}_{\mathrm{v}}{\tilde{\mathit{v}}}_k$$

(14)

Using Equation (14) and introducing the weighting matrix $\mathbf{W}$, which is used for each time step and gives a vector of matrices $\tilde{\mathbf{W}}$, the vector of slack variables $\mathit{ϵ}$ and the constraint violation penalty $\rho$, the cost function for the optimal control problem over the prediction horizon can defined as:

$$J_k\left({\tilde{\mathit{u}}}_k,\mathbf{ϵ}\right)={\left(\mathbf{\Omega }{\mathit{x}}_k+{\Gamma }_{\mathrm{u}}{\tilde{\mathit{u}}}_k+{\mathbf{\Gamma }}_{\mathrm{v}}{\tilde{\mathit{v}}}_k\right)}^{\mathrm{T}}\tilde{\mathbf{W}}\left(\mathbf{\Omega }{\mathit{x}}_k+{\Gamma }_{\mathrm{u}}{\tilde{\mathit{u}}}_k+{\mathbf{\Gamma }}_{\mathrm{v}}{\tilde{\mathit{v}}}_k\right)+{\mathsf{ϵ}}^{\mathrm{T}}\mathsf{\rho }\mathsf{ϵ}$$

(15)

and rewritten as

$$\begin{array}{c}{J}_k\left({\tilde{\mathit{u}}}_k,\mathsf{ϵ}\right)=\frac{1}{2}{\tilde{\mathit{u}}}_k^{\mathrm{T}}{\mathbf{\Gamma }}_{\mathrm{u}}^{\mathrm{T}}\tilde{\mathbf{W}}{\mathbf{\Gamma }}_{\mathrm{u}}{\tilde{\mathit{u}}}_k+\left({\mathit{x}}_k^{\mathrm{T}}{\tilde{\mathbf{\Omega }}}^{\mathrm{T}}\tilde{\mathbf{W}}{\mathbf{\Gamma }}_{\mathrm{u}}+{\tilde{\mathit{v}}}_k^{\mathrm{T}}{\mathbf{\Gamma }}_{\mathrm{v}}^{\mathrm{T}}\tilde{\mathbf{W}}{\mathbf{\Gamma }}_{\mathrm{u}}\right){\tilde{\mathit{u}}}_k+\frac{1}{2}{\mathsf{ϵ}}^{\mathrm{T}}\mathsf{\rho }\mathsf{ϵ}+\dots \hfill \\ \dots +\frac{1}{2}{\tilde{\mathit{v}}}_k^{\mathrm{T}}{\mathbf{\Gamma }}_{\mathrm{v}}^{\mathrm{T}}\tilde{\mathbf{W}}{\mathbf{\Gamma }}_{\mathrm{v}}{\tilde{\mathit{v}}}_k+{\mathit{x}}_k^{\mathrm{T}}{\mathbf{\Omega }}^{\mathbf{T}}{\tilde{\mathbf{T}}}^2{\mathbf{\Gamma }}_{\mathrm{v}}{\tilde{\mathit{v}}}_k+{\mathit{x}}_k^{\mathrm{T}}{\mathbf{\Omega }}^{\mathrm{T}}\tilde{\mathbf{W}}\mathbf{\Omega }{\mathit{x}}_k\hfill \end{array}$$

(16)

The cost function is only a function of the input and the slack variable, and all terms that do not include these are constant and can be ignored. The final cost function is shown in Equation (17).

$$J_k\left({\tilde{\mathit{u}}}_k,\mathsf{ϵ}\right)=\frac{1}{2}{\tilde{\mathit{u}}}_k^{\mathrm{T}}{\mathbf{\Gamma }}_{\mathrm{u}}^{\mathrm{T}}\tilde{\mathbf{W}}{\mathbf{\Gamma }}_{\mathrm{u}}{\tilde{\mathit{u}}}_k+\left({\mathit{x}}_k^{\mathrm{T}}{\tilde{\mathbf{\Omega }}}^{\mathrm{T}}\tilde{\mathbf{W}}{\mathbf{\Gamma }}_{\mathrm{u}}+{\tilde{\mathit{v}}}_k^{\mathrm{T}}{\mathbf{\Gamma }}_{\mathrm{v}}^{\mathrm{T}}\tilde{\mathbf{W}}{\mathbf{\Gamma }}_{\mathrm{u}}\right){\tilde{\mathit{u}}}_k+\frac{1}{2}{\mathsf{ϵ}}^{\mathrm{T}}\mathsf{\rho }\mathsf{ϵ}$$

(17)

The QP KWIK solver [100], which is based on the dual algorithm of Goldfarb and Idnani [101], is used to solve the quadratic problem in order to obtain the optimal control inputs for each time step. This algorithm is pre-implemented in the MATLAB© Model Predictive Control Toolbox, but can also be implemented in non-proprietary software such as GNU Octave or other programming languages. Another option is to use a solver from the YALMIP toolset [102] if the Model Predictive Control Toolbox and Optimization Toolbox are not available in MATLAB©.

3.1.4. Actuator Limitations, System Constraints and Weights

The optimization problem formulated in Equation (17) needs to be solved with consideration for the constraints of the control input (actuator force and rate) and output constraints on the dynamic wheel-load and suspension travel. This is necessary to ensure that the problem solution stays within the valid range of the simulation model. The optimization problem is therefore subject to the following linear inequality input constraints

$$\begin{array}{cc}\hfill F_{\mathrm{Ac},\min }\le & F_{\mathrm{Ac},k}\le F_{\mathrm{Ac},\max }\hfill \end{array}$$

(18a)

$$\begin{array}{cc}\hfill {\dot{F}}_{\mathrm{Ac},\min }\le & {\dot{F}}_{\mathrm{Ac},k}\le {\dot{F}}_{\mathrm{Ac},\max }\hfill \end{array}$$

(18b)

and output constraints

$$\begin{array}{cc}\hfill F_{\mathrm{dyn},\min }-{\mathsf{ϵ}}_k\le & F_{\mathrm{dyn},k}\le F_{\mathrm{dyn},\max }+{\mathsf{ϵ}}_k\hfill \end{array}$$

(19a)

$$\begin{array}{cc}\hfill z_{\mathrm{sus},\min }-{\mathsf{ϵ}}_k\le & z_{\mathrm{sus},k}\le z_{\mathrm{sus},\max }+{\mathsf{ϵ}}_k\hfill \end{array}$$

(19b)

$F_{\mathrm{Ac},k}$ is the actuator force (control input) at time step $\mathrm{k}$, $F_{\mathrm{dyn},k}$ and $z_{\mathrm{sus},k}$ are the discrete output values of the dynamic wheel-load and the suspension-deflection. The full output vector including the body acceleration $a_{\mathrm{Bo},k}$ is shown in Appendix A. The values of the constraints on the actuator force are shown in Section 3.2. For all simulations, the output constraints were set according to

Table 2. Output constraints.

Parameter	Value	Unit
$F_{\mathrm{dyn},\min }$	−5400	$\mathrm{N}$
$F_{\mathrm{dyn},\max }$	inf	$\mathrm{N}$
$z_{\mathrm{sus},\min }$	−0.05	$\mathrm{m}$
$z_{\mathrm{sus},\max }$	0.05	$\mathrm{m}$

.

The constraint on the dynamic wheel-load ensures that the actual wheel-load (static + dynamic) is always positive and lift off does not occur. Lift off is non-linear and can therefore not be incorporated in this model. The maximum value of the dynamic wheel load in pressure is unconstrained. The suspension travel is limited to $10\mathrm{c}\mathrm{m}$ to ensure that it stays within the linear range. This value is set conservatively and may be increased depending on the vehicle. Unlike the input constrains, the output constraints are soft constraints, which allows for a minor violation of the constraint. It is done to ensure feasibility of the optimization problem. The validity of the results was checked in post-processing, see Section 3.3.

Figure 4 shows an example of the implementation of the actuator constraints in the controller. The dashed red line shows the desired force of a non-restricted controller without delay. The blue line shows the force of a controller, in which the maximum force is limited to $2000\mathrm{N}$ and the maximum rate to $10\mathrm{k}\mathrm{N}\mathrm{s}$⁻¹. The restricted force also contains the delay of one time step from controller output to actuator input. The dotted orange line shows the highest possible change of the actuator force according to the rate constraint. At around $0.04\mathrm{s}$, the unrestricted force changes with a step to $2500\mathrm{N}$. The restricted force starts to rise accordingly at the maximum rate of change, delayed by one time step, and saturates at the maximum force limit. At around $0.45\mathrm{s}$, the desired force descends below the maximum force limit at a rate smaller than the maximum rate restriction. The restricted force follows the unrestricted, delayed by one time step. Finally, the unrestricted force drops from $1000\mathrm{N}$ to $0\mathrm{N}$, and the restricted force follows in compliance with the rate limitation.

3.2. Test Scenarios

Different test scenarios were used in this research, defined by the road-profile, the vehicle speed and the optimization objective (

Table 3. Simulation scenarios.

Road-Profile	Speed	Objective	Time
Bump	$36\mathrm{k}\mathrm{m}\mathrm{h}$−1	${\ddot{z}}_{\mathrm{Bo}}$	$5\mathrm{s}$
ISO A-B	$45\mathrm{k}\mathrm{m}\mathrm{h}$−1 $90\mathrm{k}\mathrm{m}\mathrm{h}$−1 $180\mathrm{k}\mathrm{m}\mathrm{h}$−1	${\ddot{z}}_{\mathrm{Bo}}$	$40\mathrm{s}$ $20\mathrm{s}$ $10\mathrm{s}$
ISO B-C	$90\mathrm{k}\mathrm{m}\mathrm{h}$−1	${\ddot{z}}_{\mathrm{Bo}}$ $F_{\mathrm{dyn}}$	$20\mathrm{s}$

). A simulation was performed on the single-obstacle profile at $36\mathrm{k}\mathrm{m}\mathrm{h}$⁻¹, with the objective of minimizing the vehicle body acceleration for comfort. The single-obstacle profile is shown in Figure 5. The figure only shows the section of the bump; the full profile is extended as a flat road before and after the bump.

On the stochastic road-profile A-B, three simulations were carried out at speeds of $45\mathrm{k}\mathrm{m}\mathrm{h}$⁻¹, $90\mathrm{k}\mathrm{m}\mathrm{h}$⁻¹ and $180\mathrm{k}\mathrm{m}\mathrm{h}$⁻¹. In all three cases, the objective was to minimize the vehicle body acceleration. On the B-C profile, two simulations were performed at $90\mathrm{k}\mathrm{m}\mathrm{h}$⁻¹; in one case, the objective was to minimize the body acceleration, while in the other case, the objective was to minimize the dynamic wheel-load variation. Figure 6a shows the stochastic road-profiles. The spatial power spectral density of the stochastic road-profiles can be seen in Figure 6b.

All simulations were carried out at a sampling frequency of $100\mathrm{Hz}$ and a preview horizon of $2\mathrm{s}$. The control horizon was set to its maximum of $1.99\mathrm{s}$, which is one sample less than the preview horizon. A summary of these settings is given in

Table 4. Simulation settings.

Name	Parameter	Value	Unit
Sampling time	$t_{\mathrm{s}}$	0.01	$\mathrm{s}$
Preview horizon	$t_{\mathrm{p}}$	2	$\mathrm{s}$
Control horizon	$t_{\mathrm{c}}$	1.99	$\mathrm{s}$

. The quarter-car model is parameterized according to

Table 5. Model parameters.

Name	Parameter	Value	Unit
Sprung mass	$m_{\mathrm{Bo}}$	485	$\mathrm{k}\mathrm{g}$
Unsprung mass	$m_{\mathrm{Wh}}$	65	$\mathrm{k}\mathrm{g}$
Suspension stiffness	$c_{\mathrm{Bo}}$	24,000	$\mathrm{N}\mathrm{m}$−1
Tire stiffness	$c_{\mathrm{Wh}}$	360,000	$\mathrm{N}\mathrm{m}$−1
Suspension damping	$d_{\mathrm{Bo}}$	1500	$\mathrm{N}\mathrm{s}\mathrm{m}$−1
Tire damping	$d_{\mathrm{Wh}}$	80	$\mathrm{N}\mathrm{s}\mathrm{m}$−1
Contact patch length	$l_{\mathrm{Cp}}$	0.1	$\mathrm{m}$
Actuator delay	$t_{\mathrm{Ac}}$	0.01	$\mathrm{s}$

. These parameters correspond to a vehicle with a total mass of $2200\mathrm{k}\mathrm{g}$, a body eigenfrequency of $1.1\mathrm{Hz}$, a wheel eigenfrequency of $12.2\mathrm{Hz}$ and $22$% critical damping. These are exemplary values for a modern, premium class luxury vehicle, for example a BMW 7-Series or a Mercedes S-Class. The formulae for calculating the full vehicle parameters can be found in Appendix B.

The actuator limitations of maximum force and maximum rate of change were varied from $0\mathrm{N}$ to 10,000 N and $0\mathrm{N}\mathrm{s}$⁻¹ to 10,000 N s¹ respectively. A passive vehicle and a vehicle with an MPC controller without input constraints were also simulated according to the scenario. The passive vehicle was the baseline and the ideal MPC vehicle was the reference for what could be achieved in the best case.

Table 6. Actuator limitations (input constraints) for the single-obstacle scenario.

${\mathit{F}}_{\mathbf{Ac},\mathbf{lim}}$in kN	0.25	0.5	0.7	1	1.25	1.5	1.75	2	2.25	2.5	3	3.5	4	5	10
${\dot{F}}_{\mathrm{Ac},\lim }$ in $\mathrm{k}\mathrm{N}\mathrm{s}$−1	1	2	3	4	5	7.5	10	12.5	15	17.5	20	22.5	25	50	100

shows the different values of the input constraints in the single-obstacle scenario. They were defined in an iterative process to achieve good resolution in the areas of interest; these are shown in Section 3.3. All combinations of force and rate constraints were simulated.

Table 7. Actuator limitations (input constraints) for the stochastic road-profile scenario.

${\mathit{F}}_{\mathbf{Ac},\mathbf{lim}}$ in kN	0.1	0.2	0.3	0.4	0.5	0.75	1	1.25	1.5	1.75	2	2.5	3	4	5	10
${\dot{F}}_{\mathrm{Ac},\lim }$ in $\mathrm{k}\mathrm{N}\mathrm{s}$−1	1	5	10	15	20	25	30	35	40	45	50	60	70	80	90	100

shows the different input constraints for the stochastic road-profile scenarios. All possible combinations were simulated. Compared with the limitations for the bump-profile, the values were chosen to have better resolution at lower maximum force limits and higher maximum rate of change limits. Again, this was determined in preliminary experiments and will be explained in Section 3.3.

Table 8. Weighting factors for the two optimization objectives.

	${\mathbf{w}}_{{\ddot{\mathit{z}}}_{\mathbf{Bo}}}$	${\mathbf{w}}_{{\mathit{z}}_{\mathbf{sus}}}$	${\mathbf{w}}_{{\mathit{F}}_{\mathbf{dyn}}}$	${\mathbf{w}}_{{\mathit{F}}_{\mathbf{act}}}$	$\mathsf{\rho }$
$\mathbf{min}\left({\ddot{\mathit{z}}}_{\mathbf{Bo}}\right)$	10×${\left(1\right)}^2$	$\mathbf{0}.\mathbf{1}$×${\left(100\right)}^2$	$\mathbf{0}.\mathbf{1}$×${\left({10}^{-3}\right)}^2$	-	$\mathbf{1}$×${\mathbf{10}}^{\mathbf{10}}$
$\min \left(F_{\mathrm{dyn}}\right)$	$\mathbf{0}.\mathbf{1}$×${\left(1\right)}^2$	$\mathbf{1}$×${\left(100\right)}^2$	$\mathbf{10}$×${\left({10}^{-3}\right)}^2$	$\mathbf{4}.\mathbf{1}$×${\left({10}^{-3}\right)}^2$	$\mathbf{1}$×${\mathbf{10}}^{\mathbf{10}}$

shows the weighting factors, which were used for the two different optimization objectives. Each one consists of the relative weighting factor in bold font and a normalization term in brackets. Normalization is used to even out the order of the different variables, which is necessary to apply relative weighting. The weighting factor for the slack variable ($\mathsf{\rho }$) was always set to 1 × 10${}^{10}$. For optimizing minimum body acceleration, the relative weighting for body acceleration ${\mathrm{w}}_{{\ddot{z}}_{\mathrm{Bo}}}$ was set to 10, while relative weighting factors for suspension-deflection ${\mathrm{w}}_{z_{\mathrm{sus}}}$ and dynamic wheel-load ${\mathrm{w}}_{F_{\mathrm{dyn}}}$ were set to $0.1$. These three weighting factors form the trace of the diagonal weighting matrix $\mathbf{W}$, which is concatenated to the vector of weighting matrices $\tilde{\mathbf{W}}$ in the cost function (Equation (17)). Suspension-deflection and dynamic wheel-load weighting were set to non-zero, because otherwise the controller would not consider these in the optimization, and run into the respective constraint, resulting in poorer minimization of body acceleration. For optimization with respect to the dynamic wheel-load, the weighting factor ${\mathrm{w}}_{z_{\mathrm{sus}}}$ was increased to 1, because the controller was more prone to hitting the suspension-deflection constraint. An additional term to consider the actuator force was also added to the cost function for stability reasons. Without this term, the controller would become unstable with this objective. It was added by defining another output equation for the input force and adding the weighting factor ${\mathrm{w}}_{F_{\mathrm{act}}}$ to the weighting Matrix $\mathbf{W}$. It avoids the input force increasing beyond all limits and become unstable. The value was determined by performing a sensitivity analysis with respect to the minimal achievable dynamic wheel-load variation.

3.3. Data Processing and Analysis

A number of steps were carried out to obtain the results presented in this research paper. The process is summarized in Figure 7. As shown in Section 3.2, each scenario was simulated with a variety of actuator constraints and therefore consists of multiple simulations. A single simulation outputs the trajectories of the body acceleration ${\ddot{z}}_{\mathrm{Bo}}$, suspension travel $s_{\mathrm{sus}}$ and dynamic wheel-load $F_{\mathrm{dyn}}$, along with information about the MPC solution in each time step. This information contains the slack variable, the number of iterations of the solver, the computation time and whether an optimal solution could be obtained.

3.3.1. Evaluation of the Single Simulation Results

In this investigation, the root mean square (RMS) value is used to quantify the outputs of a simulation (${\ddot{z}}_{\mathrm{Bo}}$, $s_{\mathrm{sus}}$ and $F_{\mathrm{dyn}}$). The ISO 2631 cannot be used here, for two reasons. First, the model only outputs the vertical body acceleration, so, the directional weighting of the ISO 2631 cannot be applied and the comfort value would not show a benefit compared to an analyses of the RMS value. Second, the frequency weighting of the ISO 2631 is a non-linear operation and therefore cannot be directly implemented in the MPC controller. In the case of the quarter-car model, the reduction of the body acceleration directly correlates to the reduction of the ISO 2631 comfort value, as can be seen in Figure A1 in Appendix C. The plausibility and validity of the results were also checked, as were violated constraints and reasonable trajectories of the state and output variables. Minor violations of the dynamic wheel-load and suspension travel limits were considered admissible because the suspension travel limit was set conservatively to $\pm 5\mathrm{c}\mathrm{m}$. In reality, the wheel could lift off for a short time and the inertia of the wheel could be used to generate a force larger than the static wheel-load. The maximum and mean values of the violations are shown in Appendix E, along with the number of violation time steps.

3.3.2. Processing of the Scenario Data

In order to enable a condensed overview of the results in each scenario, the results of the single simulations are subjected to further processing. First, the RMS results are combined into a surface as a function of the maximum force and the maximum rate constraint. This surface is displayed as surface plots for ${\ddot{z}}_{\mathrm{Bo}}$, $s_{\mathrm{sus}}$ and $F_{\mathrm{dyn}}$ and as contour plots for ${\ddot{z}}_{\mathrm{Bo}}$ or $F_{\mathrm{dyn}}$ respectively. Additionally, the gradient of the surface of ${\ddot{z}}_{\mathrm{Bo}}$ or $F_{\mathrm{dyn}}$ is calculated and used to display the trajectory of gradient-descent from the origin. This trajectory line can be used to analyze the optimal combination of force and rate limitations. In other words, it makes it possible to see how much rate of change an actuator should be able to supply for a given target of maximum force, and vice versa.

4. Results

The following section presents and describes the results for the different simulation scenarios. These include RMS surface plots of the body acceleration ${\ddot{z}}_{\mathrm{Bo}}$, the dynamic wheel-load $F_{\mathrm{dyn}}$, the suspension-deflection $s_{\mathrm{sus}}$ and the actuator force $F_{\mathrm{act}}$, as well as a contour plot of the body acceleration ${\ddot{z}}_{\mathrm{Bo}}$ or dynamic wheel-load $F_{\mathrm{dyn}}$ with the gradient-descent line. FFT gain plots of the body acceleration ${\ddot{z}}_{\mathrm{Bo}}$ and the dynamic wheel-load $F_{\mathrm{dyn}}$ are given in Appendix D and Appendix E along with heat-map plots showing possible constraint violations. The results of each scenario in the passive simulation and the unrestricted MPC actuator simulation are summarized in

Table 9. RMS values for the passive and unrestricted MPC simulations.

	Bump 36 km h${}^{-1}$		AB 90 km h${}^{-1}$		AB 45 km h${}^{-1}$		AB 180 km h${}^{-1}$		BC 90 km h${}^{-1}$
	PAS	MPC	PAS	MPC	PAS	MPC	PAS	MPC	PAS	MPC
${\ddot{z}}_{\mathrm{Bo}},\mathrm{rms}$	$1.29\mathrm{m}\mathrm{s}$−2	$0.05\mathrm{m}\mathrm{s}$−2	$0.71\mathrm{m}\mathrm{s}$−2	$0.01\mathrm{m}\mathrm{s}$−2	$0.38\mathrm{m}\mathrm{s}$−2	$0.004\mathrm{m}\mathrm{s}$−2	$1.54\mathrm{m}\mathrm{s}$−2	$0.48\mathrm{m}\mathrm{s}$−2	$1.51\mathrm{m}\mathrm{s}$−2	$0.54/1.81\mathrm{m}\mathrm{s}$−2
$z_{\mathrm{sus}},\mathrm{rms}$	$2.2\mathrm{c}\mathrm{m}$	$2.0\mathrm{c}\mathrm{m}$	$0.7\mathrm{c}\mathrm{m}$	$1.0\mathrm{c}\mathrm{m}$	$0.4\mathrm{c}\mathrm{m}$	$0.6\mathrm{c}\mathrm{m}$	$1.8\mathrm{c}\mathrm{m}$	$1.9\mathrm{c}\mathrm{m}$	$1.6\mathrm{c}\mathrm{m}$	$1.9/0.5\mathrm{c}\mathrm{m}$
$F_{\mathrm{dyn}},\mathrm{rms}$	$611\mathrm{N}$	$211\mathrm{N}$	$987\mathrm{N}$	$1840\mathrm{N}$	$503\mathrm{N}$	$1160\mathrm{N}$	$1960\mathrm{N}$	$2626\mathrm{N}$	$2024\mathrm{N}$	$2762/1412\mathrm{N}$
$F_{\mathrm{act}},\mathrm{rms}$	-	$528\mathrm{N}$	-	$580\mathrm{N}$	-	$362\mathrm{N}$	-	$846\mathrm{N}$	-	$884/737\mathrm{N}$

.

In the single-obstacle bump scenario, the MPC controller can reduce body acceleration from $1.29\mathrm{m}{\mathrm{s}}^{-2}$ to $0.05\mathrm{m}{\mathrm{s}}^{-2}$ and dynamic wheel-load variation from $611\mathrm{N}$ to $211\mathrm{N}$. Although the objective was to minimize body acceleration, in the case of the single-obstacle with a smooth surface, the reduction of body acceleration also entailed a reduction in dynamic wheel-load variation, because the bump mainly excited the body eigenmode. Suspension-deflection RMS was slightly reduced from $2.2\mathrm{c}\mathrm{m}$ to $2.0\mathrm{c}\mathrm{m}$, and the RMS actuator force was $528\mathrm{N}$. On the AB profile at $90\mathrm{km}{\mathrm{h}}^{-1}$, the body acceleration RMS was reduced from $0.71\mathrm{m}{\mathrm{s}}^{-2}$ to $0.01\mathrm{m}{\mathrm{s}}^{-2}$, at the cost of an increase in dynamic wheel-load variation from $987\mathrm{N}$ to $1840\mathrm{N}$. Suspension travel is also slightly increased from $0.7\mathrm{c}\mathrm{m}$ to $1.0\mathrm{c}\mathrm{m}$ and the RMS actuator force was $580\mathrm{N}$. At $45\mathrm{km}{\mathrm{h}}^{-1}$ on this profile, the body acceleration was reduced from $0.38\mathrm{m}{\mathrm{s}}^{-2}$ to $0.004\mathrm{m}{\mathrm{s}}^{-2}$, with an increase from $503\mathrm{N}$ to $1160\mathrm{N}$ in dynamic wheel-load variation. Suspension travel is increased from $0.4\mathrm{c}\mathrm{m}$ to $0.6\mathrm{c}\mathrm{m}$ with an RMS actuator force of $362\mathrm{N}$. At $180\mathrm{km}{\mathrm{h}}^{-1}$, the body acceleration decreased from $1.54\mathrm{m}{\mathrm{s}}^{-2}$ to $0.48\mathrm{m}{\mathrm{s}}^{-2}$. Due to the constraints on dynamic wheel-load and suspension-deflection, the MPC controller was not able to decrease the body acceleration any further. The dynamic wheel-load variation rose from $1960\mathrm{N}$ to $2626\mathrm{N}$, suspension-deflection RMS showed a small increase from $1.8\mathrm{c}\mathrm{m}$ to $1.9\mathrm{c}\mathrm{m}$, and the actuator RMS force was $846\mathrm{N}$. On the BC profile where the objective was to minimize body acceleration, it was reduced from $1.51\mathrm{m}{\mathrm{s}}^{-2}$ to $0.54\mathrm{m}{\mathrm{s}}^{-2}$. Again, the controller was not able to reduce the body acceleration RMS any further due to the output constraints. Dynamic wheel-load and suspension-deflection increased from $2024\mathrm{N}$ to $2762\mathrm{N}$ and from $1.6\mathrm{c}\mathrm{m}$ to $1.9\mathrm{c}\mathrm{m}$ respectively. The RMS actuator force was $884\mathrm{N}$. When the objective was to minimize dynamic wheel-load variation, its RMS value was decreased to $1412\mathrm{N}$ at the cost of an increase in body acceleration RMS to $1.81\mathrm{m}{\mathrm{s}}^{-2}$. Here, the suspension-deflection RMS value fell to $0.5\mathrm{c}\mathrm{m}$ and the actuator RMS force was $737\mathrm{N}$.

4.1. Single-Obstacle Bump

In the single-obstacle scenario, the vehicle drove over a bump (see Section 3.2) at a speed of $36\mathrm{km}{\mathrm{h}}^{-1}$. This resembles a situation in a restricted traffic zone in a residential area. Figure 8 shows the RMS surface plots for this scenario. In each plot, the lower left axis denotes the maximum actuator force $F_{\mathrm{Ac}},\max$ in $\mathrm{k}\mathrm{N}$, and the lower right axis represents the maximum actuator rate of change ${\dot{F}}_{\mathrm{Ac}},\max$ in $\mathrm{kN}{\mathrm{s}}^{-1}$. The z-axis shows the respective RMS values. In each plot, the intersection of the surface with the xz and xy plane corresponds to the RMS value of the passive vehicle. The transparent grey plane shows the result for the MPC controller with a non-restricted actuator. In this scenario, the reduction in body acceleration also led to a reduction in dynamic wheel-load variation, because the bump mainly excited the vehicle’s body eigenmode. Both decrease with increasing maximum force and rate limit. The suspension-deflection RMS decreases for low maximum forces, and peaks for low maximum force and high rate of change limits. The RMS actuator force saturates around a maximum force limit of $1.5\mathrm{k}\mathrm{N}$ and a maximum rate limit of $10\mathrm{kN}{\mathrm{s}}^{-1}$.

Figure 9 shows the contour plot of the RMS body acceleration with isolines from $0.9\mathrm{m}{\mathrm{s}}^{-2}$ to $0.05\mathrm{m}{\mathrm{s}}^{-2}$ in $0.2\mathrm{m}{\mathrm{s}}^{-2}$ steps. The closer the isolines are to each other, the steeper the gradient of the surface in this area. It can be seen that the RMS body acceleration depends greatly on the maximum rate limit and that a relatively small maximum force is sufficient, provided the actuator can supply a high enough rate of change. The largest decrease in body acceleration can be seen up to a maximum force limit of $1500\mathrm{N}$, and no further decrease in body acceleration can be achieved with maximum forces larger than $2250\mathrm{N}$. This can be deduced by following the gradient descent line, which points completely in the direction of the increased rate of change limit, after $2250\mathrm{N}$ at $26\mathrm{kN}{\mathrm{s}}^{-1}$. The perpendicular shape of the isolines also shows that it is not possible to compensate a low maximum rate of change with a high maximum force limit and vice versa. For a target value of $0.1\mathrm{m}{\mathrm{s}}^{-2}$ the ideal actuator (minimum maximum force and maximum rate) for the chosen output constraints would have a maximum force of around $2000\mathrm{N}$ and a maximum rate of change of around $21.5\mathrm{kN}{\mathrm{s}}^{-1}$.

The FFT plots of body acceleration (Figure A2) and dynamic wheel-load (Figure A3), which display the amplitude of the signals in the frequency domain for different actuator configurations, can be found in Appendix D. Here, it can be seen that the vehicle is mainly excited in the body eigenfrequency range and that the controller lowers the peak for increasing actuator potential. No violation of the dynamic wheel-load constraint occurred in this scenario. Suspension-deflection constraint violations larger than $0.5\mathrm{m}\mathrm{m}$ only occurred in the simulations with a maximum force limit of $250\mathrm{N}$ and $500\mathrm{N}$, with rate constraints of up to $3\mathrm{kN}{\mathrm{s}}^{-1}$.

4.2. Road-Profile AB

Three different scenarios were investigated on the AB road-profile (see Section 3.2). The objective was to minimize body acceleration while travelling at $90\mathrm{km}{\mathrm{h}}^{-1}$, $45\mathrm{km}{\mathrm{h}}^{-1}$ and $180\mathrm{km}{\mathrm{h}}^{-1}$. The $90\mathrm{km}{\mathrm{h}}^{-1}$ scenario corresponds to travelling on a good country road, the $45\mathrm{km}{\mathrm{h}}^{-1}$ scenario to inner city driving and the $180\mathrm{km}{\mathrm{h}}^{-1}$ scenario to driving on a German motorway.

4.2.1. Good Country Road at $90\mathrm{km}{\mathrm{h}}^{-1}$

In the first scenario with a speed of $90\mathrm{km}{\mathrm{h}}^{-1}$, the body acceleration RMS was reduced to nearly zero at the cost of an almost doubled dynamic wheel-load variation, as shown in Figure 10. Suspension travel RMS shows a moderate increase with a maximum for low maximum rate limits and high maximum force limits. Suspension-deflection is smallest for maximum force limits smaller than $1000\mathrm{N}$. Actuator force RMS increases with increasing maximum force and rate limits up to the value of the unrestricted actuator.

Figure 11 shows the contour plot of the body acceleration RMS values, with isolines from $0.5\mathrm{m}{\mathrm{s}}^{-2}$ to $0.1\mathrm{m}{\mathrm{s}}^{-2}$ in $0.1\mathrm{m}{\mathrm{s}}^{-2}$ steps. Again, the influence of the actuator rate of change is dominant, and the perpendicular contour lines indicate that an increase in one limit cannot compensate for the restriction due to the other one. The gradient is steepest for actuator maximum force limits smaller than $1000\mathrm{N}$. Above this limit, the actuator rate of change is predominant. For a target value of $0.1\mathrm{m}{\mathrm{s}}^{-2}$ RMS body acceleration, the ideal actuator for the chosen output constraints would need to have a maximum force of roughly $1400\mathrm{N}$ and a maximum rate of change of about $62\mathrm{kN}{\mathrm{s}}^{-1}$ in this scenario. Above a maximum force limit of $2000\mathrm{N}$ and a maximum rate limit of $100\mathrm{kN}{\mathrm{s}}^{-1}$, no significant decrease in body acceleration RMS value could be achieved.

The FFT plots of body acceleration (Figure A4) and dynamic wheel-load for selected limits can be found in Appendix D. They show how a decrease in body acceleration leads to an increase in dynamic wheel-load variation around the wheel eigenfrequency. Even though the MPC controller reduces the dynamic wheel-load variation at the body eigenmode and can hold the variation at the level of the passive suspension between $3\mathrm{Hz}$ and $8\mathrm{Hz}$ and above $14\mathrm{Hz}$, it cannot avoid the increase around the wheel eigenfrequency, where the objective is to minimize body acceleration.

Contrary to the single-obstacle scenario, there was no violation of the suspension-deflection limit. Instead, for maximum force limits higher than $1000\mathrm{N}$ combined with rate limits above $70\mathrm{kN}{\mathrm{s}}^{-1}$, a small violation of the dynamic wheel-load constraint appeared (Appendix E, Figure A15). The worst case is a violation in three time steps with a mean value of $589\mathrm{N}$, which is around $10\%$ of the static wheel-load for $0.03\mathrm{s}$.

4.2.2. Inner City at $45\mathrm{km}{\mathrm{h}}^{-1}$

The RMS performance indicators in Figure 12 show qualitatively the same characteristics in the inner city scenario as for the good country road scenario, but at a quantitatively lower level. The body acceleration RMS value saturates early at the level of the unrestricted actuator. This is also the case for suspension-deflection, dynamic wheel-load and actuator force RMS values.

Qualitatively, the contour plot of the body acceleration in this scenario is also comparable to the good country road scenario, as shown in Figure 13. However, the gradient descent line is shifted slightly downwards, and for a target value of $0.1\mathrm{m}{\mathrm{s}}^{-2}$, the ideal actuator would only need to supply a maximum force of around $600\mathrm{N}$ and a maximum rate of around $24\mathrm{kN}{\mathrm{s}}^{-1}$ under the chosen output constraints. Above a $1600\mathrm{N}$ maximum force limit and $82\mathrm{kN}{\mathrm{s}}^{-1}$, no further improvement in body acceleration RMS can be observed. The FFT plots of this scenario are shown in Appendix D, Figure A6 and Figure A7. No constraints were violated in this scenario.

4.2.3. Motorway at $180\mathrm{km}{\mathrm{h}}^{-1}$

The final scenario on the AB profile is driving on the motorway at $180\mathrm{km}{\mathrm{h}}^{-1}$. Again, qualitatively, the RMS surface plots in Figure 14 show the same behavior as in the $45\mathrm{km}{\mathrm{h}}^{-1}$ and $90\mathrm{km}{\mathrm{h}}^{-1}$ scenarios, but this time shifted to higher values. Furthermore, the RMS values only approach the values of the unrestricted actuator for the highest maximum rate limitations and maximum force limits above $1500\mathrm{N}$, and the body acceleration RMS value can only be reduced to $0.48\mathrm{m}{\mathrm{s}}^{-1}$.

In the contour plot in Figure 15, the gradient descent line is shifted upwards compared to the other two scenarios with the AB profile. However, as with the other scenarios, one limit cannot compensate for the other one. In this case, the target value of $0.1\mathrm{m}{\mathrm{s}}^{-1}$ cannot be reached. Instead, $0.6\mathrm{m}{\mathrm{s}}^{-1}$ is used as a target, both here and for the BC profile. To achieve this RMS body acceleration, the ideal actuator would need to have a maximum force of $1600\mathrm{N}$ and a maximum rate of $52.5\mathrm{kN}{\mathrm{s}}^{-1}$.

The FFT plots in Appendix D show that an actuator with a $1000\mathrm{N}$ maximum force and a maximum rate of $52.5\mathrm{kN}{\mathrm{s}}^{-1}$ can significantly lower the body acceleration and dynamic wheel-load variation up to $3\mathrm{Hz}$, while maintaining the wheel-load variation of the passive suspensions at $3\mathrm{Hz}$ and above. A further decrease in body acceleration, achieved with the unrestricted actuator, entails much higher dynamic wheel-load variation around the wheel eigenfrequency. In this scenario, the dynamic wheel-load constraint was the restrictive one. In Figure A16 in Appendix E, it can be seen that the most violations occurred in the simulations where the maximum force was higher than $1000\mathrm{N}$ and the rate limits were above $50\mathrm{kN}{\mathrm{s}}^{-1}$. However, their maximum and mean value is only $5\mathrm{N}$, which is negligible. The only relevant violations with mean values between $200\mathrm{N}$ and $300\mathrm{N}$ occurred in the simulation with the lowest rate or lowest force constraint. These only occurred in two to three time steps, which is admissible.

4.3. Road-Profile BC

Two scenarios were carried out on road-profile BC, which were both driven at $90\mathrm{km}{\mathrm{h}}^{-1}$. This corresponds to travelling on a substandard country road in Germany. In one scenario, the objective was again to minimize body acceleration, while in the other, it was to minimize the dynamic wheel-load variation.

4.3.1. Bad Country Road at $90\mathrm{km}{\mathrm{h}}^{-1}$

The RMS values in this scenario are comparable to the RMS values of the scenario on the AB profile at $180\mathrm{km}{\mathrm{h}}^{-1}$, as can be seen in in Figure 16. The only difference is that the body acceleration cannot be minimized so well for low maximum force and rate limitations, and the lowest possible body acceleration RMS value is slightly higher.

In the contour plot of RMS body acceleration in Figure 17, it can be seen that the gradient for the maximum force is steepest up to $1000\mathrm{N}$. Here too, the rate limitation has the biggest influence on the reduction in body acceleration. No relevant further improvement can be seen for maximum force values above $2600\mathrm{N}$. To achieve the target value of $0.6\mathrm{m}{\mathrm{s}}^{-2}$, the ideal actuator would need to posses a maximum force of $2000\mathrm{N}$ and a maximum rate of $67\mathrm{kN}{\mathrm{s}}^{-1}$ under the chosen output constraints. >From the FFT plots of body acceleration (Figure A10) and dynamic wheel-load (Figure A11) in Appendix D, it ca be seen that an actuator with a maximum force of $500\mathrm{N}$ and a maximum rate of $10\mathrm{kN}{\mathrm{s}}^{-1}$ is hardly able to reduce the body eigenfrequency peak compared to the passive vehicle, and it is this that constitutes the main difference to the $180\mathrm{km}{\mathrm{h}}^{-1}$ AB profile scenario. It is thus the main reason why actuators with low maximum force and rate limits have higher body acceleration RMS values on the BC profile. In terms of dynamic wheel-load constraint violations, this scenario is also similar to the $180\mathrm{km}{\mathrm{h}}^{-1}$ AB profile scenario, which means there are some minor violations for high force and rate limits and a few larger violations for low maximum force and rate limits, as can be seen in Figure A17. But unlike the other scenario, some violations of the suspension travel limits occurred here too (see Figure A18). This violation was the largest in the simulation with a maximum force of $300\mathrm{N}$ and a maximum rate of $5\mathrm{kN}{\mathrm{s}}^{-1}$, with a maximum value of $1.2\mathrm{cm}$ and a mean value of $0.5\mathrm{cm}$ over $0.19\mathrm{s}$, which is still acceptable. In general, the most suspension-deflection violations arose in the simulations with this maximum force or rate limit.

4.3.2. Minimization of Dynamic Wheel-Load on a Bad Country Road

This scenario investigates the requirements placed on the actuator with the objective of minimizing dynamic wheel-load variation. For this investigation, it was necessary to change the weighting factors in the cost function (see Table 8 in Section 3.2). The RMS plots in Figure 18 show that reduction in dynamic wheel-load goes hand in hand with an increase in body acceleration, which is the exact opposite compared to minimizing the body acceleration and clearly indicates that the conflict of objectives, as described in Section 2.3, also exists for active systems. Compared to minimizing body acceleration, the suspension-deflection is also reduced in this case, influenced by the weighting factor for suspension-deflection, which was increased for stability reasons. The actuator force RMS value is lower compared to the other scenario, which is influenced by the addition of the actuator force weighting.

Figure 19 shows the RMS contour plot of the dynamic wheel-load variation. This result indicates that reduction in RMS dynamic wheel-load variation to $1520\mathrm{N}$, which already represents $82\%$ of the possible reduction with an unrestricted actuator, is possible with a maximum force of about $750\mathrm{N}$, provided the maximum rate limit is higher than $42\mathrm{kN}{\mathrm{s}}^{-1}$. The ideal actuator for achieving an RMS value of $1420\mathrm{N}$ under the given constraints would need to have a maximum force of $2250\mathrm{N}$ and a maximum rate of $72\mathrm{kN}{\mathrm{s}}^{-1}$. No noteworthy improvement in dynamic wheel-load variation is seen for maximum force limits higher than $3000\mathrm{N}$ and rate limits higher than $100\mathrm{kN}{\mathrm{s}}^{-1}$. Any limits higher than $2250\mathrm{N}$ and $72\mathrm{kN}{\mathrm{s}}^{-1}$ can only reduce the dynamic wheel-load RMS value by a further $8\mathrm{N}$ at best.

The FFT plots of body acceleration (Figure A12) and dynamic wheel-load (Figure A13) in Appendix D show how the reduction in dynamic wheel-load around the wheel eigenfrequency leads to an increase in body acceleration between $4\mathrm{Hz}$ and $11\mathrm{Hz}$. It can also be observed, that the body eigenfrequency peak is shifted towards $2.5\mathrm{Hz}$ and that the dynamic wheel-load variation actually worsens between $0\mathrm{Hz}$ and $4\mathrm{Hz}$. The dynamic wheel-load (Figure A19) and suspension-deflection (Figure A20) constraint violations are comparable to those when optimizing for body acceleration. Dynamic wheel-load constraint violation was present in simulations with low maximum force or rate limitations. Suspension-deflection violation was worse again for simulations with a maximum force limit of $300\mathrm{N}$ or a maximum rate limit of $5\mathrm{kN}{\mathrm{s}}^{-1}$. Contrary to optimization for body acceleration, dynamic wheel-load constraint violations become smaller and disappear for high maximum force and rate limitations, due to the optimization target. Furthermore, the violations of the suspension-deflection constraint are slightly reduced for higher limits, due to the change in its weighting factor.

5. Discussion

The results of this research are discussed in Section 5.1. Points of criticism and potential for future work are outlined in Section 5.2.

5.1. Interpretation of the Results

In general, the results obtained in this research are plausible and were also found to be robust in repeated simulations. If a scenario was repeated multiple times with the same weighting factors and constraints, the results of the individual simulations did not change, which indicates that the optimal trajectory of control inputs for each simulation was unique. Of course the results depend on the chosen weighting factors and constraints. Here, the weighting factors were chose for either maximum comfort or maximum road-holding, so the results would of course change with different objectives. The other influential factors are the output constraints. They were set to maintain a positive wheel-load and to keep the suspension travel within $10\mathrm{c}\mathrm{m}$. These could also be changed to allow more suspension travel, where possible on the simulated vehicle, or to keep the minimum wheel-load above $0\mathrm{N}$, which again would change the obtained results. The purpose of this research article is to present a tool set and method for investigating these dependencies, and as with any simulation tool, it is the responsibility of the user to define reasonable boundary conditions and understand the validity and limitations of the results.

Another point worthy of discussion are the results obtained for low maximum force and maximum rate limits in the bump scenario, the $180\mathrm{km}{\mathrm{h}}^{-1}$ AB profile scenario, and the two scenarios on the BC profile. In these simulations, the controller had problems staying within the output constraints on dynamic wheel-load and/or suspension-deflection. This is not a problem of the controller itself, because when the higher limits on the actuator where higher, the controller remained within these constraints. The problem is that not even the passive vehicle stays within the chosen constraints, as can be seen in in Appendix F

Table A1. Constraint violations in the passive simulations.

Violation	Bump 36 km h${}^{-1}$	AB 90 km h${}^{-1}$	AB 45 km h${}^{-1}$	AB 180 km h${}^{-1}$	BC 90 km h${}^{-1}$
Max. $F_{\mathrm{dyn}}$	$0\mathrm{N}$	$0\mathrm{N}$	$0\mathrm{N}$	$906\mathrm{N}$	$2293\mathrm{N}$
Mean $F_{\mathrm{dyn}}$	$0\mathrm{N}$	$0\mathrm{N}$	$0\mathrm{N}$	$537\mathrm{N}$	$610\mathrm{N}$
Samples $F_{\mathrm{dyn}}$	$0$	$0$	$0$	$2$	$8$
Max $z_{\mathrm{sus}}$	$4.1\mathrm{c}\mathrm{m}$	$0\mathrm{c}\mathrm{m}$	$0\mathrm{c}\mathrm{m}$	$0\mathrm{c}\mathrm{m}$	$0\mathrm{c}\mathrm{m}$
Mean $z_{\mathrm{sus}}$	$2.1\mathrm{c}\mathrm{m}$	$0\mathrm{c}\mathrm{m}$	$0\mathrm{c}\mathrm{m}$	$0\mathrm{c}\mathrm{m}$	$0\mathrm{c}\mathrm{m}$
Samples $z_{\mathrm{sus}}$	$28$	$0$	$0$	$0$	$0$

. The controller has to balance the objective conflict of trying to stay within the constraints while minimizing body acceleration or dynamic wheel-load. The absolute values in these simulations are very dependent on the output constraint settings. This only occurred for the lowest rate and force limitations on the actuator and does not affect the general results in these scenarios.

Although the focus of this research is on the tool set and the method for investigating the influence of actuator limitations, the results can also be interpreted for the chosen output constraints and weighting factors. In general it could be seen that a relatively low maximum actuator force of around 1000 to $2000\mathrm{N}$ is sufficient, provided the actuator can supply a high enough rate of change. In most scenarios, a threshold for the maximum force was apparent, after which it was only possible to improve the result with higher rate limits.

Table A2. Optimal actuator limits for the chosen target values.

	Bump 36 km h${}^{-1}$	AB 90 km h${}^{-1}$	AB 45 km h${}^{-1}$	AB 180 km h${}^{-1}$	BC 90 km h${}^{-1}$	BC 90 km h${}^{-1}$ $\left({\mathit{F}}_{\mathit{dyn}}\right)$
Target	$0.1\mathrm{m}\mathrm{s}$−2	$0.1\mathrm{m}\mathrm{s}$−2	$0.1\mathrm{m}\mathrm{s}$−2	$0.6\mathrm{m}\mathrm{s}$−2	$0.6\mathrm{m}\mathrm{s}$−2	$1520\mathrm{N}$
$F_{\mathrm{Ac}},\max$	$2000\mathrm{N}$	$1400\mathrm{N}$	$600\mathrm{N}$	$1600\mathrm{N}$	$2000\mathrm{N}$	$2250\mathrm{N}$
${\dot{F}}_{\mathrm{Ac}},\max$	$21.5\mathrm{k}\mathrm{N}\mathrm{s}$−1	$62.0\mathrm{k}\mathrm{N}\mathrm{s}$−1	$24.0\mathrm{k}\mathrm{N}\mathrm{s}$−1	$52.5\mathrm{k}\mathrm{N}\mathrm{s}$−1	$67.0\mathrm{k}\mathrm{N}\mathrm{s}$−1	$72.0\mathrm{k}\mathrm{N}\mathrm{s}$−1

in Appendix F gives an overview of the maximum actuator force and rate limits obtained by choosing arbitrary target values for the body acceleration and dynamic wheel-load RMS. Again, the target value needs to be determined by the user of this method. For the chosen target values of $0.1\mathrm{m}{\mathrm{s}}^{-2}$, $0.6\mathrm{m}{\mathrm{s}}^{-2}$ and $1520\mathrm{N}$, the actuator needs to supply a maximum force within $600\mathrm{N}$ to $2250\mathrm{N}$ and a maximum rate within $21.5\mathrm{kN}{\mathrm{s}}^{-1}$ to $72.0\mathrm{kN}{\mathrm{s}}^{-1}$, where the higher force and rate requirements are due to optimization of dynamic wheel-load variation on the BC profile, and the lowest stem from the inner city and single-obstacle scenarios. An actuator with a maximum force of $2000\mathrm{N}$ and a maximum rate of $60\mathrm{kN}{\mathrm{s}}^{-1}$ would therefore be a suitable compromise for these scenarios. Comparing this to the eABC system recently introduced by Daimler [19], which has maximum force limits of $6000\mathrm{N}$ and $7000\mathrm{N}$, and maximum rate limits of $24\mathrm{kN}{\mathrm{s}}^{-1}$ and $20\mathrm{kN}{\mathrm{s}}^{-1}$ for the front and rear axles of a Mercedes GLE, it can be seen that in terms of ride comfort due to road irregularities, it might have been beneficial to design a system with a lower maximum force and higher rate of change capabilities. However, for requirements on roll and pitch compensation under lateral and longitudinal acceleration, high static forces need to be supplied and the actuator dynamics are of less importance, but with regard to autonomous driving, where vehicles are expected to travel with lower longitudinal and lateral accelerations, the focus might need to shift towards systems with higher rate capabilities and lower maximum forces.

5.2. Critisism and Future Work

This research is the first to present a method for the investigation of actuator limitation influences on ride comfort. There are some things to keep in mind when applying this method. First, it is a simplification to assume that an actuator is perfectly characterized by its maximum force and a constant maximum rate limit. Usually, actuators in suspension control systems show PT-1 or PT-2 behavior as a step response, because the actuator itself has a controller that regulates the output force and is tuned to display stable behavior. This means that the actual rate limitation is a function of the desired step height and the actuator force output is also subject to inaccuracies. Neither are currently considered, but the method can be used to analyze the general sensitivity to these constraints. Another thing to keep in mind is that if the vertical vehicle dynamics controller is of the standard feedback type, such as a sky-hook controller, the actuator rate also influences its stability. The rate of change is even more important, because it is the main reason for the phase delay from the desired to the actual force, and the greater the phase delay, the lower the possible gain of the controller. The aim is to extend this method in the future, so as to encompass actuator accuracy, power limitation, and dynamic actuator behavior using a full vehicle model; this also including roll and pitch dynamics, and the possibility of an active anti-roll bar control system. It is also planned to extend the method to analyze semi-active suspension systems and to compare these results with those from active suspension systems.