集成参数辨识的车辆动力学自适应模型预测控制

Adaptive MPC with integrated parameter identification for vehicle dynamics

摘要

针对载重和路面变化导致的车辆控制性能下降问题，提出了一种集成参数辨识的优化控制方法。首先，分析了模型参数变化对传统模型预测控制器的状态预测和约束的影响；之后，基于李雅普诺夫理论设计了自适应律，在线修正模型参数并补偿控制输入，使真实系统与名义模型响应一致；同时，提出一种基于Brush轮胎模型和递归最小二乘法的路面附着系数在线辨识方法，确保控制约束的准确性。在此基础上构建了分层控制器：上层为集成参数自适应的模型预测控制，下层为轮胎力分配器。通过CarSim/Simulink联合仿真验证，该方法能有效修正参数误差并保持控制性能，提升车辆的安全性、稳定性和鲁棒性。

Abstract

This paper addresses the issue of vehicle control performance degradation caused by variations in load and road conditions， and proposes an optimization control method integrated with parameter identification. First， the impact of model parameter changes on the state prediction and constraints of traditional Model Predictive Controllers（MPC） is analyzed. Then， an adaptive law based on Lyapunov theory is designed to online adjust model parameters and compensate for control inputs， ensuring that the real system's response aligns with the nominal model. Additionally， an online road adhesion coefficient identification method based on the Brush tire model and Recursive Least Squares algorithm is proposed to ensure the accuracy of control constraints. Building upon this， a hierarchical controller is constructed， with the upper layer being the MPC integrated with parameter adaptation， and the lower layer being the tire force allocator. Simulation results from CarSim/Simulink co-simulation demonstrate that this method effectively corrects parameter errors and maintains control performance， thereby improving vehicle safety， stability， and robustness.

Graphical abstract

关键词

控制理论与控制工程 / 车辆控制 / 分层控制 / 自适应控制 / 参数辨识 / 模型预测控制

Key words

control theory and control engineering / vehicle control / hierarchical control / adaptive control / parameter identification / model predictive control

0 引 言

车辆行驶过程中，由于载重的不同，车辆参数（如整车质量和转动惯量等）会随之发生改变，车辆行驶在不同路面行驶时附着系数也会改变，两种情况会影响车辆的动态响应和驱动能力^［1］。需要针对性地进行控制算法设计来应对车辆参数的变化带来的影响，这方面的研究值得关注与重视。

现有研究针对模型预测控制（MPC）在车辆控制中的模型参数不确定性问题提出了多种解决方案：Jiang等^［2］采用极限学习机预测误差并设计PID前馈补偿器；Mata等^［3］提出鲁棒MPC将预测偏差约束为有限集；Yamauchi等^［4］利用神经网络实时估计建模误差；Kim等^［5］通过基于模型的强化学习在线更新预测模型；Lee等^［6］在鲁棒MPC框架下集成转向与制动控制以减小滞后误差；Su等^［7］则采用高斯回归过程双补偿函数修正规划器与控制器模型参数。这些方法均通过不同途径提升了参数不确定性下的路径跟踪性能。车辆控制中，汽车的驱动能力与路面附着系数有关，即与MPC的约束条件有关，为准确描述其约束条件，相关研究值得注意，刘建泽等^［8］基于改进的最小二乘法并采集悬架响应数据进行辨识；李杰等^［9］基于车辆响应的NARX神经网络进行辨识，但他们的估计方法都依靠整车质量。

现有研究在提升控制性能方面虽取得进展，但仍存在局限^［10］。对于基于数据的方法，局限性在于对数据的强依赖性。这类方法通过历史数据建立输入输出映射关系，其性能受制于训练数据的完备性和代表性。当系统运行于数据分布之外的工况时，模型的泛化能力将下降。鲁棒控制方法能保证最坏情况的系统稳定性，但控制器牺牲了在常规工况下的最优性能；其次，当实际扰动超出预设的不确定性边界时，系统的性能退化往往呈现非线性突变的特征。参数自适应方法通过在线调整模型参数展现出更好的动态适应能力，但其有效性依赖于参数辨识的准确性。现有研究在解决参数耦合问题方面仍存在明显不足，特别是当多个关键参数需要同时辨识时，系统往往面临可辨识性降低的问题。这种参数间的相互干扰不仅影响收敛速度，可能导致辨识结果偏离真实值。目前，对于路面附着系数的识别方法大多需要依靠整车质量，而整车质量根据载重的不同会变化，影响辨识结果。

当被控对象的参数发生变化时，控制器会基于参数变化前的模型进行预测。然而，此时系统的实际状态可能与预测状态存在差异，从而影响控制效果，这给基于模型的控制器的应用带来巨大的影响。例如，在汽车在运输过程中，由于车内增加了货物或人员，车辆的转动惯量和质量发生变化，这些变化会影响车辆的状态响应，使汽车的实际运动状态与预测状态产生偏差，进而偏离行驶轨道，造成安全隐患。车辆运行过程中，地面与轮胎的附着系数发生变化，控制器无法准确描述驱动能力，导致汽车无法充分发挥其驱动能力或使得控制输入超出车辆的实际驱动能力，影响行驶安全。本文在MPC的基础上，提出了一种自适应律，以抑制模型参数变化引起的控制性能下降。首先，基于4WD4WS动力学模型，分析预测偏差的来源和其对状态响应造成的影响。随后，构造状态偏差方程，并应用李雅普诺夫稳定性理论分析模型参数收敛于真实值的条件，根据此条件确定自适应律调整控制器参数。另一方面，考虑到路面附着系数对轮胎力约束的影响，使用Brush轮胎模型与递归最小二乘法在线辨识路面附着系数，确保控制器的约束与实际驱动能力相符。本文设计分层控制器，其中改进的控制器作为上层控制器。下层控制器根据上层控制器的控制输入进行轮胎力的分配。主要贡献点在于，基于李雅普诺夫稳定性设计自适应律应用于控制器中提升控制精度，并基于Brush轮胎模型设计不依靠整车重量的附着系数辨识方法。最后，通过物理仿真实验，验证本文设计方法的有效性。

1 车辆动力学模型及其最大驱动能力

1.1　4WD4WS车辆动力学模型

本文所建立的车辆模型为考虑车辆纵向，横向和横摆运动，四轮转向和驱动的车辆动力学模型。为简化计算，假设车辆在水平路面上行驶，不考虑空气阻力和悬架的影响，如图1所示^［11］。

根据牛顿第二定律，车辆行驶状态与整车合力的状态方程为：

\left\{\begin{array}{l}{\dot{V}}_x=V_y{w}_{\mathrm{r}}+\frac{\sum F_x}{m}\\ {\dot{V}}_y=-V_x{w}_{\mathrm{r}}+\frac{\sum F_y}{m}\\ {\dot{w}}_{\mathrm{r}}=\frac{\sum M_z}{I_z}\end{array}\right.（1）

式中：V_x\mathrm{、}{V}_y\mathrm{、}{w}_{\mathrm{r}}为车辆的纵向速度、侧向速度以及横摆角速度；\sum F_x\mathrm{、}\sum F_y\mathrm{、}\sum M_z为纵向、侧向合力与横摆力矩；I_{\mathrm{z}}\mathrm{、}m为车辆的转动惯量和质量。

整车合力与和合力矩由轮胎提供，下式给出了轮胎力与整车合力与合力矩的关系：

\left\{\begin{array}{l}\sum F_x={\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{4}}}{F}_{xi}\mathrm{c}\mathrm{o}\mathrm{s}{\delta }_i-{\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{4}}}{F}_{yi}\mathrm{s}\mathrm{i}\mathrm{n}{\delta }_i\\ \sum F_y={\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{4}}}{F}_{xi}\mathrm{s}\mathrm{i}\mathrm{n}{\delta }_i+{\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{4}}}{F}_{yi}\mathrm{c}\mathrm{o}\mathrm{s}{\delta }_i\\ \sum M_z={\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{4}}}{\left(-\mathrm{1}\right)}^{i-\mathrm{1}}{l}_{\mathrm{s}}\left(F_{xi}\mathrm{c}\mathrm{o}\mathrm{s}{\delta }_i-F_{yi}\mathrm{s}\mathrm{i}\mathrm{n}{\delta }_i\right)+\\ {\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{2}}}{l}_{\mathrm{f}}\left(F_{xi}\mathrm{s}\mathrm{i}\mathrm{n}{\delta }_i+F_{yi}\mathrm{c}\mathrm{o}\mathrm{s}{\delta }_i\right){\displaystyle \sum _{i=\mathrm{3}}^{\mathrm{4}}}{l}_{\mathrm{r}}\left(F_{xi}\mathrm{s}\mathrm{i}\mathrm{n}{\delta }_i+F_{yi}\mathrm{c}\mathrm{o}\mathrm{s}{\delta }_i\right)\end{array}\right.（2）

