求解大規(guī)模線性時不變系統(tǒng)的最優(yōu)2模型降階問題的共軛梯度法*

曾泰山，魯春元，陳 劍

(1.中山大學(xué)數(shù)學(xué)與計(jì)算科學(xué)學(xué)院，廣東廣州510275;2.廣東藥學(xué)院醫(yī)藥信息工程學(xué)院，廣東廣州510006;3.佛山科學(xué)技術(shù)學(xué)院數(shù)學(xué)系，廣東佛山528000)

曾泰山1，魯春元2，陳 劍3

(1.中山大學(xué)數(shù)學(xué)與計(jì)算科學(xué)學(xué)院，廣東廣州510275;2.廣東藥學(xué)院醫(yī)藥信息工程學(xué)院，廣東廣州510006;3.佛山科學(xué)技術(shù)學(xué)院數(shù)學(xué)系，廣東佛山528000)

針對最優(yōu)2模型降階問題，提出了適合大規(guī)模多輸入多輸出系統(tǒng)的共軛梯度法。該方法僅需利用一階導(dǎo)數(shù)信息，存儲量少，計(jì)算復(fù)雜度低，且具有超線性收斂性。實(shí)驗(yàn)結(jié)果顯示了算法的有效性。

模型降階;共軛梯度法;Grassmann流形;線性時+不變系統(tǒng)

近年來，模型降階越來越受到重視。它可以大大降低大規(guī)模系統(tǒng)模擬和控制中的時間復(fù)雜度，在超大規(guī)模集成電路設(shè)計(jì)，實(shí)時控制，天氣預(yù)報(bào)等領(lǐng)域有著重要應(yīng)用。關(guān)于模型降階的綜述，參見文獻(xiàn)[1-2]。

給定矩陣 A∈Rn×n，B∈Rn×p和 C∈Rq×n，線性時不變系統(tǒng)可以表示成:

其中t≥0表示時間變量，u(t)∈Rp表示輸入，y(t)∈Rq表示輸出，x(t)∈Rn表示系統(tǒng)的狀態(tài)。n是系統(tǒng)的階。p和q是系統(tǒng)的輸入和輸出個數(shù)。當(dāng)系統(tǒng)的階n很大的時候，由于需要大規(guī)模計(jì)算和存儲，這使得系統(tǒng)的數(shù)值模擬和控制非常困難。解決這一問題的一種關(guān)鍵方法是構(gòu)造一個低階的系統(tǒng)去逼近原始的高階系統(tǒng)，并保持相應(yīng)的系統(tǒng)特性。

由于正交投影能更好地保持系統(tǒng)的穩(wěn)定性，本文用正交投影的辦法構(gòu)造降階系統(tǒng)。假設(shè)m=n?；谡煌队暗哪Ｐ徒惦A是選取一個正交矩陣U∈Rn×m，使得UTU=I。我們記A^=UTAU，B^=UTB，C^=CU。通過矩陣U的投影，我們可以獲得如下的降階系統(tǒng):

定義全階系統(tǒng)(A，B，C)的傳遞函數(shù)為

方程(3)和(4)描述的降階系統(tǒng)的傳遞函數(shù)為傳遞函數(shù)G(s)的2范數(shù)的平方定義為矩陣積分的跡[9]

定義代價函數(shù)

其中，Om表示所有m×m正交矩陣構(gòu)成的正交群，U∈Rn×m滿足UTU=I。由定義可知，等價類[U]中的矩陣的列向量張成的線性子空間相同。在實(shí)際計(jì)算中，我們將選擇其中一個正交矩陣U∈Rn×m來代表整個等價類。關(guān)于Grassmann流形的幾何性質(zhì)，參看文獻(xiàn)[10-11]。

由方程(6)以及代價函數(shù)J(U)的定義 (7)可以看出

也就是說，代價函數(shù)J(U)只依賴于U的列向量張成的線性子空間。這意味著代價函數(shù)J(U)實(shí)際上是Grassmann流形上的光滑函數(shù)。因此，最優(yōu)2模型降階問題可以改寫為Grassmann流形上的最小化問題

通過利用Grassmann流形的幾何性質(zhì)，我們將提出求解問題(9)的數(shù)值算法。

在(6)定義的誤差傳遞函數(shù)Ge(s)可改寫為

