模型预测控制(model predictive control，MPC)基于系统模型的预测求解有限时域最优控制问题，并将最优控制序列第一元素作用于系统，滚动向前^[1-2]. 由于滚动优化机制，MPC具有本质鲁棒性^[3]. 为了提高鲁棒性，近年来通过在最优控制问题中显式考虑扰动边界，提出一系列鲁棒MPC (robust MPC，RMPC)策略^[3-8]. 由于RMPC直接考虑扰动所有可能的实现，存在严重的保守性质.

在实际中，扰动通常以典型“实现”作用于系统，引入机会约束和场景优化可以降低RMPC保守性^{[4, 9-15]}. 场景优化将时域内扰动的实现视作“场景”，构建RMPC随机最优控制问题^[13]. 若控制律对所有场景满足约束，则实际系统以较高的置信度满足预设机会约束. Calafiore等^[14]针对线性变参系统，提出在线抽取场景的RMPC策略. Lorenzen等^[15]离线抽取场景以保证机会约束满足，给出系统稳定所需的场景数量边界. 通常，场景优化由随机凸优化理论保证机会约束概率，较高概率需要在扰动实现分布中抽取大量场景^[16-18]，但实际中难以获取扰动实现的完整概率分布，这限制了场景优化的应用范围.

为了建立少量数据下的约束满足，可以使用学习方法离线建模^[19-21]，再设计RMPC策略^[22]，如集员^[23]与Kinky推断方法^[24]. 这些方法能够从数据中学习扰动边界，但无法得到扰动随机分布，实现机会约束. 高斯过程(Gaussian process，GP)^[25-26]能够同时提供最优预测值与预测方差，已用于随机MPC设计. Li等^[26]基于GP提出解析近似MPC算法，处理加性不确定性问题. Grancharova等^[27]使用GP设计随机MPC，用于燃烧控制. Wang等^[28]基于GP给出针对饮用水网络的MPC算法. Hewing等^[19-20]结合GP和MPC，实现车辆的自动驾驶控制. 现有算法在理论上缺少稳定性与可行性结果^[27-28]，且假定扰动为独立同分布，限制了该算法的理论拓展与应用范围.

本文考虑约束不确定线性系统，提出基于高斯回归的场景优化鲁棒模型预测控制策略. 基于先验数据集进行高斯回归学习，获得扰动的完整概率模型. 在线抽取大量场景构建松弛机会约束下的随机凸优化问题，通过开环优化、闭环反馈镇定系统. 建立保证闭环系统鲁棒稳定性和满足松弛机会约束的充分条件. 应用DC-DC转换器和车辆巡航仿真，验证了本文策略的有效性.

考虑离散时间不确定线性系统：

(1) $ {{\boldsymbol{x}}_{k+1}} = {\boldsymbol{A}}{{\boldsymbol{x}}_k}+{\boldsymbol{B}}{{\boldsymbol{u}}_k}+{{\boldsymbol{w}}_k};\;k \in {N_0}. $

式中：x_k、u_k和w_k分别为k时刻的系统状态变量、控制变量和未知不确定扰动，x_k∈Rⁿ，u_k∈R^m，w_k∈Rⁿ；A和B为已知矩阵；N_i为所有不小于i的整数集合. 系统状态和控制输入满足如下约束：

(2) $ {{\boldsymbol{x}}_k} \in X,\;{{\boldsymbol{u}}_k} \in U,\;\;k \in {N_0}. $

式中：X = {x∈Rⁿ| f_X(x)≤0}和U = {u∈R^m| f_U(u)≤0}为紧凸集，f_X(x)和f_U(u)分别为x和u的凸函数. 假设系统状态可测，矩阵对(A, B)可镇定，即存在增益矩阵K_f使A+BK_f稳定. 系统存在给定的先验数据集D，其中包含长度为L的s条历史运行轨迹. 先验数据集是学习型控制律的基本条件^{[19, 24, 26]}，可由给定的标准测试数据、专家知识或系统运行历史数据产生^[29-31]，更多关于先验数据集获取的知识可以参见文献[32].

假设1 　系统(1)的不确定扰动w未知但有界，即w_k∈W, $k \in N_{\rm{0}} $，其中集合W $\subset$Rⁿ为紧集.

假设2 　系统(1)存在包含原点的一个凸邻域X_f及其局部控制律u = K_f x，使X_f为闭环系统的鲁棒正不变集(robust positive invariant set，RPI)，即X_f={x∈X| (A+BK_f)x+w∈X_f, K_f x∈U, $\forall $x∈X_f, w∈W}.

