求解線性橢圓型優(yōu)化控制問題的多重網(wǎng)格方法

邱 岳 東

(暨南大學 信息科學技術(shù)學院，廣州 510632)

有限差分法是計算機數(shù)值模擬最早采用的方法，該方法將求解域劃分為差分網(wǎng)格，用有限個網(wǎng)格節(jié)點代替連續(xù)的求解域．有限差分法主要包括三個過程，即網(wǎng)格生成、有限差分離散方程的形成以及離散方程的求解．其中，有限差分離散方程的求解能夠影響到數(shù)值計算結(jié)果的精確性和數(shù)值計算時間的高效性[1-2]．

我們通?？梢酝ㄟ^加密網(wǎng)格來提高計算精度．但是，有限差分法迭代計算的收斂速度隨著計算網(wǎng)格的增加而越來越慢[3-4]，計算精度無法繼續(xù)提高，因此，我們需要尋找一種高效的求解方法．

多重網(wǎng)格法是求解大規(guī)?？茖W計算問題的有效方法，其基本原理在于一定的網(wǎng)格最容易消除與網(wǎng)格步長相對應的誤差分量，對于離散方程，它在細網(wǎng)格上利用通常的迭代法消去殘量的高頻部分，然后把低頻部分轉(zhuǎn)移到粗網(wǎng)格上進行校正，經(jīng)過若干次循環(huán)以后，獲得滿足精度的解[5]．該算法主要包括以下幾個部分：網(wǎng)格剖分、網(wǎng)格粗化方法、模型問題的離散格式、松弛方法、限制與延拓以及迭代技術(shù)等．

本文利用集體平滑多重網(wǎng)格方法和多重網(wǎng)格最優(yōu)化方法求解線性橢圓型優(yōu)化控制問題，研究分析了多重網(wǎng)格法的網(wǎng)格剖分稀疏程度和最大循環(huán)次數(shù)對計算速度及精度的影響．文獻[6]通過局部傅里葉分析證明了集體平滑多重網(wǎng)格方法的收斂性，而文獻[7]則利用了利普希茨連續(xù)性證明了多重網(wǎng)格最優(yōu)化方法的收斂性．此外，數(shù)值試驗說明，這兩種方法的收斂因子和迭代步數(shù)不隨網(wǎng)格步長變化而增加，并且，集體平滑多重網(wǎng)格方法比多重網(wǎng)格最優(yōu)化方法的收斂速度更快和更好的穩(wěn)定性．

1 模型問題

1.1 優(yōu)化控制問題的模型

優(yōu)化控制是指在給定的約束條件下，尋求一個控制系統(tǒng)，使給定的被控系統(tǒng)性能指標取得最大值或最小值的控制．為了求解線性橢圓型優(yōu)化控制問題，我們選取如下目標泛函：

(1)

其中:Ω為R2上的一個凸集，v>0為控制目標泛函的權(quán)重，z∈L2(Ω)是目標函數(shù)，而u是在Ω上的控制變量的分布函數(shù)．

此外，我們選取的目標泛函還需滿足一定的偏微分方程約束優(yōu)化條件，因此，在凸集Ω上的線性橢圓型優(yōu)化控制問題的模型可以表示為：

minu∈UH(y,u),

-Δy-u=f,y|?Ω=0

-Δp+y=z,p|?Ω=0

vu-p=0.

(2)

其中:f∈L2(Ω)和U=L2(Ω)，p為拉格朗日算子．

1.2 優(yōu)化控制問題的離散形式

本文考慮線性橢圓型控制問題偏微分方程約束條件基于有限差分法的離散化形式．首先，用Ωk來表示區(qū)域Ω離散化的參數(shù)形式，網(wǎng)格大小為hk(k=1,2,…L)，令hk-1=2hk，內(nèi)部網(wǎng)格點的數(shù)目為nk，任意函數(shù)在區(qū)域Ωk里進行離散化之后，都可用一個nk維的向量表示．

那么，偏微分方程約束條件的離散形式如下：

-Δkyk-uk=fk,