實(shí)際上，矩陣Ae，Be和Ce定義了一個系統(tǒng)實(shí)現(xiàn)為(Ae，Be，Ce)的誤差系統(tǒng)。Ge稱為是誤差系統(tǒng)的傳遞函數(shù)。誤差系統(tǒng)的可控性Gramian矩陣Ec和可觀測性Gramian矩陣Eo由如下的Lyapunov方程所定義

我們將可控性Gramian矩陣Ec和可觀測性Gramian矩陣Eo做如下分劃co

Lyapunov方程(10)和(11)可以轉(zhuǎn)化為

代價函數(shù) J(U)可以用 Gramian矩陣 Ec表示如下[9]

我們定義n×m的矩陣JU為代價函數(shù)J(U)相對于矩陣U∈Rn×m里面的元素的偏導(dǎo)數(shù):i=1，2，…，n，j=1，2，…，m 。由文獻(xiàn) [7]，偏導(dǎo)數(shù)JU可以由以下方式給出:假設(shè)P，Q，X和Y是方程(14)、(15)、(16)和 (17)的解。假設(shè)R定義為

代價函數(shù)J(U)在點(diǎn)U的偏導(dǎo)數(shù)為

[11]，在點(diǎn)[U]∈Gr(n，m)處的切空間 T[U]Gr(n，m) 可以表示為

代價函數(shù) J的梯度 ▽J∈ T[U]Gr(n，m) 為

2 共軛梯度法

下面我們給出共軛梯度算法的概要。假設(shè)(A，B，C)是全階系統(tǒng)的狀態(tài)空間實(shí)現(xiàn)。記Uk為第k步獲得的模型降階投影矩陣，k=0，1，…。在共軛梯度算法中，通過以Uk為起點(diǎn)，以Fk為方向的測地線上搜索得到下一個投影矩陣Uk+1。流形上的測地線流形上兩點(diǎn)最短路徑。以Uk為起點(diǎn)，方向?yàn)镕k的Grassmann流形的測地線為

其中 Fk=WkΛkVTk是 Fk奇異值分解，Wk∈ Rn×m，Vk∈ Rm×m，Λk=diag(λ1k，λ2k，…，λmk) 。選取一個合適的步長tk≥0，我們可以構(gòu)造下一步的模型降階投影矩陣Uk+1，

為了計(jì)算下一步的搜索方向Fk+1，首先計(jì)算在點(diǎn)Uk+1處的梯度Gk+1=▽J(Uk+1)。通過求解以下的Sylvester方程來計(jì)算矩陣Pk+1，Qk+1，Xk+1和Yk+1，

矩陣Rk+1為

在點(diǎn)Uk+1偏導(dǎo)數(shù)JUk+1為JUk+1=2Rk+1。代價函數(shù)J在點(diǎn)[Uk+1]∈Gr(n，m)的梯度為

下面將介紹共軛梯度法的搜索方向構(gòu)造。計(jì)算舊搜索方向Fk的平行移動

選取新的搜索方向?yàn)楣曹椞荻确较?，即舊搜索方向的平行移動和新的梯度的線性組合

在式子 (28)中，τFk為切向量Fk的沿著測地線的平行移動。切向量在流形上的平行移動是歐氏空間中向量平移的推廣。類似于歐氏空間中共軛梯度法，組合系數(shù)γk可以通過下面的式子 (Polak-Ribiere法則)得到

其中τGk是梯度Gk的平行移動

值得指出的是，如果令γk=0，則Fk+1=－Gk+1，即搜索方向?yàn)樨?fù)梯度方向，此時共軛梯度法就退化為文獻(xiàn)[7]中的梯度流算法。

下面是本文提出的共軛梯度法 (CGA)的主要框架。

Algorithm 1:共軛梯度法 (CGA)

1) 初始化:選取U0∈Rn×m，UT0U0=I，計(jì)算G0=▽J(U0)。置F0=－G0;

2)For k=0，1，2，…，N － 1

a)計(jì)算Fk的奇異值分解Fk=WkΛkVTk;

b)求解極小化問題

其中:Uk(t) =UkVkcos(tΛk)VTk+Wksin(tΛk)VTk;

令 tk=tmin，Uk+1=Uk(tk);

c)計(jì)算梯度Gk+1=▽J(Uk+1);

d)平行移動切向量Fk和Gk到點(diǎn)[Uk+1]:

