Hybrid Model Predictive Control of Semiactive Suspension in Electric Vehicle with Hub-Motor

: In hub-motor electric vehicles (HM-EVs), the unbalanced electromagnetic force generated by the HM will further deteriorate the dynamic performance of the electric vehicle. In this paper, a semiactive suspension control method is proposed for HM-EVs. A quarter HM-EV model with an electromechanical coupling effect is established . The model consists of three parts: a motor model, road excitation model and vehicle model. A hybrid model predictive controller (HMPC) is designed based on the developed model, taking into account the nonlinear constraints of damping force. The focus is on improving the vertical performance of the HM-EV. Then, a Kalman filter is designed to provide the required state variables for the controller. The proposed control algorithm and con-strained optimal control (COC) algorithm are simulation compared under random road excitation and bump road excitation, and the results show that the proposed control algorithm can improve ride comfort, reduce motor vibration, and improve handling stability more substantially.


Introduction
Due to increased pressure for environmental protection and energy shortages, electric vehicles' development is being highly valued by many research institutions [1]. Electric vehicles are environmentally friendly and excel in terms of sporting performance due to the fast and precise torque of the electromotor [2,3]. Hub-motor electric vehicles (HM-EVs) have the advantages mentioned above and eliminate the need for a powertrain, improving transmission efficiency and space utilization, making power control easier to achieve [4,5]. Thanks to its unique technical advantages, HM-EV has received unprecedented attention and development. It has become one of the focuses and hotspots of electric vehicle technology research and has become a unique direction for electric vehicles' development. The industry also refers to the HMs as the final drive form of electric vehicles [6].
However, because the HM-EV highly integrates the motor, reduction mechanism, and brake in the wheel, it increases the unsprung mass of the vehicle and deteriorates the vertical performance of the vehicle [7]. Additionally, tire bounce, support shaft bending, bearing wear caused by the excitation of the road will cause the stator and rotor to be eccentric and cause uneven distribution of the motor air gap along the circumference [8,9]. The unbalanced electromagnetic force generated by the uneven air gap will be transmitted to the wheels through the reduction mechanism or directly, which is equivalent to adding a certain external load to the wheels, thus affecting the dynamic characteristics of the suspension [10]. To solve these problems, many scholars have suppressed various adverse effects of HMs on vehicles from both structural and control perspectives. Nagaya et al. [7] used special springs and dampers to connect the motor to the unsprung mass, to turn the motor into a dynamic vibration absorber (DVA) for the unsprung mass, reducing the vibration of the sprung mass and improving the smoothness of the vehicle ride, as well as reducing the input to the motor vibration. Luo et al. [11] proposed a new type of HM structure with rubber bushing between the stator and rotor. Through the comparison and simulation of two HM schemes with bushing and without bushing, it was concluded that the vertical performance of the motor with bushing was significantly improved. Qin et al. [12] proposed a new method of vibration damping based on the DVA structure, and the simulation results show that the ride comfort and handling stability were effectively improved after adopting the new DVA structure. Finally, a multibody simulation was carried out using muti-body dynamic software (LMS motion) to verify the feasibility of the proposed DVA structure. There are many scholars already working on suppressing the adverse effects caused by HM from an active suspension perspective. Liu et al. [13] proposed an optimal control strategy for electric vehicle wheels, which optimizes the parameters of the DVA structure according to the typical working conditions, and then adopts the fuzzy control method and linear quadratic regulator (LQR) method to control the actuator force of the DVA structure and suspension force, respectively, to reduce the impact on the HM and improve the ride comfort. Shao et al. [14] present an HM-EV active suspension fault-tolerant fuzzy H∞ controller that takes into account sprung mass variation, actuator faults, and control input constraints, and verifies the effectiveness of the controller by simulation. He et al. [15] proposed a fuzzy optimal sliding mode control method using normalization and hierarchical analysis to select weighting coefficients, and verified the superiority of using this control method for electric vehicles by simulation. Li et al. [16] proposed a multiobjective optimization control method for active suspension, using a particle swarm algorithm to generate the optimal parameters of the active suspension and the simulation results show that the optimized active suspension can effectively reduce the unbalanced electromagnetic force while simultaneously turning the motor eccentricity so it is kept within a reasonable range. Wu et al. [17] proposed an H2 active suspension control scheme with good robustness, and the simulation results show that the suspension system has a satisfactory control effect under complex electromagnetic force excitation and dynamic perturbation. Shao et al. [18] proposed a hybrid controller consisting of an H∞ suspension controller and a switched reluctance motor controller, which consists of current chopping control and pulse width modulation control, and simulation results show that the proposed control method can effectively reduce the unbalanced electromagnetic force and air gap eccentricity. Liu et al. [19] proposed a two-stage optimization method based on the mathematical model of a quarter vehicle with active suspension equipped with a dynamic vibration absorber. Firstly, an LQR controller was designed to suppress the vibration of the wheel motor. Then, a controller with finite frequency H was designed to improve the ride comfort of the vehicle. Although the active suspension exhibits an excellent performance, its structure is complex and expensive, limiting the promotion of its application. The emergence and mature development of controllable dampers has led to a greater focus on semiactive suspension. Xu et al. [20] proposed a hybrid controller based on the hybrid acceleration drive damping algorithm and used the multiobjective optimization method to determine the parameters of the controller, and finally the simulation was carried out under the random road excitation and bump excitation simulation; the results show that the proposed control method can effectively reduce the vibration of the sprung mass and motor. Mauricio et al. [21] used four different semiactive suspension controllers in an HM-EV and evaluated and compared the whole vehicle to show that it is feasible to use four different semiactive suspension controllers for the vehicles. From the available research, it can be seen that there are few studies on semiactive suspension in HM-EVs.
This research proposes a damping control method for the HM-EV semiactive suspension to improve the ride comfort, handling stability and reduce the eccentricity between the stator and rotor. The primary contributions of this study are as follows: (1) developed a hub motor direct drive-air suspension cooperative system; (2) designed a hybrid model predictive controller to suppress the vibration of the entire system; (3) designed a Kalman filter to provide state variables that cannot be measured for the hybrid model predictive controller.
The rest of this article is organized as follows: In the second section, the establishment of the 1/4 cooperative system model is introduced. Then, the control problem is proposed in the third section. To address this problem, a hybrid model predictive controller and a Kalman filter are designed. In the fifth section, a numerical simulation is performed to verify the effectiveness of the proposed method and summarize the conclusions.

