四轮独立驱动电动汽车轨迹跟踪及稳定性协调控制

摘要

针对紧急避障轨迹跟踪和车辆稳定性多目标控制问题，以四轮独立驱动智能电动汽车为研究对象，提出一种分层结构的轨迹跟踪及稳定性控制方法，上层控制器采用线性时变模型预测控制（LTV MPC）生成期望前轮转角和附加横摆力矩，并通过PID速度跟踪将车速变化考虑到模型预测优化求解中生成总需求转矩，下层控制器采用二次规划将上层计算的广义力最优分配给四个车轮。其中，预测模型为8自由度车辆模型，被控对象为14自由度车辆模型，并建立联合工况刷子轮胎模型。不同车速，路面附着条件以及有无稳定性控制的双移线工况仿真结果表明，所提控制方法具有较好的鲁棒性和轨迹跟踪性能，极限工况下在保证车辆稳定性的同时提高了轨迹跟踪精度。

关键词

四轮独立驱动电动汽车; 轨迹跟踪; 横摆稳定性; 模型预测控制; 控制分配

自动驾驶车辆是集成了环境感知、决策规划和运动控制等先进智能技术的移动装备，目前大部分自动驾驶采用线控底盘电动车辆^［1］。近年，自动驾驶紧急避撞的安全性受到极大关注，对于紧急避撞的研究，大多集中在提高路径跟踪的精度，然而，在一些复杂驾驶场景，车辆在避撞时，紧急转向可能造成车辆失稳，发生侧滑或侧翻。因此，需要将车辆稳定性控制考虑到紧急避撞控制策略中。分布式驱动电车四轮转矩独立可控，通过两侧车轮转矩差生成附加横摆力矩保证车辆稳定性，成为解决该问题的一种有效途径^［2］。

自动驾驶车辆路径跟踪的目标是以最小横向偏差和航向偏差跟踪期望路径，其本质是求解多约束条件下的优化问题。目前提出了许多算法，如纯跟踪方法，Stanley方法，最优预瞄控制方法，PID控制，LQR方法，SMC方法以及$H_{\mathrm{\infty }}$鲁棒控制方法等^［3-5］，然而轨迹跟踪和车辆稳定性之间可能会因控制目标冲突而使系统性能恶化，影响车辆的稳定性。因此，需要将轨迹跟踪和稳定性考虑到统一控制框架中，实现无人驾驶车辆轨迹稳定跟踪^［6］。文献［5］提出分层架构的模型预测轨迹跟踪与横摆稳定性的协调控制策略，上层为模型预测运动跟踪控制，下层为纵向力优化分配控制，仿真验证了所提控制器在不同车速不同附着条件下的鲁棒性。文献［7］针对路径跟踪性能和车辆稳定性多目标控制问题，提出分层的无人车轨迹跟踪和稳定性协调控制方法，上层采用滑模控制减小路径跟踪航向偏差和横向偏差，同时保证车辆稳定性，下层根据上层需求设计四轮轮胎力优化分配方法。文献［8］提出综合横向和纵向控制的分层控制结构，上层控制采用非线性模型预测控制实现参考轨迹和纵向速度的跟踪，下层控制考虑电机力矩输出和路面附着约束，将轮胎力分配问题转化为带约束的非线性规划问题，实现四轮轮胎力的最优分配。直接横摆力矩控制（DYC）通过精确快速调节电机力矩，产生附加横摆力矩，提高车辆的操纵稳定性，因此，融合DYC考虑横摆稳定性的路径跟踪控制也逐步展开^［9］。文献［10］提出一种分层的路径跟踪控制框架，考虑车辆参数和外部干扰的不确定性，上层采用约束跟随的自适应鲁棒控制方法，下层采用二次规划控制分配方法，实现前轮转向与直接横摆力矩控制的协同控制，提高了路径跟踪的精度和鲁棒性。

上述方法大多基于简化单轨车辆模型设计上层控制器，较少考虑车辆侧倾与横摆运动的耦合，此外，在下层控制分配中考虑的目标函数大多也比较单一，对于多目标协同优化较少考虑。为进一步探讨线性时变模型预测控制在考虑整车操纵稳定性的无人驾驶车辆紧急避障轨迹跟踪控制上的应用，本文以四轮独立驱动电动车辆为研究对象，提出一种分层架构的模型预测轨迹跟踪及稳定性控制方法，上层控制基于LTV MPC生成期望前轮转角及附加横摆力矩，其中，预测模型为8自由度车辆模型，被控对象为14自由度车辆模型，轮胎模型采用纵滑侧偏联合工况下的刷子轮胎模型，以适应复杂极限工况的车辆状态预测和仿真，纵向控制采用PID速度跟踪方法，通过嵌入MPC框架将纵向速度变化考虑到模型预测优化求解中，下层控制采用二次规划将上层控制计算的总需求转矩和附加横摆力矩最优分配给四个车轮，结合前轮转角输入，实现无人驾驶车辆的紧急避障轨迹跟踪及横摆稳定性控制。

1 车辆动力学建模

1.1　轮胎模型

极限工况下轮胎存在较大滑动又存在侧偏的情况，其纵向力和侧向力相互耦合，并满足附着椭圆约束。刷子轮胎模型将轮胎表示为可变形刷子包围的刚性环，利用物理参数将轮胎的特征表达，常用来表征轮胎的非线性特性，其表达式相比魔术公式简洁，且有一定的精度，纵滑侧偏联合工况下刷子轮胎模型如下^［11］：

$$F_{x,i}=\frac{C_x\left(\frac{{\kappa }_i}{\mathrm{1}+{\kappa }_i}\right)}{f_i}{F}_i$$

（1）

$$F_{y,i}=-\frac{C_{\alpha }\left(\frac{\mathrm{t}\mathrm{a}\mathrm{n}{\alpha }_i}{\mathrm{1}+{\kappa }_i}\right)}{f_i}{F}_i$$

（2）

其中：

