一类用于一般线性椭圆最优控制问题的共轭梯度算法
A Conjugate Gradient Algorithm for General Linear Elliptic Optimal Control Problems
摘要: 本文针对一类一般线性椭圆最优控制问题,应用了一种带有强Wolfe-Powell准则的共轭梯度算法对其进行求解,并从理论上分析了带有强Wolfe-Powell准则的共轭梯度算法的可行性以及全局收敛性。最后,我们做了相应的数值实验去验证算法的有效性。
Abstract: In this paper, we investigate a class of general linear elliptic optimal control problems. The conjugate gradient method with strong Wolfe-Powell criterion is used to solve this research problem. The feasibility and global convergence of conjugate gradient method with strong Wolfe-Powell criterion is theoretically proved. Finally, we do the corresponding numerical experiments to verify the feasibility of algorithm.
文章引用:林炯浩. 一类用于一般线性椭圆最优控制问题的共轭梯度算法[J]. 应用数学进展, 2025, 14(4): 355-365. https://doi.org/10.12677/aam.2025.144168

1. 引言

近年来,偏微分方程约束的最优控制问题在实际问题中得到了广泛的应用。例如:在生物医学领域,文献[1]研究了一个带有交叉扩散算子的耦合非线性抛物型最优控制问题,该问题描述了肿瘤细胞的密度、效应免疫细胞、循环淋巴细胞群和化疗药物浓度,以求达到控制药物的注射量并避免药物的副作用。在公共交通领域,文献[2]描述了多个高速列车运行的分布式控制,该模型克服了通信约束,实现了有效的速度控制。在生物研究领域,文献[3]用分布式最优控制问题描述生物模式中新的机械化学模型。在化学研究领域,文献[4]探究邻苯二甲酸酐的最大合成产率,建立了一个带有终端约束和控制变量约束的最优控制问题,并提出了一种基于罚函数方法和遗传算法的数值方法用于求解,最后得到了“最优反应温度”及“最佳反应物浓度”。

在本文中,我们主要关注以下的一般线性椭圆最优控制问题:

min u J( u )= 1 2 y y d L 2 ( Ω ) 2 + 1 2 u L 2 ( Ω ) 2 (1)

