离散耦合Riccati矩阵方程的无反演牛顿迭代法
Newton Iteration Method without Inversion for Discrete Coupled Riccati Matrix Equations
摘要: 本文研究了离散耦合Riccati矩阵方程的求解问题。我们利用矩阵反演等式得到离散耦合Riccati矩阵方程的等价形式,并通过牛顿求根方法构造无反演的迭代算法求解离散耦合Riccati矩阵方程。此外我们证明了无反演牛顿迭代算法的收敛性。最后,利用数值例子来验证提出的理论结果。
Abstract: In this paper, the discrete coupled Riccati matrix equation is considered. We obtain the equivalent form of the discrete coupled Riccati matrix equation by using the matrix inversion equation. And an iterative algorithm without inversion is constructed to solve the discrete coupled Riccati matrix equation by Newton’s method of finding roots. In addition, we prove the convergence of Newtonian iterative algorithm without inversion. Finally, a numerical example demonstrates the theoretical results.
文章引用:庄柔柔. 离散耦合Riccati矩阵方程的无反演牛顿迭代法[J]. 应用数学进展, 2025, 14(4): 342-354. https://doi.org/10.12677/aam.2025.144167

1. 引言

在控制理论的发展中,求解Riccati矩阵方程是非常重要的课题。从六十年代到如今,大量的研究成果设计出的最优控制器[1] [2]都与Riccati矩阵方程存在唯一正定解等价,因此,得到耦合Riccati矩阵方程的解是非常重要的。

Riccati矩阵方程是非线性矩阵方程,其常用的求解方法是迭代求解。文献[3]研究了Markov系统所得到的耦合Riccati矩阵方程的迭代算法,并提出将该方程变成Lyapunov矩阵方程来进行迭代求解。文献[4]使用牛顿法来获得Riccati矩阵方程的唯一正定解,给出在任意初始条件下所求的解都收敛于矩阵方程的唯一解。文献[5]根据矩阵特征值的性质,得到解的存在性并给出一种固定点迭代算法求解矩阵方程。文献[6]利用最新估计信息,给出带有权重因子的并行迭代算法,该方法对算法的收敛性能有显著的提高。考虑到离散耦合Riccati矩阵方程的特殊性,文献[7]利用一种保持结构的倍速算法对矩阵方程进行求解,讨论了保持结构倍速算法的收敛性,并通过引入最新估计信息对算法进一步改进。文献[8]通过处理一类非线性矩阵方程中的矩阵求逆运算,提出了一种近似替换的迭代算法。

基于以上研究,本文在文献[4]算法基础上提出求解离散耦合Riccati矩阵方程的无反演牛顿迭代算法,并对算法进行收敛性分析,最后通过数值例子来验证所提算法是有效的。

2. 主要问题和预备知识

本文主要研究如下的离散耦合Riccati矩阵方程:

P i = A i T G i A i A i T G i B i ( R i + B i T G i B i ) 1 B i T G i A i + Q i , (1)

其中 iS S={ 1,2,,N } 是有限集合, G i = ε i ( P ) ε i ( P )= j=1 N π ij P j P=( P 1 , P 2 ,, P N ) π ij 是转移概率。

A i R n×n R i R m×m B i R n×m 以及 Q i R n×n 是对称半正定矩阵, P i R n×n 是方程(1)的对称半正定解。

我们的目标是通过牛顿方法构造无逆的迭代算法,在介绍无反演牛顿迭代算法之前,有必要陈述一些需要使用到的结论。

引理2.1 [9] A C 是两个具有适当维数的非奇异矩阵, B D 是另外两个具有适当维数的矩阵,那么有

( A+BCD ) 1 = A 1 A 1 B ( D A 1 B+ C 1 ) 1 D A 1 . (2)

那么利用引理2.1中的矩阵反演等式,可以将方程(1)等价地表示为:

P i = A i T ( G i 1 + S i ) 1 A i + Q i ,iS, (3)

其中 S i = B i R i 1 B i T

引理2.2 [10]假设 T M R n×n 是具有相同阶的对称正定矩阵,那么

TM M 1 T 1 .

引理2.3 [11]如果 E F R n×n 是具有相同阶的对称正定矩阵,且 F>0 ,那么

EFE+ F 1 2E.