$$F_i=\left\{\begin{array}{c}{f}_i-\frac{\mathrm{1}}{\mathrm{3}\mu F_{z,i}}{f}_i^{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{27}{\mu }^{\mathrm{2}}{F}_{z,i}^{\mathrm{2}}}{f}_i^{\mathrm{3}}\mathrm{i}\mathrm{f}{f}_i\le \mathrm{3}\mu F_{z,i}\\ \mu F_{z,i}\mathrm{其他}\end{array}\right.$$

（3）

$$f_i=\sqrt[]{C_x^{\mathrm{2}}{\left(\frac{{\kappa }_i}{\mathrm{1}+{\kappa }_i}\right)}^{\mathrm{2}}+C_{\alpha }^{\mathrm{2}}{\left(\frac{\mathrm{t}\mathrm{a}\mathrm{n}{\alpha }_i}{\mathrm{1}+{\kappa }_i}\right)}^{\mathrm{2}}}$$

（4）

式中： ${\kappa }_i$为车轮的纵向滑移率。考虑到非线性轮胎模型计算量大，不利于控制实时性，当滑移率和侧偏角都限制在较小的值时，常采用线性轮胎模型模拟线性范围内的轮胎侧向力和纵向力。在小侧偏角和滑移率的假设下，轮胎侧向力可表示为：$F_{yt}=C_{\alpha }(\mu ,F_z)\alpha$，其中$C_{\alpha }$称为轮胎的侧偏刚度，与路面附着系数和垂向力有关；纵向力可表示为：$F_{xt}=C_x(\mu ,F_z)s_{\mathrm{f},\mathrm{r}}$，其中$C_x$是纵向刚度，也与路面附着系数和垂向力有关^［12］。

1.2　8自由度车辆模型

8自由度的车辆模型常作为简化的低阶模型，研究无大幅度纵向加速度的车辆操纵。8自由度的车辆模型考虑了车体纵向，横向，横摆和侧倾动力学以及四个车轮的转动动力学，俯仰运动未建模，前后悬架由等效侧倾刚度和侧倾阻尼系数表示，其受力分析如图1所示^［12-14］：

$$m\left(\dot{u}-{\omega }_{z}v\right)=\sum F_{x\mathrm{g}\mathrm{i}\mathrm{j}}+\left(m_{\mathrm{u}\mathrm{f}}a-m_{\mathrm{u}\mathrm{r}}b\right){\omega }_z^{\mathrm{2}}-\mathrm{2}{h}_{\mathrm{r}\mathrm{c}}{m}_{\mathrm{s}}{\omega }_z{\omega }_x$$

（5）

$$m\left(\dot{v}+{\omega }_{z}u\right)=\sum F_{y\mathrm{g}\mathrm{i}\mathrm{j}}+\left(m_{\mathrm{u}\mathrm{r}}b-m_{\mathrm{u}\mathrm{f}}a\right){\dot{\omega }}_z+h_{\mathrm{r}\mathrm{c}}{m}_{\mathrm{s}}{\dot{\omega }}_x$$

（6）

$$I_z{\dot{\omega }}_z+I_{xz}{\dot{\omega }}_x=\left(F_{y\mathrm{g}\mathrm{l}\mathrm{f}}+F_{y\mathrm{g}\mathrm{r}\mathrm{f}}\right)a-\left(F_{y\mathrm{g}\mathrm{l}\mathrm{r}}+F_{y\mathrm{g}\mathrm{r}\mathrm{r}}\right)b+\frac{\left(F_{x\mathrm{g}\mathrm{r}\mathrm{f}}-F_{x\mathrm{g}\mathrm{l}\mathrm{f}}\right)c_{\mathrm{f}}}{\mathrm{2}}+\frac{\left(F_{x\mathrm{g}\mathrm{r}\mathrm{r}}-F_{x\mathrm{g}\mathrm{l}\mathrm{r}}\right)c_{\mathrm{r}}}{\mathrm{2}}+\left(m_{\mathrm{u}\mathrm{r}}b-m_{\mathrm{u}\mathrm{f}}a\right)\left(\dot{v}+{\omega }_{z}u\right)$$

（7）

$$\left(I_x+m_{\mathrm{s}}{h}_{\mathrm{r}\mathrm{c}}^{\mathrm{2}}\right){\dot{\omega }}_x+I_{xz}{\dot{\omega }}_z=m_{\mathrm{s}}g{h}_{\mathrm{r}\mathrm{c}}\varphi -\left(k_{\varphi \mathrm{f}}+k_{\varphi \mathrm{r}}\right)\varphi -\left(b_{\varphi \mathrm{f}}+b_{\varphi \mathrm{r}}\right)\dot{\varphi }+h_{\mathrm{r}\mathrm{c}}{m}_{\mathrm{s}}\left(\dot{v}+{\omega }_{z}u\right)$$

（8）

式中： $h_{\mathrm{r}\mathrm{c}}=\frac{h_{\mathrm{r}\mathrm{c}\mathrm{f}}b+h_{\mathrm{r}\mathrm{c}\mathrm{r}}a}{a+b}$； $u$、$v$分别为纵向车速和横向车速； ${\omega }_z$、${\omega }_x$分别为横摆角速度和侧倾角速度； $\varphi$为车身侧倾角； $I_x$、$I_z$分别为车辆绕X轴和绕Z轴转动惯量； $I_{xz}$为侧倾横摆的惯性积； $m$为总质量； $m_{\mathrm{s}}$为簧载质量； $m_{\mathrm{u}\mathrm{f}}$、$m_{\mathrm{u}\mathrm{r}}$分别为前后非簧载质量； $F_{xgij}$、$F_{ygij}$分别为车辆坐标系下轮胎纵向力和侧向力，下标ij组合表示前右（rf）、前左（lf）、后右（rr）、后左（lr）； $a$、$b$分别为重心距前后轴距离； $c_{\mathrm{f}}$、$c_{\mathrm{r}}$分别为前后轮距； $h_{\mathrm{r}\mathrm{c}\mathrm{f}}$、$h_{\mathrm{r}\mathrm{c}\mathrm{r}}$分别为前后侧倾中心距簧载质量的距离； $k_{\varphi \mathrm{f}}$、$k_{\varphi \mathrm{r}}$分别为前后悬架等效侧倾刚度； $b_{\varphi \mathrm{f}}$、$b_{\varphi \mathrm{r}}$分别为前后悬架等效侧倾阻尼； $\mathrm{g}$为重力加速度。附加横摆力矩$∆M_z$来自4个车轮纵向力差异，车轮侧向力和纵向力经坐标转化后，式（7）可表示为：

$$I_z{\dot{\omega }}_z+I_{xz}{\dot{\omega }}_x=\left(F_{y\mathrm{t}\mathrm{l}\mathrm{f}}\mathrm{c}\mathrm{o}\mathrm{s}\delta +F_{y\mathrm{t}\mathrm{r}\mathrm{f}}\mathrm{c}\mathrm{o}\mathrm{s}\delta \right)a-\left(F_{y\mathrm{t}\mathrm{l}\mathrm{r}}+F_{y\mathrm{t}\mathrm{r}\mathrm{r}}\right)b+\frac{\left(F_{y\mathrm{t}\mathrm{l}\mathrm{f}}\mathrm{s}\mathrm{i}\mathrm{n}\delta -F_{y\mathrm{t}\mathrm{r}\mathrm{f}}\mathrm{s}\mathrm{i}\mathrm{n}\delta \right)c_{\mathrm{f}}}{\mathrm{2}}+\left(m_{\mathrm{u}\mathrm{r}}b-m_{\mathrm{u}\mathrm{f}}a\right)\left(\dot{v}+{\omega }_{z}u\right)+∆M_z$$

（9）

式中： $∆M_z=\left(F_{x\mathrm{t}\mathrm{l}\mathrm{f}}\mathrm{s}\mathrm{i}\mathrm{n}\delta +F_{x\mathrm{t}\mathrm{r}\mathrm{f}}\mathrm{s}\mathrm{i}\mathrm{n}\delta \right)a+\frac{\left(F_{x\mathrm{t}\mathrm{r}\mathrm{f}}\mathrm{c}\mathrm{o}\mathrm{s}\delta -F_{x\mathrm{t}\mathrm{l}\mathrm{f}}\mathrm{c}\mathrm{o}\mathrm{s}\delta \right)c_{\mathrm{f}}}{\mathrm{2}}+\frac{\left(F_{x\mathrm{t}\mathrm{r}\mathrm{r}}-F_{x\mathrm{t}\mathrm{l}\mathrm{r}}\right)c_{\mathrm{r}}}{\mathrm{2}}$； $F_{x\mathrm{t}\mathrm{r}\mathrm{f}}$、$F_{y\mathrm{t}\mathrm{i}\mathrm{j}}$ 分别表示车轮坐标系下的纵向力和侧向力。此外， 考虑前轮驱动车辆，以前后右侧车轮为例，车轮的旋转动力学方程表示为^［13］：

$$\begin{array}{c}\begin{array}{l}{J}_w{\dot{\omega }}_{\mathrm{r}\mathrm{f}}=T_{d\mathrm{r}\mathrm{f}}-T_{b\mathrm{r}\mathrm{f}}-r_{\mathrm{r}\mathrm{f}}{F}_{xt\mathrm{r}\mathrm{f}}\\ J_w{\dot{\omega }}_{\mathrm{r}\mathrm{r}}=-T_{b\mathrm{r}\mathrm{r}}-r_{\mathrm{r}\mathrm{r}}{F}_{xt\mathrm{r}\mathrm{r}}\end{array}\end{array}$$

（10）

式中： $T_{d\mathrm{r}\mathrm{f}}$是传递给车轮的驱动力矩； $T_{b\mathrm{r}\mathrm{f}}$是车轮的制动力矩； $r_{\mathrm{r}\mathrm{f}}$和$r_{\mathrm{r}\mathrm{r}}$分别是前后轮的滚动半径； $J_w$是车轮的转动惯量。

图1  8自由度车辆动力学模型

Fig.1  Schematic of the 8-DOF vehicle model

1.3　14自由度车辆模型

14自由度车辆模型在8自由度车辆模型的基础上考虑了一维悬架运动，还可预测车辆的俯仰运动，提供了建立非线性弹簧和减振器的可能性，因此可以用于测试防侧翻控制策略的有效性。也就是说，虽然14自由度模型相比更复杂度多体动力学模型，自由度较少且有一定限制，但该模型可以足够表达在大部分控制系统设计中重要的车辆运动。14自由度车辆动力学模型如图2所示，该模型在之前的研究工作和相关文献中已建立，本文不再详述，见文献^［12-15］。

图2  14自由度车辆动力学模型

Fig.2  Schematic of 14-DOF vehicle model

2 轨迹跟踪及稳定性控制器设计

本文所提的分层轨迹跟踪及稳定性控制架构如图3所示，上层控制采用线性时变模型预测控制实现对期望轨迹和期望运动的跟踪，考虑控制输入和系统状态输出的约束，将轨迹跟踪和稳定性控制统一考虑在模型预测控制框架中；下层控制为转矩分配层，根据上层计算的附加横摆力矩和总的需求转矩，以最小化转矩跟踪误差、轮胎利用率以及纵向滑移率为优化目标，并考虑路面附着极限和电机最大峰值转矩，通过二次规划方法将上层广义力最优分配给4个车轮。

图3  分层的轨迹跟踪及稳定性统一控制架构

Fig.3  Hierarchical coordinated control strategy structure of trajectory tracking and yaw stability

2.1　车辆参考运动状态

本文基于2自由度单轨车辆模型生成理想横摆角速度和质心侧偏角，基于侧向力和侧偏角之间的线性关系假设，忽略轮胎侧向力的非线性特性，可得^［16-17］：

$${\dot{\psi }}_{\mathrm{s}}=\frac{u/L}{\mathrm{1}+K{u}^{\mathrm{2}}}{\delta }_{\mathrm{f}}$$

（11）

$${\beta }_{\mathrm{s}}=\frac{b+ma{u}^{\mathrm{2}}/(C_{\alpha r}L)}{L\left(\mathrm{1}+K{u}^{\mathrm{2}}\right)}{\delta }_{\mathrm{f}}$$

（12）

其中：$K=\frac{m}{L^{\mathrm{2}}}\left(\frac{a}{C_{\alpha r}}-\frac{b}{C_{\alpha f}}\right)$为稳定性因子。考虑路面附着系数约束，则理想横摆角速度上限约束为^［16-17］：

$${\dot{\psi }}_{\mathrm{r}\mathrm{e}\mathrm{f}}=\mathrm{m}\mathrm{i}\mathrm{n}\left\{\left|{\dot{\psi }}_{\mathrm{s}}\right|,\left|\mathrm{0.85}\frac{\mu g}{u}\right|\right\}·\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\left(\delta \right)$$

（13）

参考质心侧偏角需限制其上限防止该值过大，因为在大质心侧偏角时，轮胎将进入非线性，并接近附着极限，对上限的经验值一般设置为arctan （0.02$\mu g$）。因此，为考虑路面附着极限并补偿非线性特性，对名义质心侧偏角进行修正，得到^［17］：

$${\beta }_{\mathrm{r}\mathrm{e}\mathrm{f}}=\mathrm{m}\mathrm{i}\mathrm{n}\left\{\begin{array}{c}\left|{\beta }_{\mathrm{s}}\right|,\left|\left(\frac{b}{u^{\mathrm{2}}}+\frac{ma}{C_{\alpha r}L}\right)\mu g\right|\\ ,\left|\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}\left(\mathrm{0.02}\mu g\right)\right|\end{array}\right\}·\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\left(\delta \right)$$