注1 　假设1的扰动可以描述系统参数摄动、外界附加扰动、建模误差等^[3]. 对于假设2中的终端域X_f，Kouvaritakis等^[33]给出求解最大RPI(maximal RPI, MRPI)作为X_f 的计算方法，其中K_f通过求解系统(1)的Riccati方程得到.

考虑系统(1)，用K_f参数化控制输入^[26]：

(3) $ {\boldsymbol{u}} = {{\boldsymbol{K}}_{\rm{f}}}{\boldsymbol{x}}+{\boldsymbol{v}}. $

式中：v为控制摄动变量，通过在线滚动优化计算得到. 参数化控制律通过将MPC开环优化问题转化为可计算的闭环控制律优化问题，可以降低RMPC控制器的保守性，使得算法可以用于开环不稳定系统^{[2-4, 33]}. 令A_K = A+BK_f，则系统(1)等价于

(4) $ {{\boldsymbol{x}}_{k+1}} = {{\boldsymbol{A}}{\rm{_K}}}{{\boldsymbol{x}}_k}+{{\boldsymbol{Bv}}_k}+{{\boldsymbol{w}}_k};\;k \in {N_0}. $

本文的目标是计算系统(1)的控制律(3)，使得闭环系统在任意有界扰动w作用下，状态进入并始终保持在集合X_f内，同时系统状态和控制输入满足预设概率为P_p的机会约束：

(5) $ {P}({{\boldsymbol{u}}_k} \in U,{{\boldsymbol{x}}_k} \in X) \geqslant P_{\rm{p}};\;k \in {N_1}. $

采用高斯过程回归和模型预测控制方法计算控制律(3).

通过GP学习时间指标k到扰动w_k的映射，基于预测时域内扰动实现的联合高斯分布抽取场景，构建RMPC优化问题，设计新的场景优化RMPC算法.

考虑系统(1)及先验数据集D，有

(6) $ {{\boldsymbol{w}}_k} = {{\boldsymbol{x}}_{k+1}} - ({{\boldsymbol{Ax}}_k}+{{\boldsymbol{Bu}}_k}). $

先验数据集D中包含任意时刻k的s个扰动值w_k，其中 $k \in I_0^{L - 1} $( $I_0^{L - 1} $表示集合{0, 1,···, L−1}).

GP由均值与协方差函数唯一确定. 令向量k_D = [0, 1, 2,···, L]^T，使用先验数据集D对GP进行训练学习. 在学习完毕后，对于待预测的输入数据向量k^*，有如下的联合高斯分布：

(7) $ \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{k}}_D}} \\ {{{\boldsymbol{k}}^*}} \end{array}} \right]\sim {\rm{N}}\left( {\left[ {\begin{array}{*{20}{c}} {{\boldsymbol{E}}({{\boldsymbol{k}}_D})} \\ {{\boldsymbol{E}}({{\boldsymbol{k}}^*})} \end{array}} \right],\left[ {\begin{array}{*{20}{c}} {{\boldsymbol{K}}({{\boldsymbol{k}}_D},{{\boldsymbol{k}}_D})}&{{\boldsymbol{K}}{{({{\boldsymbol{k}}^*},{{\boldsymbol{k}}_D})}^{\rm{T}}}} \\ {{\boldsymbol{K}}({{\boldsymbol{k}}^*},{{\boldsymbol{k}}_D})}&{{\boldsymbol{K}}({{\boldsymbol{k}}^*},{{\boldsymbol{k}}^*})} \end{array}} \right]} \right). $

GP是概率模型，均值向量的行数与协方差阵的维度等于输入时间指标的数目.

由高斯分布的性质可知，k^*相关输出的条件分布服从高斯分布^[5]，均值向量与协方差矩阵分别为

(8) $ {\boldsymbol{E}}({{\boldsymbol{w}}_{{k^*}}}) = {\boldsymbol{K}}({{\boldsymbol{k}}^*},{{\boldsymbol{k}}_D}){\boldsymbol{K}}{({{\boldsymbol{k}}_D},{{\boldsymbol{k}}_D})^{ - 1}}{{\boldsymbol{k}}_D}, $

(9) $ \begin{gathered} {\boldsymbol{K}}({{\boldsymbol{k}}^*},{{\boldsymbol{k}}^*}) = {\boldsymbol{K}}({{\boldsymbol{k}}_D},{{\boldsymbol{k}}^*}) - {\boldsymbol{K}}({{\boldsymbol{k}}^*},{{\boldsymbol{k}}_D}){\boldsymbol{K}}{({{\boldsymbol{k}}_D},{{\boldsymbol{k}}_D})^{ - 1}}{\boldsymbol{K}}{({{\boldsymbol{k}}^*},{{\boldsymbol{k}}_D})^{\rm{T}}}. \end{gathered} $

