金輝 ,肖志華 ,祁振中

(1.長江大學信息與數(shù)學學院,湖北 荊州 434023;2.西北大學數(shù)學學院,陜西 西安 710127)

1.引言

在工程應用領域中,許多物理系統(tǒng)或現(xiàn)象都可以用數(shù)學模型來描述,然而大多數(shù)的數(shù)學模型規(guī)模都極其復雜,使得對這些系統(tǒng)的直接仿真模擬和理論分析變得十分困難.因此,有效減小系統(tǒng)規(guī)模和縮短計算仿真時間的需求日益上升.模型降階正是基于此思想,將一個較大規(guī)模的復雜系統(tǒng)轉化為一個近似的較小系統(tǒng),同時保留了原始系統(tǒng)的一些重要性能,例如穩(wěn)定性、結構性和無源性等[1].模型降階方法已經(jīng)普遍用于超大規(guī)模集成電路模擬、計算電磁學、微機電系統(tǒng)和流體力學等領域[2].

模型降階自提出以來,已經(jīng)形成了一系列成熟的降階理論和方法,例如Krylov子空間方法、平衡截斷方法、本征正交分解方法和正交多項式方法等.其中平衡截斷是模型降階最著名的方法之一[3],其基本思想是計算出控制系統(tǒng)的可控Gram矩陣和可觀Gram矩陣,進一步構造合適的映射子空間獲得降階模型.該方法可以構造出保持原始系統(tǒng)穩(wěn)定性的降階模型,并且給出近似先驗誤差.然而其需求解兩個大規(guī)模Lyapunov方程,該過程計算量大且計算時間長,這導致了其應用的局限性.在這種情況下,許多學者作出相應改進,提出了近似平衡截斷方法[4].該方法以一種數(shù)值高效的方式獲得近似平衡系統(tǒng),被廣泛應用于線性系統(tǒng)及非線性系統(tǒng)的模型降階中.

對于線性系統(tǒng),相應的模型降階方法已經(jīng)日趨完善.而實際的物理系統(tǒng)本質是非線性的,許多物理系統(tǒng)都可以通過雙線性系統(tǒng)模型來描述.雙線性系統(tǒng)是一類特殊的非線性系統(tǒng),是連接線性系統(tǒng)和非線性系統(tǒng)之間的橋梁,該系統(tǒng)經(jīng)常出現(xiàn)在非線性電路、熱過程和生態(tài)系統(tǒng)等科學領域[5].關于雙線性系統(tǒng)的模型降階,已經(jīng)發(fā)展了多種方法,例如平衡相關方法、Krylov子空間方法和插值法等[6-8].值得注意的是,Shaker和Tahavori在文[9]中定義了雙線性系統(tǒng)的交叉Gram矩陣.該交叉Gram矩陣同時反映了系統(tǒng)的可控性和可觀性,并且當雙線性系統(tǒng)有界輸入有界輸出(BIBO)穩(wěn)定時,其為廣義Sylvester方程的解.作為雙線性系統(tǒng)的一個特殊子類,K-power雙線性系統(tǒng)的應用也很廣泛,涉及液壓驅動建模、多項式系統(tǒng)、非線性系統(tǒng)辨識等領域[10-11].關于K-power雙線性系統(tǒng)的模型降階方法,在文[10]中,Al-Baiyat和Bettayeb首先證明了K-power雙線性系統(tǒng)具有塊對角可控Gram矩陣和可觀Gram矩陣,并提出了一種簡化求解Lyapunov方程計算的保結構平衡截斷方法.XIE和Syrmos[11]證明了K-power雙系統(tǒng)的雙線性廣義奇異攝動近似(GSPA)可以轉化為σ-交互系統(tǒng)的直接截斷,基于此,提出了一種基于雙線性σ-交互系統(tǒng)分析的方法.WANG和JIANG[12]還成功地將矩匹配方法和最優(yōu)H2方法應用于K-power雙線性系統(tǒng)的模型降階,進一步得到保持原始系統(tǒng)結構特征的降階模型.

本文基于Laguerre函數(shù)和雙線性系統(tǒng)的交叉Gram矩陣,提出了一種K-power雙線性系統(tǒng)的保結構模型降階方法.該方法首先將其等價雙線性系統(tǒng)的交叉Gram矩陣用Laguerre函數(shù)進行展開,并結合其正交性,進一步給出交叉Gram矩陣的低秩因子,最后通過構造K-power雙線性系統(tǒng)各子系統(tǒng)的投影變換得到降階模型.相比于經(jīng)典的平衡截斷方法,此方法不需要直接求解Sylvester方程來獲得交叉Gram矩陣,具有較高的計算效率和靈活度,且具有一定的自適應性.