（14）

2.2　基于LTV MPC的轨迹跟踪及稳定性控制

上层控制实现对期望运动状态和期望轨迹的跟踪，生成前轮转角控制输入和附加横摆力矩输入。前轮转角和附加横摆力矩通过模型预测控制算法求解，并将PID速度跟踪控制考虑到MPC的优化求解中，在每次优化求解时，更新预测模型的纵向速度。

2.2.1　车辆模型线性化

由于控制器内部预测模型为非线性，采用非线性模型预测控制（NMPC）的预测控制算法求解困难且实时性差，因此对非线性预测模型进行线性化处理，采用线性模型预测控制进行轨迹跟踪。

基于非线性8自由度车辆模型，设系统状态量为$\mathit{\chi }={\left[\dot{y},\dot{x},\varphi ,\psi ,\dot{\varphi },\dot{\psi },Y,X\right]}^{\mathrm{T}}$，控制量为车轮转角和附加横摆力矩$\mathit{u}={\left[\delta ,∆M_z\right]}^{\mathrm{T}}$，选取系统的输出量为车辆的航向角、横向位置、纵向位置、横摆角速度和质心侧偏角，即$\mathit{\eta }={\left[\psi ,Y,X,\dot{\psi },\beta \right]}^{\mathrm{T}}$。

车辆动力学模型的一般表示形式如下：

$$\dot{\mathit{\chi }}=f(\mathit{\chi },\mathit{u})$$

（15）

设系统在工作点$\left({\mathit{\chi }}_{\mathrm{o}},{\mathit{u}}_{\mathrm{o}}\right)$处高阶可微分，通过在工作点处泰勒级数展开，并忽略高阶项，则^［12-14，18］

$$\dot{\mathit{\chi }}=f\left({\mathit{\chi }}_{\mathrm{o}},{\mathit{u}}_{\mathrm{o}}\right)+{\left.\frac{\partial f(\mathit{\chi },\mathit{u})}{\partial \mathit{\chi }}\right|}_{\begin{array}{c}\mathit{\chi }={\mathit{\chi }}_{\mathrm{o}}\\ \mathit{u}={\mathit{u}}_{\mathrm{o}}\end{array}}\left(\mathit{\chi }-{\mathit{\chi }}_{\mathrm{o}}\right)+{\left.\frac{\partial f(\mathit{\chi },\mathit{u})}{\partial \mathit{u}}\right|}_{\begin{array}{c}\mathit{\chi }={\mathit{\chi }}_{\mathrm{o}}\\ \mathit{u}={\mathit{u}}_{\mathrm{o}}\end{array}}\left(\mathit{u}-{\mathit{u}}_{\mathrm{o}}\right)$$

（16）

从而得到线性误差模型^［12-14，18］：

$$\dot{\tilde{\mathit{\chi }}}=\mathit{A}\tilde{\mathit{\chi }}+\mathit{B}\tilde{\mathit{u}}$$

（17）

式中： $\dot{\tilde{\mathit{\chi }}}=\dot{\mathit{\chi }}-{\dot{\mathit{\chi }}}_{\mathrm{o}}$， $\tilde{\mathit{\chi }}=\mathit{\chi }-{\mathit{\chi }}_{\mathrm{o}}$， $\tilde{\mathit{u}}=\mathit{u}-{\mathit{u}}_{\mathrm{o}}$； $\mathit{A}$和$\mathit{B}$分别为系统矩阵和输入矩阵。

为将线性化模型应用于MPC控制器的设计，将其转化为如下离散状态空间形式^［12-14，18］：

$$\mathit{\chi }(k+\mathrm{1})={\mathit{A}}_d\mathit{\chi }(k)+{\mathit{B}}_d\mathit{u}(k)+{\mathit{d}}_k(k)$$

（18）

式中：${\mathit{A}}_d=\mathit{I}+T\mathit{A}$；${\mathit{B}}_d=T\mathit{B}$；${\mathit{d}}_k(k)=f({\mathit{\chi }}_{\mathrm{o}}(k),{\mathit{u}}_{\mathrm{o}}(k))-\left({\mathit{A}}_d{\mathit{\chi }}_{\mathrm{o}}(k)+{\mathit{B}}_d{\mathit{u}}_{\mathrm{o}}(k)\right)$； $T$为采样时间。

