張凱院，寧倩芝，牛婷婷

(西北工業(yè)大學(xué)應(yīng)用數(shù)學(xué)系，陜西 西安 710072)

1 引言

Riccati矩陣方程及其特殊解研究在線性二次優(yōu)化、衛(wèi)星編隊(duì)重構(gòu)與位置保持以及振動(dòng)控制等問題中有重要的應(yīng)用[1～4]。近年來，人們對(duì)離散時(shí)間代數(shù)Riccati矩陣方程DTARME(Discrete Time Algebraic Riccati Matrix Equation)的特殊解問題進(jìn)行了許多研究，并建立了一些有效的算法。例如，Dragan V等人[1]研究了一類DTARME的最大解與穩(wěn)定解以及最小半正定解的存在性問題，并給出了解的極限表達(dá)式；Davies R等人[3]給出了DTARME解矩陣新的上界；El-Sayed S M等人[5]建立了求特殊DTARME最大對(duì)稱正定解的無逆迭代算法，并給出了算法的收斂性結(jié)論。因此，研究DTARME的數(shù)值計(jì)算方法是有意義的。

考慮來源于線性二次優(yōu)化問題的DTARME[1]:

γ1(X)+[γ2(X)+CTD][γ3(X)+

DTD]-1[γ2(X)+CTD]T=CTC

(1)

其中Ai、C∈Rn×n，Bi、D∈Rn×m。

研究DTARME(1)對(duì)稱解的數(shù)值算法，關(guān)鍵是如何轉(zhuǎn)化涉及未知矩陣的求逆計(jì)算問題。本文首先對(duì)DTARME(1)中的逆矩陣采用矩陣級(jí)數(shù)方法進(jìn)行等價(jià)轉(zhuǎn)化，然后運(yùn)用牛頓算法求對(duì)應(yīng)的非線性矩陣方程的對(duì)稱解，并采用修正共軛梯度算法MCG(Modified Conjugate Gradient)求由牛頓算法每一步迭代計(jì)算導(dǎo)出的線性矩陣方程LME(Linear Matrix Equation)的對(duì)稱解或者對(duì)稱最小二乘解(Ls解)，建立求DTARME(1)的對(duì)稱解的雙迭代算法。文中采用的MCG算法不同于通常的共軛梯度法，它不要求涉及的等價(jià)線性代數(shù)方程組的系數(shù)矩陣對(duì)稱正定、可逆或者列滿秩，因此總是可行的。本文建立的雙迭代算法僅要求DTARME(1)有對(duì)稱解，不要求DTARME(1)的對(duì)稱解唯一，也不要求由牛頓算法每一步迭代計(jì)算導(dǎo)出的LME有對(duì)稱解。

2 求DTARME(1)對(duì)稱解的牛頓算法

引進(jìn)矩陣函數(shù)f(X)=[γ3(X)+DTD]-1，γ4(X)=γ2(X)+CTD，令:

ψ(X)=γ1(X)+γ4(X)f(X)[γ4(X)]T-CTC

當(dāng)Y接近于零矩陣時(shí)，可使f(X)γ3(Y)的譜半徑小于1，由此導(dǎo)出:

(2)

□

第1步給定初始矩陣X(1)∈SRn×n，置k:=1。

第2步如果ψ(X(k))=O，停止；否則，求Y(k)∈SRn×n，使得：

(3)

第3步計(jì)算X(k+1)=X(k)+Y(k)，置k:=k+1，轉(zhuǎn)第2步。

需要指出：對(duì)于牛頓算法中的某個(gè)k，當(dāng)LME(3)沒有對(duì)稱解Y(k)時(shí)，可用它的對(duì)稱Ls解Y(k)來代替，這也是本文研究的雙迭代算法的一個(gè)特點(diǎn)。對(duì)于牛頓算法有如下收斂性結(jié)論[6,7]：假設(shè)X*∈SRn×n是DTARME(1)的單根，且初始矩陣X(1)充分接近于X*，那么由牛頓算法確定的矩陣序列{X(k)}收斂于X*。

3 求LME(3)對(duì)稱解與對(duì)稱Ls解的MCG算法

下面建立求LME(3)對(duì)稱解與對(duì)稱Ls解的MCG算法。令Ci=Ai-Bif(X)[γ4(X)]T，當(dāng)Y∈SRn×n時(shí)，LME(3)的一般形式可寫為：

(4)

問題Ⅰ 當(dāng)LME(4)有對(duì)稱解時(shí)，求Y∈SRn×n，使?jié)M足LME(4)。

3.1 求解問題Ⅰ的MCG算法

(5)

記LME(4)和LME(5)在Y=Y(k)處的殘量依次為:

借鑒文獻(xiàn)[8]的基本原理，通過修改有關(guān)矩陣的類型或者算法，建立求解問題Ⅰ的MCG算法(算法1)如下:

第1步給定初始矩陣Y(1)∈SRn×n，置k:=1，計(jì)算:

Rk=F-w(Y(k))

第2步若Rk=O，或者Rk≠O而Zk=O，停止；否則，計(jì)算:

