随着移动机器人技术的突破性进展^[1-2]，类车机器人凭借突出的载重能力与高速运动特性^[3-4]，在工业物流、智能制造、现代农业等领域获得了广泛应用. 虽然类车机器人的转向形式和系统布局与无人驾驶车辆相似，但其主要面向低速工况，无须考虑高速或低附着工况下的高动态问题，故而常采用非标轮胎、低功率转向电机、简化差动机构等低成本部件^[5-7]. 因此，直接将面向无人驾驶车辆的路径跟踪控制方法应用于类车机器人，往往会出现难以应对曲率突变、控制精度不足的问题. 在此背景下，专门针对类车机器人路径跟踪控制的研究工作日益增多^[5-6,8-17]. 在众多控制方法中，模型预测控制（model predictive control, MPC）因具有能够显式处理系统约束的优势^[18]，成为当前应用前景较好的一种技术方案^[19].

现有MPC控制方法在工程应用中面临两难困境. 非线性模型预测控制（nonlinear MPC, NMPC）具有较高的控制精度，但其计算复杂度较高，导致实时性不足. 尽管通过优化求解算法^[20-22]、优化时域参数^[23]以及降低搜索维度^[24]等手段可以提升计算效率，但是NMPC所需的计算资源仍然较多，导致该控制方法在应用于类车机器人时面临着较高的成本压力. 线性模型预测控制（linear MPC, LMPC）可以确保实时性，在前轮转向车辆^[25-28]、四轮转向车辆^[29-30]、拖挂式农用机械^[31]的路径跟踪控制中得到了广泛应用. 但是，这种控制方法一般采用单参考点（single reference point, SRP）控制架构，仅能依据当前位姿偏差进行预测，无法有效利用前方参考路径信息. 特别是在跟踪曲率突变路径时，SRP架构由于缺乏前瞻性路径信息，导致路径跟踪控制系统无法预判即将到来的路径变化^[32]. 当类车机器人以较高速度通过U型弯之类的弯道时，这种响应特性会引发显著的位移误差累积，甚至造成系统失控. 因此，虽然SRP-LMPC在平直路径场景下表现良好，但在仓库、车间、农田等存在大幅曲率突变参考路径的类车机器人典型工况中，控制精度明显不足.

为了同时提高类车机器人路径跟踪控制的精度和实时性，首先相对SRP-LMPC重新设定线性展开点，并对差分模型进行修正，以更准确地预测系统状态. 随后，采用非线性与线性航迹推算方法获得非线性补偿量，并据此设计多参考点的优化目标函数，构建出新型的路径跟踪控制方法，即基于非线性补偿的多参考点线性模型预测控制（multiple reference points LMPC, MRP-LMPC）. 通过MATLAB与Carsim联合仿真、硬件在环仿真进行测试验证，还与NMPC、SRP-LMPC、PID（proportional integral derivative）等其他控制方法进行比较研究.

目前基于LMPC的路径跟踪控制方法通常采用差分线性模型^[25-26,31]或误差模型^[27-30]作为预测模型，两者虽然最终形式有所不同，理论基础也有所差异，但是控制原理基本相同. 如图1所示，这些控制方法的整体控制框架均是首先读取被控对象当前位置与参考路径之间的误差，如位移误差和航向误差，然后利用差分线性模型或误差模型，基于当前误差值预测未来的误差值，从而实现误差值在未来一定时域内的收敛. 这种控制方法在参考路径曲率平滑且变化速度较慢时能够实现较精确的路径跟踪控制，但是由于未能引入被控对象前方的参考路径信息，在参考路径曲率大幅突变时容易出现被控对象无法及时响应，导致误差大幅增加的问题.

多参考点模型预测控制可以通过预测模型预测类车机器人未来的位置、航向状态，如图2所示. 通过使这些状态尽可能地靠近预测时域内的参考路径点，可以实现对具有大幅曲率突变的参考路径的精确跟踪. 不过，目前多参考点模型预测控制方法多为NMPC方法，其非线性优化求解过程计算成本较高，因此实时性难以得到保障. 考虑到LMPC可以将优化目标转化为标准二次型，其计算成本远低于NMPC，因此本研究拟解决的问题可以归结为提出能够引入多个参考点的LMPC控制方法的问题.

所提出的多参考点线性MPC控制器总体架构如图3所示.

类车机器人一般使用非标准轮胎，这使得其轮胎侧偏刚度参数难以精确测量. 此类机器人行驶速度较低，通常不会达到轮胎的非线性响应区域，因此可以采用运动学模型来设计控制器. 此外，考虑到左右侧的运动特性相同，模型可以进一步简化为如图4所示的单轨模型. 其中，$ \delta $表示前轮转角；$ \theta $表示航向角；$ v $为车辆速度；$ P $为控制点，为后轴中心点，其横纵坐标分别为$ x $和$ y $；$ XOY $为全局坐标系.