本文首先介紹Laguerre函數(shù)的相關重要性質,接著詳細給出保結構的K-power雙線性系統(tǒng)降階過程,并給出了相應的降階算法.最后,通過數(shù)值實驗有效驗證了算法的有效性.

2.矩陣指數(shù)函數(shù)的Laguerre函數(shù)展開

矩陣指數(shù)函數(shù)eAt又稱狀態(tài)轉移矩陣,在現(xiàn)代控制理論中,可以將其用于狀態(tài)方程的求解.本節(jié)的目標是利用Laguerre函數(shù)[13]對矩陣指數(shù)函數(shù)進行展開.

類似地,對于含有一個穩(wěn)定矩陣A(若Cn×n的所有特征值嚴格位于復平面的左半平面,則稱A為穩(wěn)定矩陣) 的矩陣指數(shù)函數(shù)eAt,可作下列Laguerre展開[16]

其中Laguerre系數(shù)矩陣Ak滿足

I為單位矩陣.

3.K-power雙線性系統(tǒng)的保結構模型降階算法

考慮如下K-power雙線性系統(tǒng)[1,17]：

雙線性系統(tǒng)(3.3)的交叉Gram矩陣[9]為

其滿足如下廣義的Sylvester方程:

i1,2,···,l,j1,2,···,N ?1.

其中{α1,α2,···,αN}為CF-ADI參數(shù),且Re(αj)<0.假設所有CF-ADI參數(shù)為?α,我們有

因此,除了符號的正負性,低秩因子F與CF-ADI迭代得到的低秩因子Zn相同,即,當所有參數(shù)αj與α相同時,上述方法與ADI方法具有等價性。

由式(3.5)得到交叉Gram矩陣Rl的低秩因子,可以將其應用到等價的雙線性系統(tǒng)來構造投影變換矩陣,進而生成降階模型.但該方法沒有保留K-power雙線性系統(tǒng)自身的結構,而在工程設計和控制應用中,通常要求保留原系統(tǒng)的結構.K-power雙線性系統(tǒng)包含k個子系統(tǒng),且這些子系統(tǒng)依次耦合,其具有保留該系統(tǒng)耦合結構的必要性.為此,通過構造該系統(tǒng)各子系統(tǒng)對應的投影矩陣的思想[19],首先對F和G作如下分塊:

可得系統(tǒng)(3.1)的降階系統(tǒng)如下

上述降階過程可由如下算法1描述：

算法1K-power雙線性系統(tǒng)基于Laguerre函數(shù)的保結構模型降階算法

其中σ(A)為A的特征值.對于更一般的情況,如何選取最優(yōu)的Laguerre參數(shù)α需要進一步研究.數(shù)值結果表明,該最優(yōu)問題對我們算法中參數(shù)α的選取具有指導作用.

降階系統(tǒng)(3.6)可以寫為如下等價的雙線性系統(tǒng)形式

定理3.1[20]對于雙線性系統(tǒng)Σ,如果A是穩(wěn)定的,并且可控Gram矩陣P和可觀Gram矩陣Q均存在,則雙線性系統(tǒng)Σ的H2范數(shù)表示為

對于精確求解雙線性系統(tǒng)可控Gram矩陣P和可觀Gram矩陣Q的平衡截斷方法[6]和某些H2優(yōu)化模型降階方法[17],可以得到誤差系統(tǒng)的Hankel范數(shù)或H2范數(shù).一般而言,對于近似平衡截斷方法,即,通過近似求解或分解系統(tǒng)可控Gram矩陣P和可觀Gram矩陣Q,很難得到與精確平衡截斷方法類似的誤差估計式.盡管如此,值得注意的是,XIANG在文[21]中證明了矩陣指數(shù)函數(shù)正交多項式近似的指數(shù)收斂性,受此啟發(fā),我們可以將該結論推廣應用到(雙)線性系統(tǒng)Gram矩陣的低秩分解的收斂性分析,并由此估計降階系統(tǒng)與原系統(tǒng)的誤差,該問題有待進一步的研究.

4.數(shù)值實驗

