基于最小殘差思想的一類代數(shù)Riccati方程的牛頓迭代求解

周立平

（湖南科技學(xué)院 理學(xué)院，湖南 永州 425199）

基于最小殘差思想的一類代數(shù)Riccati方程的牛頓迭代求解

周立平

（湖南科技學(xué)院 理學(xué)院，湖南 永州 425199）

本文在牛頓迭代法框架下，借用最小殘差法思想，結(jié)合預(yù)條件共軛最小二乘法，提出了一種新的算法對一類代數(shù)Riccati方程的數(shù)值解進(jìn)行了研究，并給出了具體的數(shù)值算例。算例結(jié)果驗(yàn)證了該方法具有良好的收斂性。

代數(shù)Riccati方程（ARE）；最小殘差法；Lyapunov方程；預(yù)條件共軛梯度最小二乘法（PCGLS）；牛頓迭代法

1　引 言

代數(shù) Riccati 方程源于控制理論中的數(shù)值問題，它在非穩(wěn)定系統(tǒng)模型降階，線性二次校正問題和H2優(yōu)化控制中有十分重要的作用。這些問題的解決在一定條件下都可以轉(zhuǎn)化為求（廣義）代數(shù) Riccati方程的穩(wěn)定解。本文旨在研究如下一類（廣義）代數(shù) Riccati 方程（ARE）的數(shù)值解：

求解ARE（1.1）的直接法有schur 向量法[1]和符號(hào)函數(shù)法[2]等，它們主要是基于系數(shù)矩陣的特征分解而提出的。假設(shè)ARE（1.1）的系數(shù)矩陣具有良好的結(jié)構(gòu)，如它們是低秩或稀疏的，若仍采有直接法，則系數(shù)矩矩陣的良好結(jié)構(gòu)特點(diǎn)在計(jì)算中就被白白浪費(fèi)。因此，直接法并不適合用于大規(guī)模稀疏系統(tǒng)的計(jì)算。由于方程（1.1）是非線性的二次矩陣方程，這使得在牛頓迭代法框架下對該類方程進(jìn)行數(shù)值求解成為研究ARE求解的主流方法。

應(yīng)用牛頓迭代法求解ARE，已經(jīng)產(chǎn)生了一系列的求解解法，如Chandrasekhar迭代法[3，4]和NM-Kleinman 迭代法[5，6]等，其主要思想是將二次方程線性化后再求解。Peter Benner等人在文[7]中給出了低秩牛頓Cholesky因子法（LRCF-NM）和低秩牛頓 Galekin-ADI 方法（LRCF-NM-GP-ADI），并對 ARE方程的低秩形式解的理論也進(jìn)行了討論。

這些方法的共同關(guān)健處是求解一系列（或降階）Lyapunov方程。這使得對方程ARE（1.1）的求解是否快速穩(wěn)定主要取決于這些Lyapunov方程的求解實(shí)現(xiàn)。Lyapunov方程的求解一直是數(shù)值代數(shù)專家研究的熱點(diǎn)問題，已經(jīng)取得一系列研究成果[8，9]。最近，Lin與 V.Simoncini 在文[10]中提出了基于最小殘差法（MINRES）求解大規(guī)模Lyapunov 方程的方法，并做了一定的理論分析，取得了很好的效果。

在文[10]的啟發(fā)下，本文旨在牛頓迭代法框架下，借用最小殘差法的思想，利用預(yù)條件共軛最小二乘法（PCGLS）求解Lyapunov方程，進(jìn)而得到ARE（1.1）的數(shù)值逼近解。為方便起見，對將用到的符號(hào)作如下約定：表示矩陣A的Frobenius范數(shù)；I表示相應(yīng)階數(shù)的單位矩陣；若A與B為同類型實(shí)矩陣，則為定義的矩陣A與B內(nèi)積。

一般說來，方程（1.1）不一定可解，即使有解也可能并不唯一。為方便研究，我們對系數(shù)矩陣給出一定的限制條件使其能具有唯一半正定解。幸運(yùn)的是，許多源于實(shí)際問題的ARE 方程其系數(shù)矩陣大多也能滿足該限制條件。下面的引理給出了ARE（1.1）有唯一半正定解的存在條件。

