一类Maxwell方程最优控制问题的差分格式

Difference Schemes for a Class of Maxwell Equation Optimal Control Problems

摘要: 针对一类时域Maxwell方程最优控制问题，借助Lagrange乘子法和Helmholtz分解理论，推导出由状态方程、对偶状态方程、一阶最优性条件构成的最优性系统，分别利用时域有限差分(FDTD)格式和四阶差分格式离散最优性系统，得到最优性系统的离散格式。最后，建立数值算例分别求解基于两种差分格式的Maxwell方程最优控制问题。为了避免求解大型耦合代数方程，采用迭代法进行计算，数值结果验证了理论分析结论的正确性。

Abstract: Aiming at a class of time-domain Maxwell equation optimal control problems, an optimality system composed of equation of state, dual equation of state and first-order optimality condition is derived by using Lagransssge multiplier method and Helmholtz decomposition theory. The finite-difference time-domain (FDTD) scheme and fourth-order difference scheme are used to discretization the optimality system, respectively. The discrete scheme of the optimality system is obtained. Finally, numerical examples are established to solve the Maxwell equation optimal control problem based on the two difference schemes. In order to avoid solving large coupled algebraic equations, the iterative method is used for calculation, and the numerical results verify the correctness of the theoretical analysis.

1. 引言

电磁学在许多现代应用和技术中起着重要的作用，它们包括在能源科学、纳米技术、生命科学、磁约束聚变、磁悬浮技术、微波加热、传感器技术等领域的应用。近年来，在物理或者工程等实际领域中的很多问题都可以划归为Maxwell方程约束的最优控制问题。这种复杂的电磁过程的最优控制是具有挑战性的，需要仔细的数学和数值研究。对于Maxwell方程最优控制问题，通常采用有限元法(FEM)进行求解，即建立离散的目标泛函，通过有限元方法推导出状态变量和对偶状态变量的离散格式及其最优性条件进行求解。众所周知，求解Maxwell方程的方法除了FEM之外，还有矩量法(MoM)、边界元法(BEM)、差分法等。本文主要考虑Maxwell方程最优控制问题的FDTD法、四阶差分法，其中FDTD法在时间和空间上均为二阶收敛，四阶差分法在时间上二阶收敛，空间上四阶收敛。FDTD法是由K.S.Yee [1] 在1966年首次提出的一种计算电磁场的新方法，许多研究(如 [2] - [7])表明，在简单介质中，麦克斯韦方程的高阶FDTD法比Yee格式 [8] [9] 更准确和有效。特别是高阶空间差分格式降低了色散误差和相位速度各向异性误差。

随着Maxwell方程的数值计算得到了深入的研究，截止目前，Maxwell方程最优控制问题已有一些研究成果，如抛物型涡流方程 [10] [11] [12] [13] 的最优控制，这与电磁流测量中的应用问题有关(另见 [14])。对Maxwell方程的边界可控性的早期贡献可以在 [15] [16] 中找到。我们还提到了 [17] - [23] 对线性时谐波涡流方程的最优控制 [24]，而对于非线性的情况，2013年，I.Yousept等人 [25] 对拟线性椭圆型电磁场偏微分方程的最优控制展开了研究。2017年，S.Nicaise [26] 和V.Bommer [27] 等人研究了拟线性抛物型Maxwell方程的最优控制。对于文献 [27] 中研究的时域Maxwell方程最优控制问题，在求解时域Maxwell方程最优控制问题时采用FEM进行求解。根据文献检索可知，对于Maxwell方程最优控制问题，用差分法进行数值求解的研究工作很少，大多都是基于差分法进行Maxwell方程的数值求解。因此，本文讨论时域Maxwell方程最优控制问题的差分法。

本文的结构如下：第二节给出模型问题，采用先优化后离散的方案，利用偏微分方程最优控制理论给出对偶状态方程和一阶最优性条件，得到最优性系统；第三节给出最优性系统的标量形式，利用FDTD法和四阶差分法对最优性系统建立离散格式；第四节采用迭代法计算Maxwell方程最优控制问题，并给出了基于FDTD法和四阶差分法的数值实验结果，验证了两种差分格式的有效性；第五节是总结部分。

2. 模型描述

设 $\Omega \subset {ℝ}^2$ 是一个有连通Lipschitz边界 $\Gamma =\partial \Omega$ 的有界域。此外，我们考虑了一个具有连通的Lipschitz边界 ${\Gamma }_c=\partial {\Omega }_c$，其中控制域 ${\Omega }_c\subseteq \Omega$。我们的目标是在 ${\Omega }_c$ 中找到一个最优的电流密度和随时间变化的振幅，从而驱动电场和磁场达到所需的水平。

本文考虑如下Maxwell方程最优控制问题：

$\text{Minimize}\frac{1}{2}{\Vert E\left(T\right)-E_d\Vert }_{L_{ϵ}^2\left(\Omega \right)}^2+{\Vert H\left(T\right)-H_d\Vert }_{L_{\mu }^2\left(\Omega \right)}^2+\frac{{\lambda }_u}{2}{\Vert curlu\Vert }_{L^2\left({\Omega }_c\right)}^2+\frac{{\lambda }_a}{2}{\Vert a\Vert }_{L^2\left(0,T\right)}^2,$ (1)

函数 $u$ 和a满足状态方程：