s.t.{ Ly+ c 0 y=u+f inΩ, y=0 onΩ,

其中 yY:= H 0 1 ( Ω ) uU:= L 2 ( Ω ) Ω R n  ( n=2,3 ) 是一个带有光滑边界 Ω 的开、凸区域,函数 y d ,f 均为已知项, α>0 为正则化参数,L为一致性椭圆微分算子并定义为:

Ly= i,j=1 n ( a ij y x i ) x j ,

其中 a ij , c 0 L ( Ω ) c 0 0 a ij = a ji 以及存在常数 θ>0 ,使得

i,j=1 n a ij ( x ) ε i ε j θ ε 2 ,a.exΩ,ε R n .

参考文献[5] [6]可知,对于任意的 1p<+ ,状态方程具有唯一解 y u H 0 1 ( Ω ) W 2,p ( Ω ) 。由于 J( y,u ) 是一致凸的,所以问题 ( 1 ) 具有唯一解。

针对线性椭圆最优控制问题的研究,目前已出现许多成果。如文献[7],陈和宋等人针对一个带有关于控制变量的“盒子约束”的一般线性椭圆最优控制问题,陈和宋等人先是使用了标准线性有限元方法对问题进行离散,再利用一种改进的多层次交替方向乘子法(mADMM)对离散化后的问题进行数值求解。陈和宋等人证明了提出的多层次交替方向乘子法(mADMM)具有良好的收敛性质以及算法复杂度 o( 1/k ) 。在文献[8]中,林同样利用有限元离散和一种修改的交替方向乘子法(mADMM)对一类具有混合控制–状态约束的线性椭圆最优控制问题进行研究,林从理论上证明了提出的修改的交替方向乘子法(mADMM)的全局收敛性和收敛速率。文献[9]则是针对一类带有 L 1 -“控制项”以及关于控制变量的逐点控制约束的线性椭圆最优控制问题,宋和陈等人讨论了多种类型的有限元先验误差估计,同时提出了一种新方案去离散 L 1 范数并获得了相应的先验误差估计,误差分析证明了提出的离散方案并不影响误差估计的精度阶。文献[10]则探讨了一种虚单元法并将其应用于一类带有逐点控制约束的线性椭圆最优控制问题,王和周分别推导了关于控制变量、状态变量以及协态变量的先验误差估计和后验误差估计。

针对定义在规则区域上的椭圆偏微分方程,五点中心有限差分方法的误差估计精度阶为 o( h 2 ) ,具有较高的收敛速度,以及强Wolfe-Powell准则在收敛速度以及稳定性方面都做了相应的改进,具有较快的收敛速度以及较强的数值稳定性。

在本文中,我们主要参考文献[11]以及文献[6]的共轭梯度算法,并将其应用于求解问题(1),并分析算法的可行性以及全局收敛性。

接下来我将简要介绍在本文中用到的两个算法,即五点中心有限差分方法以及带强Wolfe-Powell准则的共轭梯度算法。

2. 预备知识

首先,我们定义问题(1)的梯度如下:

J( y,u )=αup=g( u ),

其中 : L 2 ( Ω ) L 2 ( Ω ) 是梯度算子,p

{ Lp+ c 0 p+y= y d inΩ, p=0 onΩ, (2)

的解。y

{ Ly+ c 0 y=u+f inΩ, y=0 onΩ, (3)

的解。为了保证共轭梯度法的顺利运行,我们归纳需数值求解的方程如下:

{ L y i + c 0 y i = u i +f inΩ, y i =0 onΩ, L p i + c 0 p i + y i = y d inΩ, p i =0 onΩ, (4)

接下来,我们便利用五点中心有限差分方法对问题(4)进行离散,我们只给出有限差分离散的过程和结果,关于有限差分方法的误差估计研究部分,详细可查阅文献[12]

首先,我们将问题定义在二维情况,即 Ω=( a,b )×( c,d ) R 2 。其次,我们使用均匀笛卡尔网格对区间Ω进行均匀划分,假设“x轴方向”的划分区间为 N 1 ,“y轴方向”的划分区间为 N 2 ,进一步我们可以得到“x轴方向”的步长为 h 1 = ( ba )/ N 1 以及“y轴方向”的步长为 h 2 = ( dc )/ N 2 。因此,我们可以得到Ω的网格点集:

Ω h ={ ( x 1 k , x 2 j )=( a+k h 1 ,c+j h 2 )|k=0,, N 1 ;j=0,, N 2 }

以及 Ω h 的边界点集 Ω h

Ω h ={ ( a,c+j h 2 ),( b,c+j h 2 ),( a+k h 1 ,c ),( a+k h 1 ,d )|k=0,, N 1 ;j=0,, N 2 }.

为了表述的方便,我们定义 u k,j i y k,j i p k,j i f k,j y dk,j 分别为 u i ( x 1 k , x 2 j ) y i ( x 1 k , x 2 j ) p i ( x 1 k , x 2 j ) f( x 1 k , x 2 j ) y d ( x 1 k , x 2 j ) 的近似, u h y h p h f h y dh 分别为 [ u k,j ] n1×n1 [ y k,j ] n1×n1 [ p k,j ] n1×n1 [ f k,j ] n1×n1 [ y dk,j ] n1×n1 的近似。进而我们可以得到(4)的五点中心有限差分格式如下:

{ L h y h i + c 0 y h i = u h i + f h in Ω h , y h i =0 on Ω h (5)

{ L h p h i + c 0 p h i + y h i = y dh in Ω h , p h i =0 on Ω h (6)

其中 L h 为一致性椭圆微分算子L的离散化形式, L h y h = y h A h 1 + A h 2 y h 以及 A h m = tridiag( a mm ,2 a mm , a mm )/ h m 2 m=1,2 。同时,问题(1)的目标泛函的离散化格式也表述如下:

J h ( u h )= 1 2 y h i y dh L h 2 ( Ω h ) 2 + α 2 u h i L h 2 ( Ω h ) 2 ,

其中 L h 2 ( Ω h ) 为定义在 Ω h 上的离散 L 2 范数,并满足:对任意 N+1 阶矩阵 w h ,都有

w h L h 2 ( Ω h ) 2 = ( w h , w h ) L h 2 ( Ω h ) = h 1 h 2 4 l l T ( w h w h ),

其中 l= [ 1,2,2,,2,1 ] T N+1 维向量, ( w h w h )=[ w h ( i,j )× w h ( i,j ) ]

接下来,我们进一步阐述用于求解问题(5)~(6)的共轭梯度算法。

首先,我们定义在数值迭代求解过程中,对于所有正整数 i0 ,共轭梯度更新公式如下:

d h i = g h i + β i d h i1 , (7)

其中 g h i 是梯度 g( u h i ) 的缩写,调整系数 β i 满足

{ β 0 =0, β i = g h i g h i1 ,i=1,2, (8)

进而我们定义控制变量更新公式如下:

u h i = u h i1 + τ i d h i1 .

其中参数 τ i 满足强Wolfe-Powell准则:

{ J h ( u h i + τ i d h i ) J h ( u h i )+ρ τ i ( g h i , d h i ), | ( g h i+1 , d h i ) |σ( g h i , d h i ) (9)

以及 0<ρ<σ<0.5 。梯度更新公式定义如下:

g h i+1 = g h i + τ i g ˜ h i , (10)

其中

g ˜ h i =α d h i p 1h i .

p 1h i

{ L h p 1h i + c 0 p 1h i = y 1h i in Ω h , p 1h i =0 on Ω h (11)

的解。 y 1h i

{ L h y 1h i + c 0 y 1h i = d h i in Ω h , y 1h i =0 on Ω h

的解。我们进一步构造状态变量和协态变量的更新公式:

{ y h i = y h i1 + τ i1 y 1h i1 , p h i = p h i1 + τ i1 p 1h i1 .

综合以上信息,我们给出结合了强Wolfe-Powell准则的共轭梯度算法的流程:

输入:剖分区间数 N 1 , N 2 ,控制变量初始估计 u h 0 ,误差精度 ε

输出:控制变量数值解 u h i ,状态变量数值解 y h i

步骤 i=0 ,对

{ L h y h i + c 0 y h i = u h i + f h , y h i =0, L h p h i + c 0 p h i + y h i = y dh , p h i =0,

进行数值求解,进而得到数值解 y h i p h i

步骤二:计算梯度

g h i =α u h i p h i

以及共轭梯度

{ d h 0 = g h 0 i=0, d h i = g h i + β i d h i1 i1.

步骤三:如果 g h i L h 2 ( Ω h ) ε ,算法停止迭代,输出数值解。

否则,对

{ L h p 1h i + c 0 p 1h i = y 1h i in Ω h , p 1h i =0 on Ω h , L h y 1h i + c 0 y 1h i = d i in Ω h , y 1h i =0 on Ω h

进行数值求解,以得到 p 1h i y 1h i ,从而我们可以计算 g ˜ h i =α d h i p 1h i 。同时,根据强Wolfe-Powell准则确定最优迭代步长 τ i 。从而进一步更新控制变量和状态变量

{ u h i+1 = u h i + τ i d i , y h i+1 = y h i + τ i y 1h i .

i=i+1 ,返回步骤二。

3. 性质分析

接下来我们将讨论带强Wolfe-Powell准则的共轭梯度算法的可行性以及全局收敛性。为了方便,我

们把 L h 2 ( Ω h ) 简写为 。首先,我们先讨论带强Wolfe-Powell准则的共轭梯度算法的可行性,即在每一

步迭代中,均存在 τ i 满足强Wolfe-Powell准则。在证明强Wolfe-Powell准则在每一步迭代中均成立之前,我们先给出以下引理:

引理3.1假定共轭梯度 d h i 由(6)~(7)定义以及 τ i 满足强Wolfe-Powell准则(8),则对于所有 i0 ,有

1 σ1 g h i 2 ( d h i , g h i ) 2σ1 1σ g h i 2 .

证明: i=0 时, d h 0 = g h 0 ,此时 ( g h 0 , d h 0 )= g h 0 2 成立,显然

1 σ1 g h 0 2 ( d h 0 , g h 0 ) 2σ1 1σ g h 0 2

成立。

i1 时,我们使用归纳法进行证明。

假设当 1in 时,

1 σ1 g h i 2 ( d h i , g h i ) 2σ1 1σ g h i 2

成立。

i=n+1 时,根据强Wolfe-Powell准则(8)第二个不等式,我们可以得到

( g h n+1 , d h n+1 )=( g h n+1 , g h n+1 + β n+1 d h n ) = g h n+1 2 + g h n+1 2 g h n 2 ( g h n+1 , d h n ) g h n+1 2 + g h n+1 2 g h n 2 ( σ )( g h n , d h n ).

由于当 1in 时, ( g h i , d h i ) 1 σ1 g h i 2

因此,

g h n+1 2 + g h n+1 2 g h n 2 ( σ )( g h n , d h n ) g h n+1 2 + σ 1σ g h n+1 2 = 2σ1 1σ g h n+1 2 .

综合以上结果,我们可以得到

( g h n+1 , d h n+1 ) 2σ1 1σ g h n+1 2 .

同理可得

( g h n+1 , d h n+1 )= g h n+1 2 + g h n+1 2 g h n 2 ( g h n+1 , d h n ) g h n+1 2 + g h n+1 2 g h n 2 σ( g h n , d h n ) g h n+1 2 + g h n+1 2 g h n 2 σ 1 σ1 g h n 2 = 1 σ1 g h n+1 2 .

故引理得证。 □

定理3.1对任意 i0 ,都存在 τ i 满足强Wolfe-Powell准则(8)。

证明:先分析(8)中第二个不等式,将(9)代入(8)中第二个不等式,可以得到

( σ1 )( g h i , d h i ) τ i ( g ˜ h i , d h i )( 1+σ )( g h i , d h i ).

显然,当

( σ1 ) ( g h i , d h i ) ( g ˜ h i , d h i ) τ i ( 1+σ ) ( g h i , d h i ) ( g ˜ h i , d h i )

时, | ( g h i+1 , d h i ) |σ( g h i , d h i ) 成立。

接下来分析(8)中第一个不等式,易知

J h ( u h i+1 ) J h ( u h i )= 1 2 y h i+1 y dh 2 α 2 u h i+1 2 ( 1 2 y h i y dh 2 α 2 u h i 2 ) =( y h i y dh , τ i y 1h i )+ ( τ i ) 2 2 y 1h i 2 +α τ i ( u h i , d h i )+ α ( τ i ) 2 2 d h i 2 . (12)

根据 ( , ) 的定义,可以得到

( y h i y dh , y 1h i )=( L h p h i + c 0 p h i , y 1h i )=( L h y 1h i + c 0 y 1h i , p h i )=( d h i , p h i ). (13)

同理可以得到

y 1h i 2 =( y 1h i , y 1h i )=( L h p 1h i + c 0 p 1h i , y 1h i )=( L h y 1h i + c 0 y 1h i , p 1h i )=( d h i , p 1h i ). (14)

把(13)、(14)代入到(12)中,以及利用(9)中第一个不等式,得到

J h ( u h i+1 ) J h ( u h i )= ( τ i ) 2 2 ( d h i , g ˜ h i )+ τ i ( d h i , g h i )ρ τ i ( d h i , g h i ).

显然,当

0 τ i 2( ρ1 ) ( d h i , g h i ) ( d h i , g ˜ h i )

时,不等式

J h ( u h i + τ i d h i ) J h ( u h i )+ρ τ i ( g h i , d h i )

成立。

由于 σ1>2( ρ1 ) ,因此,存在

τ i [ ( σ1 ) ( d h i , g h i ) ( d h i , g ˜ h i ) ,min{ 2( ρ1 ) ( d h i , g h i ) ( d h i , g ˜ h i ) ,( 1+σ ) ( d h i , g h i ) ( d h i , g ˜ h i ) } ],

使得强Wolfe-Powell准则(8)成立。 □

以上我们证明了带强Wolfe-Powell准则的共轭梯度算法的可行性。除了算法的可行性之外,算法的收敛性也是讨论的重点。共轭梯度算法的全局收敛性是其得到广泛应用的一个原因之一,接下来,我们便讨论带强Wolfe-Powell准则的共轭梯度算法的全局收敛性。首先,我们先引入一个引理。

引理3.2[6]假定 g h i d h i 分别由(10)、(7)定义,则 g h i d h i 满足:

lim n i=0 n ( w i ) 2 g h i 2 =0,

其中

w i = ( g h i , d h i ) g h i 2 d h i 2 ,i=0,1,.

定理3.2带强Wolfe-Powell准则的共轭梯度算法满足:

lim n inf in g h i 2 =0.

证明:我们将利用反证法证明以上定理。

不妨设当 i0 时,存在常数 ε 2 > ε 1 ,使得

ε 2 > g h i 2 > ε 1 . (15)

根据 d h i β i 的定义,可以得到

d h i+1 2 =( d h i+1 , d h i+1 ) =( g h i+1 + β i+1 d h i , g h i+1 + β i+1 d h i ) =( g h i+1 , g h i+1 )2( β i+1 d h i , g h i )+( β i+1 d h i , β i+1 d h i ) = g h i+1 2 2 g h i+1 2 g h i 2 ( d h i , g h i+1 )+ g h i+1 4 g h i 4 d h i 2 . (16)

根据强Wolfe-Powell准则(8)中第二个不等式以及引理3.1,(16)满足

d h i+1 2 ( 1+σ 1σ ) g h i+1 2 + g h i+1 4 g h i 4 d h i 2 .

t i = d h i+1 2 g h i+1 4

以及

α= 1+σ ( 1σ ) ε 1 .

进一步根据(15)以及强Wolfe-Powell准则(8)中第二个不等式,可得

t i α+ t i1 iα+ t 0 ,

并满足

i=0 n ( w i ) 2 g h i 2 1 ε 2 3 ( 2σ1 ) 2 ( 1σ ) 2 i=1 n ( t i1 ) 2 1 ε 2 3 ( 2σ1 ) 2 ( 1σ ) 2 i=1 n ( 1 iα+ t 0 ) 2 .

由于

lim n 1 ε 2 3 ( 2σ1 ) 2 ( 1σ ) 2 i=1 n ( 1 iα+ t 0 ) 2 >0,

lim n i=0 n ( w i ) 2 g h i 2 lim n 1 ε 2 3 ( 2σ1 ) 2 ( 1σ ) 2 i=1 n ( 1 iα+ t 0 ) 2 .

lim n i=0 n ( w i ) 2 g h i 2 >0.

这与引理3.2的结果相矛盾,故

lim n inf in g h i 2 =0.

4. 数值实验

为了验证带强Wolfe-Powell准则的共轭梯度算法的有效性,我们利用以下问题对算法进行验证。

Ω= ( 0,1 ) 2 L=Δ 为拉普拉斯算子, c 0 =1 α=1 以及

f( x 1 , x 2 )=( 12 π 2 )sinπ x 1 sinπ x 2 sin2π x 1 sin2π x 2 ,

y d ( x 1 , x 2 )=( 18 π 2 )sin2π x 1 sin2π x 2 +sinπ x 1 sinπ x 2 .

以上问题存在精确解 u * ( x 1 , x 2 )=sin2π x 1 sin2π x 2

当剖分区间长度为 N 1 = N 2 =10 ,以及“停机精度” ε= 10 5 ,我们得到如下的数值结果。在图1,我们给出了精确解 u * ( x 1 , x 2 ) 的图像。在图2,我们给出了带强Wolfe-Powell准则的共轭梯度算法计算的结果。图3则是描述了数值解与精确解之间的误差,从图3可以看出数值解与精确解之间是比较接近的。当剖分区间长度为 N 1 = N 2 =10 ,以及“停机精度” ε= 10 4 ,我们将实验结果与文献[11]的实验结果相比较。结果显示:我们的误差精度阶为104,文献[11]的误差精度阶为101,显然,我们的算法精度更高。

Figure 1. Exact control

1. 精确控制

Figure 2. Numerical control

2. 数值控制

Figure 3. Error

3. 误差

5. 结语

本文利用有限差分法对一般线性椭圆最优控制问题进行离散化,并提出了带有强Wolfe-Powell准则的共轭梯度算法来求解有限维空间上的优化问题。此外,本文还验证了带有强Wolfe-Powell准则的共轭梯度算法在求解过程中的可行性,以及从理论上分析了带有强Wolfe-Powell准则的共轭梯度算法具有全局收敛性。接下来,我们会进一步研究带有强Wolfe-Powell准则的共轭梯度算法的收敛速度,以及如何提高带有强Wolfe-Powell准则的共轭梯度算法的计算速度。

参考文献

[1] Sowndarrajan, P.T., Manimaran, J., Debbouche, A. and Shangerganesh, L. (2019) Distributed Optimal Control of a Tumor Growth Treatment Model with Cross-Diffusion Effect. The European Physical Journal Plus, 134, Article No. 463.
https://doi.org/10.1140/epjp/i2019-12866-8
[2] Li, S., Yang, L. and Gao, Z. (2020) Distributed Optimal Control for Multiple High-Speed Train Movement: An Alternating Direction Method of Multipliers. Automatica, 112, Article 108646.
https://doi.org/10.1016/j.automatica.2019.108646
[3] Liu, C. and Zhang, X. (2019) Optimal Distributed Control for a New Mechanochemical Model in Biological Patterns. Journal of Mathematical Analysis and Applications, 478, 825-863.
https://doi.org/10.1016/j.jmaa.2019.05.057
[4] Antipina, E.V., Mustafina, S.I., Antipin, A.F. and Mustafina, S.A. (2020) A Numerical Algorithm for Solving Optimal Control Problems with Terminal Constraints for Dynamical Systems. Optoelectronics, Instrumentation and Data Processing, 56, 671-678.
https://doi.org/10.3103/s8756699020060035
[5] Casas, E. (2006) Using Piecewise Linear Functions in the Numerical Approximation of Semilinear Elliptic Control Problems. Advances in Computational Mathematics, 26, 137-153.
https://doi.org/10.1007/s10444-004-4142-0
[6] Zhang, Z. and Chen, X. (2021) A Conjugate Gradient Method for Distributed Optimal Control Problems with Nonhomogeneous Helmholtz Equation. Applied Mathematics and Computation, 402, Article 126019.
https://doi.org/10.1016/j.amc.2021.126019
[7] Chen, X., Song, X., Chen, Z. and Yu, B. (2020) A Multi-Level ADMM Algorithm for Elliptic PDE-Constrained Optimization Problems. Computational and Applied Mathematics, 39, Article No. 331.
https://doi.org/10.1007/s40314-020-01379-1
[8] 林继桐. 椭圆最优控制问题的修改交替方向乘子计算方法[J]. 应用数学进展, 2022, 11(11): 7936-7945.
https://doi.org/10.12677/aam.2022.1111840
[9] Song, X., Chen, B. and Yu, B. (2017) Error Estimates for Sparse Optimal Control Problems by Piecewise Linear Finite Element Approximation. arXiv: 1709.09539.
https://doi.org/10.48550/arXiv.1709.09539
[10] Wang, Q. and Zhou, Z. (2022) A Priori and a Posteriori Error Analysis for Virtual Element Discretization of Elliptic Optimal Control Problem. Numerical Algorithms, 90, 989-1015.
https://doi.org/10.1007/s11075-021-01219-1
[11] Li, B. and Liu, S. (2007) Conjugate Gradient-Boundary Element Solution for Distributed Elliptic Optimal Control Problems. Journal of Mathematical Analysis and Applications, 335, 1219-1237.
https://doi.org/10.1016/j.jmaa.2007.02.003
[12] Chen, X. (2005) Finite Difference Smoothing Solutions of Nonsmooth Constrained Optimal Control Problems. Numerical Functional Analysis and Optimization, 26, 49-68.
https://doi.org/10.1081/nfa-200051629