唐升國,蔣耀林2,,王兆鴻

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

1 引言

近年來,在眾多工程應(yīng)用領(lǐng)域,如電子系統(tǒng)、控制系統(tǒng)、微電子系統(tǒng)等,都涉及大型復(fù)雜動力系統(tǒng)的計算機(jī)仿真、優(yōu)化與控制。這些系統(tǒng)一般用微分方程或差分方程描述。隨著人們對系統(tǒng)精確度要求的不斷提高,導(dǎo)致方程的維數(shù)也越來越大,有時甚至可以達(dá)到1011數(shù)量級。然而在實踐中由于受到計算機(jī)有限內(nèi)存空間和運算時間的限制,大規(guī)模系統(tǒng)的分析和仿真有時甚至是無法實現(xiàn)的。因此, 能否在有限計算資源下快速分析和仿真系統(tǒng)成為了制約工業(yè)前期設(shè)計的一個重要因素。模型降階方法能夠有效地解決這些問題[1]。此方法用一個低維系統(tǒng)去近似原始高維系統(tǒng),并且可以保持原始系統(tǒng)的一些重要特性,如穩(wěn)定性、可控性和可觀性等。目前學(xué)者們已經(jīng)提出了許多針對線性定常系統(tǒng)的模型降階方法,如正交多項式模型降階方法[2]、平衡截斷模型降階方法[3,4]和最優(yōu)化模型降階方法[5]等。

線性時變系統(tǒng)作為一類特殊的動力系統(tǒng)已被廣泛研究[6,7],其中關(guān)于線性時變系統(tǒng)模型降階的研究已獲得大量關(guān)注。文獻(xiàn)[8]利用TVP(time varying Padé)算法,將大規(guī)模線性時變系統(tǒng)變?yōu)榫哂邢嗨戚斎胼敵鰝鬏斕匦缘慕惦A系統(tǒng)。文獻(xiàn)[9]提出了離散和連續(xù)時間線性時變系統(tǒng)的平衡截斷模型降階方法,同時還分析了誤差界和保穩(wěn)定性。文獻(xiàn)[10]定義了線性時變對稱系統(tǒng)的交叉Gram矩陣,并且基于可控、可觀和交叉Gram矩陣之間的聯(lián)系,研究了平衡截斷模型降階方法的保對稱性。針對線性周期時變系統(tǒng),文獻(xiàn)[11]利用積分全等變換提出了基于投影的模型降階算法,其中降階系統(tǒng)能夠匹配原始系統(tǒng)一定數(shù)目的矩。文獻(xiàn)[12]借助離散Laguerre多項式,提出了線性離散周期時變系統(tǒng)基于Krylov子空間的模型降階方法。經(jīng)查閱文獻(xiàn)發(fā)現(xiàn),目前關(guān)于線性離散時變系統(tǒng)基于離散正交多項式的時域模型降階方法的研究較少。因此本文研究線性離散時變系統(tǒng)基于離散Walsh函數(shù)的時域模型降階。

針對線性離散時變系統(tǒng),本文提出兩種基于離散Walsh函數(shù)的模型降階方法:①根據(jù)遞推法得到線性離散時變系統(tǒng)狀態(tài)方程的解,借助離散Walsh矩陣的性質(zhì)得到離散Walsh函數(shù)系數(shù)矩陣,用所得系數(shù)矩陣構(gòu)造正交投影矩陣,進(jìn)而得到降階系統(tǒng)。②把系統(tǒng)的各個變量用Walsh函數(shù)展開,從而獲得對應(yīng)的矩陣方程。通過求解矩陣方程得到狀態(tài)變量的離散Walsh函數(shù)系數(shù)矩陣,再結(jié)合狀態(tài)變量的非零初值得到正交投影矩陣并構(gòu)造降階系統(tǒng)。在單側(cè)投影的理論框架下,理論分析表明這兩種方法所得降階系統(tǒng)分別能夠匹配原始系統(tǒng)一定數(shù)目的展開系數(shù)。數(shù)值算例驗證了本文所提方法的可行性和有效性。

2 預(yù)備知識

在本節(jié)中,將介紹離散Walsh函數(shù)的定義和性質(zhì)。文獻(xiàn)[13]給出了關(guān)于k的第i階離散Walsh函數(shù)φi(k)。該函數(shù)被定義在N=2l個點上,其中l(wèi)是非負(fù)整數(shù);i和k是小于N的非負(fù)整數(shù)。令

(i)decimal=(il-1il-2…i0)binary

(k)decimal=(kl-1kl-2…k0)binary

(1)

則離散Walsh函數(shù)定義為

(2)

其中

g0(i)=il-1,g1(i)=il-1+il-2, …,gl(i)=i1+i0

例如,求N=24時的φ3(k),首先將十進(jìn)制數(shù)3表示為四位的二進(jìn)制數(shù)

(3)decimal=(0011)binary,

同時求出

g0(3)=i3=0,g1(3)=i3+i2=0,