Electromechanical Coupling Model
HM-EVs are complex electromechanical coupling systems, which includes HM, suspension, and wheels. In this research, an integrated model of a quarter HM-EV is established, which contains three parts: the permanent magnet brushless direct current (PMBLDC) motor model, which provides an unbalanced electromagnetic force to the wheels and suspension; vehicle model, which emphasizes the vehicle dynamic response; road excitation model, which provides external excitation to the system.

BLDC Motor Model
In recent years, the PMBLDC motor has received more interest in automotive applications due to the higher reliability, lower maintenance and quieter operation the BLDC motor has compared to its brushed DC counterpart [22]. This paper uses an outer rotor PMBLDC motor as the hub drive motor.
The current equation of the outer rotor PMBLDC is as follows: where ua, ub, and uc represent the three-phase winding phase voltage, R is the stator phase resistance, L is the winding inductance, M is the mutual inductance, pop is the differential operator, ia, ib, and ic represent the three-phase winding phase current, and ea, eb, and ec represent the three-phase winding back electromotive force. According to Equation (1), the electromagnetic torque equation and mechanical motion equation of motor are as follows: J n dω n dt where Te is the electromagnetic torque of the motor, TL is the load torque of the motor, Jn is the moment of inertia of the motor, B is the friction coefficient, and ωn is the angular velocity of the rotor. Due to the complex installation conditions and mechanical environment of the hub motor on an electric vehicle, the motor is under vertical load during operation, which causes the motor stator-rotor to be eccentric and thus generates unbalanced electromagnetic force. The motor noneccentricity and eccentricity are shown in Figure 1.
The air gap magnetic field of a PMBLDC motor includes the permanent magnet magnetic field and the armature reaction magnetic field, which can be added together by the linear superposition method to obtain the radial air gap magnetic field and the tangential air gap magnetic field of a PMBLDC motor-the equations are as follows [23]: B er (r,α,t)=[B mr (r,α,t)+B ar (r,α,t)]ε δ (4) where r is the radius at the intended position, α is the stator angle, t is time, Ber is the radial air gap magnetic field in the eccentric state, Bmr is the radial magnetic field of the permanent magnet, Bar is the radial magnetic field of armature reaction, Bet is the tangential air gap magnetic field in the eccentric state, Bmt is the tangential magnet field of the permanent magnet, Bat is the tangential magnetic field of armature reaction, and εδ is the correction coefficient of magnetic conductivity in the eccentric state. Finally, the vertical and longitudinal unbalanced electromagnetic forces caused by the eccentricity between the stator and rotor of the motor are calculated as follows: where L is the axial length of the motor and µ0 is vacuum permeability.

