非线性随机系统的稳定性最优控制：一种基于参数梯度的方法

《Automatica》：Stabilizing optimal control for nonlinear stochastic systems: A parametric gradient-based approach

　　伊藤裕二|藤本健二
丰田中央研发实验室有限公司，日本爱知县长久手市480-1192

**摘要**
本研究提出了一种用于设计非线性随机系统稳定次优控制器的方法。这些系统包含表示动态不确定性的时不变随机参数，这在最优控制中带来了两个关键难题。首先，时不变随机性违反了最优性原理和哈密顿-雅可比方程（Hamilton–Jacobi equations），而这些是解决最优控制问题的基本工具。其次，非线性系统必须对这些随机参数具有鲁棒稳定性。为了同时克服这些困难，本研究提出了一种基于参数梯度的方法，并加入了一个惩罚函数。控制器和成本函数使用基函数进行参数化，然后通过最小化参数化的成本函数来优化控制器。这种方法的关键挑战在于适当参数化成本函数以及求出成本的梯度。本研究提供了最优参数化成本及其梯度的显式表达式。此外，还提出了一个合适的惩罚函数，以确保即使使用梯度方法也能实现鲁棒稳定性。因此，梯度方法产生的次优反馈控制器能够保证系统的鲁棒稳定性。通过数值模拟验证了所提出方法的有效性，并展示了其与其他基线方法的性能对比。

**1. 引言**
现实世界中的各种动态系统都存在不确定性和非线性。例如，在时变不确定环境中交互的机器人操作器（Liu等人，2024年）、机器人航天器（Nakka和Chung，2023年）、鲁棒悬挂系统（Bai和Wang，2022年）、温室生产系统（Svensen等人，2024年）以及使用连续搅拌罐反应器的化学过程（Wu等人，2022年）等控制输入仿射系统中都包含不确定的物理参数和/或噪声。当系统具有有界不确定性且呈仿射关系时，问题会更易于处理。车辆编队系统（Yin等人，2022年）的机械效率和阻力系数是不确定的；流行病传播模型（Scarabaggio等人，2022年）的潜伏期和恢复率也是不确定的；纳米级振荡器模型（Ito等人，2019年）的摩擦系数和弹簧常数也是不确定的。不幸的是，这种不确定性会降低系统的控制性能和稳定性。这些例子促使我们考虑不确定非线性系统的稳定性和控制问题。

尽管系统的不确定性可以通过时不变随机参数来成功表示，但这些参数使得非线性最优控制问题变得复杂。时不变参数取遵循概率分布的随机值，但在时间上是恒定的。由于时不变参数违反了系统动态的马尔可夫性质（Markov property），它们与解决非线性最优控制的强大工具（如哈密顿-雅可比-贝尔曼方程和最优性原理）不兼容（Lewis等人，2012年）。这一限制的详细信息在第3.2节的备注1中讨论。实际上，许多随机最优控制方法（如McAllister和Rawlings，2023年；Nakka和Chung，2023年）考虑的是随时间变化的随机性，而不是时不变参数。

对于时不变随机非线性系统的控制来说，满足三个关键要求是理想的。首先，本研究关注的是离线设计随机最优控制器，而不是在线自适应控制（如强化学习RL，Baird，1995年；Jin等人，2019年；Meyn，2022年）。其次，应保证随机非线性系统的鲁棒稳定性。第三，应实现比前馈控制器更鲁棒的反馈控制器。与RL相比，满足这些要求还有以下优势：无需在线进行大量计算即可轻松实现，并且能够保证鲁棒稳定性，不受在线适应误差的影响。

不幸的是，同时满足上述三个要求具有挑战性。基于多项式混沌（PC，Bhattacharya，2014年；Blanchard，2010年；Ehlers和Pieper，2011年；Fisher和Bhattacharya，2009年；Templeton，2009年）和梯度下降（Ito等人，2024年）的控制设计方法仅适用于线性系统。PC（Xiu和Karniadakis，2002年）在处理随机不确定性方面非常有效，即使对于非线性控制问题也是如此（Lefebvre等人，2020年；Listov等人，2024年）。虽然已经开发了其他针对非线性系统的技术（Cottrill，2012年；Okura和Fujimoto，2016年；Phelps等人，2016年；Walton等人，2016年），但这些非线性控制器都是前馈式的而非反馈式的。尽管随机模型预测控制（Scarabaggio等人，2022年；Svensen等人，2024年；Yin等人，2022年）在实现非线性反馈控制方面很有前景，但它们在在线控制器设计时需要巨大的计算成本，使得实时控制变得困难。最重要的是，上述控制方法并未保证非线性随机系统的稳定性，这表明同时实现性能最优性和稳定性是具有挑战性的。

本研究提出了一种方法，用于解决包含时不变随机参数的非线性系统的最优控制问题，同时满足上述三个要求。所提出的方法提供了能够保证鲁棒稳定性的离线反馈控制器。我们的主要贡献如下：
- **简单且通用的概念**：为了解决涉及时不变随机参数的非线性最优控制问题，我们提出了一种基于梯度的方法，并对成本函数和控制器进行了参数化。此外，该方法结合了平方和（SOS）算法，以确保系统对未知随机参数的鲁棒稳定性。所提出的方法提供了一个满足三个要求的次优控制器：反馈控制、离线设计和鲁棒稳定性，而现有的设计方法尚未同时实现这些目标。
- **参数梯度理论**：本研究设计了一个适当的近似成本函数，其梯度用于梯度方法中，而不是真实梯度。由于通常无法评估无限时间范围内的成本函数，因此将无限时间范围的成本近似为参数函数，从而最小化与最优性条件相关的残差。定理1和定理2分别推导出了最优近似成本函数及其梯度的显式形式。参数函数可以表示为连续基函数的线性组合，因此得到的梯度易于实现，因为它以控制器参数的显式函数形式给出。
- **鲁棒稳定性理论**：本研究保证了即使存在未知随机参数，非线性多项式系统的鲁棒稳定性。我们在梯度方法中引入了一个惩罚函数，将基于SOS的稳定性分析与成本函数的最小化相结合。定理4保证了只要通过梯度方法成功更新控制器，鲁棒稳定性就能得到保持。定理5推导出了开始梯度方法所需的初始控制器。
- **验证**：通过一个单一的不确定多项式系统验证了所提出方法的有效性。与三种基线控制器和另一种现有方法相比，平均控制性能得到了改进。

