一类随机抛物最优控制问题的集成蒙特卡洛隐式Euler法

An Ensemble Monte Carlo Implicit Euler Method for a Class of Stochastic Parabolic Optimal Control Problems

摘要: 针对扩散项带随机系数的演化方程随机最优控制问题，蒙特卡洛方法是一种非常重要的方法，针对随机量，每次抽样以后，按照以前确定问题的数值方法便可以求解，但是蒙特卡洛方法收敛速度慢，近似解误差要达到理想精度，需要大量的抽样，每次抽样需要计算一次逆矩阵，这样完成一个问题的求解，就需要计算大量的逆矩阵。我们针对一类随机抛物最优控制问题，提出一种数值求解集成蒙特卡洛隐式Euler法，即时间方向采用隐式Euler法，空间上采用线性有限元方法，随机量采用集成蒙特卡洛方法，控制量采用变分离散技巧，该方法针对所有抽样，只需要求一次逆矩阵，大大减少了计算量，通过误差分析和数值计算，该方法不影响时间方向和空间方向误差的收敛速度。

Abstract: Monte Carlo method is an important method for the stochastic optimal control problem of evolution equation with random coefficient in diffusion term. For random quantity, after each sampling, the problem can be solved using the numerical method for the corresponding deterministic problem. However, Monte Carlo method converges slowly, and the approximate solution error needs a large number of samples to achieve the desired accuracy. Each sample needs to calculate the inverse matrix. So a large number of inverse matrices need to be calculated to solve a problem. For a class of stochastic parabolic optimal control problems, we propose an ensemble Monte Carlo implicit Euler method. That is, the implicit Euler method is used for the time space, the linear finite element method is used for the space, an ensemble Monte Carlo method is used for the random coefficient, and the variational discretization technique is used for the control. This method only needs to calculate an inverse matrix once for all samples, which greatly reduces the cost. Through error analysis and numerical experiments, this method does not affect the convergence rate for time and physical space.

1. 引言

自然生活中，许多现象和工程应用都可以用微分方程来建模，自上世纪七十年代Lions的奠基性工作以来，偏微分方程及最优控制理论已经得到了广泛的研究和分析 [1] [2] [3]。实际上，在自然界中存在着大量的不确定性，比如结构力学、地震学、金融数学等领域，都存在随机的因素，因此需要用随机模型来描述这些问题。如同确定型微分方程一样，这些随机问题的准确解很难求解，用数值方法来逼近这些问题的准确解是一种必然的选择。因此，随机问题的理论分析以及数值近似成为了当前研究的热点问题。

随机最优控制问题无论在数学、金融、物理还是工程上都有广泛的应用。目前，求解随机偏微分方程最优控制问题的常用方法有蒙特卡洛(MC)方法、随机有限元方法、随机配置法等。MC方法是进行数值仿真中处理随机数据时的首选方法，此方法在求解随机偏微分方程数值解时是最简单易行的方法，它在执行过程中只需不断重复求解给定样本下的一个确定的偏微分方程。众所周知，MC方法其收敛速度为 $O\left(N^{-1/2}\right)$ (N是样本点数)，为了达到某个误差精度，该方法在求解时需要大量的抽样。尽管如此，但由于方法简单，因此该方法及其改进受到了很多学者的青睐。确定型最优控制问题的数值解，需要通过多次迭代来求解状态方程和对偶状态方程。而对于随机偏微分方程最优控制问题的数值求解，MC方法对模型中的随机量进行抽样离散时，求解每个抽样点对应的确定性偏微分方程花费的计算量比直接求解确定型最优控制问题的计算量多得多。因此，随机最优控制问题的数值求解计算量是相当惊人的。

对于椭圆随机最优控制问题已有一些研究成果 [4] [5] [6]。对于抛物随机最优控制问题，巩本学等人 [7] 采用随机Galerkin方法对控制受限的随机抛物方程最优控制问题进行了研究，得到了控制、状态、伴随变量的先验误差估计，其中随机系数的处理是通过K-L展开，该方法是将随机量展开成一个有限维数的空间，维数越高，误差越小，但是计算量随着维数的增加而呈指数增加，“维数烦恼”在此得到体现。蒙特卡洛方法是一种简单的方法，每次抽样后，只需要按照原来确定问题的方法来求解即可，只是需要求解这样的问题很多(大量)次，计算量非常大。骆艳等人在求解随机抛物偏微分方程初边值问题时采用了一种基于时间的集成算法 [8]，每个时间步迭代的矩阵是相同的，大大地提高了这种具有时间迭代的计算效率。

本文针对一类抛物随机最优控制问题，在文献 [8] 的基础上利用集成MC隐式Euler法来研究其数值近似，随机系数的处理通过定义集成平均的方式使得计算量大大减少。

2. 模型描述