第3步計(jì)算:

Rk+1=F-w(Y(k+1))

第4步置k:=k+1，轉(zhuǎn)第2步。

可以驗(yàn)證，算法1中的矩陣滿足Y(k)、Zk∈SRn×n。對(duì)于算法1有以下收斂性定理(證明過程類似文獻(xiàn)[8])。

定理1設(shè)LME(4)有對(duì)稱解，則對(duì)任意初始矩陣Y(1)∈SRn×n，算法1可在有限步迭代計(jì)算后得到問題Ⅰ的一個(gè)解，即LME(4)的一個(gè)對(duì)稱解。

3.2 求解問題Ⅱ的MCG算法

在算法1中，當(dāng)Rk≠O而Zk=O時(shí)，算法1中斷，這表明LME(4)沒有對(duì)稱解[8]。因此，需要求解問題Ⅱ，即求LME(4)的對(duì)稱Ls解。下面通過構(gòu)造等價(jià)的LME，將求LME(4)的對(duì)稱Ls解問題，轉(zhuǎn)化為求等價(jià)的LME的對(duì)稱解問題。然后參照算法1，建立求LME(4)的對(duì)稱Ls解的迭代算法。引進(jìn)記號(hào)：

其中，Tn,n表示n2階置換矩陣。

定理2求解問題Ⅱ等價(jià)于求LME:

g(Y)=Q

(6)

的對(duì)稱解，且LME(6)一定有對(duì)稱解。

證明當(dāng)Y∈SRn×n時(shí)，YT=Y。因此，求解問題Ⅱ等價(jià)于求Y∈SRn×n，使得：

(7)

下面證明求極小值問題(7)的對(duì)稱解等價(jià)于求LME(6)的對(duì)稱解。將矩陣方程組:

(8)

□

參照算法1及文獻(xiàn)[9]，建立求LME(6)的對(duì)稱解，即求解問題Ⅱ的MCG算法(算法2)如下:

第1步給定初始矩陣Y(1)∈SRn×n，置k:=1，計(jì)算:

第2步若Rk=O，停止；否則，計(jì)算:

第3步計(jì)算:

Rk+1=Q-g(Y(k+1))

第4步置k:=k+1，轉(zhuǎn)第2步。

可以驗(yàn)證，算法2中的矩陣滿足Y(k)、Zk∈SRn×n，對(duì)于算法2有以下收斂性定理(證明過程類似文獻(xiàn)[9])。

定理3LME(6)總是有對(duì)稱解的，對(duì)任意初始矩陣Y(1)∈SRn×n，算法2可在有限步計(jì)算后得到問題Ⅱ的一個(gè)解，即LME(4)的一個(gè)對(duì)稱Ls解。

4 數(shù)值算例

求DTARME(1)的對(duì)稱解，可采用以下兩種計(jì)算方案。

方案1

第1步給定初始矩陣X(1)∈SRn×n，置k:=1。

第2步如果ψ(X(k))=O，停止；否則，采用算法1求Y(k)∈SRn×n，使?jié)M足LME(3)；當(dāng)算法1中斷時(shí)(表明LME(3)無對(duì)稱解)，采用算法2求Y(k)∈SRn×n，使得:

(9)

第3步計(jì)算X(k+1)=X(k)+Y(k)，置k:=k+1，轉(zhuǎn)第2步。

方案2

第1步給定初始矩陣X(1)∈SRn×n，置k:=1。

第2步如果ψ(X(k))=O，停止；否則，采用算法2求Y(k)∈SRn×n，使?jié)M足式(9)。當(dāng)LME(3)有對(duì)稱解時(shí)，它的對(duì)稱Ls解就是它的對(duì)稱解；

第3步計(jì)算X(k+1)=X(k)+Y(k)，置k:=k+1，轉(zhuǎn)第2步。

例1考慮文獻(xiàn)[5]中的例3.2(n=3)，采用雙迭代算法和El-SayedSM算法(文獻(xiàn)[5]算法)求方程X+ATX-1A=I(DTARME(1)的特殊情形)的對(duì)稱解，其系數(shù)矩陣等如下：

根據(jù)文獻(xiàn)[5],取初始矩陣為X(1)=I∈SRn×n，迭代次數(shù)和計(jì)算時(shí)間對(duì)比見表1。

Table 1 First results of example 1

若更換方程X+ATX-1A=I的系數(shù)矩陣等如下：

取初始矩陣為X(1)=I∈SRn×n，迭代次數(shù)和計(jì)算時(shí)間對(duì)比見表2。

Table 2 Second results of example 1

計(jì)算結(jié)果表明：表1中El-SayedSM算法效率較高，而表2中雙迭代算法效率較高，其原因是表2中方程的系數(shù)矩陣不滿足文獻(xiàn)[5]中的要求。

Table 3 Results of example 2(1)

