Development and evaluation of torque distribution methods for four-wheel independent-drive electric vehicle based on parameter estimation

Abstract

A hierarchical controller including a vehicle condition parameter estimator for improving the vehicle stability is theoretically applied to a four-wheel independent-drive electric vehicle. This study is aimed at combining the parameters estimation with vehicle stability control and evaluating different optimization algorithms based on the computation time and fluctuation. The hierarchical controller consists of a vehicle condition parameter estimator, active yaw moment controller, and torque distribution controller. The vehicle condition parameter estimator is based on the Unscented Kalman Filter, which avoids the truncation error caused by Taylor expansion to convert a nonlinear system into a linear system. The active yaw moment controller is designed based on the sliding mode control. The controller employs co-control between the vehicle sideslip angle and yaw rate, which properly accounts for both factors. The objective function of the torque distribution controller is based on the tire utilization function and constraint condition. The torque distribution and wheel slip rate are accounted for simultaneously. The torque distribution methods are tested based on nonlinear programming, quadratic programming, and weighted least squares for the sake of comparison. The simulation results demonstrate the accuracy of the proposed vehicle state estimator and the effectiveness of the proposed stability controller.

Introduction

A four-wheel independent-drive (4WID) electric vehicle represents an innovative design in which four in-wheel motors are used to drive the wheels independently. The electric vehicle’s transmission efficiency is also improved because the conventional mechanical drive system is eliminated. The independent four-motor design also frees up interior space, is more conducive to wire-controlled applications, and allows a variety of novel control methods to be applied. 4WID electric vehicles can achieve electric power steering, differential drive assisted steering control, ¹ electronic differentials, ² energy regenerative braking, ³ and direct yaw moment control, ⁴ making the design extremely popular. ⁵

However, it is difficult to obtain the vehicle condition parameters directly through a sensor in a vehicle with a 4WID form; there is, therefore, currently a demand for novel vehicle state estimation algorithms suited to this purpose. Imslanda et al. ⁶ explored a new method for obtaining the longitudinal speed from the wheel speed and longitudinal acceleration in which the current state of the vehicle is initially determined, followed by a comprehensive calculation of the wheel speed and longitudinal acceleration by adjusting the weight coefficient to determine the vehicle speed. Hrgetic et al. ⁷ established another vehicle state (speed) estimation method using a nonlinear, adaptive extended Kalman filter (EKF). Li and Zhang ⁸ established another vehicle sideslip angle estimation method using a hybrid Kalman filter. However, the above researches did not apply the estimated results to the vehicle stability control. Therefore, the accuracy of the estimated state variables in vehicle control cannot be verified. On the contrary, general vehicle parameter estimation methods, such as EKF, are easy to introduce truncation errors in the estimation process.

There has been a significant amount of research to date on independent driving, electric vehicle stability control, and torque distribution. The available methods for solving the torque distribution issue include a generalized inverse method, ⁹ penalty function method, ¹⁰ chain increment method, and mathematical optimization method. Maeda et al. ¹¹ investigated various independent-drive electric vehicle torque distribution strategies and built test prototypes to verify them; Li et al. ¹² designed a yaw moment controller based on sliding mode control (SMC) and a torque distribution controller based on a penalty function. The above stability control methods of the vehicle are based on the assumption that the vehicle state has been obtained, and there are few studies comparing the computational time costs and fluctuation of different control methods.

In order to estimate vehicle parameters and verify the application of those parameters in vehicle stability control, the vehicle stability controller with parameter estimator is built. The Unscented Kalman Filter (UKF)–based estimation method avoids the effect of truncation error on the accuracy of results. The main aspects of this article are as follows: To solve the difficulty in obtaining the vehicle condition parameters of 4WID vehicles, this study establishes a hierarchical control system including a vehicle condition parameters estimator. In addition, the computational times required by several optimization algorithms in a MATLAB environment are compared and analyzed, which can provide a reference for the selection of varying optimization algorithms.

The remainder of this article is as follows: The upper level vehicle condition parameter estimator is described in “Vehicle condition parameters estimator” section. “Active yaw moment controller” section introduces the middle level active yaw moment controller. The design of the lower level torque distribution controller is presented in “Torque distribution controller” section. In “Simulation results” section, the co-simulation conditions based on CarSim and MATLAB/Simulink are described, and the results are analyzed. Finally, some concluding remarks are given in “Conclusion” section.

Vehicle condition parameters estimator

A block diagram of the proposed stability control scheme is shown in Figure 1. This hierarchical controller was designed for control of the handling stability of electric vehicles. We established a vehicle condition parameter estimator, an active yaw moment controller, and a torque distribution controller to support better handling stability. The vehicle state parameter estimator uses an UKF method to estimate the longitudinal velocity, lateral velocity, yaw rate, sideslip angle, and tire forces. The proportion-integral (PI) control was used to calculate the generalized longitudinal force required for vehicle driving. The active yaw moment control layer adopts a SMC method to obtain the additional yaw moment necessary to maintain the handling stability of the vehicle by the joint control of the sideslip angle and yaw rate. This additional yaw moment is used alongside the generalized longitudinal force demand as the constraint condition under which the tire utilization rate serves as the objective function. The torque distribution was achieved through three mathematical optimization methods, and then the independent control of the four drive wheels was achieved such that the vehicle effectively retains its handling stability even under high-speed or low road friction coefficient conditions.

Figure 1.

Block diagram of the proposed stability control scheme.

Seven degrees of freedom vehicle model

Assume that the vehicle model only has seven degrees of freedom (7DOF), namely, longitudinal displacement, lateral displacement, rotation around the z-axis, and rotational freedom of the four wheels. According to the model shown in Figure 2 and the above assumption, Newton’s theory of mechanics can be used to establish the following vehicle motion equation