设 $\left(\Omega ,\Re ,P\right)$ 是完备的概率空间， $\Omega$ 是所有可能结果的集合， $\Re$ 是事件的 $\sigma$ -代数， $P:\Re \to \left[0,1\right]$ 是一个概率测度。记 $L^2\left(D\right)$ 空间上的内积和范数分别为 $\left(\cdot ,\cdot \right)$, $\Vert \cdot \Vert$。记 $W^{s,q}\left(D\right)$ 为 $L^q\left(D\right)$ 空间上具有s阶导数的Sobolev空间，其中s是一个非负整数， $1\le q\le +\infty$。当 $q=2$ 时，常用 $H^s\left(D\right)$ 来代替 $W^{s,2}\left(D\right)$，特别地， $H_0^1\left(D\right)$ 是 $H^1\left(D\right)$ 的子空间。 $H^{-1}\left(D\right)$ 是 $H_0^1\left(D\right)$ 的对偶空间，其上的范数定义为 ${\Vert f\Vert }_{-1}={\sup }_{0\ne v\in H_0^1\left(D\right)}\left(f,v\right)/\Vert \nabla v\Vert$。对于 $1\le p,r<\infty$，空间 $L^p\left(\Omega ,W^{s,q}\left(D\right)\right)$ 包含了所有随机函数 $v:\Omega \times D\to R$，它们相对于 $\Re \otimes B\left(D\right)$ $\sigma$ -代数是可测的，其中 $B\left(D\right)$ 是D的 $\sigma$ -代数， $L^p\left(\Omega ,W^{s,q}\left(D\right)\right)$ 空间上的范数定义为

${\Vert v\Vert }_{L^p\left(\Omega ,W^{s,q}\left(D\right)\right)}={\left(E\left[{\Vert v\Vert }_{W^{s,q}\left(D\right)}^p\right]\right)}^{1/p}={\left(E{\left({\displaystyle \underset{\left|\alpha \right|\le s}{\sum }{\displaystyle {\int }_D{\left|{\partial }^{\alpha }v\right|}^q\text{d}x}}\right)}^{p/q}\right)}^{1/p}$

和

${\Vert v\Vert }_{L^{\infty }\left(\Omega ,W^{s,q}\left(D\right)\right)}=\text{ess}\underset{\Omega }{\sup }{\Vert v\Vert }_{W^{s,q}\left(D\right)}$.

空间 $L^r\left(L^p\left(\Omega ,W^{s,q}\left(D\right)\right),\left(0,T\right]\right)$ 包含了所有随机函数 $v:\Omega \times D\times \left[0,T\right]\to R$，该空间上的范数同上面定义方式一致 [7]。

考虑如下一维随机抛物方程最优控制问题：

$\underset{u\in F}{\min }J\left(y,u\right)=\frac{1}{2}E\left({\displaystyle {\int }_0^T{\Vert y-y_d\Vert }_{L^2\left(D\right)}^2}\text{d}t\right)+\frac{\alpha }{2}E\left({\displaystyle {\int }_0^T{\Vert u\Vert }_{L^2\left(D\right)}^2}\text{d}t\right)$, (1)

其中 $y,u$ 满足如下的抛物方程初边值条件：

