# Lateral Stability Control of Four-Wheel Independent Drive Electric Vehicles Based on Model Predictive Control

Authors: Bin Huang; Sen Wu; Song Huang; Xiang Fu
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

Four-wheel independent drive electric vehicle was used as the research object to discuss the lateral stability control algorithm, thus improving vehicle stability under limit conditions. After establishing hierarchical integrated control structure, we designed the yaw moment decision controller based on model predictive control (MPC) theory. Meanwhile, the wheel torque was assigned by minimizing the sum of consumption rates of adhesion coefficients of four tires according to the tire friction ellipse theory. The integrated simulation platform of Carsim and Simulink was established for simulation verification of yaw/rollover stability control algorithm. Then, we finished road experiment verification of real vehicle by integrated control algorithm. The result showed that this control method can achieve the expectation of effective vehicle tracking, significantly improving the lateral stability of vehicle.


---

## Body

## 1. Introduction

With independent controllable drive/brake moment, rapid moment response, and measurable torque and speed, four-wheel independent drive electric vehicle has advantages in improving vehicle stability. When the vehicle has the risk of instability or has lost stability, it is difficult for the ground to maintain vehicle tracking and attitude adjustment with sufficient lateral force. Meanwhile, the vehicle based on active steering has little effect on driving direction control [1, 2]. In differential braking/driving process, the longitudinal force of the tire is adjusted to control yaw stability. It still has significant control effect even when the vehicle yaw movement loses stability. Differential braking is integrated with differential driving to broaden the range of active yaw moment, thus reducing the decrease in vehicle speed and keeping the vehicle in a stable running state.Domestic and foreign scholars have made lots of researches on yaw stability control based on distributed structure electric vehicle. The developed control strategy can be divided into yaw moment decision and wheel torque distribution control layers. References [3–5] focus on achieving direct yaw control by distributing the rear wheel driving or braking force. However, a single wheel produces limited yaw response with limited adhesion during the process of steering instability. Effective use of the adhesion of multiple wheels can greatly enhance the stability. In [6], the concept of additional yaw torque control is proposed. Based on independent control of braking force, the in-wheel motor is integrated with the rear wheel drive force control. Relying on the yaw control of driving and braking forces, the yaw response decreases at least 10% of the delay, and the lateral acceleration rises by up to 40%. In [7–11], the yaw moment decision is completed by adaptive control, fuzzy logic, and so on. Braking and driving coordination control can be used to achieve faster yaw response. However, these researches only apply simple method for hierarchical control. Only controlling the current state makes the existence of improved space for system robustness and accuracy. In [12], the classical sliding mode control theory is used to calculate the yaw moment required for dynamic control of in-wheel motor driving vehicle. The sliding surface is formed by the yaw rate error and the error change rate. The torque outputs of the motors are coordinated to achieve the above yaw moment value. When the system parameters are uncertain, this control method has good robustness. However, the high-frequency switching of sliding mode controller introduces a certain degree of “chattering” phenomenon, which affects the control precision of vehicle dynamics.For the vehicles with high centroid, when they run on the high adhesion road with high speed, the emergency avoidance operation can easily cause the risk of rollover because of large lateral acceleration. The rollover stability control (RSC) is as important as the yaw stability control (YSC). Fujimto and Hori’s team of University of Tokyo have done a lot of related research on their developed small test vehicle. In [13, 14], RSC and YSC are designed by utilizing two-degree-of-freedom (2-DOF) control based on disturbance observer (DOB); the tracking capability and robustness for lateral acceleration disturbance against such as side blast are realized by the proposed DOB. And they further propose electronic stability program based on DOB, which achieves integrated three dimensional vehicle motion control [15]. The DOB is based on the model following control; it works very well if the velocity is constant and not so high. But, at high speed, the control effect is reduced. The purpose of this paper is to propose integrated yaw/rollover stability control; MPC-based four-wheel braking/driving torque coordinated control strategy was proposed to improve the driving stability of four-wheel independent electric drive off-road vehicle. The proposed MPC, which is independent of exact mathematical model, compared with the DOB can also overcome the influence of modeling errors and environmental uncertainties. MPC controls the current vehicle state while predicting the states in sampling periods. The current input is corrected to obtain more prominent control accuracy, stability and robustness. Therefore, it is suitable for solving nonlinear, multiple input, and output problems [16].In the work, MPC-based four-wheel braking/driving torque coordinated control strategy was proposed to improve the driving stability of four-wheel independent electric drive off-road vehicle. According to MPC system structure, we established predictive control-oriented vehicle dynamics model. After that, the work improved vehicle yaw and rollover stability under limit conditions by taking driver’s input as the reference and road adhesion as constraint. Vehicle braking/driving force was allocated by MPC to minimize rollover index LTR, the deviations of yaw rate and side slip angle. Then, we used Carsim and Simulink joint simulation platform including professional vehicle and road models for yaw stability control under low adhesion coefficient and double lane-change conditions as well as rollover stability control algorithm validation under high adhesion coefficient and hook conditions. At last, real vehicle test verification was performed under high speed double lane-change condition.


## 2. Lateral Stability Control Strategy

The control goal of lateral stability control strategy enables the vehicle to run in the control of the driver without the risk of rollover, achieving convenient operation and stable driving. The lateral stability of the vehicle includes both the yaw and the roll stabilities. Scholars take linear yaw characteristic deduced from two-degree-of-freedom vehicle model as the ideal characteristic of the vehicle. In the work, we also select it as the control target of ideal yaw stability of the vehicle. For rollover stability, it is impossible to obtain the ideal rollover characteristic by establishing a simplified vehicle model. The current rollover condition of the vehicle is characterized only by selecting the appropriate rollover index. Besides, a suitable predicted value is selected to measure the possibility of rollover. In the rollover stability control, the expected rollover possibility should be as small as possible.


### 2.1. Yaw Stability Control Analysis

The yaw stability is characterized by state variables including side slip angle and yaw rate. Therefore, we determine the yaw stability of vehicle, and whether to carry out yaw stability control according to the two state variables.When the side slip angle of vehicle satisfies the following equation, the vehicle is in a stable zone; otherwise, it is in an unstable zone, and the stability control must be added.(1)β+B1β˙≤B2,where B1 and B2 are the parameters relevant to vehicle characteristics and driving state. In the work, we refer to the method of determining the boundary values of side slip angle by vehicle speed and pavement friction coefficient in [17]. Then, the relation between B1 and B2 is obtained by fitting.(2)β˙+B1β≤B2,B1μ=15.62μ2-34.37μ-6.719,B2V,μ=0.0002343μ2V2-0.000516μV2-0.7498μ2+1.650μ,where μ is the tire-road friction coefficient and V is the travel speed. Equation (2) is used to determine the stability boundary of side slip angle. After that, we obtain the estimated values of side slip angle and its derivative through the state estimation procedure which has been developed and validated by the authors in [18], thus deriving actual state point of the vehicle. The driving state of vehicle can be determined and expressed as (3) according to the distance between actual state point and stable boundary.(3)ΔD=β+B1β˙-B2.According to (1)–(3), we obtain the vehicle stability judgment phase diagram based on the side slip angle (see Figure 1).Figure 1
Vehicle stability judgment phase diagram based on the side slip angle.If the side slip angle and yaw rate of the vehicle satisfy the following equation, then the vehicle is in a stable state; otherwise, the vehicle is in an unstable state.(4)γ-γ∗≤γth,β+B1β˙-B2≤0,where γth is the threshold of yaw rate determined by the change rate of steering angle.In order to prevent the control system from frequent operation, the yaw stability control module is activated only when the deviation reaches a certain value. The control threshold is subtracted to determine effective deviation of additional yaw moment.(5)eγ=Δγ-Δγ+Δγ>Δγ+Δγ-Δγ-Δγ<Δγ-0others,eβ=Δβ-Δβ+Δβ>Δβ+Δβ-Δβ-Δβ<Δβ-0others.


### 2.2. Rollover Stability Control Analysis