该模型可以通过以下数学公式进行描述：

(1)$ \left[\begin{array}{c}\dot{x}\\\dot{y}\\\dot{\theta }\end{array}\right]=\left[\begin{array}{c}v\cos\; \theta \\v\sin\; \theta \\\left(v\tan\; \delta \right)/l\end{array}\right]. $

式中：$ l $为轴距.

在路径跟踪控制过程中常假设纵向速度为恒定值，因此模型可进一步简化为向量表达形式：

(2)$ \left.\begin{split} & \dot{\boldsymbol{\xi}}=f(\boldsymbol{\xi}, \boldsymbol{\kappa}), \\& \boldsymbol{\xi}=\left[x, \;\; y, \;\; \theta\right]^{\mathrm{T}}, \\&\boldsymbol{\kappa}=[\delta].\end{split}\right\} $

式中：$ \boldsymbol{\xi } $和$ \boldsymbol{\kappa } $为状态向量和控制向量，$ f $表示映射.

在运动学模型的线性化处理过程中，由于线性化可在任意位置进行展开，可以将线性化点选定为该步展开时类车机器人的当前所在位置，然后将模型线性化展开为

(3)$ \left.\begin{split} &\boldsymbol{\dot{\tilde{\xi}}}=\boldsymbol{A} \tilde{\boldsymbol{\xi}}+\boldsymbol{B} \tilde{\boldsymbol{\kappa}}, \\&\tilde{\boldsymbol{\xi}}=\boldsymbol{\xi}-\boldsymbol{\xi}_0, \\&\tilde{\boldsymbol{\kappa}}=\boldsymbol{\kappa}-\boldsymbol{\kappa}_0, \\&\boldsymbol{A}=\left[\begin{array}{ccc}0 & 0 & -v \sin \;\theta \\0 & 0 & v \cos \;\theta \\0 & 0 & 0\end{array}\right] ,\\&\boldsymbol{B}=\left[0 ,\;\; 0 ,\;\; 1\right]^{\mathrm{T}}.\end{split}\right\}$

式中：下标0表示该步展开时类车机器人的状态，上标~表示差分量，$ \boldsymbol{A} $和$ \boldsymbol{B} $表示状态系数矩阵和输入系数矩阵.

由于控制周期较短，可以假设线性化后的差分量在控制周期内呈线性变化：

(4)$ \mathbf{\dot{\tilde{\boldsymbol{\xi }} }}\left(t+1\right)T\approx \mathbf{\tilde{\boldsymbol{\xi }}}\left(t+1\right)-\mathbf{\tilde{\boldsymbol{\xi }}}\left(t\right). $

式中：$ t $表示$ t $时刻，$ t+1 $表示$ t+1 $时刻，$ T $为控制周期.

联立式(3)、(4)，可以得到如下差分模型：

(5)$ \left.\begin{split} & \tilde{\boldsymbol{\xi}}(t+1)=\tilde{\boldsymbol{A}} \tilde{\boldsymbol{\xi}}(t)+\tilde{\boldsymbol{B}} \tilde{\boldsymbol{\kappa}}(t+1), \\ & \tilde{\boldsymbol{A}}=\boldsymbol{I}+T \boldsymbol{A}, \\&\tilde{\boldsymbol{B}}=T \boldsymbol{B}.\end{split}\right\} $

式中：$ {\tilde{\boldsymbol{A}}} $和${\tilde{\boldsymbol{B}}} $表示差分状态系数矩阵和差分输入系数矩阵，$ \boldsymbol{I} $表示单位矩阵.

差分模型仅能通过当前偏差预测未来偏差，因此须提取状态量：

(6)$ \begin{split}\boldsymbol{\xi }\left(t+1\right)=&{\tilde{\boldsymbol{A}}\boldsymbol{\xi }}\left(t\right)+{\tilde{\boldsymbol{B}}\boldsymbol{\kappa }}\left(t+1\right)+{\boldsymbol{\xi }}_{\text{0}}\left(t+1\right)-\\&{\tilde{\boldsymbol{A}}}{\boldsymbol{\xi }}_{\text{0}}\left(t\right)-{\tilde{\boldsymbol{B}}}{\boldsymbol{\kappa }}_{0}\left(t+1\right).\end{split} $

为了构建便于考虑平顺性的差分输入，状态量模型可以扩增为