引理1.1[12]設(shè)矩陣對為可穩(wěn)化的，矩陣對為可測的，則ARE（1.1）有唯一半正定解，且該解是穩(wěn)定的。

2　求解ARE（1.1）的新算法

對其求Frechét導(dǎo)數(shù)，得到

從而，得到如下算法。

算法 2.1 （ARE-NM）

Step 3.求解Lyapunov方程的解Nl

注意到算法（2.1）中需要求解一系列Lyapunov 方程（2.2），本文假設(shè)在迭代過程中始終是半正定的（這一假設(shè)在一定條件下是能成立的，可參考文[9]的第3節(jié)）。設(shè)具有分解式

得到如下一般形式的Lyapunov方程

對于方程（2.5），現(xiàn)在最常用的方法是投影空間方法，其思想主要是建立一個(gè)近似空間（搜索空間），利用投影算子將方程投影到該空間得到新的降階矩陣方程，然后對其進(jìn)行求解。如果得到的近似解不能達(dá)到精度，則擴(kuò)充該空間，并重復(fù)上述工作直至收斂。近似空間的一種標(biāo)準(zhǔn)取法是取為Krylov 子空間，如

也可取其擴(kuò)展形式[13]：

為了簡化問題，我們下面以標(biāo)準(zhǔn)Krylov子空間k為代表對算法進(jìn)行設(shè)計(jì)。設(shè)第k次擴(kuò)充時(shí)Krylov子空間為

不失一般性，引入算子

則（2.7）等價(jià)為

下面我們討論（2.8）的求解。為加快收斂速度，我們嘗試采用一些理論成熟的預(yù)條件方法，如采用預(yù)條件共軛梯度最小二乘方法（PCGLS）求（2.8）的解。利用拉直算子可將（2.8）化為矩陣-向量方程形式，為方便起見，我們直接使用其矩陣-矩陣形式。記算子à的伴隨算子為取預(yù)條件子為M，從而，求解 （2.8）的PCGLS算法可描述如下：

算法 2.2 （PCGLS）

結(jié)合算法2.1和2.2，我們得到了如下求解ARE（1.1）的基于最小殘差的預(yù)條件共軛最小二乘牛頓迭代算法。算法 2.3 （ARE-MINRES-PCGLS-NM）

Step 4.利用Arnold算法或分塊Arnold算法[11]擴(kuò)充Krylov子空間kk，得到kn和nk+1；

從算法2.3可以看出，內(nèi)循環(huán)變量k的每次增加都要求用PCGLS方法求解一個(gè)最小殘差F-范數(shù)的優(yōu)化問題。但事實(shí)上，在迭代過程中，我們真正關(guān)心的只是殘差的F-范數(shù)值用以判定算法是否繼續(xù)迭代，并不關(guān)心此時(shí)方程（2.8）的解，只有在達(dá)到收斂條件時(shí)才需要計(jì)算其解?；诖耍覀兛梢詫λ惴?.3進(jìn)行改進(jìn)，采用新的策略直接對（2.8）中的范數(shù)進(jìn)行估計(jì)，而略去計(jì)算中間步驟中的求解過程，只在范數(shù)達(dá)到條件時(shí)才計(jì)算。顯然，這能夠大大節(jié)省計(jì)算成本，加快收斂速度，尤其是在矩陣方程的階數(shù)較大的情形。為此，先介紹下面兩個(gè)引理。

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

下面我們將運(yùn)用算法2.4求解ARE（1.1）的實(shí)例。實(shí)驗(yàn)是在硬件環(huán)境CPU為 Pentium（R） Core Dual E5400，內(nèi)存為 2G ，軟件環(huán)境為 Matlab R2009a 下完成。為簡單起見， 取ARE（1.1）中 E為單位矩陣。由于空間Ek所含信息涵蓋了空間k，

因此，在實(shí)驗(yàn)中用kE代替算法2.4中k。