式中：协方差矩阵K的每一项由预先指定的核函数表述. 选择核函数为径向基函数：

(10) $ k({k_1},{k_2}) = {\sigma ^2}\exp\; [ - {{{{({k_1} - {k_2})}^2}} /( {2{l^2}})}]. $

式中：k₁与k₂为任意2个时间指标标量;待求解参数σ² 用于控制输出信号方差；l反映输入尺度放缩，可以在离线阶段基于D通过极大似然法或者最大后验估计求得. w(k^*)的预测是由式(8)、(9)决定的联合高斯分布，其维数等于向量k^*中的元素个数. 在系统运行的任一时刻，可以认为预测时域内扰动分布已知，能够在线抽取大量场景.

令x_t为当前时刻t的系统状态测量，根据GP给出的预测时域内的联合高斯分布，在时刻t抽取M个扰动场景 ${\boldsymbol{\delta}}^{(1)}_t$, ${{\boldsymbol{\delta}}}^{(2)}_t$,···, ${{\boldsymbol{\delta}}}^{(M)}_t $，其中第i个场景为 $\boldsymbol{\delta}_t^{(i)}=\left[\boldsymbol{w}_{0 \mid t}^{(i)}, \boldsymbol{w}_{1| t}^{(i)}, \cdots, \boldsymbol{w}_{N-1 \mid t}^{(i)}\right]\left(i \in I_1^M\right)$，即预测时域内扰动w的一种可能实现. 选定权重矩阵Q、R、P>0，对第i个扰动场景定义二次型目标函数^{[26, 33]}：

(11) $ J^{(i)}=\sum_{j=0}^{N-1}\left(\boldsymbol{x}_{j \mid t}^{(i) \mathrm{T}} \boldsymbol{Q} \boldsymbol{x}_{j \mid t}^{(i)}+\boldsymbol{u}_{j \mid t}^{\mathrm{T}} \boldsymbol{R} \boldsymbol{u}_{j \mid t}\right)+\boldsymbol{x}_{N \mid t}^{(i) \mathrm{T}} \boldsymbol{P} \boldsymbol{x}_{N \mid t}^{(i)} .$

构建时刻t的有限时域最优控制问题如下：

$\tag{12a} \mathop {\min }\limits_{{{\boldsymbol{V}}_t},{z_t},{q_t}} \;\;{z_t}+\gamma {q_t}. $

$\tag{12b} {\rm{s.t.}}\;\;\;{J^{(i)}} \leqslant \;{z_t}, $

$\tag{12c} {\boldsymbol{u}}_{j|t}^{(i)} = {{\boldsymbol{K}}_{\rm{f}}}{\boldsymbol{x}}_{j|t}^{(i)}+{{\boldsymbol{v}}_{j|t}}, $

$\tag{12d} {\boldsymbol{x}}_{j+1|t}^{(i)} = {{\boldsymbol{A}}_{\rm{K}}}{\boldsymbol{x}}_{j|t}^{(i)}+{{\boldsymbol{Bv}}_{j|t}}+{\boldsymbol{w}}_{j|t}^{(i)}, $

$\tag{12e} {f_X}({\boldsymbol{x}}_{j+1|t}^{(i)}) - {\bf{1}}{q_t} \leqslant {\boldsymbol{0}}, $

$\tag{12f} {f_U}({\boldsymbol{u}}_{j|t}^{(i)}) - {\bf{1}}{q_t} \leqslant {\boldsymbol{0}}, $

$\tag{12g} {f_{{X_{\rm{f}}}}}({\boldsymbol{x}}_{N|t}^{(i)}) - {\bf{1}}{q_t} \leqslant {\boldsymbol{0}},\; $

$\tag{12h} {\boldsymbol{x}}_{0|t}^{(i)} = {{\boldsymbol{x}}_t},\; $

$\tag{12i} {q_t} \geqslant 0,\;i \in I_1^M,\;j \in I_0^{N - 1}. $