(7)$ \left.\begin{split} & \hat{\boldsymbol{\xi}}(t+1)=\hat{\boldsymbol{A}} \hat{\boldsymbol{\xi}}(t)+\hat{\boldsymbol{B}} \tilde{\boldsymbol{\kappa}}(t+1)+ \\& \qquad\hat{\boldsymbol{\xi}_0}(t+1)-\hat{\boldsymbol{A}} \hat{\boldsymbol{\xi}}_0(t)-\hat{\boldsymbol{B}} \tilde{\boldsymbol{\kappa}}_0(t+1), \\&\hat{\boldsymbol{\xi}}=\left[\boldsymbol{\xi}^{\mathrm{T}} ,\;\; \boldsymbol{\kappa}\right]^{\mathrm{T}}, \\&\hat{\boldsymbol{A}}=\left[\begin{array}{cc}\hat{\boldsymbol{A}} & \tilde{\boldsymbol{B}} \\{\left[0,\;\; 0 ,\;\;0\right]} & \boldsymbol{I}\end{array}\right], \\&\hat{\boldsymbol{B}}=\left[\tilde{\boldsymbol{B}}^{\mathrm{T}} ,\;\; \boldsymbol{I}\right]^{\mathrm{T}}.\end{split}\right\} $

式中：上标^表示扩增量.

根据式(7)可知，类车机器人在下一时刻的状态由3个部分组成. 第1部分是当前时刻状态量的线性迭代结果，第2部分是下一时刻输入量的线性迭代结果，第3部分是状态量非线性迭代与状态量线性迭代之间的差值.

将迭代差值定义为非线性补偿量：

(8)$ \boldsymbol{\gamma }\left(t+1\right)={{\hat{\boldsymbol{\xi }}}}_{\text{0}}\left(t+1\right)-{\hat{\boldsymbol{A}}}{{\hat{\boldsymbol{\xi }}}}_{\text{0}}\left(t\right)-{\hat{\boldsymbol{B}}}{\mathbf{\tilde{\boldsymbol{\kappa }}}}_{0}\left(t+1\right). $

式中：$ \boldsymbol{\gamma } $表示非线性补偿量.

然后，即可获得线性化的状态量模型：

(9)$ {\hat{\boldsymbol{\xi }}}\left(t+1\right)={\hat{\boldsymbol{A}}\hat{\boldsymbol{\xi }}}\left(t\right)+{\hat{\boldsymbol{B}}\tilde{\boldsymbol{\kappa }}}\left(t+1\right)+\boldsymbol{\gamma }\left(t+1\right). $

线性化的状态量模型在每一次迭代前均须先获取非线性补偿量，而非线性补偿量是状态量非线性迭代预测值和状态量线性迭代预测值之间的差值，因此首先应计算状态量的非线性迭代预测值.

由于控制周期为较小时间值，可以采用欧拉法获取状态量的非线性迭代预测值：

(10)$ {\left.\begin{split}\hat{\boldsymbol{\xi}}_0(t+1)= & \;\hat{\boldsymbol{\xi}}_0(t)+T{f}\left(\hat{\boldsymbol{\xi}}_0(t), \tilde{\boldsymbol{\kappa}}_0(t+1)\right) ,\\[-3pt]& \vdots \\[-3pt]\hat{\boldsymbol{\xi}}_0(t + c)= & \;\hat{\boldsymbol{\xi}}_0(t + c - 1) + T{f}\left(\hat{\boldsymbol{\xi}}_0(t + c - 1), \tilde{\boldsymbol{\kappa}}_0(t + c)\right), \\[-3pt]& \vdots \\[-3pt]\hat{\boldsymbol{\xi}}_0(t + p)= &\; \hat{\boldsymbol{\xi}}_0(t + p - 1) + T{f}\left(\hat{\boldsymbol{\xi}}_0(t + p - 1), \tilde{\boldsymbol{\kappa}}_0(t + c)\right).\end{split}\right\}} $

式中：$ c $为控制步数，$ p $为预测步数.

将式(10)代入式(8)，即可以获得每一步展开时的非线性补偿量：

(11)$ {\left.\begin{split}\boldsymbol{\gamma}(t+1)= &\; \hat{\boldsymbol{\xi}}_0(t+1)- \hat{\boldsymbol{A}} \hat{\boldsymbol{\xi}}_0(t)-\hat{\boldsymbol{B}} \tilde{\boldsymbol{\kappa}}_0(t+1), \\[-3pt]& \vdots \\[-3pt]\boldsymbol{\gamma}(t+c)= & \;\hat{\boldsymbol{\xi}}_0(t+c)- \hat{\boldsymbol{A}} \hat{\boldsymbol{\xi}}_0(t+c-1)-\hat{\boldsymbol{B}} \tilde{\boldsymbol{\kappa}}_0(t+c), \\[-3pt]& \vdots \\[-3pt]\boldsymbol{\gamma}(t+p)= & \;\hat{\boldsymbol{\xi}}_0(t+p)- \hat{\boldsymbol{A}} \hat{\boldsymbol{\xi}}_0(t+p-1)-\hat{\boldsymbol{B}} \tilde{\boldsymbol{\kappa}}_0(t+c).\end{split}\right\}} $