$\{\begin{array}{ll}ϵE_t-curlH+\sigma E={\chi }_{{\Omega }_c}ua\hfill & \text{in}\Omega \times \left(0,T\right),\hfill \\ \mu H_t+curlE=0\hfill & \text{in}\Omega \times \left(0,T\right),\hfill \\ E\times n=0\hfill & \text{in}\Gamma \times \left(0,T\right),\hfill \\ E\left(·,0\right)=E^0\hfill & \text{in}\Omega ,\hfill \\ H\left(·,0\right)=H^0\hfill & \text{in}\Omega ,\hfill \end{array}$ (2)

对电流密度 $u=u\left(x\right)$ 的无散度约束(电荷守恒定律)以及对随时间变化的振幅 $a=a\left(t\right)$ 的约束为：

$\{\begin{array}{l}divu=0\text{in}{\Omega }_c,\\ a_{\min }\le a\left(t\right)\le a_{\max }\text{a}\text{.e}\text{.}\text{in}\left(0,T\right),\end{array}$ (3)

以上描述的是一个受Maxwell方程约束的最优控制问题，假设Maxwell方程的解满足完美导体(PEC)边界条件。在该问题中， $E$ 表示电场， $H$ 表示磁场， $ϵ$ 、 $\mu$ 、 $\sigma$ ： $\Omega \to {ℝ}^{2\times 2}$ 分别表示介电常数、磁导率、电导率， $n$ 为边界 $\Gamma$ 的外法向量， $n_c$ 为边界 ${\Gamma }_c$ 的外法向量，向量场 $E_d$ 、 $H_d$ 分别表示最终时刻 $T\in {ℝ}^{+}$ 时所期望的电场和磁场，两者均为给定函数。常数 $a_{\min }\le a_{\max }$ 表示过程中所允许振幅的最小值和最大值， ${\lambda }_u$ 、 ${\lambda }_a\in {ℝ}^{+}$ 为控制成本参数。目标泛函(1)中加权范数 ${\Vert \cdot \Vert }_{L_{ϵ}^2\left(\Omega \right)}^2={\left(ϵ\cdot ,\cdot \right)}_{L^2\left(\Omega \right)}$， ${\Vert \cdot \Vert }_{L_{\mu }^2\left(\Omega \right)}^2={\left(\mu \cdot ,\cdot \right)}_{L^2\left(\Omega \right)}$。

引理1.1 [17]最优控制问题(1)存在一个最优解。

接下来引入集合：

$H\left(div;{\Omega }_c\right)=\left\{q\in L^2\left({\Omega }_c\right)|divq\in L^2\left({\Omega }_c\right)\right\},$

$H_0\left(curl;{\Omega }_c\right)=\left\{q\in H\left(curl;{\Omega }_c\right)|q\times n_c=0\text{on}{\Gamma }_c\right\},$

$H\left(div=0;{\Omega }_c\right)=\left\{q\in H\left(div;{\Omega }_c\right)|divq=0\text{in}{\Omega }_c\right\},$

${\mathcal{U}}^{feas}=H_0\left(curl;{\Omega }_c\right)\cap H\left(div=0;{\Omega }_c\right),$

${\mathcal{A}}^{feas}=\left\{a\in L^2\left(0,T\right)|a_{\min }\le a\left(t\right)\le a_{\max }\text{a}\text{.e}\text{.}\text{in}\left(0,T\right)\right\},$

利用Lagrange乘子法，以及Helmholtz分解理论，有以下结论成立：

定理1.2令 $\left(u,a\right)\in {\mathcal{U}}^{feas}\times {\mathcal{A}}^{feas}$，问题(1)~(3)存在唯一的最优解 $\left\{E,H,u,a\right\}$，当且仅当存在对偶状态变量 $\left\{K,Q\right\}$，使得 $\left\{E,H,K,Q,u,a\right\}$ 满足如下的最优性系统：

$\{\begin{array}{ll}ϵE_t-curlH+\sigma E={\chi }_{{\Omega }_c}ua\hfill & \text{in}\Omega \times \left(0,T\right),\hfill \\ \mu H_t+curlE=0\hfill & \text{in}\Omega \times \left(0,T\right),\hfill \\ E\times n=0\hfill & \text{in}\Gamma \times \left(0,T\right),\hfill \\ E\left(·,0\right)=E_0\hfill & \text{in}\Omega ,\hfill \\ H\left(·,0\right)=H_0\hfill & \text{in}\Omega ,\hfill \end{array}$ (4)

$\{\begin{array}{ll}ϵK_t-curlQ-\sigma K=0\hfill & \text{in}\Omega \times \left(0,T\right),\hfill \\ \mu Q_t+curlK=0\hfill & \text{in}\Omega \times \left(0,T\right),\hfill \\ K\times n=0\hfill & \text{in}\Gamma \times \left(0,T\right),\hfill \\ K\left(T\right)=E\left(T\right)-E_d\hfill & \text{in}\Omega ,\hfill \\ Q\left(T\right)=H\left(T\right)-H_d\hfill & \text{in}\Omega ,\hfill \end{array}$ (5)

$\{\begin{array}{ll}curlcurlu+\nabla \phi =-{\lambda }_u^{-1}{\left(K,a\right)}_{L^2(0,T)}\hfill & \text{in}{\Omega }_c,\hfill \\ divu=0\hfill & \text{in}{\Omega }_c,\hfill \\ u\times n_c=0\hfill & \text{on}{\Gamma }_c\hfill \end{array}$ (6)

$a\left(t\right)={\mathbb{P}}_{\left[a_{\min },a_{\max }\right]}\left(-\frac{1}{{\lambda }_a}{\displaystyle {\int }_{{\Omega }_c}K\left(x,y,t\right)\cdot u\left(x,y\right)\text{d}x\text{d}y}\right),\forall t\in \left[0,T\right].$ (7)

证明：定义Lagrange函数

$\begin{array}{l}L\left(E,H,u,a,K,Q\right)\\ =\frac{1}{2}{\Vert E\left(T\right)-E_d\Vert }_{L_{ϵ}^2\left(\Omega \right)}^2+{\Vert H\left(T\right)-H_d\Vert }_{L_{\mu }^2\left(\Omega \right)}^2+\frac{{\lambda }_u}{2}{\Vert curlu\Vert }_{L^2\left({\Omega }_c\right)}^2+\frac{{\lambda }_a}{2}{\Vert a\Vert }_{L^2\left(0,T\right)}^2\\ -{\displaystyle {\iint }_{\Omega \times \left(0,T\right)}\left(ϵE_t-curlH+\sigma E-{\chi }_{{\Omega }_c}ua\right)\cdot K\text{d}x\text{d}y\text{d}t}-{\displaystyle {\iint }_{\Omega \times \left(0,T\right)}\left(\mu H_t+curlE\right)\cdot Q\text{d}x\text{d}y\text{d}t}\\ =\frac{1}{2}{\Vert E\left(T\right)-E_d\Vert }_{L_{ϵ}^2\left(\Omega \right)}^2+{\Vert H\left(T\right)-H_d\Vert }_{L_{\mu }^2\left(\Omega \right)}^2+\frac{{\lambda }_u}{2}{\Vert curlu\Vert }_{L^2\left({\Omega }_c\right)}^2+\frac{{\lambda }_a}{2}{\Vert a\Vert }_{L^2\left(0,T\right)}^2\\ -{\displaystyle {\iint }_{\Omega \times \left(0,T\right)}ϵE_t\cdot K\text{d}x\text{d}y\text{d}t}-{\displaystyle {\iint }_{\Omega \times \left(0,T\right)}\sigma E\cdot K\text{d}x\text{d}y\text{d}t}-{\displaystyle {\iint }_{\Omega \times \left(0,T\right)}\left(-curlH-{\chi }_{{\Omega }_c}ua\right)\cdot K\text{d}x\text{d}y\text{d}t}\end{array}$

$\begin{array}{l}-{\displaystyle {\iint }_{\Omega \times \left(0,T\right)}curlE\cdot Q\text{d}x\text{d}y\text{d}t}-{\displaystyle {\iint }_{\Omega \times \left(0,T\right)}\mu H_t\cdot Q\text{d}x\text{d}y\text{d}t}\\ =\frac{1}{2}{\Vert E\left(T\right)-E_d\Vert }_{L_{ϵ}^2\left(\Omega \right)}^2+{\Vert H\left(T\right)-H_d\Vert }_{L_{\mu }^2\left(\Omega \right)}^2+\frac{{\lambda }_u}{2}{\Vert curlu\Vert }_{L^2\left({\Omega }_c\right)}^2+\frac{{\lambda }_a}{2}{\Vert a\Vert }_{L^2\left(0,T\right)}^2\\ -{\displaystyle {\int }_{\Omega }ϵE\left(T\right)\cdot K\left(T\right)\text{d}x\text{d}y}+{\displaystyle {\iint }_{\Omega \times \left(0,T\right)}ϵE\cdot K_t\text{d}x\text{d}y\text{d}t}-{\displaystyle {\iint }_{\Omega \times \left(0,T\right)}\sigma E\cdot K\text{d}x\text{d}y\text{d}t}\\ -{\displaystyle {\iint }_{\Omega \times \left(0,T\right)}\left(-curlH-{\chi }_{{\Omega }_c}ua\right)\cdot K\text{d}x\text{d}y\text{d}t}-{\displaystyle {\iint }_{\Omega \times \left(0,T\right)}curlQ\cdot E\text{d}x\text{d}y\text{d}t}\\ -{\displaystyle {\iint }_{\Gamma \times \left(0,T\right)}E\times n\cdot Q\text{d}x\text{d}y\text{d}t}-{\displaystyle {\iint }_{\Omega \times \left(0,T\right)}\mu H_t\cdot Q\text{d}x\text{d}y\text{d}t},\end{array}$

最后一个等式由Green公式得到，且由边界条件可知，最后一个等式的最后一项为零。

$\begin{array}{l}{D}_{E}L\left(E,H,u,a,K,Q\right)E\\ =ϵ{\displaystyle {\int }_{\Omega }\left(E\left(T\right)-E_d\right)\cdot E\left(T\right)\text{d}x\text{d}y}-ϵ{\displaystyle {\int }_{\Omega }E\left(T\right)\cdot K\left(T\right)\text{d}x\text{d}y}-{\displaystyle {\iint }_{\Omega \times \left(0,T\right)}ϵK_t\cdot E\text{d}x\text{d}y\text{d}t}\\ -{\displaystyle {\iint }_{\Omega \times \left(0,T\right)}curlQ\cdot E\text{d}x\text{d}y\text{d}t}-{\displaystyle {\iint }_{\Omega \times \left(0,T\right)}\sigma K\cdot E\text{d}x\text{d}y\text{d}t}\\ ={\displaystyle {\int }_{\Omega }\left[\left(E\left(T\right)-E_d\right)-K\left(T\right)\right]\cdot ϵE\left(T\right)\text{d}x\text{d}y}-{\displaystyle {\iint }_{\Omega \times \left(0,T\right)}\left(ϵK_t-curlQ-\sigma K\right)\cdot E\text{d}x\text{d}y\text{d}t}\\ =0,\end{array}$

则有

$K\left(T\right)=E\left(T\right)-E_d\text{in}\Omega ,$

$ϵK_t-curlQ-\sigma K=0\text{in}\Omega \times \left(0,T\right).$

类似地，

$\begin{array}{l}{D}_{H}L\left(E,H,u,a,K,Q\right)H\\ ={\displaystyle {\int }_{\Omega }\left[\left(H\left(T\right)-H_d\right)-Q\left(T\right)\right]\cdot \mu H\left(T\right)\text{d}x\text{d}y}+{\displaystyle {\iint }_{\Omega \times \left(0,T\right)}\left(\mu Q_t+curlK\right)\cdot H\text{d}x\text{d}y\text{d}t}\\ +{\displaystyle {\iint }_{\Gamma \times \left(0,T\right)}K\times n\cdot H\text{d}x\text{d}y\text{d}t}\\ =0,\end{array}$

则有

$Q\left(T\right)=H\left(T\right)-H_d\text{in}\Omega ,$

$\mu Q_t+curlK=0\text{in}\Omega \times \left(0,T\right),$

$K\times n=0\text{in}\Gamma \times \left(0,T\right).$

(6)、(7)式由Lagrange函数对控制 $u$ 和a求导，以及Helmholtz理论得到，具体证明见文献 [27]。

证毕。

3. 最优性系统的离散格式

3.1. 最优性系统

在二维直角坐标系下，设所有的物理量均与z轴坐标无关，则定理1.2中最优性系统(4)~(7)式可以表示成如下标量形式：

$\{\begin{array}{ll}ϵ\frac{\partial E_x}{\partial t}+\sigma E_x=\frac{\partial H}{\partial y}+{\chi }_{{\Omega }_c}{u}_{x}a,\hfill & \left(x,y,t\right)\in \Omega \times \left(0,T\right),\hfill \\ ϵ\frac{\partial E_y}{\partial t}+\sigma E_y=-\frac{\partial H}{\partial x}+{\chi }_{{\Omega }_c}{u}_{y}a,\hfill & \left(x,y,t\right)\in \Omega \times \left(0,T\right),\hfill \\ \mu \frac{\partial H}{\partial t}=\frac{\partial E_x}{\partial y}-\frac{\partial E_y}{\partial x},\hfill & \left(x,y,t\right)\in \Omega \times \left(0,T\right),\hfill \\ E_x\left(x,c,t\right)=E_x\left(x,d,t\right)=E_y\left(a,y,t\right)=E_y\left(b,y,t\right)=0,\hfill & \left(x,y,t\right)\in \Gamma \times \left(0,T\right),\hfill \\ E_x\left(x,y,0\right)=E_x^0\left(x,y\right),E_y\left(x,y,0\right)=E_y^0\left(x,y\right),H\left(x,y,0\right)=H^0\left(x,y\right),\hfill & \left(x,y\right)\in \Omega ,\hfill \end{array}$ (8)

$\{\begin{array}{ll}ϵ\frac{\partial K_x}{\partial t}-\sigma K_x=\frac{\partial Q}{\partial y},\hfill & \left(x,y,t\right)\in \Omega \times \left(0,T\right),\hfill \\ ϵ\frac{\partial K_y}{\partial t}-\sigma K_y=-\frac{\partial Q}{\partial x},\hfill & \left(x,y,t\right)\in \Omega \times \left(0,T\right),\hfill \\ \mu \frac{\partial Q}{\partial t}=\frac{\partial K_x}{\partial y}-\frac{\partial K_y}{\partial x},\hfill & \left(x,y,t\right)\in \Omega \times \left(0,T\right),\hfill \\ K_x\left(x,c,t\right)=K_x\left(x,d,t\right)=K_y\left(a,y,t\right)=K_y\left(b,y,t\right)=0,\hfill & \left(x,y,t\right)\in \Gamma \times \left(0,T\right),\hfill \\ K_x\left(x,y,T\right)=E_x^T\left(x,y\right)-E_d\left(x,y\right),K_y\left(x,y,T\right)=E_y^T\left(x,y\right)-E_d\left(x,y\right),\hfill & \left(x,y\right)\in \Omega ,\hfill \\ Q\left(x,y,T\right)=H^T\left(x,y\right)-H_d\left(x,y\right),\hfill & \left(x,y\right)\in \Omega ,\hfill \end{array}$ (9)

$\{\begin{array}{ll}\Delta u_x+\frac{\partial \phi }{\partial x}=-{\lambda }_u^{-1}{\left(K_x,a\right)}_{L^2\left(0,T\right)},\hfill & \left(x,y\right)\in {\Omega }_c,\hfill \\ \Delta u_y+\frac{\partial \phi }{\partial y}=-{\lambda }_u{}^{-1}{\left(K_y,a\right)}_{L^2\left(0,T\right)},\hfill & \left(x,y\right)\in {\Omega }_c,\hfill \\ \frac{\partial \phi }{\partial x}+\frac{\partial \phi }{\partial y}=0,\hfill & \left(x,y\right)\in {\Omega }_c,\hfill \\ u_x\left(x,c\right)=u_x\left(x,d\right)=u_y\left(a,y\right)=u_y\left(b,y\right)=0,\hfill & (x,y)\in {\Gamma }_c,\hfill \end{array}$ (10)

$a\left(t\right)={\mathbb{P}}_{\left[a_{\min },a_{\max }\right]}\left(-\frac{1}{{\lambda }_a}{\displaystyle {\int }_{{\Omega }_c}{K}_x\left(x,y,t\right)u_x\left(x,y\right)+K_y\left(x,y,t\right)u_y\left(x,y\right)\text{d}x\text{d}y}\right),\forall t\in \left[0,T\right],$

这里 $\left(E_x\left(x,y,t\right),E_y\left(x,y,t\right)\right)=E$ 是二维电场强度， $H\left(x,y,t\right)\equiv H_z\left(x,y,t\right)$ 表示z方向的磁场强度，电流密度 $u=\left(u_x\left(x,y\right),u_y\left(x,y\right)\right)$。 $T>0$ 为总时间，取 $\Omega =\left[a,b\right]\times \left[c,d\right]$，这里的a和b为正常数。初始条件中 $E_x^0\left(x,y\right)$， $E_y^0\left(x,y\right)$， $H^0\left(x,y\right)$ 为给定函数。

对区域 $\Omega$ 进行均匀剖分， $a=x_0<x_1<\cdots <x_{N_x}=b$， $c=y_0<y_1<\cdots <y_{N_y}=d$，对时间区域 $\left[0,T\right]$ 剖分为 $N_t$ 个网格，则离散的时间 $t_k=k\Delta t$， $\Delta t=\frac{T}{N_t}$， $k=0,1,\cdots ,N_t$。x方向网格点 $x_i=i\Delta x$， $\Delta x=\frac{b-a}{N_x}$， $i=0,1,\cdots ,N_x$，y方向网格点 $y_j=j\Delta y$， $\Delta y=\frac{d-c}{N_y}$， $j=0,1,\cdots ,N_y$。这里的Δx和Δy可取不同的值。

在接下来的两小节中讨论了最优性系统的FDTD格式和四阶差分格式，即对状态方程和对偶状态方程分别采用FDTD格式和四阶差分格式，对最优性系统中的静磁鞍点方程(10)采用MAC格式，MAC格式类似FDTD格式在空间上的离散。对于未知函数 $E_x$ 、 $K_x$ 和 $u_x$ 在x方向的中点取值， $E_y$ 、 $K_y$ 和 $u_y$ 在y方向的中点取值，H、Q和 $\phi$ 在x和y方向的中点取值， $E_x$ 、 $E_y$ 、 $K_x$ 和 $K_y$ 在时间整节点上取值，H和Q在时间半节点上取值。

3.2. 最优性系统的FDTD格式

记 $v_{i,j}^n=v\left(x_i,y_j,t_n\right)$，则

${\delta }_x{v}_{i,j}^n=\frac{v_{i+\frac{1}{2},j}^n-v_{i-\frac{1}{2},j}^n}{\Delta x},{\delta }_y{v}_{i,j}^n=\frac{v_{i,j+\frac{1}{2}}^n-v_{i,j-\frac{1}{2}}^n}{\Delta y},$

${\delta }_t{v}_{i,j}^n=\frac{v_{i,j}^{n+\frac{1}{2}}-v_{i,j}^{n-\frac{1}{2}}}{\Delta t},{\bar{\delta }}_t{v}_{i,j}^n=\frac{v_{i,j}^{n+\frac{1}{2}}+v_{i,j}^{n-\frac{1}{2}}}{\Delta t},$

因此最优性系统的离散格式为：

$\{\begin{array}{l}ϵ\frac{E_{x_{i+\frac{1}{2},j}}^{n+1}-E_{x_{i+\frac{1}{2},j}}^n}{\Delta t}+\sigma \frac{E_{x_{i+\frac{1}{2},j}}^{n+1}+E_{x_{i+\frac{1}{2},j}}^n}{2}={\delta }_y{H}_{i+\frac{1}{2},j}^{n+\frac{1}{2}}+{\chi }_{{\Omega }_c}{u}_x{}_{i+\frac{1}{2},j}{a}^n,\\ ϵ\frac{E_{y_{i,j+\frac{1}{2}}}^{n+1}-E_{y_{i,j+\frac{1}{2}}}^n}{\Delta t}+\sigma \frac{E_{y_{i,j+\frac{1}{2}}}^{n+1}+E_{y_{ij+\frac{1}{2}}}^n}{2}=-{\delta }_x{H}_{i,j+\frac{1}{2}}^{n+\frac{1}{2}}+{\chi }_{{\Omega }_c}{u}_y{}_{i,j+\frac{1}{2}}{a}^n,\\ \mu \frac{H_{i+\frac{1}{2},j+\frac{1}{2}}^{n+\frac{3}{2}}-H_{i+\frac{1}{2},j+\frac{1}{2}}^{n+\frac{1}{2}}}{\Delta t}=-{\delta }_x{E}_{y_{i+\frac{1}{2},j+\frac{1}{2}}}^{n+1}+{\delta }_y{E}_{x_{i+\frac{1}{2},j+\frac{1}{2}}}^{n+1},\\ E_{x_{i+\frac{1}{2},0}}^n=E_{x_{i+\frac{1}{2},N_y}}^n=E_{y_{0,j+\frac{1}{2}}}^n=E_{y_{N_x,j+\frac{1}{2}}}^n=0,\\ E_{x_{i+\frac{1}{2},j}}^0=E_x^0\left(x_{i+\frac{1}{2}},y_j\right),E_{y_{i,j+\frac{1}{2}}}^0=E_y^0\left(x_i,y_{j+\frac{1}{2}}\right),H_{i+\frac{1}{2},j+\frac{1}{2}}^{\frac{1}{2}}=H^0\left(x_{i+\frac{1}{2}},y_{j+\frac{1}{2}}\right),\end{array}$ (12)

$\{\begin{array}{l}ϵ\frac{K_{x_{i+\frac{1}{2},j}}^{n+1}-K_{x_{i+\frac{1}{2},j}}^n}{\Delta t}-\sigma \frac{K_{x_{i+\frac{1}{2},j}}^{n+1}+K_{x_{i+\frac{1}{2},j}}^n}{2}={\delta }_y{Q}_{i+\frac{1}{2},j}^{n+\frac{1}{2}},\\ ϵ\frac{K_{y_{i,j+\frac{1}{2}}}^{n+1}-K_{y_{i,j+\frac{1}{2}}}^n}{\Delta t}-\sigma \frac{K_{y_{i,j+\frac{1}{2}}}^{n+1}+K_{y_{ij+\frac{1}{2}}}^n}{2}=-{\delta }_x{Q}_{i,j+\frac{1}{2}}^{n+\frac{1}{2}},\\ \mu \frac{Q_{i+\frac{1}{2},j+\frac{1}{2}}^{n+\frac{3}{2}}-Q_{i+\frac{1}{2},j+\frac{1}{2}}^{n+\frac{1}{2}}}{\Delta t}=-{\delta }_x{K}_{y_{i+\frac{1}{2},j+\frac{1}{2}}}^{n+1}+{\delta }_y{K}_{x_{i+\frac{1}{2},j+\frac{1}{2}}}^{n+1},\\ K_{x_{i+\frac{1}{2},0}}^n=K_{x_{i+\frac{1}{2},N_y}}^n=K_{y_{0,j+\frac{1}{2}}}^n=K_{y_{N_x,j+\frac{1}{2}}}^n=0,\\ K_{x_{i+\frac{1}{2},j}}^T=E_x^T\left(x_{i+\frac{1}{2}},y_j\right)-E_{x_d}\left(x_{i+\frac{1}{2}},y_j\right),k_{y_{i,j+\frac{1}{2}}}^T=E_y^T\left(x_i,y_{j+\frac{1}{2}}\right)-E_{y_d}\left(x_i,y_{j+\frac{1}{2}}\right),\\ Q_{i+\frac{1}{2},j+\frac{1}{2}}^T=H^T\left(x_{i+\frac{1}{2}},y_{j+\frac{1}{2}}\right)-H_d\left(x_{i+\frac{1}{2}},y_{j+\frac{1}{2}}\right),\end{array}$ (13)

$\{\begin{array}{l}\frac{4u_{x_{i+\frac{1}{2},j}}-u_{x_{i-\frac{1}{2},j}}-u_{x_{i+\frac{3}{2},j}}-u_{x_{i+\frac{1}{2},j-1}}-u_{x_{i+\frac{1}{2},j+1}}}{\Delta x\Delta y}+\frac{{\phi }_{i+\frac{1}{2},j+\frac{1}{2}}-{\phi }_{i-\frac{1}{2},j+\frac{1}{2}}}{\Delta x}=-\frac{1}{{\lambda }_u}{\left(K_{x_{i+\frac{1}{2},j}}^n,a^n\right)}_{L^2\left(0,T\right)},\\ \frac{4u_{y_{i,j+\frac{1}{2}}}-u_{y_{i-1,j+\frac{1}{2}}}-u_{y_{i+1,j+\frac{1}{2}}}-u_{y_{i,j-\frac{1}{2}}}-u_{y_{i,j+\frac{3}{2}}}}{\Delta x\Delta y}+\frac{{\phi }_{i+\frac{1}{2},j+\frac{1}{2}}-{\phi }_{i+\frac{1}{2},j-\frac{1}{2}}}{\Delta y}=-\frac{1}{{\lambda }_u}{\left(K_{y_{i,j+\frac{1}{2}}}^n,a^n\right)}_{L^2\left(0,T\right)},\\ \frac{u_{x_{i+\frac{1}{2},j}}-u_{x_{i+\frac{3}{2},j}}}{\Delta x}+\frac{u_{y_{i,j+\frac{1}{2}}}-u_{y_{i,j-\frac{1}{2}}}}{\Delta y}=0,\\ u_{x_{i+\frac{1}{2},0}}=u_{x_{i+\frac{1}{2},N_y}}=u_{y_{0,j+\frac{1}{2}}}=u_{y_{N_x,j+\frac{1}{2}}}=0,\end{array}$ (14)

$a^n={\mathbb{P}}_{\left[a_{\min },a_{\max }\right]}\left(\frac{1}{{\lambda }_a}{\displaystyle {\int }_{{\Omega }_c}{K}_{x_{i+\frac{1}{2},j}}^n{u}_{x_{i+\frac{1}{2},j}}+K_{y_{i,j+\frac{1}{2}}}^n{u}_{y_{i,j+\frac{1}{2}}}\text{d}x\text{d}y}\right).$ (15)

3.3. 最优性系统的四阶差分格式

在Maxwell方程中，我们使用四阶差分格式来近似这些偏导数 $\frac{\partial }{\partial x}$ ：

${\delta }_x^{\left(4\right)}{v}_{i,j}=\frac{27\left(v_{i+\frac{1}{2},j}-v_{i-\frac{1}{2},j}\right)-\left(v_{i+\frac{3}{2},j}-v_{i-\frac{3}{2},j}\right)}{24}={\frac{\partial v}{\partial x}|}_{i,j}+{\rm O}\left(h_x^4\right),$

这里 $v_{i,j}$ 表示v在 $\left(x_i,y_j\right)$ 处的近似，在计算靠近边界的节点时，使用单边差分近似：

${\delta }_x^{\left(4\right)}{v}_{1,j}=\frac{-22v_{\frac{1}{2},j}+17v_{\frac{3}{2},j}+9v_{\frac{5}{2},j}-5v_{\frac{7}{2},j}+v_{\frac{9}{2},j}}{24},$

${\delta }_x^{\left(4\right)}{v}_{N_x-1,j}=\frac{-v_{N_x-\frac{9}{2},j}+5v_{N_x-\frac{7}{2},j}-9v_{N_x-\frac{5}{2},j}-17v_{N_x-\frac{3}{2},j}+22v_{N_x-\frac{1}{2},j}}{24},$

相似地可由 ${\delta }_y^{\left(4\right)}$ 近似 $\frac{\partial }{\partial y}$，这里算子可以移半个网格点，即在上述公式中，i和j可以移 $\frac{1}{2}$。时间导数项仍采用二阶差分，最后得到最优性系统四阶差分的离散格式：

$\{\begin{array}{l}ϵ\frac{E_{x_{i+\frac{1}{2},j}}^{n+1}-E_{x_{i+\frac{1}{2},j}}^n}{\Delta t}+\sigma \frac{E_{x_{i+\frac{1}{2},j}}^{n+1}+E_{x_{i+\frac{1}{2},j}}^n}{2}=\frac{1}{h_y}{\delta }_y^{\left(4\right)}{H}_{i+\frac{1}{2},j}^{n+\frac{1}{2}}+{\chi }_{{\Omega }_c}{u}_x{}_{i+\frac{1}{2},j}{a}^n,\\ ϵ\frac{E_{y_{i,j+\frac{1}{2}}}^{n+1}-E_{y_{i,j+\frac{1}{2}}}^n}{\Delta t}+\sigma \frac{E_{y_{i,j+\frac{1}{2}}}^{n+1}+E_{y_{ij+\frac{1}{2}}}^n}{2}=-\frac{1}{h_x}{\delta }_x^{\left(4\right)}{H}_{i,j+\frac{1}{2}}^{n+\frac{1}{2}}+{\chi }_{{\Omega }_c}{u}_y{}_{i,j+\frac{1}{2}}{a}^n,\\ \mu \frac{H_{i+\frac{1}{2},j+\frac{1}{2}}^{n+\frac{3}{2}}-H_{i+\frac{1}{2},j+\frac{1}{2}}^{n+\frac{1}{2}}}{\Delta t}=-\frac{1}{h_x}{\delta }_x^{\left(4\right)}{E}_{y_{i+\frac{1}{2},j+\frac{1}{2}}}^{n+1}+\frac{1}{h_y}{\delta }_y^{\left(4\right)}{E}_{x_{i+\frac{1}{2},j+\frac{1}{2}}}^{n+1},\\ E_{x_{i+\frac{1}{2},0}}^n=E_{x_{i+\frac{1}{2},N_y}}^n=E_{y_{0,j+\frac{1}{2}}}^n=E_{y_{N_x,j+\frac{1}{2}}}^n=0,\\ E_{x_{i+\frac{1}{2},j}}^0=E_x^0\left(x_{i+\frac{1}{2}},y_j\right),E_{y_{i,j+\frac{1}{2}}}^0=E_y^0\left(x_i,y_{j+\frac{1}{2}}\right),H_{i+\frac{1}{2},j+\frac{1}{2}}^{\frac{1}{2}}=H^0\left(x_{i+\frac{1}{2}},y_{j+\frac{1}{2}}\right),\end{array}$ (16)

$\{\begin{array}{l}ϵ\frac{K_{x_{i+\frac{1}{2},j}}^{n+1}-K_{x_{i+\frac{1}{2},j}}^n}{\Delta t}-\sigma \frac{K_{x_{i+\frac{1}{2},j}}^{n+1}+K_{x_{i+\frac{1}{2},j}}^n}{2}=\frac{1}{h_y}{\delta }_y^{\left(4\right)}{Q}_{i+\frac{1}{2},j}^{n+\frac{1}{2}},\\ ϵ\frac{K_{y_{i,j+\frac{1}{2}}}^{n+1}-K_{y_{i,j+\frac{1}{2}}}^n}{\Delta t}-\sigma \frac{K_{y_{i,j+\frac{1}{2}}}^{n+1}+K_{y_{ij+\frac{1}{2}}}^n}{2}=-\frac{1}{h_x}{\delta }_x^{\left(4\right)}{Q}_{i,j+\frac{1}{2}}^{n+\frac{1}{2}},\\ \mu \frac{Q_{i+\frac{1}{2},j+\frac{1}{2}}^{n+\frac{3}{2}}-Q_{i+\frac{1}{2},j+\frac{1}{2}}^{n+\frac{1}{2}}}{\Delta t}=-\frac{1}{h_x}{\delta }_x^{\left(4\right)}{K}_{y_{i+\frac{1}{2},j+\frac{1}{2}}}^{n+1}+\frac{1}{h_y}{\delta }_y^{\left(4\right)}{K}_{x_{i+\frac{1}{2},j+\frac{1}{2}}}^{n+1},\\ K_{x_{i+\frac{1}{2},0}}^n=K_{x_{i+\frac{1}{2},N_y}}^n=K_{y_{0,j+\frac{1}{2}}}^n=K_{y_{N_x,j+\frac{1}{2}}}^n=0,\\ K_{x_{i+\frac{1}{2},j}}^0=K_x^0\left(x_{i+\frac{1}{2}},y_j\right),K_{y_{i,j+\frac{1}{2}}}^0=K_y^0\left(x_i,y_{j+\frac{1}{2}}\right),Q_{i+\frac{1}{2},j+\frac{1}{2}}^{\frac{1}{2}}=Q^0\left(x_{i+\frac{1}{2}},y_{j+\frac{1}{2}}\right),\end{array}$ (17)

$\{\begin{array}{l}\frac{4u_{x_{i+\frac{1}{2},j}}-u_{x_{i-\frac{1}{2},j}}-u_{x_{i+\frac{3}{2},j}}-u_{x_{i+\frac{1}{2},j-1}}-u_{x_{i+\frac{1}{2},j+1}}}{\Delta x\Delta y}+\frac{{\phi }_{i+\frac{1}{2},j+\frac{1}{2}}-{\phi }_{i-\frac{1}{2},j+\frac{1}{2}}}{\Delta x}=-\frac{1}{{\lambda }_u}{\left(K_{x_{i+\frac{1}{2},j}}^n,a^n\right)}_{L^2\left(0,T\right)},\\ \frac{4u_{y_{i,j+\frac{1}{2}}}-u_{y_{i-1,j+\frac{1}{2}}}-u_{y_{i+1,j+\frac{1}{2}}}-u_{y_{i,j-\frac{1}{2}}}-u_{y_{i,j+\frac{3}{2}}}}{\Delta x\Delta y}+\frac{{\phi }_{i+\frac{1}{2},j+\frac{1}{2}}-{\phi }_{i+\frac{1}{2},j-\frac{1}{2}}}{\Delta y}=-\frac{1}{{\lambda }_u}{\left(K_{y_{i,j+\frac{1}{2}}}^n,a^n\right)}_{L^2\left(0,T\right)},\\ \frac{u_{x_{i+\frac{1}{2},j}}-u_{x_{i+\frac{3}{2},j}}}{\Delta x}+\frac{u_{y_{i,j+\frac{1}{2}}}-u_{y_{i,j-\frac{1}{2}}}}{\Delta y}=0,\\ u_{x_{i+\frac{1}{2},0}}=u_{x_{i+\frac{1}{2},N_y}}=u_{y_{0,j+\frac{1}{2}}}=u_{y_{N_x,j+\frac{1}{2}}}=0,\end{array}$ (18)

$a^n={\mathbb{P}}_{\left[a_{\min },a_{\max }\right]}\left(\frac{1}{{\lambda }_a}{\displaystyle {\int }_{{\Omega }_c}{K}_{x_{i+\frac{1}{2},j}}^n{u}_{x_{i+\frac{1}{2},j}}+K_{y_{i,j+\frac{1}{2}}}^n{u}_{y_{i,j+\frac{1}{2}}}\text{d}x\text{d}y}\right).$ (19)

注：在以上两种格式中， $u_{x_{-\frac{1}{2},j}}$ 、 $u_{x_{N_x+\frac{1}{2},j}}$ 、 $u_{y_{i,-\frac{1}{2}}}$ 、 $u_{y_{i,N_y+\frac{1}{2}}}$ 为当x、y变量越过边界时取的“虚拟”函数值。为处理越界问题，采用线性插值，用越界点前(或后)两个点的线性表示来近似代替该点的值。

4. 数值实验

考虑区域 $\Omega ={\Omega }_c=\left[0,1\right]\times \left[0,1\right]$，时间域取 $\left[0,1\right]$。空间步长和时间步长分别取 $\Delta x=\Delta y=2^{-6}$， $\Delta t=\Delta x^2$。令 $\epsilon =1,\mu =10,\sigma =0,{\lambda }_u={\lambda }_a=0.1,a_{\min }=-9,a_{\max }=20$。选择电磁场

$E_{x_d}=-2\pi \sin \left(2\pi y\right)\cos \left(2\pi x\right),$

$E_{y_d}=2\pi \cos \left(2\pi y\right)\sin \left(2\pi x\right),$

$H_d=\cos \left(2\pi y\right)\cos \left(2\pi x\right),$

$E_x^0=-\frac{1}{10}\sin \left(\pi y\right)\cos \left(\pi x\right),$

$E_y^0=\frac{1}{10}\cos \left(\pi y\right)\sin \left(\pi x\right),$

$H^0=0,$

设置初始常函数 $a_0=4,u_{x_0}=E_{x_d},u_{y_0}=E_{y_d}$。在下述迭代算法中取 $eps={10}^{-10}$。

迭代算法为：

4.1. FDTD格式的数值模拟

以下给出了状态变量、对偶状态变量以及控制变量的数值模拟，具体的计算结果如图1所示，其中图1(a)描述的是最优电流密度函数 $u_x$ 的数值解，图1(b)描述的是最优电流密度函数 $u_y$ 的数值解，图1(c)描述的是 $t=1$ 时最优电场强度 $E_x$ 的数值解，图1(d)描述的是 $t=1$ 时最优电场强度 $E_y$ 的数值解，图1(e)描述的是 $t=1$ 时最优磁场强度H的数值解。

(a) 最优电流密度函数 $u_x$ 的数值解(b) 最优电流密度函数 $u_y$ 的数值解

(c) $t=1$ 时最优电场强度 $E_x$ 的数值解 (d) $t=1$ 时最优电场强度 $E_y$ 的数值(e) $t=1$ 时最优磁场强度H的数值解

Figure 1. Calculation results of the optimality system based on FDTD scheme

图1. 最优性系统基于FDTD格式的计算结果

(a) 最优电流密度函数 $u_x$ 的数值解(b) 最优电流密度函数 $u_y$ 的数值解

(c) $t=1$ 时最优电场强度 $E_x$ 的数值解 (d) $t=1$ 时最优电场强度 $E_y$ 的数值(e) $t=1$ 时最优磁场强度H的数值解

Figure 2. Results of the optimality system based on the fourth difference scheme

图2. 最优性系统基于四阶差分格式的计算结果

4.2. 四阶差分格式的数值模拟

以下给出了状态变量、对偶状态变量以及控制变量的数值模拟，具体的计算结果如图2所示，其中图2(a)描述的是最优电流密度函数 $u_x$ 的数值解，图2(b)描述的是最优电流密度函数 $u_y$ 的数值解，图2(c)描述的是 $t=1$ 时最优电场强度 $E_x$ 的数值解，图2(d)描述的是 $t=1$ 时最优电场强度 $E_y$ 的数值解，图2(e)描述的是 $t=1$ 时最优磁场强度H的数值解。

在选取相同大小的参数下，FDTD格式迭代三次便能满足 ${\Vert u_k-u_{k+1}\Vert }_{L^2\left({\Omega }_c\right)}+{\Vert a_k-a_{k+1}\Vert }_{L^2\left(0,T\right)}\le {10}^{-10}$，而四阶差分格式迭代四次达到了同样的约束条件。由于问题本身没有解析解，无法从误差大小的角度来比较精度，但是从图1(a)、图1(b)和图2(a)、图2(b)来看，四阶差分格式的图像更加平滑，因此四阶差分格式更加精确。

5. 总结

本文主要研究受Maxwell方程约束的最优控制问题，求解难度主要来源于Maxwell方程的复杂性。该最优控制问题带有两个控制 $u$ 和a，且是控制受控的。本文的工作是基于Lagrange乘子法和Helmholtz分解理论，将Maxwell方程最优控制问题转化为由状态方程、对偶状态方程、最优性条件三者联立形成的最优性系统，再分别利用FDTD格式、四阶差分格式建立最优性系统的离散格式，在此基础上分别求解最优性系统中的状态量、对偶状态量，再带入最优性条件中求解控制 $u$ 和a。为了避免求解大型耦合代数方程，采用迭代法进行计算，数值结果验证理论分析结论的正确性。

基金项目

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