$\{\begin{array}{l}{\partial }_{t}y-{\partial }_x\left(a\left(\omega ,x,t\right){\partial }_{x}y\left(\omega ,x,t\right)\right)=u\left(\omega ,x,t\right),\left(\omega ,x,t\right)\in \Omega \times D\times \left[0,T\right],\\ y\left(\omega ,a,t\right)=y\left(\omega ,b,t\right)=0,\left(\omega ,t\right)\in \Omega \times \left[0,T\right],\\ y\left(\omega ,x,0\right)=g^0\left(x\right),\left(\omega ,x\right)\in \Omega \times D.\end{array}$ (2)

其中 $D=\left[a,b\right]$ 表示空间区域， $\omega \in \Omega$ 为随机变量，E表示对随机变量取期望， ${\partial }_t$ 和 ${\partial }_x$ 分别表示关于时间和空间的偏导数， $y_d:D\times \left[0,T\right]\to R$ 给定， $\alpha >0$ 和 $T>0$ 分别是一个给定的常数， $F\subset L^2\left(L^2\left(D\right),\left(0,T\right]\right)$ 为凸集。

对任意 $t\in \left[0,T\right]$，为了保证问题(1)~(2)解的存在唯一性，假设系数 $a\left(\omega ,x,t\right)$ 是有界且一致强制的，即存在 $a_{\min },a_{\max }\in R$ 使得在 $\Omega \times D\times \left[0,T\right]$ 上几乎处处成立 $a_{\min }\le a\left(\omega ,x,t\right)\le a_{\max }$。

引入下列记号：

$\left({\partial }_{t}y,v\right)={\displaystyle {\int }_D{\partial }_{t}yv}\text{d}x$, $b\left(y,v\right)={\displaystyle {\int }_{D}a{\partial }_{x}y{\partial }_{x}v}\text{d}x$, $\left(u,v\right)={\displaystyle {\int }_{D}uv\text{d}x}$.

其中双线性形式 $b:H_0^1\left(D\right)\times H_0^1\left(D\right)\to R$，则状态方程(2)的弱形式可以表述为：对于给定的 $t>0,\omega \in \Omega$，求 $y\left(\omega ,\cdot ,t\right)\in H_0^1\left(D\right)$，使得

$\left({\partial }_{t}y,v\right)+b\left(y,v\right)=\left(u,v\right),\forall v\in H_0^1\left(D\right)$.

利用Lagrange乘数法可得求解原问题(1)~(2)等价于求 $y\left(\omega ,\cdot ,t\right)\in H_0^1\left(D\right)$, $u\left(\omega ,\cdot ,t\right)\in F$, $p\left(\omega ,\cdot ,t\right)\in H_0^1\left(D\right)$ 满足如下最优性系统：

$\{\begin{array}{l}\left({\partial }_{t}y,v\right)+b\left(y,v\right)=\left(u,v\right),\forall v\in H_0^1,t\in \left(0,T\right],\hfill \\ y\left(\omega ,a,t\right)=y\left(\omega ,b,t\right)=0,\left(\omega ,t\right)\in \Omega \times \left[0,T\right],\hfill \\ y\left(\omega ,x,0\right)=g^0\left(x\right),\left(\omega ,t\right)\in \Omega \times D,\hfill \\ \begin{array}{l}-\left({\partial }_{t}p,q\right)+b\left(q,p\right)=\left(y-y_d,q\right),\forall q\in H_0^1\left(D\right),t\in \left[0,T\right],\hfill \\ p\left(\omega ,a,t\right)=p\left(\omega ,b,t\right)=0,\left(\omega ,t\right)\in \Omega \times \left[0,T\right],\hfill \\ p\left(\omega ,x,T\right)=0,\left(\omega ,x,t\right)\in \Omega \times D\times \left[0,T\right],\hfill \\ {\displaystyle {\int }_0^T\left(p+\alpha u,w-u\right)\text{d}t}\ge 0,\forall w\in F.\hfill \end{array}\hfill \end{array}$ (3)

由于F为凸集，因此数值求解具有随机系数的抛物型最优控制问题(1)~(2)，就等价于求解最优性系统(3)。问题(3)的求解需要涉及到随机抛物方程的数值模拟，以下首先介绍数值模拟抛物方程定解问题的方法，然后再介绍(3)的计算方法。

3. 集成MC隐式Euler法

3.1. 随机抛物方程初边值问题的集成MC隐式Euler法

此处仅考虑随机抛物方程初边值问题的数值方法，因此此处将式(2)中的 $u\left(\omega ,x,t\right)$ 替换成 $f\left(\omega ,x,t\right)$ 并假设已知。在抛物方程初边值问题(2)的数值求解方法中，MC方法是一种非常重要的方法，在给定每个抽样(给定 $\omega$ )后，利用原来确定问题的数值方法及程序便可求解。这样为了得到式(2)的统计解，就需要将其数值求解K次，然后平均，其中第k次就相当于求解下述问题

$\{\begin{array}{l}{\partial }_t{y}_k-{\partial }_x\left(a\left({\omega }_k,x,t\right){\partial }_x{y}_k\left({\omega }_k,x,t\right)\right)=f\left({\omega }_k,x,t\right),\left({\omega }_k,x,t\right)\in \Omega \times D\times \left[0,T\right],\\ y_k\left({\omega }_k,a,t\right)=y_k\left({\omega }_k,b,t\right)=0,\left({\omega }_k,t\right)\in \Omega \times \left[0,T\right],\\ y_k\left({\omega }_k,x,0\right)=g^0\left(x\right),\left({\omega }_k,x\right)\in \Omega \times D.\end{array}$ (4)

在对问题(4)进行离散时，考虑到离散格式的稳定性，常用隐式Euler法。因此对不同的k，对应不同的系数矩阵，就要算一次矩阵的逆。集成方法得到的离散格式对所有抽样只需要算一次逆矩阵，这就大大减少了计算量。下面介绍处理问题(4)的集成MC隐式Euler法。

选取一组独立同分布的随机样本 $\left\{{\omega }_k\right\}$，对应于不同的k，有不同的 解 $y_k=y\left({\omega }_k,\cdot ,\cdot \right)$，其中 $k=1,2,\cdots ,K$。对时间区域 $\left[0,T\right]$ 进行均匀划分， $0=t_0<t_1<\cdots <t_{N-1}<t_N=T$，时间节点 $t_i=i\tau ,i=0,1,\cdots ,N$，其中步长为 $\tau =T/N$。定义随机系数的集成均值如下：

${\bar{a}}^i\left(x\right):=\bar{a}\left(x,t_i\right)=\frac{1}{K}{\displaystyle \underset{k=1}{\overset{K}{\sum }}{a}_k\left(\omega ,x,t_i\right)}=\frac{1}{K}{\displaystyle \underset{k=1}{\overset{K}{\sum }}a\left({\omega }_k,x,t_i\right)}$. (5)

则基于时间的状态方程可用如下的离散格式来近似，即

$\frac{y_k^{i+1}-y_k^i}{\tau }-{\partial }_x\left({\bar{a}}^{i+1}{\partial }_x{y}_k^{i+1}\right)-{\partial }_x\left(\left(a_k^{i+1}-{\bar{a}}^{i+1}\right){\partial }_x{y}_k^i\right)=f_k^{i+1}$, (6)

其中 $a_k^{i+1}=a\left({\omega }_k,x,t_{i+1}\right)$, $y_k^i=y\left({\omega }_k,x,t_i\right)$, $f_k^{i+1}=f\left({\omega }_k,x,t_{i+1}\right)$。将式(6)变形得

$\frac{y_k^{i+1}}{\tau }-{\partial }_x\left({\bar{a}}^{i+1}{\partial }_x{y}_k^{i+1}\right)=\frac{y_k^i}{\tau }+{\partial }_x\left(\left(a_k^{i+1}-{\bar{a}}^{i+1}\right){\partial }_x{y}_k^i\right)+f_k^{i+1}$.

对空间 $\left[a,b\right]$ 采用等距剖分 $T_h$，设步长为 $h=\left(b-a\right)/M$，节点坐标为 $x_j=a+jh$, $j=0,1,\cdots ,M$。取有限元空间 $V_h=\left\{v_h|v_h\in P_1\left(x\right),K\in T_h\right\}$，此处的K表示剖分 $T_h$ 的某个区间 $\left[x_j,x_{j+1}\right]$， $P_1\left(x\right)$ 表示一次多项式。用 $y_{k,h}^i\in V_h$ 表示 $a_k\left(\omega ,x,t_i\right)=a\left({\omega }_k,x,t_i\right)$ 时式(6)的近似解在 $t=t_i$ 时刻的值，即 $y_{k,h}^i=y_h\left({\omega }_k,x,t_i\right)$。则式(6)的集成MC隐式Euler格式为：求 $y_{k,h}^i\in V_h,i=1,\cdots ,N-\text{1}$，使得对任意 $v_h\in V_h$ 有

$\left(\frac{y_{k,h}^{i+1}-y_{k,h}^i}{\tau },v_h\right)+\left({\bar{a}}^{i+1}{\partial }_x{y}_{k,h}^{i+1},{\partial }_x{v}_h\right)+\left(\left(a_k^{i+1}-{\bar{a}}^{i+1}\right){\partial }_x{y}_{k,h}^i,{\partial }_x{v}_h\right)=\left(f_k^{i+1},v_h\right).$ (7)

其中，初始条件 $y_{k,h}^0\left({\omega }_k,x,0\right)=g_k^0\left(x\right)$。

接下来考虑其矩阵形式，在 $V_h$ 上，设 ${\phi }_j\left(x\right)$ 为节点基函数，则 $y_{k,h}^i\left({\omega }_k,x,t_i\right)$ 可表示为

${\displaystyle {\sum }_{j=1}^{M-1}{y}_{k,h}^i\left({\omega }_k,t_i\right){\phi }_j\left(x\right)}$，在式(7)中取 $v_h={\phi }_l\left(x\right)$ 便可得：

$\begin{array}{l}{\displaystyle {\int }_D{\displaystyle \underset{j=1}{\overset{M-1}{\sum }}{y}_{k,h}^{i+1}}{\phi }_j\left(x\right){\phi }_l\text{d}x}+\tau {\displaystyle {\int }_D{\displaystyle \underset{j=1}{\overset{M-1}{\sum }}{y}_{k,h}^{i+1}}{\bar{a}}^{i+1}{{\phi }^{\prime }}_j\left(x\right){{\phi }^{\prime }}_l\text{d}x}+\tau {\displaystyle {\int }_D{\displaystyle \underset{j=1}{\overset{M-1}{\sum }}{y}_{k,h}^i}\left(a_k^{i+1}-{\bar{a}}^{i+1}\right){{\phi }^{\prime }}_j\left(x\right){{\phi }^{\prime }}_l\text{d}x}\\ ={\displaystyle {\int }_D{\displaystyle \underset{j=1}{\overset{M-1}{\sum }}{y}_{k,h}^i}{\phi }_j\left(x\right){\phi }_l\text{d}x}+\tau {\displaystyle {\int }_D{f}_k^{i+1}}{\phi }_l\text{d}x\end{array}$

即

$\begin{array}{l}{\displaystyle \underset{j=1}{\overset{M-1}{\sum }}{y}_{k,h}^{i+1}}{\displaystyle {\int }_D{\phi }_j\left(x\right){\phi }_l\text{d}x}+\tau {\displaystyle \underset{j=1}{\overset{M-1}{\sum }}{y}_{k,h}^{i+1}}{\displaystyle {\int }_D{\bar{a}}^{i+1}{{\phi }^{\prime }}_j\left(x\right){{\phi }^{\prime }}_l\text{d}x}\\ =\tau {\displaystyle \underset{j=1}{\overset{M-1}{\sum }}{y}_{k,h}^i}{\displaystyle {\int }_D\left({\bar{a}}^{i+1}-a_k^{i+1}\right){{\phi }^{\prime }}_j\left(x\right){{\phi }^{\prime }}_l\text{d}x}+{\displaystyle \underset{j=1}{\overset{M-1}{\sum }}{y}_{k,h}^i}{\displaystyle {\int }_D{\phi }_j\left(x\right){\phi }_l\text{d}x}+\tau {\displaystyle {\int }_D{f}_k^{i+1}}{\phi }_l\text{d}x\end{array}$

记

$m_{j,l}={\displaystyle {\int }_D{\phi }_j\left(x\right){\phi }_l\text{d}x}$, ${\bar{a}}_{l,j}^{i+1}={\displaystyle {\int }_D{\bar{a}}^{i+1}{{\phi }^{\prime }}_j\left(x\right){{\phi }^{\prime }}_l\text{d}x}$, $a_{k,l,j}^{i+1}={\displaystyle {\int }_D\left({\bar{a}}^{i+1}-a_k^{i+1}\right){{\phi }^{\prime }}_j\left(x\right){{\phi }^{\prime }}_l\text{d}x}$, $f_{k,l}^{i+1}={\displaystyle {\int }_D{f}_k^{i+1}}{\phi }_l\text{d}x$.

则式(7)的矩阵格式为

$\left(M+\tau {\bar{A}}^{i+1}\right)Y_k^{i+1}=\left(M+\tau A_k^{i+1}\right)Y_k^i+\tau F_k^{i+1}$. (8)

为了便于比较和分析计算效果，我们同时给出一般的，没有采用集成方法MC隐式Euler格式：求 $y_{k,h}^i\in V_h,i=0,1,\cdots ,N-\text{1}$，使得对任意 $v_h\in V_h$ 有

$\left(\frac{y_{k,h}^{i+1}-y_{k,h}^i}{\tau },v_h\right)+\left(a_k^{i+1}{\partial }_x{y}_{k,h}^{i+1},{\partial }_x{v}_h\right)=\left(f_k^{i+1},v_h\right)$.

其中，初始条件 $y_{k,h}^0\left({\omega }_k,x,0\right)=g_k^0\left(x\right)$。类似的方法可得其矩阵形式为：

$\left(M+\tau A_k^{i+1}\right)Y_k^{i+1}=M{Y}_k^i+\tau F_k^{i+1}$. (9)

将式(8)和式(9)相比较，为了得到 $Y_k^{i+1}$，集成方法针对不同的抽样 ${\omega }_k$ 只需要求一次矩阵的逆 ${\left(M+\tau {\bar{A}}^{i+1}\right)}^{-1}$，而非集成方法，针对不同的抽样 ${\omega }_k$，每迭代一次都要算一次矩阵的逆 ${\left(M+\tau A_k^{i+1}\right)}^{-1}$，这样，抽样次数K越大，计算逆矩阵的次数就越多，造成计算量也越大。

3.2. 随机抛物方程初边值问题的集成MC隐式Euler格式的稳定性及误差分析

接下来先讨论集成格式(7)的稳定性，假设下述两个条件是成立的：

i) 存在一个正常数 $\theta$，使得对于任意的 $t\in \left[0,T\right]$，有

$P\left\{\omega \in \Omega ;\underset{x\in \bar{D}}{\min }a\left(\omega ,x,t\right)>\theta \right\}=1$.

ii) 存在正常数 ${\theta }_{-}$ 和 ${\theta }_{+}$，使得对于任意的 $t\in \left[0,T\right]$，有

$P\left\{{\omega }_k\in \Omega ;{\theta }_{-}\le \left|a\left({\omega }_k,x,t\right)-\bar{a}\left(x,t\right)\right|\le {\theta }_{+}\right\}=1$.

其中，条件(i)保证了问题的一致强制性，条件(ii)表示系数 $a\left({\omega }_k,x,t\right)$ 到式(5)定义的集成均值 $\bar{a}\left(x,t\right)$ 的距离是一致有界的。

引理 [8] 3.1. 设 $f_k=f\left({\omega }_k,\cdot ,\cdot \right)\in L^2\left(\left(0,T\right],H^{-1}\left(D\right)\right)$ 并且成立假设条件(i) (ii)及 $\theta >{\theta }_{+}$，则集成格式(7)是稳定的，并且式(7)的数值解满足

$\begin{array}{l}E\left[{\Vert y_{k,h}^N\Vert }^2\right]+{\theta }_{-}\tau E\left[{\Vert y_{k,h}^N\Vert }^2\right]+\left(\theta -{\theta }_{+}\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}E\left[{\Vert {\partial }_x{y}_{k,h}^i\Vert }^2\right]}\\ \le C\tau {\displaystyle \underset{i=1}{\overset{N}{\sum }}E\left[{\Vert f_k^i\Vert }_{-1}^2\right]}+C\tau E\left[{\Vert {\partial }_x{y}_{k,h}^0\Vert }^2\right]+E\left[{\Vert y_{k,h}^0\Vert }^2\right]\end{array}$

