邢秩源，王 麗，蔣耀林，2

(1.新疆大學(xué)數(shù)學(xué)與系統(tǒng)科學(xué)學(xué)院，新疆 烏魯木齊 830046；2.西安交通大學(xué)數(shù)學(xué)與統(tǒng)計學(xué)院，陜西 西安 710049)

1 引言

在眾多諸如控制、力學(xué)、流體機械等工程領(lǐng)域，都涉及復(fù)雜系統(tǒng)的計算機仿真、優(yōu)化與控制，而這類系統(tǒng)的模型大多通過偏微分方程來建立。通過對相應(yīng)的偏微分問題進(jìn)行相當(dāng)精細(xì)的模擬求解，可以為實際工程提供安全、高效的數(shù)據(jù)。隨著高性能計算機的快速發(fā)展，偏微分方程的數(shù)值求解方法，如有限差分法，有限元法、有限體積法、譜方法等得到了普遍應(yīng)用。實際工程模型的復(fù)雜性決定了這些偏微分方程的數(shù)值離散系統(tǒng)的規(guī)模十分龐大，這使得單純采取數(shù)值方法求解復(fù)雜和長時間尺度的偏微分方程具有較大局限性，在有限計算資源條件下往往得不到理想的數(shù)值模擬效果。模型降階方法[1]能夠利用低維模型去逼近原高維模型，獲得較低計算代價下的可接受降階模型，極大地節(jié)約計算資源和提高模擬效率。因此偏微分方程的模型降階方法是很重要的研究課題。

模型降階方法在很好地近似高維系統(tǒng)解的同時還能保持原高維系統(tǒng)的一些重要特征，如穩(wěn)定性、結(jié)構(gòu)特性等。隨著模型降階方法理論的成熟，其在許多工程領(lǐng)域中發(fā)揮著越來越重要的作用，也因其具有高效快速的計算性能被逐漸應(yīng)用到偏微分方程的求解中，如本征正交分解模型降階方法[2，3](POD， Proper Orthogonal Decomposition)、Krylov子空間模型降階方法[4]、平衡截斷模型降階方法[5，6]和黎曼優(yōu)化模型降階方法[7，8]等。POD方法是用低維數(shù)據(jù)信息來近似高維數(shù)據(jù)信息的一類常用的降階基方法，該方法已廣泛應(yīng)用于各類非線性系統(tǒng)。對于含參數(shù)偏微分方程，由于方程系數(shù)及相應(yīng)的解隨著參數(shù)的變化而不斷改變，直接用傳統(tǒng)模型降階方法(如POD)求解并不能達(dá)到期望的求解效果，故對參數(shù)進(jìn)行預(yù)處理是很有必要的。Peter Binev等人通過Greedy算法及Kolmogorov寬度對含參數(shù)偏微分方程進(jìn)行降階求解，并驗證了該方法的收斂速率[9]。Haasdonk等人提出了含參數(shù)偏微分方程的Greedy-POD降階基方法，分析了該算法的收斂性[10]。文獻(xiàn)[11]進(jìn)一步闡述了降階基空間的構(gòu)造，分析了先驗和后驗誤差，并驗證了該方法對含參數(shù)偏微分方程求解的有效性。

本文進(jìn)一步提出含參數(shù)偏微分方程的Greedy-KPOD模型降階方法。該方法建立在Krylov enhanced POD (KPOD)方法[12]和Greedy 參數(shù)優(yōu)化算法取樣的基礎(chǔ)上，結(jié)合Galerkin變分理論及參數(shù)分離方法而提出。同時，給出單邊及雙邊Greedy-KPOD方法的具體算法格式。所提出方法具有原始問題參數(shù)結(jié)構(gòu)和較高的計算精度，且相比傳統(tǒng)Greedy-POD方法具有更少的計算復(fù)雜度。

2 含參數(shù)的Galerkin變分離散系統(tǒng)

設(shè)T>0，Ω∈Rd(d=1，2，3)為連通區(qū)域，邊界?Ω為分段光滑的，μ∈P，其中P?R是參數(shù)空間?？紤]如下含參數(shù)偏微分方程

(1)

(2)

成立。選取標(biāo)準(zhǔn)連續(xù)有限元空間Xh：= span{φ1(x)，φ2(x)，…，φN(x)}，則方程(2)的空間有限元離散的Galerkin逼近為：對?t∈[0，T]，求uh∈Xh?X有

(3)