-Δkpk+yk=zk,

vuk-pk=0.

(3)

這里，-Δk表示負五點拉普拉斯算子，x∈Ωk且x=(ihk,jhk)，i,j是網(wǎng)格點的按字典序的下標指數(shù)．因此，我們由Taylor展式可得

(4)

其中:在變量yi,j和pi,j的更新過程中，假定Ai,j、Bi,j的值為常量，并且可得

Ai,j=-(yi-1,j+yi+1,j+yi,j-1+yi,j+1)-h2fi,j,

Bi,j=-(pi-1,j+pi+1,j+pi,j-1+pi,j+1)-h2zi,j

(5)

將式(4)經(jīng)過化簡變形，并且定義yi,j=yi,j(ui,j),pi,j=pi,j(ui,j)，可以得到關(guān)于ui,j的四次多項式方程，并且根據(jù)Cardano-Tartaglia公式求得四次多項式方程的最小實數(shù)值：

(6)

該四次多項式方程的最小實數(shù)根即為橢圓型優(yōu)化控制問題目標泛函的有限差分離散化問題的解，有了這一條件，我們可以得到一個完整且有效的平滑迭代．

2 集體平滑多重網(wǎng)格方法和多重網(wǎng)格最優(yōu)化方法

集體平滑多重網(wǎng)格方法是基于一種應用到具有集體平滑性質(zhì)的非線性多重網(wǎng)格全近似存儲方法，主要是通過解決相配的偏微分方程在處理約束條件中的優(yōu)化變量來解決優(yōu)化控制問題，而且這種方法需要對問題進行集體平滑處理與內(nèi)部網(wǎng)格算子轉(zhuǎn)移．文獻[8-9]都較詳細地介紹了集體平滑多重網(wǎng)格方法在求解狀態(tài)和控制約束條件下線性優(yōu)化控制問題的應用．

Nash[10]和Lewis[11]首次將多重網(wǎng)格最優(yōu)化方法應用到約束優(yōu)化控制問題的求解．在應用多重網(wǎng)格最優(yōu)化方法求解控制問題過程中，我們假設(shè)其外循環(huán)的控制函數(shù)為惟一獨立的變量，而內(nèi)循環(huán)則由一個單重網(wǎng)格和其他多重網(wǎng)格法構(gòu)成的，而且，多重網(wǎng)格最優(yōu)化方法具有更簡易、更廣泛的應用性[12]．

2.1 集體平滑多重網(wǎng)格方法的結(jié)構(gòu)

為了闡明集體平滑多重網(wǎng)格方法，我們首先考慮如下方程組

Ak(ωk)=fk,

(7)

前光滑：

(8)

粗網(wǎng)格問題：

1) 計算細網(wǎng)格虧損量：

4)虧損量校正(由細網(wǎng)格到粗網(wǎng)格)

(9)

對粗網(wǎng)格問題Ak-1(ωk-1)=fk-1應用集體平滑算法進行r次循環(huán)(r=r1+r2)求得wk-1．

粗網(wǎng)格校正：

后光滑：

在Ωk上再對方程Ak(ωk)=fk應用γ2次平滑算法迭代，得到

其中:rk為殘差因子，ωk-1表示一個相對于ωk的粗網(wǎng)格估計，同時，Ak-1(·)表示相對于Ak(·)較粗網(wǎng)格的散方程組，而且τk-1表示由細網(wǎng)格到粗網(wǎng)格的殘差校正．

2.2 多重網(wǎng)格最優(yōu)化方法的結(jié)構(gòu)

為了能夠更好地了解多重網(wǎng)格最優(yōu)化方法，我們考慮一個局部離散的凸規(guī)劃問題

(10)

而求解這個凸規(guī)劃問題就等同于求解

▽Hk(uk)=fk．

(11)

其中:Hk(·)表示目標泛函，那么，在優(yōu)化控制問題上，為了計算Hk(·)在uk的值，我們應計算yk(uk)的值，并且，還需要yk(uk)、pk(uk)來求解梯度函數(shù)▽Hk(·)的值．假設(shè)Sk為在向量空間Vk上的一個優(yōu)化迭代算法，使得