3.1收斂準(zhǔn)則

在算法2.4中我們將采用如下關(guān)于相對殘差范數(shù)取值作為外迭代收斂準(zhǔn)則

利用（2.11）計(jì)算得到的殘差范數(shù)minRes來判斷內(nèi)迭代是否完成，即是否要計(jì)算相應(yīng)Lyapunov方程的解。為使求解ARE獲得較高精度，這要求相應(yīng)的Lyapunov方程求解具有一個(gè)較高的精度。因此，我們直接取 minRes

3.2預(yù)條件子的取法

由于預(yù)條件子M的主要作用是使M能盡量與à*à接近，又

因此，可取算子M滿足

3.3實(shí)驗(yàn)結(jié)果

我們以一個(gè)鋼軌冷卻系統(tǒng)模型[13]為研究對象，其冷卻過程可以用一個(gè)不穩(wěn)定的二維熱方程進(jìn)行描述，對其進(jìn)行有限元半離散后最終該模型轉(zhuǎn)化為對形如（1.1）的代數(shù)Riccati方程求解。在稍作處理后，取定s=t=6，表3.1給出了階 n=821

表3.1 數(shù)值實(shí)驗(yàn)結(jié)果

從表3.1可以看出，算法2.4對ARE（1.1）的求解是非常有效的。

[1]Laub A.A schur method for solving algebraic Riccati equations[J].IEEE Trans. Automat. Control.,1979,24:913-921.

[2]Byers R.Solving the algebraic Riccati equation with the matrix sign function[J].Linear Algebra Appl.,1987,85:267-279.

[3]H.T.Banks,K.Ito.A numerical algorithm for optimal feedback gains in high dimensional linear quadratic regulator problems[J]. SIAM J.Control.Optim.,1991,29(3):499-515.

[4]J.A.Burns,K.P.Hulsing.Numerical methods for approximating functional gains in LQR boundary control problems[J]. Mathematical and Computer Modelling,2001,33(1):89-100.

[5]D.Kleinman.On an Iterative Technique for Riccati Equation Computations[J].IEEE Transactions on Automat.Control.,1968,13: 114-115.

[6]I.G.Rosen and C. Wang.A multilevel technique for the approximate solution of operator Lyapunov and algebraic Riccati equations [J].SIAM J.Numer.Anal.,1995,32(2):514-541.[7] P. Benner and J. Saak. A Galerkin-Newton-ADI method for solving large-scale algebraic Riccati equations [J]. DFG Priority Programme 1253 "Optimization with Partial Differential Equations", 2010. Preprint SPP1253-090.

[8] Gajic Z, Qureshi M. Lyapunov Matrix equation in system stability and control [M]. Academic Press: San Diego, 1995.

[9] P. Benner, J. R. Li, T. Penzl. Numerical solution of large-scale Lyapunov equations, Riccati equations, and linear-quadratic optimal control problems [J]. Numer. Linear Algebra Appl., 2008, 15(9):755–777.

[10] Lin Y, V. Simoncini. Minimal residual methods for large scale Lyapunov equations [J]. Appl. Numer. Math., 2013, 72:52-71.

[11] Y. Saad. Iterative methods for sparse linear systems [M]. 2nded., Society for Industrial and Applied Mathematics, Philadelphia, PA, 2003.

[12] Kemin Zhou, John C.Doyle, Keith Glover. Robust and optimal control [M]. New Jersey: Prentice Hall, 2007.

[13] V.Simoncini. A new iterative method for solving large-scale Lyapunov matrix equations [J]. SIAM J. Sci. Comput., 2007, 29(3): 1268–1288.

（責(zé)任編校：何俊華）

O241.6； O151.21

A

1673-2219（2015）05-0005-06

2014－12－05

湖南省教育廳資助科研項(xiàng)目（12C0688），湖南省自然科學(xué)基金資助項(xiàng)目（12JJ3077）。

周立平（1978－），男，湖南永州人，副教授，碩士，研究方向?yàn)閿?shù)值代數(shù)和偏微分方程數(shù)值解。