劉 萌,苗 真,蔣耀林,

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

1 引言

隨著現(xiàn)代集成電路技術(shù)的快速發(fā)展,各類器件大規(guī)模集成化以及信號(hào)頻率不斷提高,大型或復(fù)雜動(dòng)力系統(tǒng)仿真給工程技術(shù)人員帶來了巨大的挑戰(zhàn)。為了解決這個(gè)問題,必須給出新型的仿真方法來降低系統(tǒng)的復(fù)雜程度和分析難度,進(jìn)而減少計(jì)算量,提高計(jì)算效率。這一背景下模型降階方法應(yīng)運(yùn)而生,其主要思想是用一個(gè)小規(guī)模的系統(tǒng)去近似原始大規(guī)模系統(tǒng),降階系統(tǒng)在允許誤差范圍內(nèi)保持系統(tǒng)輸入輸出特性的同時(shí)還要保持原始系統(tǒng)的穩(wěn)定性、無源性等主要性質(zhì)[1,2]。

在模型降階方法成體系出現(xiàn)之前,Padé逼近作為一種近似傳遞函數(shù)的方法,在數(shù)學(xué)領(lǐng)域已經(jīng)被廣泛應(yīng)用。以Padé有理逼近思想為基礎(chǔ),進(jìn)一步產(chǎn)生了漸近波形估計(jì)方法(AWE),它是一種顯式矩匹配方法[3],通過包含相對(duì)少量的主極點(diǎn)和零點(diǎn)的降階模型來近似傳遞函數(shù)。雖然電路可以產(chǎn)生大量的極點(diǎn)和留數(shù),但在特定的頻率范圍內(nèi),所有的極點(diǎn)和留數(shù)的貢獻(xiàn)并不一定相等。瞬態(tài)響應(yīng)通?？梢杂奢^少的主極點(diǎn)和留數(shù)控制,AWE方法可以提取這些主極點(diǎn)和留數(shù),從而構(gòu)建低階的傳遞函數(shù),在合理精度內(nèi)逼近電路的原始響應(yīng)。這種方法只適用于較小規(guī)模的系統(tǒng),其主要考慮降階系統(tǒng)如何保持系統(tǒng)的矩。文獻(xiàn)[4]提出了面向脈沖響應(yīng)矩的計(jì)算方法,對(duì)大規(guī)模系統(tǒng)進(jìn)行矩匹配過程提供了理論基礎(chǔ)。矩匹配類方法的基本思想是使降階后的傳遞函數(shù)匹配原始系統(tǒng)傳遞函數(shù)的前若干矩。針對(duì)單輸入單輸出系統(tǒng),文獻(xiàn)[5]利用矩的概念研究了線性和非線性系統(tǒng)的模型降階問題。文獻(xiàn)[6]提出了一種獲得線性系統(tǒng)傳遞函數(shù)最優(yōu)降階模型的新方法。它以迭代的方式使用多點(diǎn)Padé逼近來高效地生成最優(yōu)模型。該方法通過將多個(gè)點(diǎn)的Padé逼近簡化為等價(jià)的Taylor級(jí)數(shù)來計(jì)算近似值。一維AWE的研究逐漸成熟,文獻(xiàn)[7]提出了M-D AWE方法,該方法將一維的AWE方法推廣至任意m維。一階極點(diǎn)AWE方法無法處理重復(fù)極點(diǎn)問題,文獻(xiàn)[8]將傳統(tǒng)的一階極點(diǎn)方法推廣到高階極點(diǎn)情形,提出了廣義AWE模型降階方法,對(duì)AWE方法進(jìn)行了擴(kuò)充。將有理逼近相關(guān)理論應(yīng)用到模型降階領(lǐng)域當(dāng)中,有助于提高模擬精度,例如,與經(jīng)典Padé逼近相對(duì)照,Chebyshev-Padé逼近具有良好的整體逼近性質(zhì),且非常接近最佳一致有理逼近[9]。

本文在傳統(tǒng)AWE算法的基礎(chǔ)上提出Chebyshev-Padé型極點(diǎn)降階算法,使之能夠處理含重復(fù)極點(diǎn)的復(fù)雜問題,并進(jìn)一步將其應(yīng)用至傳輸線方程中。所提出的基本方法能夠在很少的降階階數(shù)下提高模型的近似精度。

2 Chebyshev-Padé型極點(diǎn)降階算法

本節(jié)主要介紹Chebyshev-Padé型極點(diǎn)降階算法,分別討論單重極點(diǎn)和多重極點(diǎn)情形。

考慮線性時(shí)不變系統(tǒng)

(1)