$$\{\begin{array}{c}M{a}_x=(F_{\mathrm{xfl}}+F_{\mathrm{xfr}})\cos \delta -(F_{\mathrm{yfl}}+F_{\mathrm{yfr}})\sin \delta +F_{\mathrm{xrl}}+F_{\mathrm{xrr}}\hfill \\ \multicolumn{1}{c}{M{a}_y=\sum _{\frac{i=f,l}{j=f,r}}(F_{\mathrm{xij}}\sin \delta +F_{\mathrm{yin}}\cos \delta )}\\ \multicolumn{1}{c}{I_z\stackrel{·}{\gamma }=\sum _{\frac{i=f,l}{j=f,r}}(F_{\mathrm{xij}}\sin \delta +F_{\mathrm{yin}}\cos \delta )a_i+\sum _{\frac{i=f,l}{j=f,r}}(F_{\mathrm{xij}}\cos \delta +F_{\mathrm{yin}}\sin \delta )\frac{l_i}{2}}\end{array}$$

(1)

where $a_f=a$ and $a_r=b$ denote the distance from the center of gravity (CG) to the front and rear axle, respectively.

Figure 2.

7DOF nonlinear vehicle model.

The four-wheel rotation equations are

$$\begin{array}{c}{J}_{\omega }{\stackrel{·}{\omega }}_{\mathrm{fl}}=T_{\mathrm{dfl}}-F_{\mathrm{xfl}}R-M_{\mathrm{ffl}}\hfill \\ J_{\omega }{\stackrel{·}{\omega }}_{\mathrm{fr}}=T_{\mathrm{dfr}}-F_{\mathrm{xfr}}R-M_{\mathrm{ffr}}\hfill \\ J_{\omega }{\stackrel{·}{\omega }}_{\mathrm{rl}}=T_{\mathrm{drl}}-F_{\mathrm{xrl}}R-M_{\mathrm{frl}}\hfill \\ J_{\omega }{\stackrel{·}{\omega }}_{\mathrm{rr}}=T_{\mathrm{drr}}-F_{\mathrm{xrr}}R-M_{\mathrm{frr}}\hfill \end{array}$$

(2)

where the subscripts fl, fr, rl, and rr denote the front-left wheel, front-right wheel, rear-left wheel, and rear-right wheel, respectively, and have the same meaning in the following sections.

Vehicle condition parameters estimation

It is more straightforward to estimate the probability density distribution of nonlinear functions than the nonlinear function itself, ¹³ and thus we designed our vehicle condition parameter estimation using a UKF to approximate a nonlinear function. The core algorithm of UKF is unscented transform (UT), that is, in the original state distribution, the set of points is selected according to certain rules such that the mean and covariance of the set are equal to the mean and covariance of the original state distribution. The point set is then plugged into the nonlinear function to obtain the corresponding value of the nonlinear function. Unlike an EKF estimation method, UKF does not require calculating the Jacobian matrix complex, thus preventing truncation errors owing to the Taylor expansion of the nonlinear system.

The general nonlinear system can be described using a state space equation encompassing the state equation and measurement equation. Its discrete form is as follows

$$\begin{array}{c}{\mathit{x}}_{k+1}=f({\mathit{x}}_k,{\mathit{u}}_k)+{\mathit{w}}_k\hfill \\ {\mathit{y}}_{k+1}=h({\mathit{x}}_{k+1},{\mathit{u}}_{k+1})+{\mathit{v}}_{k+1}\hfill \end{array}$$

(3)

where x is the system state variable, y is the system measurement variable, u is the import variable, w is the system process noise with variance Q , and v is the system measurement noise with variance $\widehat{\mathit{R}}$ .

The necessary vehicle state parameters defined for observation are as follows

$$\begin{array}{c}\multicolumn{1}{c}{\mathit{x}={\left[\begin{array}{ccccccccccc}{x}_1\hfill & x_2\hfill & x_3\hfill & x_4\hfill & x_5\hfill & x_6\hfill & x_7\hfill & x_8\hfill & x_9\hfill & x_{10}\hfill & x_{11}\hfill \end{array}\right]}^T}\\ \multicolumn{1}{c}{={\left[\begin{array}{ccccccccccc}{v}_x\hfill & v_y\hfill & \gamma \hfill & F_{\mathrm{yfl}}\hfill & F_{\mathrm{yfr}}\hfill & F_{\mathrm{yrl}}\hfill & F_{\mathrm{yrr}}\hfill & F_{\mathrm{zfl}}\hfill & F_{\mathrm{zfr}}\hfill & F_{\mathrm{zrl}}\hfill & F_{\mathrm{zrr}}\hfill \end{array}\right]}^T}\end{array}$$

(4)

The measurement variables are

$$\mathit{y}={\left[\begin{array}{ccc}{y}_1\hfill & y_2\hfill & y_3\hfill \end{array}\right]}^T={\left[\begin{array}{ccc}{a}_x\hfill & a_y\hfill & \gamma \hfill \end{array}\right]}^T$$

(5)

The longitudinal tire force can be obtained using the motor torque and equation (2), and $\delta$ , $a_x$ , and $a_y$ can be measured using sensors, the input variable of which is as follows

$$\begin{array}{l}\mathit{u}={\left[\begin{array}{ccccccc}{u}_1& u_2& u_3& u_4& u_5& u_6& u_7\end{array}\right]}^T\\ ={\left[\begin{array}{ccccccc}\delta & a_x& a_y& F_{xfl}& F_{xfr}& F_{xrl}& F_{xrr}\end{array}\right]}^T\end{array}$$

(6)

According to the vehicle model, the magic formula tire model ¹⁴ and tire vertical force formula comprise the estimator state equation (7), and the measurement equation (8), shown below

$$\{\begin{array}{cc}\multicolumn{1}{c}{{\stackrel{·}{v}}_x=}& \frac{v_y\gamma +[(F_{\mathrm{xfl}}+F_{\mathrm{xfr}})\cos \delta -F_f-(F_{\mathrm{yfl}}+F_{\mathrm{yfr}})\sin \delta ]}{M}+w_1\hfill \\ \multicolumn{1}{c}{{\stackrel{·}{v}}_y=}& \frac{-v_x\gamma +[(F_{\mathrm{xfl}}+F_{\mathrm{xfr}})\sin \delta +(F_{\mathrm{yfl}}+F_{\mathrm{yfr}})\cos \delta +F_{\mathrm{yrl}}+F_{\mathrm{yrr}}]}{M}+w_2\hfill \\ \multicolumn{1}{c}{\stackrel{·}{\gamma }=}& \frac{[a(F_{\mathrm{xfl}}\sin \delta +F_{\mathrm{yfl}}\cos \delta +F_{\mathrm{xfr}}\sin \delta -F_{\mathrm{yfr}}\cos \delta )-\frac{l_f}{2}(F_{\mathrm{xfr}}\cos \delta -F_{\mathrm{yfr}}\sin \delta -F_{\mathrm{xfl}}\cos \delta +F_{\mathrm{yfl}}\sin \delta )-b(F_{\mathrm{xrl}}+F_{\mathrm{xrr}})+\frac{l_r}{2}(F_{\mathrm{xrr}}-F_{\mathrm{xrl}})]}{I_z}+w_3\hfill \\ \multicolumn{1}{c}{F_{\mathrm{yi}}=}& D\sin \left(C\arctan (B{X}_1-E(B{X}_1-\arctan \left(B{X}_1\right)))\right)+w_j(i=\mathrm{fl},\mathrm{fr},\mathrm{rl},\mathrm{rr};j=4,5,6,7)\hfill \\ \multicolumn{1}{c}{F_{\mathrm{zi}}=}& \frac{1}{2}\mathrm{Mg}·\frac{b}{l}-\frac{1}{2}M{a}_x·\frac{h_g}{l}\pm \frac{1}{2}M{a}_y·\frac{b}{l}·\frac{h_g}{l_f}+w_j(i=\mathrm{fl},\mathrm{fr};j=8,9)\hfill \\ \multicolumn{1}{c}{F_{\mathrm{zi}}=}& \frac{1}{2}\mathrm{Mg}·\frac{a}{l}+\frac{1}{2}M{a}_x·\frac{h_g}{l}\pm \frac{1}{2}M{a}_y·\frac{a}{l}·\frac{h_g}{l_r}+w_j(i=\mathrm{rl},\mathrm{rr};j=10,11)\hfill \end{array}$$

(7)

$$\{\begin{array}{c}{a}_x=\frac{[(F_{\mathrm{xfl}}+F_{\mathrm{xfr}})\cos \delta -F_f-(F_{\mathrm{yfl}}+F_{\mathrm{yfr}})\sin \delta ]}{M}+v_1\hfill \\ a_y=\frac{[(F_{\mathrm{xfl}}+F_{\mathrm{xfr}})\sin \delta +(F_{\mathrm{yfl}}+F_{\mathrm{yfr}})\cos \delta +F_{\mathrm{yrl}}+F_{\mathrm{yrr}}]}{M}+v_2\hfill \\ \gamma =\gamma +v_3\hfill \end{array}$$

(8)

The vehicle sideslip angle is

$$\beta =\arctan \frac{v_y}{v_x}$$

(9)

where B, C, D, and E are the fitting coefficient for the magic tire formula. ¹⁵

UKF estimation process

Filter initialization

The nonlinear system has system and measurement noises, and thus the extended vector is defined as ${\mathit{x}}_k^a=[x_k{\omega }_k{v}_k{]}^T$ . The initial preconditions are as follows

$$\{\begin{array}{cc}\multicolumn{1}{c}{{\overline{\mathit{x}}}_0=}& E\left({\mathit{x}}_0\right)\hfill \\ \multicolumn{1}{c}{{\mathit{P}}_0=}& E\left[({\mathit{x}}_0-{\overline{\mathit{x}}}_0){({\mathit{x}}_0-{\overline{\mathit{x}}}_0)}^T\right]\hfill \\ \multicolumn{1}{c}{{\overline{\mathit{x}}}_0^a=}& E\left[{\mathit{x}}^a\right]=\left[\begin{array}{ccc}{\overline{\mathit{x}}}_0^T\hfill & 0\hfill & 0\hfill \end{array}\right]\hfill \\ \multicolumn{1}{c}{{\mathit{P}}_0^a=}& E\left[({\mathit{x}}_0^a-{\overline{\mathit{x}}}_0^a){({\mathit{x}}_0^a-{\overline{\mathit{x}}}_0^a)}^T\right]\hfill \end{array}$$

(10)

The 2n + 1 Sigma sampling point

$$\{\begin{array}{c}{\chi }_{0,k|k}^a={\overline{x}}_{k|k}^a\\ {\chi }_{i,k|k}^a={\overline{x}}_{k|k}^a+\sqrt{\left(n_a+\lambda \right)P_{k|k}^a}i=1,\dots ,n\\ {\chi }_{i,k|k}^a={\overline{x}}_{k|k}^a-\sqrt{\left(n_a+\lambda \right)P_{k|k}^a}i=n+1,\dots ,2n\end{array}$$

(11)

The Sigma point weight

$$\{\begin{array}{c}{W}_0^{(m)}=\frac{\lambda }{(n_a+\lambda )}\hfill \\ W_0^{(c)}=\frac{\lambda }{(n_a+\lambda )}+(1-{\alpha }^2+{\beta }_t)\hfill \\ W_i^{(m)}=W_i^{(c)}=\frac{1}{\left[2(n_a+\lambda )\right]}i=1,\dots ,2n\hfill \end{array}$$

(12)

where $n_a=2n+n_v$ , $\lambda ={\alpha }^2(n_a+e)-n_a$ , n is the dimension of the state vector, and $n_v$ is the measured noise dimension. Here, e is usually set to zero to ensure a positive definite of the covariance matrix, and $\alpha$ determines the intensity of the sigma point near ${\overline{\mathit{x}}}_0^a$ , usually ranging from ${10}^{-4}\le \alpha \le 1$ . In this article, $\alpha ={10}^{-3}$ and ${\beta }_t=2$ .

UKF calculation

The time update functions

$$\{\begin{array}{cc}\multicolumn{1}{c}{{\chi }_{i,k+1|k}^a=}& f({\chi }_{i,k|k}^a,u_k)\hfill \\ \multicolumn{1}{c}{{\overline{x}}_{k+1|k}^a=}& \sum _{i=0}^{2n_a}{W}_i^{\left(m\right)}{\chi }_{i,k+1|k}^a\hfill \\ \multicolumn{1}{c}{y_{i,k+1|k}^a=}& h({\chi }_{i,k+1|k}^a,u_k)\hfill \\ \multicolumn{1}{c}{{\overline{y}}_{k+1|k}^a=}& \sum _{i=0}^{2n_a}{W}_i^{\left(m\right)}{y}_{i,k+1|k}^a\hfill \\ \multicolumn{1}{c}{P_{k+1|k}^a=}& \sum _{i=0}^{2n_a}{W}_i^{\left(c\right)}({\chi }_{i,k+1|k}^a-{\overline{\chi }}_{k+1|k}^a){({\chi }_{i,k+1|k}^a-{\overline{\chi }}_{k+1|k}^a)}^T\hfill \end{array}$$

(13)

The measurement update functions

$$\{\begin{array}{c}{P}_{(y,y),k+1|k}^a=\sum _{i=0}^{2n_a}{W}_i^{\left(c\right)}(y_{k+1|k}^a-{\overline{y}}_{k+1|k}^a){(y_{k+1|k}^a-{\overline{y}}_{k+1|k}^a)}^T\hfill \\ \multicolumn{1}{c}{P_{(x,y),k+1|k}^a=\sum _{i=0}^{2n_a}{W}_i^{\left(c\right)}(x_{k+1|k}^a-{\overline{x}}_{k+1|k}^a){(x_{k+1|k}^a-{\overline{x}}_{k+1|k}^a)}^T}\\ \multicolumn{1}{c}{K_{k+1}=P_{(x,y),k+1|k}^a{\left(P_{(y,y),k+1|k}^a\right)}^{-1}}\\ \multicolumn{1}{c}{P_{k+1|k+1}^a=P_{k+1|k}^a-K_{k+1}{P}_{(y,y),k+1|k}^a{K}_{k+1}^T}\\ \multicolumn{1}{c}{{\overline{x}}_{k+1|k+1}^a={\overline{x}}_{k+1|k}^a+K_{k+1}(y_k^a-{\overline{y}}_{k+1|k}^a)}\end{array}$$

(14)

Active yaw moment controller

Two degrees of freedom vehicle model

The vehicle model can be considered as having two degrees of freedom (2DOF), namely, lateral motion along the y-axis and yaw motion about the z-axis, when the vehicle is driving at a constant velocity. The following assumptions are made for a vehicle model: (1) The lateral acceleration is below 0.4 g, (2) the lateral characteristics of the tire are within a linear range, (3) the suspension effects are ignored, (4) the aerodynamics are ignored, and (5) under a steady state, ${\stackrel{·}{v}}_x=0$ and $\stackrel{·}{\gamma }=0$ . Based on the above assumptions, the established 2DOF vehicle model can be considered linear. The linear 2DOF vehicle model can reflect the linear relationship between the driver’s intention and the vehicle’s stability. The proposed SMC controller was established according to the 2DOF vehicle model, and the yaw moment was determined based on the corresponding longitudinal force. The 2DOF vehicle motion differential equations are

$$M{v}_x(\stackrel{·}{\beta }+\gamma )=Y_f+Y_r$$

(15)

$$I_z\stackrel{·}{\gamma }=a{Y}_f-b{Y}_r$$

(16)

The following equation can be obtained from the vehicle linear 2DOF model shown in Figure 3

$$\begin{array}{c}(k_1+k_2)\beta +\frac{1}{v_x}(a{k}_1-b{k}_2)\gamma -k_1\delta =M{v}_x(\stackrel{·}{\beta }+\gamma )\hfill \\ \multicolumn{1}{c}{(a{k}_1-b{k}_2)\beta +\frac{1}{v_x}(a^2{k}_1+b^2{k}_2)\gamma -a{k}_1\delta =I_z\stackrel{·}{\gamma }}\end{array}\}$$

(17)

Figure 3.

2DOF vehicle model.

Considering the state of stability of the vehicle, the steady-state yaw rate gain and steady-state sideslip angle gain are obtained from equation (17).

The steady-state yaw rate gain is as follows

$${\left(\frac{{\gamma }_d}{\delta }\right)}_S=\frac{\frac{v_x}{l}}{1+{\mathrm{Kv}}_x^2}$$

(18)

The steady-state sideslip angle gain is the following

$${\left(\frac{{\beta }_d}{\delta }\right)}_S=\frac{\frac{b}{l}-\frac{{\mathrm{amv}}_x^2}{k_2{l}^2}}{1+{\mathrm{Kv}}_x^2}$$

(19)

where K is a stability factor

$$K=\frac{m}{l^2}(\frac{a}{k_2}-\frac{b}{k_1})$$

Because the 2DOF vehicle model analysis does not include the longitudinal motion of the vehicle, the effects of the longitudinal force on the yaw moment are not reflected in equations (15) and (16); in other words, they only include information about the yaw moment caused by the lateral tire force. However, for certain practical application, the control of the yaw moment is only achieved by controlling the torque of the in-wheel motor, that is, only the longitudinal tire force is controllable. Thus, to remedy this, it is necessary to added the following to equation (16)

$$I_z\stackrel{·}{\gamma }=a{Y}_f-b{Y}_r+M_{\mathrm{zx}}$$

(20)

where $M_{\mathrm{zx}}$ denotes the yaw moment caused by the longitudinal tire force, and in the next section, it also denotes the output of the active yaw moment controller.

Active yaw moment controller based on SMC

The sideslip angle and yaw rate are two important state variables that measure the vehicle stability. The goal of the stability control is to keep $\beta$ as small as possible while maintaining the expectations of $\gamma$ tracking. According to the coupling relationship between $\beta$ and $\gamma$ , as shown in equations (15) and (16), we propose a method of joint control of $\beta$ and $\gamma$ ; the proposed method resolves the issues wherein $\gamma$ cannot track the ideal value during independent control of $\beta$ , or when $\beta$ is too large for an independent control of $\gamma$ . ⁴ The sliding surface is defined as follows

$$S=\gamma -{\gamma }_d+c(\beta -{\beta }_d),c>0$$

(21)

where c is the joint control parameter.

According to the exponential reaching law obtained in this the arrival conditions ¹⁶

$$\stackrel{·}{S}+\mathrm{kS}=-\epsilon \mathrm{sgn}\left(S\right)k>0,\epsilon >0$$

(22)

According to the Lyapunov function $V=S^2/2$ , it is easy to prove

$$\stackrel{·}{V}=\stackrel{·}{S}S=-(\epsilon \left|S\right|+k{S}^2)\le 0$$

(23)

and therefore, the system is guaranteed to be stable.

The following equation holds

$$\stackrel{·}{\gamma }={\stackrel{·}{\gamma }}_d-c(\stackrel{·}{\beta }-{\stackrel{·}{\beta }}_d)-k[\gamma -{\gamma }_d+c(\beta -{\beta }_d)]-\epsilon \mathrm{sgn}\left(S\right)$$

(24)

where k is the coefficient of the defined sliding surface, and $k=0.2$ in this study.

Substituting equation (24) into equation (20) yields the additional yaw moment in which $\beta$ and $\gamma$ are combined

$$M_{z,(\beta +\gamma )}=-\left(Y_{f}a-Y_{r}b\right)+I_z\left\{{\stackrel{.}{\gamma }}_d-c\left(\stackrel{.}{\beta }-{\stackrel{.}{\beta }}_d\right)-k\left[\gamma -{\gamma }_d+c\left(\beta -{\beta }_d\right)\right]-\epsilon \mathrm{sgn}\left(S\right)\right\}$$

(25)

Torque distribution controller

Objective function

We used the minimum vehicle tire utilization as the objective function. ¹⁷ This allowed us to reduce the possibility of sliding or locking the wheel while reasonably allocating the load of each tire, ¹⁸ providing a sufficient tire adhesion margin to maintain the level of stability. If the all tire utilization reaches a maximum of 1, the tire meets the limits of adhesion and external influences threaten the vehicle stability.

The tire utilization formula is

$${\eta }_{\mathrm{ij}}=\frac{\sqrt{F_{\mathrm{xij}}^2+F_{\mathrm{yij}}^2}}{{\mu }_{\mathrm{ij}}{F}_{\mathrm{zij}}},\mathrm{ij}=\mathrm{fl},\mathrm{fr},\mathrm{rl},\mathrm{rr}$$

(26)

The minimum sum of the four tire utilization squares serves as the objective function. The 4WID electric vehicle has a longitudinal drive force that is independently controllable, whereas the lateral tire force cannot be directly controlled. Therefore, we revised the objective function accordingly to control the longitudinal tire force and minimize its utilization to avoid longitudinal tire force saturation to affect vehicle stability. The objective function obtained after rewriting is as follows

$$\underset{u}{min}J=\sum _{\mathrm{ij}=\mathrm{fl},\mathrm{fr},\mathrm{rl},\mathrm{rr}}\frac{F_{\mathrm{xij}}^2}{{\left({\mu }_{\mathrm{ij}}{F}_{\mathrm{zij}}\right)}^2}$$

(27)

The objective function is expressed in matrix form

$$\underset{u}{min}J={\mathit{U}}^T\mathit{W}\mathit{U}$$

(28)

where

$$\mathit{U}={\left[F_{\mathrm{xfl}}{F}_{\mathrm{xfr}}{F}_{\mathrm{xrl}}{F}_{\mathrm{xrr}}\right]}^T$$

and

$$\mathit{W}=\left[\begin{array}{cccc}\frac{1}{{\left({\mu }_{\mathrm{fl}}{F}_{\mathrm{zfl}}\right)}^2}\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& \frac{1}{{\left({\mu }_{\mathrm{fr}}{F}_{\mathrm{zfr}}\right)}^2}\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & \frac{1}{{\left({\mu }_{\mathrm{rl}}{F}_{\mathrm{zrl}}\right)}^2}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & \frac{1}{{\left({\mu }_{\mathrm{rr}}{F}_{\mathrm{zrr}}\right)}^2}\hfill \end{array}\right]$$

Objective function constraint

The active yaw moment controller is used to calculate the additional yaw moment required to maintain the vehicle stability during transit. The additional yaw moment is achieved by controlling the longitudinal force of the tire, and thus, it is necessary to distribute the wheel torque; the distributed torque must satisfy the driving force of the vehicle. We designed two control inputs for this purpose: the additional yaw moment $M_{z,(\beta +\gamma )}$ and generalized longitudinal force $F_{\mathrm{xG}}$ .

The additional yaw moment was calculated via equation (25) and the generalized longitudinal force via PI control. It is necessary to ensure that the actual longitudinal acceleration is close to the driver’s intended acceleration, which is expressed through the following equations

$$F_{\mathrm{xG}}=k_P(a_d-a_x)+k_I\underset{0}{\overset{t}{\int }}(a_d-a_x)\mathrm{dt}$$

(29)

where $F_{\mathrm{xG}}=F_{\mathrm{xfl}}\cos \delta +F_{\mathrm{xfr}}\cos \delta +F_{\mathrm{xrl}}+F_{\mathrm{xrr}}$ , $k_P$ is the proportional coefficient, for example, $k_P=20$ , and $k_I$ is the integral coefficient, for example, $k_I=10$ .

When the vehicle is operating, the power and stability are maintained as long as equation (30) is satisfied

$$\begin{array}{c}{F}_{\mathrm{xG}}=F_{\mathrm{xfl}}\cos \delta +F_{\mathrm{xfr}}\cos \delta +F_{\mathrm{xrl}}+F_{\mathrm{xrr}}\hfill \\ \multicolumn{1}{c}{M_{z,(\beta +\gamma )}=(F_{\mathrm{xfr}}-F_{\mathrm{xfl}})(\frac{l_f}{2}\cos \delta +a\sin \delta )+(F_{\mathrm{xrr}}-F_{\mathrm{xrl}})\frac{l_r}{2}}\end{array}\}$$

(30)

Equation (30) also reveals the relationship between the control and output parameters. Therefore, the equality constraint of the objective function can be obtained

$$\mathit{V}=\mathit{BU}$$

(31)

where $\mathit{V}=[\begin{array}{cc}{F}_{\mathrm{xG}}\hfill & M_{z,(\beta +\gamma )}\hfill \end{array}{]}^T$ for the control input matrix; the nonlinear weight coefficient matrix

$$\mathit{B}=\left[\begin{array}{cccc}\cos \delta \hfill & \cos \delta \hfill & 1\hfill & 1\hfill \\ \multicolumn{1}{c}{-\frac{l_f}{2}\cos \delta +a\sin \delta }& \frac{l_f}{2}\cos \delta +a\sin \delta \hfill & -\frac{l_r}{2}\hfill & \frac{l_r}{2}\hfill \end{array}\right]$$

The torque provided by the drive motor is limited by the maximum torque of the motor. According to the design requirements, $T_{\mathrm{ij}\max }$ is chosen as 280 $\mathrm{N}·\mathrm{m}$ . Thus, the inequality constraint of the objective function can be derived

$$-\frac{T_{\mathrm{ij}\max }}{R}\le F_{\mathrm{xij}}\le \frac{T_{\mathrm{ij}\max }}{R}$$

(32)

Objective function solution

Nonlinear programming method

The electric vehicle torque distribution problem was solved based on the nonlinear programming (NP). The standard form of the NP is

$$\begin{array}{cc}\multicolumn{1}{c}{min}& J=\sum _{\mathrm{ij}=\mathrm{fl},\mathrm{fr},\mathrm{rl},\mathrm{rr}}\frac{F_{\mathrm{xij}}^2}{{\left({\mu }_{\mathrm{ij}}{F}_{\mathrm{zij}}\right)}^2}\hfill \\ \multicolumn{1}{c}{s.t.}& \mathit{BU}=\mathit{V}\hfill \\ \multicolumn{1}{c}{}& {\mathit{U}}_{min}\le \mathit{U}\le {\mathit{U}}_{max}\hfill \end{array}$$

(33)

This problem can be solved by using fmincon tool in a MATLAB optimization toolbox. With this method, the solution method can be autonomously chosen according to the calculation applied, such as a sequence quadratic programming (SQP) algorithm, an interior point penalty function method, or a feasible set method, among others.

Quadratic programming method

According to the above optimization objective function and the constraint condition, the torque distribution problem can be solved via quadratic programming (QP). The standard form of the QP method is

$$\begin{array}{cc}\multicolumn{1}{c}{minJ}& ={\mathit{U}}^T\mathit{WU}\hfill \\ \multicolumn{1}{c}{s.t.}& \mathit{BU}=\mathit{V}\hfill \\ \multicolumn{1}{c}{}& {\mathit{U}}_{min}\le \mathit{U}\le {\mathit{U}}_{max}\hfill \end{array}$$

(34)

This problem can be solved using the quadprog tool in the MATLAB optimization toolbox.

By comparing equations (33) and (34), the objective function of the NP is similar to the QP. As the difference between the two, the objective function of the NP is not required to form the objective function, whereas QP requires transforming the objective function, for example, $({\mathit{x}}^T\mathit{Hx})/2+{\mathit{f}}^T\mathit{x}$ , as the standard form. In addition, the first algorithm is a combined algorithm that can select the appropriate algorithm based on the computational scale.

Weighted least squares method

The independent driving electric vehicle torque distribution problem can also be formulated as a weighted least squares (WLS) problem. Equation (34) is the optimal solution to the objective function within the range of equation (30), and thus, the objective function and constraint condition can be transformed into the SQP form

$$\{\begin{array}{c}\mathit{u}=\arg \underset{\mathit{U}\in \mathrm{\Omega }}{min}\Vert {\mathit{W}}_{\mathit{U}}(\mathit{U}-{\mathit{U}}_d){\Vert }_2\hfill \\ \multicolumn{1}{c}{\mathbf{\Omega }=\arg \underset{{\mathit{u}}_{min}\le \mathit{u}\le {\mathit{u}}_{max}}{min}\Vert {\mathit{W}}_{\mathit{V}}(\mathit{BU}-\mathit{V}){\Vert }_2}\end{array}$$

(35)

The feasible solution field is $\mathit{\Omega }$ , where ${\mathit{W}}_{\mathit{U}}=\mathrm{diag}(1/F_{\mathrm{zij}})$ is the weight matrix that controls the output vector u . The greater the vertical tire force, the greater the tire distribution torque.

Here, ${\mathit{W}}_{\mathit{V}}=\mathbf{diag}(1,10)$ specifies that the demand for the yaw moment is preferentially satisfied if the demand for a longitudinal force cannot be satisfied, that is, to reduce the vehicle power and thus ensure the stability of the vehicle.

The SQP solution steps are more cumbersome. First, the feasible solution domain is guaranteed as a precondition for the optimal solution. This degrades the real-time performance of the algorithm. We introduce the weight coefficient $\epsilon$ to transform the problem into a more efficient WLS form, as shown in equation (36). ¹⁹ An accurate solution can be guaranteed when the penalty factor $\epsilon$ value is sufficiently large (e.g. $\epsilon =1000$ )

$$\mathit{U}=\arg \underset{{\mathit{u}}_{min}\le \mathit{u}\le {\mathit{u}}_{max}}{min}(\Vert {\mathit{W}}_{\mathit{U}}(\mathit{U}-{\mathit{U}}_d){\Vert }_2^2+\epsilon \Vert {\mathit{W}}_{\mathit{V}}(\mathit{BU}-\mathit{V}){\Vert }_2^2)$$

(36)

Equation (36) can be transformed into a standard least squares form using the deformation shown in equation (37), which can be solved through the least squares tool in the MATLAB optimization toolbox

$$\begin{array}{l}{\Vert {\mathit{W}}_{\mathit{U}}\left(\mathit{U}-{\mathit{U}}_d\right)\Vert }_2^2+\epsilon {\Vert {\mathit{W}}_{\mathit{V}}\left(\mathit{B}\mathit{U}-\mathit{V}\right)\Vert }_2^2\\ ={\Vert \underset{\mathit{A}}{\underbrace{\left(\begin{array}{c}\epsilon {\mathit{W}}_{\mathit{V}}\mathit{B}\\ {\mathit{W}}_{\mathit{U}}\end{array}\right)}}\mathit{U}-\underset{\mathit{b}}{\underbrace{\left(\begin{array}{c}\epsilon {\mathit{W}}_{\mathit{V}}\mathit{V}\\ {\mathit{W}}_{\mathit{U}}{\mathit{U}}_d\end{array}\right)}}\Vert }^2={\Vert \mathit{A}\mathit{U}-\mathit{b}\Vert }^2\end{array}$$

(37)

Simulation results

To validate the estimation results of the vehicle state estimator, we conducted a double line change (DLC) simulation test under a vehicle speed of 100 km/h and road adhesion coefficient of 0.8. The driver model for the steering wheel angle input was established in CarSim software using a preview time set to 1 s. The vehicle parameters were selected according to CarSim, as shown in Table 1. CarSim is a widely used vehicle simulation program and has high credibility, ensuring the reliability of the simulation test results. ²⁰ The DLC path was set as shown in Figure 4, and the results of the UKF are shown in Figures 5–7.

Table 1.

Simulation parameters of the vehicle.

Parameters	Description	Values
M	Vehicle mass	1480 kg
$J_{\omega }$	Wheel spin inertia	0.8 kg m2
$I_z$	Vehicle spin inertia	2200 kg m2
l	Wheelbase	2.5 m
a	Distance from CG to the front axle	1.3 m
b	Distance from CG to the rear axle	1.2 m
$l_f$	Front track width	1.6 m
$l_r$	Rear track width	1.5 m
$h_g$	Height of CG	0.5 m
R	Wheel rolling radius	0.32 m

CG: center of gravity.

Figure 4.

DLC test path traffic cones setting.

Figure 5.

Estimation results of vehicle condition parameter: (a) Vehicle longitudinal velocity, (b) vehicle lateral velocity, (c) Yaw rate and (d) sideslip angle.

Figure 6.

Estimation results of vehicle condition parameter: (a) Left front wheel lateral force, (b) Right front wheel lateral force, (c) Left rear wheel lateral force and (d) Right rear wheel lateral force.

Figure 7.

Estimation results of vehicle condition parameter: (a) Left front wheel vertical force, (b) Right front wheel vertical force, (c) Left rear wheel vertical force and (d) Right rear wheel vertical force.

Figure 5(a)–(d) show the vehicle speed, yaw rate, and sideslip angle estimation results of the vehicle state estimator, all of which are accurate. Figure 6(a)–(d) show the lateral tire force estimation results in which an error is shown at the peak, because the tires in this stage of force change violently and are highly nonlinear; the magic formula tire model for fitting the lateral force contains an error under these conditions. Figure 7(a)–(d) show the vertical force estimation results of the tire, where the proposed estimator is again highly accurate.

To verify the control effect of the active yaw moment controller and torque distribution controller, we conducted a DLC simulation test with a vehicle speed of 100 km/h and road adhesion coefficient of 0.4. The DLC path was set as shown in Figure 4. Three methods of torque distribution were tested. The active yaw moment controller output results are shown in Figure 8, the results of three methods of vehicle stability control are shown in Figure 9, and their resulting torque distribution are shown in Figure 10. “Method 1” indicates the NP, “Method 2” represents the QP, and “Method 3” is the WLS.

Figure 8.

Additional yaw moment controller outputs.

Figure 9.

Three optimization methods for stability control: (a) Control before and after vehicle sideslip angle. (b) Control before and after yaw rate. (c) Control before and after $\beta -\gamma$ phase plane.

Figure 10.

Three optimization methods torque distribution results: (a) Torque outputs under NP control. (b) Torque outputs under QP control. (c) Torque outputs under WLS control.

Figure 8 shows the output of the active yaw moment controller based on the SMC, which is one of the constraints of the torque distribution controller. Figure 9 shows the control results of the three vehicle drive torque distribution methods tested. We used two important parameters, yaw rate $\gamma$ and sideslip angle $\beta$ , to characterize the vehicle stability based on the $\beta -\gamma$ phase plane. The phase plane analysis method is commonly applied for similar purposes, such as the $\beta -\stackrel{·}{\beta }$ phase plane representing the sideslip angle and its changing trend, ²¹ or the $\beta -\gamma$ phase plane representing the sideslip angle and yaw rate.^22,23 The second method is particularly well applicable to investigating vehicle stability, because it includes two distinct parameters that intuitively characterize the stability. Figure 9(a) shows the results of the sideslip angle before and after the control, where the sideslip angle is larger and the stability is poorer when the vehicle is steering, whereas both are improved when using the control scheme. The three control methods discussed in this article perform well and are fairly similar. Figure 9(b) shows the yaw rate before and after the control, where such control has little effect on the steering intention of the driver, whereas the vehicle shows slight understeering characteristics indicative of enhanced stability. Figure 9(c) shows the results of the $\beta -\gamma$ phase plane before and after the control; again, the plane demonstrates an enhanced stability after the control scheme is imposed.

In the following, we will analyze the characteristics of three methods. As we all know, variance is a useful method to represent data fluctuations. However, the mean of the three methods has a large difference, so it is not advisable to directly compare the variance of the data. Thus, we have adopted another method to represent the fluctuation of data. First is the pairwise combination of the data of three methods, followed by the calculation of the differences of these combinations and, finally, calculation of the variances or means of the above differences. The calculation method is as follows

$$\psi \left(T\right)=\frac{\sum _{m=1}^4\left\{\mathrm{var}[T_{\mathrm{dm},q}-T_{\mathrm{dm},p}]\right\}}{4}$$

(38)

$$\{\begin{array}{c}\varphi \left(\gamma \right)=\mathrm{mean}\left(|{\gamma }_q-{\gamma }_p|\right)\hfill \\ \multicolumn{1}{c}{\xi \left(\beta \right)=\mathrm{mean}\left(|{\beta }_q-{\beta }_p|\right)}\end{array}$$

(39)

where $\psi (T)$ denotes the variance of motor torque error, $\varphi (\gamma )$ represents the mean of yaw rate error, $\xi (\beta )$ denotes the mean of sideslip angle error, the subscript $m=\mathrm{fl},\mathrm{fr},\mathrm{rl},\mathrm{rr}$ , and var denotes variance. The setting of subscripts q and p is shown in Table 2.

Table 2.

Comparison of three optimization methods.

	q = QP, p = NP	q = WLS, p = NP	q = WLS, p = QP
$\psi \left(T\right)$	0.2532	4.5197	4.5799
$\varphi \left(\gamma \right)$	2.30E-07	1.48E-07	3.79E-07
$\xi \left(\beta \right)$	2.76E-06	6.88E-06	4.12E-06

QP: quadratic programming; NP: nonlinear programming; WLS: weighted least squares.

All three methods perform similarly in providing effective vehicle stability control as shown in Figure 9 and $\varphi (\gamma ),\xi (\beta )$ in Table 2. The $\varphi (\gamma )$ and $\xi (\beta )$ are between 10⁻⁷ and 10⁻⁵; this represents a small difference. Figure 10 shows the results of the three torque distribution optimization methods. Among them, the fluctuation levels of the NP and QP algorithms are smaller as shown in Figure 10 and $\psi (T)$ in Table 2. These two algorithms reflect well the trend of the additional yaw moment and yield with relatively stable curves. The WLS method is less stable and fluctuates significantly when disturbed, because it is largely different from the other two methods. This can be clearly seen from Figure 10(c) and $\psi (T)$ in Table 2. The motor will be shocked too much. However, this method also comes with a lower calculation time cost. The calculation time of three solution methods in MATLAB, as shown in Table 3, is as follows: the calculation time of WLS is 432.10 s, QP algorithm is 530.74 s, and NP method is 1036.86 s.

Table 3.

Run time of three methods.

Method	Run time ( $\mathrm{s}$ )
NP	1036.86
QP	530.74
WLS	432.10

QP: quadratic programming; NP: nonlinear programming; WLS: weighted least squares.

Conclusion

1. This article presented a novel vehicle stability controller comprised of a vehicle state estimator, active yaw moment controller, and torque distribution controller. The vehicle state estimator was designed based on a UKF, the active yaw moment controller was designed based on the SMC theory, and the torque distribution controller was based on three mathematical optimization algorithms.

2. The UKF-based vehicle state estimator was shown to accurately estimate such parameters as the vehicle speed, sideslip angle, yaw rate, and tire force. This lays the foundation for acquisition of the active yaw moment control and torque distribution control information.

3. We simulated a DLC using a high-speed and low adhesion coefficient under which the proposed controller performed extremely well.

4. The control output results of the NP, QP algorithm, and WLS method were compared, and it was determined that all three are feasible but with different advantages and disadvantages. The NP and QP showed the best stability, whereas the WLS method showed the fastest calculation time.

In a future work, we will apply hardware-in-the-loop testing on the algorithm to verify the real-time performance of the control method in the C programming language environment.