成立，其中uh0表示初始條件u0(x；μ)在離散空間Xh中的插值。對于任意μ∈，t∈[0，T]，方程(3)的解uh(t，x；μ)可以用有限元基函數(shù)φi(x)， (i=1，2，…，N)來逼近，從而得到方程(3)的N個含參數(shù)微分方程，其簡化的矩陣形式為

(4)

3 參數(shù)系統(tǒng)的Greedy-KPOD模型降階方法

對于參數(shù)微分方程，執(zhí)行Greedy-POD方法時需要生成大量快照點并進(jìn)行矩陣分解，其計算規(guī)模是很大的。為此，提出單邊及雙邊Greedy-KPOD模型降階方法來降低大規(guī)模矩陣分解的代價。

3.1 Greedy 參數(shù)優(yōu)化算法

Greedy算法于70年代早期被引入以來，已被廣泛地應(yīng)用于優(yōu)化等問題[13]。該算法是一種通過迭代來選擇最優(yōu)樣本點的參數(shù)優(yōu)化算法，即在滿足一定條件的情況下，逐步迭代選出最優(yōu)參數(shù)。具體地，在迭代過程中假設(shè)給定參數(shù)集

H={μ1，μ2，…，μJ}，H?P

在H中，近似誤差‖u-Vuh‖Xh與后驗誤差Δn(μ)滿足如下關(guān)系

‖u-Vuh‖Xh≤Δn(μ)，?μ∈H

其中u為原始系統(tǒng)的解，uh為有限元離散解，V為離散函數(shù)空間到連續(xù)函數(shù)空間的投影。每次迭代都在剩余參數(shù)集中選擇使后驗誤差最大時對應(yīng)的參數(shù)，故對?μ∈H，迭代過程中每步都需要計算后驗誤差。

假設(shè)已選出n個參數(shù)μi，i=1，2，…，n，其中n

且

其中0<εg≤1為給定的正定充分小的臨界值。

3.2 參數(shù)系統(tǒng)的KPOD方法

隨著模型降階方法的應(yīng)用越來越廣泛，其理論因?qū)嶋H問題的需求而不斷完善。其中KPOD模型降階方法的核心是基于系統(tǒng)傳遞函數(shù)的矩匹配特性構(gòu)造Krylov子空間，進(jìn)而通過與POD方法結(jié)合構(gòu)造降階所需的變換矩陣。該方法比單純POD方法有更低的計算復(fù)雜度。

為便于討論，僅考慮零初始狀態(tài)及齊次的情況，即u0(x；μ)=0及F(μ)=0。定義系統(tǒng)(4)的輸出Y(t；μ)∈Rp及輸出矩陣C∈Rp×N。則有如下輸入輸出參數(shù)系統(tǒng)

(5)

對于參數(shù)系統(tǒng)(5)，根據(jù)Greedy算法選定的參數(shù)μi，i=1，2，…，n，得到對應(yīng)的原始系統(tǒng)的傳遞函數(shù)

H(s；μi)=C(sM-(μi))-1R(μi)

設(shè)S(μi)非奇異，將H(s；μi)在s=0處Taylor展開為

故原始系統(tǒng)的傳遞函數(shù)H(s；μi)在s=0處的第j階矩為

ml(μi)=-C((μi)-1M)l(μi)-1R(μi)

(6)

類似地，投影系統(tǒng)的傳遞函數(shù)在s=0處的第l階矩為

(7)

下面針對參數(shù)系統(tǒng)(5)，提出如下單邊及雙邊KPOD模型降階方法。

3.2.1 單邊KPOD模型降階方法

對于單邊KPOD模型降階方法，首先利用系統(tǒng)的輸入信息，構(gòu)造輸入Krylov子空間Kk((μi)-1M，(μi)-1R(μi))。其次，通過塊 Arnoldi過程得到一組正交基矩陣=[(μ1)，(μ2)，…，(μn)]∈RN×(pkn)，利用可得如下投影系統(tǒng)

(8)

證明：因為矩陣

(μi)-1M∈Kk((μi)-1M，(μi)-1R(μi))

根據(jù)式(7)和式(8)可以得到

m0(μi)=-CT(T(μi))-1TR(μi)

=-CTξ0=-CT(μi)-1R(μi)=M0(μi)

((μi)-1M)l(μi)-1R(μi)=ξl

成立。為證明此結(jié)論，用歸納法可證明下式成立

((T(μi))-1TM)l(T(μi))-1TR(μi)=ξl

(9)

