一种基于对角化的抛物型最优控制问题的预处理子

A Preconditioner for Parabolic Optimal Control Problems Based on Diagonalization

摘要: 本文研究了求解抛物型偏微分方程约束的最优控制问题，利用对角化技巧，提出了一个新的基于对角化的预处理子，用于快速求解大型稀疏方程组。数值实验说明了预处理子的良好加速效果和稳定性。

Abstract: In this paper, we study the solution of the optimal control problem constrained by parabolic partial differential equations. By using the diagonalization technique, we propose a new preconditioner based on diagonalization to quickly calculate the large sparse equation system. Numerical experiments demonstrate the good acceleration effect and stability of the preconditioner.

1. 引言

最优控制是指在给定的约束条件下，寻求一个控制，使给定的系统性能指标达到极大值(或极小值)。这类问题广泛存在于社会问题或技术领域中。由于在绝大多数情况下，不能通过解析的方式求解最优控制问题，因此构造有效的数值计算方法成为必要。通常有两种不同的方法来处理偏微分方程约束的最优控制问题。一种是先优化后离散法，另一种是先离散后优化法，参见 [1] [2] [3] [4] 。本文选择先优化后离散的方法，来求解偏微分方程最优控制问题。

对于包含时间项的偏微分方程最优控制问题，状态方程随时间正向发展，伴随方程随时间反向发展。这意味着需要一次性求解正倒向两个方程耦合的大型代数方程组，因此导致计算成本很高。求解这类问题经常用到的方法是一种所谓的ParaDiag算法。该算法的核心步骤通常是将时间离散矩阵近似对角化，得到预处理子，然后进行分解，使得系统可以被一次性求解，以达到提高计算效率的目的，近年来该领域有许多相关的文献。这些文献中，ParaDiag算法大致被分为两类：① ParaDiag-I：使用直接法进行求解。这类方法最早由Maday [5] 等人提出，为了保证得到时间离散化矩阵可以很好地对角化，该方法会采取变化的时间步长或是使用混合的时间离散格式 [6] 。② ParaDiag-II：使用迭代法进行求解。该类方法允许通过强制修改系数矩阵来使得它可以对角化，近年来被许多文献所研究 [7] [8] [9] [10] 。

本文结构如下：第2节引入研究的抛物型方程约束的最优控制问题，对该问题进行时空间方向的离散。然后提出一种预处理子，用于加速计算大型稀疏线性方程组；第3节是数值实验，说明提出的预处理子的良好性能；第4节为全文的总结。

2. 离散格式与预处理

2.1. 全离散格式

考虑如下跟踪型分布式控制问题：

$\min J\left(y,u\right):=\frac{1}{2}{\Vert y-g\Vert }_{L^2\left(\Omega \times \left(0,T\right)\right)}^2+\frac{\gamma }{2}{\Vert u\Vert }_{L^2\left(\Omega \times \left(0,T\right)\right)}^2,$ (1)

使得