計(jì)算新的搜索方向

3) 獲得投影矩陣U=UN，計(jì)算

在上述算法中，需要求解極小化問 (30)，這可以通過不完全線性搜索算法來求得，例如Armijo搜索方法，可參見文獻(xiàn)[10]。

對于Grassmann流形上的最優(yōu)化問題，共軛梯度法在滿足一定的條件下達(dá)到超線性收斂[10-11]。因此，本文提出的共軛梯度算法也能達(dá)到超線性收斂。

另外，由于本文所提出的共軛梯度算法不需要求解大規(guī)模的 Lyapunov方程，因此計(jì)算量大大減少。下面，我們來分析共軛梯度算法的計(jì)算復(fù)雜性。在本文中，將以乘法次數(shù)來衡量計(jì)算復(fù)雜性。記Ns極小化問題 (30)的最大搜索步數(shù)。記N為最大迭代次數(shù)。

定理1假設(shè)A∈Rn×n是個具有N(A)個非零元素的稀疏矩陣。假設(shè)存在一個只需要O(rn+N(A))次乘法運(yùn)算的線性方程求解方法求解方程其中r是個固定的整數(shù)。則共軛梯度法CGA計(jì)算復(fù)雜度為證明 在第k步迭代中，k=0，1，…，算法CGA的計(jì)算花費(fèi)主要來自以下三個方面:

第一部分是Pk、Qk、Xk和Yk的計(jì)算，需要求解Sylvester方程 (23)-(26)。在文獻(xiàn) [7]中指出，如果存在一個只需要O(rn+N(A))次乘法運(yùn)算的方法求解方程(31)，則利用文獻(xiàn) [12]中的方法求解Sylvester方程 (23)-(26)的計(jì)算復(fù)雜度只需第二部分是搜索方向Fk的計(jì)算。對給定的A，B，C和Uk，需要O(mN(A)+nm2+nmp+nmq)次乘法運(yùn)算來計(jì)算Rk。因此，它需要O(mN(A)+nm2+nmp+nmq)次乘法運(yùn)算來計(jì)算Fk。

第三部分是非精確搜索方法求解最小化問題(30)以獲得步長tmin。對極小化問題(30)的每一個搜索步，我們需要計(jì)算測地線方程(21)。易知，測地線的計(jì)算復(fù)雜度為O(nm2)。因?yàn)樽畲蟮乃阉鞑綌?shù)為Ns，所以非精確搜索方法求解最小化問題(30)的計(jì)算復(fù)雜度為O(Nsnm2)。

總而言之，算法CGA每一步迭代需要O(nmr+mN(A)+Nsnm2+nmp+nmq)次乘法次數(shù)。由于算法CGA的最大迭代次數(shù)為N，因此算法的總的復(fù)雜度為如果降階系統(tǒng)的階m遠(yuǎn)遠(yuǎn)小于原始大規(guī)模系統(tǒng)的階n，且系統(tǒng)的輸入輸出個數(shù)p和q也遠(yuǎn)遠(yuǎn)小于n，此時共軛梯度法的計(jì)算復(fù)雜度為線性O(shè)(Nn)，其中N為迭代次數(shù)。

3 數(shù)值算例

在本節(jié)，我們測試比較了本文提出的共軛梯度法的有效性。

例1考慮在區(qū)域Ω =(0，1)2上的熱傳導(dǎo)方程。熱傳導(dǎo)方程可以有如下形式

其中 u=u(t，x，y) ，(x，y) ∈ Ω ，t∈[0，∞) 。假設(shè)微分方程在空間域上等距剖分，格點(diǎn)數(shù)為d×d。導(dǎo)出的剛度矩陣A∈Rd2×d2是稀疏的、穩(wěn)定的，帶寬為d。系統(tǒng)的階為n=d2。假設(shè)b1∈Rn是一個所有元素都為1的向量。b2∈Rn是個隨機(jī)向量。假設(shè)B=[b1，b2]，C=BT。此時，構(gòu)造的系統(tǒng)(A，B，C)是個多輸入多輸出系統(tǒng)。