其中，C是与 $K,h,\tau$ 无关的常数。

下面对式(4)进行误差估计，设 $y_k^i$ 是在 $\omega ={\omega }_k$ 和 $t=t_i$ 时式(2)的解， $y_{k,h}^i$ 是式(7)的解，记 $y_{k,h}^i=y_h\left({\omega }_k,x,t_i\right)$，定义全离散集成MC隐式Euler格式(7)的近似为：

${\Phi }_h^i=\frac{1}{K}{\displaystyle \underset{k=1}{\overset{K}{\sum }}{y}_{k,h}^i}$.

接下来，将对 $E\left[y_k^i\right]-{\Phi }_h^i$ 进行估计，考虑到数值求解时我们采用的是MC方法进行随机抽样，因此，该式误差来源由有限元离散误差和抽样误差两部分构成，即可以写成如下

$E\left[y_k^i\right]-{\Phi }_h^i=\left(E\left[y_k^i\right]-E\left[y_{k,h}^i\right]\right)+\left(E\left[y_{k,h}^i\right]-{\Phi }_h^i\right)$.

注意到，右边第一部分为有限元离散产生的误差，第二部分是MC抽样产生的误差。

引理 [8] 3.2 对于给定的 $f_k\in L^2\left(\left(0,T\right],H^{-1}\left(D\right)\right)$ 和 $y_k^0\in L^2\left(H_0^1\left(D\right),H^2\left(D\right)\right)$，在假设(i) (ii)及 $\theta >{\theta }_{+}$ 成立的条件下，有下面式子成立