理论上，当每一步迭代时的输入与实际输入完全相同时，获得的非线性补偿量最为准确. 然而，实际输入是一个未知数，需要控制器求解后才能获取. 不过，根据式(11)可知，输出的补偿量是非线性迭代结果和线性迭代结果的差分量. 为了保障实时性，控制步数常应取较小值^[23]，且控制周期普遍较短，根据式(5)可以明显看出，线性迭代过程中状态量起主导作用，其中控制增量的系数即为迭代周期. 另外，式(11)输出的是差分量，微小数值经差分运算后数值进一步减小，故而完全可以忽略不计. 因此，可以假设每一步迭代时的输入均等于前一个控制周期的实际输入，即各时刻的控制增量均为零：

(12)$ {{\tilde{\boldsymbol{\kappa }}}}_{{0}}\left(t+i\right)={\bf{0}};\; i=1,2,\cdots ,c. $

根据式(9)的迭代形式，可得预测时域内的所有预测状态量：

(13)$ \left.\begin{split}\hat{\boldsymbol{\xi}}(t+1)= & \hat{\boldsymbol{A}} \hat{\boldsymbol{\xi}}(t)+\hat{\boldsymbol{B}} \tilde{\kappa}(t+1)+{\boldsymbol{\gamma}}(t+1), \\& \vdots \\\hat{\boldsymbol{\xi}}(t+c)= & \hat{\boldsymbol{A}} \hat{\boldsymbol{\xi}}(t+c-1)+ \hat{\boldsymbol{B}} \tilde{\boldsymbol{\kappa}}(t+c)+{\boldsymbol{\gamma}}(t+c) ,\\& \vdots \\\hat{\boldsymbol{\xi}}(t+p)= & \hat{\boldsymbol{A}} \hat{\boldsymbol{\xi}}(t+p-1)+ \hat{\boldsymbol{B}} \tilde{{\boldsymbol { k }}}(t+c)+{\boldsymbol{\gamma}}(t+p).\end{split}\right\} $

为了便于后续构建标准二次型形式的优化目标函数，可将式(13)改为矩阵形式：

(14)$ \left.\begin{split} &\boldsymbol{Y}(t)=\boldsymbol{\varPsi} \hat{\xi}(t)+\boldsymbol{\varTheta} \boldsymbol{K}(t)+\boldsymbol{\varOmega} \boldsymbol{\varGamma}(t),\\& \boldsymbol{Y}(t)=\left[\hat{{\boldsymbol{\xi}}}(t+1),\; \cdots ,\; \hat{{\boldsymbol{\xi}}}(t+p)\right]^{\mathrm{T}} ,\\&\boldsymbol{\varPsi}=\left[\hat{\boldsymbol{A}} ,\; \hat{\boldsymbol{A}}^2 ,\; \cdots ,\; \hat{\boldsymbol{A}}^p\right]^{\mathrm{T}}, \\&\boldsymbol{\varTheta}=\left[\begin{array}{cccc}\hat{\boldsymbol{B}} &{\boldsymbol{ 0}} & \cdots & {\bf{0}} \\\hat{\boldsymbol{A}} \hat{\boldsymbol{B}} & \hat{\boldsymbol{B}} & \cdots & {\bf{0}} \\\vdots & \vdots & \; & \vdots \\\hat{\boldsymbol{A}}^{c-1} \hat{\boldsymbol{B}} & \hat{\boldsymbol{A}}^{c-2} \hat{\boldsymbol{B}} & \cdots & \hat{\boldsymbol{B}} \\\hat{\boldsymbol{A}}^c \hat{\boldsymbol{B}} & \hat{\boldsymbol{A}}^{c-1} \hat{\boldsymbol{B}} & \cdots & \hat{\boldsymbol{A}} \hat{\boldsymbol{B}} \\\vdots & \vdots & \; & \vdots \\\hat{\boldsymbol{A}}^{p-1} \hat{\boldsymbol{B}} & \hat{\boldsymbol{A}}^{p-2} \hat{\boldsymbol{B}} & \cdots & \hat{\boldsymbol{A}}^{p-c} \hat{\boldsymbol{B}}\end{array}\right] , \\&\boldsymbol{K}(t)=\left[\tilde{\boldsymbol{\kappa}}(t+1) ,\; \cdots ,\; \tilde{\boldsymbol{\kappa}}(t+c)\right]^{\mathrm{T}}, \\&\boldsymbol{\varOmega}=\left[\begin{array}{cccc}\boldsymbol{I} & {\bf{0}} & \cdots & {\bf{0}} \\\hat{\boldsymbol{A}} & \boldsymbol{I} & \cdots & {\bf{0}} \\\vdots & \vdots & \; & \vdots \\\hat{\boldsymbol{A}}^{p-1} & \hat{\boldsymbol{A}}^{p-2} & \cdots & \boldsymbol{I}\end{array}\right] ,\\&\boldsymbol{\varGamma}(t)=\left[\boldsymbol{\gamma}(t+1) ,\; \cdots ,\; \boldsymbol{\gamma}(t+p)\right]^{\mathrm{T}}.\end{split}\right\} $

