基于快速模型预测控制的重型车辆侧翻稳定性研究

Research on Heavy Vehicle Rollover Stability Based on Fast Model Predictive Control

摘要

为解决重型车辆易发生侧翻的问题,提出一种基于快速模型预测的主动安全控制方法.首先,建立重型车辆非线性动力学模型,利用T-S模糊方法将非线性系统转化为线性系统;采用分层控制策略,上层设计模型预测控制（model predictive control,MPC）防侧翻控制器获得最优横摆力矩,下层利用单侧车轮控制方式进行制动力分配.然后,为减少控制器的计算负担,引入Laguerre函数对控制输入序列进行重构,采用原始对偶内点法处理约束条件,实时高效地求解二次规划函数.最后,在典型工况下进行MATLAB/Simulink和TruckSim联合仿真.结果表明,相较于传统MPC控制器,所提的快速MPC控制器能够在保证性能的同时有效提高计算效率,可以及时避免车辆侧翻失稳.

Abstract

In order to solve the problem that heavy vehicles are prone to rollover， an active safety control method based on fast model prediction was proposed. Firstly， a nonlinear dynamic model of the heavy vehicle was established， and the nonlinear system was transformed into a linear system by the T-S fuzzy method. A hierarchical control strategy was adopted， where the upper layer was designed to obtain the optimal yaw moment by the model predictive control （MPC） anti-rollover controller， and the lower layer utilized the unilateral wheel control method for the allocation of the braking force. Then， in order to reduce the computational burden of the controller， the Laguerre function was introduced to reconstruct the control input sequence， and the constraints were handled by the primal-dual interior point method to solve the quadratic programming function efficiently in real time. Finally， the co-simulation of MATLAB/Simulink and TruckSim was carried out under typical working conditions， and the results show that compared with the traditional MPC controller， the proposed fast MPC controller can effectively improve the computational efficiency while guaranteeing the performance and can avoid the vehicle rollover instability.

Graphical abstract

关键词

重型车辆 / 模型预测控制 / Laguerre函数 / 原始对偶内点法 / 侧翻稳定性

Key words

heavy vehicle / model predictive control / Laguerre function / primal-dual interior-point method / rollover stability

重型特种车辆在军事国防方面发挥了重要作用，但由于其质心高、质量大的结构特点，导致其在行驶过程中更容易发生侧翻^［1］.侧翻是一种极其严重的事故，这会使得车辆瞬态失稳，一旦发生就可能导致严重人员伤亡和财产损失.因此，在车辆运行过程中，对侧翻稳定性进行控制具有十分重要的意义^［2］.

近年来，已有多种控制算法应用在重型车辆防侧翻领域^［3-4］，例如PID控制、滑模控制、鲁棒控制、模糊控制等.虽然上述控制算法可以成功防止车辆侧翻，但难以处理多约束问题.MPC作为一种可以有效处理上述问题的算法，近年来在车辆稳定性控制方面有了越来越多的应用.文献［5］设计了一种分层防侧翻控制策略，采用MPC算法有效防止车辆侧翻.但是传统MPC需要通过滚动优化来确定当前的最优控制策略，计算负担大，减少MPC计算量的问题亟待解决.文献［6-7］提出了一种基于Laguerre函数的模型预测控制，重构控制输入信号，减小系统矩阵维度，提高MPC的计算效率.文献［8］为了提高MPC控制器的动态响应速度，在实际车辆中采用了Laguerre函数，快速求得最优控制输入，及时对车辆的横摆稳定性进行控制.虽然引入Laguerre函数可以提高传统MPC控制器的响应速度，但是二次规划的求解部分仍然比较复杂.文献［9-10］针对MPC控制方法在线计算量大的问题，采用原始对偶内点（primal-dual interior-point，PDIP）法来计算二次规划部分，实时高效地处理了约束条件，通过求解搜索路径和搜索步长快速迭代，求得最优解.

本文建立了纵向车速时变的三自由度车辆侧翻模型，采用T-S模糊方法进行线性化处理.设计了分层控制策略：上层采用MPC控制器，获得维持车辆侧翻稳定性所需的最优横摆力矩，引入Laguerre函数和原始对偶内点法，并将两者有机结合，既对输入序列进行重构，减小了系统矩阵的维度，又实时高效地处理约束条件，加快了MPC控制器的响应速度；下层采用单侧分配方式对附加横摆力矩进行分配.选取双移线和鱼钩工况，通过MATLAB/Simulink与TruckSim联合仿真，验证所提快速MPC算法的有效性.

1 系统模型

1.1 车辆动力学模型

忽略空气动力学的影响，将模型进行简化，建立包含侧向、侧倾、横摆运动的三自由度车辆动力学模型，如图1所示^［11］.根据上述模型，车辆的动力学方程表示形式如下：

\left.\begin{array}{c}{I}_x\ddot{\varphi }=m{a}_{y}h\mathrm{c}\mathrm{o}\mathrm{s}\varphi +m_{\mathrm{s}}gh\mathrm{s}\mathrm{i}\mathrm{n}\varphi -c_{\mathrm{\phi }}\dot{\varphi }-k_{\mathrm{\phi }}\varphi ,\\ m{a}_y-m_{\mathrm{s}}h\ddot{\varphi }=F_{\mathrm{f}}\mathrm{c}\mathrm{o}\mathrm{s}\delta +F_{\mathrm{r}},\\ I_z\dot{\gamma }=F_{\mathrm{f}}a\mathrm{c}\mathrm{o}\mathrm{s}\delta -F_{\mathrm{r}}b+M_z.\end{array}\right\}(1)

式中：a_y={\dot{v}}_y+v_x\gamma为车辆的侧向加速度；I_x为整车绕x轴的转动惯量；I_z为整车绕z轴的转动惯量；v_x和v_y分别为车辆系统的纵向速度和横向速度；m为整车质量；h为车辆质心到侧倾中心的距离；m_{\mathrm{s}}为簧载质量；k_{\mathrm{\phi }}和c_{\mathrm{\phi }}分别为车辆的等效侧倾刚度和等效阻尼系数；F_{\mathrm{f}}和F_{\mathrm{r}}分别为前、后轮胎的侧偏力；\gamma为横摆角速度；a和b分别为车辆质心到前、后轴的距离；\varphi为侧翻角；\delta为前轮转角；M_z为横摆力矩.

前、后轮胎的侧向力表示为

\left.\begin{array}{c}{F}_{\mathrm{f}}=k_{\mathrm{f}}\cdot {\alpha }_{\mathrm{f}},\\ F_{\mathrm{r}}=k_{\mathrm{r}}\cdot {\alpha }_{\mathrm{r}}.\end{array}\right\}(2)

式中：k_{\mathrm{f}}，k_{\mathrm{r}}分别代表前轮与后轮的侧偏刚度；{\alpha }_{\mathrm{f}}，{\alpha }_{\mathrm{r}}分别代表前轮与后轮的侧偏角，对应的几何关系如下：