When the off-road vehicle with high centroid runs in the high adhesion road with high speed, the emergency avoidance operation can easily cause the risk of rollover because of large lateral acceleration. The tire load transfer rate is the most direct indicator to describe vehicle rollover condition, which can be expressed by(6)LTR=Fzo-FziFzo+Fzi,where Fzo, Fzi are outside and inside wheel vertical load.The vehicle rollover state identified by the LTR value reflects the current vehicle state. For the driver, the vehicle rollover is possible to happen before the controller identifies the wheel off the ground by calculating LTR = 1. In the vehicle rollover stability analysis, it is important to predict the rollover state in the coming time.The system warns the driver in advance to take an efficient operation within sufficient time or enables the rollover control system to take appropriate action by a command, effectively improving the rollover stability of vehicle [19].In 2001, Bo-Chiuan Chen from University of Michigan firstly proposed the concept of Time-To-Rollover (TTR), which characterizes the time interval from the current moment to the rollover time [15]. TTR is a good index to measure the risk degree of vehicle rollover. The TTR is calculated as follows. The vehicle prediction model can predict vehicle movement state within the longest N steps (Ts represents a step) according to the current vehicle state parameter and the steering wheel angle input. If LTR does not exceed its threshold within N steps, then TTR=TTRth (TTRth represents the threshold), indicating that there is no risk of rollover in prediction time; if LTR exceeds the threshold in the nth step (n<N), then TTR=n·Ts, indicating that there is a certain risk of rollover. The closer the TTR is to 0, the greater the risk of rollover is. At this moment, the prediction model stops the forecast and takes TTR value as the output. Meanwhile, the rollover warning system issues a warning and commands the rollover control system to intervene. Figure 2 shows the specific process.Figure 2
TTR rollover warning algorithm.Wherein, the rollover prediction model takes the current vehicle speed, the steering wheel angle, the additional yaw moment calculated at the previous time, and the vehicle sensor measurement value or the vehicle state estimation value as input. During the prediction process, it is assumed that the vehicle speed, steering wheel angle, and additional yaw moment remain unchanged.The problem of vehicle yaw stability can occur in the road surface with any adhesion condition, and the rollover stability problem mainly occurs in the high adhesion road. Therefore, the yaw stability control takes effect in low adhesion road. In the high adhesion road, both yaw and rollover stability controls are possibly triggered synchronously or asynchronously.


### 2.3. Yaw Moment Decision Controller Based on MPC

The integrated controller determines the total driving/braking torque demand based on the throttle/brake pedal opening or the difference between the target and actual vehicle speed. YSC yaw moment decision module calculates the required yaw moment ensuring vehicle yaw stability by taking the difference between the actual and ideal values of centroid slip/yaw angle velocity as the input. RSC yaw moment decision module calculates the required yaw moment for rollover stabilization control by taking the LTR value output by the rollover prediction system as the input. The integrated control module of yaw/rollover stability finally determines the total yaw moment based on yaw and rollover state. When the sum of four-wheel target driving/braking torque is equal to the total demand torque, then the total demand torque is assigned to each wheel according to a certain distribution rule. Thus, the longitudinal forces of wheels produce the desired yaw moments. Wherein, the model prediction theory is used for YSC and RSC yaw moment decision modules.


#### 2.3.1. Vehicle Reference Model

With little computation and good real-time control, the linear two-DOF vehicle model can better describe the driver’s driving intention. The two degrees of freedom contain lateral and yaw directions. Figure3 shows the specific reference vehicle model.Figure 3
Linear two-DOF vehicle model.In Figure3, β is the side slip angle; γ is the yaw rate; δ is the front wheel rotor angle; V is the centroid velocity; α1 and α2 are the slip angles of front and rear wheels; vx and vy are the components of centroid velocity in axes x and y; u1 and u2 are the velocities at the midpoints of front and rear axles; Fy1 and Fy2 are the slip forces of front and rear wheels.(7)Fy1=k1α1,Fy2=k2α2,where k1 and k2 are the equivalent cornering stiffness of front and rear axles. According to the provisions of the tire coordinate system, the slip angles of front and rear wheels are expressed as follows:(8)α1=β+aγvx-δ,α2=β-bγvx,where a, b are the distance from centroid to front and rear axle.The differential equations of motion in lateral and yaw directions are expressed as follows:(9)k1+k2β+ak1-bk2vxγ-k1δ=mv˙y+vxγ,ak1-bk2β+a2k1+b2k2vxγ-ak1δ=Izγ˙.In stable state, the yaw angle velocity of vehicle is fixed. Here, γ˙=0; vy=0. Equation (9) is used to obtain the reference values of centroid slip and yaw angle velocities.(10)βd=b1+Kvx2L+amvx21+Kvx2k2L2δ,γd=vx/L1+Kvx2δ,where L is the wheelbase, K is the stability factor, and K=ma/k2-b/k1/L2.It is assumed that the tire of linear two-DOF vehicle model is linear. When the vehicle is in the limit state, the steady-state response of reference model is not suitable as a reference value. Therefore, the reference value is replaced by boundary value [20]. Equations (11) and (12) show the boundary values of the side slip angle and yaw rate.The boundary values of the side slip angle are as follows:(11)βd_b=arctan⁡0.02μg.The boundary values of the yaw rate are as follows:(12)γd_b=0.85μgvx,where μ is the pavement friction coefficient.


#### 2.3.2. Vehicle Prediction Model

In Figure4, m is the mass of the vehicle; ms is the sprung mass; vx, vy are the component of the centroid velocity on the x- and y-axes; γ is the yaw rate; ϕv is the roll angle of the vehicle; ϕr is the lateral slope of the road; δ is the front wheel angle; ay is the lateral acceleration of the vehicle; e is the distance from the center of mass to the roll axis; Kϕ is the roll stiffness of the suspension; Cϕ is the roll damping of the suspension; Ix, Iz are the rotational inertia of the axle around the x- and z-axes of the vehicle; Ixeq is the moment of inertia of the sprung mass around the roll axis: Ixeq=Ix+mse2.Figure 4
Three-DOF vehicle model.In general, the more accurate prediction model leads to better control effect of the controller. However, the precise model weakens the real-time performance of controller. At present, most controllers in real vehicle use two or three-DOF model [21]. Therefore, we established a three-DOF vehicle model with additional yaw moment Mz for prediction (see Figure 4). The dynamic equations in lateral, yaw, and roll directions are expressed as follows:(13)mvxβ˙+γ-mseϕ¨=Fy1+Fy2,Izγ˙=aFy1-bFy2+Mz,Ixeqϕ¨=msegϕ+ay-Kϕϕ-Cϕϕ˙.According toϕ˙ = ϕ˙, (7), (8), and (13), the differential equation of three-DOF vehicle model can be rewritten into a matrix form:(14)mvx00-mse0Iz00-msevx00Ixeq0010β˙γ˙ϕ˙ϕ¨=CDvx-mvx00DEvx000msevxmseg-Kϕ-Cϕ0001βγϕϕ˙+0100Mz+-k1-ak100δ,where C=k1+k2;D=ak1-bk2;E=a2k1+b2k2. xt=βγϕϕ˙T is selected as the state vector; u(t)=Mz is the control input; w(t)=δ is the disturbance variable. Equation (14) can also be expressed as(15)Fx˙t=Gxt+Put+Qwt.The nonlinear equation (14) is linearized by Taylor series expansion and simplifies the system design by reducing the higher-order derivative as a small amount. Equation (16) is the state space equation after linearization.(16)x˙t=Axt+B1ut+B2wt,yt=Cxt,where A=F-1G;B1=F-1P;B2=F-1Q;C=I4. I4 is the four-order unit matrix.The two-stage fold line tire model is used for vehicle prediction (see Figure5). When the tire works in the linear region, the cornering stiffness can be linear. After the tire enters the nonlinear region, the side angle increases. However, the concerning force remains unchanged; namely, the slope of the curve Fy-α is zero.Figure 5
Two-stage fold line tire model.


#### 2.3.3. Design of MPC Control Algorithm