類似于文獻(xiàn)[3]中的證明，有

ml(μi)=-CT((T(μi))-1TM)l(T(μi))-1TR(μi)

=-CTξl=-CT((μi)-1M)l(μi)-1R(μi)=Ml(μi)故ml(μi)=Ml(μi)，l=0，1，…，k-1。證畢。

3.2.2 雙邊KPOD模型降階方法

雙邊KPOD模型降階方法是利用系統(tǒng)的輸入及輸出信息分別構(gòu)造輸入Krylov子空間Kk((μi)-1M，(μi)-1R(μi))及輸出Krylov子空間Kk((μi)-TMT，(μi)-TCT)。進(jìn)而分別得到輸入、輸出Krylov子空間的正交基矩陣=[(μ1)，(μ2)，…，(μn)]及=[(μ1)，(μ2) …，(μn)]，令V=[，]∈RN×2pkn。

證明：由于矩陣

因此存在一個矩陣ζl∈2kn×p滿足((μi)-TMT)l(μi)-TC=Vζl，l=0，1，…，2k-1等價于進(jìn)而，可以得到

CTV((VT(μi)V)-1VTMV)l

=(CT(μi)-1)(μi)V((VT(μi)V)-1VTMV)l

=CT((μi)-1M)l-1(μi)-1MV

=CT((μi)-1M)lV

此外，根據(jù)定理1可得

V((VT(μi)V)-1VTMV)k-1(VT(μi)V)-1VTR(μi)

因此，當(dāng)l=k，k+1，…，2k-1時，由定理1的證明過程可得ml(μi)=Ml(μi)，l=0，1，…，2k-1，即匹配原始系統(tǒng)的前2k階矩。證畢。

3.3 Greedy-KPOD模型降階方法

在Greedy算法的迭代過程中，每次都要計算剩余參數(shù)集中參數(shù)所對應(yīng)的后驗誤差Δn(μ)，因此降低后驗誤差的計算復(fù)雜度是很有必要的。

在參數(shù)系統(tǒng)(5)中，由于該系統(tǒng)的系數(shù)矩陣S(μ)、L(μ)及R(μ)與參數(shù)μ相關(guān)，隨著參數(shù)的變化每次降階都要在原始大規(guī)模系數(shù)矩陣上進(jìn)行，這無疑增加了降階過程的計算復(fù)雜度。故在降階前，首先對參數(shù)進(jìn)行分離，將含參系數(shù)矩陣通過近似的方式分離為參數(shù)部分與非參數(shù)部分，即

其中Sq∈RQ×Q、Lq∈RQ×Q及Rq∈RQ×Q，q=1，2，…，Q，分別表示S(μ)、L(μ)、R(μ)在點μ0處的近似Taylor展開系數(shù)矩陣。然后將上式代入原參數(shù)系統(tǒng)(5)中，得到近似參數(shù)系統(tǒng)為

(10)

(11)

算法1：單邊Greedy-KPOD模型降階方法

輸入變量：系數(shù)矩陣M，q，Rq，C；參數(shù)樣本集H?；迭代臨界值Nm；誤差臨界值εg；降階階數(shù)r；初始參數(shù)μ1∈。

輸出變量：變換矩陣Ψr和降階系統(tǒng)的系數(shù)矩陣Mr，qr，Rqr，Cr。

1)n=0，δ0=1

2)Dowhilen≤Nm且δn>εg

3) 令n=n+1

B(μn)=-1(μn)R(μn)

5) 對輸入Krylov子空間Kk(A(μn)，B(μn))實施塊Arnoldi算法得基底矩陣∈RN×pkn，對奇異值分解得到正交陣Φ

7)Ψr=[ΨrΦ]，并對Ψr進(jìn)行QR分解得到Γ

8) 截斷Γ的前r列得到Γr，并令Ψr=Γr

9) 令H=H∪{μn}

10)fori=1：length(H)

11) 利用Ψr得到降階解

12) 計算后驗誤差Δn(μi)

13)endfor

16) end while

為了進(jìn)一步利用輸出矩陣信息，在單邊Greedy-KPOD算法基礎(chǔ)上給出下列雙邊算法。

算法2：雙邊Greedy-KPOD模型降階方法

輸入變量：系數(shù)矩陣M，q，Rq，C；參數(shù)樣本集H?P；迭代臨界值Nm；誤差臨界值εg；降階階數(shù)r；初始參數(shù)μ1∈P。