8自由度非线性车辆模型系统矩阵和输入矩阵分别为：

$$\mathit{A}={\left.\frac{\partial f(\mathit{\chi },\mathit{u})}{\partial \mathit{\chi }}\right|}_{\begin{array}{c}\mathit{\chi }={\mathit{\chi }}_{\mathrm{o}}\\ \mathit{u}={\mathit{u}}_{\mathrm{o}}\end{array}}=\left[\begin{array}{cccccccc}\frac{\partial f_{\dot{y}}}{\partial \dot{y}}& \frac{\partial f_{\dot{y}}}{\partial \dot{x}}& \frac{\partial f_{\dot{y}}}{\partial \varphi }& \mathrm{0}& \frac{\partial f_{\dot{y}}}{\partial \dot{\varphi }}& \frac{\partial f_{\dot{y}}}{\partial \dot{\psi }}& \mathrm{0}& \mathrm{0}\\ \frac{\partial f_{\dot{x}}}{\partial \dot{y}}& \frac{\partial f_{\dot{x}}}{\partial \dot{x}}& \frac{\partial f_{\dot{x}}}{\partial \varphi }& \mathrm{0}& \frac{\partial f_{\dot{x}}}{\partial \dot{\varphi }}& \frac{\partial f_{\dot{x}}}{\partial \dot{\psi }}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}\\ \frac{\partial f_{\dot{\varphi }}}{\partial \dot{y}}& \frac{\partial f_{\dot{\varphi }}}{\partial \dot{x}}& \frac{\partial f_{\dot{\varphi }}}{\partial \varphi }& \mathrm{0}& \frac{\partial f_{\dot{\varphi }}}{\partial \dot{\varphi }}& \frac{\partial f_{\dot{\varphi }}}{\partial \dot{\psi }}& \mathrm{0}& \mathrm{0}\\ \frac{\partial f_{\dot{\psi }}}{\partial \dot{y}}& \frac{\partial f_{\dot{\psi }}}{\partial \dot{x}}& \frac{\partial f_{\dot{\psi }}}{\partial \varphi }& \mathrm{0}& \frac{\partial f_{\dot{\psi }}}{\partial \dot{\varphi }}& \frac{\partial f_{\dot{\psi }}}{\partial \dot{\psi }}& \mathrm{0}& \mathrm{0}\\ \mathrm{c}\mathrm{o}\mathrm{s}\psi & \mathrm{s}\mathrm{i}\mathrm{n}\psi & \mathrm{0}& \dot{x}\mathrm{c}\mathrm{o}\mathrm{s}\psi -\dot{y}\mathrm{s}\mathrm{i}\mathrm{n}\psi & \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ -\mathrm{s}\mathrm{i}\mathrm{n}\psi & \mathrm{c}\mathrm{o}\mathrm{s}\psi & \mathrm{0}& -\dot{y}\mathrm{c}\mathrm{o}\mathrm{s}\psi -\dot{x}\mathrm{s}\mathrm{i}\mathrm{n}\psi & \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right]$$

$$\mathit{B}={\left.\frac{\partial f(\mathit{\chi },\mathit{u})}{\partial \mathit{u}}\right|}_{\begin{array}{c}\mathit{\chi }={\mathit{\chi }}_{\mathrm{o}}\\ \mathit{u}={\mathit{u}}_{\mathrm{o}}\end{array}}=\left[\begin{array}{cccccccc}\frac{\partial f_{\dot{y}}}{\partial \delta }& \frac{\partial f_{\dot{x}}}{\partial \delta }& \mathrm{0}& \mathrm{0}& \frac{\partial f_{\dot{\varphi }}}{\partial \delta }& \frac{\partial f_{\dot{\psi }}}{\partial \delta }& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \frac{\partial f_{\dot{x}}}{\partial ∆M_z}& \mathrm{0}& \mathrm{0}& \frac{\partial f_{\dot{\varphi }}}{\partial ∆M_z}& \frac{\partial f_{\dot{\psi }}}{\partial ∆M_z}& \mathrm{0}& \mathrm{0}\end{array}\right]$$

2.2.2　车辆状态预测

为了考虑控制输入增量，定义新的系统状态变量$\mathit{\xi }(k)=\left[\begin{array}{c}\mathit{\chi }(k)\\ \mathit{u}(k-\mathrm{1})\end{array}\right]$，系统的输出变量为$\mathit{\eta }(\mathrm{k})$，控制输入增量为$\mathrm{\Delta }\mathit{u}(k)=\mathit{u}(k)-\mathit{u}(k-\mathrm{1})$，则新的离散化的状态空间控制器可表示为^［12-14，18］：

$$\mathit{\xi }(k+\mathrm{1})={\tilde{\mathit{A}}}_d\mathit{\xi }(k)+{\tilde{\mathit{B}}}_d\mathrm{\Delta }\mathit{u}(k)+{\tilde{\mathit{d}}}_k(k)\mathit{\eta }(k)={\tilde{\mathit{C}}}_d\mathit{\xi }(k)$$

（19）

式中： ${\tilde{\mathit{A}}}_d=\left[\begin{array}{cc}{\mathit{A}}_d& {\mathit{B}}_d\\ {\mathrm{0}}_{m\times n}& {\mathit{I}}_m\end{array}\right];{\tilde{\mathit{B}}}_d=\left[\begin{array}{c}{\mathit{B}}_d\\ {\mathit{I}}_m\end{array}\right];{\tilde{\mathit{d}}}_k\left(k\right)=\left[\begin{array}{c}{\mathit{d}}_k\left(k\right)\\ {\mathrm{0}}_m\end{array}\right];{\tilde{\mathit{C}}}_d=\left[\begin{array}{c}{\mathit{C}}_d\\ {\mathrm{0}}_{p\times m}\end{array}\right]$（其中，m表示控制输入维度，n表示状态变量的维度，p表示输出变量的维度。）

在预测时域$N_p$上的系统状态输出矩阵可表示为^［12-14，18］：

$${\mathit{\eta }}_m(k)={\mathit{\Theta }}_m\mathit{\xi }(k)+{\mathit{\Gamma }}_m\Delta U_m+{\mathit{\Psi }}_m{D}_k$$

（20）

式中：

$${\mathit{\eta }}_m={\left[\begin{array}{ccc}\mathit{\eta }\left(k+\mathrm{1}\right)& \dots & \mathit{\eta }\left(k+N_p\right)\end{array}\right]}^{\mathrm{T}}$$

$${\mathit{\Theta }}_m={\left[\begin{array}{cccccc}{\tilde{\mathit{C}}}_d{\tilde{\mathit{A}}}_d& {\tilde{\mathit{C}}}_d{\tilde{\mathit{A}}}_d^{\mathrm{2}}& \dots & {\tilde{\mathit{C}}}_d{\tilde{\mathit{A}}}_d^{N_c}& \dots & {\tilde{\mathit{C}}}_d{\tilde{\mathit{A}}}_d^{N_p}\end{array}\right]}^{\mathrm{T}}$$

$$\mathrm{\Delta }{\mathit{U}}_m={\left[\begin{array}{ccccc}\mathrm{\Delta }\mathit{u}(k)& \dots & \mathrm{\Delta }\mathit{u}(k+m)& \dots & \mathrm{\Delta }\mathit{u}(k+N_c-\mathrm{1})\end{array}\right]}^{\mathrm{T}}$$

$${\mathit{D}}_k={\left[\begin{array}{cccc}{\tilde{\mathit{d}}}_k(k)& {\tilde{\mathit{d}}}_k(k+\mathrm{1})& \dots & {\tilde{\mathit{d}}}_k(k+N_p-\mathrm{1})\end{array}\right]}^{\mathrm{T}}$$