\left.\begin{array}{c}{\alpha }_{\mathrm{f}}=\delta -\beta -\frac{a\gamma }{v_x},\\ {\alpha }_{\mathrm{r}}=\frac{b\gamma }{v_x}-\beta .\end{array}\right\}(3)

式中，\beta为质心侧偏角.

假设侧翻角很小，\mathrm{s}\mathrm{i}\mathrm{n}\varphi =\varphi，\mathrm{c}\mathrm{o}\mathrm{s}\varphi =\mathrm{1}，将式（2）和式（3）代入式（1）中，可得车辆三自由度动力学模型的状态空间方程：

\dot{\mathit{x}}=\mathit{A}\mathit{x}+\mathit{B}\mathit{u}.(4)

式中：\mathit{x}={\left[v_y,\gamma ,\varphi ,\dot{\varphi }\right]}^{\mathrm{T}}为系统状态；\mathit{u}={\left[\delta ,M_z\right]}^{\mathrm{T}}为系统输入；系数矩阵\mathit{A}和\mathit{B}如下所示：

\left.\begin{array}{l}{A}_{\mathrm{11}}=-\frac{\left(k_{\mathrm{f}}+k_{\mathrm{r}}\right)I_x}{(m{I}_x-m_{\mathrm{s}}^{\mathrm{2}}{h}^{\mathrm{2}})v_x},A_{\mathrm{12}}=\frac{\left(k_{\mathrm{r}}b-k_{\mathrm{f}}a\right)I_x}{(m{I}_x-m_{\mathrm{s}}^{\mathrm{2}}{h}^{\mathrm{2}})v_x}-v_x,\\ A_{\mathrm{13}}=\frac{m_{\mathrm{s}}h\left(m_{\mathrm{s}}gh-k_{\mathrm{\phi }}\right)}{m{I}_x-m_{\mathrm{s}}^{\mathrm{2}}{h}^{\mathrm{2}}},A_{\mathrm{14}}=\frac{-m_{\mathrm{s}}h{c}_{\mathrm{\phi }}}{m{I}_x-m_{\mathrm{s}}^{\mathrm{2}}{h}^{\mathrm{2}}},\\ A_{\mathrm{21}}=\frac{k_{\mathrm{r}}b-k_{\mathrm{f}}a}{I_z{v}_x},A_{\mathrm{22}}=\frac{k_{\mathrm{f}}{a}^{\mathrm{2}}+k_{\mathrm{r}}{b}^{\mathrm{2}}}{I_z{v}_x},A_{\mathrm{23}}=\mathrm{0},A_{\mathrm{24}}=\mathrm{0},\\ A_{\mathrm{31}}=\mathrm{0},A_{\mathrm{32}}=\mathrm{0},A_{\mathrm{33}}=\mathrm{0},A_{\mathrm{34}}=\mathrm{1},\\ A_{\mathrm{41}}=\frac{m_{\mathrm{s}}h\left(k_{\mathrm{f}}+k_{\mathrm{r}}\right)}{(m{I}_x-m_{\mathrm{s}}^{\mathrm{2}}{h}^{\mathrm{2}})v_x},A_{\mathrm{42}}=\frac{m_{\mathrm{s}}h\left(k_{\mathrm{f}}a-k_{\mathrm{r}}b\right)}{(m{I}_x-m_{\mathrm{s}}^{\mathrm{2}}{h}^{\mathrm{2}})v_x},\\ A_{\mathrm{43}}=\frac{m(k_{\mathrm{\phi }}-m_{\mathrm{s}}gh)}{m{I}_x-m_{\mathrm{s}}^{\mathrm{2}}{h}^{\mathrm{2}}},A_{\mathrm{44}}=\frac{-m{c}_{\mathrm{\phi }}}{m{I}_x-m_{\mathrm{s}}^{\mathrm{2}}{h}^{\mathrm{2}}}.\end{array}\right\}(5)

\left.\begin{array}{l}{B}_{\mathrm{11}}=\frac{k_{\mathrm{f}}{I}_x}{(m{I}_x-m_{\mathrm{s}}^{\mathrm{2}}{h}^{\mathrm{2}})},B_{\mathrm{12}}=\mathrm{0},B_{\mathrm{21}}=\frac{k_{\mathrm{f}}a}{I_z},B_{\mathrm{22}}=\frac{\mathrm{1}}{I_z},\\ B_{\mathrm{31}}=\mathrm{0},B_{\mathrm{32}}=\mathrm{0},B_{\mathrm{41}}=\frac{m_{\mathrm{s}}h{k}_{\mathrm{f}}}{(m{I}_x-m_{\mathrm{s}}^{\mathrm{2}}{h}^{\mathrm{2}})},B_{\mathrm{42}}=\mathrm{0}.\end{array}\right\}(6)

1.2 变速度车辆模型T-S模糊线性化

状态空间方程（4）中的纵向车速是时变的，模糊规则如下：

\left.\begin{array}{l}\dot{\mathit{x}}\left(t\right)={\mathit{A}}_i\mathit{x}\left(t\right)+\mathit{B}\mathit{u}\left(t\right).\\ \mathrm{s}.\mathrm{t}.\mathit{v}(t)=\{v_{\mathrm{1}}\left(t\right),v_{\mathrm{2}}\left(t\right),\dots ,v_q\left(t\right)\}\subset F_{ij},\\ i=\mathrm{1,2},\dots ,r,j=\mathrm{1,2},\dots ,q.\end{array}\right\}(7)

式中：{\mathit{A}}_i为车辆状态矩阵；\mathit{v}\left(t\right)为车辆速度；F_{ij}为模糊集；q为车辆速度分量的数量；r为if-then规则的数量.

V_i\left(\mathit{v}\left(t\right)\right)={\displaystyle \prod _{j=\mathrm{1}}^q}{F}_{ij}\left(v_j\left(t\right)\right).(8)

式中，V_i\left(\mathit{v}\left(t\right)\right)表示v_j\left(t\right)在F_{ij}中的隶属度.

隶属度函数表示为

O_i\left(\mathit{v}\left(t\right)\right)=\frac{V_i\left(\mathit{v}\left(t\right)\right)}{{\displaystyle \sum _{i=\mathrm{1}}^r}{V}_i\left(\mathit{v}\left(t\right)\right)}.(9)

建立车辆的T-S模糊线性化模型：

\left.\begin{array}{c}\dot{\mathit{x}}={\displaystyle \sum _{i=\mathrm{1}}^r}{O}_i\left(\mathit{v}\left(t\right)\right)\left({\mathit{A}}_i\mathit{x}+\mathit{B}\mathit{u}\right),\\ \mathit{y}=\mathit{C}\mathit{x}.\end{array}\right\}(10)

式中：\mathit{C}=\left[\begin{array}{cccc}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]； y 为输出矩阵.

2 防侧翻控制系统设计