式中：$ \boldsymbol{Y} $表示预测状态量矩阵，$ \boldsymbol{\varPsi } $表示扩增差分状态系数矩阵，$ \boldsymbol{\varTheta } $表示扩增差分输入系数矩阵，$ \boldsymbol{K} $表示输入量矩阵，$ \boldsymbol{\varOmega } $表示非线性补偿系数矩阵，$ \boldsymbol{\varGamma } $表示非线性补偿量矩阵.

由于预测模型（式(14)）可以给出未来$ p $个时刻的状态，在构造优化目标函数时必须先获取$ p $个对应的参考路径点信息，构建如下参考路径信息：

(15)$ {\boldsymbol{Y}}_{\text{ref}}\left(t\right)\text=\left[\begin{array}{c}{{\hat{\boldsymbol{\xi }}}}_{\text{ref}}\left(t|1\right)\\{{\hat{\boldsymbol{\xi }}}}_{\text{ref}}\left(t|2\right)\\\vdots \\{{\hat{\boldsymbol{\xi }}}}_{\text{ref}}\left(t|p\right)\end{array}\right]. $

式中：${{\boldsymbol{\hat{\xi }}}_{\text{ref}}}\left( t|p \right) $表示t时刻的第p个参考路径点信息.

取参考路径上距离类车机器人当前位置最近的点作为基准点，$ t $时刻的第1个参考路径点应位于类车机器人当前位置前方，与参考路径基准点之间的距离应等于路径跟踪控制系统的期望速度与控制周期之积. 此外，各参考路径点之间的距离也应等于路径跟踪控制系统的期望行驶速度与控制周期之积. 在此基础上，为了兼顾路径跟踪的精确性与平顺性，可将优化目标函数设计为如下形式：

(16)$ J=\left|\left|\boldsymbol{Y}\left(t\right)-{\boldsymbol{Y}}_{\text{ref}}\left(t\right)\right|\right|_{\boldsymbol{Q}}^{2}+\left|\left|\boldsymbol{K}\left(t\right)\right|\right|_{\boldsymbol{R}}^{2}. $

式中：$ \boldsymbol{Q} $为精确性权重矩阵，$ \boldsymbol{R} $为平顺性权重矩阵.

将式(14)、(15)代入式(16)，并化为标准二次型：

(17)$ \left.\begin{split} & J=\frac{1}{2} \boldsymbol{K}(t)^{\mathrm{T}} \boldsymbol{H K}(t)+\boldsymbol{K}(t)^{\mathrm{T}} \boldsymbol{G}, \\& \boldsymbol{H}=2\left(\boldsymbol{\varTheta}^{\mathrm{T}} \boldsymbol{Q} \boldsymbol{\varTheta}+\boldsymbol{R}\right)+{\boldsymbol{\rho}}, \\& \boldsymbol{G}=2 \boldsymbol{\varTheta}^{\mathrm{T}} \boldsymbol{Q}\left(\boldsymbol{\varPsi} {\boldsymbol{\xi}}(t)+\boldsymbol{\varOmega}\boldsymbol{\varGamma}(t)-\boldsymbol{Y}_{\mathrm{ref}}(t)\right).\end{split}\right\} $

式中：$ \boldsymbol{H } $为标准二次型的二次项系数矩阵，$ \boldsymbol{G} $为标准二次型的一次项系数矩阵，${\boldsymbol{ \rho}} $为松弛因子.

接着引入系统约束条件，即可获得最终的优化目标函数：

(18)$ \left.\begin{split}&\mathop {\min }\limits_{\mathop \omega \nolimits_\gamma } {J}_{\text{m}}=\dfrac{1}{2}\boldsymbol{K}{\left(t\right)}^{\text{T}}\boldsymbol{H}\boldsymbol{K}\left(t\right)+\boldsymbol{K}{\left(t\right)}^{\text{T}}\boldsymbol{G};\\{{\mathrm{s.t.}}}& \left\{\begin{split} & -{\delta }_{\lim }\leqslant \delta \leqslant {\delta }_{\lim },\\& -\Delta {\delta }_{\lim }\leqslant \Delta \delta \leqslant \Delta {\delta }_{\lim }.\end{split}\right.\end{split}\right\} $

式中：$ {\delta }_{\lim } $为前轮转角极限，$\Delta {\delta }_{\lim } $为前轮转角在控制周期内的增量极限. 求解式(18)，并将获取的控制增量序列的第1个值与当前控制量相加，即可获得当前控制器应输出的控制信号.

根据式(18)可知，矩阵$ \boldsymbol{H } $具有带状稀疏结构，对稀疏、带状的QP（quadratic programming）问题，每次迭代的主要工作是求解某个带状线性方程组，其代价可用函数$ {C}_{\text{L}}\left(p\right) $表示，设这些求解器需要不超过$ {K}_{\text{QP}} $次迭代收敛，$ {K}_{\text{QP}} $为小常数或可预估上界，则LMPC的求解时间$ {T}_{\text{LMPC}} $满足：