式中：F_{xi}为轮胎纵向力；F_{yi}为轮胎侧向力；{\delta }_i为轮胎转向角；l_{\mathrm{s}}为半轮距；l_{\mathrm{f}}\mathrm{、}{l}_{\mathrm{r}}分别为车辆前轴与后轴到质心的距离。

控制系统通过控制驱动电机和轮胎侧偏角来控制轮胎力，式（3）给出了相应的轮胎逆模型和轮胎力矩平衡条件，从而确定下层控制器的控制输入：轮胎侧偏角和驱动电机扭矩。

（3）

式中：{\alpha }_i为轮胎侧偏角；T_{\mathrm{t}i}为驱动电机扭矩；C_{\mathrm{s}i}\mathrm{、}{C}_{\mathrm{\alpha }i}分别为轮胎纵向刚度与侧偏刚度；\mu为轮胎与地面的附着系数；F_{zi}为轮胎所受到的垂向力；w_i为轮胎旋转角速度；R_{\mathrm{e}i}为轮胎有效旋转半径；J_i为轮胎转动惯量。

1.2　最大驱动能力

整车驱动力由轮胎力提供，而由于摩擦圆的存在，即轮胎受到的摩擦力具有边界，所以整车驱动能力也具有边界^［12］：

\sum F_x\in \left[\begin{array}{cc}\underline{\varpi }& \overline{\varpi }\end{array}\right],\sum F_y\in \left[\underline{\tau },\overline{\tau }\right],\sum M_z\in \left[\begin{array}{cc}\underline{\kappa }& \overline{\kappa }\end{array}\right]（4）

其中：

\left\{\begin{array}{l}\overline{\varpi }={\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{4}}}{\overline{F}}_{mxi}\\ \underline{\varpi }={\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{4}}}{\underline{F}}_{mxi}\end{array}\right.，\left\{\begin{array}{l}\overline{\tau }={\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{4}}}{\overline{F}}_{myi}\\ \underline{\tau }={\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{4}}}{\underline{F}}_{myi}\end{array}\right.，\left\{\begin{array}{l}{\overline{\theta }}_i=\mathrm{s}\mathrm{i}{\mathrm{n}}^{-\mathrm{1}}\left(\raisebox{1ex}{$F_{yi}$}\!\left/ \!\raisebox{-1ex}{$\mu F_{zi}$}\right.\right)\\ {\underline{\theta }}_i=\mathrm{c}\mathrm{o}{\mathrm{s}}^{-\mathrm{1}}\left(\raisebox{1ex}{$F_{xi}$}\!\left/ \!\raisebox{-1ex}{$\mu F_{zi}$}\right.\right)\end{array}\right.；

\left\{\begin{array}{l}\overline{\kappa }={\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{4}}}{\left(-\mathrm{1}\right)}^{i-\mathrm{1}}{\overline{F}}_{mxi}{l}_{\mathrm{s}}+{\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{2}}}{\overline{F}}_{myi}{l}_{\mathrm{f}}-{\displaystyle \sum _{i=\mathrm{3}}^{\mathrm{4}}}{\overline{F}}_{myi}{l}_{\mathrm{r}}\\ \underline{\kappa }={\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{4}}}{\left(-\mathrm{1}\right)}^{i-\mathrm{1}}{\underline{F}}_{mxi}{l}_{\mathrm{s}}+{\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{2}}}{\underline{F}}_{myi}{l}_{\mathrm{f}}-{\displaystyle \sum _{i=\mathrm{3}}^{\mathrm{4}}}{\underline{F}}_{myi}{l}_{\mathrm{r}}\end{array}\right.，

\overline{\varpi }\mathrm{、}\underline{\varpi }\mathrm{、}\overline{\tau }\mathrm{、}\underline{\tau }\mathrm{、}\overline{\kappa }\mathrm{、}\underline{\kappa }分别表示整车所受纵向合力、侧向合力与横摆力矩的最值；{\overline{F}}_{mxi}\mathrm{、}{\underline{F}}_{mxi}\mathrm{、}{\overline{F}}_{myi}\mathrm{、}{\underline{F}}_{mxi}为单个轮胎车身坐标轴上纵向力与横向力的最值；{\overline{\theta }}_i\mathrm{、}{\underline{\theta }}_i为轮胎力导向角范围，表示在当前轮胎受力状态下，通过增加横向或纵向轮胎力至摩擦圆边界所能达到的轮胎力夹角范围。式（4）为模型预测控制器求解优化问题的约束范围，避免控制输入超出车辆驱动能力。

单个轮胎的轮胎力限制可表示为：