Quarter Vehicle Model
Despite its simplicity, the quarter vehicle model is widely used in vertical dynamics analysis and controller design [24]. Figure 2 depicts the physical structure and equivalent model of the quarter vehicle model for HM-EV. The dynamics equations can be described according to Newton's law of motion as follows: where ms is sprung mass, z s is vertical acceleration of sprung mass, z s is vertical velocity of sprung mass, zs is vertical displacement of sprung mass, mus is stator and residual unsprung mass, z us is vertical acceleration of stator and residual unsprung mass, z us is vertical velocity of stator and residual unsprung mass, z us is vertical displacement of stator and residual unsprung mass, mur is rotor, rim and tire mass, z ur is vertical acceleration of rotor, rim and tire mass, z ur is vertical velocity of rotor, rim and tire mass, zur is vertical displacement rotor, rim and tire mass, cs is damping coefficient, ks is spring stiffness, kb is spring stiffness, kt is tire stiffness, q is road excitation and udam is damping force. The air suspension has a low frequency of vibration, can change the height of the vehicle and can improve control stability and ride comfort [25]. Therefore, in this paper, air springs have been chosen as the elastic element of the suspension system. Hengjia states that the differential equation for the change in air pressure in an air spring with time can be expressed as [26]: Assuming that the air pressure, volume and air mass changes inside the air spring are much less than its initial state, the above equation can be linearized as: Xinbo indicates that the air spring model can be ultimately represented as [27]: p as A as (z as0 +z s -z us )=κRT m -κp as A as (z s -z us ) Considering that the rate of change in air pressure and suspension dynamic deflection has a large influence on the dynamic characteristics of the system and therefore is retained, and the current air pressure and current suspension dynamic deflection are replaced by the initial air pressure and initial suspension dynamic deflection, the linearized air spring model can be represented as: p as A as z as0 =κRT m -κp as0 A as (z s -z us ) (11) where pas is the air pressure of the air spring, Vas is the volume of the air spring, Vas0 is the initial volume of the air spring, mas is the air mass in the air spring, mas0 is the initial air mass of the air spring, κ is the air polytropic exponent, R is the perfect gas constant, T is the air temperature, zas0 is the initial air spring displacement, Aas is the air spring crosssection area (assuming it remains constant) and qm is mass flow of air flowing into the air spring.
The air spring is not inflated or deflated in this paper, so the air spring model can be represented as: p as A as z as0 =-κp as0 A as (z s -z us ) (12) From the above equation, it follows that p as A as =-κp as0 A as 2 (z s -z us )/V as0 (13) The gas spring stiffness can therefore be expressed as: As the main components of the suspension system, dampers have a significant effect on the smoothness and handling stability of the vehicle. The damping coefficient of traditional dampers is generally fixed, making it difficult to meet the performance requirements of vehicles under different working conditions. In this paper, a controllable dampener was used, which can achieve a constantly adjustable damping coefficient by adjusting the opening of the damping valve. The controllable damper force Fcs can be expressed as a function of the control current and the relative velocity of the suspension. The formulas are described as follows [28]: Many methods can be used to describe such a function and this paper adopts a nonparametric model [29], in which the controllable damper force is expressed as: where A(i) is the maximum damping force at a given control current, an is a polynomial parameter, k is the largest order of the polynomial, Sb is a shape function, b0 and V0 are the shape function parameter. The parameters are shown in Table 1 [28] with the force-velocity characteristics depicted in Figure 3.  Think of the damping coefficient of a system as a combination of uncontrollable cmin and controllable parts [0, cmax-cmin]. According to JASOC602:2001 [30], the damping force when the piston speed is 0.31 m/s represents the equivalent damping force, and the corresponding equivalent damping coefficient can be obtained accordingly. For a controllable damper, the nonadjustable part of the damper is always at work, so it is assumed that cs=cmin. Parameters of the quarter HM-EV model are shown in Table 2