针对重型车辆侧翻稳定性问题，提出一种分层控制策略.选择横向载荷转移率（lateral load transfer ratio， LTR）作为侧翻评价指标，设计侧翻预警系统.上层采用MPC控制器计算维持车辆侧翻稳定性所需的横摆力矩；下层采用单侧方式分配力矩.随后，引入Laguerre函数和原始对偶内点法，以减小MPC控制器的计算量.侧翻稳定性的整体控制框图如图2所示.

2.1 侧翻预警系统

LTR反映车辆侧翻风险的大小，其定义如下^［12］：

\mathrm{L}\mathrm{T}\mathrm{R}=\frac{F_{z\mathrm{r}}-F_{z\mathrm{l}}}{F_{z\mathrm{r}}+F_{z\mathrm{l}}}(11)

由于F_{z\mathrm{l}}和F_{z\mathrm{r}}难以进行实时测量，以图1中侧倾中心O为参考点列力矩平衡方程：

将LTR转换为容易测量的形式，得

\left.\begin{array}{l}{F}_{z\mathrm{l}}\frac{B_{\mathrm{w}}}{\mathrm{2}}+m_{\mathrm{s}}{a}_{y}h\mathrm{c}\mathrm{o}\mathrm{s}\varphi +m_{\mathrm{s}}gh\mathrm{s}\mathrm{i}\mathrm{n}\varphi -F_{z\mathrm{r}}\frac{B_{\mathrm{w}}}{\mathrm{2}}=\mathrm{0},\\ F_{z\mathrm{l}}+F_{z\mathrm{r}}=mg,\\ h\mathrm{c}\mathrm{o}\mathrm{s}\varphi \approx h,\\ h\mathrm{s}\mathrm{i}\mathrm{n}\varphi \approx h\varphi .\end{array}\right\}(12)

将式（12）代入式（11）中，得

\mathrm{L}\mathrm{T}\mathrm{R}=\frac{\mathrm{2}\left({\dot{v}}_y+v_x\gamma \right)h-\mathrm{2}g\varphi h}{g{B}_{\mathrm{w}}}.(13)

LTR是随速度变化的函数，因此，需要研究LTR在时变速度条件下的预警情况.当LTR大于安全阈值时，车辆具有不稳定倾向，控制器需进行防侧翻控制.

2.2 防侧翻控制器设计

选取二自由度车辆模型的横摆角速度{\gamma }_{\mathrm{s}}和质心侧偏角{\beta }_{\mathrm{s}}的稳态响应值作为参考值^［13］：

{\gamma }_{\mathrm{s}}=\frac{v_x\delta }{\left(a+b\right)(\mathrm{1}+k_{\mathrm{u}\mathrm{s}}{v}_x^{\mathrm{2}})},(14)

{\beta }_{\mathrm{s}}=\frac{l_{\mathrm{r}}\delta }{\left(a+b\right)(\mathrm{1}+k_{\mathrm{u}\mathrm{s}}{v}_x^{\mathrm{2}})}+\frac{m{l}_{\mathrm{f}}{v}_x\delta }{k_{\mathrm{f}}\left(a+b\right)(\mathrm{1}+k_{\mathrm{u}\mathrm{s}}{v}_x^{\mathrm{2}})}.(15)

式中：k_{\mathrm{u}\mathrm{s}}=\frac{m}{\left(a+b\right)}\left(\frac{a}{k_{\mathrm{f}}}-\frac{b}{k_{\mathrm{r}}}\right)为转向稳定性系数；l_{\mathrm{f}}，l_{\mathrm{r}}为前后轴到车辆质心的距离.

考虑路面最大附着系数的影响，可得二者最大值：

{\gamma }_{\mathrm{m}\mathrm{a}\mathrm{x}}=\frac{\mu g}{v_x},(16)

{\beta }_{\mathrm{m}\mathrm{a}\mathrm{x}}=\mu g\left[\frac{b}{v_x^{\mathrm{2}}}+\frac{ma}{k_{\mathrm{2}}\left(a+b\right)}\right].(17)

式中，\mu为路面附着系数.

由此可计算得横摆角速度{\gamma }_{\mathrm{d}}、质心侧偏角{\beta }_{\mathrm{d}}、横向速度v_{y\mathrm{d}}的理想值如下：

{\gamma }_{\mathrm{d}}=\mathrm{m}\mathrm{i}\mathrm{n}\left(\left|{\gamma }_{\mathrm{s}}\right|,\left|{\gamma }_{\mathrm{m}\mathrm{a}\mathrm{x}}\right|\right)\cdot \mathrm{s}\mathrm{g}\mathrm{n}\left({\gamma }_{\mathrm{s}}\right),(18)

{\beta }_{\mathrm{d}}=\mathrm{m}\mathrm{i}\mathrm{n}\left(\left|{\beta }_{\mathrm{s}}\right|,\left|{\beta }_{\mathrm{m}\mathrm{a}\mathrm{x}}\right|\right)\cdot \mathrm{s}\mathrm{g}\mathrm{n}\left({\beta }_{\mathrm{s}}\right),(19)

{v_y}_{\mathrm{d}}=v_x{\beta }_{\mathrm{d}}.(20)

将式（10）进行离散化处理，可得

\left.\begin{array}{c}\mathit{x}\left(k+\mathrm{1}\right)={\displaystyle \sum _{i=\mathrm{1}}^r}{O}_i\left(\mathit{v}\left(k\right)\right)\left({\mathit{A}}_{\mathrm{M}}\mathit{x}\left(k\right)+{\mathit{B}}_{\mathrm{M}}\mathit{u}\left(k\right)\right),\\ \mathit{y}\left(k\right)=\mathit{C}\mathit{x}\left(k\right).\end{array}\right\}(21)

式中：\mathit{x}\left(k+\mathrm{1}\right)=\left[\begin{array}{c}{v}_y\left(k+\mathrm{1}\right)\\ \gamma \left(k+\mathrm{1}\right)\\ \varphi \left(k+\mathrm{1}\right)\\ \dot{\varphi }\left(k+\mathrm{1}\right)\end{array}\right]；\mathit{x}\left(k\right)=\left[\begin{array}{c}{v}_y\left(k\right)\\ \gamma \left(k\right)\\ \varphi \left(k\right)\\ \dot{\varphi }\left(k\right)\end{array}\right]；{\mathit{A}}_{\mathrm{M}}=\mathrm{1}+{\mathit{A}}_{\mathit{i}}T为离散状态转移矩阵；{\mathit{B}}_{\mathrm{M}}=T\mathit{B}为离散化后的输入矩阵；T为采样周期.

MPC控制器的目标函数为

\mathit{J}={{\displaystyle \sum _{z=\mathrm{1}}^{N_{\mathrm{p}}}}‖\mathit{y}\left(k+z\right)-{\mathit{y}}_{\mathrm{d}}\left(k+z\right)‖}_{\mathit{Q}}^{\mathrm{2}}+{{\displaystyle \sum _{z=\mathrm{0}}^{N_{\mathrm{c}}-\mathrm{1}}}‖\mathrm{\Delta }\mathit{u}\left(k+z\right)‖}_{\mathit{R}}^{\mathrm{2}}.(22)