我們將使用本文提出的共軛梯度法 CGA將此大規(guī)模多輸入多輸出系統(tǒng)降階。我們將與以下的方法比較我們的結(jié)果:平衡截?cái)喾椒?BTM[9];快速梯度流算法 FGFA[7]。表1給出了數(shù)值結(jié)果的比較，2相對誤差定義為表中(n，m)給出了系統(tǒng)的階，其中n是原始大規(guī)模系統(tǒng)的階，m是降階系統(tǒng)的階。從表中可以看到，共軛梯度法的計(jì)算結(jié)果最好。

表1 2相對誤差比較Table 1 Comparison of Relative 2error

表1 2相對誤差比較Table 1 Comparison of Relative 2error

BTM FGFA CGA(n，m)6.08 e-3 4.05 e-3 4.03 e-3(1600，3) 1.22 e-2 5.96 e-3 5.88 e-3(3600，3)(900，3)1.17 e-2 7.25 e-3 7.10 e-3

為了觀察算法的收斂情況，圖1給出n=900的時算法FGFA和CGA的收斂曲線。在每一步迭代中，我們計(jì)算了2相對誤差。從圖1中可以看到，兩個算法FGFA，CGA的相對誤差都隨著迭代下降，而算法CGA的收斂速度達(dá)到了超線性。

圖1 收斂曲線比較Fig.1 Comparison of convergence curves

4 總 結(jié)

參考文獻(xiàn):

[1]ANTOULAS A C．Approximation of large-scale dynamical systems[M]．Adv Des Control 6 SIAM，Philadelphia，2005．

[2]GUGERCIN S，ANTOULAS A C．A survey of model reduction by balanced truncation and some new results[J]．International Journal of Control，2004，77(8):748–766．

[3]HUANG X X，YAN W Y，TEO K．Linear-optimal model reduction [J]．IEEE Trans Automat Contr，2001，46:1279–1284．

[4]HYLAND D C，BERNSTEIN D S．The optimal projection equations for model reduction and the relationships among the methods of wilson，skelton and moore[J]．IEEE Trans Automat Contr，1985，30:1201–1211．

[5]LEPSCHY A，MIAN G A，PINATO G，et al．Rational approximation:a non-gradient algorithm[C]//Proc 30th IEEE CDC，1991:2321–2323．

[6]YAN W Y，LAM J．An approximate approach to optimal model reduction[J]．IEEE Trans Automat Contr，1999，44:1341–1358．

[8]GUGERCIN S，ANTOULAS A．C，BEATTIE C．2model reduction for large-scale linear dynamical systems[J]．SIAM J Matrix Anal Appl，2008，30:609–638．

[9]ZHOU K，DOYLE J C，GLOVER K．Robust and optimal control[M]．New Jersey:Prentice-Hall，1996．

[10]ABSIL P A，MAHONY R，SEPULCHRE R．Optimization algorithms on matrix manifolds[M]．Princeton U-niversity Press，2008．

[11]EDELMAN A，ARIAS T A，SMITHS T．The geometry of algorithms with orthogonality constraints[J]．SIAM J Matrix Anal Appl，1999，20:303–353．

[12]VASILYEV D，WHITE J．A more reliable reduction algorithm for behavior model extraction[C]//Proc Int Conf on Computer，Aided Design(ICCAD)，2005:812–819．

A Conjugated Gradient Algorithm for Optimal2Model Reduction of Large Scale Dynamical Systems

ZENG Taishan1，LU Chunyuan2，CHEN Jian3
(1.School of Mathematics and Computational Science，Sun Yat-sen University，Guangzhou 510275，China;2.College of Medical Information Engineering，Guangdong Pharmaceutical University，Guangzhou 510006，China;3.Department of Mathematics，F(xiàn)oshan University，F(xiàn)oshan 528000，China)

A conjugated gradient algorithm with super-linear convergence which is suitable for the optimal2model reduction of the multi-input multi-output large scale dynamical systems is proposed．The proposed algorithm computes only first-order derivative of the cost function．It has low storage requirement and computational cost．Numerical example demonstrates the approximation accuracy and computational efficiency．

model reduction;conjugated gradient method;Grassmann manifold;linear time invariant system

O241.8

A

0529－6579(2011)02－0001－05

2010－04－30

國家自然科學(xué)基金資助項(xiàng)目 (10771224);中山大學(xué)985項(xiàng)目專項(xiàng)基金資助項(xiàng)目

曾泰山 (1981年生)，男，博士生;E-mail:ztszsu@163.com