$\begin{array}{l}E\left[{\Vert E\left[y_k^N\right]-{\Phi }_h^N\Vert }^2\right]+{\theta }_{-}\tau E\left[{\Vert {\partial }_x\left(E\left[y^N\right]-{\Phi }_h^N\right)\Vert }^2\right]\\ +\left(\theta -{\theta }_{+}\right)\tau {\displaystyle \underset{i=1}{\overset{N}{\sum }}E\left[{\Vert {\partial }_x\left(E\left[y^i\right]-{\Phi }_h^i\right)\Vert }^2\right]}\\ \le \frac{1}{K}\left(C\tau {\displaystyle \underset{i=1}{\overset{N}{\sum }}E\left[{\Vert f_k^i\Vert }_{-1}^2\right]+C\tau E\left[{\Vert {\partial }_x{y}_{k,h}^0\Vert }^2\right]}+E\left[{\Vert y_{k,h}^0\Vert }^2\right]\right)+C\left({\tau }^2+h^2\right)\end{array}$

其中，C是与 $K,h,\tau$ 无关的正常数。

4. 随机抛物最优控制问题的集成MC隐式Euler法

现在我们来考虑随机抛物最优控制问题，对问题(2)中随机系数 $a\left(\omega ,x,t\right)$ 的随机量进行随机抽样 $\left\{{\omega }_k\right\}$，便得 $a_k\left(\omega ,x,t\right)=a\left({\omega }_k,x,t\right)$，如同前面取 ${\bar{a}}^i\left(x\right)$，在随机抛物方程初边值问题集成MC隐式Euler法的基础上，根据最优性条件(3)，采用变分离散的技巧 [9]，我们可以得到问题(2)的第k次集成MC隐式Euler格式为：求 $\left(y_{k,h}^i,p_{k,h}^i,u_{k,h}^i\right)\in V_h\times V_h\times F$，使得对任意 $\left(v_h,q_h\right)\in V_h\times V_h$ 有下式成立：