MPC belongs to a discrete control. The state equation of (16) cannot be used directly as a prediction model. It needs to be transformed into incremental form by discretization.First, (16) can be transformed into(17)xt+1-xtΔt=Axt+B1ut+B2wt,yt=Cxt.Then,(18)xt+1=Acxt+B1cut+B2cwt,yt=Cxt,where Δt is sampling time, Ac=Δt·A+I4;B1c=Δt·B1;B2c=Δt·B2Lett=k, k is sampling point. Then get(19)xk+1=Acxk+B1cuk+B2cwk,yk=Cxk.Ac, B1c, and B2c are discrete matrices, in order to reduce or eliminate static errors; (19) is rewritten as an incremental model (20).(20)Δxk+1=AcΔxk+B1cΔuk+B2cΔwk,yk=CΔxk+yk-1,where Δx(k)=x(k)-x(k-1);Δu(k)=u(k)-u(k-1);Δw(k)=w(k)-w(k-1).According to the MPC theory, the most recent measured value is taken as the initial condition to predict the future dynamic based on the predictive model. We make the following assumptions in the work:(
1
) The predicted time domain is p; the control time domain is m; and p≥m (m=3, p=10).(
2
) At the current k moment, the estimated and measured values (x^(k) and y(k)) of state vector are known.(
3
) The controlled variable remains unchanged outside of control time domain. Δu(k+i)=0(i≥m).(
4
) The measurable disturbance variable remains unchanged after moment k. Δw(k+i)=0(i≥1).Within predicted time domainp, the prediction expression for the pth step state vector at moment k is as follows:(21)Δxk+1∣k=AcΔx^k+B1cΔuk+B2cΔwkΔxk+2∣k=Ac2Δx^k+AcB1cΔuk+B1cΔuk+1+AcB2cΔwk⋮Δxk+p∣k=AcpΔxk+Acp-1B1cΔuk+Acp-2B1cΔuk+1+⋯+B1cΔuk+p-1+Acp-1B2cΔwk.At moment k, we predict the output from k+1 to k+p according to the predicted value of state vector in (22)yk∣k=ykyk+1∣k=CAcΔx^k+CB1cΔuk+CB2cΔwk+ykyk+2∣k=CAc2+CAcΔx^k+CAcB1c+CB1cΔuk+CB1cΔuk+1+CAcB2c+CB2cΔwk+ykyk+3∣k=∑i=13CAciΔx^k+∑i=13CAci-1B1cΔuk+⋯+CB1cΔuk+2+∑i=13CAci-1B2cΔuk+yk⋮yk+p∣k=∑i=1pCAciΔx^k+∑i=1pCAci-1B1cΔuk+⋯+CB1cΔuk+p-1+∑i=1pCAci-1B2cΔuk+yk.The vectors Yp(k) and ΔUm(k) are defined as follows:(23)Ypk=yk+1∣kyk+2∣k⋯yk+p∣kT,ΔUmk=ΔukΔuk+1⋯Δuk+m-1T.The pth step control output of the prediction system can be expressed as(24)Ypk=SxΔx^k+Lyyk+SuΔUk+SwΔwk,where(25)Sx=CAcCAc2+CAc⋮∑i=1pCAci;Sw=CB2cCAcB2c+CB2c⋮∑i=1pCAci-1B2c;Lc=I4I4⋮I4,Su=CB1c00⋯0CAcB1cCB1c0⋯0⋮⋮⋮⋯⋮∑i=1pCAciB1c∑i=1p-1CAci-1B1c⋯⋯∑i=1p-m+1CAci-1B1c.Similarly, we can deduce thepth step predictive equation with constraint output.MPC is an algorithm of determining control strategy by optimization. Firstly, it is necessary to determinem control increments Δu(k),Δu(k+1),…,Δu(k+m-1) from moment k. In the next p moment, the predicted output values y(k+1∣k),y(k+2∣k),…,y(k+p∣k) of the controlled object can approximate to the given expected value d(k+i),i=1,2,…,p.In the control process, the control increment should not have sharp change. Based on this factor, soft constraints are added to optimization performance indexes. Thus, the optimization performance index of momentk can be expressed as(26)minJk=∑i=1pqidk+i-yk+1+∑j=1mriΔu2k+j-1,where qi and ri are the weighting coefficients, which represent error tracking and suppression of control increment. If the constraint is not considered, the above problem is described as follows.TakingΔUk=ΔukΔuk+1⋯Δuk+m-1T as optimal variable, the performance index (26) is minimized under prediction model (20). Likewise, optimal performance index (26) can also be written as the form of vector.(27)minJk=Dpk-YpkQ2+ΔUmkR2,where Dpk=d(k+1)d(k+2)⋯d(k+p)T is the expected output value in the optimization time domain, namely, side slip angle and yaw rate based on vehicle model output. It is assumed that the expected value is fixed in the optimization time domain; namely, d(k+1)=d(k+2)=⋯=d(k+i);Q=diagq1q2⋯qp; and R=diagr1r2⋯rm are diagonal matrices comprised of weighting coefficients.At the momentk, Dp(k), and Yp(k) are known. In order to minimize ΔUm(k) by J(k), the optimal control sequence is obtained according to the necessary condition of the extremum: dJ(k)/dΔUm(k)=0.(28)ΔU∗k=SuTQTQSu+QTQ-1SuTQTQEpk+1∣k,where Ep(k+1∣k) is control deviation, which can be calculated on line by the following equation:(29)Epk+1∣k=Dpk-SxΔx^k-Lyyk-SwΔwk.In each sampling period, the first element of the derived optimal solution ΔU∗(k) is the control increment of MPC.(30)Δuk=100ΔU∗k,where u(k)=u(k-1)+Δu(k) is the additional yaw moment provided for the torque distribution layer.The Lyapunov functionV∗(k) at the moment k is defined as the optimal performance index of open loop.(31)V∗k=defJ∗k.Relationship between V∗(k) and V(k+1) is established:(32)V∗k=Vk+1+Jy∗1∣k,u∗0∣k,Δu∗0∣k-Jyp∣k+1,up-1∣k+1,Δup-1∣k+1.Since V(k+1)≥V∗(k+1), we can derive(33)V∗k+1-V∗k≤y∗p+1∣k22-y∗1∣k22+u∗p-1∣k22-u∗0∣k22-Δu∗0∣k22. As y(p+1∣k)22≤y(1∣k)22 and u(p∣k)22≤u(0∣k)22, we can get(34)V∗k+1-V∗k≤0.According to the Lyapunov second stability theorem, the proposed MPC closed-loop system is stable.


### 2.4. Wheel Torque Distribution

For the tire, lower pavement adhesion consumption rate leads to greater adhesion margin and larger distance between the tire and nonlinear saturation zone, indirectly improving the stability of the vehicle. According to tire friction ellipse theory, the work distributes the wheel torque for the purpose of minimum sum of pavement adhesion consumption rates of four tires. Equation (35) shows the objective function.(35)minJμ=∑i=14CiFxi2+Fyi2μ2Fzi2,i=1,2,3,4,where Fxi is the tire longitudinal force; Fyi the tire lateral force; Fzi the vertical load on the wheel; Ci the weight coefficient within the interval (0,1). The tire longitudinal force Fxi can be obtained by (36)Fxi=Twi-Ieqω˙irw,where Twi is the drive torque on the wheel; rw the wheel radius; ω˙i the wheel rotation angle acceleration; Ieq the wheel equivalent inertia of in-wheel motor and reducer inertia.The vertical loadFzi of wheel consists of static load and dynamic load caused by longitudinal and lateral accelerations without considering the centroid change caused by vehicle rollover and the influence of suspension system. The vertical load of each wheel is calculated by (37)Fz1Fz2Fz3Fz4=mg2Lbbaa+maxh2L-1-111+mayhBL-bb-aa,where a and b are the distance from the centroid to front and rear axle; L is the wheelbase; h is the height of centroid; B is the trend.The objective function in (35) determines smallest sum of the load rates of four tires. However, it is impossible to ensure that each tire works at a low load-rate state. Therefore, the utilization factor of each tire is adjusted by the weighting factor Ci. Under different driving conditions, Ci is not a fixed value. In simulation and debugging process, it is found that the value has great relationship with tire vertical load and road adhesion condition. The following equation is selected in the work:(38)Ci=μFzimg2,i=1,2,3,4.Here, we focus on longitudinal driving torque distribution in the conventional operating conditions. It is difficult to acquire the lateral force of tire in actual driving, so the pavement adhesion consumption rate is simplified to the longitudinal adhesion consumption rate. Equation (35) can be transformed into the following equation:(39)minJμ=∑i=14CiFxiμFzi2.When the yaw moment control is not required, the sum of longitudinal forces of four wheels satisfies the total moment demand, which is obtained by analyzing the pedal input. Therefore, constraint equation is expressed as follows:(40)∑i=14Fxi=fθp,where θp is the accelerator pedal opening; f(~) the pedal analysis strategy. In addition, the longitudinal force of tire is also restricted by pavement adhesion condition and motor output capacity.(41)-μFzi≤Fxi≤min⁡μFzi,Fmax,i=1,2,3,4,where Fmax is the maximum wheel drive force provided by the motor; namely, Fmax=Twmax/rw. Therefore, the optimization problem can be expressed as follows:(42)minJμ=∑i=14CiFxiμFzi2s.t.∑i=14Fxi=fθp-μFzi≤Fxi≤min⁡μFzi,FmaxFx2-Fx1B2+Fx4-Fx3B2=Tz.The nonlinear constrained optimization problem is solved by numerical optimization method to obtain the longitudinal forces for four tires. Then, we derive the drive torques which should be allocated to four wheels by calculation.