$${\mathit{\Gamma }}_m=\left[\begin{array}{cccc}{\tilde{\mathit{C}}}_d{\tilde{\mathit{B}}}_d& \mathrm{0}& \cdots & \mathrm{0}\\ {\tilde{\mathit{C}}}_d{\tilde{\mathit{A}}}_d{\tilde{\mathit{B}}}_d& {\tilde{\mathit{C}}}_d{\tilde{\mathit{B}}}_d& \cdots & \mathrm{0}\\ ⋮& ⋮& \cdots & \mathrm{0}\\ {\tilde{\mathit{C}}}_d{\tilde{\mathit{A}}}_d^{N_c-\mathrm{1}}{\tilde{\mathit{B}}}_d& {\tilde{\mathit{C}}}_d{\tilde{\mathit{A}}}_d^{N_c-\mathrm{2}}{\tilde{\mathit{B}}}_d& \cdots & {\tilde{\mathit{C}}}_d{\tilde{\mathit{B}}}_d\\ ⋮& ⋮& \cdots & ⋮\\ {\tilde{\mathit{C}}}_d{\tilde{\mathit{A}}}_d^{N_p-\mathrm{1}}{\tilde{\mathit{B}}}_d& {\tilde{\mathit{C}}}_d{\tilde{\mathit{A}}}_d^{N_p-\mathrm{2}}{\tilde{\mathit{B}}}_d& \cdots & {\tilde{\mathit{C}}}_d{\tilde{\mathit{A}}}_d^{N_p-N_c}{\tilde{\mathit{B}}}_d\end{array}\right]$$

$${\mathit{\Psi }}_m=\left[\begin{array}{ccccc}{\tilde{\mathit{C}}}_d& \mathrm{0}& \mathrm{0}& \cdots & \mathrm{0}\\ {\tilde{\mathit{C}}}_d{\tilde{\mathit{A}}}_d& {\tilde{\mathit{C}}}_d& \mathrm{0}& \cdots & \mathrm{0}\\ {\tilde{\mathit{C}}}_d{\tilde{\mathit{A}}}_d^{\mathrm{2}}& {\tilde{\mathit{C}}}_d{\tilde{\mathit{A}}}_d& {\tilde{\mathit{C}}}_d& \cdots & \mathrm{0}\\ ⋮& ⋮& ⋮& \cdots & ⋮\\ {\tilde{\mathit{C}}}_d{\tilde{\mathit{A}}}_d^{N_p-\mathrm{1}}& {\tilde{\mathit{C}}}_d{\tilde{\mathit{A}}}_d^{N_p-\mathrm{2}}& {\tilde{\mathit{C}}}_d{\tilde{\mathit{A}}}_d^{N_p-\mathrm{3}}& \cdots & {\tilde{\mathit{C}}}_d\end{array}\right]$$

2.2.3　目标函数设计

为使无人驾驶车辆准确跟踪期望轨迹和期望运动状态，保持车辆的横摆稳定性，控制器在每个控制周期求解滚动时域优化问题。采用二次型目标函数，在每个有限时域中，求解最优控制输入增量，进而得到总的需求转矩和期望附加横摆力矩，设计如下目标函数^［12-14，18］：

$$\underset{\mathrm{\Delta }{U}_t,\epsilon }{\mathrm{m}\mathrm{i}\mathrm{n}}J\left(k\right)={\displaystyle \sum _{i=\mathrm{1}}^{N_p}}\parallel \eta \left(k+i|t\right)-{\eta }_{\mathrm{r}\mathrm{e}\mathrm{f}}\left(k+i|t\right){\parallel }_{\mathit{Q}}^{\mathrm{2}}+{\displaystyle \sum _{i=\mathrm{1}}^{N_c}}\parallel \mathrm{\Delta }U\left(k+i|t\right){\parallel }_{\mathit{R}}^{\mathrm{2}}+\rho {\epsilon }^{\mathrm{2}}$$

（21）

式中： $\rho$和$\epsilon$分别是权重因子和松弛因子；目标函数第一项是使系统状态跟随期望车辆状态，参考变量包含车辆的参考横向和纵向位置，航向角，参考横摆角速度和质心侧偏角；第二项是对控制输入增量的约束，使车辆控制输入更加平顺；$\mathit{Q}$和$\mathit{R}$分别代表目标函数状态输入和控制输出的权重矩阵。

2.2.4　约束条件及优化求解

控制输入约束：

$$\begin{array}{c}{\mathit{u}}_{\mathrm{m}\mathrm{i}\mathrm{n}}\left(t+k\right)\le \mathit{u}\left(t+k\right)\le {\mathit{u}}_{\mathrm{m}\mathrm{a}\mathrm{x}}\left(t+k\right),\\ k=\mathrm{0,1},\cdots ,N_c-\mathrm{1}\end{array}$$

（22）

控制输入增量约束：

$$\begin{array}{c}\mathrm{\Delta }{\mathit{u}}_{\mathrm{m}\mathrm{i}\mathrm{n}}\left(t+k\right)\le \mathrm{\Delta }\mathit{u}\left(t+k\right)\le \mathrm{\Delta }{\mathit{u}}_{\mathrm{m}\mathrm{a}\mathrm{x}}\left(t+k\right),\\ k=\mathrm{0,1},\cdots ,N_c-\mathrm{1}\end{array}$$

（23）

系统输出约束：

$$\begin{array}{c}{\mathit{y}}_{\mathrm{m}\mathrm{i}\mathrm{n}}\left(t+k\right)\le \mathit{y}\left(t+k\right)\le {\mathit{y}}_{\mathrm{m}\mathrm{a}\mathrm{x}}\left(t+k\right),\\ k=\mathrm{0,1},\cdots ,N_c-\mathrm{1}\end{array}$$

（24）

目标函数求解的变量为控制时域内的控制增量，约束条件需要以控制增量或控制增量乘转换矩阵的形式表达。因此，对式（22）进行转换，求出相应的转换矩阵^［17］。设${\mathit{U}}_{\mathrm{t}}={\mathrm{1}}_{N_c}\otimes \mathit{u}\left(k-\mathrm{1}\right)$，

$$\mathit{A}=\underset{N_c\times N_c}{\underbrace{\left[\begin{array}{ccccc}\mathrm{1}& \mathrm{0}& \cdots & \cdots & \mathrm{0}\\ \mathrm{1}& \mathrm{1}& \mathrm{0}& \cdots & \mathrm{0}\\ \mathrm{1}& \mathrm{1}& \mathrm{1}& \ddots & \mathrm{0}\\ ⋮& ⋮& \ddots & \ddots & \mathrm{0}\\ \mathrm{1}& \mathrm{1}& \cdots & \mathrm{1}& \mathrm{1}\end{array}\right]}}\otimes {\mathit{I}}_m$$

（25）

式中： ${\mathrm{1}}_{{\mathrm{N}}_c}$表示行数为$N_c$的列向量； ${\mathit{I}}_m$表示维度为m的单位矩阵； $\otimes$表示克罗内克积；$\mathit{u}\left(k-\mathrm{1}\right)$表示上一时刻实际控制输入。结合式（22、23、25），可得^［18］：

$${\mathit{U}}_{\mathrm{m}\mathrm{i}\mathrm{n}}\le \mathit{A}\mathrm{\Delta }{\mathit{U}}_t+{\mathit{U}}_t\le {\mathit{U}}_{\mathrm{m}\mathrm{a}\mathrm{x}}$$

（26）

式中，$U_{\mathrm{m}\mathrm{i}\mathrm{n}}$和$U_{\mathrm{m}\mathrm{a}\mathrm{x}}$分别为控制时域内的控制量最小值和最大值的集合。

将预测方程（20）代入目标函数（21），经过矩阵计算，设预测时域内的跟踪误差${\mathit{e}}_t=({\mathit{\Theta }}_m\mathit{\xi }(k)+{\mathit{\Psi }}_m{\tilde{\mathit{d}}}_k-{\mathit{\eta }}_{\mathrm{r}\mathrm{e}\mathrm{f}})$，模型预测控制在每一步求解带约束的优化问题等价于求解以下标准二次规划问题^［12-14，18］：