$\begin{array}{l}\left(\frac{y_{k,h}^{i+1}-y_{k,h}^i}{\tau },v_h\right)+\left({\bar{a}}^{i+1}{\partial }_x{y}_{k,h}^{i+1},{\partial }_x{v}_h\right)+\left(\left(a_k^{i+1}-{\bar{a}}^{i+1}\right){\partial }_x{y}_{k,h}^i,{\partial }_x{v}_h\right)=\left(f_k^{i+1},v_h\right),\\ y_{k,h}^0\left({\omega }_k,x,0\right)=g_k^0\left(x\right),i=0,1,\cdots ,N-1\end{array}$ (10)

$\begin{array}{l}-\left(\frac{p_{k,h}^{i+1}-p_{k,h}^i}{\tau },q_h\right)+\left({\bar{a}}^i{\partial }_x{p}_{k,h}^i,{\partial }_x{q}_h\right)+\left(\left(a_k^i-{\bar{a}}^i\right){\partial }_x{p}_{k,h}^i,{\partial }_x{q}_h\right)=\left(y_{k,h}^i-y_d^i,q_h\right),\\ p_{k,h}^N=0,i=N-1,\cdots ,1,0\end{array}$ (11)

$\left(p_{k,h}^i+\alpha u_{k,h}^i,w-u_{k,h}^i\right)\ge 0,\forall w\in F,i=0,1,\cdots ,N-1,N$ (12)