## 2.1. Yaw Stability Control Analysis

The yaw stability is characterized by state variables including side slip angle and yaw rate. Therefore, we determine the yaw stability of vehicle, and whether to carry out yaw stability control according to the two state variables.When the side slip angle of vehicle satisfies the following equation, the vehicle is in a stable zone; otherwise, it is in an unstable zone, and the stability control must be added.(1)β+B1β˙≤B2,where B1 and B2 are the parameters relevant to vehicle characteristics and driving state. In the work, we refer to the method of determining the boundary values of side slip angle by vehicle speed and pavement friction coefficient in [17]. Then, the relation between B1 and B2 is obtained by fitting.(2)β˙+B1β≤B2,B1μ=15.62μ2-34.37μ-6.719,B2V,μ=0.0002343μ2V2-0.000516μV2-0.7498μ2+1.650μ,where μ is the tire-road friction coefficient and V is the travel speed. Equation (2) is used to determine the stability boundary of side slip angle. After that, we obtain the estimated values of side slip angle and its derivative through the state estimation procedure which has been developed and validated by the authors in [18], thus deriving actual state point of the vehicle. The driving state of vehicle can be determined and expressed as (3) according to the distance between actual state point and stable boundary.(3)ΔD=β+B1β˙-B2.According to (1)–(3), we obtain the vehicle stability judgment phase diagram based on the side slip angle (see Figure 1).Figure 1
Vehicle stability judgment phase diagram based on the side slip angle.If the side slip angle and yaw rate of the vehicle satisfy the following equation, then the vehicle is in a stable state; otherwise, the vehicle is in an unstable state.(4)γ-γ∗≤γth,β+B1β˙-B2≤0,where γth is the threshold of yaw rate determined by the change rate of steering angle.In order to prevent the control system from frequent operation, the yaw stability control module is activated only when the deviation reaches a certain value. The control threshold is subtracted to determine effective deviation of additional yaw moment.(5)eγ=Δγ-Δγ+Δγ>Δγ+Δγ-Δγ-Δγ<Δγ-0others,eβ=Δβ-Δβ+Δβ>Δβ+Δβ-Δβ-Δβ<Δβ-0others.


## 2.2. Rollover Stability Control Analysis

When the off-road vehicle with high centroid runs in the high adhesion road with high speed, the emergency avoidance operation can easily cause the risk of rollover because of large lateral acceleration. The tire load transfer rate is the most direct indicator to describe vehicle rollover condition, which can be expressed by(6)LTR=Fzo-FziFzo+Fzi,where Fzo, Fzi are outside and inside wheel vertical load.The vehicle rollover state identified by the LTR value reflects the current vehicle state. For the driver, the vehicle rollover is possible to happen before the controller identifies the wheel off the ground by calculating LTR = 1. In the vehicle rollover stability analysis, it is important to predict the rollover state in the coming time.The system warns the driver in advance to take an efficient operation within sufficient time or enables the rollover control system to take appropriate action by a command, effectively improving the rollover stability of vehicle [19].In 2001, Bo-Chiuan Chen from University of Michigan firstly proposed the concept of Time-To-Rollover (TTR), which characterizes the time interval from the current moment to the rollover time [15]. TTR is a good index to measure the risk degree of vehicle rollover. The TTR is calculated as follows. The vehicle prediction model can predict vehicle movement state within the longest N steps (Ts represents a step) according to the current vehicle state parameter and the steering wheel angle input. If LTR does not exceed its threshold within N steps, then TTR=TTRth (TTRth represents the threshold), indicating that there is no risk of rollover in prediction time; if LTR exceeds the threshold in the nth step (n<N), then TTR=n·Ts, indicating that there is a certain risk of rollover. The closer the TTR is to 0, the greater the risk of rollover is. At this moment, the prediction model stops the forecast and takes TTR value as the output. Meanwhile, the rollover warning system issues a warning and commands the rollover control system to intervene. Figure 2 shows the specific process.Figure 2
TTR rollover warning algorithm.Wherein, the rollover prediction model takes the current vehicle speed, the steering wheel angle, the additional yaw moment calculated at the previous time, and the vehicle sensor measurement value or the vehicle state estimation value as input. During the prediction process, it is assumed that the vehicle speed, steering wheel angle, and additional yaw moment remain unchanged.The problem of vehicle yaw stability can occur in the road surface with any adhesion condition, and the rollover stability problem mainly occurs in the high adhesion road. Therefore, the yaw stability control takes effect in low adhesion road. In the high adhesion road, both yaw and rollover stability controls are possibly triggered synchronously or asynchronously.


## 2.3. Yaw Moment Decision Controller Based on MPC

The integrated controller determines the total driving/braking torque demand based on the throttle/brake pedal opening or the difference between the target and actual vehicle speed. YSC yaw moment decision module calculates the required yaw moment ensuring vehicle yaw stability by taking the difference between the actual and ideal values of centroid slip/yaw angle velocity as the input. RSC yaw moment decision module calculates the required yaw moment for rollover stabilization control by taking the LTR value output by the rollover prediction system as the input. The integrated control module of yaw/rollover stability finally determines the total yaw moment based on yaw and rollover state. When the sum of four-wheel target driving/braking torque is equal to the total demand torque, then the total demand torque is assigned to each wheel according to a certain distribution rule. Thus, the longitudinal forces of wheels produce the desired yaw moments. Wherein, the model prediction theory is used for YSC and RSC yaw moment decision modules.


### 2.3.1. Vehicle Reference Model

With little computation and good real-time control, the linear two-DOF vehicle model can better describe the driver’s driving intention. The two degrees of freedom contain lateral and yaw directions. Figure3 shows the specific reference vehicle model.Figure 3
Linear two-DOF vehicle model.In Figure3, β is the side slip angle; γ is the yaw rate; δ is the front wheel rotor angle; V is the centroid velocity; α1 and α2 are the slip angles of front and rear wheels; vx and vy are the components of centroid velocity in axes x and y; u1 and u2 are the velocities at the midpoints of front and rear axles; Fy1 and Fy2 are the slip forces of front and rear wheels.(7)Fy1=k1α1,Fy2=k2α2,where k1 and k2 are the equivalent cornering stiffness of front and rear axles. According to the provisions of the tire coordinate system, the slip angles of front and rear wheels are expressed as follows:(8)α1=β+aγvx-δ,α2=β-bγvx,where a, b are the distance from centroid to front and rear axle.The differential equations of motion in lateral and yaw directions are expressed as follows:(9)k1+k2β+ak1-bk2vxγ-k1δ=mv˙y+vxγ,ak1-bk2β+a2k1+b2k2vxγ-ak1δ=Izγ˙.In stable state, the yaw angle velocity of vehicle is fixed. Here, γ˙=0; vy=0. Equation (9) is used to obtain the reference values of centroid slip and yaw angle velocities.(10)βd=b1+Kvx2L+amvx21+Kvx2k2L2δ,γd=vx/L1+Kvx2δ,where L is the wheelbase, K is the stability factor, and K=ma/k2-b/k1/L2.It is assumed that the tire of linear two-DOF vehicle model is linear. When the vehicle is in the limit state, the steady-state response of reference model is not suitable as a reference value. Therefore, the reference value is replaced by boundary value [20]. Equations (11) and (12) show the boundary values of the side slip angle and yaw rate.The boundary values of the side slip angle are as follows:(11)βd_b=arctan⁡0.02μg.The boundary values of the yaw rate are as follows:(12)γd_b=0.85μgvx,where μ is the pavement friction coefficient.


