陳世軍

(陽光學(xué)院，福建 福州 350015)

0 引言

在時(shí)變廣義系統(tǒng)的穩(wěn)定性和具有對稱循環(huán)結(jié)構(gòu)的廣義大系統(tǒng)的最優(yōu)控制中存在求解Riccati矩陣方程(RME)解的問題，尤其是求解廣義Riccati矩陣方程的約束解問題[1-2].RME的一般形式可表示為

(1)

近年來，許多學(xué)者對非線性矩陣方程的數(shù)值解進(jìn)行了研究，并給出了一些有效的算法.例如：Long等[3]證明了求解非對稱Riccati方程迭代算法的收斂性；張凱院等[4]用牛頓算法和修正共軛梯度算法研究了Riccati矩陣方程的對稱解問題；Y.Q.Lin等[5]分析了非線性代數(shù)Riccati方程迭代算法的收斂性；Y.F.Ke等[6]利用交替方向法求解了矩陣方程的非負(fù)解；周富照等[7]等給出了子矩陣約束下矩陣方程的正交投影迭代解.本文借鑒非精確牛頓算法(In-Newton算法)、修正共軛梯度算法(MCG算法)和正交投影迭代算法(OPA算法)原理，建立了一種求廣義方程(1)的對稱和反對稱約束解的迭代算法，并通過數(shù)值算例比較了兩類算法的計(jì)算效率.

1 建立求方程(1)Ω1-2解的非精確Newton算法

為了計(jì)算方便，引入如下記號(hào)：

根據(jù)ψ(X)、φX(Y)和γ(Y)的表達(dá)式和矩陣乘法的性質(zhì)，易推導(dǎo)出

ψ(X+Y)=ψ(X1+Y1,X2+Y2)=ψ(X)+φX(Y)+γ(Y)，

(2)

其中φX(Y)是ψ(X)在點(diǎn)X沿Y方向的Frechet導(dǎo)數(shù).

由Newton算法的基本原理可知，當(dāng)Y1和Y2的范數(shù)較小時(shí)，可以舍去式(2)右端關(guān)于Y1和Y2的高次項(xiàng)γ(Y)，即用式(2)右端的線性部分近似可得ψ(X+Y)≈ψ(X)+φX(Y).設(shè)(X1,X2)∈Ω1-2是方程(1)的近似解，則求方程(1)的一組解(X1,X2)等價(jià)于求校正值(Y1,Y2)∈Ω1-2使ψ(X+Y)=O，并可以把求方程(1)的解X近似為求方程(3)的解(Y1,Y2).方程(3)的表達(dá)式為:

φX(Y)=-ψ(X).

(3)

本文利用Newton算法求方程(1)Ω1-2解的步驟如下：

Step 3 計(jì)算X(k +1)=X(k)+Y(k)，令k∶=k+1，轉(zhuǎn)Step 2.

參照Newton算法，本文給出求方程(1)Ω1-2解的非精確Newton算法如下：

(4)

Step 3 計(jì)算X(k +1)=X(k)+Y(k)，令k∶=k+1，轉(zhuǎn)step 2.

2 建立求方程(3)Ω1-2解和最小二乘解的MCG算法

(5)

其中Ai,Bi,Ci,Di,F∈Rn ×n(i=1,2,3,4,5).本文主要解決下面兩個(gè)問題：

問題Ⅰ設(shè)方程(5)有Ω1-2解，求(Y1,Y2)∈Ω1-2，使(Y1,Y2)滿足方程(5).

2.1 求解問題Ⅰ的MCG算法

為建立求解問題Ⅰ的MCG算法，首先引入如下記號(hào)：

求解問題Ⅰ的MCG算法1的步驟如下:

Step 4 令k∶=k+1，轉(zhuǎn)Step 2.

2.2 求解問題Ⅱ的MCG算法

當(dāng)MCG算法1中斷時(shí)，需要求解問題Ⅱ.首先構(gòu)造與方程(5)等價(jià)的約束正規(guī)矩陣方程，并且把求方程(5)的Ω1-2最小二乘解轉(zhuǎn)化為求其約束正規(guī)矩陣方程的Ω1-2解；然后參照MCG算法1，建立求方程(5)Ω1-2最小二乘解的MCG算法，即求問題Ⅱ的解.引入如下記號(hào)：

定理2[9]求解問題Ⅱ等價(jià)于求矩陣方程f(Y)=N的Ω1-2解，且該矩陣方程一定有Ω1-2解.