Road Excitation
The road excitation is the main excitation of the suspension system, so a suitable road excitation model is needed for control effect verification. The change in the height of the pavement relative to the reference surface along the length of the roadway is commonly referred to as the pavement unevenness function. In this paper, two types of road excitation were used: random road excitation and bump road excitation.
Random road surface unevenness can be represented by the power spectral density (PSD) function [31]. The PSD of velocity in the spatial frequency domain of the road surface is shown in the following equation where n is the spatial frequency, n0 is the reference spatial frequency, Gq(n0) is the power spectral density at the reference spatial frequency and v is the vehicle velocity. The road speed signal is white noise with a finite bandwidth, so the white noise can be passed through a filter to generate a time-domain signal of the road elevation. The filter is shown in the following equation where is the cutting frequency and = 0.01 HZ. Road excitation at 20 m/s on a Class C road and on a Class B road is shown in Figure 4. Then, the road time-domain signal can be expressed as follows: In order to study the dynamic response characteristics of the HM-EV in this paper introduces bump road excitation as shown in the following equation.
where a represents the height of the bump, l represents the length of the bump and v represents the vehicle velocity. In this paper, two different heights of bump road were selected as excitations: a = 0.1 m and a = 0.02 m. For other parameters, v = 1.5 m/s, l = 0.85 m was chosen. Bump road excitation is shown in Figure 5.

Hybrid Model Predictive Controller Design
Suspension performance directly affects passenger comfort and vehicle driving safety. Conventional suspensions have difficulty adapting to changes in the road surface, speed and load due to their stiffness and damping, making it difficult to meet the high demands of vehicle performance in changing conditions. This paper proposes a hybrid model predictive controller for the nonlinear constraint of damping forces to improve the overall system vertical performance deterioration due to the HM.
For traditional control methods, the optimization problem with state constraints has been a research challenge [32]. Model predictive control (MPC) can handle optimization problems with constraints in a finite time domain. At each sampling time interval, MPC transforms the system control problem into an open-loop optimal control problem solving a finite time domain with a step length of N according to the amount of state of the system x(k) at the current moment, obtains the optimal control variables {u(k), u(k + 1), ..., u(k + N -1)} within N steps by an online calculation method, applies the first step control amount u(k) to the system and corrects the future output at the next moment according to the error between the system output and the predicted output to realize the closed-loop of the overall rolling optimization. The control procedure is shown in the Figure 6. Subject to where Q and R are the corresponding weighting matrixes. The optimization problem described above can be translated into a standard form of linear quadratic programming.

Subject to
A e x=b e A j x≤b j Equation (7) can be expressed in the form of a state-space: Since MPC solvs the problem in the discrete time domain, it is necessary to discretize the existing continuous time state equations. The equation of state after discretization can be written as follows: where For the quarter vehicle model for HM-EVs, the evaluation indicators are usually the acceleration of vertical vibration of the sprung mass z s , the eccentricity of the motor z usz ur and tire dynamic load k t (z us -q). The acceleration of the vertical vibration of the sprung mass represents the comfort of the vehicle; the eccentricity of the motor represents the safety of the motor-too large an eccentricity will damage the motor, but also increase the electromagnetic torque fluctuation and deteriorate the comfort and handling stability of the vehicle; the dynamic load of the tire represents the grounding of the tire, and its value is smaller, indicating that the vehicle handling stability is better.
Excessive suspension dynamic deflection can lead to an increase in the probability of hitting a limit during its movement, so it is necessary to constrain the suspension dynamic deflection. Suspension limit deflection is typically ±60 mm, and the probability of the suspension hitting the limit is 0.3% when the suspension limit deflection is three times larger than the suspension dynamic travel. Therefore, it is necessary to limit the suspension dynamic deflection to ±20 mm. For the damping force, the direction of the damping force must be the same as the direction of the relative speed of the suspension and must also meet the maximum damping force constraint. The constraints are shown in the following equation: -20 mm≤x 4 (k)≤20 mm 0≤u(k)(x 1 (k)-x 2 (k))≤(c max -c min )(x 1 (k)-x 2 (k)) 2 (25) Since both pavement excitation and unbalanced electromagnetic forces are unmeasurable, the chosen predictive equations are as follows: The optimization problem can be expressed as the following equation: In the above constraint, the computational equation used for the constraint is related to the logical relationship between state variables, so not only the magnitude of the constraint has to be satisfied, but the correct constraint also has to be switched according to the logical relationship between the state variables. It is clear that traditional linear MPC does not handle these types of constraints. In order to better deal with the above constraints, a hybrid model in the damping adjustment process needs to be established.
Due to the complex engineering background, there is no unified way to construct hybrid system models for different research fields. The main models currently used to describe hybrid systems are: stepwise structural hybrid system models, hybrid automata models, piecewise affine models, mixed logical dynamical (MLD) models and linear complementary models [33].
The MLD model is a modeling framework for hybrid systems proposed by Bemporad and Morari in 1999, which can effectively solve the problem of modeling continuous dynamic processes coupled with discrete events in hybrid systems, and at the same time, with certain theoretical derivations, the MLD model can also be transformed into other model structures of hybrid systems, thus it has a strong generality. The MLD-based model can effectively handle the integrated problem of optimizing control laws for hybrid systems. Therefore, in this paper, the MLD modeling approach was used to model the mixing dynamics in the damping regulation process.
The standard form of the MLD model is as follows [34]: where x is the state variable, y is the output variable, u is the system input and δ is the auxiliary variable. x,y,u and δ can be continuous and logical variables. Two logical and continuous auxiliary variables need to be introduced in order to accurately describe the nonlinear constraints.
where δ v , δ u are auxiliary continuous variables.
The hybrid system description language (HYSDEL) can be used to translate the above problems into a standard MLD model.
The optimization problem can then be transformed into the following form Subject to

