邱 岳 東

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

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

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

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

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

1 模型問題

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

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

(1)

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

此外，我們選取的目標(biāo)泛函還需滿足一定的偏微分方程約束優(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)格點(diǎn)的數(shù)目為nk，任意函數(shù)在區(qū)域Ωk里進(jìn)行離散化之后，都可用一個(gè)nk維的向量表示．

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

-Δkyk-uk=fk,

-Δkpk+yk=zk,

vuk-pk=0.

(3)

這里，-Δk表示負(fù)五點(diǎn)拉普拉斯算子，x∈Ωk且x=(ihk,jhk)，i,j是網(wǎng)格點(diǎn)的按字典序的下標(biāo)指數(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的四次多項(xiàng)式方程，并且根據(jù)Cardano-Tartaglia公式求得四次多項(xiàng)式方程的最小實(shí)數(shù)值：

(6)

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

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

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

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

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

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

Ak(ωk)=fk,

(7)

前光滑：

(8)

粗網(wǎng)格問題：

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

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

(9)

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

粗網(wǎng)格校正：

后光滑：

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

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

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

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

(10)

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

▽Hk(uk)=fk．

(11)

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

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

前光滑：

(12)

粗網(wǎng)格問題：

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

(13)

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

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

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

后光滑：

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

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

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

算例 測試函數(shù)f(x,y)，選擇目標(biāo)函數(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 目標(biāo)函數(shù) z

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

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

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

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

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

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

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

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

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

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

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

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

3)在表2中的這幾種情形，多重網(wǎng)格最優(yōu)化方法都極大地提高了單重網(wǎng)格最優(yōu)化格式的獨(dú)立性，并且根據(jù)表1相應(yīng)的數(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ù)值試驗(yàn)結(jié)果表明，相對于多重網(wǎng)格最優(yōu)化方法，集體平滑多重網(wǎng)格方法的求解速度更快和具有更好的穩(wěn)定性，同時(shí)，這兩種方法都說明了網(wǎng)格疏密程度的獨(dú)立性和參數(shù)收斂的獨(dú)立性，此外，對于單重網(wǎng)格優(yōu)化問題，應(yīng)用最優(yōu)化問題多重網(wǎng)格方法可得到更為準(zhǔn)確的數(shù)值結(jié)果；應(yīng)用集體平滑多重網(wǎng)格方法，則需要對問題設(shè)計(jì)一個(gè)具體的集體平滑格式，再根據(jù)相應(yīng)的迭代法進(jìn)行運(yùn)算，而最優(yōu)化問題多重網(wǎng)格方法則只需要利用相應(yīng)的迭代法進(jìn)行操作．

參考文獻(xiàn)：

[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.