(19)$ {T}_{\text{LMPC}}\leqslant {K}_{\text{QP}}{C}_{\text{L}}\left(p\right). $

$ {K}_{\text{QP}} $和$ {C}_{\text{L}}\left(p\right) $的行为较为可控，可以给出较稳定的最坏情况上界，因此能为LMPC给出硬性时间保证.

对于NMPC，求解方法须进行外层迭代，在每次外层迭代中须计算关于当前轨迹的雅可比矩阵以及可能的Hessian近似，并解一个线性化的QP系统. 若将计算雅可比的成本记为$ {C}_{\text{J}}\left(p\right) $，在每次外层迭代中又须解一类与LMPC类似的线性子问题，则总时间$ {T}_{\text{NMPC}} $可以估计为

(20)$ {T}_{\text{NMPC}}\approx {K}_{\text{NLP}}\left({C}_{\text{L}}\left(p\right)+{C}_{\text{J}}\left(p\right)\right). $

式中：$ {K}_{\text{NLP}} $为外层迭代次数，对于路径跟踪控制这类没有进一步结构化假设的问题，其难以保证存在小的常数上界.

对比两者可得：

(21)$ \frac{{T}_{\text{NMPC}}}{{T}_{\text{LMPC}}}\approx \frac{{K}_{\text{NLP}}}{{K}_{\text{QP}}}\left(1+\frac{{C}_{\text{J}}\left(p\right)}{{C}_{\text{L}}\left(p\right)}\right). $

$ {K}_{\text{NLP}} $通常显著大于$ {K}_{\text{QP}} $，$ {C}_{\text{J}}\left(p\right) $大于0，因此多参考点LMPC的求解时间理论上必然小于多参考点NMPC的求解时间.

类车机器人的路径跟踪控制方法通常通过MATLAB仿真^[9,11-13]或MATLAB与CarSim联合仿真^[5-6]进行验证. 鉴于CarSim在建模精度上的优势，采用基于MATLAB与CarSim搭建的联合仿真平台评估所提出的控制算法. 仿真软件的版本为MATLAB R2022a与CarSim 2017.1，运行平台为搭载Intel(R) Core(TM) i5-10200H @ 2.40 GHz处理器的Dell G15 5510笔记本电脑. 为了保证在读取实时性指标时不包含Windows系统其他任务的时间，采用MATLAB的Real-Time Synchronization模块，优先保证控制系统实时运行，使用MATLAB的tic与toc指令读取控制算法的求解时间. 被控类车机器人轴距为2 m，前轮转角范围为[−0.524, 0.524] rad，转角变化速率范围为[−0.262, 0.262] rad/s. 仿真设置的行驶速度为6 m/s，控制周期为0.02 s. 仿真设计涵盖2种典型路径场景，分别是U型弯道与单移线道路. 在仿真结果中，提出的路径跟踪控制方法记为MRP-LMPC，另外设置了3种对照控制方法，分别为NMPC、PID（proportional integral derivative）与SRP-LMPC.

U型弯是类车机器人的常见工作场景，本组仿真采用的U型弯由直线和圆弧组成，圆弧的半径为20 m. 在仿真中，为了保障实时性^[23]，预测步数取值为15，控制步数取值为2. 考虑到若航向误差和平顺性惩罚项相对位移误差的权重较大，虽然可以更好地保障稳定控制和平顺性，但是可能会出现直线行驶时残差难以消除的问题，因此为了平衡精确性和平顺性，精确性权重取值为diag {10,10,1,0}，平顺性权重取值为diag {1}，松弛因子取值为10.

如图5所示为U型弯仿真结果. 其中，图5(a)为类车机器人实际行驶路径与参考路径的对照，图5(b)为图5(a)的局部放大图. 由图可知，MRP-LMPC控制器、NMPC控制器和PID控制器均能完成路径跟踪，而SRP-LMPC控制器能够完成直线段的路径跟踪，在入弯后出现了偏离参考路径的现象. 由于在实际运行中偏离参考路径过远会导致碰撞风险，设置类车机器人偏离参考路径1 m为控制失败判定依据，因此可以判定SRP-LMPC控制器失败. 图5(c)显示了位移误差$ {e}_{\text{d}} $. 位移误差是路径跟踪控制的主要精确性指标，其定义是路径跟踪控制过程中，类车机器人偏离参考路径的距离值. 由图可知，SRP-LMPC控制器的位移误差发散，MRP-LMPC控制器的位移误差幅值虽然略大于NMPC的，但远小于PID控制器的. 在本组仿真中，排除人为设定的起始误差外，MRP-LMPC控制器的位移误差最大幅值仅为0.0971 m，整体仿真过程的平均幅值仅为0.0503 m. 图5(d)显示了航向误差$ {e}_{\text{h}} $. 航向误差是路径跟踪控制的次要精确性指标，其定义是路径跟踪控制过程中，类车机器人的航向与参考路径上与其距离最近点的航向之间的偏差值. 由图可知，SRP-LMPC控制器的航向误差发散，MRP-LMPC控制器的航向误差幅值虽然略大于NMPC的，但是远小于PID控制器的，其最大幅值仅为0.0448 rad，平均幅值仅为0.0086 rad. 图5(e)为被控对象执行的控制信号，即前轮转角，根据前轮转角幅值及其变化趋势可知，MRP-LMPC控制器和NMPC控制器一样，均未因前轮转角约束、前轮转角增量约束的影响出现大幅振荡现象. 图5(f)为实时性指标，即控制器在每个控制周期内消耗的求解时间t_s. 由图可知，MRP-LMPC控制器在实时性方面显著优于NMPC控制器，在每个控制周期内的求解时间最大值为5.95 ms，平均值为3.53 ms.