g2(3)=i2+i1=1,g3(3)=i1+i0=2。

于是由式(1)和(2)可以得到

φ3(k)=(-1)[g0(3)k0+g1(3)k1+g2(3)k2+g3(3)k3]

即

φ3(k)=[1 1 1 1-1-1-1-1 1 1 1 1-1-1-1-1]。

離散多項式{φi(k),i=0, 1, …,N-1}構(gòu)成一個完備集并滿足正交性

(3)

其中δij是克羅內(nèi)克δ函數(shù)(Kronecker delta)。

給定任意有界離散序列x(k)∈Rn,x(k)可以被展成離散Walsh級數(shù)

(4)

其中

X=[x0x1…xN-1],xi∈Rn,

Φ(k)=[φ0(k)φ1(k) …φN-1(k)]T,

上標(biāo)“T”表示向量或矩陣的轉(zhuǎn)置。利用正交性(3),xi可以被計算出,即

離散Walsh矩陣W定義為[14]

W=[Φ(0)Φ(1) …Φ(N-1)]

離散Walsh函數(shù)具有移位性質(zhì)[15],即存在唯一的矩陣P∈RN×N使得

PΦ(0)=0,PΦ(k+1)=Φ(k)

(5)

成立,其中k=0, 1, …,N-2。事實上,式(5)還滿足PW=WM,其中M被定義為

離散Walsh函數(shù)還具有乘法性質(zhì),它適用于適當(dāng)維數(shù)的向量c,即

Φ(k)Φ(k)Tc=CΦ(k)

(6)

其中矩陣C是一個與向量c相關(guān)的矩陣[16]。進(jìn)一步,當(dāng)k=0, 1, …,N-1時,由式(6)可得[Φ(0)Φ(0)TcΦ(1)Φ(1)Tc…Φ(N-1)Φ(N-1)Tc]=CW,

即

…Φ(N-1)Φ(N-1)Tc]W。

3 線性離散時變系統(tǒng)的離散Walsh函數(shù)模型降階方法

在這一節(jié),將提出兩種基于離散Walsh函數(shù)的模型降階方法?？紤]如下線性離散時變系統(tǒng)

(7)

其中A(k)∈Rn×n,B(k)∈Rn×p,C(k)∈Rq×n是系數(shù)矩陣;x(k)∈Rn是n維狀態(tài)變量,x(0)是初值;u(k)∈Rp和y(k)∈Rq分別是p維輸入變量和q維輸出變量。

本文的主要目的是構(gòu)造正交投影矩陣V∈Rn×(?n),然后在單側(cè)投影框架下構(gòu)造系統(tǒng)(7)的降階系統(tǒng)

(8)

3.1 基于系統(tǒng)狀態(tài)方程解的模型降階方法

采用遞推法,由系統(tǒng)(7)的第一個方程可知

x(k)=A(k-1)A(k-2)…A(0)x(0)+

×A(k-2)…A(j+1)B(j)u(j),

其中k≥2。根據(jù)式(4)可得

[x(0)x(1) …x(N-1)]

=X[Φ(0)Φ(1) …Φ(N-1)],

(9)

從而

(10)

其中W是離散Walsh矩陣。

接下來,利用離散Walsh系數(shù)矩陣X構(gòu)造正交投影矩陣V。需要執(zhí)行Gram-Schmidt(GS)過程來構(gòu)造V的列vi(i=0, 1, …, N-1)。此外,QR分解和改進(jìn)的Gram-Schmidt(MGS)過程也可以用來構(gòu)造正交投影矩陣V。上述降階過程可以總結(jié)為如下算法。

算法1:基于系統(tǒng)狀態(tài)方程解的模型降階方法

輸入變量:離散時變系統(tǒng){A(k),B(k),C(k)}及初值x(0)

1)通過式(10)計算出離散Walsh系數(shù)矩陣X;

2)構(gòu)造投影矩陣V∈Rn×N,它滿足

colspan{X}?colspan{V},VTV=IN;

注 1:在式(10)中, 由于離散Walsh矩陣W∈RN×N非奇異,所以有X∈Rn×N(N?n)和V∈Rn×N(N?n), 因此降階系統(tǒng)的階數(shù)為=N=2l。