式中：N_{\mathrm{p}}表示预测时域；N_{\mathrm{c}}表示控制时域；\mathit{y}为输出量预测值；{\mathit{y}}_{\mathrm{d}}为输出量期望值；\mathrm{\Delta }\mathit{u}为控制量增量；\mathit{Q}为输出权重矩阵；\mathit{R}为控制量增量权重矩阵；z为当前时刻到未来时刻的步数.

控制量\mathit{u}，控制量增量\mathrm{\Delta }\mathit{u}及输出量\mathit{y}的约束条件为

\left.\begin{array}{c}\mathrm{\Delta }{\mathit{u}}_{\mathrm{m}\mathrm{i}\mathrm{n}}\left(k+z\right)\le \mathrm{\Delta }\mathit{u}\left(k+z\right)\le \mathrm{\Delta }{\mathit{u}}_{\mathrm{m}\mathrm{a}\mathrm{x}}\left(k+z\right),\\ -{\mathit{u}}_{\mathrm{m}\mathrm{i}\mathrm{n}}\left(k+z\right)+\mathit{u}\left(k-\mathrm{1}\right)\le \mathit{u}\left(k+z\right)\le \\ -{\mathit{u}}_{\mathrm{m}\mathrm{a}\mathrm{x}}\left(k+z\right)+\mathit{u}\left(k-\mathrm{1}\right),\\ {\mathit{y}}_{\mathrm{m}\mathrm{i}\mathrm{n}}\left(k+z\right)\le \mathit{y}\left(k+z\right)\le {\mathit{y}}_{\mathrm{m}\mathrm{a}\mathrm{x}}\left(k+z\right).\end{array}\right\}(23)

将目标函数转换为易于求解的二次规划形式：

\left.\begin{array}{l}\mathrm{m}\mathrm{i}\mathrm{n}\mathit{J}=\frac{\mathrm{1}}{\mathrm{2}}\mathrm{\Delta }{\mathit{u}}^{\mathrm{T}}{\mathit{H}}_{\mathrm{1}}\mathrm{\Delta }\mathit{u}+{{\mathit{f}}_{\mathrm{1}}}^{\mathrm{T}}\mathrm{\Delta }\mathit{u},\\ \mathrm{s}.\mathrm{t}.{\mathit{A}}_{\mathrm{c}}\mathrm{\Delta }\mathit{u}\le {\mathit{b}}_{\mathrm{c}}.\end{array}\right\}(24)

式中：{\mathit{A}}_{\mathrm{c}}为二次规划约束矩阵；{\mathit{b}}_{\mathrm{c}}为约束向量；{\mathit{H}}_{\mathrm{1}}=\mathrm{2}{\mathit{S}}^{\mathrm{T}}\mathit{Q}\mathit{S}+\mathit{R}；{{\mathit{f}}_{\mathrm{1}}}^{\mathrm{T}}=\mathrm{2}{\mathit{S}}^{\mathrm{T}}\mathit{Q}{\displaystyle \sum _{i=\mathrm{1}}^r}{O}_{\mathrm{i}}\left(\mathit{v}\left(k\right)\right)A_{\mathrm{M}}X-\mathrm{2}{\mathit{S}}^{\mathrm{T}}\mathit{Q}{\mathit{Y}}_{\mathrm{d}}；

\mathit{S}={\displaystyle \sum _{i=\mathrm{1}}^r}{O}_i\left(\mathit{v}\left(k\right)\right){\mathit{B}}_{\mathrm{M}}；\mathit{X}为预测状态序列；{\mathit{Y}}_{\mathrm{d}}为期望输出序列.

上层MPC控制器获得横摆力矩以确保车辆稳定运行，而下层控制器的主要目的是对力矩进行合理分配.采用单侧分配的方法将横摆力矩分配给每个车轮，这样可以避免由于单个车轮制动压力过大，而出现抱死现象，具体力矩分配情况如下：

T_{\mathrm{l}\mathrm{f}}=\frac{\mathrm{4}{F}_{z\mathrm{l}\mathrm{f}}{M}_{z}R}{\left(F_{z\mathrm{l}\mathrm{f}}+F_{z\mathrm{l}\mathrm{r}}\right)B_{\mathrm{w}}},(25)

T_{\mathrm{l}\mathrm{r}}=\frac{\mathrm{4}{F}_{z\mathrm{l}\mathrm{r}}{M}_{z}R}{\left(F_{z\mathrm{l}\mathrm{f}}+F_{z\mathrm{l}\mathrm{r}}\right)B_{\mathrm{w}}},(26)

T_{\mathrm{r}\mathrm{f}}=\frac{\mathrm{4}{F}_{z\mathrm{r}\mathrm{f}}{M}_{z}R}{\left(F_{z\mathrm{r}\mathrm{f}}+F_{z\mathrm{r}\mathrm{r}}\right)B_{\mathrm{w}}},(27)

T_{\mathrm{r}\mathrm{r}}=\frac{\mathrm{4}{F}_{z\mathrm{r}\mathrm{r}}{M}_{z}R}{\left(F_{z\mathrm{r}\mathrm{f}}+F_{z\mathrm{r}\mathrm{r}}\right)B_{\mathrm{w}}}.(28)

2.3 MPC控制算法改进

由于MPC控制器中N_{\mathrm{p}}和N_{\mathrm{c}}的取值影响控制效果，当N_{\mathrm{p}}和N_{\mathrm{c}}取值较大时，控制器性能较好，但预测模型内部的参数矩阵维度随之增加，导致控制器计算量增大.因此引入Laguerre函数，在保证控制器性能的前提下，减少计算量，提高系统的运行速度.

Laguerre函数通过逼近控制输入\mathrm{\Delta }\mathit{u}来减少计算负担，其正交变换如下所示：

{\Gamma }_N\left(z_{\mathrm{a}}\right)=\frac{\sqrt[]{\mathrm{1}-p^{\mathrm{2}}}}{\mathrm{1}-pz{\mathrm{a}}^{-\mathrm{1}}}{\left(\frac{z{\mathrm{a}}^{-\mathrm{1}}-p}{\mathrm{1}-pz{\mathrm{a}}^{-\mathrm{1}}}\right)}^{N-\mathrm{1}}.(29)

式中：N是Laguerre函数的级数；p是Laguerre函数的极点，；z\mathrm{a}是正交变换的复频域变量.Laguerre函数的序列表示为向量形式：

\mathit{L}\left(k\right)={\left[{\mathit{l}}_{\mathrm{1}}\left(k\right),{\mathit{l}}_{\mathrm{2}}\left(k\right),\dots ,{\mathit{l}}_N\left(k\right)\right]}^{\mathrm{T}}.(30)

相邻时刻离散Laguerre函数之间的递推关系如下所示：

\mathit{L}\left(k+\mathrm{1}\right)={\mathit{A}}_{\mathrm{L}}\mathit{L}\left(k\right).(31)