$$J(\mathit{\xi }(t),\mathit{u}(t-\mathrm{1}),\mathrm{\Delta }\mathit{U}(t))=[\mathrm{\Delta }\mathit{U}{(t)}^{\mathrm{T}}{,\epsilon ]}^{\mathrm{T}}{\mathit{H}}_t[\mathrm{\Delta }\mathit{U}{(t)}^{\mathrm{T}},\epsilon ]+{\mathit{G}}_t[\mathrm{\Delta }\mathit{U}{(t)}^{\mathrm{T}},\epsilon ]$$

（27）

s.t. $\mathrm{\Delta }{\mathit{U}}_{\mathrm{m}\mathrm{i}\mathrm{n}}\le \mathrm{\Delta }{\mathit{U}}_t\le \mathrm{\Delta }{\mathit{U}}_{\mathrm{m}\mathrm{a}\mathrm{x}}$

$${\mathit{U}}_{\mathrm{m}\mathrm{i}\mathrm{n}}\le \mathit{A}\mathrm{\Delta }{\mathit{U}}_t+{\mathit{U}}_t\le {\mathit{U}}_{\mathrm{m}\mathrm{a}\mathrm{x}}$$

$${\mathit{Y}}_{\mathrm{m}\mathrm{i}\mathrm{n}}\le {\mathit{\Theta }}_m\mathit{\xi }\left(k\right)+{\mathit{\Gamma }}_m\mathrm{\Delta }{\mathit{U}}_m+{\mathit{\Psi }}_m{\mathit{D}}_k\le {\mathit{Y}}_{\mathrm{m}\mathrm{a}\mathrm{x}}$$

式中：${\mathit{H}}_t=\left[\begin{array}{cc}{\mathit{\Gamma }}_m^{\mathrm{T}}\mathit{Q}{\mathit{\Gamma }}_m+\mathit{R}& \mathrm{0}\\ \mathrm{0}& \rho \end{array}\right]$； ${\mathit{G}}_t=\left[\begin{array}{cc}\mathrm{2}{\mathit{e}}_t^{\mathrm{T}}\mathit{Q}{\mathit{\Gamma }}_m& \mathrm{0}\end{array}\right]$

通过在每个控制周期内求解上述带约束的目标函数，得到一系列控制时域的控制输入增量：

$$\mathrm{\Delta }{\mathit{U}}_t^{*}={\left[\begin{array}{c}\mathrm{\Delta }{\mathit{u}}_t^{*},\mathrm{\Delta }{\mathit{u}}_{t+\mathrm{1}}^{*},\dots ,\mathrm{\Delta }{\mathit{u}}_{k+N_c-\mathrm{1}}^{*}\end{array}\right]}^{\mathrm{T}}$$

（28）

将控制序列的第一个元素作为实际控制输入，作用于被控对象：

$$\mathit{u}(t)=\mathit{u}(t-\mathrm{1})+\mathrm{\Delta }{\mathit{u}}_t^{*}$$

（29）

2.3　基于PID的纵向速度跟踪控制

为了使被控车辆跟踪目标车速且减少计算量，本文采用PID控制生成总的驱动或制动转矩实现速度跟踪。

速度跟踪的误差定义为：

$$\begin{array}{c}{e}_u=u_d-u\\ e_{\dot{u}}={\dot{u}}_d-\dot{u}\\ e_{s_x}={\int }_{t_{\mathrm{1}}}^{t_{\mathrm{2}}}{e}_{u}dt=e_u{T}_s\end{array}$$

（30）

式中： $u_d$表示期望速度； $u$表示被控车辆模型的速度；${\dot{u}}_d$表示期望的纵向加速度； $\dot{u}$表示被控车辆模型的纵向加速度； $T_s$表示采样时间。

在PID控制中，当被控对象的车速低于期望车速时，总的驱动制动转矩分别为$T_a=K_P{e}_u+K_I{e}_{s_x}+K_D{e}_{\dot{u}}$， $T_b=\mathrm{0}$；当被控对象的车速高于期望车速时， $T_a=\mathrm{0}$，$T_b=K_P{e}_u+K_I{e}_{s_x}+K_D{e}_{\dot{u}}$；当被控对象的车速与期望车速相等时，$T_a=\mathrm{0}$， $T_b=\mathrm{0}$。

3 下层控制分配方法

下层转矩分配控制根据车辆稳定行驶的横摆力矩及总转矩需求，考虑电机、路面等约束条件，对分布式驱动四轮转矩进行动态分配，优化目标包括最小化车轮转矩跟踪误差，最小化各个轮胎利用率以及最小化纵向滑移能量损失，此外，还考虑轮胎路面附着极限和各电机输出转矩的约束。

（1） 转矩跟踪误差、期望横摆力矩、总转矩和下层各车轮转矩的动力学关系为如下^［19］：

$$\left\{\begin{array}{c}{T}_x=T_{\mathrm{f}\mathrm{l}}+T_{\mathrm{f}\mathrm{r}}+T_{\mathrm{r}\mathrm{l}}+T_{\mathrm{r}\mathrm{r}}\\ ∆M=\frac{\mathrm{1}}{R}\left(\left(a\mathrm{s}\mathrm{i}\mathrm{n}\delta -\frac{c_{\mathrm{f}}\mathrm{c}\mathrm{o}\mathrm{s}\delta }{\mathrm{2}}\right)·T_{\mathrm{f}\mathrm{l}}+\left(a\mathrm{s}\mathrm{i}\mathrm{n}\delta +\frac{c_{\mathrm{f}}\mathrm{c}\mathrm{o}\mathrm{s}\delta }{\mathrm{2}}\right)·T_{\mathrm{f}\mathrm{r}}-\frac{c_{\mathrm{r}}}{\mathrm{2}}·T_{\mathrm{r}\mathrm{l}}+\frac{c_{\mathrm{r}}}{\mathrm{2}}·T_{\mathrm{r}\mathrm{r}}\right)\end{array}\right.$$

（31）

方程组（31）可用式${\mathit{B}}_d{\mathit{u}}_d={\mathit{v}}_d$表示。其中： ${\mathit{v}}_d=[T_x{∆M]}^{\mathrm{T}};{\mathit{u}}_d=[T_{\mathrm{f}\mathrm{l}}{T}_{\mathrm{f}\mathrm{r}}{T}_{\mathrm{r}\mathrm{l}}{T}_{\mathrm{r}\mathrm{r}}{]}^{\mathrm{T}};{\mathit{B}}_d=\left[\begin{array}{cccc}\mathrm{1}& \mathrm{1}& \mathrm{1}& \mathrm{1}\\ \left(-\frac{c_{\mathrm{f}}\mathrm{c}\mathrm{o}\mathrm{s}\delta }{\mathrm{2}R}+\frac{a\mathrm{s}\mathrm{i}\mathrm{n}\delta }{R}\right)& \left(\frac{c_{\mathrm{f}}\mathrm{c}\mathrm{o}\mathrm{s}\delta }{\mathrm{2}R}+\frac{a\mathrm{s}\mathrm{i}\mathrm{n}\delta }{R}\right)& -\frac{c_{\mathrm{r}}}{\mathrm{2}R}& \frac{c_{\mathrm{r}}}{\mathrm{2}R}\end{array}\right]$，$T_{\mathrm{f}\mathrm{l}}\mathrm{、}{T}_{\mathrm{f}\mathrm{r}}\mathrm{、}{T}_{\mathrm{r}\mathrm{l}}\mathrm{、}{T}_{\mathrm{r}\mathrm{r}}$为各车轮转矩，$d$为轮距，$R$为各轮的有效滚动半径。这样，车轮转矩误差目标函数可表示为^［19］：

$${\mathit{J}}_{\mathrm{1}}=\frac{\mathrm{1}}{\mathrm{2}}{\left({\mathit{B}}_d{\mathit{u}}_d-{\mathit{v}}_d\right)}^{\mathrm{T}}{\mathit{W}}_{\mathrm{1}}\left({\mathit{B}}_d{\mathit{u}}_d-{\mathit{v}}_d\right)$$

（32）

（2）　轮胎利用率，该指标表达了轮胎稳定性裕量，轮胎利用率越高，则其稳定性裕量越小^［20］，目标函数如下：

$$J_{\mathrm{2}}={\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{4}}}\frac{F_{xi}^{\mathrm{2}}+F_{yi}^{\mathrm{2}}}{{\left(\mu F_{zi}\right)}^{\mathrm{2}}},i=\mathrm{1,2},\mathrm{3,4}$$