单移线是类车机器人的另一种常见工作场景，本组仿真采用的单移线由直线和圆弧组成，圆弧的半径为30 m. 仿真中控制器参数与上一组中的保持一致. 如图6所示为单移线仿真的主要结果. 其中，图6(a)为类车机器人实际行驶路径与参考路径的对照，由图可知，MRP-LMPC控制器和NMPC控制器均能完成路径跟踪，而PID控制器和SRP-LMPC控制器能够完成直线段的路径跟踪，在换道后分别出现了大幅振荡和控制失败的现象. 图6(b)显示了位移误差，由图可知，PID控制器和SRP-LMPC控制器的误差出现了发散，MRP-LMPC控制器的位移误差幅值虽然略大于NMPC，但是其最大幅值仅为0.1185 m，平均幅值仅为0.0286 m，表现出了较佳的精确性. 图6(c)为前轮转角，根据前轮转角幅值及其变化趋势可知，MRP-LMPC控制器和NMPC控制器相同，均未因前轮转角约束、前轮转角增量约束的影响出现大幅振荡现象. 图6(d)为控制器在每个控制周期内消耗的求解时间，由图可知，MRP-LMPC控制器在实时性方面显著优于NMPC控制器.

为了系统评估不同路径跟踪控制器的性能，选取主要精确性指标和实时性指标展开讨论. 精确性指标以位移误差的最大幅值衡量. 由于最大求解时间具有一定的随机性，采用平均求解时间作为实时性分析指标. 各控制器的性能指标如表1所示. 其中，e_d,max为最大位移误差，t_s,avg为平均求解时间. 可以看出，2种工况的最大位移误差存在明显差异. 单移线路径的误差普遍高于U型弯路径，这是由于单移线路径的曲率变化幅度更大，达到了0. 066 m⁻¹，使类车机器人跟踪参考路径的难度更大. 在此工况下，MRP-LMPC的误差虽大于NMPC，但最大增加幅度仅为0.0949 m. 此外，MRP-LMPC的精度下降程度远低于PID. 相比NMPC，PID的误差增加幅度达到0.3708 m. 同时，在SRP-LMPC仿真中，类车机器人行驶至弯道路径时出现误差发散，无法完成有效跟踪，这些现象表明MRP-LMPC在精确性方面明显优于PID和SRP-LMPC. 在实时性方面，仿真结果表明路径形状对平均求解时间基本没有影响. MRP-LMPC相较于NMPC能够将求解时间降低超过39.14 %，虽然仍不及PID的极低求解时间，但其计算负担已显著减轻. 上述结果表明，尽管MRP-LMPC尚无法在超低成本的单片机级控制器上直接部署，但随着工业控制计算机成本降低，该控制器在仿真条件下已能满足类车机器人路径跟踪控制的实际应用需求，同时兼顾精确性与实时性优势.

为了进一步检验控制方法在真实硬件条件下的工程可行性与可实施性，并保证实验的安全性，采用由搭载机器人操作系统（robot operating system, ROS）的工业控制计算机、MATLAB和CarSim组成的硬件在环测试平台进行实验验证. 实验平台如图7所示，ROS系统搭载于一台超恩智能EXC-1000-2G型工业计算机上，配备Intel(R) Core(TM) i7-8700 @ 3.20 GHz处理器.

控制器部署于ROS系统，通过高速局域网接收类车机器人仿真模型反馈的位姿信号，并输出控制信号. MATLAB和CarSim构成被控的类车机器人仿真模型，部署于一台Dell笔记本电脑. ROS与Simulink之间的通信与控制频率为20 Hz，路径跟踪控制器的控制频率为50 Hz，控制器的求解时间通过ROS系统的rosrun plotjuggler plotjuggler工具读取. 硬件在环实验使用的参考路径、控制器参数均与联合仿真中的一致. 另外为了使测试条件与真实环境更加一致，还在类车机器人仿真模型反馈的横坐标、纵坐标信号中分别加入了幅值为1 cm的随机误差，以模拟真实环境中的定位系统误差，在实验结果中标记为MRP-LMPC-D.