式中：{\mathit{A}}_{\mathrm{L}}为Laguerre函数状态转移矩阵，

{\mathit{A}}_{\mathrm{L}}=\left[\begin{array}{cccc}p& \mathrm{0}& \cdots & \mathrm{0}\\ \theta & p& \cdots & \mathrm{0}\\ ⋮& ⋮& & ⋮\\ {\left(-\mathrm{1}\right)}^{N-\mathrm{2}}{p}^{N-\mathrm{2}}\theta & {\left(-\mathrm{1}\right)}^{N-\mathrm{3}}{p}^{N-\mathrm{3}}\theta & \cdots & p\end{array}\right]，\mathit{L}\left(\mathrm{0}\right)=\sqrt[]{\theta }{\left[\mathrm{1},-p,\dots ,{\left(-\mathrm{1}\right)}^{N-\mathrm{1}}{p}^{N-\mathrm{1}}\right]}^{\mathrm{T}}，\theta =\left(\mathrm{1}-p^{\mathrm{2}}\right).\theta为Laguerre函数极点p推导的归一化系数.

因此，k+z时刻的控制增量由Laguerre项表示如下：

\mathrm{\Delta }\mathit{u}(k+z)=\mathit{L}{\left(z\right)}^{\mathrm{T}}\mathit{\eta },(32)

\mathit{u}\left(k+z\right)=\mathit{u}\left(k-\mathrm{1}\right)+{\displaystyle \sum _{i=\mathrm{0}}^z}L{\left(i\right)}^{\mathrm{T}}\eta.(33)

式中：\mathit{L}\left(z\right)={\left[{\mathit{l}}_{\mathrm{1}}\left(z\right)，{\mathit{l}}_{\mathrm{2}}\left(z\right)，\dots ，{\mathit{l}}_N\left(z\right)\right]}^{\mathrm{T}}，\mathit{\eta }={\left[{\mathit{c}}_{\mathrm{1}},{\mathit{c}}_{\mathrm{2}},\cdots ,{\mathit{c}}_N\right]}^{\mathrm{T}}表示Laguerre函数的系数矩阵.

预测模型表示为

\begin{array}{c}\mathit{x}\left(k+z\left|k\right.\right)={\displaystyle \sum _{i=\mathrm{1}}^r}{O}_i\left(\mathit{v}\left(k\right)\right)\left[{\mathit{A}}_{\mathrm{M}}^z\mathit{x}\left(k\right)\right.+\\ \left.{\displaystyle \sum _{i=\mathrm{0}}^{z-\mathrm{1}}}{\mathit{A}}_{\mathrm{M}}^{z-i-\mathrm{1}}{\mathit{B}}_{\mathrm{M}}\mathit{L}{\left(i\right)}^{\mathrm{T}}\mathit{\eta }\right].\end{array}(34)

具有约束的目标函数可以重构为

\mathit{J}={{\displaystyle \sum _{z=\mathrm{1}}^{N_{\mathrm{p}}}}‖\mathit{y}\left(k+z\left|k\right.\right)-{\mathit{y}}_{\mathrm{d}}\left(k+z\right)‖}_{\mathit{Q}}^{\mathrm{2}}+{{\displaystyle \sum _{z=\mathrm{0}}^{N_{\mathrm{c}}-\mathrm{1}}}‖\mathit{L}{\left(z\right)}^{\mathrm{T}}\mathit{\eta }‖}_{\mathit{R}}^{\mathrm{2}}.(35)

用Laguerre函数的正交性简化式（35）的第二项，目标函数可以改写为

\mathit{J}={{\displaystyle \sum _{z=\mathrm{1}}^{N_{\mathrm{p}}}}‖\mathit{y}\left(k+z\left|k\right.\right)-{\mathit{y}}_{\mathrm{d}}\left(k+z\right)‖}_{\mathit{Q}}^{\mathrm{2}}+{\mathit{\eta }}^{\mathrm{T}}\mathit{R}\mathit{\eta }.(36)

将采用Laguerre函数逼近后的控制输入代入约束条件可得

\left.\begin{array}{c}\mathrm{\Delta }{\mathit{u}}_{\mathrm{m}\mathrm{i}\mathrm{n}}\left(k+z\right)\le \mathit{L}{\left(z\right)}^{\mathrm{T}}\mathit{\eta }\le \mathrm{\Delta }{\mathit{u}}_{\mathrm{m}\mathrm{a}\mathrm{x}}\left(k+z\right),\\ -{\mathit{u}}_{\mathrm{m}\mathrm{i}\mathrm{n}}\left(k+z\right)+\mathit{u}\left(k-\mathrm{1}\right)\le \mathit{u}\left(k-\mathrm{1}\right)+{\displaystyle \sum _{i=\mathrm{0}}^z}\mathit{L}{\left(i\right)}^{\mathrm{T}}\mathit{\eta }\le \\ -{\mathit{u}}_{\mathrm{m}\mathrm{a}\mathrm{x}}\left(k+z\right)+\mathit{u}\left(k-\mathrm{1}\right),\\ {\mathit{y}}_{\mathrm{m}\mathrm{i}\mathrm{n}}\left(k+z\right)\le \mathit{y}\left(k+z\right)\le {\mathit{y}}_{\mathrm{m}\mathrm{a}\mathrm{x}}\left(k+z\right).\end{array}\right\}(37)

将简化后的目标函数表示为二次型形式：

\left.\begin{array}{c}\mathrm{m}\mathrm{i}\mathrm{n}\mathit{J}=\frac{\mathrm{1}}{\mathrm{2}}{\mathit{\eta }}^{\mathrm{T}}{\mathit{H}}_{\mathrm{2}}\mathit{\eta }+{{\mathit{f}}_{\mathrm{2}}}^{\mathrm{T}}\mathit{\eta },\\ \mathrm{s}.\mathrm{t}.{\mathit{A}}_{\mathrm{c}}\mathit{\eta }\le {\mathit{b}}_{\mathrm{c}}.\end{array}\right\}(38)

式中：{\mathit{H}}_{\mathrm{2}}为二次项系数矩阵，{\mathit{H}}_{\mathrm{2}}=\mathrm{2}{\displaystyle \sum _{z=\mathrm{1}}^{N_{\mathrm{p}}}}\mathit{\sigma }\left(z\right)\mathit{Q}\mathit{\sigma }{\left(z\right)}^{\mathrm{T}}+\mathit{R}；{\mathit{f}}_{\mathrm{2}}为一次性系数向量，{\mathit{f}}_{\mathrm{2}}=\mathrm{2}{\displaystyle \sum _{z=\mathrm{1}}^{N_{\mathrm{p}}}}\sigma \left(z\right)Q\cdot \left({\displaystyle \sum _{i=\mathrm{1}}^r}{O}_i\left(\mathit{v}\left(k\right)\right){\mathit{A}}_{\mathrm{M}}^k\mathit{x}\left(k\right)-{\mathit{y}}_{\mathrm{d}}\left(k+z\right)\right)；\mathit{\sigma }\left(z\right)为权重系数矩阵，\mathit{\sigma }{\left(z\right)}^{\mathrm{T}}={\displaystyle \sum _{i=\mathrm{1}}^r}{O}_i\left(\mathit{v}\left(k\right)\right)\cdot{\displaystyle \sum _{i=\mathrm{1}}^{z-\mathrm{1}}}{\mathit{A}}_{\mathrm{M}}^{z-i-\mathrm{1}}{\mathit{B}}_{\mathrm{M}}\mathit{L}\left(i\right).

