参数化最优控制问题的初值校正法

An Initial Value Correction Method for the Parameterized Optimal Control Problem

摘要: 针对带随机系数椭圆方程约束的最优控制问题，我们提出一种结合共轭梯度法的初值校正新方法。该方法通过少量的样本求解最优控制，然后得到的解作为抽取其它样本的初始迭代值，从而减少整个系统求解最优控制平均值的成本。我们理论上说明该方法的合理性。两个数值例子表明我们的方法是充分优于未初值校正方法。

Abstract: For optimal control problems with random elliptic equation constraints, we propose an initial value correction new method combined with the conjugate gradient method. This method solves the optimal control through a small number of samples, and then the solution is used as the initial iteration value of extracting other samples, thus reducing the cost of solving the optimal control average for the whole system. We theoretically illustrate the rationality of the method. Two numerical examples show that our method is sufficiently superior to the no initial value correction method.

1. 引言

最优控制问题广泛应用于工程和应用科学领域。对于确定性最优控制问题，存在着大量的数学理论和计算方法 [1] [2] [3] [4] [5]。然而，边界条件、材料参数往往不能被精确地测量，所以随机模型更现实。与确定性最优控制问题一样，随机最优控制问题很难得到解析解，数值解是应用这些问题的一个很好的选择。在过去的几十年里，随机最优控制问题的数值研究已经引起了人们的关注 [6] [7] [8] [9]。

对于随机系数的每一个固定实现，我们求解一个确定的最优控制系统。然后，估计最优控制的统计量(如：期望值、方差)。这使我们能够研究控制对偏微分方程约束中的随机波动的敏感性，并为受不确定输入影响的控制装置的设计提供实际信息。迄今为止，具有随机系数的椭圆偏微分方程已经得到了很好的研究，并且有大量有关这方面数值方法 [10] [11]。具有随机系数的最优控制问题的研究相对少一些。 [9] 研究了受随机椭圆型偏微分方程约束的随机最优控制问题的最优解可测性。在传统有限元方法的框架下， [12] 研究了随机最优控制问题最优解的存在唯一性。

众所周知，偏微分方程约束优化问题的数值解在计算上是昂贵的，因为它需要数值求解由最优性条件产生的偏微分方程系统的解：状态问题、伴随问题。特别是偏微分方程的参数随机化时，这使计算变得更具挑战性，并可能导致高维随机空间中的“维数诅咒”。目前已有一些求解随机最优控制问题的有效方法。

A. Borzi和G. Winckel将随机配置法与多重网格结合求解具有随机系数抛物型和椭圆型最优控制问题 [3]。O. Ghattas and A. Alexanderian，针对均值方差风险规避最优控制问题，运用泰勒近似，简化原本最优控制模型 [13]。蒙特卡罗方法(MC)是在随机空间中得到近似最优解的最直接的方案之一。然而，它的收敛率很低，并且抽取每个样本都必须解决一个完全确定性的最优控制问题，导致沉重的计算成本。为克服MC的不足，M. Hinze，F. Kuo，S. Vandewalle等人使用多水平蒙特卡洛方法和拟蒙特卡洛方法模拟最优控制的均值，降低MC的计算成本 [9] [14] [15]。但是他们没有降低单个样本成本，而是从方差减少的角度降低计算整个最优控制的样本量。

很多时候我们使用MC抽样时，没有充分利用已抽样本的信息，这也是导致计算花费增加一个原因。为更好利用已抽样本的信息，我们提出初值校正的方法，该方法通过少量的样本求解最优控制平均值，然后得到的平均值作为抽取其它样本的初始迭代值，从而减少求解最优控制统计量的成本。结合理论分析和数值实验，充分说明了方法实用性。

文章主要分为四部分，第一部分给出模型及一些基本的假设。第二部分给出最优控制问题的梯度及有限元离散模型。第三部分分析求解最优控制问题的误差，并提出基于初值校正的算法，分析算法的可靠性。第四部分给出两个数值算例，比较了初值校正与初值未校正算法的差异。

2. 模型及假设

我们用 $\left(\Omega ,\mathcal{A},P\right)$ 表示一个概率空间，其中 $\Omega$ 是一个样本空间， $\mathcal{A}\subset 2^{\Omega }$ 是一个 $\sigma$ 代数， $P:\mathcal{A}\to \left[0,1\right]$ 是一个概率空间，给定巴拿赫空间 $\left(X,{\Vert \cdot \Vert }_X\right)$， $u\in L^p\left(\Omega ,X\right)$ 是强可测函数，并满足 $u:\Omega \to X$， ${\Vert u\Vert }_{L^p\left(\Omega ,X\right)}<\infty$，其中：