下面通過一個數(shù)值算例來驗證上述保結構模型降階算法的有效性.該數(shù)值實驗在Intel(R)Core (TM) i5-10210U CPU (1.60 GHZ) PC上進行,內(nèi)存為16.00GB.同時,選用ode15s函數(shù)在MATLAB R2018b中求解系統(tǒng).

考慮由兩個子系統(tǒng)組成的單輸入單輸出(SISO)K-power雙線性系統(tǒng)[12,17,19],并選取每個子系統(tǒng)的階數(shù)均為3000.該系統(tǒng)的系數(shù)矩陣如下:

其中B1,R3000.通過本文給出的算法1(Algorithm 1)對該系統(tǒng)進行降階,該系統(tǒng)優(yōu)化問題(3.7)的解為α ≈10.數(shù)值結果表明,該α值是有效的.取Laguerre函數(shù)展開項數(shù)N4,誤差限tol10-10,由此自適應得到的兩個降階子系統(tǒng)的階數(shù)均為4.同時,分別以文[22]中的Krylov子空間方法(Krylov),文[19]中的Laguerre正交多項式方法(Laguerre)和文[10]中的平衡截斷方法(BT)對該系統(tǒng)進行降階.在Krylov子空間方法中,取q19匹配該系統(tǒng)第一個子系統(tǒng)的前9階矩,同時選取p22,q24匹配第二個子系統(tǒng)的前4階矩[19],由此得到的兩個降階子系統(tǒng)的階數(shù)為9和8.由Laguerre多項式方法和BT方法分別得到的兩個降階子系統(tǒng)的階數(shù)也為9和8.

表1分別給出了該系統(tǒng)的輸入變量u1(t)10sin(10t+5)和u2(t)sin(4t)e-0.2t時各模型降階方法的CPU運行時間(包括降階的時間和模擬降階系統(tǒng)的時間)、最大相對誤差∥y(t)?yr(t)∥2/∥y(t)∥2和加速比,其中yr(t)為降階系統(tǒng)的輸出響應.

表1 各模型降階方法的CPU運行時間、最大相對誤差和加速比

圖1和圖2分別給出了輸入變量為u1(t)時各降階系統(tǒng)和原始系統(tǒng)的瞬態(tài)響應及其對應的相對誤差.

圖1 輸入為u1時各降階系統(tǒng)的瞬態(tài)響應

圖2 輸入為u1(t)時瞬態(tài)響應的相對誤差

圖3和圖4分別給出了輸入變量為u2(t)時各降階系統(tǒng)和原始系統(tǒng)的瞬態(tài)響應及其對應的相對誤差.

圖4 輸入為u2(t)時瞬態(tài)響應的相對誤差

從模擬結果可以看出,對于該SISO K-power雙線性系統(tǒng)模型,由算法1得到的降階系統(tǒng)對原始系統(tǒng)有很好地瞬態(tài)響應近似效果.在相同的近似精度下,算法1構造的降階系統(tǒng)比BT方法構造的降階系統(tǒng)的階數(shù)小.算法1構造的降階系統(tǒng)比Krylov子空間方法和Laguerre多項式方法構造的降階系統(tǒng)有略好的近似精度,并且算法1得到的降階系統(tǒng)的階數(shù)小于這兩種方法(Krylov,Laguerre)得到的降階系統(tǒng)的階數(shù).此外,算法1構造和模擬降階系統(tǒng)的時間少于其他三種方法(Krylov,Laguerre,BT)所用的時間.這些結果表明,本文提出的算法對于該K-power雙線性系統(tǒng)是有效的.

5.結論

本文提出了一種K-power雙線性系統(tǒng)基于Laguerre函數(shù)的保結構模型降階方法.該方法首先將K-power雙線性系統(tǒng)轉化為等價的雙線性系統(tǒng),再結合Laguerre函數(shù)的正交性,將矩陣指數(shù)函數(shù)通過Laguerre函數(shù)展開,然后計算出雙線性系統(tǒng)交叉Gram矩陣的低秩分解因子,從而構造K-power雙線性系統(tǒng)各子系統(tǒng)相應的投影變換,最后得到保結構的降階模型.通過所提算法得到的降階系統(tǒng)能很好地近似原始系統(tǒng)的瞬態(tài)相應,并且該方法計算高效,具有一定的自適應性.最后,數(shù)值實驗有效驗證了該算法的有效性.