成立．那么，多重網(wǎng)格最優(yōu)化方法的結(jié)構(gòu)為：

前光滑：

(12)

粗網(wǎng)格問題：

2)計算由細網(wǎng)格到粗網(wǎng)格的梯度校正，

(13)

對粗網(wǎng)格最小化問題應用MGOPT算法的一次循環(huán)，

粗網(wǎng)格到細網(wǎng)格的最小化：

2)在方向e上進行一次線性查找來求得步長的值αk；

后光滑：

在Ωk上再對方程▽Hk(uk)=fk應用γ2次平滑算法迭代，得到

3 數(shù)值試驗與分析

為驗證集體平滑多重網(wǎng)格方法和多重網(wǎng)格最優(yōu)化方法的收斂速度及其穩(wěn)定性，給出了如下數(shù)值算例．見圖1、2.

算例 測試函數(shù)f(x,y)，選擇目標函數(shù)z，使其滿足模型問題：

-Δy-u=f,

y|?Ω=0.

假設(shè)Ω=(0,1)×(0,1)和f,z∈L2(Ω)，

f(x,y)=sin(x2+y2),

z(x,y)=

圖1 測試函數(shù)f(x,y)

圖2 目標函數(shù) z

數(shù)值試驗中，前光滑次數(shù)和后光滑次數(shù)統(tǒng)一選擇為2次，即γ1=γ2=2，以Gauss-Seidel迭代法作為光滑算子，M為細層網(wǎng)格剖分數(shù)．取迭代終止條件為：tol=10-8，Iter為集體平滑多重網(wǎng)格方法的迭代次數(shù)，γ次迭代后，集體平滑多重網(wǎng)格方法的收斂因子記為

線性橢圓型的控制問題的數(shù)值試驗結(jié)果可由表1，2給出．

表1 利用集體平滑方法求解得到的數(shù)值結(jié)果

根據(jù)表1中的數(shù)據(jù)，時間單位為s，可以總結(jié)出以下的結(jié)論：

1)在多重網(wǎng)格循環(huán)次數(shù)相同的前提下，多重網(wǎng)格剖分的數(shù)目對運行的時間的影響相對較大，即影響運行時間的主要因素是求解區(qū)域剖分疏密程度，網(wǎng)格剖分越密，運行的時間越長；

2) 當參數(shù)取不同的值時，集體平滑多重網(wǎng)格方法都是在經(jīng)過7次迭代之內(nèi)收斂的；

3)集體平滑多重網(wǎng)格方法的迭代次數(shù)跟參數(shù)的取值和網(wǎng)格的疏密程度是互不影響的；

4)在相同的網(wǎng)格疏密程度的前提下，集體平滑多重網(wǎng)格的循環(huán)次數(shù)對運行時間的影響相對較小，即增大多重網(wǎng)格的循環(huán)次數(shù)時，運行的時間并沒有大幅度變化，而是在一個很小的范圍內(nèi)浮動．

其中GM是梯度法，MGM是梯度格式的多重網(wǎng)格最優(yōu)化方法，NCG表示非線性共軛梯度法，MNCG表示非線性共軛梯度格式的多重網(wǎng)格最優(yōu)化方法，表2表示利用以上幾種求解線性橢圓優(yōu)化控制問題所耗用的CPU的時間，時間單位為s．

對于單重網(wǎng)格最優(yōu)化方法，可利用梯度法和非線性共軛梯度法來進行求解，并且，根據(jù)表2中的數(shù)據(jù)，我們可以得到以下的結(jié)論：

1)當參數(shù)v的取值減小時，在對應相同的網(wǎng)格疏密程度的前提下，梯度法和梯度格式的多重網(wǎng)格最優(yōu)化方法運行所需的時間是隨著網(wǎng)格規(guī)模的增加而增加，并且，同時使用非線性共軛梯度法和非線性共軛梯度格式的多重網(wǎng)格最優(yōu)化算法也具有這一特性；