引入Laguerre函数后，利用级数N和极点p代替控制时域N_{\mathrm{c}}.当p=\mathrm{0}时，N=N_{\mathrm{c}}；当时，N明显小于N_{\mathrm{c}}，系统的矩阵维度减小.引入Laguerre函数可以在保证控制效果的同时，还能降低MPC控制器的计算量.然而，控制器的二次规划求解过程仍存在计算复杂度高的问题.为此采用原始对偶内点法处理约束条件，在保证数值稳定性的前提下快速迭代求解出维持车辆稳定所需的最优横摆力矩.

将松弛向量\mathit{\epsilon }引入式（38）中，将不等式约束转换为等式约束，可得原问题如下：

\left.\begin{array}{l}\mathrm{m}\mathrm{i}\mathrm{n}{\mathit{J}}_{\mathrm{P}}=\frac{\mathrm{1}}{\mathrm{2}}{\mathit{\eta }}^{\mathrm{T}}{\mathit{H}}_{\mathrm{2}}\mathit{\eta }+{{\mathit{f}}_{\mathrm{2}}}^{\mathrm{T}}\mathit{\eta },\\ \begin{array}{c}\mathrm{s}.\mathrm{t}.{\mathit{A}}_{\mathrm{c}}\mathit{\eta }+\mathit{\epsilon }={\mathit{b}}_{\mathrm{c}},\\ \mathit{\epsilon }\ge \mathrm{0}.\end{array}\end{array}\right\}(39)

式中，{\mathit{J}}_{\mathrm{P}}为原问题的目标函数.

对偶问题为

\left.\begin{array}{c}\mathrm{m}\mathrm{a}\mathrm{x}{\mathit{J}}_{\mathrm{D}}=-\frac{\mathrm{1}}{\mathrm{2}}{\mathit{\eta }}^{\mathrm{T}}{\mathit{H}}_{\mathrm{2}}\mathit{\eta }+{\mathit{\lambda }}^{\mathrm{T}}{\mathit{b}}_{\mathrm{c}},\\ \mathrm{s}.\mathrm{t}.\mathit{\lambda }\ge \mathrm{0}.\end{array}\right\}(40)

式中：{\mathit{J}}_{\mathrm{D}}为对偶问题的目标函数；\mathit{\lambda }为拉格朗日乘子.

扰动KKT条件表示为

\left.\begin{array}{c}{\mathit{H}}_{\mathrm{2}}\mathit{\eta }+{\mathit{f}}_{\mathrm{2}}-{\mathit{A}}_{\mathrm{c}}^{\mathrm{T}}\mathit{\lambda }=\mathrm{0},\mathit{\lambda }\ge \mathrm{0},\\ {\mathit{A}}_{\mathrm{c}}\mathit{\eta }+\mathit{\epsilon }={\mathit{b}}_{\mathrm{c}},\mathit{\epsilon }\ge \mathrm{0},\\ \mathit{M}\mathit{Y}\mathit{e}=\mathit{\tau }\mathit{e}.\end{array}\right\}(41)

式中：矩阵\mathit{M}为对偶间隙的权重矩阵；对角矩阵\mathit{Y}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left({\epsilon }_{\mathrm{1}},{\epsilon }_{\mathrm{2}},\cdots ,{\epsilon }_{m_{\mathrm{1}}}\right)；\mathit{e}\in {\mathbf{R}}^{m_{\mathrm{1}}}；m_{\mathrm{1}}为向量\mathit{e}的维度；\mathit{\tau }为对偶间隙参数；对偶间隔表示为

{\mathit{\delta }}_{\mathrm{P}\mathrm{D}}={\mathit{\mu }}^{\mathrm{T}}\mathit{\epsilon }=m_{\mathrm{1}}\mathit{\tau }.(42)

当\mathit{\tau }\to \mathrm{0}，{\mathit{\delta }}_{\mathrm{P}\mathrm{D}}\to \mathrm{0}，{\mathit{Z}}^{\left(k\right)}=\left[{\mathit{\eta }}^{\left(j\right)},{\mathit{\epsilon }}^{\left(j\right)},{\mathit{\lambda }}^{\left(j\right)}\right]为当前的迭代点；k,j为不同序列的迭代次数；{\mathit{\delta }}^{\left(k\right)}=\left[\mathrm{\Delta }{\mathit{\eta }}^{\left(j\right)},\mathrm{\Delta }{\mathit{\epsilon }}^{\left(j\right)},\mathrm{\Delta }{\mathit{\lambda }}^{\left(j\right)}\right]为当前的搜索方向；\mathrm{\Delta }\mathit{\eta }，\mathrm{\Delta }\mathit{\epsilon }，\mathrm{\Delta }\mathit{\lambda }分别为\mathit{\eta }，\mathit{\epsilon }，\mathit{\lambda }的迭代增量.

将{\mathit{Z}}^{\left(j+\mathrm{1}\right)}代入扰动KKT条件中，通过牛顿迭代法求解搜索方向：

\left.\begin{array}{c}{\mathit{H}}_{\mathrm{2}}\left({\mathit{\eta }}^{\left(j\right)}+\mathrm{\Delta }{\mathit{\eta }}^{\left(j\right)}\right)+{\mathit{f}}_{\mathrm{2}}-{\mathit{A}}_{\mathrm{c}}^{\mathrm{T}}\left({\mathit{\lambda }}^{\left(j\right)}+\mathrm{\Delta }{\mathit{\lambda }}^{\left(j\right)}\right)=\mathrm{0},\\ {\mathit{M}}^{\left(j\right)}{\mathit{Y}}^{\left(j\right)}\mathit{e}+{\mathit{M}}^{\left(j\right)}\mathrm{\Delta }{\mathit{\epsilon }}^{\left(j\right)}+{\mathit{Y}}^{\left(j\right)}\mathrm{\Delta }{\mathit{\lambda }}^{\left(j\right)}={\mathit{\tau }}^{\left(j+\mathrm{1}\right)}\mathit{e},\\ {\mathit{A}}_{\mathrm{c}}\left({\mathit{\eta }}^{\left(j\right)}+\mathrm{\Delta }{\mathit{\eta }}^{\left(j\right)}\right)+\left({\mathit{\epsilon }}^{\left(j\right)}+\mathrm{\Delta }{\mathit{\epsilon }}^{\left(j\right)}\right)={\mathit{b}}_{\mathrm{c}}.\end{array}\right\}(43)