### 2.3.2. Vehicle Prediction Model

In Figure4, m is the mass of the vehicle; ms is the sprung mass; vx, vy are the component of the centroid velocity on the x- and y-axes; γ is the yaw rate; ϕv is the roll angle of the vehicle; ϕr is the lateral slope of the road; δ is the front wheel angle; ay is the lateral acceleration of the vehicle; e is the distance from the center of mass to the roll axis; Kϕ is the roll stiffness of the suspension; Cϕ is the roll damping of the suspension; Ix, Iz are the rotational inertia of the axle around the x- and z-axes of the vehicle; Ixeq is the moment of inertia of the sprung mass around the roll axis: Ixeq=Ix+mse2.Figure 4
Three-DOF vehicle model.In general, the more accurate prediction model leads to better control effect of the controller. However, the precise model weakens the real-time performance of controller. At present, most controllers in real vehicle use two or three-DOF model [21]. Therefore, we established a three-DOF vehicle model with additional yaw moment Mz for prediction (see Figure 4). The dynamic equations in lateral, yaw, and roll directions are expressed as follows:(13)mvxβ˙+γ-mseϕ¨=Fy1+Fy2,Izγ˙=aFy1-bFy2+Mz,Ixeqϕ¨=msegϕ+ay-Kϕϕ-Cϕϕ˙.According toϕ˙ = ϕ˙, (7), (8), and (13), the differential equation of three-DOF vehicle model can be rewritten into a matrix form:(14)mvx00-mse0Iz00-msevx00Ixeq0010β˙γ˙ϕ˙ϕ¨=CDvx-mvx00DEvx000msevxmseg-Kϕ-Cϕ0001βγϕϕ˙+0100Mz+-k1-ak100δ,where C=k1+k2;D=ak1-bk2;E=a2k1+b2k2. xt=βγϕϕ˙T is selected as the state vector; u(t)=Mz is the control input; w(t)=δ is the disturbance variable. Equation (14) can also be expressed as(15)Fx˙t=Gxt+Put+Qwt.The nonlinear equation (14) is linearized by Taylor series expansion and simplifies the system design by reducing the higher-order derivative as a small amount. Equation (16) is the state space equation after linearization.(16)x˙t=Axt+B1ut+B2wt,yt=Cxt,where A=F-1G;B1=F-1P;B2=F-1Q;C=I4. I4 is the four-order unit matrix.The two-stage fold line tire model is used for vehicle prediction (see Figure5). When the tire works in the linear region, the cornering stiffness can be linear. After the tire enters the nonlinear region, the side angle increases. However, the concerning force remains unchanged; namely, the slope of the curve Fy-α is zero.Figure 5
Two-stage fold line tire model.


### 2.3.3. Design of MPC Control Algorithm

MPC belongs to a discrete control. The state equation of (16) cannot be used directly as a prediction model. It needs to be transformed into incremental form by discretization.First, (16) can be transformed into(17)xt+1-xtΔt=Axt+B1ut+B2wt,yt=Cxt.Then,(18)xt+1=Acxt+B1cut+B2cwt,yt=Cxt,where Δt is sampling time, Ac=Δt·A+I4;B1c=Δt·B1;B2c=Δt·B2Lett=k, k is sampling point. Then get(19)xk+1=Acxk+B1cuk+B2cwk,yk=Cxk.Ac, B1c, and B2c are discrete matrices, in order to reduce or eliminate static errors; (19) is rewritten as an incremental model (20).(20)Δxk+1=AcΔxk+B1cΔuk+B2cΔwk,yk=CΔxk+yk-1,where Δx(k)=x(k)-x(k-1);Δu(k)=u(k)-u(k-1);Δw(k)=w(k)-w(k-1).According to the MPC theory, the most recent measured value is taken as the initial condition to predict the future dynamic based on the predictive model. We make the following assumptions in the work:(
1
) The predicted time domain is p; the control time domain is m; and p≥m (m=3, p=10).(
2
) At the current k moment, the estimated and measured values (x^(k) and y(k)) of state vector are known.(
3
) The controlled variable remains unchanged outside of control time domain. Δu(k+i)=0(i≥m).(
4
) The measurable disturbance variable remains unchanged after moment k. Δw(k+i)=0(i≥1).Within predicted time domainp, the prediction expression for the pth step state vector at moment k is as follows:(21)Δxk+1∣k=AcΔx^k+B1cΔuk+B2cΔwkΔxk+2∣k=Ac2Δx^k+AcB1cΔuk+B1cΔuk+1+AcB2cΔwk⋮Δxk+p∣k=AcpΔxk+Acp-1B1cΔuk+Acp-2B1cΔuk+1+⋯+B1cΔuk+p-1+Acp-1B2cΔwk.At moment k, we predict the output from k+1 to k+p according to the predicted value of state vector in (22)yk∣k=ykyk+1∣k=CAcΔx^k+CB1cΔuk+CB2cΔwk+ykyk+2∣k=CAc2+CAcΔx^k+CAcB1c+CB1cΔuk+CB1cΔuk+1+CAcB2c+CB2cΔwk+ykyk+3∣k=∑i=13CAciΔx^k+∑i=13CAci-1B1cΔuk+⋯+CB1cΔuk+2+∑i=13CAci-1B2cΔuk+yk⋮yk+p∣k=∑i=1pCAciΔx^k+∑i=1pCAci-1B1cΔuk+⋯+CB1cΔuk+p-1+∑i=1pCAci-1B2cΔuk+yk.The vectors Yp(k) and ΔUm(k) are defined as follows:(23)Ypk=yk+1∣kyk+2∣k⋯yk+p∣kT,ΔUmk=ΔukΔuk+1⋯Δuk+m-1T.The pth step control output of the prediction system can be expressed as(24)Ypk=SxΔx^k+Lyyk+SuΔUk+SwΔwk,where(25)Sx=CAcCAc2+CAc⋮∑i=1pCAci;Sw=CB2cCAcB2c+CB2c⋮∑i=1pCAci-1B2c;Lc=I4I4⋮I4,Su=CB1c00⋯0CAcB1cCB1c0⋯0⋮⋮⋮⋯⋮∑i=1pCAciB1c∑i=1p-1CAci-1B1c⋯⋯∑i=1p-m+1CAci-1B1c.Similarly, we can deduce thepth step predictive equation with constraint output.MPC is an algorithm of determining control strategy by optimization. Firstly, it is necessary to determinem control increments Δu(k),Δu(k+1),…,Δu(k+m-1) from moment k. In the next p moment, the predicted output values y(k+1∣k),y(k+2∣k),…,y(k+p∣k) of the controlled object can approximate to the given expected value d(k+i),i=1,2,…,p.In the control process, the control increment should not have sharp change. Based on this factor, soft constraints are added to optimization performance indexes. Thus, the optimization performance index of momentk can be expressed as(26)minJk=∑i=1pqidk+i-yk+1+∑j=1mriΔu2k+j-1,where qi and ri are the weighting coefficients, which represent error tracking and suppression of control increment. If the constraint is not considered, the above problem is described as follows.TakingΔUk=ΔukΔuk+1⋯Δuk+m-1T as optimal variable, the performance index (26) is minimized under prediction model (20). Likewise, optimal performance index (26) can also be written as the form of vector.(27)minJk=Dpk-YpkQ2+ΔUmkR2,where Dpk=d(k+1)d(k+2)⋯d(k+p)T is the expected output value in the optimization time domain, namely, side slip angle and yaw rate based on vehicle model output. It is assumed that the expected value is fixed in the optimization time domain; namely, d(k+1)=d(k+2)=⋯=d(k+i);Q=diagq1q2⋯qp; and R=diagr1r2⋯rm are diagonal matrices comprised of weighting coefficients.At the momentk, Dp(k), and Yp(k) are known. In order to minimize ΔUm(k) by J(k), the optimal control sequence is obtained according to the necessary condition of the extremum: dJ(k)/dΔUm(k)=0.(28)ΔU∗k=SuTQTQSu+QTQ-1SuTQTQEpk+1∣k,where Ep(k+1∣k) is control deviation, which can be calculated on line by the following equation:(29)Epk+1∣k=Dpk-SxΔx^k-Lyyk-SwΔwk.In each sampling period, the first element of the derived optimal solution ΔU∗(k) is the control increment of MPC.(30)Δuk=100ΔU∗k,where u(k)=u(k-1)+Δu(k) is the additional yaw moment provided for the torque distribution layer.The Lyapunov functionV∗(k) at the moment k is defined as the optimal performance index of open loop.(31)V∗k=defJ∗k.Relationship between V∗(k) and V(k+1) is established:(32)V∗k=Vk+1+Jy∗1∣k,u∗0∣k,Δu∗0∣k-Jyp∣k+1,up-1∣k+1,Δup-1∣k+1.Since V(k+1)≥V∗(k+1), we can derive(33)V∗k+1-V∗k≤y∗p+1∣k22-y∗1∣k22+u∗p-1∣k22-u∗0∣k22-Δu∗0∣k22. As y(p+1∣k)22≤y(1∣k)22 and u(p∣k)22≤u(0∣k)22, we can get(34)V∗k+1-V∗k≤0.According to the Lyapunov second stability theorem, the proposed MPC closed-loop system is stable.