2)當參數(shù)v的值確定時，相對應的網(wǎng)格剖分的疏密程度不相同，使用非線性共軛梯度格式的多重網(wǎng)格最優(yōu)化算法都比使用梯度格式的多重網(wǎng)格最優(yōu)化算法的收斂速度要快；

3)在表2中的這幾種情形，多重網(wǎng)格最優(yōu)化方法都極大地提高了單重網(wǎng)格最優(yōu)化格式的獨立性，并且根據(jù)表1相應的數(shù)據(jù)結(jié)果，多重網(wǎng)格最優(yōu)化方法在網(wǎng)格剖分的疏密程度、收斂速度等因素上，跟集體平滑多重網(wǎng)格方法的數(shù)據(jù)結(jié)果有一定的可比性．

表2 利用多重網(wǎng)格最優(yōu)化方法得到的數(shù)值結(jié)果

4 結(jié) 語

本文主要探討了求解線性橢圓型優(yōu)化控制問題的兩種多重網(wǎng)格法：集體平滑多重網(wǎng)格方法和多重網(wǎng)格最優(yōu)化方法．由數(shù)值試驗結(jié)果表明，相對于多重網(wǎng)格最優(yōu)化方法，集體平滑多重網(wǎng)格方法的求解速度更快和具有更好的穩(wěn)定性，同時，這兩種方法都說明了網(wǎng)格疏密程度的獨立性和參數(shù)收斂的獨立性，此外，對于單重網(wǎng)格優(yōu)化問題，應用最優(yōu)化問題多重網(wǎng)格方法可得到更為準確的數(shù)值結(jié)果；應用集體平滑多重網(wǎng)格方法，則需要對問題設(shè)計一個具體的集體平滑格式，再根據(jù)相應的迭代法進行運算，而最優(yōu)化問題多重網(wǎng)格方法則只需要利用相應的迭代法進行操作．

參考文獻：

[1] BORZ`1 A. Multigrid methods for parabolic distributed optimal control problems [J]. J Comput Appl Math , 2003, 157 (2): 365-382.

[2] BORZ`1 A. A multigrid scheme for elliptic constrained optimal control problems [J]. Comput Optim Appl, 2005, 31(3): 309-333.

[3] NOCEDAL J, WRIGHT S J. Numerical optimization [M]. New York: Springe, 1999.

[4] BORZ`1 A. High-order discretization andmultigrid solution of elliptic nonlinear constrained optimal control problems [J]. J Comput Appl Math, 2007, 200(1): 67-85.

[5] BORZ`1 A, KUNISCH K, KWAK D. Accuracy and convergence properties of the finite difference multigrid solution of an optimal control optimality system [J]. SIAM J Control Optim, 2003, 41(5): 1477-1497.

[6] BORZ`1 A. On the convergence of the MG/OPT method [J]. PAMM, 2005, 5(1): 735-736.

[7] OH S, BOUMAN C,WEBB K. Multigrid tomographic inversion with variable resolution data and image spaces [J]. IEEE Trans Image Proc, 2006, 15(9): 2805-2819.

[8] OH S, MILSTEIN A, BOUMAN C,etal. A general framework for nonlinearmultigrid inversion [J]. IEEE Trans Image Proc, 2005. 14(1): 125-140.

[9] LEWIS R M, NASH S G. A multigrid approach to the optimization of systems governed by differential equations [C]//Long Beach: Symposium on Multidisciplinary Analysis and Optimization, CA, 2000.

[10] LEWIS R M, NASH, S G. Model problems for the multigrid optimization of systems governed by differential equations [J]. SIAM J Sci Comput, 2005, 26(6): 1811-1837.

[11] NASH S G. A multigrid approach to discretized optimization problems [J]. Optim Methods Softw, 2000, 14(2): 99-116.

[12] OH S, MILSTEIN A, BOUMAN C,etal. Multigrid algorithms for optimization and inverse problems [J]. IEEE Trans Image Proc, 2005, 14(2): 125-140.