其中E,A∈Rn×n是常數(shù)矩陣,b∈Rn×1,c∈R1×n,x(t)∈Rn是系統(tǒng)的狀態(tài)變量,u(t)∈R是系統(tǒng)的輸入變量,y(t)∈R是系統(tǒng)的輸出變量。系統(tǒng)(1)經(jīng)Laplace變換后滿足Y(s)=H(s)U(s)=c(sI-A)-1bU(s),其中Y(s)是y(t)的Laplace變換,U(s)是u(t) 的Laplace變換,H(s)表示系統(tǒng)的輸入輸出關(guān)系,即為傳遞函數(shù)?？珊営浵到y(tǒng)(1)為{E,A,b,c}。

2.1 Chebyshev多項(xiàng)式及其性質(zhì)

定義1[9]多項(xiàng)式

Tn(x)=cos(narccosx),-1≤x≤1,n=0,1,…

稱為Chebyshev多項(xiàng)式。

性質(zhì)1[9]Tn(x)=2xTn-1(x)-Tn-2(x),n=2,3,…

T0(x)=1,T1(x)=x

性質(zhì)2[9]Tn是首項(xiàng)系數(shù)為2n-1的n次代數(shù)多項(xiàng)式,且T2k(x) 只含x的偶次冪,T2k-1(x)只含x的偶次冪。

GT(x)=X

(2)

其中

(3)

2.2 Chebyshev-Padé型極點(diǎn)降階

對(duì)于線性時(shí)不變系統(tǒng)(1),將傳遞函數(shù)H(s)在s=0處做Taylor級(jí)數(shù)展開,可得

(4)

將傳遞函數(shù)H(s)在s=∞處進(jìn)行Laurent展開,可得

(5)

(6)

(7)

記Ci=miGi,C-i=m-iG(i=1,2,…)為廣義矩,其中Gi表示G的第i行。

用下述函數(shù)對(duì)原始系統(tǒng)進(jìn)行有理逼近

(T(s))

(8)

(9)

(T(s))=1+b1T1(s)+…+bqTq(s)

(10)

(11)

(12)

上述降階過程可用算法1來描述。

算法1:Chebyshev-Padé型極點(diǎn)降階算法

輸入:線性系統(tǒng){E,A,b,c},降階次數(shù)q

2) 將H(s)轉(zhuǎn)化為Chebyshev多項(xiàng)式

2.3 Chebyshev-Padé型逼近定理

為便于描述,下面簡記矩陣

(13)

(14)

存在的充分必要條件是

rank{HC(L,M-1,M-1)}=rank{HC(L+1,M,M-1)}

(15)

證明:由于

(16)

成立的充分必要條件是方程

QM(T(s))H(T(s))-PL(T(s))=O(s(L+M+1))

(17)

有解,令(17)式兩邊對(duì)應(yīng)系數(shù)相等可得

(18)

進(jìn)一步,方程(18)可改寫為

(19)

方程(18)和(19)的系數(shù)矩陣分別為HC(L+1,M,M-1)和HC(L,M-1,M-1),因此(17)式有解的充要條件為方程組系數(shù)矩陣和增廣矩陣具有相同的秩,即定理結(jié)論成立。

推論1 若矩陣HC(L,M-1,M-1)是非奇異的,則函數(shù)H(T(s))存在有理逼近式

(20)

3 傳輸線系統(tǒng)的降階

本節(jié)將上述算法應(yīng)用于傳輸線系統(tǒng)。傳輸線是一種導(dǎo)行電磁波結(jié)構(gòu),多導(dǎo)體傳輸線(MTL)的應(yīng)用涵蓋了很寬的頻域范圍和各種傳輸線類型,從工頻輸電線路到微波電路。對(duì)于n+1個(gè)平行于z軸的導(dǎo)體構(gòu)成的MTL,其傳輸線方程是一組耦合的2n個(gè)關(guān)于n個(gè)電壓Vi(z,t)和電流Ii(z,t) (i=1,2,…,n)的一階偏微分矩陣方程。大量的導(dǎo)體會(huì)產(chǎn)生巨大數(shù)目的耦合MTL系統(tǒng),為了方便計(jì)算,需要降低結(jié)構(gòu)的表征階數(shù),用較少的有限個(gè)主導(dǎo)極點(diǎn)表示MTL傳輸函數(shù)[10]。

3.1 傳輸線模型

假設(shè)導(dǎo)體為理想導(dǎo)體,包圍在導(dǎo)體周圍的介質(zhì)是均勻的,考慮傳輸線上的一個(gè)微分長度單元,可由單位長度電阻R,單位長度電感L,單位長度電導(dǎo)G,單位長度電容C對(duì)其進(jìn)行描述。對(duì)傳輸線做出以下兩個(gè)假設(shè):1)傳輸線的橫向尺寸在其長度方向上是連續(xù)的;2)傳輸線上的電磁波為準(zhǔn)TEM波,即傳輸線上沿波傳輸方向的電磁分量遠(yuǎn)遠(yuǎn)小于其橫向分量。設(shè)i(x,t)與v(x,t)分別為傳輸線上x處的分布電流和分布電壓,根據(jù)基爾霍夫電壓定律KVL和基爾霍夫電流定律KCL可得如下傳輸線系統(tǒng)方程