x(k+1)=Ax(k)+B 1 u(k)+B 2 δ(k)+B 3 z(k) y(k)=Cx(k)+D 1 u(k)+D 2 δ(k)+D 3 z(k) E 2 δ(k)+E 3 z(k)≤E 1 u(k)+E 3 x(k)+E5
where Q = and N is the prediction step. The following vectors are defined: Then, the prediction equation of the system can be written as: The performance metrics of the system can be further expressed as: where The constraints can be further translated as follows: By removing the constant term, the above problem can be transformed into a mixed integer quadratic programming problem, which is mathematically described as follows: For the mixed-integer quadratic planning problem, the main methods are: the branch delimitation method, genetic algorithm, etc. In this paper, we use the branch delimitation method and use the Gurobi solver to improve the solution speed.

HM-EV Semiactive Suspension Kalman Filter Design
In the hybrid model predictive controller design, we assumed that all state variables were measurable, but, in reality, not all state variables are measurable. Some noise and interference is inevitable in the input and output measurement data of the system. It is then necessary to reconstruct the system using a state observer and correct the state and output values of the reconstructed system by the error between the measured output and the observer output, so that the state and output values of the observer approximate the real state and output of the system.
The equation of state for a multivariate linear system in a stochastic environment, as described by the Kalman filter, is as follows: ω,v are a white process noise and measurement noise, which are subjected to a multivariate Gaussian distribution and are independent of each other, as follows: where Q is the process noise covariance matrix and R is the measurement noise covariance matrix.
The principle of Kalman filtering is to use the Kalman gain to correct the state prediction so that it approximates the actual value.
The state variable can be estimated by the following equation: where the optimal observer gain L can be calculated by the discrete Riccati equation.

M(k)=P(k)C T (CP(k)C T +R) -1 Z(k)=(I-M(k)C)P(k)(I-M(k)C) T +M(k)RM T (k)
The Kalman filter minimizes the state estimation error covariance by solving the discrete algebraic Riccati equation to obtain the optimal gain matrix L(k).
For an HM-EV semiactive suspension system, the measurable outputs are the acceleration of vertical vibration of the sprung mass and the suspension dynamic deflection. Thus, using these two measurable via the Kalman filter to estimate the state of the entire system. Under the disturbance of white noise, the observer results under C-class road excitation are shown in the figure below.
As shown in Figures 7-16, for each state variable, the difference between the estimated and actual values of the Kalman filter is very small under different noise disturbances and can meet the control requirements.