控制量由摄动序列 ${\boldsymbol{V}}_t=\left[{\boldsymbol{v}}_{0 \mid t}, {\boldsymbol{v}}_{1 \mid t}, \cdots, {\boldsymbol{v}}_{N-1 \mid t}\right]$、约束松弛变量q_t和所有场景对应的目标函数上界z_t构成. 式(12b)为目标函数的上界约束，式(12c)为预镇定状态反馈律，式(12d)为基于抽取场景的模型预测，式(12e)~(12g)分别为状态、输入和终端域X_f约束松弛，式(12h)为优化问题的初始条件，式(12i)保证约束被松弛.1是维度适当的全1列向量，γ>0为q_t对应的权重，当q_t无穷大或γ = 0时，约束默认满足；当γ无穷大或q_t=0时，该优化问题对应于无松弛RMPC问题.q_t反映系统在当前状态下对约束的违反程度. 设定γ为一个较大的常数，以惩罚约束违反，保证优化问题对任意抽取的场景均可行，实现控制器鲁棒性和约束松弛的折中. 优化问题(12)的解是当前状态x_t的隐函数，受模型参数(A, B, f_X, f_U)和设计参数(N, M, Q, R, γ)的综合影响. 由于扰动场景是随机抽取的，对于相同的初始条件，抽取的场景不同，会得到不同的最优解.

在线求解优化问题(12)，得到最优解 $q_t^* $、 $z_t^* $、 ${\boldsymbol{V}}_t^*=\left[{\boldsymbol{v}}_{0 \mid t}^*, {\boldsymbol{v}}_{1 \mid t}^*, \cdots, {\boldsymbol{v}}_{N-1 \mid t}^*\right] $. 根据MPC滚动优化原理，计算场景优化RMPC控制量为

(13) $ {\boldsymbol{u}}_t^{{\rm{MPC}}} = {{\boldsymbol{K}}_{\rm{f}}}{{\boldsymbol{x}}_t}+{\boldsymbol{v}}_{0|t}^*, $

并作用于系统(1). 在下一时刻，测量系统状态，更新优化问题(12)的初始条件，再重复整个计算过程，实现场景优化RMPC的滚动时域优化控制.

引理1^[18] 　考虑如下优化问题：

(14) $ \mathop {\min }\limits_{\boldsymbol{s}} \;\;\{ {{\boldsymbol{c}}^{\rm{T}}}{\boldsymbol{s}};\;\;{h_{\boldsymbol{\delta}} }({\boldsymbol{s}}) \leqslant {\boldsymbol{0}}\} . $

式中：c∈R^n×1为常数向量， ${h_{\boldsymbol{\delta}} }\left( \cdot \right) $为取决于随机变量 ${\boldsymbol{\delta}} $的凸约束. 设定常数β∈(0, 1.0)与机会约束的预设概率P_p，若随机抽取 ${\boldsymbol{\delta}} ^{(i)} $的场景数目M使得下式成立：

(15) $ \sum\limits_{j = 0}^{d - 1} {\left( {\begin{array}{*{20}{c}} M \\ j \end{array}} \right)} {(1 - P_{\rm{p}})^j}{P_{\rm{p}}^{M - j}} \leqslant \beta , $

其中d为优化变量s的维度，则优化问题

(16) $ \mathop {\min }\limits_{\boldsymbol{s}} \;\;\{ {{\boldsymbol{c}}^{\rm{T}}}{\boldsymbol{s}};\quad {h_{{\delta ^{(i)}}}}({\boldsymbol{s}}) \leqslant {\boldsymbol{0}},\;i \in I_1^M\} $

的解s^*以1−β的置信度满足P(h_δ(s)≤0) ≥ P_p.

注2 　式(15)的等号仅对全支撑(fully-supported)问题成立^[18]. 在GP学习中，先验数据可能不严格满足高斯过程，或者由于先验数据的不完备使得GP预测存在偏差，式(15)的保守性质缓解了由先验数据带来的控制器性能下降或闭环系统失控问题.

对于优化问题(12)，应用引理1的结果，可得关于场景优化非滚动优化控制量的性质.

引理2 　考虑优化问题(12)及最优摄动序列解 ${\boldsymbol{V}}_t^*=\left[{\boldsymbol{v}}_{0 \mid t}^*, \cdots, {\boldsymbol{v}}_{N-1 \mid t}^*\right] $. 若假设1和2成立，且P_p、β和M满足式(15)，则在t, t+1, ···, t+N−1时刻依次调用V_t^*的元素形成控制量(3)，序列 ${\boldsymbol{V}}_t^* $以1 − β置信度，不低于P_p的概率和约束违反q_t满足：

a）将系统(1)在N步内驱动到终端域X_f内，即