若 $F=L^2\left(L^2\left(D\right),\left(0,T\right]\right)$，则式(12)等价于

$u_{k,h}^i=-\frac{p_{k,h}^i}{\alpha },i=0,1,\cdots ,N-1,N$.

若 $F=U_{ad}=\left\{u|u_a\le u\left(\omega ,x,t\right)\le u_b,a.e.\text{in}\Omega \times D\times \left[0,T\right],u_a<u_b\right\}$，则式(12)等价于

$u_{k,h}^i=P_{ad}\left(-\frac{p_{k,h}^i}{\alpha }\right),i=0,1,\cdots ,N-1,N$.

其中 $P_{ad}\left(v\right)=\max \left\{u_a,\min \left\{v,u_b\right\}\right\}$。

以下考虑(10)~(12)的迭代算法，为此对任意 $v\in L^2\left(L^{\text{2}}\left(D\right),\left(0,T\right]\right)$，定义

$\left|\right|\left|v\right|\left|\right|={\displaystyle \underset{i=1}{\overset{N}{\sum }}\tau {\displaystyle \underset{j=1}{\overset{M}{\sum }}h{\left(v\left(x_j,t_i\right)\right)}^2}}$.

根据求解最优控制问题的思路，给出如下迭代算法

算法I (Algorithm I)

为了计算出目标函数的近似值，将式(1)进行离散得到：

$J\left(y,u\right)\approx \frac{1}{K}{\displaystyle \underset{k=1}{\overset{K}{\sum }}{J}_k\left(y_{k,h},u_{k,h}\right)}$,

其中 $J_k\left(y_{k,h},u_{k,h}\right)$ 的表达式为：

$J_k\left(y_{k,h},u_{k,h}\right)=\frac{1}{2}{\displaystyle \underset{i=\text{1}}{\overset{N}{\sum }}\tau {\displaystyle \underset{j=1}{\overset{M}{\sum }}h{\left(y_{k,h}^i\left(x_j\right)-y_d\left({\omega }_k,x_j,t_i\right)\right)}^2}}+\frac{\alpha }{2}{\displaystyle \underset{i=1}{\overset{N}{\sum }}\tau {\displaystyle \underset{j=1}{\overset{M}{\sum }}h}{\left(u_{k,h}^i\left(x_j\right)\right)}^2}$.

这里我们采用的是复化右矩形公式离散第k次抽样计算得到的目标函数值的近似。

5. 误差估计

考虑到数值求解(10)-(12)时我们采用的是MC方法进行随机抽样，因此，误差来源由离散误差和抽样误差两部分构成。

$\left(u,y,p\right)$ 和 $\left(u_{k,h}^i,y_{k,h}^i,p_{k,h}^i\right)$ 分别是问题(3)和(10)-(12)的解， $y\in L^2\left(L^{\text{2}}\left(\Omega ,H_0^1\left(D\right)\right),\left(0,T\right]\right)$，用样本均值 $E_K\left[y\right]$ 来近似 $E\left[y\right]$，在 $\omega ={\omega }_k,t=t_i$ 时， $E_K\left[y\right]$ 定义为

$E_K\left[y\right]=\frac{1}{K}{\displaystyle \underset{k=1}{\overset{K}{\sum }}y\left({\omega }_k,\cdot ,t_i\right)}$,

我们使用符号 ${\Vert \cdot \Vert }_V:={\Vert \cdot \Vert }_{L^2\left(L^{\text{2}}\left(\Omega ,H_0^1\left(D\right)\right),\left(0,T\right]\right)}$，则有如下误差