${\Vert u\Vert }_{L^p\left(\Omega ,X\right)}=\{\begin{array}{l}{\left({\displaystyle {\int }_{\Omega }{\Vert u\Vert }_X^p\text{d}P\left(\omega \right)}\right)}^{\frac{1}{p}}1\le p<\infty ,\\ {\text{esssup}}_{\omega \in \Omega }{\Vert u\Vert }_{X}p=\infty .\end{array}$

随机变量v的期望算子，方差算子表示如下：

$E\left[u\right]={\displaystyle {\int }_{\Omega }u}\text{d}P\left(\omega \right),Var\left(u\right)=E\left[{\left(u-E\left[u\right]\right)}^2\right]=E\left[u^2\right]-E{\left[u\right]}^2.$

2.1. 状态方程与目标函数

对于我们提出的方法并不假定任何特定的随机偏微分方程，然而，为了阐明方法的思想，我们考虑状态方程是随机椭圆方程约束下的最优控制问题。

空间区域 $D={\left[0,1\right]}^2\subset R^2$， $\partial D$ 表示边界。在D上温度分布由状态y描述。控制u是D上的一个热源，这里我们假设u与时间无关。热传导系数是一个随机场k： $D\times \Omega \to R:\left(x,\omega \right)\to k\left(x,\omega \right)$。整个系统用下面带随机系数的PDE描述：

$\begin{array}{l}-\nabla \cdot \left(k\left(x,\omega \right)\nabla y\left(x,\omega \right)\right)=u\left(x\right),\left(x,\omega \right)\in D\times \Omega \\ y\left(x,\omega \right)=0,\left(x,\omega \right)\in \partial D\times \Omega \end{array}$ (1)

$k_{\min }\left(\omega \right)={\min }_{x\in D}k\left(x,\omega \right),k_{\max }\left(\omega \right)={\max }_{x\in D}k\left(x,\omega \right).$

我们给出以下假设：

A1: $k_{\min }\left(\omega \right)>0$，a.s.， $k_{\min }^{-1}\in L^p\left(\Omega \right)$，其中 $p\in \left[1,\infty \right)$，

A2: $k\in L^p\left(\Omega ,C^1\left(\bar{D}\right)\right)$，其中 $p\in \left[1,\infty \right)$，

A3: $u\in L^2\left(D\right)$。

进一步我们写出(1)的变分形式

$\left(k\left(x,\omega \right)\nabla y\left(x,\omega \right),\nabla v\right)=\left(u,v\right),v\in H_0^1\left(D\right).$

由文献 [11] [13] 知，上述方程几乎确定有解，并且 $y\in L^2\left(\Omega ,H^2\left(D\right)\right)\cap L^2\left(\Omega ,H_0^1\left(D\right)\right)$。

我们研究下面参数化的二次型目标函数

$J=\frac{1}{2}{\Vert y\left(x,\omega \right)-y_d\left(x\right)\Vert }_{L^2\left(D\right)}^2+\frac{\lambda }{2}{\Vert u\left(x\right)\Vert }_{L^2\left(D\right)}^2,$

其他不同类型的目标函数可参考 [14] [16]，注意到这里控制与参数 $\omega$ 无关。 $\lambda$ 表示正则化参数，对每一个固定的参数 $\omega$，状态y与u满足方程(1)。

2.2. 随机场 $k\left(x,\omega \right)$

我们假设k是对数随机场，即： $k\left(x,\omega \right)=\exp \left(z\left(x,\omega \right)\right)$，其中z是高斯随机场， $E\left[z\left(x,\omega \right)\right]=0$，根据文献 [17]，高斯随机场的协方差定义为：

$C\left(x_1,x_2\right)=Cov\left(z\left(x_1,\omega \right),z\left(x_2,\omega \right)\right)={\sigma }^2\exp \left(-\frac{{\Vert x_1-x_2\Vert }_1}{a}\right),$

其中 ${\sigma }^2$ 是随机场的方差，a是相关长度。为得到z的样本，通常方法有循环抽样 [18]，K-L展开 [19]。本文运用K-L开展

$z\left(x,\omega \right)=E\left[z\left(x,\omega \right)\right]+{\displaystyle \underset{n=1}{\overset{\infty }{\sum }}\sqrt{{\theta }_n}}{\xi }_n\left(\omega \right)f_n\left(x\right).$ 有关K-L展开的详细内容可参考 [18]，我们选择相关长度 $a=0.3$，截断数 $n_{KL}=500$。关于K-L展开抽样产生的误差不是文章关注的重点。图1给出不同方差下随机场的实现。

Figure 1. Realization of random field

图1. 随机场的实现

3. 最优控制的梯度及离散最优控制

为估计最优控制期望，我们采用MC进行抽样。固定随机参数 $\omega$，得到一个参数化的最优控制问题：

$\underset{u\in L^2\left(D\right)}{\min }{J}^{\omega }\left(u\right)=\frac{1}{2}{\Vert y^{\omega }-y_d\Vert }_{L^2\left(D\right)}^2+\frac{\lambda }{2}{\Vert u\Vert }^2,$

$-\left(k^{\omega }\nabla y^{\omega },\nabla v\right)=\left(u,v\right),v\in H_0^1\left(D\right),$

其中 $y^{\omega }$ 表示依赖随机参数 $\omega$，标准的证明方法 [1] [9]，我们可以知道最优控制有唯一解 $u^{*,\omega }$，同时可以得到最优控制的梯度：

$\begin{array}{l}\left(k^{\omega }{y}^{\omega },\nabla v\right)=\left(u,v\right),v\in H_0^1\left(D\right),\\ \left(k^{\omega }\nabla p^{\omega },\nabla v\right)=\left(y^{\omega }-y_d,v\right),v\in H_0^1\left(D\right)\\ \nabla J^{\omega }\left(u\right)=\lambda u+p^{\omega }.\end{array}$

下面我们离散这个系统，设 ${\mathcal{T}}_h$ 是D的三角形剖分，最大网格尺寸 $h=\underset{\mathcal{T}\in {\mathcal{T}}_h}{\max }diam\left(\mathcal{T}\right)$，并且 $\bar{D}={\displaystyle \underset{\mathcal{T}\in {\mathcal{T}}_h}{\cup }\bar{\mathcal{T}}}$。假定三角形剖分是拟一致的，我们定义有限元空间：

$V_h=\{v_h\in C^0\left(\bar{D}\right):v_h$ 在每个三角形单元 $\mathcal{T}\in {\mathcal{T}}_h$ 是线性多项式，且 ${v_h|}_{\partial D}=0\}$。

得到下面有限元离散形式

$\begin{array}{l}\left(k^{\omega }\nabla y_h^{\omega },\nabla v\right)=\left(u,v\right),v\in V_h,\\ \left(k^{\omega }\nabla p_h^{\omega },\nabla v\right)=\left(y_h^{\omega }-y_d,v\right),v\in V_h\\ \nabla J^{\omega }\left(u\right)=\lambda u+p_h^{\omega }.\end{array}$

进一步我们得到矩阵形式：

$\begin{array}{l}{A}_h{y}^{\omega }=M_{h}u,\\ A_h{p}^{\omega }=M_h\left(y^{\omega }-y_d\right),\\ \nabla J^{\omega }\left(u\right)=\lambda u+p^{\omega },\end{array}$

其中h表示离散参数， $A_h$ 表示刚度矩阵依赖于 $\omega$， $M_h$ 是对称的质量矩阵独立于 $\omega$。

4. 数值方法

对于上述最优控制问题，常见的数值优化方法是梯度下降法，而在实际应用中，如果我们需要较高精度数值，需要更细的步长，及很小的正则化参数，这使得整个离散最优控制系统很病态，增加最优控制的数值求解花费。

下面，我们统计上采用蒙特卡洛方法(MC)，MC方法它的收敛率不受维数的影响，但是抽样成本昂贵，因为每次抽样都要求解一个最优控制系统，因此我们希望通过一些技巧减少MC抽样成本。优化算法采用共轭梯度法，分析获取控制期望所产生的成本。

对于每一个固定样本 $\omega$，离散最优控制问题的数值解我们用 $u_{k,h}^{*,\omega }$ 表示，其中k优化算法的迭代次数，因此：

$E\left[u_{k,h}^{*}\right]\approx E\left[u^{*}\right].$

我们定义

$E_N\left[u_{k,h}^{*}\right]=\frac{u_{k,h}^{*,{\omega }_1}+\cdots +u_{k,h}^{*,{\omega }_N}}{N},$

其中N表示样本量， ${\left\{u_{k,h}^{*,{\omega }_i}\right\}}_{i=1}^N$ 独立同分布， $E_N\left[u_{k,h}^{*}\right]$ 是 $E\left[u_{k,h}^{*}\right]$ 的近似。进一步。我们给出求解最优控制 $E\left[u^{*}\right]$ 的均方误差：

$\begin{array}{c}{\Vert E_N\left[u_{k,h}^{*}\right]-E\left[u^{*}\right]\Vert }_{L^2\left(\Omega ;L^2\left(D\right)\right)}^2={\Vert E_N\left[u_{k,h}^{*}\right]-E\left[u_{k,h}^{*}\right]+E\left[u_{k,h}^{*}\right]-E\left[u^{*}\right]\Vert }_{L^2\left(\Omega ;L^2\left(D\right)\right)}^2\\ ={\Vert E_N\left[u_{k,h}^{*}\right]-E\left[u_{k,h}^{*}\right]\Vert }_{L^2\left(\Omega ;L^2\left(D\right)\right)}^2+{\Vert E\left[u_{k,h}^{*}\right]-E\left[u^{*}\right]\Vert }_{L^2\left(\Omega ;L^2\left(D\right)\right)}^2.\end{array}$

这里误差分为两部分，前面部分表示MC抽样产生的统计误差，而后面部分表示有限元误差和共轭梯度迭代产生的误差。我们先估计前面部分：

${\Vert E_N\left[u_{k,h}^{*}\right]-E\left[u_{k,h}^{*}\right]\Vert }_{L^2\left(\Omega ;L^2\left(D\right)\right)}^2=\frac{1}{N}{\displaystyle {\int }_D}Var\left(u_{kh}^{*}\right)\text{d}x.$

后面部分，我们简单运用三角不等式将有限元离散和共轭梯度法误差分开：

${\Vert E\left[u_{k,h}^{*}\right]-E\left[u^{*}\right]\Vert }_{L^2\left(\Omega ;L^2\left(D\right)\right)}^2\le 2\left({\Vert u_h^{*}-u^{*}\Vert }_{L^2\left(\Omega ;L^2\left(D\right)\right)}^2+{\Vert u_{k,h}^{*}-u_h^{*}\Vert }_{L^2\left(\Omega ;L^2\left(D\right)\right)}^2\right)$

根据有限元误差估计 [1] [9]，我们知道 ${\Vert u_h^{*}-u^{*}\Vert }_{L^2\left(\Omega ;L^2\left(D\right)\right)}\le C_1{h}^2$，同时我们知道共轭梯度法是线性收敛，所以 ${\Vert u_{k,h}^{*}-u_h^{*}\Vert }_{L^2\left(\Omega ;L^2\left(D\right)\right)}\le C_2\exp \left(-k\right)$，因此我们得到求解最优控制的误差估计：

${\Vert E_N\left[u_{k,h}^{*}\right]-E\left[u^{*}\right]\Vert }_{L^2\left(\Omega ;L^2\left(D\right)\right)}\le \frac{1}{\sqrt{N}}{\left({\displaystyle {\int }_D}Var\left(u_{kh}^{*}\right)\text{d}x\right)}^{\frac{1}{2}}+C_1{h}^2+C_2\exp \left(-k\right).$

为总的误差小于 $ϵ$ ( $ϵ<1$ )，我们需要选择

$N=O\left({ϵ}^{-2}\right),h=O\left(\sqrt{ϵ}\right),k=O\left(\log \frac{1}{ϵ}\right).$

因此总的花费为 $O\left({ϵ}^{-2}\log \frac{1}{ϵ}\right)$。仔细地讲我们至少需要求解 $O\left({ϵ}^{-2}\log \frac{1}{ϵ}\right)$ 次维度为 $O\left(\sqrt{ϵ}\right)$ 的线性方程组(非线性方程组)，成本是非常高的，尤其整个离散最优控制系统比较病态时( $\lambda ,h$ 比较小时)。

为减少优化方法迭代次数，下面提出一种后验方法选取优化方法的初始迭代值，减少整个最优控制的花费。

首先在给出算法之前，我们给出一些准备，记

${\left(u,v\right)}_{M_h}=u^{\text{T}}{M}_{h}v,{\Vert u\Vert }_{M_h}=\sqrt{{\left(u,u\right)}_{M_h}}.$

第k次迭代用到的变量我们用(k)上标表示， $\nabla J^{\left(k\right)}$ 表示迭代k次计算的梯度，共轭梯度法的搜索方向

$d^{\left(k\right)}=-\nabla J^{\left(k\right)}+{\alpha }^{\left(k\right)}{d}^{\left(k-1\right)},$ 其中 $d^{\left(0\right)}=\nabla J^{\left(0\right)},$ (2)

关于 ${\alpha }^{\left(k\right)}$，我们使用戴袁公式

${\alpha }^{\left(k\right)}=\frac{{\Vert \nabla J^{\left(k\right)}\Vert }_{M_h}^2}{{\left(d^{\left(k-1\right)},\nabla J^{\left(k\right)}-\nabla J^{\left(k-1\right)}\right)}_{M_h}}.$ (3)

关于戴袁公式更多内容可参考文献 [19]。控制更新公式为 $u^{\left(k+1\right)}=u^{\left(k\right)}+s^{\left(k\right)}{d}^{\left(k\right)}$，其中 $s^{\left(k\right)}$ 表示步长。下面是共轭梯度算法：

我们使用MC估计控制期望时，每次抽样时都需要调用算法1，抽取所有的样本所用的初值往往都是一样的。这没有充分使用已抽样的样本信息。对此我们提出一种初值校正的方法。

通过少量的样本求解最优控制，然后将迭代算法的初值用平均值代替，这个平均值或许能够靠近每个样本的精确值，从而减少样本的迭代次数。

更精确的讲，根据切比雪夫不等式

$P\left({\Vert u_{kh}^{*}-E\left[u_{kh}^{*}\right]\Vert }_{L^2\left(D\right)}>ϵ\right)\le \frac{E\left[{\Vert u_{kh}^{*}-E\left[u_{kh}^{*}\right]\Vert }_{L^2\left(D\right)}^2\right]}{{ϵ}^2}=\frac{{\displaystyle {\int }_D}Var\left(u_{kh}^{*}\right)\text{d}x}{{ϵ}^2}.$

很容易地看出，当方差小时，大多数的样本的精确值是靠近平均值。使得大多数的样本的迭代次数减少。当方差比较大时，算法2的或许效果不是那么明显，可以从后面的实验看出这一点。但是这样校正初值总是比没有校正初值的效果好，因为许多情况下，样本总是在平均值附近扰动。下面我们给出了两个数值例子，表明算法2的实用性。

5. 数值实验

这部分，我们给出两个例子，并比较了初值校正(aNCG-MC)和初值未校正算法(NCG-MC)。设期望状态

$y_d\left(x\right)=\{\begin{array}{l}1x\in \left[0.25,0.75\right]\times \left[0.25,0.75\right],\\ 0其他.\end{array}$

关于随机场的参数： $a=0.3$， $n_{KL}=500$，我们数值实验在MATLAB2016a上进行(Intel(R)Core(TM)i7-6700HQ CPU@2.60GHz 2.59GHz)，关于求解线性代数方程，使用matlab中的反斜杠算子，两个算例的数值结果见图2、图3。

5.1. 算例1

这里我们考虑第3节的最优控制问题，设置空间离散步长 $h=\frac{1}{2^5}$，迭代误差 $ϵ=1e-5$，初始迭代 $u^{\left(0\right)}=1$，抽样数 $N=200$，校正样本量 $N_0=20$，设置不同的正则化参数 $\lambda =1e-4,1e-6,1e-8$，选取步长

$s^{\left(k\right)}=\frac{-{\left(A_h^{-1}{M}_h{d}^{\left(k\right)},y^{\left(k-1\right)}-y_d\right)}_{M_h}-\lambda {\left(d^{\left(k\right)},u^{\left(k-1\right)}\right)}_{M_h}}{{\Vert A_h^{-1}{M}_h{d}^{\left(k\right)}\Vert }_{M_h}^2+\lambda {\Vert d^{\left(k\right)}\Vert }_{M_h}^2},$

其中 $A_h^{-1}$ 表示刚度矩阵 $A_h$ 的逆，通过算法2求解，结果见表1，表里的数据表示抽取200个样本总的迭代数. 我们固定 $\lambda =1e-6$，设置不同的迭代误差 $ϵ=1e-3,1e-4,1e-5$，见表2的结果。

Table 1. Example 1: aNCG-MCand NCG-MC ( ϵ = 1 e − 5 )

表1. 算例1：aNCG-MC与NCG-MC ( $ϵ=1e-5$ )

从实验可以看出，初值校正是优于初值未校正，尤其是当 $\lambda$ 很小时。但当方差较大时，样本差异大时，初值校正的方法效果稍差一些，这时我们可以考虑分类的思想，将样本相似的聚在一起，减少整体样本的差异性。

Figure 2. Example 1 Numerical results ( $\lambda =1e-6,ϵ=1e-5,{\sigma }^2=0.5,h=\frac{1}{2^6}$ )

图2. 算例1数值结果( $\lambda =1e-6,ϵ=1e-5,{\sigma }^2=0.5,h=\frac{1}{2^6}$ )

Table 2. Example 1: aNCG-MCand NCG-MC ( λ = 1 e − 6 )

表2. 算例1：aNCG-MC与NCG-MC ( $\lambda =1e-6$ )

5.2. 算例2

我们考虑下面的半线性状态方程

$\begin{array}{l}-\nabla \cdot \left(k\left(x,\omega \right)\nabla y\left(x,\omega \right)\right)+f\left(y\left(x,\omega \right)\right)=u\left(x\right),\left(x,\omega \right)\in D\times \Omega ,\\ y\left(x,\omega \right)=0,\left(x,\omega \right)\in \partial D\times \Omega ,\end{array}$

其中f表示非线性反应项。设置 $f\left(x\right)=5+x^3$，关于上述方程解的存在唯一性可参考文献 [3] [20]。下面我们给出这个问题的梯度

$\begin{array}{l}\left(k^{\omega }\nabla y^{\omega },\nabla v\right)+\left(f\left(y^{\omega }\right),v\right)=\left(u,v\right),v\in H_0^1\left(D\right),\\ \left(k^{\omega }\nabla p^{\omega },\nabla v\right)+\left(f^{\prime }\left(y^{\omega }\right)p^{\omega },v\right)=\left(y^{\omega }-y_d,v\right),v\in H_0^1\left(D\right),\\ \nabla J^{\omega }\left(u\right)=\lambda u+p^{\omega },\end{array}$ (4)

其中 $f^{\prime }$ 表示f的导数。用 $\mathcal{N}$ 表示求解状态变量的算子，即 $y_{\omega }=\mathcal{N}\left(u\right)$，离散(4)可得

$\begin{array}{l}{y}^{\omega }=\mathcal{N}\left(u\right)\\ C_h{p}^{\omega }=M_h\left(y^{\omega }-y_d\right),\\ \nabla J^{\omega }\left(u\right)=\lambda u+p^{\omega },\end{array}$

其中 $C_h=A_h+B_h$， $B_h$ 表示伴随方程左端第二项的离散有限元矩阵。取迭代步长

$s^{\left(k\right)}=\frac{-{\left(\nabla J^{\left(k\right)},d^{\left(k\right)}\right)}_{M_h}}{{\left(1-p^{\left(k\right)}⊛f^{″}\left(y^{(k-1)}\right),\left(C_h^{-1}{d}^{\left(k\right)}\right)⊛\left(C_h^{-1}{d}^{\left(k\right)}\right)\right)}_{M_h}+\lambda {\Vert d^{\left(k\right)}\Vert }_{M_h}^2},$

其中 $a⊛b$ 表示 $a,b$ 对应元素相乘组成的列向量， $C_h^{-1}$ 表示 $C_h$ 的逆， $f^{″}$ 表示f的二阶导。我们采用牛顿法求解状态方程，其他的参数设置均与问题一相同，结果见表3、表4。

Figure 3. Example2 Numerical results ( $\lambda =1e-6,ϵ=1e-5,{\sigma }^2=0.5,h=\frac{1}{2^6}$ )

图3. 算例2数值结果( $\lambda =1e-6,ϵ=1e-5,{\sigma }^2=0.5,h=\frac{1}{2^6}$ )

Table 3. Example 2: aNCG-MCandNCG-MC ( ϵ = 1 e − 5 )

表3. 算例2：aNCG-MC与NCG-MC ( $ϵ=1e-5$ )

Table 4. Example 2: aNCG-MCand CG-MC ( λ = 1 e − 6 )

表4. 算例2: aNCG-MC与CG-MC ( $\lambda =1e-6$ )

6. 结论

本文通过初值校正的方法减少每次抽样的迭代次数，从而减少整个MC抽样成本。两个数值例子表明我们的方法是充分优于未初值校正方法。初值校正方法可结合其他的优化方法(原始对偶有效集、牛顿法等等)，或许有更好的效果。也可以结合多水平蒙特卡洛方法或拟蒙特卡洛方法，进一步减少计算最优控制的成本。

基金项目

项目名称：随机最优控制问题的高效Monte Carlo有限元法；合同编号：国家自然科学基金项目(11961008)。