将式（43）写成矩阵形式：

\left[\begin{array}{ccc}{\mathit{H}}_{\mathrm{2}}& \mathrm{0}& -{\mathit{A}}_{\mathrm{c}}^{\mathrm{T}}\\ \mathrm{0}& {\mathit{M}}^{\left(j\right)}& {\mathit{Y}}^{\left(j\right)}\\ {\mathit{A}}_{\mathrm{c}}& \mathit{I}& \mathrm{0}\end{array}\right]\left[\begin{array}{c}\mathrm{\Delta }{\mathit{\eta }}^{\left(j\right)}\\ \mathrm{\Delta }{\mathit{\epsilon }}^{\left(j\right)}\\ \mathrm{\Delta }{\mathit{\lambda }}^{\left(j\right)}\end{array}\right]=\left[\begin{array}{c}{\mathit{A}}_{\mathrm{c}}^{\mathrm{T}}{\mathit{\lambda }}^{\left(j\right)}-{\mathit{f}}_{\mathrm{2}}-{\mathit{H}}_{\mathrm{2}}{\mathit{\eta }}^{\left(j\right)}\\ {\mathit{\tau }}^{\left(j+\mathrm{1}\right)}\mathit{e}-{\mathit{M}}^{\left(j\right)}{\mathit{Y}}^{\left(j\right)}\mathit{e}\\ {\mathit{b}}_{\mathrm{c}}-{\mathit{\epsilon }}^{\left(j\right)}-{\mathit{A}}_{\mathrm{c}}{\mathit{\eta }}^{\left(j\right)}\end{array}\right].(44)

对式（44）进行变换，得到如下形式：

\left[\begin{array}{c}{\mathit{H}}_{\mathrm{2}}{\mathit{\eta }}^{\left(j\right)}+{\mathit{f}}_{\mathrm{2}}-{\mathit{A}}_{\mathit{c}}^{\mathrm{T}}{\mathit{\lambda }}^{\left(j\right)}\\ {\left({\mathit{Y}}^{\left(j\right)}\right)}^{-\mathrm{1}}{\mathit{\tau }}^{\left(j+\mathrm{1}\right)}\mathit{e}-{\left({\mathit{Y}}^{\left(j\right)}\right)}^{-\mathrm{1}}{\mathit{M}}^{\left(j\right)}{\mathit{Y}}^{\left(j\right)}\mathit{e}\\ {\mathit{b}}_{\mathrm{c}}-{\mathit{\epsilon }}^{\left(j\right)}-{\mathit{A}}_{\mathrm{c}}{\mathit{\eta }}^{\left(j\right)}\end{array}\right]=\left[\begin{array}{c}{\mathit{v}}_{\eta }^{\left(j\right)}\\ -{\mathit{v}}_{\epsilon }^{\left(j\right)}\\ {\mathit{v}}_{\lambda }^{\left(j\right)}\end{array}\right].(45)

式中，{\mathit{v}}_{\eta }^{\left(j\right)}，{\mathit{v}}_{\epsilon }^{\left(j\right)}，{\mathit{v}}_{\lambda }^{\left(j\right)}分别为左边矩阵元素的等效代替.

解决以下线性规划问题：

\mathrm{m}\mathrm{i}\mathrm{n}\frac{\mathrm{1}}{\mathrm{2}}{\left[\mathrm{\Delta }{\mathit{\eta }}^{\left(j\right)}\right]}^{\mathrm{T}}{{\mathit{H}}_{\mathrm{2}}}^{\left(j\right)}\mathrm{\Delta }{\mathit{\eta }}^{\left(j\right)}+{\left(\mathrm{\Delta }{\mathit{\eta }}^{\left(j\right)}\right)}^{\mathrm{T}}{{\mathit{f}}_{\mathrm{2}}}^{\left(j\right)}.(46)

式中：{{\mathit{H}}_{\mathrm{2}}}^{\left(j\right)}={\mathit{H}}_{\mathrm{2}}+{\mathit{A}}_{\mathrm{c}}^{\mathrm{T}}{\left({\mathit{Y}}^{\left(j\right)}\right)}^{-\mathrm{1}}{\mathit{M}}^{\left(j\right)}{\mathit{A}}_{\mathrm{c}}；{{\mathit{f}}_{\mathrm{2}}}^{\left(j\right)}={\mathit{v}}_{\eta }^{\left(j\right)}-

{\mathit{A}}_{\mathrm{c}}^{\mathrm{T}}\left[{\mathit{v}}_{\epsilon }^{\left(j\right)}+{\left({\mathit{Y}}^{\left(j\right)}\right)}^{-\mathrm{1}}{\mathit{M}}^{\left(j\right)}{\mathit{v}}_{\lambda }^{\left(j\right)}\right].

求得搜索方向为

\left.\begin{array}{c}\mathrm{\Delta }{\mathit{\eta }}^{\left(j\right)}=-{\left({{\mathit{H}}_{\mathrm{2}}}^{\left(j\right)}\right)}^{-\mathrm{1}}{{\mathit{f}}_{\mathrm{2}}}^{\left(j\right)},\\ \mathrm{\Delta }{\mathit{\epsilon }}^{\left(j\right)}=-{\mathit{v}}_{\lambda }^{\left(j\right)}-{\mathit{A}}_{\mathrm{c}}\mathrm{\Delta }{\mathit{\eta }}^{\left(j\right)},\\ \mathrm{\Delta }{\mathit{\lambda }}^{\left(j\right)}={\mathit{v}}_{\epsilon }^{\left(j\right)}-{\left({\mathit{Y}}^{\left(j\right)}\right)}^{-\mathrm{1}}{\mathit{M}}^{\left(j\right)}\mathrm{\Delta }{\mathit{\epsilon }}^{\left(j\right)}.\end{array}\right\}(47)

计算搜索步长{\mathit{\alpha }}^{\left(j\right)}为

(48)

式中：{\mathit{\alpha }}_{\mathrm{P},\mathrm{m}\mathrm{i}\mathrm{n}}^{\left(j\right)}为最小原始步长；{\mathit{\alpha }}_{\mathrm{D},\mathrm{m}\mathrm{i}\mathrm{n}}^{\left(j\right)}为最小对偶步长.

更新迭代点：

{\mathit{Z}}^{\left(j+\mathrm{1}\right)}={\mathit{Z}}^{\left(j\right)}+{\left[\mathrm{\Delta }{\mathit{\eta }}^{\left(j\right)}，\mathrm{\Delta }{\mathit{\epsilon }}^{\left(j\right)}，\mathrm{\Delta }{\mathit{\lambda }}^{\left(j\right)}\right]}^{\mathrm{T}}{\mathit{\alpha }}^{\left(j\right)}.(49)