需要指出，本文建立的雙迭代算法不對(duì)DTARME(1)的系數(shù)矩陣做附加限定，不必選取特定的初始矩陣，但求出的只是DTARME(1)的對(duì)稱解。而文獻(xiàn)[5]中El-SayedSM算法對(duì)DTARME(1)的系數(shù)矩陣有要求，選取特定的初始矩陣時(shí)，可求出DTARME(1)的最大對(duì)稱正定解。

(1) 取牛頓算法的初始矩陣為X(1)=4ones(n)-I∈SRn×n，兩種方案的迭代次數(shù)和計(jì)算時(shí)間對(duì)比見表3。

兩種方案求出的對(duì)稱解均為(未知矩陣階數(shù)n=100情形)：

X(5)=0.0399ones(100)-2.4582I+0.2791J

(2) 取牛頓算法的初始矩陣為X(1)=6ones(n)∈SRn×n，兩種方案的迭代次數(shù)和計(jì)算時(shí)間對(duì)比見表4。

兩種方案求出的對(duì)稱解均為(未知矩陣階數(shù)n=100情形)X(3)=0.0181ones(100)。對(duì)比(1)中的計(jì)算結(jié)果可見，該DTARME(1)的對(duì)稱解不唯一。

計(jì)算結(jié)果表明：當(dāng)某一步的LME(3)有對(duì)稱解時(shí)，方案1一般比方案2的效率高，其原因是當(dāng)LME(3)有對(duì)稱解時(shí)，方案2總是采用求LME(3)的對(duì)稱Ls解的算法2進(jìn)行計(jì)算，而算法2比算法1的計(jì)算耗時(shí)多；反之，當(dāng)每一步的LME(3)都沒有對(duì)稱解時(shí)，方案2比方案1(a)的效率高。

Table 4 Results of example 2(2)

5 結(jié)束語

運(yùn)用牛頓算法求DTARME(1)的對(duì)稱解，并采用MCG算法求由牛頓算法每一步迭代計(jì)算導(dǎo)出的LME的對(duì)稱解或者對(duì)稱最小二乘解，建立了求DTARME(1)的對(duì)稱解的雙迭代算法，數(shù)值算例表明雙迭代算法是有效的。修改算法中初始矩陣的類型，算法1中涉及的矩陣Z1與Zk+1的計(jì)算公式，以及算法2中涉及的LME(6)的構(gòu)造方式，還可建立求DTARME(1)的其它特殊解的雙迭代算法。

[1]DraganV,MorozanT.AclassofdiscretetimegeneralizedRiccatiequations[J].JournalofDifferenceEquationsandApplications, 2010, 16(4):291-320.

[2]ParkHE,ParkSY,ChoiKH.Satelliteformationreconfigurationandstation-keepingusingstate-dependentRiccatiequation technique[J]. Aerospace Science and Technology, 2011, 15(6):440-452.

[3] Davies R, Shi P, Wiltshire R. New upper solution bounds of the discrete algebraic Riccati matrix equation[J]. Journal of Computational and Applied Mathematics, 2008, 213(2):307-315.

[4] Wei Yan-ding, Lou Jun-qiang, Lü Yong-gui, et al. Research on linear quadratic optimal control problem in active vibration control[J]. Journal of Zhejiang University (Engineering Science), 2009, 43(3):420-424. (in Chinese)

[5] El-Sayed S M, Al-Dbiban A M. A new inversion free iteration for solving the equationX+AHX-1A=Q[J]. Journal of Computational and Applied Mathematics, 2005, 181(1):148-156.

[6] Higham N J, Kim H M. Solving a quadratic matrix equation by Newton’s method with exact line searches[J]. SIAM Journal on Matrix Analysis and Applications, 2001, 23(2):303-316.

[7] Long J H, Hu X Y, Zhang L. Improved Newton’s method with exact line searches to solve quadratic matrix equation[J]. Journal of Computational and Applied Mathematics, 2008, 22(2):645-654.

[8] Zhang Kai-yuan, Yuan Fei. The modified conjugate gradient method for the symmetric solution of the general linear matrix equation[J]. Numerical Mathematics:a Journal of Chinese Universities, 2011, 33(3):215-224. (in Chinese)

[9] Liu Xiao-min, Zhang Kai-yuan. The MCG method for a different constrained least square solution of two-variable matrix equations[J]. Acta Mathematicae Applicatae Sinica (Chinese Series), 2011, 34(5):938-948. (in Chinese)

附中文參考文獻(xiàn)

[4] 魏燕定, 婁軍強(qiáng), 呂永桂, 等. 振動(dòng)主動(dòng)控制中線性二次型最優(yōu)控制問題研究[J]. 浙江大學(xué)學(xué)報(bào)(工學(xué)版), 2009, 43(3):420-424.

[8] 張凱院, 袁飛.求一般線性矩陣方程對(duì)稱解的修正共軛梯度法[J]. 高等學(xué)校計(jì)算數(shù)學(xué)學(xué)報(bào), 2011,33(3):215-224.

[9] 劉曉敏, 張凱院.雙變量LMEs一種異類約束最小二乘解的MCG算法[J]. 應(yīng)用數(shù)學(xué)學(xué)報(bào), 2011,34(5):938-948.