(21)

其中R為電阻,L為電感,G為電導(dǎo),C為電容。關(guān)于時(shí)間t做Laplace變換,可得

(22)

其中V(x,s),I(x,s),U(s),Y(s)分別表示,v(x,t),u(t),y(t)經(jīng)過Laplace變換后的函數(shù)。可以將系統(tǒng)寫為如下矩陣形式

(23)

對(duì)傳輸線方程使用有限差分進(jìn)行空間離散,將其轉(zhuǎn)化為狀態(tài)方程形式

(24)

其中

3.2 傳輸線模型的降階方法

對(duì)H(s)在s=0處做Taylor級(jí)數(shù)展開得

(25)

即Mi為系統(tǒng)的矩。對(duì)H(s)在s=∞處做Laurent級(jí)數(shù)展開得

(26)

(27)

其中F(0)=F(s)|s=0是直流響應(yīng),且

(T(s))=a0+a1T(s)+…+aq-1Tq-1(s)

(T(s))=1+b1T(s)+…+bqTq(s)

通過比較F(T(s))(T(s))=F(0)(T(s)) 中的項(xiàng)Ti(s)(i=q,q+1,…,2q-1)的系數(shù),可得到如下線性方程組

(28)

若pi≠pj(i≠j),將pj代入下式,得到關(guān)于kj(j=1,2,…,q)的線性方程組

(29)

否則通過如下線性方程組求解出留數(shù)

4 數(shù)值實(shí)驗(yàn)

為驗(yàn)證本文算法的有效性,下面給出兩個(gè)數(shù)值算例。

例1 考慮如下的線性時(shí)不變系統(tǒng)

(30)

圖2 n=300,q=2時(shí)的系統(tǒng)輸出響應(yīng)對(duì)比

由輸出響應(yīng)對(duì)比圖可以看出n=300時(shí),q=2的降階系統(tǒng)的輸出響應(yīng)與原始系統(tǒng)的輸出響應(yīng)曲線擬合程度更好。為了進(jìn)一步表示算法的降階效果,在圖3、圖4中分別給出了基于Chebyshev-Padé型極點(diǎn)降階算法與傳統(tǒng)AWE的原始系統(tǒng)與降階系統(tǒng)的輸出相對(duì)誤差對(duì)比。

圖3 n=300,q=3時(shí)的系統(tǒng)輸出相對(duì)誤差

圖4 n=300,q=2時(shí)的系統(tǒng)輸出相對(duì)誤差

由相對(duì)誤差響應(yīng)圖可以看出系統(tǒng)在較長時(shí)間模擬時(shí),兩種算法均在q=2時(shí)的相對(duì)誤差更小。相比之下,Chebyshev-Padé型極點(diǎn)降階算法的相對(duì)誤差比傳統(tǒng)AWE更小。

例2 考慮一根半徑為rw=10mils,長l=1 m的導(dǎo)線位于無限大接地平面以上,距地高度h=2 cm。導(dǎo)線的電阻為5Ω,電感為L=3*10-4μH,電容為C=2.5*10-3μF。通過空間離散可得600階的傳輸線系統(tǒng)。分別采用Chebyshev-Padé型極點(diǎn)降階算法與傳統(tǒng)AWE,將上述傳輸線系統(tǒng)降至2階,得到如下結(jié)果。

圖5和圖6分別顯示了該傳輸線系統(tǒng)的傳遞函數(shù)在頻域上的幅值和相對(duì)誤差,可以看出將600階的傳輸線系統(tǒng)降至2階時(shí),Chebyshev-Padé型極點(diǎn)降階算法的傳遞函數(shù)頻域幅值擬合效果良好,且擬合度相較AWE更高,具有明顯優(yōu)勢(shì)。

圖5 n=600,q=2時(shí)的傳遞函數(shù)幅值對(duì)比

圖6 n=600,q=2時(shí)的傳遞函數(shù)幅值相對(duì)誤差

5 結(jié)論

本文提出了Chebyshev-Padé型極點(diǎn)降階算法,通過將傳遞函數(shù)展開為Chebyshev級(jí)數(shù),可以在匹配較少的廣義矩的條件下得到更為精確的降階系統(tǒng)。此外,根據(jù)傳遞函數(shù)的極點(diǎn)重?cái)?shù)的不同情形可以通過Chebyshev-Padé型極點(diǎn)降階算法得到對(duì)應(yīng)的留數(shù),并進(jìn)一步獲得傳遞函數(shù)和輸出響應(yīng)。將該算法應(yīng)用于傳輸線模型中,其數(shù)值結(jié)果表明降階系統(tǒng)能夠較好地近似原始系統(tǒng)。