**2. 符号说明**
- \(I_n\)：\(n\times n\)的单位矩阵
- \([v]_i\)：向量 \(v \in \mathbb{R}^n\) 的第 \(i\) 个分量
- \([M]_{i,j}\)：矩阵 \(M \in \mathbb{R}^{n\times m}\) 的第 \(i\) 行第 \(j\) 列的元素
- \(\vec{M} \equiv \begin{bmatrix} [M]_{1,1}, & \ldots, [M]_{n,1}, & \ldots, [M]_{1,m}, & \ldots, [M]_{n,m} \end{bmatrix} \in \mathbb{R}^{nm}\)：矩阵 \(M\) 的向量化形式
- \(\vec{e}_n\)：\(\begin{bmatrix} 1, & 0, & \ldots, 0 \end{bmatrix}^{n(n+1)/2}\)：对称矩阵 \(M \in \mathbb{R}^{n\times n}\) 的半向量化形式
- \(Ma \otimes Mb \in \mathbb{R}^{na\times mb}\)：矩阵 \(Ma \in \mathbb{R}^{na\times ma}\) 和 \(Mb \in \mathbb{R}^{nb\times mb}\) 的克罗内克积（Kronecker product），满足 \([Ma \otimes Mb]_{nb}(ia^{-1} + ib, mb(ja^{-1} + jb) = [Ma]_{ia, ja}[Mb]_{ib, jb}\)
- \(\lambda_{\text{min}}(M)\)：对称矩阵 \(M \in \mathbb{R}^{n\times n}\) 的最小实特征值
- \(\frac{\partial \Psi(v)}{\partial v}^T \in \mathbb{R}^{n\times m}\)：函数 \(\Psi(v) \in \mathbb{R}^n\) 关于 \(v \in \mathbb{R}^m\) 的偏导数，其中 \(\frac{\partial \Psi(v)}{\partial v}^T_{i,j} = \frac{\partial \Psi(v)}{i}{\partial v}_j\)
- \(E[\Psi(\theta)]\)：函数 \(\Psi(\theta)\) 关于随机变量 \(\theta\) 的期望值，其中 \(E[\Psi(\theta)|v]\) 表示在给定 \(v\) 的条件下的期望值

**3. 问题设定**
**3.1. 目标系统**
我们关注以下具有不确定参数 \(\theta \in \mathbb{R}^d\) 的输入仿射非线性反馈系统：
\[
x'(t, \theta) = f(x(t, \theta), \theta) + G(x(t, \theta), \theta)u(x(t, \theta),
\]
其中 \(x(t, \theta) \in \mathbb{R}^n\) 是依赖于时间 \(t\) 和 \(\theta\) 的状态。设 \(\S_\theta \equiv \{ \theta \in \mathbb{R}^d \mid p(\theta) > 0 \}\) 表示 \(\theta\) 的概率分布函数（PDF）。对于每个 \(\theta \in \S_\theta\)，设状态反馈控制器 \(u: \mathbb{R}^n \to \mathbb{R}^m\)、漂移项 \(f: \mathbb{R}^n \to \mathbb{R}^n\) 和输入向量场 \(G: \mathbb{R}^n \to \mathbb{R}^{n\times m}\) 在 \(x(t, \theta)\) 上局部利普希茨连续，并且满足 \(f(0, \theta) = 0\)。初始状态 \(x(0, \theta) \equiv x_0\) 是随机的，并假设它遵循一个概率分布 \(p(x_0)\)，使得 \(x_0\) 的每个单项式的期望值都是有限的。公式（1）中的输入仿射结构对于在第4.2节实现成本函数的最优近似、第4.3节近似成本函数的梯度以及在第4.4节将稳定性条件表示为基于SOS的约束是必要的。虽然在描述车辆动力学中的转向控制等情况下输入仿射性质是限制性的（Listov等人，2024年），但一阶泰勒展开等近似方法有助于克服这一限制。我们在本文中始终遵守假设1和假设2。

**假设1：多态非线性动力学**
存在一个正整数 \(K\)，连续函数 \(h_k: \mathbb{R}^n \to [0,1]\)，局部利普希茨连续函数 \(f_k: \mathbb{R}^n \to \mathbb{R}^n\)，以及局部利普希茨连续函数 \(G_k: \mathbb{R}^n \to \mathbb{R}^{n\times m}\)（\(k \in \{1,\ldots, K\}\)，它们满足：
\[
f(x, \theta) = \sum_{k=1}^K h_k(x),
G(x, \theta) = \sum_{k=1}^K h_k(G(x)),
\(\sum_{k=1}^K h_k(x) = 1,
\(\forall k, f_k(0) = 0,
\]
其中 \(K\)、\(h_k\)、\(f_k\) 和 \(G_k\) 对于 \(k \in \{1,\ldots, K\}\) 是已知的。假设1是为了在第4.4节保证鲁棒稳定性而必需的，因为这个假设没有考虑某些类型的不确定性（如不确定物理参数的三角函数，Listov等人，2024年）。由公式（2）、（3）给出的系统不同于线性马尔可夫决策过程（MDPs，Jin等人，2019年），因为这两种系统都是将状态表示为 \(x, u\) 的函数。我们的系统在不确定性 \(h_k(\theta)\) 上是线性的，这使得能够对每个 \(\theta\) 鲁棒地稳定系统。同时，线性MDPs在未来状态的分布上是线性的，因此动作价值函数也可以在特征映射上是线性的。

**假设2：随机参数**
- \(\theta\) 的值是未知的，并且在时间 \(t\) 上是恒定的。
- \(\theta\) 的PDF \(p(\theta)\) 是已知的。
假设2指出了本研究的关键动机：为具有时不变不确定性的系统设计控制器。在实践中，可以使用预先收集的 \(\theta\) 数据样本来估计PDF \(p(\theta)\)（Ito等人，2019年）。

**3.2. 具有鲁棒稳定性的随机最优控制**
我们的控制目标是在 \(\theta\) 的不确定性下最小化无限时间范围的成本函数 \(J(u, x_0, \theta)\)：
\[
J(u, x_0, \theta) \equiv \lim_{T \to \infty} \int_{0}^T q(x(t, \theta)) \otimes u(x(t, \theta)) \otimes R(x(t, \theta)) dt,
\]
其中 \(q: \mathbb{R}^n \to \mathbb{R}\) 是连续的且全局正定的，即 \(\lim_{|x| \to \infty} q(x) \neq 0\)。函数 \(R: \mathbb{R}^n \to \mathbb{R}^{m\times m}\) 对于每个 \(x \in \mathbb{R}^n\) 是连续的且正定的对称函数。

我们关注反馈控制器 \(u\) 设计的几个要求：由于假设2（i），控制器 \(u\) 不能直接使用未知的参数 \(\theta\)，因此控制器应该与 \(\theta\) 无关，并旨在最小化 \(J(u, x_0, \theta)\) 关于 \(\theta\) 的期望值。因为 \(E[J(u, x_0, \theta)]\) 不依赖于 \(\theta\) 的每个值，我们可以通过将 \(x(t, \theta)\) 替换为 \(x\) 来实现相应的控制器 \(u(x)\)。期望值 \(E[\cdot]\) 可以使用 \(p(\theta)\) 和 \(p(x_0)\) 来计算。此外，离线设计控制器是可取的，以便反馈系统可以实时执行。而且，在 \(\theta \in \S_\theta\) 的不确定性下，确保反馈系统的稳定性至关重要。根据这些要求，本研究提出了主要问题：
**主要问题：** 设计一个反馈控制器 \(u\)，以最小化期望成本函数并保证鲁棒的全局渐近稳定性：
\[
\min_{u} E[J(u, x_0, \theta)] \quad \text{满足条件} \quad \forall \theta \in \S_\theta, \forall x_0 \in \mathbb{R}^n, \lim_{t \to \infty} x(t, \theta) = 0.
\]
解决这个问题有两个优势：无需在线进行大量计算即可轻松实现，并且能够保证鲁棒稳定性，不受在线适应误差的影响。此外，该解决方案的控制性能比仅关注 \(\theta\) 最坏情况的鲁棒控制 \(min_{\text{max}} \theta J(u, x_0, \theta)\) 更好。类似于Bhattacharya（2014年）、Fisher和Bhattacharya（2009年）、Lewis等人（2012年）、Phelps等人（2016年）的方法，使用期望值、输入二次成本、正定状态成本和无限时间范围作为目标函数（6）、（7）是标准的。

**备注1：主要问题的难点**
主要问题涉及以下三个难点：
- 由于系统（1）的非线性，成本函数 \(J\) 和最优控制器 \(u\) 没有明确给出（没有积分形式）。
- \(\theta\) 的未知和时不变特性不满足马尔可夫性质；即存在 \(t\) 和 \(\Delta t\) 使得 \(p(x(t+\Delta t, \theta)|x(t, \theta)) \neq p(x(t+\Delta t, \theta)|x(t, \theta), x(t-\Delta t, \theta)\)。这个违反可以通过简单示例来说明：\(x(t, \theta) = [cost, \theta] \sin t^t\) 对于 \(t \in \{2\pi, 4\pi, 6\pi, \ldots\}\) 和 \(\Delta t = \pi/2\)。有效的最优控制工具（如哈密顿-雅可比-贝尔曼方程和最优性原理）不适用，因为它们基于马尔可夫性质。尽管对于增广状态 \([x^\top, \theta^\top]\) 来说马尔可夫性质成立，但由于 \(\theta\) 是未知的，因此难以实现增广状态反馈控制器。从观测值 \(x\) 估计 \(\theta\) 需要大量的在线计算，且估计误差会破坏鲁棒稳定性。此外，仅使用x而不使用θ的反馈策略对于通过相同的控制器同时控制多个不同的系统来说很有前景。(iii)控制器不仅需要实现最优性，还需要实现非线性系统(1)的鲁棒全局渐近稳定性。由于非线性和不确定性，确保稳定性是具有挑战性的。4. 提出的方法：具有鲁棒稳定性的随机最优控制4.1. 提出方法的概述为了克服与主要问题相关的三个困难(i)–(iii)，我们基于三种技术建立了提出的方法：参数近似、梯度方法和SOS。主要问题(7)被放宽为以下问题：(8)minw,P,rg(w,P,r),(9)?g(w,P,r)?E[J?(x0,v(w,θ))]+ρ(w,P,r)，其中J?表示真实成本J的近似值，它由参数函数v(w,θ)表征。符号w是表征控制器u的决策变量。其他决策变量(P,r)和惩罚函数ρ稍后引入以保证稳定性。由于以下近似，解决问题(8)会得到一个次优控制器。原始成本J被近似为J?。控制器在参数函数u?(x,w)的类别内进行优化。我们将(7)中的鲁棒稳定性约束放宽为添加ρ(w,P,r)。(8)的非凸性也导致了控制器的次优性。这些近似技术的详细贡献在下面描述。首先，使用参数近似来以显式形式表达控制器和成本函数。我们解决困难(i)的方法是在参数函数类别中设计一个次优控制器。考虑到灵活的参数函数，预期设计的控制器可以实现高控制性能。参数化提供了控制器和成本函数之间的显式关系。这种关系有助于在后续段落中介绍的控制器优化。控制器u(x)和成本函数J(u,x,θ)在(7)中被分别替换为参数函数u?(x,w)和J?(x,v(w,θ))，如下所示：(10)u(x)=u?(x,w),(11)J(u?(•,w),x,θ)=J?(x,v(w,θ))+?(x,θ,w)，其中?(x,θ,w)表示近似残差。函数J?(x,v(w,θ))和u?(x,w)分别由v(w,θ)∈Rdv和w∈Rdw参数化，以x为基函数。应该注意的是，v(w,θ)依赖于w和θ，因为J(u?(•,w),x,θ)依赖于它们。这种近似的细节在第4.2节中介绍。其次，基于梯度的方法对于克服困难(ii)是有效的，因为不需要马尔可夫性质，这与最优性原理和HJB方程不同。梯度方法用于获得(8)的次优解：(12)w{?+1}=w{?}?α???wg(w{?},P{?},r{?}),(13)P{?+1}=P{?}?α???Pρ(w{?},P{?},r{?}),(14)r{?+1}=r{?}?α???rρ(w{?},P{?},r{?})，其中上标{?}表示第?次迭代中的决策变量。经过迭代后，可以得到一个次优解和相应的反馈控制器u?(x,w{?})。梯度推导的细节在第4.3节中介绍。虽然提出的方法可以归类为强化学习领域的策略梯度方法，但该方法的优点是可以得到显式的梯度。不需要使用时间差分学习或基于观察实验/模拟的奖励之类的评估部分。第三，(9)中的惩罚函数ρ(w,P,r)帮助控制器确保鲁棒全局渐近稳定性，这与困难(iii)相关。惩罚函数是基于SOS制定的，SOS是保证多项式非线性系统稳定性的一个有前途的工具。惩罚函数的细节在第4.4节中介绍。注2提出的方法是我们在之前的方法(Ito等人，2024年)的改进版本，由于困难(i)和(iii)之前没有应用于非线性系统。4.2. 参数近似本小节介绍了控制器u?(x,w)和近似成本函数J?(x,v(w,θ))的参数化细节。它们的定义如下：(15)?u?(x,w)?Φ(x)?w,(16)?J?(x,v(w,θ))??(x)?v(w,θ)，其中Φ:Rn→Rdw×m和?:Rn→Rdv分别是局部利普希茨连续基函数和C1连续基函数，满足?(0)=0。控制器参数w∈Rdw如图(12)所示进行优化。对于每个w和θ，需要确定参数v(w,θ)∈Rdv，使得(16)中的函数J?(x,v(w,θ))近似真实成本函数J(u?(•,w),x,θ)，即最小化两个函数之间的距离。然而，由于无法直接评估无限视界成本J(u?(•,w),x,θ)，因此很难获得这个距离。为了克服上述困难，我们使用贝尔曼残差来描述J(u?(•,w),x,θ)和J?(x,v?)之间的距离，对于每个参数v?。我们定义最优参数v(w,θ)为基于贝尔曼残差B(x,v?,w,θ)的最小化器：(17)v(w,θ)∈argminv?∈Rdv(Π(B(•,v?,w,θ)2)+η‖v?‖2),(18)?B(x,v?,w,θ)?q(x)+12u?(x,w)?R(x)u?(x,w)+?J?(x,v?)?x?(f(x,θ)+G(x,θ)u?(x,w))，其中η≥0表示一个自由参数。线性算子Π(•)由(19)给出：?Π(Ψ(•))?∫Rnπ(x)Ψ(x)dx，对于任何定义在Rn上的矩阵值连续函数Ψ，其中π:Rn→[0,∞]是一个预定义的权重函数，使得Π(•)是线性的、有界的，并且对于每个连续的Ψ都是良定义的。虽然(17)中的Π(B(•,v?,w,θ)2)可以被视为具有固定策略的均方贝尔曼误差(Baird, 1995, Feng等人，2019年，Meyn, 2022年第9.3.1节)，但发现了几处不同之处。首先，我们将残差B(x,v?,w,θ)导出为依赖于θ的函数。其次，最小化器可以在定理2中明确获得。接下来，我们介绍定理1和注3，以证明(17)适用于近似J?(x,v(w,θ))。然后，我们推导定理2以获得(17)的显式解v(w,θ)。定理1 贝尔曼残差给定参数v?∈dv, w∈dw, 和 θ∈Sθ，假设带有u(x)=u?(x,w)的反馈系统(1)是全局渐近稳定的，并且对于每个x∈Rn，都有J(u?(•,w),x,θ)<∞。如果存在β∈[0,1)使得(20)?x∈Rn,B(x,v?,w,θ)2≤β2(q(x)+12u?(x,w)?R(x)u?(x,w))2，则我们有(21)?x∈Rn,J?(x,v?)1+β≤J(u?(•,w),x,θ)≤J?(x,v?)1?β。证明在附录A中描述。注3 定理1的贡献定理1表明，将近似成本J?(x,v?)除以1±β可以提供J(u?(•,w),x,θ)的下界和上界。较小的非负β值可以提供更紧密的界限。这种关系证明了最小化贝尔曼残差Π(B(•,v?,w,θ)2的必要性，这是(17)中的第一项。(17)中的第二项确保了(17)的解v(w,θ)是唯一确定的。定理1与逆向动态规划(IDP) (Meyn, 2022年，命题3.8)不同。定理1指出了给定控制器的J?(x,v?)的近似误差，用于找到最优的J?(x,v?)。同时，IDP比较了最优控制器和近似控制器之间的成本函数值。定理2 最优近似的解假设假设1成立。假设η>0，Π(‖ψ0‖2)<∞，并且对于k∈{1,…,K}，有Π(‖ψk‖2)<∞，其中(22)?ψ0(x)?[q(x),vec(Φ(x)R(x)Φ(x)?/2)?]?，(23)?ψk(x)?vec(??(x)?x?[fk(x),Gk(x)Φ(x)?])。对于任何w∈Rdw和任何θ∈Sθ，(17)的解v(w,θ)由(24)v(w,θ)=?(∑k=1K∑k′=1Khk(θ)hk′(θ)Lk,k′(w))?1∑k=1Khk(θ)lk(w)唯一确定，其中(25)?Lk,k′(w)?(1w?Idv)?Π(ψkψk′?)(1w?Idv)+ηIdv，(26)?lk(w)?(1w?Idv)?Π(ψkψ0?)1w?w。此外，如果∑k=1K∑k′=1Khk(θ)hk′(θ)Lk,k′(w)?0，则不需要假设η>0。证明在附录B中描述。注4 定理2的贡献定理2给出了(17)的显式解v(w,θ)。将此解代入(16)可以得到近似成本函数J?(x0,v(w,θ))=?(x)?v(w,θ)的显式形式。这个成本是真实成本J(u?(•,w),x,θ)的最优近似，因为(17)中的基于贝尔曼的距离被最小化了。如果η>0或∑k=1K∑k′=1Khk(θ)hk′(θ)Π(ψkψk′?)?0成立，那么对于每个w，我们有?L?(w,θ)?∑k=1K∑k′=1Khk(θ)hk′(θ)Lk,k′(w)?0。如果(19)中的π(x)定义得当，即π(x)>0对于各种x都成立，则期望正定性成立。第5.1节中给出了一个简单的例子。在L?(w,θ)?1的计算数值稳定的范围内，η的值应该是小的。此外，我们不直接测量B(x,v?,w,θ)来获得v(w,θ)。最小化器v(w,θ)由(24)中已知的(w,θ)函数给出，因为除了θ之外，系统动态是已知的。4.3. 近似成本函数的梯度为了实现(12)中提出的梯度方法，需要近似成本函数J?(x0,v(w,θ))的梯度。我们在下面推导梯度。定理3 成本函数的显式梯度假设定理2中假设的所有条件都成立。对于任何w∈Rdw和任何θ∈Sθ，带有(24)的近似成本J?(x0,v(w,θ))的梯度由(27)给出：??[w]jE[J?(x0,v(w,θ))|θ]=?E[?(x0)]?(∑k=1K∑k′=1Khk(θ)hk′(θ)Lk,k′(w))?1×∑k=1Khk(θ)(?lk(w)?[w]j+∑k′=1Khk′(θ)?Lk,k′(w)?[w]jv(w,θ))。此外，如果∑k=1K∑k′=1Khk(θ)hk′(θ)Lk,k′(w)?0在w的邻域内成立，则不需要假设η>0。证明在附录C中描述。注5 定理3的贡献定理3提供了近似成本函数J?(x0,v(w,θ))的显式梯度。这个梯度使得可以在(12)中实现梯度方法来解决最小化问题(8)。4.4. 基于SOS的稳定性与梯度方法在本小节中，确保了带有u(x)=u?(x,w)的反馈系统(1)的鲁棒全局渐近稳定性，同时梯度方法优化控制器。我们保证稳定性的策略基于SOS。对于单项式基函数z(x)∈Rdz，定义了Lyapunov函数V(x)及其时间导数V?(x,θ)的以下SOS形式：(28)?V(x)?z(x)?Pz(x)，(29)V?(x(t,θ),θ)=?2∑k=1Khk(θ)z(x(t,θ))?Uk(x(t,θ),P,w)z(x(t,θ))，其中P∈Rdz×dz和Uk(x,P,w)∈Rdz×dz分别是最小化(8)中的决策变量和稍后要澄清的函数。如果对于每个θ，V(x)和?V?(x,θ)都是全局正定的，并且V(x)是径向无界的，即V(x)>0且V?(x,θ)<0对于x≠0，V(0)=V?(0,θ)=0，lim‖x‖→∞V(x)=∞，以及lim‖x‖→∞V?(x,θ)≠0，则满足鲁棒全局渐近稳定性。多项式Lyapunov函数(28)适用于多项式非线性系统，因为使用基于SOS的方法可以提供作为凸不等式的系统稳定性条件，这些条件易于验证。虽然这种策略基于开创性的结果(Xu等人，2009年，定理1)和Prajna等人(2004年，定理6)，但提出的方法保证了系统关于θ的鲁棒稳定性。此外，惩罚函数使得可以将稳定性保证与通过梯度方法减少的预期成本函数结合起来。本小节的目标是设计惩罚函数ρ(w,P,r)和初始决策变量(w{0},P{0},r{0})，以便在(12)–(14)中的梯度方法中实现前一段中介绍的基于SOS的策略。我们在下面介绍关键定理来证明设计。基于引理1，我们介绍了惩罚函数ρ(w,P,r)和决策变量r的细节。定理4表明所提出的惩罚函数保证了鲁棒全局渐近稳定性。定理5提供了在梯度方法开始时使用的适当初始值(w{0},P{0},r{0})。以下定义和单项式基的假设被使用。定义1 单项式基如果zˉ:Rn→Rd1ˉ×d2ˉ的每个分量都是原始单项式∏i=1n[x]i[a]i对于某些a∈{0,1,2,…}n，则称zˉ(x)为单项式基。如果满足zˉ(x)=0?x=0，则称单项式基zˉ(x)为严格的。如果zˉ(x)的所有分量都是zˉ(x)的分量并且是不同的单项式，即对于每个i≠j，[zˉ′(x)]i≠[zˉ′(x)]j对于某些x成立，则称单项式基zˉ′(x)为zˉ(x)的非冗余形式。注6 引入基函数设z:Rn→Rdz和ζ:Rn→Rdζ为单项式基，使得z(x)是严格的并且[ζ(x)]1=1。设ξ:Rn→Rdξ为ζ(x)?z(x)∈Rdzdζ的非冗余形式。设Z:Rn→Rm×dZ为连续函数，使得Gk(x)Z(x)是多项式。设计这些基的进一步指导原则对应于以下假设，这些假设对于基于SOS保证稳定性是必要的。虽然目标反馈系统限于多项式，但使用高阶单项式可以提高系统的表达能力。假设3 SOS设置基z(x)，ζ(x)和Z(x)被设计为满足以下条件：(i)对于每个k∈{1,2,…,K}，存在已知的多项式矩阵Fk:Rn→Rn×dz满足(30)fk(x)=Fk(x)z(x)。(ii)在控制器u?(x,w)中使用的基函数Φ(x)在(15)中被设置为(31)Φ(x)=z(x)?Z(x)?。(iii)对于任何P∈Rdz×dz，任何参数w∈Rdw，以及任何k∈{1,2,…,K}，以下函数：(32)?Uk(x,P,w)??P+P?2?z(x)?x?(Fk(x)+Gk(x)Z(x)W)是由基ζ(x)?ζ(x)组成的多项式矩阵，即每个(i,j)th分量满足[Uk(x此外，该解决方案在(w,P)上是双线性的，在r上是线性的。证明详见附录D。备注7：引理1的贡献由于引理1中Tk(w,P,r)的唯一性，(34)中的惩罚函数ρ(w,P,r)是明确定义的。确定常数(Cr,cr)非常直接，因为(35)–(37)简化为与x无关的Tk的简单线性方程，如第5.1节所示。使用明确定义的惩罚函数ρ(w,P,r)，我们得出了以下结果：梯度方法保证了稳定性。定理4：适当的梯度方法假设假设1和假设3成立。给定初始决策变量(w{0},P{0},r{0})和常数ρˉ∈R(34)，假设满足以下条件：(38)P{0}+P{0}??0，(39)Tk(w{0},P{0},r{0})?0。对于任何?∈{0,1,2,…}，如果梯度方法(12)–(14)提供的(w{?},P{?},r{?})满足(40)g(w{?},P{?},r{?})≤g(w{0},P{0},r{0})，(41)E[J?(x0,v(w{?},θ))]>g(w{0},P{0},r{0})?ρˉ，那么具有u(x)=u?(x,w{?})的反馈系统(1)的鲁棒全局渐近稳定性得到满足：(42)?θ∈Sθ,?x0∈Rn,limt→∞x(t,θ)=0。证明详见附录E。备注8：定理4的贡献定理4意味着设计的控制器满足鲁棒全局渐近稳定性(42)，而梯度方法优化控制器参数w。梯度方法可以自然满足条件(40)。通过将ρˉ设置为足够大的值，很容易满足(41)。如果定理4中的所有条件都成立，ρˉ的值不会影响控制器设计，因为ρ(w{?},P{?},r{?})的梯度与ρˉ无关。由于定理4要求初始决策变量(w{0},P{0},r{0})满足(38), (39)，因此它们是这样设计的。我们设计初始决策变量作为以下两个优化问题的解：(43), (48)：(43)maxQ,H,S1,…,SK,?1s.t.?k,?x,Yk(x,Q,H)+Yk(x,Q,H)?=(ζ(x)?Idz)?Sk(ζ(x)?Idz)，?k,Sk=Sk???1Idζdz,Q=Q???1Idz,?1>0。其中(44)?Yk(x,Q,H)???z(x)?x?(Fk(x)Q+Gk(x)Z(x)H)。如果获得了问题(43)的可行解(Q?,H?)，则初始决策变量由(45)?w{0}?vec(H?Q??1)，(46)?P{0}?Q??1，(47)?r{0}?r?定义，其中r?是以下问题的可行解：(48)maxr,?2?2s.t.?k,Tk(w{0},P{0},r)??2Idξ,?2>0。注意Tk(w{0},P{0},r)是满足(w,P)=(w{0},P{0})的(35)–(37)的唯一解。我们推导出以下定理以保证这些优化问题的可行性以及初始决策变量的合理性。定理5：初始控制器设计假设假设1和假设3成立，我们有以下结论：(i)存在问题(43)的可行解(Q?,H?,S1?,…,SK?,?1?)。(ii)给定问题(43)的可行解，存在问题(48)的可行解r?，对于某些Cr∈Rdr×dξ(dξ+1)/2和cr∈{0,1}在(36)中。(iii)在(45)–(47)中的初始决策变量(w{0},P{0},r{0})满足条件(38), (39)。证明详见附录F。备注9：定理5的贡献定理5确保在(45)–(47)中定义的初始决策变量总是可以获得的，并且满足定理4所需的条件(38), (39)。因此，所提出的梯度方法产生了稳定的次优控制器。此外，优化问题(43), (48)很容易解决，因为它们通过将SOS约束转换为线性等式约束而简化为凸半定规划(SDP)。具体来说，(43)中第一个约束的每个(i,j)分量是一个多项式，表示为zˉ(x)?bk,i,j(Q,H,Sk)=0，且?bk,i,j(Q,H,Sk)?Ak,i,j[vec(Q)?,vec(H)?,vec(Sk)?,1]?，其中zˉ(x)是一个单项式基，Ak,i,j是一个常数矩阵。因为这是一个零多项式，SOS简化为线性等式：bk,i,j(Q,H,Sk)=0。4.5：次优性分析在本小节中，我们分析了两个方面：所提出的梯度方法的次优性以及原始主问题(7)与其松弛(8)之间的差距。在前一方面，所提出的梯度方法适当地减少了目标函数g(w,P,r)，原因有二。第一个原因是g(w,P,r)以标准的基于梯度的方式减少(40)。这种单调非递增的方式至少产生局部最优和/或稳定解。第二个原因是，在鲁棒稳定性约束下，最小化g(w,P,r)几乎等同于最小化E[J?(x0,v(w,θ))]，因为g(w,P,r)是基于对数障碍方法(Boyd & Vandenberghe, 2004第11.2节)建立的。这种最小化差距的细节如下推导。设?SwPr?{(w,P,r)|P+P??0,?k,Tk(w,P,r)?0}表示保证稳定性条件的(w,P,r)集合。SwPr的一个内部子集由?λλλλSwPro?{(w,P,r)|λˉIdz?(P+P?)/2?λ?Idz,?k,λˉIdξ?Tk(w,P,r)?λ?Idξ}定义，对于自由参数λλλˉ>λ?>0。SwPr中可行的w集合由?Sw?{w|?P,?r,(w,P,r)∈SwPr}定义。然后，以下结果表明，最优性差距δa(κ)与ρ(w,P,r)中包含的系数κ≥0成正比，并且可以通过使用较小的κ来减少。定理6：假设E[J?(x0,v(w,θ))]在Sw上根据(24)是明确定义的。假设对于每个κ≥0，最小化器(w??,P??,r??)∈argmin(w,P,r)∈SwPrg(w,P,r)包含在SwPro中。那么，我们有δa(κ)-次优性：(49)E[J?(x0,v(w??,θ))]≤minw∈SwE[J?(x0,v(w,θ))+δa(κ)，(50)?λλδa(κ)?κ(dz+Kdξ)ln(λˉ/λ?)。证明详见附录G。接下来，我们分析原始主问题(7)与其松弛(8)之间的差距。以下定理表明，所提出的控制器u?(•,w??)是δb(κ,β1,β2)-次优的，其中δb(κ,β1,β2)是根据参数(κ,β1,β2)确定的，用于描述所提出方法的精度。较小的非负值(κ,β1,β2)会产生较小的差距δb(κ,β1,β2)。定理7：假设定理6中使用的所有假设都成立。假设问题(7)的解由u?(•,w?)给出，对于某些w?∈Sw。对于每个w∈Sw和每个θ∈Sθ，假设系统(1)在u(•)=u?(•,w)下是全局渐近稳定的，每个x∈Rn满足J(u?(•,w),x,θ)<∞，并且期望E[J(u?(•,w),x0,θ)]是明确定义的。如果存在β1∈[0,1)和β2∈[0,1)，使得对于每个x∈Rn和每个θ∈Sθ，我们有(51)B(x,v(w?,θ),w?,θ)2≤β12(q(x)+12u?(x,w?)?R(x)u?(x,w?))2，(52)B(x,v(w??,θ),w??,θ)2≤β22(q(x)+12u?(x,w??)?R(x)u?(x,w??))2，那么控制器u?(•,w??)是δb(κ,β1,β2)-次优的：(53)E[J(u?(•,w??),x0,θ)]≤E[J(u?(•,w?),x0,θ)]+δb(κ,β1,β2)，(54)?δb(κ,β1,β2)?β1+β2(1+β1)(1?β2)E[?(x0)]?E[v(w??,θ)]+δa(κ)1+β1。证明详见附录H。4.6：实现下载：下载高分辨率图像（232KB）下载：下载全尺寸图像最后，所提出的方法在算法1中进行了总结，该算法包括两个主要部分：初始决策变量的设计和控制器的优化。以下备注详细说明了算法的几个技术方面。备注10：期望值的近似我们旨在近似梯度方法(12)中使用的关于θ的期望E[…]，因为对于具有无限集合Sθ的许多类型的PDF p(θ)，期望值是不可计算的，而对于有限集合Sθ，期望值是可以计算的。可以使用几种方法来近似期望值，例如使用梯形法则的数值积分和蒙特卡洛近似。使用这些技术，(12)中的梯度可以近似为?E[J?(x0,v(w{?},θ))]/?w≈(1/S)(∑s=1S?v(w{?},θs,?)?/?w)E[?(x0)]，其中θs,?是根据近似技术采样的参数，S表示样本数量。选择S的值时，需要考虑近似精度和计算成本之间的权衡。较大的S可以提高精度，但计算量也会大大增加。因为所提出的方法基于通用梯度方法，所以它与其他算法兼容。例如，随机梯度下降方法（Sun等人，2020年第III节）、Adam（Kingma & Ba，2015年）和概率线搜索（Mahsereci & Hennig，2017年）都适用于执行梯度方法。预计使用这些技术可以减少控制器设计中的计算时间。备注11：步长在第7行中，建议确定(12)–(14)中使用的步长α?，以满足条件(40)。一种合适的方法是使用满足Wolfe条件（Nocedal & Wright，2006年第3.1节）的步长α?：(55)g?(y{?+1})≤g?(y{?})?χα?‖?g?(y{?})/?y‖2，其中y和g?(y{?})分别表示(w,P,r)的向量值集合和目标函数g(w{?},P{?},r{?})，为了简洁起见。符号χ∈(0,1)是一个自由参数。可以使用回溯方法（Nocedal & Wright，2006年第3.1节算法）来找到这样的步长α?。在每次迭代?时，我们首先设置α?=αini，其中αini是一个较大的正值。然后我们重复α?←γαα?，直到满足(55)，其中γα∈(0,1)是一个自由参数。5. 数值示例5.1：植物系统和设置让我们考虑以下具有n=2的系统(1)的示例：(56)?f(x,θ)?[x]1+[x]12+([θ]1?2)[x]13?[x]1[x]22/2+[x]2([θ]2+1)[x]1+[x]22，(57)?G(θ)?[θ]1[θ]2+1。不确定参数θ遵循二维有限集?Sθ?{0,0.1,0.9,1.0}2上的均匀分布，其中Sθ的基数相对较小，因此我们简化了各种类型控制器和参数的模拟计算负担。这个示例通过以下设置满足假设1，其中K=4：fk(x)=f(x,θk)，Gk=G(θk)，h1(θ)=(1?[θ]1)(1?[θ]2)，h2(θ)=(1?[θ]1)[θ]2，h3(θ)=[θ]1(1?[θ]2)，h4(θ)=[θ]1[θ]2，其中?θ1?[0,0]?，?θ2?[0,1]?，?θ3?[1,0]?，?θ4?[1,1]?。我们根据假设3在备注6中定义单项式基z(x)，ζ(x)，Z(x)和ξ(x)为?z(x)?[[x]1,[x]2]?，?ζ(x)?[1,[x]1,[x]2]?，?Z(x)?[1,[x]1,[x]2,[x]12,[x]1[x]2,[x]22]，以及?ξ(x)?[[x]1,[x]2,[x]12,[x]1[x]2,[x]22]?。因为(56)中的f(x,θ)分解为f(x,θ)=F(x,θ)z(x)，并且(58)?F(x,θ)?1+[x]1+([θ]1?2)[x]12?[x]22/21[θ]2+1[x]2，所以满足fk(x)=Fk(x)z(x)的矩阵Fk(x)由Fk(x)=F(x,θk)给出。在这些设置下，约束(36)简化为三个方程：[Tk]5,3=[rk]1，[Tk]3,2=[rk]2，[Tk]4,2=[rk]3。选择较差的基函数Z(x)会导致控制器设计失败，因为不满足稳定性条件(43)和/或(48)。如果Z(x)从[1,[x]1,[x]2]，[1,[x]1,[x]2,[x]12]，[1,[x]1,[x]2,[x]1,[x]2],[1,[x]1,[x]2,[x]12,[x]1[x]2]中选择，那么控制器设计就会失败。这表明使用高阶和高维的Z(x)是合适的。在随机最优控制器的设计中，(6)中的成本函数定义为q(x)=x?x，R=10。PDF p(x0)设置为在{?3,0,3}2?{[0,0]?}上的均匀分布。基?(x)包括所有度数大于0且小于或等于6的单项式，即[x]1d1[x]2d2对于所有(d1,d2)∈{0,1,2,…,6}2且0在另一方面，所提出的方法在理论上保证了稳定性，而DMPPI即使在建模中的状态在数值上是稳定的也无法保证这一点。此外，这种比较对于评估所提出方法的最优性非常有启发性。当N?≥2000时，所提出方法的预期成本低于DMPPI的成本。从这个意义上说，所提出的控制器已经得到了充分的优化。

表1. 预期成本E[J(u,x,θ)]表示了关于θ和x0的平均控制性能，以及在使用MATLAB 2024b（CPU：Intel(R) Xeon(R) Platinum 8488C 2.40 GHz；RAM：256 GB）时执行的计算时间。

| 控制器 | 预期成本 | 计算时间（秒） |
|-------|---------|---------|
| DMPPI(Δtc=5) | 发散 | 0.099 |
| DMPPI(Δtc=0.5) | 发散 | 0.921 |
| DMPPI(Δtc=0.05) | 941 | 0.973 |
| DMPPI(Δtc=0.005) | 781 | 105.079 |
| DMPPI(Δtc=0.0005) | 771 | 109.365 |
| 未优化 | 1359 | 0.088 |
| 无最优性（N?=10000） | 1130 | 0.095 |
| 无稳定性（N?=10000） | 发散 | 0.096 |
| 所提出的方法（N?=100） | 760 | 0.088 |
| 所提出的方法（N?=500） | 1000 | 0.090 |
| 所提出的方法（N?=1000） | 860 | 0.090 |
| 所提出的方法（N?=2000） | 730 | 0.090 |
| 所提出的方法（N?=8000） | 600 | 0.088 |
| 空单元（N?=10000） | 590 | 0.090 |

6. 结论

本研究提出了一种用于设计时不变非线性随机系统的次优控制器的方法，该方法满足了三个基本要求：反馈控制、离线设计和鲁棒稳定性。我们通过提出成本函数的参数近似、基于梯度的控制器优化以及基于SOS的惩罚函数，克服了与非线性最优控制、时不变随机特性和鲁棒稳定性相关的挑战。

未来的研究将致力于将所提出的方法扩展到其他类别的非线性随机系统，例如随机多面体系统，其线性情况已在Hosoe等人（2020年）的研究中进行分析。我们预期，关注确定性强化学习（RL）方法与随机RL方法之间的差距（Meyn，2022年）可以为将我们的研究扩展到一般随机过程提供有价值的思路。一个关键的挑战是使用非多项式基函数来保证鲁棒稳定性。预计使用核函数和神经网络等非多项式基函数可以进一步提高控制性能。一种可能的方法是将控制器分解为多项式项和非多项式项。根据所提出的方法设计出多项式项后，再调整非多项式项以增强控制性能。其他有前景的研究方向包括扩展到其他控制问题，例如风险敏感控制（Jacobson，1973年），以及应用于实际系统，包括具有制造变异的纳米级设备（Ito等人，2019年）。