回顾一下牛顿法求解 X 1 A=0 的根的方法:

X r+1 = X r [ ( 1+τ )IτA X r ],r=0,1,2, (4)

其中 τ>0

根据等式(4),为了避免方程(3)的矩阵反演,本文提出以下迭代格式:

{ G i ( m+1 )= j=1 i1 π ij P j ( m+1 )+ j=i N π ij P j (m) X i ( m+1 )= X i ( m )[ ( γ+1 )Iγ G i ( m+1 ) X i ( m ) ],γ>0 Y i ( m+1 )= Y i ( m )[ ( η+1 )Iη( X i ( m+1 )+ S i ) Y i ( m ) ],η( 0,1 ] P i ( m+1 )= A i T Y i ( m+1 ) A i + Q i ,iS. (5)

算法1:无反演牛顿迭代算法

步骤1 设置容许误差 ε ,输入 A i B i Q i R i ,选择初始迭代矩阵 P i ( 0 )= Q i X i ( 0 )= G i 1 ( 0 ) Y i ( 0 )= ( G i 1 ( 0 )+ S i ) 1 ,其中 G i ( 0 )= ε i ( P i ( 0 ) ) S i = B i R i 1 B i T iS

步骤2 m:=0 ,计算

log 10 δ( m )= log 10 i=1 N Δ i ( m ) F 2

其中

Δ i ( m )= P i ( m ) A i T G i ( m ) A i + A i T G i ( m ) B i ( R i + B i T G i ( m ) B i ) 1 B i T G i ( m ) A i Q i G i ( m )= j=1 N π ij P j ( m )

步骤3 如果 log 10 δ( m )<ε ,则停止计算,否则

G i ( m+1 )= j=1 i1 π ij P j ( m+1 )+ j=i N π ij P j ( m ),

X i ( m+1 )= X i ( m )[ ( γ+1 )Iγ G i ( m+1 ) X i ( m ) ],γ>0

Y i ( m+1 )= Y i ( m )[ ( η+1 )Iη( X i ( m+1 )+ S i ) Y i ( m ) ],η( 0,1 ],

P i ( m+1 )= A i T Y i ( m+1 ) A i + Q i ,i=1,2,,N,

log 10 δ( m+1 )= log 10 i=1 N Δ i ( m+1 ) F 2 .

步骤4 m:=m+1 ,返回步骤3。

注释2.1 在算法1中,每次迭代中的多次矩阵乘法运算构成了主要的计算量。假设矩阵维度为 n×n ,根据矩阵运算复杂度理论,一次矩阵乘法的计算复杂度为 O( n 3 ) 。设每次迭代执行矩阵乘法的次数为 m ,那么每次迭代的计算复杂度为 O( m n 3 )

注释2.2 在算法1中,参数 γ η 对算法收敛性产生了影响,理论上,可以通过构建误差方程对迭代格式(5)进行敏感性分析,通过分析误差方程组的解随 γ η 变化的情况,推导 γ η 在不同取值下对矩阵序列收敛性的影响。

3. 无反演牛顿迭代算法的收敛性

定理3.1 在算法1中,设 P i = P i* iS 是方程(1)的唯一正定解,记 P * =( P 1* , P 2* ,, P N* ) G i* = ε i ( P * ) ,初始条件如下:

{ P i ( 0 )= Q i X i ( 0 )= ( j=1 N π ij P j ( 0 ) ) 1 Y i ( 0 )= ( X i ( 0 )+ S i ) 1 . (6)

如果对所有 m0 ,当 γ η0 时,有

η( 1γ )0 ( 1+γ ) X i ( m )γ X i ( m ) G i ( m+1 ) X i ( m ) G i 1 ( m+1 ) i=1,2,,N (7)

那么由算法1生成的序列 P i ( m ) X i ( m ) Y i ( m ) 是有界的,即对所有 m0 有:

{ P i ( m ) P i* ,iS X i ( m ) G i* 1 ,iS Y i ( m ) ( G i* 1 + S i ) 1 ,iS (8)

证明 我们使用数学归纳法来证明这个定理,首先证明 m=0 的情况。由 P i = P i* iS ,有

P i* = A i T ( G i* 1 + S i ) 1 A i + Q i iS

可得到

P i ( 0 ) P i* iS (9)

G i ( 0 )= j=1 N π ij P j ( 0 ) iS (10)

因此,对 iS

G i ( 0 ) j=1 N π ij ( A j T ( G j* 1 + S j ) 1 A j + Q j ) = j=1 N π ij P j* = G i* .

则有

X i ( 0 )= G i 1 ( 0 ) G i* 1 iS (11)

另外,对 iS

Y i ( 0 )= ( X i ( 0 )+ S i ) 1 ( G i* 1 + S i ) 1 (12)

因此,当 m=0 iS (8)成立。假设当 m=ω 时(8)成立,即

{ P i ( ω ) P i* ,iS X i ( ω ) G i* 1 ,iS Y i ( ω ) ( G i* 1 + S i ) 1 ,iS. (13)

根据(13)得到

G 1 ( ω+1 )= j=1 N π ij P j ( ω ) j=1 N π ij P j* = G 1* .

由条件(7)有

X 1 ( ω+1 )= X 1 ( ω )( ( γ+1 )Iγ G 1 ( ω+1 ) X 1 ( ω ) ) = X 1 ( ω )( γ+1 )γ X 1 ( ω ) G 1 ( ω+1 ) X 1 ( ω ) G 1 1 ( ω+1 ) G 1* 1 . (14)

根据引理2.3和(14)有

Y 1 ( ω+1 )= Y 1 ( ω )[ ( η+1 )Iη( X 1 ( ω+1 )+ S 1 ) Y 1 ( ω ) ] = Y 1 ( ω )[ ( 1η )I+η( 2I( X 1 ( ω+1 )+ S 1 ) Y 1 ( ω ) ) ] =( 1η ) Y 1 ( ω )+η[ 2 Y 1 ( ω ) Y 1 ( ω )( X 1 ( ω+1 )+ S 1 ) Y 1 ( ω ) ] ( 1η ) Y 1 ( ω )+η ( X 1 ( ω+1 )+ S 1 ) 1 ( 1η ) ( G i* 1 + S 1 ) 1 +η ( G i* 1 + S 1 ) 1 = ( G 1* 1 + S 1 ) 1 . (15)

则有

P 1 ( ω+1 )= A 1 T Y 1 ( ω+1 ) A 1 + Q 1 A 1 T ( G 1* 1 + S 1 ) 1 + Q 1 = P 1* . (16)

故此,已证明

{ P 1 ( ω+1 ) P 1* X 1 ( ω+1 ) G 1* 1 Y 1 ( ω+1 ) ( G 1* 1 + S 1 ) 1 .

假设对 il ,有

{ P i ( ω+1 ) P i* X i ( ω+1 ) G i* 1 Y i ( ω+1 ) ( G i* 1 + S i ) 1 (17)

成立。根据(13)和(17)得到

G l+1 ( ω+1 )= j=1 l π l+1,j P j ( ω+1 )+ j=l+1 N π l+1,j P j ( ω ) j=1 l π l+1,j P j* + j=l+1 N π l+1,j P j* = G ( l+1 )* .

可得到

X l+1 ( ω+1 )= X l+1 ( ω )[ ( γ+1 )Iγ G l+1 ( ω+1 ) X l+1 ( ω ) ] =( γ+1 ) X l+1 ( ω )γ X l+1 ( ω ) G l+1 ( ω+1 ) X l+1 ( ω ) G l+1 1 ( ω+1 ) G ( l+1 )* 1 (18)

Y l+1 ( ω+1 )= Y l+1 ( ω )[ ( η+1 )Iη( X l+1 ( ω+1 )+ S l+1 ) Y l+1 ( ω ) ] = Y l+1 ( ω )[ ( 1η )I+η( 2I( X l+1 ( ω+1 )+ S l+1 ) Y l+1 ( ω ) ) ] =( 1η ) Y l+1 ( ω )+η[ 2 Y l+1 ( ω ) Y l+1 ( ω )( X l+1 ( ω+1 )+ S l+1 ) Y l+1 ( ω ) ] ( 1η ) Y l+1 ( ω )+η ( X l+1 ( ω+1 )+ S l+1 ) 1 ( 1η ) ( G ( l+1 )* 1 + S l+1 ) 1 +η ( G ( l+1 )* 1 + S l+1 ) 1 = ( G ( l+1 )* + S l+1 ) 1 . (19)

P l+1 ( ω+1 )= A l+1 T Y l+1 ( ω+1 ) A l+1 + Q l+1 A l+1 T ( G ( l+1 )* 1 + S l+1 ) 1 A l+1 + Q l+1 = P ( l+1 )* . (20)

根据(17)~(20)和(17)对 il 成立,因此(17)对所有 iS 成立。此外,根据(9) (11) (12)和归纳假设(13)可知(8)对任意 ω>0 成立,则定理得证。

定理3.1揭示了算法1产生的矩阵序列 P i ( m ) 有上界 P i* X i ( m ) 有下界 G i* 1 Y i ( m ) 有上界 ( G i* 1 + S i ) 1

定理3.2 在算法1中,设 P i = P i* iS 是方程(1)的唯一正定解,记 P * =( P 1* , P 2* ,, P N* ) G i* = ε i ( P * ) ,如果对所有 m0 记,当 γ η0 时满足条件(7),那么由算法1生成的序列 P i ( m ) X i ( m ) Y i ( m ) ,对所有 m0 满足以下结论:

{ P i ( m+1 ) P i ( m ),iS X i ( m+1 ) X i ( m ),iS Y i ( m+1 ) Y i ( m ),iS. (21)

证明 我们使用数学归纳法来证明这个定理。根据迭代格式(5)有

G 1 ( 1 )= j=1 N π 1j P j ( 0 ).

X 1 ( 1 )= X 1 ( 0 )[ ( γ+1 )Iγ G 1 ( 1 ) X 1 ( 0 ) ] =( γ+1 ) X 1 ( 0 )γ X 1 ( 0 ) G 1 ( 0 ) X 1 ( 0 ) = X 1 ( 0 ),

Y 1 ( 1 )= Y 1 ( 0 )[ ( η+1 )Iη( X 1 ( 1 )+ S 1 ) Y 1 ( 0 ) ] =( η+1 ) Y 1 ( 0 )η Y 1 ( 0 )( X 1 ( 0 )+ S 1 ) ( X 1 ( 0 )+ S 1 ) 1 = Y 1 ( 0 )>0,

P 1 ( 1 )= A 1 T Y 1 ( 1 ) A 1 + Q 1 Q 1 = P 1 ( 0 ).

这表明当 i=1 m=0 时(21)成立。我们要证当 m=0 时,对所有 iS

{ P i ( 1 ) P i ( 0 ) X i ( 1 ) X i ( 0 ) Y i ( 1 ) Y i ( 0 ) (22)

成立。当 i=1 时已证,假设 il 时(22)成立。则可得到

G l+1 ( 1 )= j=1 l π l+1,j P j ( 1 )+ j=l+1 N π l+1,j P j ( 0 ) j=1 l π l+1,j P j ( 0 ) + j=l+1 N π l+1,j P j ( 0 ) = G l+1 ( 0 ). (23)

X l+1 ( 1 )= X l+1 ( 0 )[ ( γ+1 )Iγ G l+1 ( 1 ) X l+1 ( 0 ) ] X l+1 ( 0 )[ ( γ+1 )Iγ G l+1 ( 0 ) X l+1 ( 0 ) ] =( γ+1 ) X l+1 ( 0 )γ X l+1 ( 0 ) G l+1 ( 0 ) X l+1 ( 0 ) = X l+1 ( 0 ). (24)

Y l+1 ( 1 )= Y l+1 ( 0 )[ ( η+1 )Iη( X l+1 ( 1 )+ S l+1 ) Y l+1 ( 0 ) ] Y l+1 ( 0 )[ ( η+1 )Iη( X l+1 ( 0 )+ S l+1 ) Y l+1 ( 0 ) ] =( η+1 ) Y l+1 ( 0 )η Y l+1 ( 0 )( X l+1 ( 0 )+ S l+1 ) Y l+1 ( 0 ) =( η+1 ) Y l+1 ( 0 )η Y l+1 ( 0 ) = Y l+1 ( 0 ). (25)

P l+1 ( 1 )= A l+1 T Y l+1 ( 1 ) A l+1 + Q l+1 A l+1 T Y l+1 ( 0 ) A l+1 + Q l+1 = P l+1 ( 0 ). (26)

因此,(22)对所有 iS 均成立。假设(21)对 m=ω 成立,即对所有 iS

{ P i ( ω+1 ) P i ( ω ) X i ( ω+1 ) X i ( ω ) Y i ( ω+1 ) Y i ( ω ). (27)

根据(27)可得

G 1 ( ω+2 )= j=1 N π 1j P j ( ω+1 ) j=1 N π 1j P j ( ω )= G 1 ( ω+1 ).

由条件(7)可知 X 1 ( ω+1 ) G 1 1 ( ω+1 ) ,则

X 1 ( ω+2 ) X 1 ( ω+1 )= X 1 ( ω+1 )[ ( γ+1 )Iγ G 1 ( ω+2 ) X 1 ( ω+1 ) ] X 1 ( ω+1 ) =γ X 1 ( ω+1 )γ X 1 ( ω+1 ) G 1 ( ω+2 ) X 1 ( ω+1 ) γ X 1 ( ω+1 )γ X 1 ( ω+1 ) G 1 ( ω+1 ) X 1 ( ω+1 ) =γ X 1 ( ω+1 )( X 1 1 ( ω+1 ) G 1 ( ω+1 ) ) X 1 ( ω+1 )0,

因此

X 1 ( ω+2 ) X 1 ( ω+1 ). (28)

根据引理3.3可知 2 X 1 ( ω+1 ) X 1 ( ω+1 ) G 1 ( ω+2 ) X 1 ( ω+1 ) G 1 1 ( ω+2 ) ,且由(8)可知 Y 1 ( ω+1 ) ( G 1 1 ( ω+1 )+ S 1 ) 1 以及 η( 1γ )0 ,则有

Y 1 ( ω+2 )= Y 1 ( ω+1 )[ ( η+1 )Iη( X 1 ( ω+2 )+ S 1 ) Y 1 ( ω+1 ) ] = Y 1 ( ω+1 )[ ( η+1 )Iη( ( ( γ+1 ) X 1 ( ω+1 )γ X 1 ( ω+1 ) G 1 ( ω+2 ) X 1 ( ω+1 ) )+ S 1 ) Y 1 ( ω+1 ) ] = Y 1 ( ω+1 ) [ ( η+1 )Iη( ( 1γ ) X 1 ( ω+1 ) + γ( 2 X 1 ( ω+1 ) X 1 ( ω+1 ) G 1 ( ω+2 ) X 1 ( ω+1 ) )+ S 1 ) Y 1 ( ω+1 ) ] Y 1 ( ω+1 )[ ( η+1 )Iη( ( 1γ ) X 1 ( ω+1 )+γ G 1 1 ( ω+2 )+ S 1 ) Y 1 ( ω+1 ) ] Y 1 ( ω+1 )[ ( η+1 )Iη( ( 1γ ) G 1 1 ( ω+1 )+γ G 1 1 ( ω+1 )+ S 1 ) Y 1 ( ω+1 ) ] = Y 1 ( ω+1 )[ ( η+1 )Iη( G 1 1 ( ω+1 )+ S 1 ) Y 1 ( ω+1 ) ] Y 1 ( ω+1 )[ ( η+1 )Iη( G 1 1 ( ω+1 )+ S 1 ) ( G 1 1 ( ω+1 )+ S 1 ) 1 ] =( η+1 ) Y 1 ( ω+1 )η Y 1 ( ω+1 )= Y 1 ( ω+1 ). (29)

P 1 ( ω+2 )= A 1 T Y 1 ( ω+2 ) A 1 + Q 1 A 1 T Y 1 ( ω+1 ) A 1 + Q 1 = P 1 ( ω+1 ). (30)

根据(28)~(30)可知

{ P i ( ω+2 ) P i ( ω+1 ) X i ( ω+2 ) X i ( ω+1 ) Y i ( ω+2 ) Y i ( ω+1 ) (31)

i=1 成立。假设(31)对所有 il 均成立,根据(27)可得到

G l+1 ( ω+2 )= j=1 l π l+1,j P j ( ω+2 )+ j=l+1 N π l+1,j P j ( ω+1 ) j=1 l π l+1,j P j ( ω+1 )+ j=l+1 N π l+1,j P j ( ω ) = G l+1 ( ω+1 ).

由条件(7)可知 X l+1 ( ω+1 ) G l+1 1 ( ω+1 ) ,则

X l+1 ( ω+2 ) X l+1 ( ω+1 ) = X l+1 ( ω+1 )[ ( γ+1 )Iγ G l+1 ( ω+2 ) X l+1 ( ω+1 ) ] X l+1 ( ω+1 ) =γ X l+1 ( ω+1 )γ X l+1 ( ω+1 ) G l+1 ( ω+2 ) X l+1 ( ω+1 ) γ X l+1 ( ω+1 )γ X l+1 ( ω+1 ) G l+1 ( ω+1 ) X l+1 ( ω+1 ) =γ X l+1 ( ω+1 )( X l+1 1 ( ω+1 ) G l+1 ( ω+1 ) ) X l+1 ( ω+1 )0,

因此

X l+1 ( ω+2 ) X l+1 ( ω+1 ). (32)

根据引理2.3可知 2 X l+1 ( ω+1 ) X l+1 ( ω+1 ) G l+1 ( ω+2 ) X l+1 ( ω+1 ) G l+1 1 ( ω+2 ) ,且由(8)可知 Y l+1 ( ω+1 ) ( G l+1 1 ( ω+1 )+ S l+1 ) 1 以及 η( 1γ )0 ,则有

Y l+1 ( ω+2 )= Y l+1 ( ω+1 )[ ( η+1 )Iη( X l+1 ( ω+2 )+ S l+1 ) Y l+1 ( ω+1 ) ] = Y l+1 ( ω+1 ) [ ( η+1 )Iη( ( ( γ+1 ) X l+1 ( ω+1 ) γ X l+1 ( ω+1 ) G l+1 ( ω+2 ) X l+1 ( ω+1 ) )+ S l+1 ) Y l+1 ( ω+1 ) ] = Y l+1 ( ω+1 ) [ ( η+1 )Iη( ( 1γ ) X l+1 ( ω+1 )+γ( 2 X l+1 ( ω+1 ) X l+1 ( ω+1 ) G l+1 ( ω+2 ) X l+1 ( ω+1 ) )+ S l+1 ) Y l+1 ( ω+1 ) ] Y l+1 ( ω+1 )[ ( η+1 )Iη( ( 1γ ) X l+1 ( ω+1 )+γ G l+1 1 ( ω+2 )+ S l+1 ) Y l+1 ( ω+1 ) ] Y l+1 ( ω+1 )[ ( η+1 )Iη( ( 1γ ) G l+1 1 ( ω+1 )+γ G l+1 1 ( ω+1 )+ S l+1 ) Y l+1 ( ω+1 ) ] = Y l+1 ( ω+1 )[ ( η+1 )Iη( G l+1 1 ( ω+1 )+ S l+1 ) Y l+1 ( ω+1 ) ] Y l+1 ( ω+1 )[ ( η+1 )Iη( G l+1 1 ( ω+1 )+ S l+1 ) ( G l+1 1 ( ω+1 )+ S l+1 ) 1 ] =( η+1 ) Y l+1 ( ω+1 )η Y l+1 ( ω+1 )= Y l+1 ( ω+1 ). (33)

P l+1 ( ω+2 )= A l+1 T Y l+1 ( ω+2 ) A l+1 + Q l+1 A l+1 T Y l+1 ( ω+1 ) A l+1 + Q l+1 = P l+1 ( ω+1 ). (34)

根据(32)~(34)以及归纳假设可知,(31)对所有 iS 成立。此外,结合(22) (27) (31)以及归纳假设,通过数学归纳可知(21)对任意 m0 成立,则定理得证。

定理3.2揭示了算法1产生的 P i ( m ) 是单调递增序列, X i ( m ) 是单调递减序列, Y i ( m ) 是单调递增序列。

定理3.3 P i = P i* iS 为方程(3)的唯一正定解,记 P * =( P 1* , P 2* ,, P N* ) 。此外,对所有 iS ,记 G i* = ε i ( P * ) 。如果选择 γ η 满足条件(7),那么算法1生成的矩阵序列 P i ( m ) iS 收敛于方程(3)的唯一解,即

lim m P i ( m )= P i* ,iS.

证明 根据定理3.1和3.2,由算法1生成的矩阵序列 P i ( m ) Y i ( m ) iS 是单调递增且有上界,矩

阵序列 X i ( m ) 是单调递减且有下界,因此,矩阵序列 P i ( m ) Y i ( m ) X i ( m ) 是收敛的,则有 lim m P i ( m )= P i * lim m X i ( m )= X i * lim m Y i ( m )= Y i * iS 。记 lim m G i ( m )= G i *

通过对(5)中第一个表达式的两边都取极限,可得到

G i * = lim m G i ( m )= j=1 N π ij P j * iS (35)

通过对(5)中第二个表达式的两边都取极限,可得到

X i * = X i * [ ( γ+1 )Iγ G i * X i * ] (36)

由(36)得到 X i * = ( G i * ) 1

另外,对(5)中第三个表达式的两边都取极限,可得到

Y i * = Y i * [ ( η+1 )Iη( X i * + S i ) Y i * ] (37)

由(37)有

Y i * = ( X i * + S i ) 1 (38)

最后,对(5)中第四个表达式的两边都取极限,得到

P i * = A i T Y i * A i + Q i (39)

X i * = ( G i * ) 1 和(38)代入(39)得到

P i * = A i T ( ( G i * ) 1 + S i ) 1 A i + Q i iS (40)

根据(35)中的形式,可以发现(40)相较于与方程(3)具有相同的形式。由于方程(3)有一个唯一正定解,因此 P i = P i * iS 是方程(3)的正定解,则定理得证。

注释3.1 本文算法的收敛性依赖条件(7),对于算法中条件(7)的普适性分析,可基于矩阵性质推导研究矩阵 A B C D 的特征值、特征向量与条件(7)的关系,探讨其如何改变条件(7)成立的可能性。

4. 数值例子

我们给出1个数值实例来辅助验证本文的结论。定义迭代误差为 log 10 δ( m ) δ( m ) 定义如下:

δ( m )= i=1 N Δ i ( m ) F 2

其中

Δ i ( m )= P i ( m ) A i T G i ( m ) A i + A i T G i ( m ) B i ( R i + B i T G i ( m ) B i ) 1 B i T G i ( m ) A i Q i

G i ( m )= j=1 N π ij P j ( m )

考虑以下矩阵方程组:

A 1 =( 0.0325 0.1778 0.1101 0.0083 0.0567 0.1276 0.1319 0.0100 0.1358 0.1358 0.1959 0.0139 0.0990 0.0630 0.0252 0.1225 0.1659 0.0382 0.1026 0.0400 0.0215 0.0757 0.1617 0.0804 0.2155 )

A 2 =( 0.2221 0.1400 0.0089 0.2177 0.0762 0.0757 0.1935 0.0973 0.0387 0.1502 0.2046 0.1509 0.0928 0.1328 0.1390 0.0120 0.0673 0.2262 0.1191 0.1461 0.1095 0.0223 0.1211 0.1158 0.1074 )

A 3 =( 0.1506 0.1796 0.2282 0.1311 0.0937 0.1664 0.9690 0.6038 0.8946 0.8855 0.1404 0.0893 0.0707 0.0865 0.1863 0.1651 0.1361 0.0509 0.1648 0.1770 0.4160 0.7876 0.9033 0.1543 0.3971 )

B 1 =( 0.1225 0.0205 0.4479 0.2369 0.1088 0.3594 0.5310 0.2926 0.3102 0.5904 0.4187 0.6748 0.6263 0.5094 0.7620 )

B 2 =( 0.0033 0.5025 0.9186 0.9195 0.4898 0.0944 0.3192 0.0119 0.7228 0.1705 0.3982 0.7014 0.5524 0.7818 0.4850 )

B 3 =( 0.6814 0.3000 0.2967 0.4914 0.6068 0.6617 0.6286 0.4978 0.1705 0.5130 0.2321 0.0994 0.8585 0.0534 0.8344 )

转移概率矩阵为:

π=( 0.15 0.60 0.25 0.45 0.25 0.30 0.60 0.20 0.20 )

D 1 =[ 0.1067  0.9619  0.0046 ]

D 2 =[ 0.7749  0.8173  0.8687 ]

D 3 =[ 0.0844  0.3998  0.2599 ]

I=diag( [ 1.0,1.0,1.0 ] ) R 1 =I+ D 1 T D 1 R 2 =I+ D 2 T D 2 R 3 =I+ D 3 T D 3

S 1 = B 1 R 1 1 B 1 T S 2 = B 2 R 2 1 B 2 T S 3 = B 3 R 3 1 B 3 T

Q 1 =diag( [ 0.5629, 0.3843, 0.3128, 0.7344, 0.5809 ] )

Q 2 =diag( [ 0.5629, 0.3843, 0.3128, 0.7344, 0.5809 ] )

Q 3 =diag( [ 0.5629, 0.3843, 0.3128, 0.7344, 0.5809 ] )

图1表明,无反演牛顿迭代算法在选取不同 γ η 下能收敛到方程(1)的唯一解。因此,能有效说明本文提出算法是收敛的。与文献[8]相比,本文研究的问题更有针对性,提出的算法避免了直接矩阵求逆运算,是对传统牛顿迭代法的改进。

Figure 1. The convergence curve of algorithm 1 under different values with γ , η

1. 算法1在不同 γ η 取值下的收敛曲线

5. 结论

本文利用矩阵反演等式得到离散耦合Riccati矩阵方程的等价形式,并利用牛顿求根方法构造无反演的牛顿算法来求解方程。此外,验证了无反演牛顿迭代算法具有收敛性。值得注意的是,初始矩阵的选取会影响到牛顿求根方法,如何选择初始矩阵将会是我们继续的研究对象。

参考文献

[1] Shang, Y., Liu, K., Cui, N., Wang, N., Li, K. and Zhang, C. (2020) A Compact Resonant Switched-Capacitor Heater for Lithium-Ion Battery Self-Heating at Low Temperatures. IEEE Transactions on Power Electronics, 35, 7134-7144.
https://doi.org/10.1109/tpel.2019.2954703
[2] Liu, K., Hu, X., Wei, Z., Li, Y. and Jiang, Y. (2019) Modified Gaussian Process Regression Models for Cyclic Capacity Prediction of Lithium-Ion Batteries. IEEE Transactions on Transportation Electrification, 5, 1225-1236.
https://doi.org/10.1109/tte.2019.2944802
[3] Costa, O.L.V. (1995) Discrete-Time Coupled Riccati Equations for Systems with Markov Switching Parameters. Journal of Mathematical Analysis and Applications, 194, 197-216.
https://doi.org/10.1006/jmaa.1995.1294
[4] Erfanifar, R., Sayevand, K. and Hajarian, M. (2022) Convergence Analysis of Newton Method without Inversion for Solving Discrete Algebraic Riccati Equations. Journal of the Franklin Institute, 359, 7540-7561.
https://doi.org/10.1016/j.jfranklin.2022.07.048
[5] Liu, J. and Zhang, J. (2011) The Existence Uniqueness and the Fixed Iterative Algorithm of the Solution for the Discrete Coupled Algebraic Riccati Equation. International Journal of Control, 84, 1430-1441.
https://doi.org/10.1080/00207179.2011.604794
[6] Ai-Guo Wu, and Guang-Ren Duan, (2015) New Iterative Algorithms for Solving Coupled Markovian Jump Lyapunov Equations. IEEE Transactions on Automatic Control, 60, 289-294.
https://doi.org/10.1109/tac.2014.2326273
[7] Zhang, J. and Li, S. (2020) The Structure-Preserving Doubling Numerical Algorithm of the Continuous Coupled Algebraic Riccati Equation. International Journal of Control, Automation and Systems, 18, 1641-1650.
https://doi.org/10.1007/s12555-019-0368-y
[8] Li, W. and Li, Z. (2010) A Family of Iterative Methods for Computing the Approximate Inverse of a Square Matrix and Inner Inverse of a Non-Square Matrix. Applied Mathematics and Computation, 215, 3433-3442.
https://doi.org/10.1016/j.amc.2009.10.038
[9] Kailath, T. (1980) Linear Systems (Vol. 156). Prentice-Hall.
[10] Horn, R.A. and Johnson, C.R. (2012) Matrix Analysis. 2nd Edition, Cambridge University Press.
https://doi.org/10.1017/cbo9781139020411
[11] Zhan, X. (1996) Computing the Extremal Positive Definite Solutions of a Matrix Equation. SIAM Journal on Scientific Computing, 17, 1167-1174.
https://doi.org/10.1137/s1064827594277041