參照MCG算法1的原理，建立求方程f(Y)=N的Ω1-2解，即建立求解問題Ⅱ的MCG算法2.MCG算法2的計(jì)算步驟如下：

Step 4 令k∶=k+1，轉(zhuǎn)Step 2.

3 建立求方程(3)Ω1-2解和最小二乘解的OPA算法

MCG算法雖然給出了有限步的收斂結(jié)果，但是沒有給出收斂率的估計(jì).OPA算法也是一種求解矩陣方程約束解的有效算法，其被廣泛運(yùn)用于信號(hào)處理、矩陣方程的約束解求解等問題中.參照OPA算法的基本原理，本文建立的求問題Ⅰ的OPA算法1的步驟如下：

Step 4 令k∶=k+1，轉(zhuǎn)Step 2.

由定理4可知問題Ⅰ不相容，矩陣方程(5)無Ω1-2解.根據(jù)定理2并參照OPA算法1，本文建立的求問題Ⅱ的OPA算法2如下：

Step 4 令k∶=k+1，轉(zhuǎn)Step 2.

4 數(shù)值算例

分別用兩個(gè)數(shù)值算例比較文中給出的In-Newton-MCG算法和In-Newton-OPA算法.參照文獻(xiàn)[4]的計(jì)算方案，本文將In-Newton算法作為外迭代算法，將MCG算法和OPA算法作為內(nèi)迭代算法.In-Newton-MCG算法2和In-Newton-OPA算法2的初始矩陣Y(1)取算法1終止時(shí)的Y(k).本文設(shè)計(jì)的求方程(1)Ω1-2解的計(jì)算方法如下：

In-Newton-MCG算法：

(6)

Step 3 計(jì)算X(k +1)=X(k)+Y(k)，令k∶=k+1，轉(zhuǎn)Step 2.

In-Newton-OPA算法：

Step 3 計(jì)算X(k +1)=X(k)+Y(k)，令k∶=k+1，轉(zhuǎn)Step 2.

例1用本文建立的In-Newton-MCG算法和In-Newton-OPA算法求方程(1)的一組Ω1-2解和最小二乘解，其中系數(shù)矩陣為：

Ei=Fi=Ci j=I(i,j=1,2)，Mi=Ni=I(i=1,2,3,4)，

數(shù)值實(shí)驗(yàn)利用Matlab軟件進(jìn)行.實(shí)驗(yàn)時(shí)取η=0.1，系數(shù)矩陣階數(shù)取n=4，外迭代Newton算法的終止準(zhǔn)則為10-7，MCG算法和OPA算法的終止準(zhǔn)則為10-8.MCG算法和OPA算法的最大迭代次數(shù)為n0=4 999，非精確Newton算法的初始矩陣為X1=4In和X2=On，MCG算法和OPA算法的初始矩陣均為Y1=Y2=On.用t、k0、k1、k2和error分別表示計(jì)算時(shí)間、非精確Newton算法的迭代次數(shù)、MCG算法1和OPA算法1的迭代次數(shù)、MCG算法2和OPA算法2的迭代次數(shù)和誤差范數(shù).兩種算法的計(jì)算結(jié)果如表1所示.

表1 兩種算法求方程(1)Ω1-2解的結(jié)果

當(dāng)MCG算法和OPA算法的初始矩陣取Y1=Y2=On時(shí)，In-Newton-MC G算法和In-Newton-OPA算法均可得到方程(1)的一組解：

上述兩種算法雖然都能求得方程(1)的解，但從表1可以看出，在相同的條件下In-Newton-MCG算法的效率高于In-Newton-OPA算法的效率.

例2修改例1中的部分系數(shù)矩陣，即使Mi=Ni=O(i=2,3,4)，其余系數(shù)矩陣同例1.求方程(1)的一組Ω1-2解.計(jì)算取η=0.9，非精確Newton算法的初始矩陣為X1=X2=On，MCG算法和OPA算法的初始矩陣為Y1=Y2=On，且在計(jì)算實(shí)驗(yàn)過程中依次增加系數(shù)矩陣的階數(shù).兩種算法求方程(1)Ω1-2解的結(jié)果見表2.

表2 兩種算法求方程(1)Ω1-2解的結(jié)果

由表2可知，當(dāng)選擇不同的非精確Newton算法的初始矩陣時(shí)，In-Newton-OPA算法失效，而In-Newton-MCG算法始終有效.因此，在今后的研究中我們將探討In-Newton-OPA算法在什么條件下失效，以及在什么條件下In-Newton-OPA算法比In-Newton-MCG算法有效.