$\begin{array}{c}E\left[{\Vert E\left[y\right]-E_K\left[y\right]\Vert }_V^2\right]=E\left[{\Vert E\left[y\right]-\frac{1}{K}{\displaystyle \underset{k=1}{\overset{K}{\sum }}y\left({\omega }_k,\cdot ,t_i\right)}\Vert }_V^2\right]\\ =\frac{1}{K^2}E\left[{\Vert {\displaystyle \underset{k=1}{\overset{K}{\sum }}\left(E\left[y\right]-y\left({\omega }_k,\cdot ,t_i\right)\right)}\Vert }_V^2\right]\\ =\frac{1}{K^2}{\displaystyle \underset{k=1}{\overset{K}{\sum }}E}\left[{\Vert \left(E\left[y\right]-y\left({\omega }_k,\cdot ,t_i\right)\right)\Vert }_V^2\right]\\ =\frac{1}{K^2}\left(E\left[{\Vert y\Vert }_V^2\right]-{\Vert E\left[y\right]\Vert }_V^2\right)\\ \le \frac{1}{K}[\left[{\Vert y\Vert }_V^2\right]\end{array}$

对上述式子两边开平方根，即

${\Vert E\left[y\right]-E_K\left[y\right]\Vert }_V\le K^{-1/2}{\Vert y\Vert }_V$.

接下来考虑 $E\left[y\right]$ 的有限元近似 $E_K\left[y_h\right]$，用同样方式有如下定义

$E_K\left[y_h\right]=\frac{1}{K}{\displaystyle \underset{k=1}{\overset{K}{\sum }}{y}_{k,h}^i\left({\omega }_k,\cdot ,t_i\right)}$.

对于任意的 $y,y_h\in L^2\left(L^{\text{2}}\left(\Omega ,H_0^1\left(D\right)\right),\left(0,T\right]\right)$，由三角不等式可得

${\Vert E\left[y\right]-E_K\left[y_h\right]\Vert }_V\le {\Vert E\left[y\right]-E\left[y_h\right]\Vert }_V\text{+}{\Vert E\left[y_h\right]-E_K\left[y_h\right]\Vert }_V$. (13)

可以看到，式(13)右边有两部分，第一部分是有限元离散产生的估计，第二部分是MC抽样产生的估计。

首先分析第一项，根据Cauchy-Schwarz不等式可以得到

${\Vert E\left[y\right]-E\left[y_h\right]\Vert }_V={\Vert E\left[y-y_h\right]\Vert }_V\le E\left[{\Vert y-y_h\Vert }_V\right]\le C_1{\Vert y-y_h\Vert }_V$.

根据文献 [5] 有下式成立

${\Vert y-y_h\Vert }_V\le \frac{1}{\alpha }\sqrt{\alpha {\Vert y\Vert }_V^2+\max {\left(1,a_{\min }^{-1}\right)}^2{\left({\Vert y\Vert }_V+{\Vert y_d\Vert }_V\right)}^2}\left(h+\tau \right)$.

即

${\Vert E\left[y\right]-E\left[y_h\right]\Vert }_V\le Ch$.

再根据

${\Vert y_h\Vert }_V\le C_2:=\sqrt{\left(\frac{2}{\alpha a_{\min }^2}+1\right){\Vert y\Vert }_V^2+\frac{2}{\alpha }{\Vert y_d\Vert }_V^2}$,

得

${\Vert E\left[y_h\right]-E_K\left[y_h\right]\Vert }_V\le C{K}^{-1/2}$.

其中，C是一个与h无关的正常数，仅取决于 $F,\alpha ,y_d$。综上所述，根据引理3.2整理可得

${\Vert E\left[y\right]-E_K\left[y_h\right]\Vert }_V\le C\left(K^{-1/2}+h+\tau \right)$。该结论可归结为如下定理：

定理4.1 $y,y_h\in L^2\left(L^{\text{2}}\left(\Omega ,H_0^1\left(D\right)\right),\left(0,T\right]\right)$，对 $\forall K\in \mathbb{N}$，有

${\Vert E\left[y\right]-E_K\left[y_h\right]\Vert }_V\le C\left(F,\alpha ,y_d\right)\left(K^{-1/2}+h+\tau \right)$

成立，其中 $C\left(F,\alpha ,y_d\right)>0$ 是一个与h无关的正常数，仅取决于 $F,\alpha ,y_d$。

6. 数值算例

我们给出一个数值算例，通过数值模拟来验证所提方法的有效性。该算例无精确解，取空间区间为 $D=\left[0,1\right]$，时间区间为 $\left(0,1\right]$，目标泛函中的 $\alpha =1$, $y_d=t{\left(t-1\right)}^2\sin \left(2\pi x\right)$，随机系数a表示为如下形式：

$a\left(\omega ,x\right)=8+\left(1+\omega \right)\cos \left(\pi x\right)$.

用Matlab软件进行数值模拟，以下只列出 $h=\tau =1/32$ 时的模拟情况。计算时间情况如下表1，目标泛函 $J\left(\bar{y},\bar{u}\right)$ 的近似计算值见表2。

Table 1. When h = τ = 1 / 32 , the number of samples and the time of numerical simulation (s: seconds)

表1. $h=\tau =1/32$，抽样数与模拟时间情况(s：秒)

Table 2. When h = τ = 1 / 32 , the calculation results of the objective function value under different sampling

表2. $h=\tau =1/32$，不同抽样下目标函数计算结果

7. 结论

本文针对一类随机抛物最优控制问题提出的集成蒙特卡洛隐式Euler法，大大减少了计算量，抽样数越大，误差越小，越节省计算量。该方法是针对处理随机系数而提出的一种减少计算量的方法，它不影响空间方向、时间方向的误差精度。该集成格式可用于其它带随机扩散系数的演化方程随机最优控制问题的数值求解。

基金项目

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