\left\{\begin{array}{l}{\overline{F}}_{mxi}=\mu F_{\mathrm{z}i}\mathrm{c}\mathrm{o}\mathrm{s}{\overline{\theta }}_i\mathrm{c}\mathrm{o}\mathrm{s}{\delta }_i-\mu F_{zi}\mathrm{s}\mathrm{i}\mathrm{n}{\overline{\theta }}_i\mathrm{s}\mathrm{i}\mathrm{n}{\delta }_i\\ {\overline{F}}_{myi}=\mu F_{zi}\mathrm{c}\mathrm{o}\mathrm{s}{\overline{\theta }}_i\mathrm{s}\mathrm{i}\mathrm{n}{\delta }_i+\mu F_{zi}\mathrm{s}\mathrm{i}\mathrm{n}{\overline{\theta }}_i\mathrm{c}\mathrm{o}\mathrm{s}{\delta }_i\\ {\underline{F}}_{myi}=\mu F_{zi}\mathrm{c}\mathrm{o}\mathrm{s}{\underline{\theta }}_i\mathrm{c}\mathrm{o}\mathrm{s}{\delta }_i-\mu F_{zi}\mathrm{s}\mathrm{i}\mathrm{n}{\underline{\theta }}_i\mathrm{s}\mathrm{i}\mathrm{n}{\delta }_i\\ {\underline{F}}_{mxi}=\mu F_{zi}\mathrm{c}\mathrm{o}\mathrm{s}{\underline{\theta }}_i\mathrm{s}\mathrm{i}\mathrm{n}{\delta }_i+\mu F_{zi}\mathrm{s}\mathrm{i}\mathrm{n}{\underline{\theta }}_i\mathrm{c}\mathrm{o}\mathrm{s}{\delta }_i\end{array}\right.（5）

式（5）为下层分配控制器对整车合力进行分配时对各轮胎上的轮胎力的约束。

2 控制系统结构描述

本文通过控制车辆的速度与横摆角速度，达到控制汽车按照规划路线行驶的目标。将模型预测控制器作为上层控制器，控制整车驱动力，下层控制器根据整车所需驱动力情况将上层控制器的控制输入分配给各个轮胎，并尽可能使轮胎力远离摩擦圆边界。同时，设计自适应律在线修正模型预测控制器参数。控制系统结构如图2所示。

3 控制器设计

为了使汽车的行驶状态能准确跟随期望值，并尽可能减小控制输入的波动，且确保控制输入始终保持在车辆驱动能力范围内，本节基于名义模型式（1）设计了具有控制输入约束的模型预测控制器作为上层控制器。通过对系统状态进行预测并进行滚动优化，使得控制输入能使汽车快速且平稳地实现状态跟踪^{［13，14］}。

式（1）在工作点处进行离散线性化，可得：

\mathit{x}(k+\mathrm{1})={\mathit{A}}_k\mathit{x}(k)+{\mathit{B}}_k\tilde{\mathit{u}}(k)（6）

式中：\mathit{x}(k)为当前时刻系统状态与工作点的偏差；\tilde{\mathit{u}}(k)为相对于工作点稳态值的偏差控制输入，控制器根据状态偏差输入\tilde{\mathit{u}}(k)来调整系统状态。

为了对控制增量施加约束，将k时刻控制增量与k-\mathrm{1}时刻的控制输入拆分，构造新的状态向量\mathit{\xi }(k)={\left[\mathit{x}\left(k\right),\tilde{\mathit{u}}\left(k-\mathrm{1}\right)\right]}^{\mathrm{T}}。\mathrm{\Delta }\tilde{\mathit{u}}(k)=\tilde{\mathit{u}}\left(k\right)-\tilde{\mathit{u}}\left(k-\mathrm{1}\right)为此时刻相对于上一时刻的控制增量。根据式（6）得到新状态空间表达式：

\left\{\begin{array}{l}\mathit{\xi }(k+\mathrm{1})={\mathit{A}}_{\xi }\mathit{\xi }(k)+{\mathit{B}}_{\xi }\mathrm{\Delta }\tilde{\mathit{u}}(k)\\ \mathit{\eta }(k)={\mathit{C}}_{\xi }\mathit{\xi }(k)\end{array}\right.（7）

定义预测、控制时域为N_{\mathrm{p}}和N_{\mathrm{c}}，k时刻的预测状态序列和待求输入增量序列为\mathit{Y}\left(k\right)和\mathrm{\Delta }\mathit{U}\left(k\right)，分别扩展到预测、控制时域。根据\mathit{\xi }(k)和\mathrm{\Delta }\mathit{U}\left(k\right)对系统输出进行预测的方程为：

\mathit{Y}(k)={\mathit{\psi }}_k\mathit{\xi }(k)+{\mathit{\theta }}_k\mathrm{\Delta }\mathit{U}(k)（8）

目标函数式（9）分为两部分，分别表示预测偏差和控制增量的大小：

\begin{array}{c}J(\mathit{\xi }(k),\tilde{\mathit{u}}(k-\mathrm{1}),\mathrm{\Delta }\mathit{U}(k))=\\ {\displaystyle \sum _{i=\mathrm{1}}^{N_{\mathrm{p}}}}||\mathit{\eta }(k+i|k)|{|}_{\mathit{Q}}^{\mathrm{2}}+{\displaystyle \sum _{i=\mathrm{0}}^{N_{\mathrm{c}}-\mathrm{1}}}||\mathrm{\Delta }\tilde{\mathit{u}}(k+i|k)|{|}_{\mathit{R}}^{\mathrm{2}}\end{array}（9）

式中：\mathit{Q}和\mathit{R}分别为预测偏差和控制增量的权重矩阵，为正定对角矩阵；\mathit{\eta }(k+i|k)为在k时刻时，对系统在k+i时刻输出\mathit{x}的预测；\mathrm{\Delta }\tilde{\mathit{u}}(k+i|k)为k+i时刻相对于上一时刻的控制输入增量。

在每个周期内求解优化问题：

\underset{\mathrm{\Delta }U}{\mathrm{m}\mathrm{i}\mathrm{n}}J\left(\mathit{\xi }(k),\tilde{\mathit{u}}(k-\mathrm{1}),\mathrm{\Delta }\mathit{U}(k)\right)（10）

\begin{array}{l}\mathrm{s}.\mathrm{t}.\mathrm{\Delta }{\mathit{U}}_{\mathrm{m}\mathrm{i}\mathrm{n}}\le \mathrm{\Delta }\mathit{U}\left(k\right)\le \mathrm{\Delta }{\mathit{U}}_{\mathrm{m}\mathrm{a}\mathrm{x}}\\ {\mathit{U}}_{\mathrm{m}\mathrm{i}\mathrm{n}}\le \mathit{I}\mathrm{\Delta }\mathit{U}\left(k\right)+\mathit{U}\left(k-\mathrm{1}\right)\le {\mathit{U}}_{\mathrm{m}\mathrm{a}\mathrm{x}}\end{array}

式中：\mathrm{\Delta }{\mathit{U}}_{\mathrm{m}\mathrm{i}\mathrm{n}}\mathrm{、}\mathrm{\Delta }{\mathit{U}}_{\mathrm{m}\mathrm{a}\mathrm{x}}为控制增量约束扩展到控制时域的增广形式，以防止控制增量过大。

由于系统线性化过程中引入了前馈控制，则{\mathit{U}}_{\mathrm{m}\mathrm{i}\mathrm{n}}\mathrm{、}{\mathit{U}}_{\mathrm{m}\mathrm{a}\mathrm{x}}为式（4）表示的约束减去前馈控制所得到控制器输入约束的控制时域扩展形式。\mathit{U}\left(k-\mathrm{1}\right)为\tilde{\mathit{u}}(k-\mathrm{1})扩展到控制时域的增广形式。\mathit{I}为一个N_{\mathrm{c}}维方阵，其中下三角部分的元素全为1。

求解结果为控制增量序列，将该序列的第1个元素作用于系统。下一个周期重复上述过程，实现被控车辆对参考信号的跟踪。本文的下层控制器^［12］将上层控制器控制输入转化为各轮胎转向角和驱动转矩。

首先，求解式（11）得到各轮胎应分配的轮胎力为：

\underset{u}{\mathrm{m}\mathrm{i}\mathrm{n}}J(\mathit{u})={\mathit{u}}^{\mathrm{T}}\left[\begin{array}{cc}\mathit{W}& \mathrm{0}\\ \mathrm{0}& \mathit{W}\end{array}\right]\mathit{u}\mathrm{s}.\mathrm{t}.(\mathrm{2}),(\mathrm{4})（11）

式中：\mathit{u}=\left[F_{x\mathrm{1}}\cdots F_{x\mathrm{4}},F_{y\mathrm{1}}\cdots F_{y\mathrm{4}}\right]表示应分配的轮胎力；\mathit{W}=diag\{[\raisebox{1ex}{$\mathrm{1}$}\!\left/ \!\raisebox{-1ex}{$(\mu F_{z\mathrm{1}}{)}^{\mathrm{2}}$}\right.\cdots \raisebox{1ex}{$\mathrm{1}$}\!\left/ \!\raisebox{-1ex}{$(\mu F_{z\mathrm{4}}{)}^{\mathrm{2}}$}\right.]\}。

之后根据式（3）确定各轮胎侧偏角和驱动力矩后控制各轮胎转向角和驱动电机。

4 模型不确定性分析

模型预测控制方法是一种基于系统模型设计的控制方法，预测结果会受到车辆动力学模型中质量和转动惯量等参数误差的影响，影响控制效果。接下来，将针对这一问题进行详细说明。

使用名义参数进行建模的预测方程为：

\mathit{x}(k+\mathrm{1})={\mathit{A}}_k\mathit{x}(k)+{\mathit{B}}_k\tilde{\mathit{u}}(k-\mathrm{1})+{\mathit{B}}_k\mathrm{\Delta }\tilde{\mathit{u}}(k)（12）

将式（12）中的名义值替换为真实值，可得：

\begin{array}{c}{\mathit{x}}_{\mathrm{r}\mathrm{e}\mathrm{a}}(k+\mathrm{1})={\mathit{A}}_k\mathit{x}(k)+\\ {\mathit{B}}_{k\mathrm{r}\mathrm{e}\mathrm{a}}\tilde{\mathit{u}}(k-\mathrm{1})+{\mathit{B}}_{k\mathrm{r}\mathrm{e}\mathrm{a}}\mathrm{\Delta }\tilde{\mathit{u}}(k)\end{array}（13）

式中：{\mathit{B}}_{k\mathrm{r}\mathrm{e}\mathrm{a}}=diag\left\{\left[\begin{array}{ccc}\raisebox{1ex}{$T$}\!\left/ \!\raisebox{-1ex}{$m_{\mathrm{r}\mathrm{e}\mathrm{a}}$}\right.& \raisebox{1ex}{$T$}\!\left/ \!\raisebox{-1ex}{$m_{\mathrm{r}\mathrm{e}\mathrm{a}}$}\right.& \raisebox{1ex}{$T$}\!\left/ \!\raisebox{-1ex}{$I_{z\mathrm{r}\mathrm{e}\mathrm{a}}$}\right.\end{array}\right]\right\}，T为采样周期；m_{\mathrm{r}\mathrm{e}\mathrm{a}}\mathrm{、}{I}_{z\mathrm{r}\mathrm{e}\mathrm{a}}分别为质量与转动惯量的真实值；{\mathit{x}}_{\mathrm{r}\mathrm{e}\mathrm{a}}(k+\mathrm{1})和\mathit{x}(k+\mathrm{1})分别表示在相同控制输入作用下，真实模型和名义模型在当前状态\mathit{x}(k)下的下一时刻状态。

假定一种情景，控制系统使用基于名义模型设计的模型预测控制器，则\mathrm{\Delta }\tilde{\mathit{u}}(k)为名义模型预测控制器的控制输入增量，{\mathit{x}}_{\mathrm{r}\mathrm{e}\mathrm{a}}(k+\mathrm{1})为下一时刻的基于真实参数的预测模型的预测状态，也为系统在下一时刻的实际状态。为了得到\mathrm{\Delta }\tilde{\mathit{u}}(k)与实际状态和名义预测模型的预测状态的关系，将式（12）变形为：

\begin{array}{c}\mathrm{\Delta }\tilde{\mathit{u}}(k)=-{\mathit{B}}_k^{-\mathrm{1}}{\mathit{A}}_{\mathrm{k}}\mathit{x}(k)-\\ \tilde{\mathit{u}}(k-\mathrm{1})+{\mathit{B}}_k^{-\mathrm{1}}\mathit{x}(k+\mathrm{1})\end{array}（14）

为了得到在此控制增量下系统下一时刻的状态量，将式（14）代入式（13），有：

\begin{array}{c}{\mathit{x}}_{\mathrm{r}\mathrm{e}\mathrm{a}}(k+\mathrm{1})={\mathit{A}}_k\mathit{x}(k)+{\mathit{B}}_{k\mathrm{r}\mathrm{e}\mathrm{a}}\tilde{\mathit{u}}(k-\mathrm{1})+\\ {\mathit{B}}_{k\mathrm{r}\mathrm{e}\mathrm{a}}[-{\mathit{B}}_k^{-\mathrm{1}}{\mathit{A}}_k\mathit{x}(k)-\tilde{\mathit{u}}(k-\mathrm{1})+{\mathit{B}}_k^{-\mathrm{1}}\mathit{x}(k+\mathrm{1})]=\\ (\mathit{E}-{\mathit{B}}_{k\mathrm{r}\mathrm{e}\mathrm{a}}{\mathit{B}}_k^{-\mathrm{1}}){\mathit{A}}_k\mathit{x}(k)+{\mathit{B}}_{k\mathrm{r}\mathrm{e}\mathrm{a}}{\mathit{B}}_k^{-\mathrm{1}}\mathit{x}(k+\mathrm{1})\end{array}（15）

式（15）减去式（12）得到参数误差对状态预测的影响为：

\begin{array}{c}{\mathit{x}}_{\mathrm{r}\mathrm{e}\mathrm{a}}(k+\mathrm{1})-\mathit{x}(k+\mathrm{1})=(\mathit{E}-{\mathit{B}}_{k\mathrm{r}\mathrm{e}\mathrm{a}}{\mathit{B}}_k^{-\mathrm{1}}){\mathit{A}}_k\mathit{x}(k)-\\ (\mathit{E}-{\mathit{B}}_{k\mathrm{r}\mathrm{e}\mathrm{a}}{\mathit{B}}_k^{-\mathrm{1}})\mathit{x}(k+\mathrm{1})\end{array}（16）

式（16）说明，在模型参数误差的影响下，传统模型预测控制方法的预测结果与实际状态存在偏差，虽然这种影响可通过滚动优化的方式在下一时刻有得到补偿趋势，但对于当前时刻会导致控制偏差扩大的结果，从而影响系统的控制性能。而当模型参数得到修正后，即{\mathit{B}}_{k\mathrm{r}\mathrm{e}\mathrm{a}}={\mathit{B}}_k时，根据式（16），名义预测模型中的下一时刻状态量和实际系统下一时刻的状态量的差值为零，即{\mathit{x}}_{\mathrm{r}\mathrm{e}\mathrm{a}}(k+\mathrm{1})-\mathit{x}(k+\mathrm{1})=\mathrm{0}时，预测状态与真实状态一致，控制效果得到改善。

5 自适应律设计

为解决由于参数误差导致的控制器预测状态与实际状态不符的问题，根据名义和假设的真实参数的车辆动力学方程构造状态偏差方程，之后应用李雅普诺夫稳定性理论^［15］分析模型参数与预测偏差收敛于真实值的条件来构建自适应律，从而在线调整控制器参数并对控制输入进行补偿。

设真实模型和名义模型的状态方程为：

{\mathit{x}}_{\mathrm{r}\mathrm{e}\mathrm{a}}(k+\mathrm{1})={\mathit{A}}_{\mathrm{k}}{\mathit{x}}_{\mathrm{r}\mathrm{e}\mathrm{a}}(k)+{\mathit{B}}_{\mathrm{k}\mathrm{r}\mathrm{e}\mathrm{a}}\tilde{\mathit{c}}(k)（17）

\mathit{x}(k+\mathrm{1})={\mathit{A}}_{\mathrm{k}}\mathit{x}(k)+{\mathit{B}}_{\mathrm{k}}\tilde{\mathit{u}}(k)（18）

定义名义模型与真实模型的状态偏差为：

\mathit{e}(k)=\mathit{x}(k)-{\mathit{x}}_{\mathrm{r}\mathrm{e}\mathrm{a}}(k)（19）

式中：{\mathit{x}}_{\mathrm{r}\mathrm{e}\mathrm{a}}为真实模型状态量；\mathit{x}为名义模型在上一时刻预测此时刻的状态量；\tilde{\mathit{c}}为补偿控制输入。

本文设计的核心思想是通过补偿控制输入，使上层控制器的输入作用于被控对象时，系统的状态响应与预测一致，设\tilde{\mathit{c}}与其他可观测参数的关系为：

\tilde{\mathit{c}}(k)=\mathit{\alpha }\tilde{\mathit{u}}(k)+\mathit{b}{\mathit{x}}_{\mathrm{r}\mathrm{e}\mathrm{a}}(k)+de(k)（20）

式中：\mathit{\alpha }\mathrm{、}\mathit{b}\mathrm{、}\mathit{d}为实对角矩阵，矩阵维数为3×3。

将式（17）~（20）联立，得到状态偏差\mathit{e}的状态方程为：

\begin{array}{c}\mathit{e}(k+\mathrm{1})=\mathit{x}(k+\mathrm{1})-{\mathit{x}}_{\mathrm{r}\mathrm{e}\mathrm{a}}(k+\mathrm{1})=\\ ({\mathit{A}}_k-{\mathit{B}}_{k\mathrm{r}\mathrm{e}\mathrm{a}}\mathit{d})\mathit{e}(k)-{\mathit{B}}_{k\mathrm{r}\mathrm{e}\mathrm{a}}\mathit{b}{\mathit{x}}_{\mathrm{r}\mathrm{e}\mathrm{a}}(k)+\\ ({\mathit{B}}_k-{\mathit{B}}_{k\mathrm{r}\mathrm{e}\mathrm{a}}\mathit{\alpha })\tilde{\mathit{u}}(k)\end{array}（21）

根据李雅普诺夫稳定性定理，当({\mathit{A}}_k-{\mathit{B}}_{k\mathrm{r}\mathrm{e}\mathrm{a}}\mathit{d})的特征值的绝对值都小于1，且-{\mathit{B}}_{k\mathrm{r}\mathrm{e}\mathrm{a}}\mathit{b}与({\mathit{B}}_k-{\mathit{B}}_{k\mathrm{r}\mathrm{e}\mathrm{a}}\mathit{\alpha })为元素全为0的矩阵时，\mathit{e}在原点处稳定。

为了使-{\mathit{B}}_{k\mathrm{r}\mathrm{e}\mathrm{a}}\mathit{b}与{\mathit{B}}_k-{\mathit{B}}_{k\mathrm{r}\mathrm{e}\mathrm{a}}\mathit{\alpha }能稳定于原点，对于第一项，要求\mathit{b}始终为0矩阵。对于第二项，设状态量{\mathit{E}}_{\mathrm{2}}={\mathit{B}}_k-{\mathit{B}}_{k\mathrm{r}\mathrm{e}\mathrm{a}}\mathit{\alpha }，通过对时间求导并离散，得到：

{\mathit{E}}_{\mathrm{2}}(k+\mathrm{1})={\mathit{E}}_{\mathrm{2}}(k)-{\mathit{B}}_{k\mathrm{r}\mathrm{e}\mathrm{a}}(\mathit{\alpha }(k+\mathrm{1})-\mathit{\alpha }(k))（22）

综上所述，为使\mathit{e}稳定于原点，需要使\mathit{e},{\mathit{E}}_{\mathrm{2}}稳定于原点，则构造李雅普诺夫函数为：

\begin{array}{c}V(\mathit{e}(k),{\mathit{e}}_{\mathrm{2}}(k))={\mathit{e}}_{\mathrm{2}}^{\mathrm{T}}(k){\mathit{\psi }}_{\mathrm{2}}{\mathit{B}}_{\mathrm{k}\mathrm{r}\mathrm{e}\mathrm{a}}^{-\mathrm{1}}{\mathit{e}}_{\mathrm{2}}^{}(k)+\\ {\mathit{e}}^{\mathrm{T}}(k){\mathit{\psi }}_{\mathrm{1}}({\mathit{A}}_k-{\mathit{B}}_{k\mathrm{r}\mathrm{e}\mathrm{a}}{\mathit{d})}^{-\mathrm{1}}\mathit{e}(k)\end{array}（23）

式中：{\mathit{\psi }}_{\mathrm{1}}\mathrm{、}{\mathit{\psi }}_{\mathrm{2}}为正实数对角矩阵，可记为{\mathit{\psi }}_{\mathrm{1}}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\{[{\psi }_{\mathrm{11}},{\psi }_{\mathrm{12}},{\psi }_{\mathrm{13}}]\}，{\mathit{\psi }}_{\mathrm{2}}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\{[{\psi }_{\mathrm{21}},{\psi }_{\mathrm{22}},{\psi }_{\mathrm{23}}]\}。

对矩阵{\mathit{E}}_{\mathrm{2}}和\mathit{\alpha }进行处理，得到由对角元素组成的列矩阵，以方便构造李雅普诺夫函数，则设{\mathit{e}}_{\mathrm{2}}={\mathit{E}}_{\mathrm{2}}{\mathit{I}}_{\mathrm{3}\times \mathrm{1}},\mathit{a}=\mathit{\alpha }{\mathit{I}}_{\mathrm{3}\times \mathrm{1}}，{\mathit{I}}_{\mathrm{3}\times \mathrm{1}}为元素全为1的三阶列矩阵，须选择合适的\mathit{d}使({\mathit{A}}_k-{\mathit{B}}_{k\mathrm{r}\mathrm{e}\mathrm{a}}{\mathit{d})}^{-\mathrm{1}}正定。

将\mathit{e}(k+\mathrm{1})\mathrm{、}{\mathit{e}}_{\mathrm{2}}(k+\mathrm{1})代入式（23）并与式（21）联立，得到k+\mathrm{1}时刻的李雅普诺夫函数与\mathit{e}(k)\mathrm{、}{\mathit{e}}_{\mathrm{2}}(k)的关系为：

\begin{array}{c}V(\mathit{e}(k+\mathrm{1}),{\mathit{e}}_{\mathrm{2}}(k+\mathrm{1}))={\mathit{e}}_{\mathrm{2}}^{\mathrm{T}}(k){\mathit{\psi }}_{\mathrm{2}}{\mathit{B}}_{k\mathrm{r}\mathrm{e}\mathrm{a}}^{-\mathrm{1}}{\mathit{e}}_{\mathrm{2}}^{}(k)+\\ {\mathit{e}}^{\mathrm{T}}(k)({\mathit{A}}_k-{\mathit{B}}_{k\mathrm{r}\mathrm{e}\mathrm{a}}{\mathit{d})}^{\mathrm{T}}{\mathit{\psi }}_{\mathrm{1}}\mathit{e}(k)+\\ \mathrm{2}({\mathit{\psi }}_{\mathrm{1}}{\mathit{I}}_{\mathrm{3}\times \mathrm{1}}\times \tilde{\mathit{u}}(k)\times \mathit{e}{(k))}^{\mathrm{T}}{\mathit{e}}_{\mathrm{2}}(k)-\\ \mathrm{2}({\mathit{\psi }}_{\mathrm{2}}{\mathit{I}}_{\mathrm{3}\times \mathrm{1}}\times (\mathit{a}(k+\mathrm{1})-\mathit{a}{(k)))}^{\mathrm{T}}{\mathit{e}}_{\mathrm{2}}(k)+\\ (\tilde{\mathit{u}}(k)\times {\mathit{e}}_{\mathrm{2}}{(k))}^{\mathrm{T}}{\mathit{\psi }}_{\mathrm{1}}({\mathit{A}}_k-{\mathit{B}}_{k\mathrm{r}\mathrm{e}\mathrm{a}}{\mathit{d})}^{-\mathrm{1}}\tilde{\mathit{u}}(k)\times {\mathit{e}}_{\mathrm{2}}(k)+\\ ({\mathit{B}}_{k\mathrm{r}\mathrm{e}\mathrm{a}}(\mathit{a}(k+\mathrm{1})-\mathit{a}{(k)))}^{\mathrm{T}}{\mathit{\psi }}_{\mathrm{2}}(\mathit{a}(k+\mathrm{1})-\mathit{a}(k))\end{array}（24）

将式（24）与式（23）相减得到相邻两时刻的离散李雅普诺夫函数的差值：

\begin{array}{c}\mathrm{\Delta }V=-{\mathit{e}}^{\mathrm{T}}(k){\mathit{\psi }}_{\mathrm{1}}({\mathit{A}}_k-{\mathit{B}}_{k\mathrm{r}\mathrm{e}\mathrm{a}}{\mathit{d})}^{-\mathrm{1}}\mathit{e}(k)+\\ (\tilde{\mathit{u}}(k)\times {\mathit{e}}_{\mathrm{2}}{(k))}^{\mathrm{T}}{\mathit{\psi }}_{\mathrm{1}}({\mathit{A}}_k-{\mathit{B}}_{k\mathrm{r}\mathrm{e}\mathrm{a}}{\mathit{d})}^{-\mathrm{1}}\tilde{\mathit{u}}(k)\times {\mathit{e}}_{\mathrm{2}}(k)+\\ {\mathit{e}}^{\mathrm{T}}(k)({\mathit{A}}_k-{\mathit{B}}_{k\mathrm{r}\mathrm{e}\mathrm{a}}{\mathit{d})}^{\mathrm{T}}{\mathit{\psi }}_{\mathrm{1}}\mathit{e}(k)+\\ \mathrm{2}({\mathit{\psi }}_{\mathrm{1}}{\mathit{I}}_{\mathrm{3}\times \mathrm{1}}\times \tilde{\mathit{u}}(k)\times \mathit{e}{(k))}^{\mathrm{T}}{\mathit{e}}_{\mathrm{2}}(k)-\\ \mathrm{2}({\mathit{\psi }}_{\mathrm{2}}{\mathit{I}}_{\mathrm{3}\times \mathrm{1}}\times (\mathit{a}(k+\mathrm{1})-\mathit{a}{(k)))}^{\mathrm{T}}{\mathit{e}}_{\mathrm{2}}(k)+\\ ({\mathit{B}}_{k\mathrm{r}\mathrm{e}\mathrm{a}}(\mathit{a}(k+\mathrm{1})-\mathit{a}{(k)))}^{\mathrm{T}}{\mathit{\psi }}_{\mathrm{2}}(\mathit{a}(k+\mathrm{1})-\mathit{a}(k))\end{array}（25）

为使\mathit{e}\mathrm{、}{\mathit{e}}_{\mathrm{2}}稳定于原点，需要寻求\mathrm{\Delta }V恒不大于0的条件，接下来对此进行说明。

首先，消除式（25）中的奇数项，得到：

\mathit{a}(k+\mathrm{1})=\mathit{a}(k)+{\mathit{\psi }}_{\mathrm{1}}{\mathit{\psi }}_{\mathrm{2}}^{-\mathrm{1}}{\mathit{I}}_{\mathrm{3}\times \mathrm{1}}\times \tilde{\mathit{u}}(k)\times \mathit{e}(k)（26）

假设由于\mathit{a}的离散积分作用，完成了参数整定，此时{\mathit{e}}_{\mathrm{2}}在原点附近，则将{\mathit{e}}_{\mathrm{2}}替换为关于式（26）的增量形式，将式（22）与（26）联立得到式（27）：

{\mathit{e}}_{\mathrm{2}}(k)={\mathit{B}}_{k\mathrm{r}\mathrm{e}\mathrm{a}}{\mathit{\psi }}_{\mathrm{1}}{\mathit{\psi }}_{\mathrm{2}}^{-\mathrm{1}}{\mathit{I}}_{\mathrm{3}\times \mathrm{1}}\times \tilde{\mathit{u}}(k)\times \mathit{e}(k)（27）

将式（26）（27）代入式（25）得到式（28），讨论\mathit{d}满足系统稳定所所需要的条件。

\begin{array}{c}\mathrm{\Delta }V=(\tilde{\mathit{u}}(k)\times \mathit{e}{(k))}^{\mathrm{T}}{\mathit{M}}_{\mathrm{1}}\tilde{\mathit{u}}(k)\times \mathit{e}(k)+(\tilde{\mathit{u}}(k)\times \\ \tilde{\mathit{u}}(k)\times \mathit{e}{(k))}^{\mathrm{T}}{\mathit{M}}_{\mathrm{2}}(\tilde{\mathit{u}}(k)\times \tilde{\mathit{u}}(k)\times \mathit{e}(k))+{\mathit{e}}^{\mathrm{T}}(k)·\\ (({\mathit{A}}_k-{\mathit{B}}_{k\mathrm{r}\mathrm{e}\mathrm{a}}{\mathit{d})}^{\mathrm{T}}{\mathit{\psi }}_{\mathrm{1}}-{\mathit{\psi }}_{\mathrm{1}}({\mathit{A}}_k-{\mathit{B}}_{k\mathrm{r}\mathrm{e}\mathrm{a}}{\mathit{d})}^{-\mathrm{1}})\mathit{e}(k)\end{array}（28）

其中，{\mathit{M}}_{\mathrm{1}}={\mathit{\psi }}_{\mathrm{1}}^{\mathrm{2}}({\mathit{A}}_k-{\mathit{B}}_{k\mathrm{r}\mathrm{e}\mathrm{a}}{\mathit{d})}^{-\mathrm{2}}{\mathit{\psi }}_{\mathrm{2}}^{-\mathrm{1}}{\mathit{B}}_{k\mathrm{r}\mathrm{e}\mathrm{a}}；{\mathit{M}}_{\mathrm{2}}={\mathit{\psi }}_{\mathrm{1}}^{\mathrm{3}}({\mathit{A}}_k-{\mathit{B}}_{k\mathrm{r}\mathrm{e}\mathrm{a}}{\mathit{d})}^{-\mathrm{3}}{\mathit{\psi }}_{\mathrm{2}}^{-\mathrm{2}}{\mathit{B}}_{k\mathrm{r}\mathrm{e}\mathrm{a}}^{\mathrm{2}}。

式（28）中，\mathit{d}要使关于\mathit{e}的二次型矩阵非正定，\tilde{\mathit{u}}在完成了参数整定后其中元素的值很小，选择具有较小元素的{\mathit{\psi }}_{\mathrm{1}}和较大元素的{\mathit{\psi }}_{\mathrm{2}}有利于稳定。根据式（21）稳定条件和参数范围可找到\mathit{d}，并能存在一定裕量，使存在\tilde{\mathit{u}}的项不会导致式（28）正定，从而说明参数辨识稳定在真值上。以上说明了以式（26）即\mathit{a}的表达式为基础设计的自适应律具有稳定性。

\mathit{d}与模型参数关系不明显，而\mathit{a}与模型参数关系明显，下面将讨论\mathit{a}与模型参数的关系。

当系统完成参数自整定时有：

{\mathit{e}}_{\mathrm{2}}={\mathit{B}}_k{\mathit{I}}_{\mathrm{3}\times \mathrm{1}}-{\mathit{B}}_{k\mathrm{r}\mathrm{e}\mathrm{a}}\mathit{a}=\mathrm{0}（29）

在系统启动时，自适应律还没有发挥作用，此时预测模型的参数为名义参数，所以矩阵\mathit{a}(\mathrm{0})取各元素为1的对角矩阵。

联立式（29）与（26），得到真实值与名义值的关系为：

{\mathit{B}}_k{\mathit{I}}_{\mathrm{3}\times \mathrm{1}}={\mathit{B}}_{k\mathrm{r}\mathrm{e}\mathrm{a}}({\mathit{I}}_{\mathrm{3}\times \mathrm{1}}+{\displaystyle \sum _{i=\mathrm{1}}^k}{\mathit{\psi }}_{\mathrm{1}}{\mathit{\psi }}_{\mathrm{2}}^{-\mathrm{1}}{\mathit{I}}_{\mathrm{3}\times \mathrm{1}}\times \tilde{\mathit{u}}(i)\times \mathit{e}(i))（30）

式（30）中，质量的自适应律由两个因素决定，将式（30）展开并改写为式（31），从而保证这两个因素能根据相应偏差的大小在自适应律中发挥相应比重的作用。

\left\{\begin{array}{l}{I}_{z\mathrm{r}\mathrm{e}\mathrm{a}}=I_z[\mathrm{1}+{\displaystyle \sum _{i=\mathrm{1}}^k}{\psi }_{\mathrm{13}}{\psi }_{\mathrm{23}}^{-\mathrm{1}}{e}_{\mathrm{3}}(i){\tilde{u}}_{\mathrm{3}}(i)]\\ m_{\mathrm{r}\mathrm{e}\mathrm{a}}=\frac{e_{\mathrm{1}}^{\mathrm{2}}(k)m[\mathrm{1}+{\displaystyle \sum _{i=\mathrm{1}}^k}{\psi }_{\mathrm{11}}{\psi }_{\mathrm{21}}^{-\mathrm{1}}{e}_{\mathrm{1}}(i){\tilde{u}}_{\mathrm{1}}(i)]}{e_{\mathrm{1}}^{\mathrm{2}}(k)+e_{\mathrm{2}}^{\mathrm{2}}(k)}+\\ \frac{e_{\mathrm{2}}^{\mathrm{2}}(k)m[\mathrm{1}+{\displaystyle \sum _{i=\mathrm{1}}^k}{\psi }_{\mathrm{12}}{\psi }_{\mathrm{22}}^{-\mathrm{1}}{e}_{\mathrm{2}}(i){\tilde{u}}_{\mathrm{2}}(i)]}{e_{\mathrm{1}}^{\mathrm{2}}(k)+e_{\mathrm{2}}^{\mathrm{2}}(k)}\end{array}\right.（31）

式中：{\tilde{u}}_{\mathrm{1}}\mathrm{、}{\tilde{u}}_{\mathrm{2}}\mathrm{、}{\tilde{u}}_{\mathrm{3}}为模型预测控制器的控制输入；e_{\mathrm{1}}\mathrm{、}{e}_{\mathrm{2}}\mathrm{、}{e}_{\mathrm{3}}分别为纵向速度、横向速度、横摆角速度的上一时刻预测状态和实际的偏差。

式（31）为本节设计的自适应律。在此自适应律作用下，随着时间的推移，式（30）得以满足，从而修正名义模型的参数。修正后，依据式（16）可知，模型预测控制方法的预测状态与实际被控对象的状态一致，从而优化控制效果。此外，本自适应律在一个采样周期内计算一轮，相较于模型预测控制中使用的内点法，复杂度没有显著增加。

6 地面附着系数辨识

模型预测控制器中的约束和地面附着系数\mu关系密切，且\mu的值难以通过直接测量确定，见式（5），若不获得准确的附着系数，则无法正确描述控制输入约束，本文将基于Brush轮胎模型^［14］和递归最小二乘法对\mu进行估算。

Brush模型描述了轮胎纵向力、侧向力、滑移率和侧偏角之间的复杂关系，如式（32）所示：

\left\{\begin{array}{l}{F}_{\mathrm{f}i}=\left\{\begin{array}{l}\begin{array}{c}{C}_{\mathrm{\alpha }i}{\sigma }_i-\frac{\mathrm{1}}{\mathrm{3}}\frac{C_{\mathrm{\alpha }i}^{\mathrm{2}}}{u{F}_{zi}}+\frac{\mathrm{1}}{\mathrm{27}}\frac{C_{\mathrm{\alpha }i}^{\mathrm{3}}{\sigma }_i^{\mathrm{3}}}{u^{\mathrm{2}}{F}_{zi}^{\mathrm{2}}},{\sigma }_i\le \frac{\mathrm{3}u{F}_{zi}}{C_{\mathrm{\alpha }i}}\end{array}\\ \begin{array}{c}u{F}_{zi},{\sigma }_i\le \frac{\mathrm{3}u{F}_{zi}}{C_{\mathrm{\alpha }i}}\end{array}\end{array}\right.\\ F_{xi}=\frac{k_{xi}{\sigma }_{xi}}{{\sigma }_i}{F}_{\mathrm{f}i}\\ F_{yi}=\frac{{\sigma }_{yi}}{{\sigma }_i}{F}_{\mathrm{f}i}\end{array}\right.（32）

式（32）中的各个参数可通过式（33）获得：

\left\{\begin{array}{l}{\sigma }_{yi}=\frac{V_{\mathrm{w}xi}}{R_{\mathrm{e}i}{w}_i}\mathrm{t}\mathrm{a}\mathrm{n}{\alpha }_i,{\sigma }_i=\sqrt[]{{\sigma }_{xi}^{\mathrm{2}}+{\sigma }_{yi}^{\mathrm{2}}}\\ k_i=\frac{C_{\mathrm{s}i}}{C_{\mathrm{\alpha }i}}=\frac{{\sigma }_{yi}{F}_{xi}}{{\sigma }_{xi}{F}_{yi}},{\sigma }_{xi}=s_i\end{array}\right.（33）

式中：V_{wxi}为轮胎的平移速度在车身坐标轴横向方向上的分量；s_i为轮胎纵向滑移率，实际应用中，轮胎滑移率已经可以通过ABS获得。

当轮胎力没有饱和时，将式（32）整理为观测模型：

F_{\mathrm{f}i}=\sqrt[]{\frac{F_{xi}^{\mathrm{2}}}{k_i^{\mathrm{2}}}+F_{yi}^{\mathrm{2}}}=\mathit{H}\mathit{X}（34）

其中，状态观测矩阵\mathit{H}=\left[\begin{array}{ccc}{\sigma }_i& -\raisebox{1ex}{${\sigma }_i^{\mathrm{2}}$}\!\left/ \!\raisebox{-1ex}{$\mathrm{3}{F}_{zi}$}\right.& \raisebox{1ex}{${\sigma }_i^{\mathrm{3}}$}\!\left/ \!\raisebox{-1ex}{$\mathrm{27}{F}_{zi}^{\mathrm{2}}$}\right.\end{array}\right]，待估计参数矩阵\mathit{X}={\left[\begin{array}{ccc}{C}_{\mathrm{\alpha }i}& \raisebox{1ex}{$C_{\mathrm{\alpha }i}^{\mathrm{2}}$}\!\left/ \!\raisebox{-1ex}{$\mu $}\right.& \raisebox{1ex}{$C_{\mathrm{\alpha }i}^{\mathrm{3}}$}\!\left/ \!\raisebox{-1ex}{${\mu }^{\mathrm{2}}$}\right.\end{array}\right]}^{\mathrm{T}}，\mathit{H}中相关参数和测量值{\widehat{F}}_{\mathrm{f}\mathrm{l}}可通过观测车辆系统得到，其中的轮胎横纵向力可通过轮胎位移加速度方向判断摩擦力方向，再将摩擦力分解至轮胎坐标系得到。

根据观测模型式（34），使用递归最小二乘法进行参数辨识，计算过程涉及式（35）：

\left\{\begin{array}{l}\mathit{K}\left(k\right)=\frac{\mathit{P}(k-\mathrm{1})\mathit{H}(k)}{\mathrm{1}+\mathit{H}(k)\mathit{P}(k-\mathrm{1}){\mathit{H}}^{\mathrm{T}}(k)}\\ \mathit{X}(k)=\mathit{X}(k-\mathrm{1})+\mathit{K}(k)({\widehat{F}}_{\mathrm{f}l}(k)-\mathit{H}(k)\mathit{X}(k))\\ \mathit{P}(k)=(\mathit{I}-\mathit{K}(k)\mathit{H}(k))\mathit{P}(k-\mathrm{1})\end{array}\right.（35）

式中：\mathit{I}为三维单位矩阵；记\mathit{X}={\left[\begin{array}{ccc}{X}_{\mathrm{1}}& X_{\mathrm{2}}& X_{\mathrm{3}}\end{array}\right]}^{\mathrm{T}}。

递归最小二乘法计算过程如下：

（1）初始化初始估计值{\mathit{X}}_{\mathrm{0}}和协方差矩阵{\mathit{P}}_{\mathrm{0}}；

（2）根据式（35）更新增益矩阵\mathit{K}(k)；

（3）根据式（35）更新估计参数矩阵\mathit{X}(k)；

（4）根据式（35）更新协方差矩阵\mathit{P}(k)；

（5）通过{\mu }^{\mathrm{*}}=\frac{X_{\mathrm{1}}{X}_{\mathrm{2}}}{X_{\mathrm{3}}}获得估计附着系数{\mu }^{\mathrm{*}}。下一采样时刻返回第（2）步。

通过估计\mu的数值，可准确估计汽车驱动能力，完善对控制器约束条件式（5）的描述，递归最小二乘法在一个采样周期内只进行一轮计算，相较于模型预测控制中进行求解最优问题使用的内点法，复杂度没有显著增加。

7 仿真验证及结果分析

本文基于Simulink-Carsim联合仿真平台，使用传统模型预测控制器和参数自适应模型预测控制器控制车辆行驶进行对比仿真验证，验证自适应模型预测控制器在模型参数存在误差情况下的控制效果。仿真采用双移线工况，路面附着系数取0.8。

表1中的名义质量与转动惯量为车辆空车时的参数，并作为名义模型预测控制器的参数。名义路面附着系数与侧偏刚度的关系，见式（34），作为最小二乘法的初值{\mathit{X}}_{\mathrm{0}}。表2中的整车质量和转动惯量则使用载重情况下的参数。控制器的权重矩阵和其他参数取值如下：

\left\{\begin{array}{l}\mathit{Q}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\{[\begin{array}{ccc}\mathrm{1}\mathrm{000}& \mathrm{60}\mathrm{000}& \mathrm{70}\mathrm{000}\end{array}]\}\\ \mathit{R}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\{[\begin{array}{ccc}\mathrm{10}& \mathrm{10}& \mathrm{10}\end{array}]\}\\ {\mathit{\psi }}_{\mathrm{1}}={\mathit{\psi }}_{\mathrm{2}}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\{[\begin{array}{ccc}\mathrm{1}& \mathrm{1}& \mathrm{1}\end{array}]\}\\ {\mathit{P}}_{\mathrm{0}}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\{[\begin{array}{ccc}\mathrm{50}\mathrm{000}& \mathrm{50}\mathrm{000}& \mathrm{50}\mathrm{000}\end{array}]\}\\ \mathit{d}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\{[\begin{array}{ccc}m& m& I_z\end{array}]\}\times \mathrm{0.8}/T\end{array}\right.

一般的汽车的质量在1 000~1 500 kg，本文取控制器参数（名义整车质量）为1 000 kg。载重车辆质量为1 400 kg，为真实模型参数。初始名义附着系数选择0.2，使辨识附着系数前的控制策略相对保守。

将选取的参数{\mathit{\psi }}_{\mathrm{1}}\mathrm{、}{\mathit{\psi }}_{\mathrm{2}}\mathrm{、}\mathit{d}及设定的行车状态纵向速度80 km/h、横向速度0、横摆角速度0.3 rad/s、反馈控制量\tilde{\mathit{u}}={[\begin{array}{ccc}\mathrm{500}& \mathrm{1}\mathrm{000}& \mathrm{1}\mathrm{500}\end{array}]}^{\mathrm{T}}，汽车空载与满载参数代入式（28）得到：

\left\{\begin{array}{l}\mathrm{\Delta }V\approx {\mathit{e}}^{\mathrm{T}}(k)\left[\begin{array}{ccc}-\mathrm{0.45}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& -\mathrm{0.43}& -\mathrm{0.17}\\ \mathrm{0}& -\mathrm{0.11}& -\mathrm{0.47}\end{array}\right]\mathit{e}(k)\\ \mathrm{\Delta }V\approx {\mathit{e}}^{\mathrm{T}}(k)\left[\begin{array}{ccc}-\mathrm{1.53}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& -\mathrm{1.50}& -\mathrm{0.51}\\ \mathrm{0}& -\mathrm{0.11}& -\mathrm{1.56}\end{array}\right]\mathit{e}(k)\end{array}\right.（36）

式（36）分别为向式（28）代入表中m\mathrm{、}{I}_z与m_{\mathrm{r}\mathrm{e}\mathrm{a}}\mathrm{、}{I}_{z\mathrm{r}\mathrm{e}\mathrm{a}}作为真值的结果，其值恒为负。说明对于不具有载重和具有载重的汽车，此参数都能使辨识参数稳定在真值附近。

将期望纵向速度\横向速度和横摆角速度作为系统的输入信号，验证参数自适应模型预测控制器效果相比于传统模型预测控制器，在模型参数不确定的情况下具有更好的控制性能。

匀加速工况下的两种控制器的控制效果如图3所示。此工况下，汽车在4 s内进行匀加速运动从静止到80 km/h，之后维持匀速运动。“a-mpc”为本文设计的自适应模型预测控制器的控制效果。“mpc”与“lqr”为对应控制器的控制效果。

从图3可以看出，由于参数辨识需要时间，系统启动后的0.5 s内的三者控制器效果基本一致。MPC和LQR控制器的控制效果存在静差，直到4 s后静差才逐渐消失，而自适应MPC在匀加速阶段即可消除静差，在4 s时的超调也较低。

双移线工况下的控制效果如图4所示。在开始的1 s内，在参数自适应基本完成之前，这段时间内由MPC、LQR和自适应MPC控制得到的控制效果几乎相同。1 s后，自适应律对参数修正基本完成，偏差减少。自适应MPC使质心侧偏角更小，保证车辆更好的操纵性。原因是由于名义质量参数相比于真实质量较小，传统模型预测会计算出为了维持转弯不发生偏移的情况所需的侧向力较小，但实际上需要的侧向力比控制输入大，从而导致了侧滑的发生。从行驶路径来看，MPC和LQR下路径的最大偏移量超过2 m，考虑到一般情况下，车身宽度在1.8 m左右，高速公路单车道宽为3.75 m，假设行驶路径在道路中线，则车身边缘距离单车道边缘约1 m，在最大偏差达到2 m时，会导致车辆的行驶轨迹超出单车道从而影响安全行驶。使用本文提出的参数自适应模型预测控制方法后，车辆的行驶状态与期望基本一致，能够保证行车安全。表3为各工况下各控制期望与实际状态的平均绝对误差，可以看到自适应MPC的性能相较于MPC和LQR更好。

附着系数辨识效果图如图5所示，\mu表示附着系数真值，{\mu }_i表示各轮胎处的估计附着系数，i表示相应轮胎。可以看出，在系统启动时，估计附着系数向真值迅速靠近，其中左前方轮胎的估计附着系数经过大约1 s的上升与超调过程后最终收敛于真值，其他轮胎的估计附着系数经过小于1 s的上升过程后也接近于真值。1 s的时间经过了20个采样周期，即迭代20次可得到较为准确的估计值，实验结果表明，本文设计的参数辨识方法是有效的。

8 结束语

汽车行驶过程中，由于载重变化和路面情况的不同，车辆参数的变化可能导致控制器性能下降。本文针对这一问题，一方面提出了一种模型参数自适应方法，减少预测状态与实际状态之间的偏差。该方法通过补偿控制输入，使预测模型能更准确反映真实模型的状态响应。另一方面，设计参数辨识方法在线估计轮胎与路面的附着系数，并修正控制输入约束。将两者与模型预测控制方法结合，可以改善参数误差情况下的控制效果并获取到相应参数。通过物理仿真实验验证了本文设计的有效性，结果显示本文设计的方法能够在模型参数不准确时保证控制的性能，从而保证行车安全。

参考文献

[1] Zhu S, Wei B, Liu D, et al. A dynamics coordinated control system for 4WD-4WS electric vehicles[J]. Electronics, 2022, 11(22): 410-411.

[2] Jiang C, Tian H, Hu J, et al. Learning based predictive error estimation and compensator design for autonomous vehicle path tracking[C]∥15th IEEE Conference on Industrial Electronics and Applications, Singapore, 2020: 1496-1500.

[3] Mata S, Zubizarreta A, Pinto C. Robust tube-based model predictive control for lateral path tracking[J]. IEEE Transactions on Intelligent Vehicles, 2019, 4(4): 569-577.

[4] Yamauchi Y, Saito M, Ono T. Adaptive identification method of vehicle modeling according to the fluctuation of road and running situation in autonomous driving[C]∥58th Annual Conference of the Society of Instrument and Control Engineers of Japan, Tokyo, Japan, 2019: 1255-1260.

[5] Kim T, Kim H J. Path tracking control and identification of tire parameters using on-line model-based reinforcement learning[C]∥16th International Conference on Control, Automation and Systems, Seoul, South Korea, 2016: 215-219.

[6] Lee J, Hwang Y, Choi S B. Robust tube-MPC based steering and braking control for path tracking at high-speed driving[J]. IEEE Transactions on Vehicular Technology, 2023, 72: 15301-15316.

[7] Su S, Hao C, Weaver C, et al. Double-iterative Gaussian process regression for modeling error compensation in autonomous racing[J]. IFAC-PapersOnLine, 2023, 56(2): 7940-7947.

[8] 刘建泽, 柳江, 李敏, 等. 基于最小二乘的车速解耦路面辨识方法[J]. 吉林大学学报: 工学版, 2024, 54(7): 1821-1830.

[9] Liu Jian-ze, Liu Jiang, Li Min, et al. Vehicle speed decoupling road identification method based on least squares[J]. Journal of Jilin University (Engineering and Technology Edition), 2024, 54(7): 1821-1830.

[10] 李杰, 郭文翠, 赵旗, 等. 基于车辆响应的路面不平度识别方法[J]. 吉林大学学报: 工学版, 2019, 49(6): 1810-1817.

[11] Li Jie, Guo Wen-cui, Zhao Qi, et al. Road roughness identification method based on vehicle response[J]. Journal of Jilin University (Engineering and Technology Edition), 2019, 49(6): 1810-1817.

[12] Duan G. High-order fully actuated system approaches: part V. robust adaptive control[J]. International Journal of Systems Science, 2021, 52(10): 2129-2143.

[13] Pius P E, Selekwa M F. Dynamic modelling of a 4WD/4WS ground vehicle by using Gibbs-Appell approach[C]∥ASME International Mechanical Engineering Congress and Exposition, New York, USA, 2022: V005T07A071.

[14] 顾添骠. 基于驱动可行域的电动汽车执行器容错控制[D]. 长春: 吉林大学通信工程学院, 2024.

[15] Gu Tian-biao. Fault-tolerant control of electric vehicle actuator based on drive feasible domain[D]. Changchun: College of Communication Engineering, Jilin University, 2024.

[16] Schwenzer M, Ay M, Bergs T, et al. Review on model predictive control: an engineering perspective[J]. The International Journal of Advanced Manufacturing Technology, 2021, 117(5): 1327-1349.

[17] Agarwal R P, Bohner M, Özbekler A. Lyapunov Inequalities and Applications[M]. Cham: Springer, 2021.

[18] Jelkić Z, Pikula B. The establishment of an advanced brush model for simulation of vehicle dynamics[C]∥Proceedings of the International Conference "New Technologies, Development and Applications", Springer Cham, 2022: 260-273.

基金资助

吉林省科技发展计划项目(20220301014GX)