## 2.3.1. Vehicle Reference Model

With little computation and good real-time control, the linear two-DOF vehicle model can better describe the driver’s driving intention. The two degrees of freedom contain lateral and yaw directions. Figure3 shows the specific reference vehicle model.Figure 3
Linear two-DOF vehicle model.In Figure3, β is the side slip angle; γ is the yaw rate; δ is the front wheel rotor angle; V is the centroid velocity; α1 and α2 are the slip angles of front and rear wheels; vx and vy are the components of centroid velocity in axes x and y; u1 and u2 are the velocities at the midpoints of front and rear axles; Fy1 and Fy2 are the slip forces of front and rear wheels.(7)Fy1=k1α1,Fy2=k2α2,where k1 and k2 are the equivalent cornering stiffness of front and rear axles. According to the provisions of the tire coordinate system, the slip angles of front and rear wheels are expressed as follows:(8)α1=β+aγvx-δ,α2=β-bγvx,where a, b are the distance from centroid to front and rear axle.The differential equations of motion in lateral and yaw directions are expressed as follows:(9)k1+k2β+ak1-bk2vxγ-k1δ=mv˙y+vxγ,ak1-bk2β+a2k1+b2k2vxγ-ak1δ=Izγ˙.In stable state, the yaw angle velocity of vehicle is fixed. Here, γ˙=0; vy=0. Equation (9) is used to obtain the reference values of centroid slip and yaw angle velocities.(10)βd=b1+Kvx2L+amvx21+Kvx2k2L2δ,γd=vx/L1+Kvx2δ,where L is the wheelbase, K is the stability factor, and K=ma/k2-b/k1/L2.It is assumed that the tire of linear two-DOF vehicle model is linear. When the vehicle is in the limit state, the steady-state response of reference model is not suitable as a reference value. Therefore, the reference value is replaced by boundary value [20]. Equations (11) and (12) show the boundary values of the side slip angle and yaw rate.The boundary values of the side slip angle are as follows:(11)βd_b=arctan⁡0.02μg.The boundary values of the yaw rate are as follows:(12)γd_b=0.85μgvx,where μ is the pavement friction coefficient.


## 2.3.2. Vehicle Prediction Model

In Figure4, m is the mass of the vehicle; ms is the sprung mass; vx, vy are the component of the centroid velocity on the x- and y-axes; γ is the yaw rate; ϕv is the roll angle of the vehicle; ϕr is the lateral slope of the road; δ is the front wheel angle; ay is the lateral acceleration of the vehicle; e is the distance from the center of mass to the roll axis; Kϕ is the roll stiffness of the suspension; Cϕ is the roll damping of the suspension; Ix, Iz are the rotational inertia of the axle around the x- and z-axes of the vehicle; Ixeq is the moment of inertia of the sprung mass around the roll axis: Ixeq=Ix+mse2.Figure 4
Three-DOF vehicle model.In general, the more accurate prediction model leads to better control effect of the controller. However, the precise model weakens the real-time performance of controller. At present, most controllers in real vehicle use two or three-DOF model [21]. Therefore, we established a three-DOF vehicle model with additional yaw moment Mz for prediction (see Figure 4). The dynamic equations in lateral, yaw, and roll directions are expressed as follows:(13)mvxβ˙+γ-mseϕ¨=Fy1+Fy2,Izγ˙=aFy1-bFy2+Mz,Ixeqϕ¨=msegϕ+ay-Kϕϕ-Cϕϕ˙.According toϕ˙ = ϕ˙, (7), (8), and (13), the differential equation of three-DOF vehicle model can be rewritten into a matrix form:(14)mvx00-mse0Iz00-msevx00Ixeq0010β˙γ˙ϕ˙ϕ¨=CDvx-mvx00DEvx000msevxmseg-Kϕ-Cϕ0001βγϕϕ˙+0100Mz+-k1-ak100δ,where C=k1+k2;D=ak1-bk2;E=a2k1+b2k2. xt=βγϕϕ˙T is selected as the state vector; u(t)=Mz is the control input; w(t)=δ is the disturbance variable. Equation (14) can also be expressed as(15)Fx˙t=Gxt+Put+Qwt.The nonlinear equation (14) is linearized by Taylor series expansion and simplifies the system design by reducing the higher-order derivative as a small amount. Equation (16) is the state space equation after linearization.(16)x˙t=Axt+B1ut+B2wt,yt=Cxt,where A=F-1G;B1=F-1P;B2=F-1Q;C=I4. I4 is the four-order unit matrix.The two-stage fold line tire model is used for vehicle prediction (see Figure5). When the tire works in the linear region, the cornering stiffness can be linear. After the tire enters the nonlinear region, the side angle increases. However, the concerning force remains unchanged; namely, the slope of the curve Fy-α is zero.Figure 5
Two-stage fold line tire model.


## 2.3.3. Design of MPC Control Algorithm

