一類離散時間代數Riccati矩陣方程對稱解的雙迭代算法*

張凱院，寧倩芝，牛婷婷

(西北工業大學應用數學系，陜西 西安 710072)

1 引言

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

考慮來源于線性二次優化問題的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)對稱解的數值算法，關鍵是如何轉化涉及未知矩陣的求逆計算問題。本文首先對DTARME(1)中的逆矩陣采用矩陣級數方法進行等價轉化，然后運用牛頓算法求對應的非線性矩陣方程的對稱解，并采用修正共軛梯度算法MCG(Modified Conjugate Gradient)求由牛頓算法每一步迭代計算導出的線性矩陣方程LME(Linear Matrix Equation)的對稱解或者對稱最小二乘解(Ls解)，建立求DTARME(1)的對稱解的雙迭代算法。文中采用的MCG算法不同于通常的共軛梯度法，它不要求涉及的等價線性代數方程組的系數矩陣對稱正定、可逆或者列滿秩，因此總是可行的。本文建立的雙迭代算法僅要求DTARME(1)有對稱解，不要求DTARME(1)的對稱解唯一，也不要求由牛頓算法每一步迭代計算導出的LME有對稱解。

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

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

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

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

(2)

□

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

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

(3)

第3步計算X(k+1)=X(k)+Y(k)，置k:=k+1，轉第2步。

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

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

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

(4)

問題Ⅰ 當LME(4)有對稱解時，求Y∈SRn×n，使滿足LME(4)。

3.1 求解問題Ⅰ的MCG算法

(5)

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

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

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

Rk=F-w(Y(k))

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

第3步計算:

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

第4步置k:=k+1，轉第2步。

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

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

3.2 求解問題Ⅱ的MCG算法

在算法1中，當Rk≠O而Zk=O時，算法1中斷，這表明LME(4)沒有對稱解[8]。因此，需要求解問題Ⅱ，即求LME(4)的對稱Ls解。下面通過構造等價的LME，將求LME(4)的對稱Ls解問題，轉化為求等價的LME的對稱解問題。然后參照算法1，建立求LME(4)的對稱Ls解的迭代算法。引進記號：

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

定理2求解問題Ⅱ等價于求LME:

g(Y)=Q

(6)

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

證明當Y∈SRn×n時，YT=Y。因此，求解問題Ⅱ等價于求Y∈SRn×n，使得：

(7)

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

(8)

□

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

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

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

第3步計算:

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

第4步置k:=k+1，轉第2步。

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

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

4 數值算例

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

方案1

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

第2步如果ψ(X(k))=O，停止；否則，采用算法1求Y(k)∈SRn×n，使滿足LME(3)；當算法1中斷時(表明LME(3)無對稱解)，采用算法2求Y(k)∈SRn×n，使得:

(9)

第3步計算X(k+1)=X(k)+Y(k)，置k:=k+1，轉第2步。

方案2

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

第2步如果ψ(X(k))=O，停止；否則，采用算法2求Y(k)∈SRn×n，使滿足式(9)。當LME(3)有對稱解時，它的對稱Ls解就是它的對稱解；

第3步計算X(k+1)=X(k)+Y(k)，置k:=k+1，轉第2步。

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

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

Table 1 First results of example 1

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

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

Table 2 Second results of example 1

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

Table 3 Results of example 2(1)

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

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

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

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

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

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

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

Table 4 Results of example 2(2)

5 結束語

運用牛頓算法求DTARME(1)的對稱解，并采用MCG算法求由牛頓算法每一步迭代計算導出的LME的對稱解或者對稱最小二乘解，建立了求DTARME(1)的對稱解的雙迭代算法，數值算例表明雙迭代算法是有效的。修改算法中初始矩陣的類型，算法1中涉及的矩陣Z1與Zk+1的計算公式，以及算法2中涉及的LME(6)的構造方式，還可建立求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)

附中文參考文獻

[4] 魏燕定, 婁軍強, 呂永桂, 等. 振動主動控制中線性二次型最優控制問題研究[J]. 浙江大學學報(工學版), 2009, 43(3):420-424.

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

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