（33）

下层控制器尽量使各轮胎纵向力最小化，以保留侧向力裕度，提高车辆稳定性。为方便计算，忽略轮胎侧向力因素，只考虑轮胎纵向力，则轮胎利用率目标函数可简化为^［20］：

$$J_{\mathrm{2}}={\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{4}}}\frac{F_{xi}^{\mathrm{2}}}{{\left(\mu F_{zi}\right)}^{\mathrm{2}}}={{\mathit{u}}_d}^{\mathrm{T}}{W}_{\mathrm{2}}{\mathit{u}}_d$$

（34）

式中：$F_{xi}=\frac{T_i}{R};W_{\mathrm{2}}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}{\left(\frac{\mathrm{1}}{{\left({\mu }_i{F}_{zi}R\right)}^{\mathrm{2}}}\right)}_{\mathrm{4}\times \mathrm{4}}$。

（3） 纵向滑移能量损失，定义纵向滑移能量损失向量为：

${\mathit{P}}_{\mathrm{l}\mathrm{o}\mathrm{s}\mathrm{s}}={\left[F_{x\mathrm{1}}{\kappa }_{\mathrm{1}}{v}_{x\mathrm{1}}{F}_{x\mathrm{2}}{\kappa }_{\mathrm{2}}{v}_{x\mathrm{2}}{F}_{x\mathrm{3}}{\kappa }_{\mathrm{3}}{v}_{x\mathrm{3}}{F}_{x\mathrm{4}}{\kappa }_{\mathrm{4}}{v}_{x\mathrm{4}}\right]}^{\mathrm{T}}$，则其目标函数为^［21］：

$$J_{\mathrm{3}}={\mathit{P}}_{\mathrm{l}\mathrm{o}\mathrm{s}\mathrm{s}}{\left(\tau \right)}^{\mathrm{T}}{W}_{\tau }{\mathit{P}}_{\mathrm{l}\mathrm{o}\mathrm{s}\mathrm{s}}(\tau )={\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{4}}}{q}_i{\left(v_{xi}\left(\tau \right){\kappa }_i(\tau )\right)}^{\mathrm{2}}{\left(\frac{T_i(\tau )-I_s\dot{\omega }(\tau )}{r_{\omega }}\right)}^{\mathrm{2}}$$

（35）

式中： $W_{\tau }=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(q_{\mathrm{1}}{q}_{\mathrm{2}}{q}_{\mathrm{3}}{q}_{\mathrm{4}})$。因$\dot{\omega }(\tau )$在大多数情况下较小，$q_i{\left(v_{xi}\left(\tau \right){\kappa }_i(\tau )\right)}^{\mathrm{2}}{\left(\frac{T_i(\tau )-I_s\dot{\omega }(\tau )}{r_{\omega }}\right)}^{\mathrm{2}}$可近似表示为$q_i{\left(v_{xi}\left(\tau \right){\kappa }_i(\tau )\right)}^{\mathrm{2}}{T}_i{(\tau )}^{\mathrm{2}}$。当高速行驶出现较大滑移率时，时变的权重项$q_i{\left(v_{xi}\left(\tau \right){\kappa }_i(\tau )\right)}^{\mathrm{2}}$可有效约束转矩$T_i(\tau )$的增大，从而降低滑移率${\kappa }_i\left(\tau \right)$，即通过优化纵向滑移能量损失，各车轮的滑移率均能被自适应约束^［21］。

综上，下层转矩最优分配的综合目标函数为：

$$\mathrm{m}\mathrm{i}\mathrm{n}{J}_{\mathrm{C}\mathrm{A}}={\xi }_{\mathrm{1}}·J_{\mathrm{1}}+{\xi }_{\mathrm{2}}·J_{\mathrm{2}}+{\xi }_{\mathrm{3}}·J_{\mathrm{3}}$$

（36）

式中：第一项表示各轮输出转矩与期望转矩的误差；第二项表征轮胎利用率；第三项表征纵向滑移能量损失；${\xi }_{\mathrm{1}}$、${\xi }_{\mathrm{2}}$和${\xi }_{\mathrm{3}}$是各目标函数的权重因子。

将目标函数（36）转化为标准二次规划问题，同时考虑实际工程应用中电机所能提供的最大峰值转矩和路面附着条件约束：

$$\begin{array}{l}\mathrm{m}\mathrm{i}\mathrm{n}J=\frac{\mathrm{1}}{\mathrm{2}}{{\mathit{u}}_d}^{\mathrm{T}}H{\mathit{u}}_d+{\mathit{f}}^{\mathrm{T}}{\mathit{u}}_d\\ \mathrm{s}.\mathrm{t}.{\mathit{u}}_{ij\mathrm{m}\mathrm{i}\mathrm{n}}\le {\mathit{u}}_d\le {\mathit{u}}_{ij\mathrm{m}\mathrm{a}\mathrm{x}}\end{array}$$

（37）

式中： ${\mathit{u}}_{ij\mathrm{m}\mathrm{a}\mathrm{x}}=\mathrm{m}\mathrm{i}\mathrm{n}\left(T_{i\mathrm{m}\mathrm{a}\mathrm{x}},\mu R{F}_{zi}\right)$； ${\mathit{u}}_{ij\mathrm{m}\mathrm{i}\mathrm{n}}=\mathrm{m}\mathrm{a}\mathrm{x}\left(-T_{i\mathrm{m}\mathrm{a}\mathrm{x}},-\mu R{F}_{zi}\right)$； $T_{i\mathrm{m}\mathrm{a}\mathrm{x}}$为各电机输出峰值转矩； $H={{\mathit{B}}_d}^{\mathrm{T}}{\xi }_{\mathrm{1}}{W}_{\mathrm{1}}{\mathit{B}}_d+{\xi }_{\mathrm{2}}{W}_{\mathrm{2}}+{\xi }_{\mathrm{3}}{W}_{\mathrm{3}}$； $\mathit{f}={{\mathit{B}}_d}^{\mathrm{T}}{\xi }_{\mathrm{1}}{W}_{\mathrm{1}}{\mathit{v}}_d^{\mathrm{T}}$； $W_{\mathrm{3}}=W_{\tau }·\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}{\left({\left(v_{xi}\left(\tau \right){\kappa }_i(\tau )\right)}^{\mathrm{2}}\right)}_{\mathrm{4}\times \mathrm{4}}$。通过求解上述二次规划问题获得4个车轮最优分配转矩。

4 仿真分析

为验证所提控制算法在紧急避障场景下的有效性和鲁棒性，本文以双移线轨迹作为参考路径，该工况可充分验证高速、低附路面以及紧急避障等极限工况的控制性能。车辆建模参数如表1所示。

表1  车辆建模参数

Tab.1  Vehicle modeling parameters