MPC belongs to a discrete control. The state equation of (16) cannot be used directly as a prediction model. It needs to be transformed into incremental form by discretization.First, (16) can be transformed into(17)xt+1-xtΔt=Axt+B1ut+B2wt,yt=Cxt.Then,(18)xt+1=Acxt+B1cut+B2cwt,yt=Cxt,where Δt is sampling time, Ac=Δt·A+I4;B1c=Δt·B1;B2c=Δt·B2Lett=k, k is sampling point. Then get(19)xk+1=Acxk+B1cuk+B2cwk,yk=Cxk.Ac, B1c, and B2c are discrete matrices, in order to reduce or eliminate static errors; (19) is rewritten as an incremental model (20).(20)Δxk+1=AcΔxk+B1cΔuk+B2cΔwk,yk=CΔxk+yk-1,where Δx(k)=x(k)-x(k-1);Δu(k)=u(k)-u(k-1);Δw(k)=w(k)-w(k-1).According to the MPC theory, the most recent measured value is taken as the initial condition to predict the future dynamic based on the predictive model. We make the following assumptions in the work:(
1
) The predicted time domain is p; the control time domain is m; and p≥m (m=3, p=10).(
2
) At the current k moment, the estimated and measured values (x^(k) and y(k)) of state vector are known.(
3
) The controlled variable remains unchanged outside of control time domain. Δu(k+i)=0(i≥m).(
4
) The measurable disturbance variable remains unchanged after moment k. Δw(k+i)=0(i≥1).Within predicted time domainp, the prediction expression for the pth step state vector at moment k is as follows:(21)Δxk+1∣k=AcΔx^k+B1cΔuk+B2cΔwkΔxk+2∣k=Ac2Δx^k+AcB1cΔuk+B1cΔuk+1+AcB2cΔwk⋮Δxk+p∣k=AcpΔxk+Acp-1B1cΔuk+Acp-2B1cΔuk+1+⋯+B1cΔuk+p-1+Acp-1B2cΔwk.At moment k, we predict the output from k+1 to k+p according to the predicted value of state vector in (22)yk∣k=ykyk+1∣k=CAcΔx^k+CB1cΔuk+CB2cΔwk+ykyk+2∣k=CAc2+CAcΔx^k+CAcB1c+CB1cΔuk+CB1cΔuk+1+CAcB2c+CB2cΔwk+ykyk+3∣k=∑i=13CAciΔx^k+∑i=13CAci-1B1cΔuk+⋯+CB1cΔuk+2+∑i=13CAci-1B2cΔuk+yk⋮yk+p∣k=∑i=1pCAciΔx^k+∑i=1pCAci-1B1cΔuk+⋯+CB1cΔuk+p-1+∑i=1pCAci-1B2cΔuk+yk.The vectors Yp(k) and ΔUm(k) are defined as follows:(23)Ypk=yk+1∣kyk+2∣k⋯yk+p∣kT,ΔUmk=ΔukΔuk+1⋯Δuk+m-1T.The pth step control output of the prediction system can be expressed as(24)Ypk=SxΔx^k+Lyyk+SuΔUk+SwΔwk,where(25)Sx=CAcCAc2+CAc⋮∑i=1pCAci;Sw=CB2cCAcB2c+CB2c⋮∑i=1pCAci-1B2c;Lc=I4I4⋮I4,Su=CB1c00⋯0CAcB1cCB1c0⋯0⋮⋮⋮⋯⋮∑i=1pCAciB1c∑i=1p-1CAci-1B1c⋯⋯∑i=1p-m+1CAci-1B1c.Similarly, we can deduce thepth step predictive equation with constraint output.MPC is an algorithm of determining control strategy by optimization. Firstly, it is necessary to determinem control increments Δu(k),Δu(k+1),…,Δu(k+m-1) from moment k. In the next p moment, the predicted output values y(k+1∣k),y(k+2∣k),…,y(k+p∣k) of the controlled object can approximate to the given expected value d(k+i),i=1,2,…,p.In the control process, the control increment should not have sharp change. Based on this factor, soft constraints are added to optimization performance indexes. Thus, the optimization performance index of momentk can be expressed as(26)minJk=∑i=1pqidk+i-yk+1+∑j=1mriΔu2k+j-1,where qi and ri are the weighting coefficients, which represent error tracking and suppression of control increment. If the constraint is not considered, the above problem is described as follows.TakingΔUk=ΔukΔuk+1⋯Δuk+m-1T as optimal variable, the performance index (26) is minimized under prediction model (20). Likewise, optimal performance index (26) can also be written as the form of vector.(27)minJk=Dpk-YpkQ2+ΔUmkR2,where Dpk=d(k+1)d(k+2)⋯d(k+p)T is the expected output value in the optimization time domain, namely, side slip angle and yaw rate based on vehicle model output. It is assumed that the expected value is fixed in the optimization time domain; namely, d(k+1)=d(k+2)=⋯=d(k+i);Q=diagq1q2⋯qp; and R=diagr1r2⋯rm are diagonal matrices comprised of weighting coefficients.At the momentk, Dp(k), and Yp(k) are known. In order to minimize ΔUm(k) by J(k), the optimal control sequence is obtained according to the necessary condition of the extremum: dJ(k)/dΔUm(k)=0.(28)ΔU∗k=SuTQTQSu+QTQ-1SuTQTQEpk+1∣k,where Ep(k+1∣k) is control deviation, which can be calculated on line by the following equation:(29)Epk+1∣k=Dpk-SxΔx^k-Lyyk-SwΔwk.In each sampling period, the first element of the derived optimal solution ΔU∗(k) is the control increment of MPC.(30)Δuk=100ΔU∗k,where u(k)=u(k-1)+Δu(k) is the additional yaw moment provided for the torque distribution layer.The Lyapunov functionV∗(k) at the moment k is defined as the optimal performance index of open loop.(31)V∗k=defJ∗k.Relationship between V∗(k) and V(k+1) is established:(32)V∗k=Vk+1+Jy∗1∣k,u∗0∣k,Δu∗0∣k-Jyp∣k+1,up-1∣k+1,Δup-1∣k+1.Since V(k+1)≥V∗(k+1), we can derive(33)V∗k+1-V∗k≤y∗p+1∣k22-y∗1∣k22+u∗p-1∣k22-u∗0∣k22-Δu∗0∣k22. As y(p+1∣k)22≤y(1∣k)22 and u(p∣k)22≤u(0∣k)22, we can get(34)V∗k+1-V∗k≤0.According to the Lyapunov second stability theorem, the proposed MPC closed-loop system is stable.


## 2.4. Wheel Torque Distribution

For the tire, lower pavement adhesion consumption rate leads to greater adhesion margin and larger distance between the tire and nonlinear saturation zone, indirectly improving the stability of the vehicle. According to tire friction ellipse theory, the work distributes the wheel torque for the purpose of minimum sum of pavement adhesion consumption rates of four tires. Equation (35) shows the objective function.(35)minJμ=∑i=14CiFxi2+Fyi2μ2Fzi2,i=1,2,3,4,where Fxi is the tire longitudinal force; Fyi the tire lateral force; Fzi the vertical load on the wheel; Ci the weight coefficient within the interval (0,1). The tire longitudinal force Fxi can be obtained by (36)Fxi=Twi-Ieqω˙irw,where Twi is the drive torque on the wheel; rw the wheel radius; ω˙i the wheel rotation angle acceleration; Ieq the wheel equivalent inertia of in-wheel motor and reducer inertia.The vertical loadFzi of wheel consists of static load and dynamic load caused by longitudinal and lateral accelerations without considering the centroid change caused by vehicle rollover and the influence of suspension system. The vertical load of each wheel is calculated by (37)Fz1Fz2Fz3Fz4=mg2Lbbaa+maxh2L-1-111+mayhBL-bb-aa,where a and b are the distance from the centroid to front and rear axle; L is the wheelbase; h is the height of centroid; B is the trend.The objective function in (35) determines smallest sum of the load rates of four tires. However, it is impossible to ensure that each tire works at a low load-rate state. Therefore, the utilization factor of each tire is adjusted by the weighting factor Ci. Under different driving conditions, Ci is not a fixed value. In simulation and debugging process, it is found that the value has great relationship with tire vertical load and road adhesion condition. The following equation is selected in the work:(38)Ci=μFzimg2,i=1,2,3,4.Here, we focus on longitudinal driving torque distribution in the conventional operating conditions. It is difficult to acquire the lateral force of tire in actual driving, so the pavement adhesion consumption rate is simplified to the longitudinal adhesion consumption rate. Equation (35) can be transformed into the following equation:(39)minJμ=∑i=14CiFxiμFzi2.When the yaw moment control is not required, the sum of longitudinal forces of four wheels satisfies the total moment demand, which is obtained by analyzing the pedal input. Therefore, constraint equation is expressed as follows:(40)∑i=14Fxi=fθp,where θp is the accelerator pedal opening; f(~) the pedal analysis strategy. In addition, the longitudinal force of tire is also restricted by pavement adhesion condition and motor output capacity.(41)-μFzi≤Fxi≤min⁡μFzi,Fmax,i=1,2,3,4,where Fmax is the maximum wheel drive force provided by the motor; namely, Fmax=Twmax/rw. Therefore, the optimization problem can be expressed as follows:(42)minJμ=∑i=14CiFxiμFzi2s.t.∑i=14Fxi=fθp-μFzi≤Fxi≤min⁡μFzi,FmaxFx2-Fx1B2+Fx4-Fx3B2=Tz.The nonlinear constrained optimization problem is solved by numerical optimization method to obtain the longitudinal forces for four tires. Then, we derive the drive torques which should be allocated to four wheels by calculation.


## 3. Control Algorithm Simulation Validation

The proposed controller was simulated on a 4WID vehicle in the Carsim and MATLAB environment. The vehicle under investigation is the light-duty off-road vehicle with the specifications which are listed in Table1.Table 1
Specifications of the experimental electric vehicle.
Parameter
Value
Mass
3565 kg
Yaw moment of inertia
7016 kg⋅m2
Distance from CG to front axle
1.593 m
Distance from CG to rear axle
1.707 m
Tire radius
0.435 m
Maximum motor torque
250 Nm
Maximum motor power
55.5 kW
Maximum motor speed
6000 rpm
Planetary reducer ratio
9.15The proposed controller (MPC) in this paper and its results are compared with the proposed controller (DOB) by [15]. Furthermore, the uncontrolled state is simulated, too.


### 3.1. Yaw Stability Control Algorithm Validation