当对偶间隔{\mathit{\delta }}_{\mathrm{P}\mathrm{D}}^{\left(j+\mathrm{1}\right)}满足精度要求时，迭代过程停止，此时原始对偶内点法的解为{\mathit{Z}}^{\mathrm{*}}=\left[{\mathit{\eta }}^{\mathrm{*}},{\mathit{\epsilon }}^{\mathrm{*}},{\mathit{\lambda }}^{\mathrm{*}}\right]，其中{\mathit{\eta }}^{\mathrm{*}}，{\mathit{\epsilon }}^{\mathrm{*}}，{\mathit{\lambda }}^{\mathrm{*}}分别为\mathit{\eta }，\mathit{\epsilon }，\mathit{\lambda }的最优解.将{\mathit{\eta }}^{\mathrm{*}}代入式（33），可计算维持车辆稳定行驶所需的横摆力矩.

3 仿真验证

为了验证本文所提的快速MPC控制器的实际性能与有效性，采用MATLAB/Simulink和TruckSim联合仿真进行验证.选择双移线和鱼钩两种典型工况，分别针对传统MPC、基于Laguerre函数的MPC（LMPC）以及原始对偶内点法改进的LMPC（PDIP-LMPC）进行对比分析.仿真中所使用的车辆参数如表1所示^［14］：

3.1 双移线工况仿真

车辆的初始速度设置为80 km/h，路面附着系数设置为0.85，侧翻阈值设置为0.8.在双移线工况下，LTR指标、横摆角速度和侧翻角的仿真结果分别如图3~图5所示.

由图3可知，在双移线工况下，LTR指标始终处于设定的阈值内，并且最终稳定于0，车辆未发生侧翻，3种控制策略均能维持车辆稳定.由图4和图5可知，当车辆转角发生变化时，均为PDIP-LMPC先响应，表明引入Laguerre函数和PDIP算法后提高了控制器的响应速度，控制器能快速对车辆进行控制，且车辆的各个性能指标幅值均有小幅度减小，车辆稳定性得到进一步提高.

3.2 鱼钩转向工况仿真

在鱼钩转向工况下，车辆的转角如图6所示.

车辆的初始速度和路面附着系数与双移线工况相同.LTR指标、横摆角速度和侧翻角的仿真结果分别如图7~图9所示.

由图7可知，在鱼钩转向工况下，LTR指标始终处于可控范围内，说明车辆未发生侧翻.由图8和图9可知，当车辆转角变化时，PDIP-LMPC表现出最快的响应速度.结果表明，在引入Laguerre函数和PDIP算法对传统MPC进行改进后，控制器的响应速度加快，能够快速控制车辆的侧翻稳定性.同时PDIP-LMPC控制器使车辆的横摆角速度和侧翻角的峰值减小，并且在2 s以后，各性能指标的平滑性也变得更好，车辆稳定性更优.

4 结 语

本文针对重型车辆易发生侧翻的问题，建立T-S模糊纵向速度时变线性化模型，设计基于MPC的防侧翻控制策略，通过引入Laguerre函数和PDIP算法减轻了传统MPC控制器的计算负担.提出的分层控制策略能够快速高效地维持车辆侧翻稳定性.通过MATLAB/Simulink与TruckSim联合仿真验证，结果表明，在极端转向工况下，提出的分层控制策略能够有效维持车辆的侧翻稳定性；相较于传统MPC和LMPC，PDIP-LMPC具有最快的响应速度，且各性能指标的平滑性均得到改善.该控制策略显著降低了传统MPC控制器的计算量，使控制器能够快速高效地防止车辆侧翻.

参考文献

[1] Yang X W， Wu C Z， He Y， et al. A dynamic rollover prediction index of heavy-duty vehicles with a real-time parameter estimation algorithm using NLMS method［J］. IEEE Transactions on Vehicular Technology， 2022， 71（3）： 2734-2748.

[2] Jin Z L， Zhang L， Zhang J L， et al. Stability and optimised H_∞ control of tripped and untripped vehicle rollover［J］. Vehicle System Dynamics， 2016，54（10）： 1405-1427.

[3] Song Z S， Wang L L， Liu Y M， et al. Actively steering a wheeled tractor against potential rollover using a sliding-mode control algorithm： scaled physical test［J］. Biosystems Engineering， 2022， 213： 13-29.

[4] Hu C， Wang Z F， Qin Y C， et al. Lane keeping control of autonomous vehicles with prescribed performance considering the rollover prevention and input saturation［J］. IEEE Transactions on Intelligent Transportation Systems， 2020， 21（7）： 3091-3103.

[5] Zhang Y， Lu J W， Xia G， et al. Enhancing roll stability of counterbalance forklift trucks with a novel rollover index and hierarchical adaptive model predictive control［J］. IEEE Transactions on Vehicular Technology， 2024， 73（6）： 7925-7938.

[6] Wang L P. Discrete model predictive controller design using Laguerre functions［J］. Journal of Process Control， 2004， 14（2）： 131-142.

[7] Wang W R， Yan J H， Wang H， et al. Adaptive MPC trajectory tracking for AUV based on Laguerre function［J］. Ocean Engineering， 2022， 261： 111870.

[8] Jeong D， Choi S B. Tracking control based on model predictive control using laguerre functions with pole optimization［J］. IEEE Transactions on Intelligent Transportation Systems， 2022， 23（11）： 20652-20663.

[9] Panahi S， Kashani A， Danielson C. Primal-dual interior-point algorithm for symmetric model predictive control［J］. Automatica， 2023， 155： 111157.

[10] Hu W J， Zhou Y H， Zhang Z L， et al. Model predictive control for hybrid levitation systems of maglev trains with state constraints［J］. IEEE Transactions on Vehicular Technology， 2021， 70（10）：9972-9985.

[11] Chu D F， Li Z L， Wang J M， et al. Rollover speed prediction on curves for heavy vehicles using mobile smartphone［J］. Measurement， 2018， 130：404-411.

[12] Larish C， Piyabongkarn D， Tsourapas V， et al. A new predictive lateral load transfer ratio for rollover prevention systems［J］. IEEE Transactions on Vehicular Technology， 2013， 62（7）：2928-2936.

[13] 李寿涛，魏玉博，李秋媛，等.考虑车辆侧偏刚度变化的MPC稳定性控制方法［J］.东北大学学报（自然科学版）， 2023， 44（2）： 162-167.

[14] Li Shou-tao， Wei Yu-bo， Li Qiu-yuan， et al. MPC stability control method considering the variation of vehicle’s cornering stiffness［J］. Journal of Northeastern University （Natural Science）， 2023， 44（2）： 162-167.

[15] 朱天军，麻威，王振峰，等.基于AdaBoost算法的重型车辆侧翻预警研究［J］.重庆交通大学学报（自然科学版）， 2021， 40（8）： 25-33.

[16] Zhu Tian-jun， Ma Wei， Wang Zhen-feng， et al. Rollover warning study of heavy vehicle based on AdaBoost algorithm［J］. Journal of Chongqing Jiaotong University （Natural Science）， 2021， 40（8）： 25-33.

基金资助

中央高校基本科研业务费专项资金资助项目(N2223029)