注 2:在式(10)中,如果x(0)=0或X不是列滿秩的,則可以選取X的列向量中最大線性無關(guān)組來構(gòu)造正交投影矩陣V∈Rn×N1(N1

3.2 具有系統(tǒng)非零初值的模型降階方法

事實上,在3.1節(jié)中已經(jīng)提出了具有系統(tǒng)非零初值的模型降階方法。在這一節(jié)將研究另外一種具有系統(tǒng)非零初值的模型降階方法。

根據(jù)離散Walsh級數(shù)的性質(zhì),x(k+1)∈Rn可以被展開為

于是有

[x(1)x(2) …x(N)]

(11)

又由于

[x(0)x(1) …x(N-1)]

=[x(1)x(2) …x(N)]M+[x(0) 0 … 0],

所以將式(9)和(11)帶入上式并化簡可得

(12)

其中

X0=x(0)Φ(0)T=[x(0)x(0) …x(0)]∈Rn×N。

根據(jù)離散Walsh級數(shù)的性質(zhì),系統(tǒng)(7)中的各個變量還可以表示為

x(k)=XΦ(k)=(Φ(k)T?In)vec(X),

u(k)=UΦ(k)=(Φ(k)T?Ip)vec(U),U∈Rp×N,

其中?表示Kronecker積,vec(·)表示矩陣的按列拉直算子。用類似得到式(10)的方法, 還可以得到

于是, 將上述各式代入系統(tǒng)(7)的第一個方程中可得

接著應(yīng)用性質(zhì)(6)可得

由于上述各式對k=0,1,…,N-1都成立,所以利用離散Walsh矩陣W,有

從而

(13)

成立。結(jié)合式(12)和(13)可得

(14)

算法2:具有系統(tǒng)非零初值的模型降階方法

輸入變量:離散時變系統(tǒng){A(k),B(k),C(k)}及初值x(0)

2)構(gòu)造正交投影矩陣V∈Rn×,它滿足關(guān)系

3.3 系數(shù)匹配

相似地,原始系統(tǒng)(7)的輸出變量y(k)也可以表示為

于是,由算法1得到的降階系統(tǒng)(8)滿足如下定理。

證明:在算法1中,由于x(k)的離散Walsh函數(shù)系數(shù)矩陣為

(15)

對系統(tǒng)(15)的第一個方程左乘矩陣VT,有

(16)

相似地,由算法2得到的降階系統(tǒng)(8)也滿足如下定理。

定理 2:若V是由算法2得到,則當(dāng)k=0, 1, …,-1時,成立。進(jìn)一步,有X=V和Y=。

所以x(k)?colspan{X}?colspan{V}成立,k=0, 1, …,-1。因為接下來的證明過程與定理1的證明過程類似,所以此處省略。證畢。

4 數(shù)值實驗

在這一節(jié)中,將應(yīng)用一個數(shù)值例子驗證算法1和算法2的可行性和有效性。

考慮如下偏微分方程[12]

(17)

用500個網(wǎng)格,對偏微分方程(17)進(jìn)行有限差分離散,根據(jù)初值條件可以得到499階的連續(xù)時變系統(tǒng)。令p=+∞和時間步長Δt=0.5,利用向后歐拉法可以得到一個n=499的線性離散時變系統(tǒng)。令ξ(t)=1+t,ξ1(t)=t+tsin(t),ξ2(t)=1,z(0,0.002)=0.1, …,z(0,0.6)=0.1,z(0,0.602)=0,…,z(0,0.998)=0。假設(shè)算法1的降階階數(shù)為=8,算法2的降階階數(shù)為=9。當(dāng)u(t)=(-1)2te0.01t時,圖1描述了原始系統(tǒng)和降階系統(tǒng)的輸出;圖2給出了原始系統(tǒng)和降階系統(tǒng)輸出的絕對誤差。從圖1可以看出由算法1和算法2得到的降階系統(tǒng)的輸出能夠很好地逼近原始系統(tǒng)的輸出。在圖2中,可以觀察到由算法1得到的降階系統(tǒng)的輸出能夠匹配原始系統(tǒng)的前8個輸出y(k),然而由算法2得到的降階系統(tǒng)的輸出能夠匹配原始系統(tǒng)的前9個輸出y(k)。因此,由算法1和算法2得到的降階系統(tǒng)分別滿足定理1和定理2的結(jié)論。

圖1 原始系統(tǒng)和降階系統(tǒng)的輸出

圖2 絕對誤差

5 結(jié)論

本文利用離散Walsh函數(shù)研究了線性離散時變系統(tǒng)的時域模型降階,并提出了兩種模型降階方法:①采用遞推法,得到線性離散時變系統(tǒng)的狀態(tài)變量x(k)(k=0, 1, …,-1),并結(jié)合離散Walsh函數(shù)的性質(zhì)得到狀態(tài)變量x(k)的離散Walsh函數(shù)系數(shù)矩陣,將得到的系數(shù)矩陣構(gòu)造正交投影矩陣,進(jìn)而得到降階系統(tǒng)。②將線性離散時變系統(tǒng)在離散Walsh函數(shù)構(gòu)成的空間中展成矩陣方程,通過求解矩陣方程得到狀態(tài)變量的離散Walsh函數(shù)系數(shù)矩陣。利用所得系數(shù)矩陣和非零初值x(0)構(gòu)造正交投影矩陣并得到降階系統(tǒng)。理論分析表明上述兩種方法得到的降階系統(tǒng)分別能夠匹配原始系統(tǒng)一定數(shù)目的展開系數(shù)。數(shù)值算例驗證了本文所提方法的可行性和有效性。