When the steering wheel reaches certain rotor angle, the vehicle will be in danger of yaw instability on low adhesion road. The effect of control algorithm on vehicle stability is verified by selecting the double lane condition on low adhesion road. The specific simulation conditions are as follows: the road adhesion coefficient is 0.3; the driving speed is 60 km/h. Figure6(a) shows the vehicle trajectory under three different control modes. Figure 6(b) shows steering angle for the lane-change maneuver. Figures 6(c)–6(e) give the response curves of the yaw rate, sideslip angle, and lateral acceleration under the above three kinds of control modes.Figure 6
Double lane-change simulation test on low adhesion road.
(a)
Vehicle travel path
(b)
Steering wheel angle input
(c)
Yaw rate
(d)
Side slip angle
(e)
Lateral acceleration
(f)
Wheel torqueFrom Figure6, it is seen that without the stability control, when the steering wheel reaches a certain angle, the vehicle loses stability to cause the failure of lane-change. The yaw rate and the side slip angle of the vehicle change in a wide range, which cannot effectively track the expected value. After the MPC-based or DOB-based stability control, the yaw rate and side slip angle of vehicle can effectively track the expected value to pass through the test road according to predetermined trajectory, maintaining high speed stability. It can be well controlled under these two kinds of controllers. Moreover, the MPC-based controller can achieve the same smooth nonflutter control effect as the DOB-based controller. Figure 6(f) shows the wheel torque calculated by MPC-based stability control algorithm.


### 3.2. Verification of Rollover Stability Control Algorithm

The hook condition is conducted on the road with an adhesion coefficient of 0.85. The initial speed is 80 km/h, and the maximum steering angle is 180°. Figure7(a) shows the hook condition steering wheel angle input. Figures 7(b)–7(e) give the response curves of the roll, roll rate, yaw rate, sideslip angle, and longitudinal velocity under the three kinds of control modes. Figure 7(g) shows the vehicle trajectory under three different control modes.Figure 7
Hook condition simulation test on high adhesion road.
(a)
Steering wheel angle input
(b)
Roll
(c)
Roll rate
(d)
Yaw rate
(e)
Side slip angle
(f)
Longitudinal velocity
(g)
Vehicle travel path
(h)
Wheel torqueIn hook condition, the vehicle will lose stability and roll over if there is no stability control. After MPC-based or DOB-based stability control, the vehicle does not roll over as the roll is in a reasonable range. It can be seen that the response curves of MPC-based controller are as smooth as DOB-based controller. In addition, it can also be concluded that the stability will achieve a litter better accuracy under the proposed MPC-based stability controller and keep a higher longitudinal velocity. Figure7(h) indicates the wheel torque calculated by the MPC-based stability control algorithm.


## 3.1. Yaw Stability Control Algorithm Validation

When the steering wheel reaches certain rotor angle, the vehicle will be in danger of yaw instability on low adhesion road. The effect of control algorithm on vehicle stability is verified by selecting the double lane condition on low adhesion road. The specific simulation conditions are as follows: the road adhesion coefficient is 0.3; the driving speed is 60 km/h. Figure6(a) shows the vehicle trajectory under three different control modes. Figure 6(b) shows steering angle for the lane-change maneuver. Figures 6(c)–6(e) give the response curves of the yaw rate, sideslip angle, and lateral acceleration under the above three kinds of control modes.Figure 6
Double lane-change simulation test on low adhesion road.
(a)
Vehicle travel path
(b)
Steering wheel angle input
(c)
Yaw rate
(d)
Side slip angle
(e)
Lateral acceleration
(f)
Wheel torqueFrom Figure6, it is seen that without the stability control, when the steering wheel reaches a certain angle, the vehicle loses stability to cause the failure of lane-change. The yaw rate and the side slip angle of the vehicle change in a wide range, which cannot effectively track the expected value. After the MPC-based or DOB-based stability control, the yaw rate and side slip angle of vehicle can effectively track the expected value to pass through the test road according to predetermined trajectory, maintaining high speed stability. It can be well controlled under these two kinds of controllers. Moreover, the MPC-based controller can achieve the same smooth nonflutter control effect as the DOB-based controller. Figure 6(f) shows the wheel torque calculated by MPC-based stability control algorithm.


## 3.2. Verification of Rollover Stability Control Algorithm

The hook condition is conducted on the road with an adhesion coefficient of 0.85. The initial speed is 80 km/h, and the maximum steering angle is 180°. Figure7(a) shows the hook condition steering wheel angle input. Figures 7(b)–7(e) give the response curves of the roll, roll rate, yaw rate, sideslip angle, and longitudinal velocity under the three kinds of control modes. Figure 7(g) shows the vehicle trajectory under three different control modes.Figure 7
Hook condition simulation test on high adhesion road.
(a)
Steering wheel angle input
(b)
Roll
(c)
Roll rate
(d)
Yaw rate
(e)
Side slip angle
(f)
Longitudinal velocity
(g)
Vehicle travel path
(h)
Wheel torqueIn hook condition, the vehicle will lose stability and roll over if there is no stability control. After MPC-based or DOB-based stability control, the vehicle does not roll over as the roll is in a reasonable range. It can be seen that the response curves of MPC-based controller are as smooth as DOB-based controller. In addition, it can also be concluded that the stability will achieve a litter better accuracy under the proposed MPC-based stability controller and keep a higher longitudinal velocity. Figure7(h) indicates the wheel torque calculated by the MPC-based stability control algorithm.


## 4. Verification of Real Vehicle Test

According to ISO/FDIS 3888-1 standard, we performed double lane-change field test of light-duty off-road vehicle. The piles were arranged on the test field (see Figure8).Figure 8
Double lane-change test.The gyroscope was installed in the vehicle for real-time acquisition of the yaw rate and lateral acceleration of the vehicle. The GPS speed sensor was used to measure the longitudinal and lateral velocities. The angular velocities of four wheels were collected by in-wheel motor speed sensor. We measured the steering wheel angle by steering wheel angle sensor. The test speed of vehicle was 80 km/h (±2 km/h).Figures9–11 show the double lane-change test results with MPC-based stability control, with DOB-based stability control, and without control, respectively. It can be seen that the vehicle without control has poor stability because of large difference between yaw rate and reference value. After MPC-based or DOB-based stability control, the difference obviously decreases to strengthen vehicle stability. Figure 12 shows double lane-change driving route. If the steering angle reaches a certain value, the vehicle without stability control will lose stability to approach the curb, thus failing to track the desired trajectory. With MPC-based or DOB-based stability control, the vehicle can pass through the test road in accordance with the desired trajectory. The control effects of the two controllers are comparable, and the MPC-based control in the direction of the moment of emergency change is relatively more stable.Figure 9
Double lane-change test results with MPC-based stability control.
(a)
Steering wheel angle
(b)
Yaw rateFigure 10
Double lane-change test results with DOB-based stability control.
(a)
Steering wheel angle
(b)
Yaw rateFigure 11
Double lane-change test results without stability control.
(a)
Steering wheel angle
(b)
Yaw rateFigure 12
Double lane-change driving route.


## 5. Conclusions

In the work, we proposed the lateral stability control method of in-wheel motor drive off-road vehicle based on MPC. The multiinput and multioutput system predictive control law was designed by establishing vehicle dynamics model for predictive control. The motor braking/driving torque was directly generated to conduct yaw stability control under low adhesion coefficient and double lane-change conditions as well as rollover stability control algorithm validation under high adhesion coefficient and hook conditions. At last, real vehicle test verification was performed under high speed double lane-change condition. This paper only carries out the test under the condition of invariable road surface and constant velocity but does not verify the proposed MPC-based stability control under the condition of changing road surface and time-varying velocity. Subsequent further research, we will verify it.Consequently, the following conclusions were obtained.(
1
) Distributed electric vehicle braking/driving force coordination control made full use of multiwheel adhesion to provide vehicle with greater yaw moment and faster response, increasing vehicle stability under limit conditions.(
2
) With high accuracy and good robustness, MPC was suitable for solving complex coordination problem for multiple actuators. It could also be applied in lateral stability control of distributed electric drive vehicle.(
3
) Simulation result and real vehicle test showed that MPC algorithm had better effect in integrated control of yaw/rollover stability, thus greatly improving the lateral stability of vehicle.



---

*Source: 6080763-2018-02-07.xml*