(17) $ {{P}}({{P}}({f}_{{X}_{{\rm{f}}}}({{\boldsymbol{x}}}_{t+N})-{\bf{1}}{q}_{t}\leqslant {\boldsymbol{0}})\geqslant P_{\rm{p}}）\geqslant 1-\beta . $

b）使系统(1)满足输入和状态约束(2)，即

(18a) $\tag{18a} {{P}}({{P}}({f}_{U}({{\boldsymbol{u}}}_{t+j})-{\bf{1}}{q}_{t}\leqslant {\boldsymbol{0}})\geqslant P_{\rm{p}})\geqslant 1-\beta , $

(18b) $\tag{18b} {{P}}({{P}}({f}_{X}({{\boldsymbol{x}}}_{t+j})-{\bf{1}}{q}_{t}\leqslant {\boldsymbol{0}})\geqslant P_{\rm{p}})\geqslant 1-\beta ,\;j\in {I}_{1}^{N}. $

为了设计新的场景优化RMPC 算法，定义位于终端域X_f中的以原点为中心、φ>0为半径的闭球B_φ={x| ||x||≤φ} $\subset $X_f. 由假设2可知，X_f为包含原点在其内部的紧集，故B_φ总是存在. 由式(11)、(12b)，有z ≥ J ≥ x^TQx ≥ 0，故B_φ外的系统状态x对应的最优值函数上界z必有下界 $\underline z $.

系统(1)的场景优化RMPC算法描述如下.

算法1　场景优化RMPC算法

离线阶段(初始化)：

1）设定GP核函数，基于先验数据集D学习GP参数. 设定P_p与β，计算M和终端域X_f ；设定正常数d，使得 $\underline z $ ≥ld，l∈N₁；令t = 0.

2）测量状态x_t，求解优化问题(12)得到最优解 $\left\{q_t^*, z_t^*, {\boldsymbol{V}}_t^*=\left[{\boldsymbol{v}}_{0 \mid t}^*, {\boldsymbol{v}}_{1 \mid t}^*, \cdots, \dot{{\boldsymbol{v}}}_{{N}-1 \mid t}^*\right]\right\} $，计算控制量u_t^MPC并作用于系统(1). 令 $q_t=q_t^*, z_t=z_t^*, \boldsymbol{V}_t=\boldsymbol{V}_t^* $.

在线阶段：

3）令t=t+1, $\hat{\boldsymbol{V}}_t=\left[\boldsymbol{v}_{1 \mid t-1}, \cdots, \boldsymbol{v}_{N-1 \mid t-1},{\boldsymbol{0}}\right]$, $\hat{q}_t=q_{t-1} $及 $\hat{z}_ t=\max\; \left\{0, z_{t-1}-d\right\}$. 求解优化问题(12)，得到最优解 $\left\{q_t^*, z_t^*, \boldsymbol{V}_t^*\right\} $.

4）判定如下条件：

a）若 $z_t^*>\hat{z}_t $且 $\hat{z}_t<\boldsymbol{x}_t^{\mathrm{T}} \boldsymbol{x}_t $，则令 $\boldsymbol{V}_t=\hat{\boldsymbol{V}}_t, q_t=\hat{q}_t, z_t=0 $.

b）若 $z_t^*>\hat{z}_ t $且 $\hat{z}_t \geqslant \boldsymbol{x}_t^{\mathrm{T}} \boldsymbol{x}_t $，则令 $\boldsymbol{V}_t=\hat{\boldsymbol{V}}_t, q_t=\hat{q}_{t}, z_t=\hat{z}_t$.

c）若 $z_t^* \leqslant \hat{z}_t$，则令 $\boldsymbol{V}_t=\boldsymbol{V}_t^*, q_t=q_t^*, z_t=z_t^* $.

5）由V_t第一元素计算控制量u_t^MPC并作用于系统(1)，返回3）.

在算法1中， $\left\{q_t^*, z_t^*, {\boldsymbol{V}}_t^*\right\} $表示时刻t实时求解获得的最优解， $\left\{\hat{q}_t, \hat{z}_t, \hat{\boldsymbol{V}}_t\right\} $表示备选解，{q_t, z_t, V_t}表示用于系统实际控制的解. 球半径φ可以通过求解如下的优化问题获得：

(19) $ \varphi = \mathop {\min }\limits_{\boldsymbol{r}} \;\{ \left\| {\boldsymbol{r}} \right\|;\;\;{\boldsymbol{r}} \in \partial {X_{\rm{f}}}\} . $

式中： $\partial {X_{\rm{f}}} $为集合X_f的边界. 进一步可以计算B_φ外目标函数下界：

(20) $ \underline{z} = \mathop {\min }\limits_{\boldsymbol{x}}\; \{ {{\boldsymbol{x}}^{\rm{T}}}{\boldsymbol{Qx}};\;\;{{\boldsymbol{x}}^{\rm{T}}}{\boldsymbol{x}}\; \geqslant \;{\varphi ^2}\} . $

注3 　传统场景优化RMPC的稳定性由离线设计的矩阵参数保证^[13-14]，本文只需要设计φ与 $\underline{z} $2个标量参数，简化了场景优化RMPC算法的分析与设计.

下面给出场景优化RMPC的机会约束满足和闭环系统的鲁棒稳定性分析.

定理1 　若假设1和2成立，且算法参数{M, P_p, β}满足式(15)，则闭环系统(1)在算法1控制下将以1−β置信度和松弛q_t满足所有概率为P_p的机会约束，即

(21a) $\tag{21a} {{P}}({{P}}({f}_{{X}_{{\rm{f}}}}({{\boldsymbol{x}}}_{{{t}}+{{N}}})-{\boldsymbol{1}}{q}_{t}\leqslant {\boldsymbol{0}})\geqslant P_{\rm{p}}）\geqslant 1-\beta , $

(21b) $\tag{21b} {{P}}({{P}}({f}_{U}({{\boldsymbol{u}}}_{t+j})-{\boldsymbol{1}}{q}_{t}\leqslant {\boldsymbol{0}})\geqslant P_{\rm{p}}）\geqslant 1-\beta , $

(21c) $\tag{21c} {{P}}({{P}}({f}_{X}({{\boldsymbol{x}}}_{t+j})-{\bf{1}}{q}_{t}\leqslant {\boldsymbol{0}})\geqslant P_{\rm{p}})\geqslant 1-\beta ,j\in {I}_{1}^{N}. $

证明：由算法1的第4步可知，在任意时刻有q_t = $\hat{q} _t $或 $q_t=q_t^* $，下面分类讨论.

1）若q_t = $\hat{q} _t $，由算法1第3步可得q_t = q_t−1及V_t = $\hat{{\boldsymbol{V}}} _t $，将引理2应用于t−1时刻，可知不等式(21)成立.

2）若 $q_t=q_t^* $，则 ${\boldsymbol{V}}_t={\boldsymbol{V}}_t^* $. 将引理2用于当前时刻t，可知不等式(21)成立，从而定理1结果成立. 证毕.

定理2 　若假设1和2成立，算法参数{M, P_p, β}满足式(15)，则闭环系统(1) 在算法1控制下以不低于1−β的置信度，以概率P_p收敛于终端域X_f或在有限时间内到达终端域.

证明： 根据算法1的第4步a）条件是否满足，分类讨论如下.

1）不满足第4步的a）条件，即控制器运行4.b）或4.c）. 根据第3步中 $\hat {z}_t $的取值，分别讨论如下.

1.1）若在第3步 $\hat {z}_t=\max \;\left\{0, z_{t-1}-d\right\}=0 $，则在4.b）时有 ${\boldsymbol{x}}_t^{\mathrm{T}} {\boldsymbol{x}}_t<\hat{z}_t=0 $，即x_t = 0；或在4.c）时有 $z_t^* \leqslant \hat{z}_t=0 $，即x_t = 0（显然在优化问题(12)中 $z_t^* $=0时必有x_t = 0）. 因此，系统收敛至终端域中的原点.

1.2）若在第3步 $\hat{z}_t=\max \;\left\{0, z_{t-1}-d\right\}=z_{t-1}-d $，则4.b）时满足递减性条件 $z_t=\hat{z}_t=z_{t-1}-d $，或4.c）时有 $z_t=z_t^* \leq \hat z_t=z_{t-1}-d $. 考虑到d的性质，只要闭环状态位于B_φ外，递减性条件max {0, z_t−1−d}=z_t−1−d成立，从而闭环系统必定渐近稳定于B_φ. 由于B_φ$ \subset $X_f，则系统状态收敛于终端域. 在该情况下，系统从初始时刻到终端域的总时间不高于 $\left\lceil {z_0^*/d} \right\rceil $，因为到达终端域的时间不高于到达闭球B_φ的时间，更不高于到达原点的时间.

2）当前时刻满足第4步的a）条件. 此时按下一时刻的情况分类考虑如下.

2.1）若下一时刻满足第4步的a）条件，则实际上使用备选解产生RMPC控制量（在极端情况下，使用初始时刻t = 0时的最优解V₀^* 产生RMPC控制量，即非滚动优化算法）. 根据定理1可知，系统以不低于1−β的置信度和松弛q_t，以P_p的概率满足约束且收敛于终端域.

2.2）若下一时刻满足4.b）或4.c），则由上述证明过程1）中的讨论可知，系统收敛于终端域. 证毕.

考虑如下DC-DC转换器的线性时不变模型^[26]：

(22) $ {{\boldsymbol{x}}_{k+1}} = \left[ {\begin{array}{*{20}{c}} 1&{0.007\;5} \\ { - 0.143}&{0.115} \end{array}} \right]{{\boldsymbol{x}}_k}+\left[ {\begin{array}{*{20}{c}} {4.798} \\ {0.115} \end{array}} \right]{{{u}}_k}+\left[ {\begin{array}{*{20}{c}} 1 \\ 1 \end{array}} \right]{{{w}}_k}. $

式中：扰动信号w_k由持续扰动项0.01sin k和不确定项e_k组成，即w_k = 0.01sin k + e_k，在任一k时刻，e的取值分布为均值为0、方差为0.005²的截断高斯分布. 系统状态、控制输入和扰动约束集合分别为X = {x∈R²| [−2, −3]^T≤x≤[2, 3]^T}, U = {u∈R|−0.2≤u≤0.2}和W = {w∈R| −0.02≤w≤0.02}.

在GP离线训练阶段，在约束范围内随机生成长度L = 16的50条经验轨迹数据，拟合结果如图1所示. 从图1可知，GP能够准确拟合非线性三角函数，且双侧0.8置信区间可正确覆盖约80%的点，说明GP能够正确拟合时间相关扰动分布.

在线控制阶段中，选择预测时域的时间尺度N = 8，置信度β = 10⁻⁵，权重矩阵

$ {\boldsymbol{Q}}=\left[ {\begin{array}{*{20}{c}} 1&0\\ 0&{10} \end{array}} \right]，$

R = 1，γ = 10⁴，根据式(19)、(20)计算得到φ = 0.345 9，z = 0.119 7. 常数d = 0.05和l = 2，初始点x₀ = [3.4, 3.5]^T.

为了验证本文算法的机会约束满足和低保守性，使用10 000次蒙特卡洛实验测定机会约束频率p_c与优化问题求解时间t_s，对比滚动场景优化RMPC（简记为算法1）与非滚动优化设计（简记为算法2）. 在算法1中，基于实时松弛量q_t 构建松弛约束，并判定其是否违反. 在算法2中，约束判定依赖于t=0时优化问题的松弛约束，在一次模拟中的任何约束违反视作该模拟失败. 2种算法的机会约束统计频率与求解时间如表1所示.

分析表1可知，当P_p = 0.7或0.9时，M > s，此时满足高概率机会约束需要GP学习. 随着P_p值上升，M增大，统计频率上升说明控制器的鲁棒性更强. 2种算法的实际机会约束概率均高于P_p，即满足机会约束，但存在一定的保守性，这是由于优化问题(12)仅在初始时刻接近全支撑问题. 算法1的保守性更低，这得益于松弛变量的实时更新. 在算法2下，由于初始点离终端域最远，q₀值过大导致后续约束违反较难. 算法2由于滚动优化实时刷新q_t，降低了保守性. 在时效性上，t_s随M的增大而变长，这是由于式(12)的约束数目随着M的增大而增加，但优化问题的严格凸性和可行性保证了该问题的快速求解.

图2~4验证了算法1的可行性与松弛变量的必要性. 图2中，x₁、x₂分别为系统状态x的第一和第二分量，由于x₀的一步可达集合X₁与状态约束集合X无交集，有必要引入松弛变量q₀保证RMPC控制律的可行性. 如图3、4所示为当P_p = 0.1时算法1下系统的闭环状态和输入轨迹，图3中的状态轨迹全部收敛于终端域，验证了本文算法的可行性. 如图4所示为图3对应的输入轨迹，算法1在t = 0时违反输入约束，松弛量典型值为q₀ = 0.083，说明了松弛变量的有效性.

考虑网联车自主巡航控制系统^[34]：

(23) $ \Delta d = d - {d_{{\rm{des}}}},\;\Delta v = {v_{\rm{p}}} - {v_{\rm{h}}},\;{\tau _{\rm{d}}}{\dot a_{\rm{h}}}+{a_{\rm{h}}} = {a_{{\rm{des}}}}. $

式中：Δd为车间距误差；d为实际车间距，d = x_p−x_h−l，其中x_p和x_h分别为前、后车的位移，l为前车车长；d_des为理想车间距，d_des = t_hv_h+d₀，其中v_h为后车的速度，t_h为车头时距，d₀为最小安全间距；Δv为相对速度，v_p为前车的速度；τ_d为巡航时间常数，a_h为后车的加速度， $\dot{a}_{\mathrm{h}} $为a_h的导数，a_des为期望加速度. 在车联网环境下，后车可以通过网络通信获取当前时刻的前车加速度和速度信息.

选择状态变量x = [Δd, Δv, a_h]^T，输入为a_des，扰动为a_p（其中a_p为前车的加速度），则巡航控制系统(24)的离散时间状态方程为

(24) $ {{\boldsymbol{x}}_{k+1}} = {{\boldsymbol{A}}_k}{{\boldsymbol{x}}_k}+{{\boldsymbol{B}}_k}{{{u}}_k}+{{\boldsymbol{G}}_k}{{{w}}_k}. $

式中：

$ {{\boldsymbol{A}}_k} = \left[ {\begin{array}{*{20}{c}} 1&{{T_{\rm{c}}}}&{ - {t_{\rm{h}}}{T_{\rm{c}}}} \\ 0&1&{ - {T_{\rm{c}}}} \\ 0&0&{1 - {T_{\rm{c}}}/{\tau _{\rm{d}}}} \end{array}} \right],\;{{\boldsymbol{B}}_k} = \left[ {\begin{array}{*{20}{c}} 0 \\ 0 \\ {{T_{\rm{c}}}/{\tau _{\rm{d}}}} \end{array}} \right],\;{{\boldsymbol{G}}_k} = \left[ {\begin{array}{*{20}{c}} 0 \\ {{T_{\rm{c}}}} \\ 0 \end{array}} \right] , $

其中T_c为采样离散时间，T_c > 0.

仿真中令t_h = 2.5 s，T_c = 0.4 s，τ_d= 0.5，γ = 1 500，d = 4，R = 60，Q为三阶单位阵，N = 8，先验轨迹数s = 15，场景数目M = 25. 令v = G_kw，系统扰动范围为W={v| −1≤||v||≤1}，且方差为0.12²，状态和控制约束为

$ X = \left\{ {{\boldsymbol{x}} \in {{\bf{R}}^3}\left| {\left[ {\begin{array}{*{20}{c}} { - 6} \\ { - 5} \\ { - 4} \end{array}} \right] \leqslant {\boldsymbol{x}} \leqslant \left[ {\begin{array}{*{20}{c}} 6 \\ 5 \\ 5 \end{array}} \right]} \right.} \right\},U = \{ - 4,5\}. $

使用算法1与算法2，在加速和减速2种工况中进行对比实验.

考虑减速停车工况，设定前车初始速度为4.26 m/s，仿真初始状态x₀ = [8, 2, 1]^T，使用离散sigmoid函数模拟前车速度，对应的前车加速度如表2所示. 由于前车最终停止，设定最终v_k为0. 分别使用算法1和算法2控制系统，绘制被控后车的加速度、速度、间距误差与系统状态轨迹曲线，如图5~8所示.

分析图5可知，算法1控制下的自主车辆加速度更接近前车. 图6、7中，算法1控制的自主车辆速度与间距误差变化更加平缓，意味着自主巡航车辆拥有更好的舒适性与跟踪性能. 图8中，2种算法作用下的车辆加速度、速度和理想间距均能渐近稳定于原点. 由于初始状态位于约束X外部且距离较远，系统初始轨迹位于约束X外部，此时须引入松弛因子q₀保证可行性.

考虑加速工况，设定系统初始状态x₀ = [−9, −3, −1]^T，前车初始速度为15.21 m/s，加速度如表2所示. 仿真结果如图9~12所示. 类似于减速工况分析可知，算法1控制下的自主车辆加速度更加接近于前车加速度. 自主车辆速度与间距误差更加平稳，意味着算法1控制下的车辆在前车加速工况下具有更好的舒适性与跟踪性能. 由于前车速度变化后自主巡航控制系统存在不确定扰动v_k，系统状态不会渐近稳定于原点，但始终位于终端域中.

本文设计基于GP学习的场景优化RMPC算法. 去除了传统场景优化分布已知或足够场景可获得的假设，通过使用GP学习扰动实现的取值分布，抽取场景构建随机凸优化问题，求解控制摄动序列. 引入松弛变量保证算法可行性，基于随机凸优化引理证明了闭环系统的鲁棒稳定性. DC-DC转换器和网联车巡航控制仿真结果验证了本文算法的有效性和优越性. 后续将进一步研究非准确学习下的RMPC算法及低场景数和高效算法设计.