Simulation Results
In order to verify the effectiveness of the proposed control strategy, a simulation analysis of the control strategy is required. The overall control framework is shown in Figure 17 as follows: First, the HM-EV model inputs its measurable sprung mass acceleration and suspension dynamic deflection into the Kalman filter, which estimates the system's state variables in real-time. Then, the MLD model generates the objective function from the values of state variables obtained from the Kalman filter in a predictive iterative process combined with a reference trajectory, which is then transformed into a standard form of a mixed integer quadratic program (MIQP) problem to be solved. Then, the damping reverse model generates a control current to the dampers based on the optimal damping force output from the controller and the relative velocity of the suspension. Finally, dampers generate the damping force provided to the model, thus closing the loop on the entire control strategy. For comparison purposes, the clipped optimal control (COC) was selected as the comparator. In active suspension systems, optimal control can be achieved by offline calculations or by updating the optimal feedback matrix online according to driving conditions. For the semiactive suspension system, the equivalent mathematics is described as a complex nonlinear problem due to its controllable force energy dissipation constraint and state-dependent control force time-varying constraint [35]. Typical COC can be obtained as follows: where Q is the semipositive matrix and R and N are the positive definite matrixes. The system dynamics model can be represented as u(t) can be expressed as linear state feedback.
where the state feedback gain matrix can be calculated by the following equation: Thus, the optimal control when unconstrained is as follows: There is a constraint on the damping force of the semiactive suspension, so the final output control force after passing the constraint is expressed as where denotes a damping force constraint. Assuming the vehicle is driving at 20 m/s on a C-class road surface and at 20 m/s on a B-class road surface, the damping coefficient of the passive suspension is 3000 N·s·m −1 . In order to make the results comparable, for the selection of the weight coefficients of the controller, the trial-and-error method was used to make an appropriate selection to ensure that the performance of the system could be optimized under the control of the controller. For the weighting coefficients of the hybrid model predictive controller (HMPC) controller, this paper chose Q = diag ([1, 10 × 10 7 , 10 × 10 4 ]) and R = 0.001. For the weighting coefficients of the COC controller, this paper chose Q = diag ([2.5 × 10 4 , 2 × 10 11 , 5 × 10 9 ]) and R = 0.005. Figures 18-20 shows vertical acceleration of sprung mass, motor eccentricity and tire dynamic load with a time under 20 m/s of C-class road excitation. Table 3 gives the root mean square (RMS) for the suspension performance index under B-class and Cclass road.    Table 3, both COC and HMPC effectively improve the system performance compared to passive suspension under, but the HMPC can improve the vertical performance of the vehicle to a greater extent in comparison. The improvement of the evaluation index by HMPC is approximately the same for both the Cclass road with a larger amplitude and the B-class road with a smaller amplitude, which proves the effectiveness of HMPC.   In order to better compare the effect of different control methods, peak-to-peak (PTP) values were defined which represent the maximum values minus the minimum values. Table 4 shows the PTP values of each evaluation index under two types of bump road excitation.  Table 4 that HMPC can reduce the amplitude of the evaluation index to a greater extent under bump excitation compared to COC. In addition, the improvement of PTP values under different amplitudes bumps excitation of HMPC so it is basically remains the same, which reflects that HMPC can perform well under different amplitude of road excitation, further verifying the effectiveness of HMPC.
Tire dynamic load (N)

Conclusions and Discussion
In this paper, we proposed a semiactive suspension control method that can consider the nonlinear constraint of damping force to suppress the vibration of the sprung mass and motor in order to address the deterioration of the vertical performance of HM-EV. A state observer was designed based on the Kalman filter algorithm for HM-EV to meet the control needs of the controller. The main conclusions are as follows: (1) The state observer designed based on the Kalman filter algorithm can effectively accomplish the observation task. Simulations were performed to verify the observation effect under white noise random roads and bump roads, and the results show that the errors between the observed and actual values are small under both road excitations.
(2) The proposed control algorithm can effectively improve the vertical performance of the vehicle. Simulation analysis was carried out and compared with the COC control algorithm under different magnitudes of random road excitation and bump road excitation. The results show that the proposed control algorithm can effectively reduce the sprung mass and motor vibration regardless of the excitation magnitude and outperforms the COC control algorithm. The effectiveness of the proposed control algorithm was verified.
The COC control method ignores the constraints in the optimization calculation and then feeds back to the system after limiting the calculated forces according to the current state, which destroys the optimality of the control. The HMPC control algorithm was able to consider the nonlinear constraints of the damping force directly in the controller, but this also greatly increased the computational effort of the controller and increased the computational time. The HMPC algorithm requires online optimization calculations which increases the actual workload and computation time of the control system and violates the real-time requirements of the control during the operation of the actual vehicle. In order to shorten the computation time, the online computation process of the control law can be transformed into a simple offline table look-up process using multiparameter planning techniques so that the proposed control algorithm can be applied to the real vehicle. This is one of the future research directions.