$\{\begin{array}{l}{y}_t-Ly=u+f,\left(x,t\right)\in \Omega \times \left[0,T\right]\\ y\left(x,t\right)=0,\left(x,t\right)\in \partial \Omega \times \left[0,T\right]\\ y\left(x,0\right)=y_0,x\in \Omega \end{array}$ (2)

其中L为空间算子，一般由拉普拉斯算子给出。y是状态变量，g为目标函数，u为控制变量， $\gamma >0$ 为正则化参数。

利用一阶最优性条件，消去控制变量u，得到如下简化的KKT系统 [11]

$\left[\begin{array}{cc}{L}_1& -{\gamma }^{-1}I\\ I& L_2\end{array}\right]\left[\begin{array}{c}y\\ p\end{array}\right]=\left[\begin{array}{c}f\\ g\end{array}\right].$ (3)

其中

$\{\begin{array}{l}{L}_{1}y:=\left({\partial }_t-L\right)y,y\left(x,0\right)=y_0,\\ L_{2}p:=\left(-{\partial }_t-L\right)p,p\left(x,T\right)=0.\end{array}$ (4)

为了求解系统(3)，对于空间离散化，采用有限差分格式(或有限元法，有限体积法等)。本文更加关注对时间方向离散化后的矩阵进行预处理，因此不论使用上述的何种方法对空间进行离散化，都可以利用空间离散矩阵 $L_h$ 代替具体的格式。即

$\{\begin{array}{l}{y}_t=f+L_{h}y-{\gamma }^{-1}p,\\ -p_t=g+L_{h}p-y.\end{array}$

时间上采用二阶精度的Crank-Nicolson离散格式，于是得到离散KKT系统

$\left[\begin{array}{cc}{B}_1\otimes I_x-\frac{\Delta t}{2}\left(B_2\otimes L_h\right)& -\frac{\Delta t}{2\sqrt{\gamma }}\left(B_2^{{\rm T}}\otimes I_x\right)\\ \frac{\Delta t}{2\sqrt{\gamma }}\left(B_2\otimes I_x\right)& B_1^{{\rm T}}\otimes I_x-\frac{\Delta t}{2}\left(B_2^{{\rm T}}\otimes L_h\right)\end{array}\right]\left[\begin{array}{c}\hat{Y}\\ P\end{array}\right]=\left[\begin{array}{c}\hat{F}\\ G\end{array}\right],$ (5)

其中

$Y=\left[\begin{array}{c}{y}_1\\ y_2\\ ⋮\\ y_{N_t}\end{array}\right],P=\left[\begin{array}{c}{p}_0\\ p_1\\ ⋮\\ p_{N_t-1}\end{array}\right],B_1=\left[\begin{array}{ccccc}1& & & & \\ -1& 1& & & \\ & -1& 1& & \\ & & \ddots & \ddots & \\ & & & -1& 1\end{array}\right],B_2=\left[\begin{array}{ccccc}1& & & & \\ 1& 1& & & \\ & 1& 1& & \\ & & \ddots & \ddots & \\ & & & 1& 1\end{array}\right],$

$F=\frac{\Delta t}{2}\left[\begin{array}{c}2y_0/\Delta t+L_h{y}_0+\left(f_0+f_1\right)\\ f_1+f_2\\ ⋮\\ f_{N_t-1}+f_{N_t}\end{array}\right],G=\frac{\Delta t}{2}\left[\begin{array}{c}-2y_0/\Delta t+\left(g_0+g_1\right)\\ g_1+g_2\\ ⋮\\ g_{N_t-1}+g_{N_t}\end{array}\right].$

其中 $I_t$ 是 $N_t\times N_t$ 的单位矩阵， $I_x$ 为 $N_x\times N_x$ 的单位矩阵， $N_x$ 表示空间离散剖分的网格数。“ $\otimes$ ”表示克罗内克积，是张量积的特殊形式。

2.2. 预处理

本节使用ParaDiag-II方法处理大型代数系统(5)，通过对矩阵 $B_1,B_2$ 进行修改，使得它们成为 $\alpha$ -循环矩阵。众所周知，循环矩阵都可以进行相似对角化，具体步骤参见 [12] ，因此在这一类方法中，基于对角化的预处理子被引入，以达到加速计算的目的。 $\alpha$ -循环矩阵结构如下：

$C_1=\left[\begin{array}{ccccc}1& & & & -\alpha \\ -1& 1& & & \\ & -1& 1& & \\ & & \ddots & \ddots & \\ & & & -1& 1\end{array}\right],C_2=\left[\begin{array}{ccccc}1& & & & \alpha \\ 1& 1& & & \\ & 1& 1& & \\ & & \ddots & \ddots & \\ & & & 1& 1\end{array}\right],$ (6)

$\alpha \in \left[-1,1\right]$ 是一个常数。有关 $\alpha$ -循环矩阵的介绍可以参见 [8] ， $\alpha$ -循环矩阵可以进行对角化分解

$\begin{array}{l}{C}_1={\Lambda }^{-1}{\mathbb{F}}^{\ast }{D}_1\mathbb{F}\Lambda ,D_1=diag\left({\lambda }_{1,1},{\lambda }_{1,2},\cdots ,{\lambda }_{1,N_t}\right),\\ C_2={\Lambda }^{-1}{\mathbb{F}}^{\ast }{D}_2\mathbb{F}\Lambda ,D_2=diag\left({\lambda }_{2,1},{\lambda }_{2,2},\cdots ,{\lambda }_{2,N_t}\right),\end{array}$

其中 $\Lambda =diag\left(1,{\alpha }^{\frac{1}{N_t}},{\alpha }^{\frac{2}{N_t}},\cdots ,{\alpha }^{\frac{N_t-1}{N_t}}\right)$ ， $\mathbb{F}$ 为离散傅立叶变换矩阵，两个对角矩阵 $D_1,D_2$ 的特征值分别为

${\lambda }_{1,n}=1-{\alpha }^{\frac{1}{N_t}}{\omega }^{n-1},{\lambda }_{2,n}=1+{\alpha }^{\frac{1}{N_t}}{\omega }^{n-1},n=1,2,\cdots ,N_t.$ (7)

其中 $\omega ={\text{e}}^{-2\pi i/N_t}$ ，容易知道 $\mathbb{F}$ 是一个酉矩阵，即 ${\mathbb{F}}^{\ast }\mathbb{F}=I$ ，这意味着 $\mathbb{F}$ 的逆矩阵与其共轭转置矩阵是相等的，观察(7)可知 ${\bar{\lambda }}_{1,n}={\lambda }_{2,n}$ 。利用 [7] 的结论，可以得到 $\tilde{P}$ 的分解

$\tilde{P}=\left(V\otimes I_x\right)\left(D_3\otimes I_x-\frac{\Delta t}{2}\left[\begin{array}{cc}{I}_t& \\ & I_t\end{array}\right]\otimes L_h\right)\left(V^{-1}\otimes I_x\right).$ (8)

其中

$V=\left[\begin{array}{cc}{\Lambda }_{}^{-1}{\mathbb{F}}^{\ast }& \\ & {\Lambda }_{}^{-1}{\mathbb{F}}^{\ast }\end{array}\right]\left[\begin{array}{cc}{I}_t& S_2\\ S_1& I_t\end{array}\right],S=\left[\begin{array}{cc}{I}_t& S_2\\ S_1& I_t\end{array}\right],$

$D_3=\left[\begin{array}{cc}{D}_2^{-1}{D}_1-\frac{\Delta t}{2\sqrt{\gamma }}{I}_t{S}_1& \\ & D_2^{}{D}_1^{-1}+\frac{\Delta t}{2\sqrt{\gamma }}{I}_t{S}_2\end{array}\right].$

$\begin{array}{l}{S}_1=-\frac{\sqrt{\gamma }}{\Delta t}{I}_t\left(D_2^{}{D}_1^{-1}-D_2^{-1}{D}_1+\sqrt{{\left(D_2^{}{D}_1^{-1}-D_2^{-1}{D}_1\right)}^2-\frac{\Delta t^{\text{2}}}{\gamma }{I}_t}\right),\\ S_2=-\frac{\sqrt{\gamma }}{\Delta t}{I}_t\left(D_2^{}{D}_1^{-1}-D_2^{-1}{D}_1+\sqrt{{\left(D_2^{}{D}_1^{-1}-D_2^{-1}{D}_1\right)}^2+\frac{\Delta t^{\text{2}}}{\gamma }{I}_t}\right)\text{.}\end{array}$

由 [13] 的分析可以知道，要成功地执行类似 [9] 的对角化分解，矩阵 $C_{1,2},C_{1,2}^{\ast }$ 必须同时都可以对角化，即

${\Lambda }_{}^{-1}={\Lambda }_{}^{\ast }.$

该等式和 $\left|\alpha \right|=1$ 是等价的，这个限制条件使得 $\tilde{P}$ 的形式有很大的局限性。和 [13] 中的选择类似，在本文后续的讨论中，取 $\alpha =-1$ ，此时 $C_{1,2}$ 又被称为负循环矩阵。

下面的步骤(1)~(3)给出了求解式(8)的方法。

为计算 $X_2={\tilde{P}}^{-1}r$ ，有

$\begin{array}{l}\begin{array}{cc}\left(\text{a}\right)& \end{array}{X}_1=\left(V^{-1}\otimes I_x\right)r,\\ \begin{array}{cc}\left(\text{b}\right)& \end{array}\left(d_n{I}_x-\frac{\Delta t}{2}{L}_h\right)w_n=X_{1,n},n=1,2,\cdots ,2N_t,\\ \begin{array}{cc}\left(\text{c}\right)& \end{array}\left(V\otimes I_x\right)w=X_2.\end{array}$ (9)

其中r表示任意的输入向量，

$X_1={\left[X_{1,1}^{{\rm T}},X_{1,2}^{{\rm T}},\cdots ,X_{1,2N_t}^{{\rm T}}\right]}^{{\rm T}},w={\left[w_1^{{\rm T}},w_2^{{\rm T}},\cdots ,w_{2N_t}^{{\rm T}}\right]}^{{\rm T}}.$

$d_n$ 表示对角矩阵 $D_3$ 的第n行元素。式(9)一般被称为三步对角化技巧 [6] 。该方法的好处在于：得益于矩阵V的特殊结构，步骤(a)和步骤(c)可以通过快速傅立叶变换进行高效计算；步骤(b)由2N个独立的线性系统构成，可以用高度并行的方式求解 [14] 。

根据 [15] 的分析，对角化步骤(a)和(c)会产生较大的舍入误差，定理2.1给出了适用于本文研究问题的舍入误差分析。

定理2.1令 $\lambda \in \sigma \left(L_h\right)$ 为空间离散矩阵 $L_h$ 的任意特征值，假设(8)的步骤(b)是通过直接法求解的(例如，LU分解法)。设机器精度为 $\epsilon$ ，那么舍入误差为

$\frac{{\Vert {\hat{\eta }}_{\lambda }-{\tilde{\eta }}_{\lambda }\Vert }_2}{{\Vert {\hat{\eta }}_{\lambda }\Vert }_2}\le \epsilon \left(4N_t+1\right){\text{Cond}}_2\left(\sum \left(\lambda \right)\right){\left({\text{Cond}}_2\left(S\right)\right)}^2+O\left({\epsilon }^2\right).$ (10)

其中 ${\hat{\eta }}_{\lambda }$ 为如下方程(11)中对任一输入向量 ${\hat{r}}_{\lambda }$ 的精确解， ${\tilde{\eta }}_{\lambda }$ 为方程(11)的数值近似解。

${\tilde{P}}_{\lambda }{\hat{\eta }}_{\lambda }:=\left(\left[\begin{array}{cc}{C}_2^{-1}{C}_1& -\frac{\Delta t}{2\sqrt{\gamma }}\\ \frac{\Delta t}{2\sqrt{\gamma }}& {\left(C_2^{-1}{C}_1\right)}^{{\rm T}}\end{array}\right]-\frac{\Delta t}{2}\left[\begin{array}{cc}{I}_t& \\ & I_t\end{array}\right]\right){\hat{\eta }}_{\lambda }={\hat{r}}_{\lambda },$ (11)

证明：对于方程(11)，根据 [7] 的向后误差分析，通过对角化技巧(2.4)得到的解满足如下的扰动系统：

$\left(V+\delta V_1\right)g_{\lambda }={\hat{r}}_{\lambda },\left(\Sigma \left(\lambda \right)+\delta \Sigma \right)w_{\lambda }=g_{\lambda },\left(V^{-1}+\delta V_2\right){\tilde{\eta }}_{\lambda }=w_{\lambda },$

其中 $\delta V_1$ ， $\delta \Sigma$ 和 $\delta V_2$ 表示矩阵V， $D_3$ 和 $V^{-1}$ 的微小扰动。由 [7] ，可得

$\begin{array}{l}{\Vert \delta V_1\Vert }_2\le 2\epsilon N{\Vert V\Vert }_2+O\left({\epsilon }^2\right),\\ {\Vert \delta V_2\Vert }_2\le 2\epsilon N{\Vert V^{-1}\Vert }_2+O\left({\epsilon }^2\right),\\ {\Vert \delta \sum \left(\lambda \right)\Vert }_2\le \epsilon {\Vert \sum \left(\lambda \right)\Vert }_2+O\left({\epsilon }^2\right).\end{array}$ (12)

最后一个不等式成立的原因是 $\sum \left(\lambda \right)$ 是一个对角矩阵。注意到通过对角化解方程(11)等价于精确求解 $\left({\tilde{P}}_{\lambda }+\delta {\tilde{P}}_{\lambda }\right){\tilde{\eta }}_{\lambda }={\hat{r}}_{\lambda }$ ，这里 $\delta {\tilde{P}}_{\lambda }$ 表示一个合适的扰动。实际上，有

$\left(V+\delta V_1\right)\left(\sum \left(\lambda \right)+\delta \sum \left(\lambda \right)\right)\left(V^{-1}+\delta V_2\right){\tilde{\eta }}_{\lambda }={\hat{r}}_{\lambda },$ (13)

结合(12)式和(13)可以对 $\delta {\tilde{P}}_{\lambda }$ 做出如下的估计：

${\Vert \delta {\tilde{P}}_{\lambda }\Vert }_2\le \epsilon \left(4N_t+1\right){\Vert V\Vert }_2{\Vert V^{-1}\Vert }_2{\Vert \sum \left(\lambda \right)\Vert }_2+O\left({\epsilon }^2\right).$ (14)

根据 [16] ，下列不等式成立

$\frac{{\Vert {\hat{\eta }}_{\lambda }-{\tilde{\eta }}_{\lambda }\Vert }_2}{{\Vert {\hat{\eta }}_{\lambda }\Vert }_2}\le {\text{Cond}}_2\left({\tilde{P}}_{\lambda }\right)\frac{{\Vert \delta {\tilde{P}}_{\lambda }\Vert }_2}{{\Vert {\tilde{P}}_{\lambda }\Vert }_2},$

结合(14)可以得到：

$\begin{array}{c}\frac{{\Vert {\hat{\eta }}_{\lambda }-{\tilde{\eta }}_{\lambda }\Vert }_2}{{\Vert {\hat{\eta }}_{\lambda }\Vert }_2}\le \epsilon \left(4N_t+1\right){\Vert V\Vert }_2{\Vert V^{-1}\Vert }_2{\Vert {\tilde{P}}^{-1}{}_{\lambda }\Vert }_2{\Vert \sum \left(\lambda \right)\Vert }_2+O\left({\epsilon }^2\right)\\ =\epsilon \left(4N_t+1\right){\text{Cond}}_2\left(V\right){\Vert {\tilde{P}}^{-1}{}_{\lambda }\Vert }_2{\Vert \sum \left(\lambda \right)\Vert }_2+O\left({\epsilon }^2\right).\end{array}$

由于 $V=\left[\begin{array}{cc}{\Lambda }_{\alpha }^{-1}{\mathbb{F}}^{\ast }& \\ & {\Lambda }_{\alpha }^{-1}{\mathbb{F}}^{\ast }\end{array}\right]S$ ，且 ${\tilde{P}}_{\lambda }=V\sum \left(\lambda \right)V^{-1}$ ，结合 $\mathbb{F}$ 是酉矩阵，可以得到

${\Vert V\Vert }_2\le {\Vert {\Lambda }_{\alpha }^{-1}\Vert }_2{\Vert {\mathbb{F}}^{*}\Vert }_2{\Vert S\Vert }_2\le \max \left\{1,{\left|\alpha \right|}^{\frac{N_t-1}{N_t}}\right\}\frac{1}{\sqrt{N_t}}{\Vert S\Vert }_2,$

同样地，有

${\Vert V^{-1}\Vert }_2\le {\Vert {\Lambda }_{\alpha }\Vert }_2{\Vert \mathbb{F}\Vert }_2{\Vert S^{-1}\Vert }_2\le \max \left\{1,{\left|\alpha \right|}^{-\frac{N_t-1}{N_t}}\right\}\frac{1}{\sqrt{N_t}}{\Vert S^{-1}\Vert }_2,$

于是

${\text{Cond}}_2\left(V\right)\le {\text{Cond}}_2\left(S\right),$

以及

${\Vert {\tilde{P}}^{-1}{}_{\lambda }\Vert }_2\le {\Vert S{\sum }^{-1}\left(\lambda \right)S^{-1}\Vert }_2\le {\text{Cond}}_2\left(S\right){\Vert {\sum }^{-1}\left(\lambda \right)\Vert }_2.$

综合上述讨论，便得到了不等式(10)。 ∎

上述定理表明，由(9)带来的舍入误差主要与三个因素有关，即时间离散的节点数 $N_t$ ，矩阵S的条件数以及不可避免的，对角矩阵 $\sum \left(\lambda \right)$ 的条件数。而(10)中 $\sum \left(\lambda \right)$ 的条件数这一因素无法进一步改进。于是本文考虑对 $\tilde{P}$ 进行修改，以使得矩阵S的条件数尽可能接近于或是等于一个网格无关的常数。

容易知道，消除 ${\text{Cond}}_2\left(S\right)$ 等价于令下列矩阵为零，

$D_2^{-1}{D}_1-D_1^{-1}{D}_2.$

用特征值来表征矩阵 $D_2^{-1}{D}_1$ 和 $D_1^{-1}{D}_2$ ，即

$\{\begin{array}{l}\lambda \left(D_2^{-1}{D}_1\right)=\frac{{\lambda }_{1,n}}{{\lambda }_{2,n}}=\frac{1-{\alpha }^{\frac{1}{N_t}}{\omega }^{n-1}}{1+{\alpha }^{\frac{1}{N_t}}{\omega }^{n-1}},\\ \lambda \left(D_1^{-1}{D}_2\right)=\frac{{\lambda }_{2,n}}{{\lambda }_{1,n}}=\frac{1+{\alpha }^{\frac{1}{N_t}}{\omega }^{n-1}}{1-{\alpha }^{\frac{1}{N_t}}{\omega }^{n-1}},\end{array}\begin{array}{cc}& n=1,2,\cdots ,N_t.\end{array}$ (15)

取

${\rho }_n=\frac{-2{\alpha }^{\frac{1}{N_t}}}{{\omega }^{-\left(n-1\right)}-{\alpha }^{\frac{2}{N_t}}{\omega }^{n-1}},\rho =diag\left({\rho }_1,{\rho }_2,\cdots ,{\rho }_{N_t}\right),$

于是有

$\left(D_2^{-1}{D}_1-\rho I_t\right)-\left(D_1^{-1}{D}_2+\rho I_t\right)=O$ . (16)

因此，可以用矩阵 $\left(D_2^{-1}{D}_1-\rho I_t\right)$ 代替 $D_2^{-1}{D}_1$ ，类似地，用矩阵 $\left(D_1^{-1}{D}_2+\rho I_t\right)$ 代替 $D_1^{-1}{D}_2$ 。如此一来，就可以得到

$\{\begin{array}{l}{\tilde{S}}_1=\frac{1}{2}\cdot \frac{-2\sqrt{\gamma }}{\Delta t}\cdot \frac{2i\Delta t}{2\sqrt{\gamma }}{I}_t=-i{I}_t,\\ {\tilde{S}}_2=-\frac{1}{2}\cdot \frac{2\sqrt{\gamma }}{\Delta t}\cdot \frac{2i\Delta t}{2\sqrt{\gamma }}{I}_t=-i{I}_t,\end{array}$ (17)

于是

$\tilde{S}{\tilde{S}}^{*}=\left[\begin{array}{cc}{I}_t& {\tilde{S}}_2\\ {\tilde{S}}_1& I_t\end{array}\right]\left[\begin{array}{cc}{I}_t& {\tilde{S}}_1^{\ast }\\ {\tilde{S}}_2^{\ast }& I_t\end{array}\right]=2\left[\begin{array}{cc}{I}_t& \\ & I_t\end{array}\right].$

矩阵 $\tilde{S}$ 的条件数为

${\text{Cond}}_2\left(\tilde{S}\right):={\Vert S\Vert }_2{\Vert S^{-1}\Vert }_2=\sqrt{2}{\Vert \frac{1}{\sqrt{2}}\tilde{S}\Vert }_2\frac{1}{\sqrt{2}}{\Vert \frac{1}{\sqrt{2}}{\tilde{S}}^{*}\Vert }_2={\Vert \frac{1}{\sqrt{2}}\tilde{S}\Vert }_2{\Vert \frac{1}{\sqrt{2}}{\tilde{S}}^{*}\Vert }_2=1.$

这样一来，便消除了和舍入误差有关的 ${\text{Cond}}_2\left(S\right)$ 这一因素，导出了预处理子 $\hat{P}$ ：

令

${\tilde{D}}_3=\left[\begin{array}{cc}{D}_2^{-1}{D}_1-\rho I_t-\frac{\Delta t}{2\sqrt{\gamma }}{I}_t{\tilde{S}}_1& \\ & D_2^{}{D}_1{}^{-1}+\rho I_t+\frac{\Delta t}{2\sqrt{\gamma }}{I}_t{\tilde{S}}_2\end{array}\right],$

预处理子 $\hat{P}$ 可以表示为

$\hat{P}=\left(V\otimes I_x\right)\left({\tilde{D}}_3\otimes I_x-\frac{\Delta t}{2}\left[\begin{array}{cc}{I}_t& \\ & I_t\end{array}\right]\otimes L_h\right)\left(V^{-1}\otimes I_x\right).$ (18)

值得注意的是，在消去舍入误差中 ${\text{Cond}}_2\left(S\right)$ 这一项时， ${\text{Cond}}_2\left(\sum \left(\lambda \right)\right)$ 同样会受到一定的影响。令

$\tilde{\sum }\left(\lambda \right)={\tilde{D}}_3-\frac{\lambda \Delta t}{2}\left[\begin{array}{cc}{I}_t& \\ & I_t\end{array}\right],$

因为 $\lambda \in \sigma \left(L_h\right)$ 是实对称矩阵的特征值，故 $\lambda =\bar{\lambda }$ ，可以使用矩阵的奇异值来估计 $\tilde{\sum }\left(\lambda \right)$ 以及 $\sum \left(\lambda \right)$ 的条件数，即

${\text{Cond}}_2\left(\sum \left(\lambda \right)\right)=\frac{{\sigma }_{\max }\left(\sum \left(\lambda \right)\right)}{{\sigma }_{\min }\left(\sum \left(\lambda \right)\right)}=\sqrt{\left|\frac{{\lambda }_{\max }\left(\sum \left(\lambda \right){\sum }^{*}\left(\lambda \right)\right)}{{\lambda }_{\min }\left(\sum \left(\lambda \right){\sum }^{*}\left(\lambda \right)\right)}\right|,}$

同样有

${\text{Cond}}_2\left(\tilde{\sum }\left(\lambda \right)\right)=\frac{{\sigma }_{\max }\left(\tilde{\sum }\left(\lambda \right)\right)}{{\sigma }_{\min }\left(\tilde{\sum }\left(\lambda \right)\right)}=\sqrt{\left|\frac{{\lambda }_{\max }\left(\tilde{\sum }\left(\lambda \right){\tilde{\sum }}^{*}\left(\lambda \right)\right)}{{\lambda }_{\min }\left(\tilde{\sum }\left(\lambda \right){\tilde{\sum }}^{*}\left(\lambda \right)\right)}\right|.}$

其中 ${\sigma }_{\max }(·),{\sigma }_{\min }(·)$ 分别表示最大奇异值和最小奇异值。

观察 $\sum \left(\lambda \right)$ 可知，它也是一个对角矩阵，因此特征值可以被对角线元素表征。考虑矩阵 $\sum \left(\lambda \right),\tilde{\sum }\left(\lambda \right)$ 的前 $N_t$ 行，有

$\begin{array}{l}{\lambda }_n\left(\sum \left(\lambda \right){\sum }^{*}\left(\lambda \right)\right)=\left({\lambda }_n\left(D_2^{-1}{D}_1\right)-\frac{\Delta t}{2\sqrt{\gamma }}{\lambda }_n\left(S_1\right)-\frac{\lambda \Delta t}{2}\right)\left({\lambda }_n\left(D_1^{-1}{D}_2\right)-\frac{\Delta t}{2\sqrt{\gamma }}{\lambda }_n^{*}\left(S_1\right)-\frac{\lambda \Delta t}{2}\right),\\ {\lambda }_n\left(\tilde{\sum }\left(\lambda \right){\tilde{\sum }}^{*}\left(\lambda \right)\right)=\left({\lambda }_n\left(D_2^{-1}{D}_1\right)+\frac{i\Delta t}{2\sqrt{\gamma }}-\frac{\lambda \Delta t}{2}-{\rho }_n\right)\left({\lambda }_n\left(D_1^{-1}{D}_2\right)-\frac{i\Delta t}{2\sqrt{\gamma }}-\frac{\lambda \Delta t}{2}-{\rho }_n^{*}\right),\end{array}$

${\lambda }_n(·)$ 表示第n个特征值， $n=1,2,\cdots ,N_t$ 。

令

$\Delta \kappa ={\text{Cond}}_2^2\left(\tilde{\sum }\left(\lambda \right)\right)-{\text{Cond}}_2^2\left(\sum \left(\lambda \right)\right),$

若

$\frac{{\lambda }_{\min }\left(\sum \left(\lambda \right){\sum }^{*}\left(\lambda \right)\right)}{{\lambda }_{\min }\left(\tilde{\sum }\left(\lambda \right){\tilde{\sum }}^{*}\left(\lambda \right)\right)}<\sqrt{\frac{{\lambda }_{\max }\left(\tilde{\sum }\left(\lambda \right){\tilde{\sum }}^{*}\left(\lambda \right)\right)}{{\lambda }_{\max }\left(\sum \left(\lambda \right){\sum }^{*}\left(\lambda \right)\right)}},$

则有如下的不等式估计

$\begin{array}{c}\Delta \kappa =\left|\frac{{\lambda }_{\max }\left(\tilde{\sum }\left(\lambda \right){\tilde{\sum }}^{*}\left(\lambda \right)\right)}{{\lambda }_{\min }\left(\tilde{\sum }\left(\lambda \right){\tilde{\sum }}^{*}\left(\lambda \right)\right)}\right|-\left|\frac{{\lambda }_{\max }\left(\sum \left(\lambda \right){\sum }^{*}\left(\lambda \right)\right)}{{\lambda }_{\min }\left(\sum \left(\lambda \right){\sum }^{*}\left(\lambda \right)\right)}\right|\\ <\frac{\left|{\lambda }_{\max }\left(\tilde{\sum }\left(\lambda \right){\tilde{\sum }}^{*}\left(\lambda \right)\right)\right|-\left|{\lambda }_{\max }\left(\sum \left(\lambda \right){\sum }^{*}\left(\lambda \right)\right)\right|}{\left|{\lambda }_{\min }\left(\tilde{\sum }\left(\lambda \right){\tilde{\sum }}^{*}\left(\lambda \right)\right)\right|-\left|{\lambda }_{\min }\left(\sum \left(\lambda \right){\sum }^{*}\left(\lambda \right)\right)\right|},\end{array}$

经过一系列冗杂但简单的运算可以得到

$\begin{array}{c}\Delta \kappa <\left|\frac{1+{\alpha }^{\frac{2}{N_t}}{\omega }^{2\left(n-1\right)}}{1-{\alpha }^{\frac{2}{N_t}}{\omega }^{2\left(n-1\right)}}-i\frac{\Delta t}{\sqrt{\gamma }}+\frac{\Delta t^2}{4\gamma }\right|\\ <\frac{\Delta t^2}{4\gamma }+c\left(\alpha ,\frac{\Delta t}{\sqrt{\gamma }}\right).\end{array}$ (19)

其中 $c\left(\alpha ,\frac{\Delta t}{\sqrt{\gamma }}\right)$ 是一个与 $\alpha$ 和 $\frac{\Delta t}{\sqrt{\gamma }}$ 有关的正数。于是可以知道，修改 $\tilde{P}$ 确实会增加 ${\text{Cond}}_2\left(\sum \left(\lambda \right)\right)$ ，但当时间方向的网格步长的平方 $\Delta t^2$ 和正则化参数 $\gamma$ 保持在同一数量级时， $\Delta \kappa$ 会有一个比较稳定的上界。

3. 数值算例

本节将使用带预处理子的GMRES方法来求解给定的抛物型方程最优控制问题。设置停机标准为 $\epsilon ={10}^{-6}$ ，采用的是右预处理GMRES方法，迭代的最大步数设置为IterMax = 400。空间离散算子 $L_h$ 由二阶中心差分格式导出，它是一个实对称矩阵。本节的实验目的是，针对具体的问题，在不同的正则化参数 $\gamma$ 以及时空间网格大小的选取下，获得程序输出的结果。输出的结果包括了GMRES求解器完成计算所需要的迭代次数，记为Iter、求解过程中计算机耗费的时间(单位为s)，记为CPU、状态变量的误差估计范数 $e_y$ 。此外，作为对比，将直接使用GMRES求解器而非使用预处理子来计算同样的算例。

例1给出了一维情况下抛物型PDE最优控制问题，数值结果如 表1 所示。 $N_t=64,N_x=64,\gamma ={10}^{-4}$ 时真解和数值解的对比图如图1和图2所示。

例1令空间域 $\Omega =\left(0,1\right)$ ，时间域 $T=\left(0,1\right)$ ，方程组的解析解等已知条件如下：

$\{\begin{array}{c}y=\sin \left(\pi x\right)\cos \left(\pi t\right),\\ p={\left(t-1\right)}^2\sin \left(\pi x\right),\\ f=-\pi \sin \left(\pi x\right)\sin \left(\pi t\right)-{\pi }^2\sin \left(\pi x\right)\cos \left(\pi t\right)-{\left(t-1\right)}^2\sin \left(\pi x\right)/\gamma ,\\ g=2\left(t-1\right)\sin \left(\pi x\right)+{\pi }^2\sin \left(\pi x\right){\left(t-1\right)}^2+\sin \left(\pi x\right)\cos \left(\pi t\right).\end{array}$ (20)

通过表1，可以看到，对于GMRES方法所需的迭代次数，无论正则化参数 $\gamma$ 取到何值，预处理过程都会大大减少迭代次数，使得GMRES求解过程的迭代次数保持在一个很少的范围。对于计算时间，如果正则化参数 $\gamma$ 的取值不过分的小，即 $\gamma \ge {10}^{-4}$ 时，不使用预处理方法的GMRES求解器所耗费的计算成本都远远大于经过预处理后的GMRES求解器，而在 $\gamma ={10}^{-6}$ 以及 $\gamma ={10}^{-8}$ 的情况下，使用预处理子与否对计算耗费的时间影响并没有显著差异，甚至网格较粗的情况下不使用预处理的GMRES方法在时间成本上稍显优势。对于此现象的解释是，正则化参数 $\gamma$ 与时空间网格步长的比 $\Delta t/\Delta x^2$ 之间的数量级有差异。如果增加网格大小，以使得 $\sqrt{\gamma }$ 和 $\Delta t/\Delta x^2$ 尽量接近，那么上述现象便不会再发生。为尝试验证这一解释，在 $\gamma ={10}^{-6}$ 和 $\gamma ={10}^{-8}$ 的情况下使用了更加细化的网格，大小为500 × 500，作为追加的实验。如表1中Add1和Add2的两行数据所示，当 $\sqrt{\gamma }$ 和 $\Delta t/\Delta x^2$ 比较接近时，施加了预处理过程的GMRES方法的数值表现明显更为优秀，即Add1中的结果；而 $\sqrt{\gamma }$ 和 $\Delta t/\Delta x^2$ 不那么接近时，两者在计算时间上的差异变得较小，即Add2中的结果。由此可以得到实验结果一定程度上与上述解释相符合。

Table 1. The number of iterations of the GMRES method in Example 1, the CPU time, and the error before and after using the preconditioner

表1. 使用预处理子前后，例1的GMRES方法迭代所需的次数、CPU耗费时间以及误差

Figure 1. Comparison of the real solution (left) and the numerical solution (right) of y

图1. y的真实解(左)和数值解(右)对比图

Figure 2. Comparison of the real solution (left) and the numerical solution (right) of p

图2. p的真实解(左)和数值解(右)对比图

4. 总结

本文研究了一种应用于初值型抛物型方程最优控制问题的预处理子，用于快速求解给定的模型问题。首先，给出了受抛物型方程约束的控制问题及初值条件，得到KKT系统；然后在时间方向上采用Crank-Nicolson格式，空间方向上使用离散拉普拉斯算子表示空间离散后的格式，得到了一个大型线性方程组。对离散后得到的大型稀疏代数系统，利用三步对角化技巧快速求解。为了消除对角化步骤中矩阵S的条件数这一因素，提出了一个新的预处理子。本文采用ParaDiag-II方案时，使用了广义极小残量法(GMRES)来作为迭代求解器，并分析了应用预处理子时，GMRES方法的收敛性。最后的数值实验，给出了一维的具体数值算例，数值结果表明提出的新的预处理子确实有较为优秀的加速效果和稳定性。

致谢

作者非常感谢各位老师们给出的宝贵意见和建议。