1) 執(zhí)行算法1得到單邊Greedy-KPOD方法的變換矩陣Ψr

2) 類似算法1第5步構(gòu)造輸出Krylov子空間Kk((μi)-TMT，(μi)-TCT)，得到基底矩陣，并對其奇異值分解得到

注意到，雙邊Greedy-KPOD降階系數(shù)矩陣的構(gòu)造將依賴當(dāng)前計算的參數(shù)點，即對不同參數(shù)點降階時需要進(jìn)行上述幾步預(yù)處理，因此計算復(fù)雜度相比單邊算法略大一些。

4 數(shù)值實驗

下面驗證Greedy-KPOD模型降階方法對含參數(shù)偏微分方程的可行性及優(yōu)勢。

設(shè)T∈[0，1]，Ω=(0，1)，P?[0.8，1]，給定40個參數(shù)?？紤]如下零初始狀態(tài)的含參數(shù)偏微分方程

(12)

設(shè)a1(μ)=2μ+1，a2(μ)=5μ+1，邊界條件為[ω1，ω2]=[100，0.1]。

對空間Ω等距剖分，離散自由度為N。通過Galerkin有限元離散，得到方程(12)的含參數(shù)微分方程組，根據(jù)輸出矩陣C，進(jìn)一步得到形如系統(tǒng)(5)的矩陣形式。

在參數(shù)集合P內(nèi)隨機選取參數(shù)μ=0.897，基于算法1和算法2，分別得到N=100，N=500，N=1000階降至r=13階的降階系數(shù)矩陣。進(jìn)而得到降階解與有限元解的相對誤差。同時采用Greedy-POD方法作為對比，在迭代中計算前1/5時間區(qū)間的解構(gòu)成快照集合。圖1、圖2、圖3分別為N=100，500，1000時，單邊Greedy-KPOD，雙邊Greedy-KPOD與Greedy-POD方法所對應(yīng)的降階解與有限元解的相對誤差。由圖1-3可以看出，當(dāng)N較大時，單邊及雙邊Greedy-KPOD方法相比傳統(tǒng)Greedy-POD方法的降階效果優(yōu)勢更明顯。此外，結(jié)合多次計算結(jié)果的精度比較，雙邊Greedy-KPOD方法在一些參數(shù)點處比單邊Greedy-KPOD方法更具精度優(yōu)勢，但單邊Greedy-KPOD方法相比雙邊Greedy-KPOD方法的降階效果更穩(wěn)定。

圖1 N=100時的相對誤差

為了說明降階矩陣生成過程的優(yōu)勢，表1給出了上述情況下三種模型降階方法生成變換矩陣所需的時間。從表1可以看出，隨著原始系統(tǒng)規(guī)模增大，在具備精度優(yōu)勢的前提下，兩種Greedy-KPOD方法相比傳統(tǒng)Greedy-POD方法生成變換矩陣所需時間少很多，且雙邊Greedy-KPOD比單邊Greedy-KPOD方法耗時略多一些。所以本文所提方法在降低計算復(fù)雜度方面更有優(yōu)勢。

圖2 N=500時的相對誤差

圖3 N=1000時的相對誤差

表1 生成變換矩陣所用時間

5 結(jié)論

本文研究了含參數(shù)偏微分方程的Greedy-KPOD模型降階方法，提出了基于塊 Arnoldi 過程的單邊及雙邊 Greedy-KPOD 模型降階算法，主要結(jié)論有：

1)所提方法使得降階系統(tǒng)能夠保持原始系統(tǒng)的參數(shù)結(jié)構(gòu)；

2)數(shù)值算例給出了所提方法對應(yīng)的降階解與有限元解的相對誤差，與傳統(tǒng) Greedy-POD 方法相比具有更好的降階效果；

3)實際應(yīng)用單邊和雙邊 Greedy-POD 方法時，如果在計算時間和穩(wěn)定性上有更高要求，則應(yīng)用單邊 Greedy-KPOD 方法；如果問題的輸出矩陣相對復(fù)雜和重要，可以應(yīng)用雙邊 Greedy-KPOD 方法，從而獲得更好的近似精度；

4)相比傳統(tǒng) Greedy-POD 方法需要生成大規(guī)?？煺站仃嚥⑦M(jìn)行分解方式，基于 KPOD 方法進(jìn)行最優(yōu)參數(shù)迭代，生成變換矩陣的用時更少，因此有望通過改善后應(yīng)用于實際的大規(guī)模對含參數(shù)偏微分方程求解問題。