如图8所示为U型弯测试的硬件在环实验结果. 其中，图8(a)为类车机器人实际行驶路径与参考路径的对照，由图可知即使在存在定位误差的情况下，MRP-LMPC控制器也能够完成路径跟踪. 图8(b)显示了位移误差. 在不存在定位误差时，位移误差的最大幅值为0.1897 m，平均幅值为0.0647 m，在存在定位误差时，位移误差的最大幅值为0.2486 m，平均幅值为0.1165 m. 总而言之，MRP-LMPC控制器在硬件在环实验中也表现出了较高的精确性. 虽然该控制器在硬件在环实验中的误差值相对其在联合仿真中的误差值较大，但是变化趋势基本一致，控制效果趋于平稳，这也验证了提出的控制方法对于非理想环境具有较好的适应性. 此外，对比分析硬件在环实验环境和联合仿真环境的异同后可以推定，硬件在环实验的误差值较大的原因应是硬件在环实验系统中存在的轻微随机时滞. 时滞处理可以参考文献[17]、[33]和[34]，笔者计划在未来的研究工作中对该问题作进一步探讨. 图8(c)所示为前轮转角，根据前轮转角幅值及其变化趋势可知，MRP-LMPC控制器未因前轮转角约束、前轮转角增量约束的影响出现大幅振荡现象. 图8(d)为实时性指标. 由图可知，未引入随机定位误差时，MRP-LMPC的实时性表现相对较好，在每个控制周期内的求解时间最大值为5.94 ms，平均值为3.85 ms，与联合仿真中的结果基本一致. 引入随机定位误差后，实时性有所降低，在每个控制周期内的求解时间最大值为8.67 ms，平均值为5.22 ms. 不过，即使实时性有所降低，MRP-LMPC控制器的解算时间最大值仍仅为控制周期的43.35%，能够满足实时控制需求.

如图9所示为单移线测试的硬件在环实验结果. 由图9(a)可知MRP-LMPC可以在硬件在环条件下完成单移线路径跟踪. 图9(b)显示了位移误差. 在不存在定位误差时，位移误差的最大幅值为0.1660 m，平均幅值为0.0399 m，在存在定位误差时，位移误差的最大幅值为0.2221 m，平均幅值为0.0584 m. 总而言之，MRP-LMPC控制器在硬件在环实验中也表现出了较高的精确性. 图9(c)所示为前轮转角，根据前轮转角幅值及其变化趋势可知，MRP-LMPC控制器未因前轮转角约束、前轮转角增量约束的影响出现大幅振荡现象. 图9(d)所示为实时性指标. 由图可知，在未引入随机定位误差时，MRP-LMPC在每个控制周期内的求解时间最大值为7.81 ms，平均值为4.68 ms，与联合仿真中的结果基本一致. 引入随机定位误差后，在每个控制周期内的求解时间最大值为8.83 ms，平均值为5.59 ms. 解算时间的最大值为控制周期的44.15%，能够满足实时控制需求.

本研究提出基于多参考点线性MPC的控制方法，旨在解决类车机器人路径跟踪中精确性与实时性之间存在矛盾的问题. 传统控制方法通常须在精确性和实时性之间做出妥协，精确性较高的算法计算量大，实时性不足，而实时性强的算法精确性较低. 通过引入非线性补偿和多参考点策略，MRP-LMPC在保障较高精确性的前提下，显著提高了实时性.

在精确性方面，联合仿真中MRP-LMPC在U型弯路径下的最大位移误差为0.0971 m，在单移线路径下为0.1185 m. 硬件在环实验显示，在无定位误差时该控制方法的最大位移误差为0.1897 m，引入定位误差后为0.2486 m. 与相同工况下的NMPC相比，MRP-LMPC的精确性稍逊，最大位移误差略大，增加了不超过0.0949 m，但其精度仍高于SRP-LMPC和PID，后两者在弯道出现误差发散时无法有效跟踪路径.

在实时性方面，MRP-LMPC在联合仿真中的平均求解时间小于3.53 ms. 硬件在环实验显示，在表现最差的情况下平均求解时间为5.59 ms，最大求解时间为8.83 ms，最大求解时间占控制周期的43.15%. 与NMPC相比，MRP-LMPC在相同工况下可将求解时间降低超过39.14%.

综上，MRP-LMPC在本研究条件下能同时兼顾精确性与实时性，与NMPC相比在计算速度上具有明显优势，与PID和SRP-LMPC相比在精度上更高. 其可行性在硬件在环条件下得到了验证. 未来工作将围绕控制器结构优化和全实物实验进行，以进一步提升其工程可用性和可靠性.