参数	参数值	参数	参数值
$m_{\mathrm{s}}$/kg	1 400	$m_{\mathrm{u}\mathrm{l}\mathrm{f}}$/kg	80
$m$/kg	1 720	$m_{\mathrm{u}\mathrm{r}\mathrm{f}}$/kg	80
$g$/（m/${\mathrm{s}}^2$）	9.80	$m_{\mathrm{u}\mathrm{l}\mathrm{r}}$/kg	80
$I_x$/（$\mathrm{k}\mathrm{g}·{\mathrm{m}}^2$）	900	$m_{\mathrm{u}\mathrm{r}\mathrm{r}}$/kg	80
$I_y$/（$\mathrm{k}\mathrm{g}·{\mathrm{m}}^2$）	2 000	$k_{\mathrm{t}\mathrm{f}}$/（$\mathrm{N}/\mathrm{m}$）	200 000
$I_z$/（$\mathrm{k}\mathrm{g}·{\mathrm{m}}^2$）	2 420	$k_{\mathrm{t}\mathrm{r}}$/（$\mathrm{N}/\mathrm{m}$）	200 000
$a$/m	1.14	$C_{\alpha \mathrm{f}}$/（$\mathrm{N}/\mathrm{r}\mathrm{a}\mathrm{d}$）	-44 000
$b$/m	1.40	$C_{\alpha \mathrm{r}}$/（$\mathrm{N}/\mathrm{r}\mathrm{a}\mathrm{d}$）	-47 000
$h$/m	0.75	$C_{x\mathrm{f}}$/N	5 000
$c_{\mathrm{f}}$/m	1.50	$C_{x\mathrm{r}}$/N	5 000
$c_{\mathrm{r}}$/m	1.50	$r_0$/m	0.285
$k_{\mathrm{s}\mathrm{f}}$/（$\mathrm{N}/\mathrm{m}$）	35 000	$I_w$/（$\mathrm{k}\mathrm{g}·{\mathrm{m}}^2$）	1
$k_{\mathrm{s}\mathrm{r}}$/（$\mathrm{N}/\mathrm{m}$）	30 000	$h_{\mathrm{r}\mathrm{c}\mathrm{f}}$/m	0.65
$b_{\mathrm{s}\mathrm{f}}$/（$\mathrm{N}·\mathrm{s}/\mathrm{m}$）	2 500	$h_{\mathrm{r}\mathrm{c}\mathrm{r}}$/m	0.60
$b_{\mathrm{s}\mathrm{r}}$/（$\mathrm{N}·\mathrm{s}/\mathrm{m}$）	2 000	$\mu$	0.85/0.80/0.40

4.1　不同车速双移线工况

选择在附着条件良好的路面（$\mu =\mathrm{0.85}$）上，车辆的纵向速度分别为36 km/h、72 km/h和90 km/h，进行仿真试验，仿真结果如图4—图11所示。由图4可以看出，不同车速下，车辆可准确跟踪期望轨迹，随着车速的提高，车辆实际轨迹和参考轨迹跟踪的误差逐渐增大且存在一定滞后，但最大横向偏差小于0.28 m，在可接受范围内。由图5—图7看出，随着车速的提高，横摆角速度和质心侧偏角峰值逐渐增大，但都维持在极限范围内，前轮转角和附加横摆力矩的幅值随着车速的提高逐渐增大且有一定的提前。图8—图11看出左右两侧车轮转矩差使车辆获得附加横摆力矩。仿真结果表明，所提控制器使无人驾驶车辆在不同车速下均能较好地跟踪期望轨迹同时保持车辆横摆稳定性，具有较好的鲁棒性。

图4  轨迹跟踪比较

Fig.4  Comparison of tracking trajectory

图5  横摆角速度比较

Fig.5  Comparison of vehicle yaw rate

图6  质心侧偏角比较

Fig.6  Comparison of vehicle sideslip angle

图7  前轮转角的比较

Fig.7  Comparison of front wheel steering angle

图8  附加横摆力矩比较

Fig.8  Comparison of additional yaw moment

图9  转矩优化分配 （车速36km/h）

Fig.9  Optimal torque allocation, speed 36 km/h

图10  转矩优化分配 （车速72km/h）

Fig.10  Optimal torque allocation, speed 72 km/h

图11  转矩优化分配 （车速90km/h）

Fig.11  Optimal torque allocation, speed 90 km/h

4.2　不同路面附着条件双移线工况

为了验证所提控制算法在不同路面附着情况下的鲁棒性，以50 km/ h车速，分别在高附（$\mu =$0.8）和低附（$\mu =$0.3）路面进行仿真，仿真结果如图12—图19所示。由图12和图13可以看出，在不同路面附着系数下，车辆均可准确跟踪期望轨迹和期望车速，最大横向偏差＜0.12 m，最大纵向速度偏差＜0.062 m/s。在低附路面，车辆轨迹跟踪误差增大且存在一定滞后，跟踪效果变差，但仍满足精度要求。由图14—图19看出，在相同车速下，高附路面的横摆角度和质心侧偏角峰值小于低附路面，说明车辆在高附路面的稳定性更好，裕度更高。整个过程中，低附路面下横摆角速度和质心侧偏角的幅值都维持在极限范围内，车辆没有失稳，前轮转角和附加横摆力矩的幅值，高附路面相比低附路面增大且有一定的提前。以上结果表明，所提控制器使无人驾驶车辆在不同路面附着条件下均能较好地跟踪期望轨迹和车速并保持稳定性，具有较好的鲁棒性。

图12  轨迹跟踪比较

Fig.12  Comparison of tracking trajectory

图13  速度跟踪比较

Fig.13  Comparison of speed tracking

图14  质心侧偏角比较

Fig.14  Comparison of vehicle sideslip angle

图15  横摆角速度比较

Fig.15  Comparison of vehicle yaw rate

图16  前轮转角的比较

Fig.16  Comparison of front wheel steering angle

图17  附加横摆力矩比较

Fig.17  Comparison of additional yaw moment

图18  转矩优化分配 （μ=0.3）

Fig.18  Optimal torque allocation, μ=0.3

图19  转矩优化分配 （μ=0.8）

Fig.19  Optimal torque allocation, μ=0.8

4.3　有无稳定性控制的双移线工况

为了与仅考虑轨迹跟踪的情况进行比较，以20 m/s车速，路面附着系数0.4的双移线工况仿真，对比仅有无稳定性控制的轨迹跟踪效果。仿真结果如图20—图27所示。由图20和图21可以看出，考虑稳定性控制的轨迹跟踪效果优于单纯的轨迹跟踪，轨迹跟踪的横向偏差更小，纵向速度跟踪的误差也更小。由图22—图27可以看出，有稳定性控制的情况，其横摆角速度和质心侧偏角的均方根误差均小于不考虑稳定性控制的情况，且方向盘转角相比仅轨迹跟踪的情况，幅值更小，说明通过附加横摆力矩使得车辆在低附路面下，在保持轨迹跟踪精度的同时，提高了车辆的稳定性。

图20  轨迹跟踪比较

Fig.20  Comparison of tracking trajectory

图21  速度跟踪比较

Fig.21  Comparison of speed tracking

图22  质心侧偏角比较

Fig.22  Comparison of vehicle sideslip angle

图23  横摆角速度比较

Fig.23  Comparison of vehicle yaw rate

图24  前轮转角的比较

Fig.24  Comparison of front wheel steering angle

图25  航向角跟踪比较

Fig.25  Comparison of vehicle heading angle

图26  附加横摆力矩（考虑控制分配）

Fig.26  Comparison of additional yaw moment in case of considering optimal torque allocation

图27  车轮转矩优化分配结果

Fig.27  Four wheel optimal torque allocation

5 结论及展望

本文针对分布式驱动电动汽车紧急避障轨迹跟踪及稳定性控制问题，提出一种基于LTV MPC的分层架构的轨迹跟踪及稳定性协调控制方法，其上层为线性时变模型预测轨迹跟踪控制器，下层为基于二次规划的转矩优化分配控制器，通过不同车速，不同路面附着条件以及有无稳定性控制的双移线工况进行仿真验证，得出以下结论：

（1） 所设计的分层轨迹跟踪和稳定性控制策略在不同车速、不同附着条件下有较好的鲁棒性，在高速和低附情况下的跟踪误差在可接受范围内，满足轨迹跟踪精度和车辆横摆稳定性的要求。

（2） 通过有无考虑稳定性控制的轨迹跟踪仿真对比，进一步验证了控制算法在极限工况下的有效性，使被控车辆在低附路面保持轨迹跟踪精度的同时提高了横摆稳定性。

后续将对极限工况下考虑横摆侧倾耦合稳定性判据，融合侧倾稳定性轨迹跟踪方法，以及转矩分配目标函数权重系数自适应优化方法等进一步研究。