状态转移矩阵计算

1、状态转移矩阵计算 状态转移矩阵计算状态转移矩阵计算(1/1) 3.2 状态转移矩阵计算 q 在状态方程求解中在状态方程求解中,关键是状态转移矩阵关键是状态转移矩阵 (t)的计算。的计算。 对于线性定常连续系统对于线性定常连续系统,该问题又归结为矩阵指数函数该问题又归结为矩阵指数函数 eAt的计算。的计算。 上一节已经介绍了基于拉氏反变换技术的矩阵指数函数上一节已经介绍了基于拉氏反变换技术的矩阵指数函数 eAt的计算方法的计算方法,下面讲述计算矩阵指数函数的下述其他下面讲述计算矩阵指数函数的下述其他 3种常用方法。种常用方法。 级数求和法 约旦规范形法 化eAt为A的有限多项式矩阵函数法 状态转

2、移矩阵计算 级数求和法级数求和法(1/3) 3.2.1 级数求和法 q 由上一节对矩阵指数函数的定义过程中可知由上一节对矩阵指数函数的定义过程中可知: . ! . ! 2 22 k tAtA AtI kk At e 矩阵指数函数eAt的计算可由上述定义式直接计算。 q 由于上述定义式是一个无穷级数,故在用此方法计算eAt时必须 考虑级数收敛性条件和计算收敛速度问题。 类似于标量指数函数eat,对所有有限的常数矩阵A和有限 的时间t来说,矩阵指数函数eAt这个无穷级数表示收敛。 状态转移矩阵计算 级数求和法级数求和法(2/3) q 显然显然,用此方法计算用此方法计算eAt一般不能写成封闭的、简洁

3、的解析形一般不能写成封闭的、简洁的解析形 式式,只能得到数值计算的近似计算结果。只能得到数值计算的近似计算结果。 其计算精度取决于矩阵级数的收敛性与计算时所取的其计算精度取决于矩阵级数的收敛性与计算时所取的 项数的多少。项数的多少。 如果级数收敛较慢如果级数收敛较慢,则需计算的级数项数多则需计算的级数项数多,人工计算是人工计算是 非常麻烦的非常麻烦的,一般只适用于计算机计算。一般只适用于计算机计算。 因此因此,该方法的缺点该方法的缺点: 计算量大计算量大 精度低精度低 非解析方法非解析方法,难以得到计算结果的简洁的解析表达难以得到计算结果的简洁的解析表达 式式 。 状态转移矩阵计算 级数求和法

4、级数求和法(3/3)例例3-4 q 例3-4 用直接计算法求下述矩阵的矩阵指数函数用直接计算法求下述矩阵的矩阵指数函数: .31.32 . 2 3 .1 . ! 232 10 32 10 10 01 . ! . ! 2 2 22 2 2 22 ttt ttt t t k tAtA AtI kk At e 32 10 A q 解 按矩阵指数函数的展开式计算如下: 状态转移矩阵计算 约旦规范形法约旦规范形法 (1/8) 3.2.2 约旦规范形法 q 上节给出了对角线矩阵、块对角矩阵和约旦块三种特殊形上节给出了对角线矩阵、块对角矩阵和约旦块三种特殊形 式矩阵的矩阵指数函数式矩阵的矩阵指数函数。 由于

5、由于任何矩阵都可经线性变换成为对角线矩阵或约旦任何矩阵都可经线性变换成为对角线矩阵或约旦 矩阵矩阵,因此因此 可通过线性变换将一般形式的矩阵变换成对角线可通过线性变换将一般形式的矩阵变换成对角线 矩阵或约旦矩阵矩阵或约旦矩阵, 再利用上述特殊形式矩阵的矩阵指数函数来快速再利用上述特殊形式矩阵的矩阵指数函数来快速 计算矩阵矩阵指数函数。计算矩阵矩阵指数函数。 下面讨论之。下面讨论之。 状态转移矩阵计算 约旦规范形法约旦规范形法(2/8) q 下面首先讨论矩阵指数函数的一条性质下面首先讨论矩阵指数函数的一条性质: 对矩阵对矩阵A,经变换矩阵经变换矩阵P作线性变换后作线性变换后,有有 则相应地有如下

6、矩阵指数函数的变换关系则相应地有如下矩阵指数函数的变换关系 PPPP AttAtAAt eeee 1 1 1 AP AP q 根据上述性质根据上述性质,对矩阵对矩阵A,可通过可通过线性变换方法得到对角线矩线性变换方法得到对角线矩 阵或约旦矩阵阵或约旦矩阵,然后利用该类特殊矩阵的矩阵指数函数然后利用该类特殊矩阵的矩阵指数函数,由矩由矩 阵指数函数的变换关系来求原矩阵阵指数函数的变换关系来求原矩阵A的矩阵指数函数。的矩阵指数函数。 状态转移矩阵计算 约旦规范形法约旦规范形法(4/8)例例3-5 q 例例3-5 试求如下系统矩阵的矩阵指数函数试求如下系统矩阵的矩阵指数函数 5116 6116 110

7、 A q 解解 1. 先求先求A的特征值。的特征值。由特征方程可求得特征值为 1=-1 2=-2 3=-3 2. 求特征值所对应的特征向量。求特征值所对应的特征向量。由前述的方法可求得特征值 1,2和3所对应的特征向量分别为 p1=1 0 1 p2=1 2 4 p3=1 6 9 状态转移矩阵计算 约旦规范形法约旦规范形法例例3-5 q 故将故将A变换成对角线矩阵的变换矩阵变换成对角线矩阵的变换矩阵P及其逆阵及其逆阵P-1为为 12/ 31 343 22/ 53 941 620 111 1 PP 3. 由系统矩阵和矩阵指数函数的变换关系,分别有 t t t tA APPA 3 2 1 e00 0

8、e0 00e e 300 020 001 状态转移矩阵计算 约旦规范形法约旦规范形法例例3-6 2/e27e16-2/5ee9e12-e3 e9e8-e6e6- 2/e3e4-2/5eee3-e3 ee 3232 3232 3232 1 tttttt tttt tttttt tAAt PP ttt tt ttt 32 32 32 e9-e122e- e6-6e e-e3e2- q 例例3-6 试求如下系统矩阵的矩阵指数函数试求如下系统矩阵的矩阵指数函数 032 100 010 A 状态转移矩阵计算 约旦规范形法约旦规范形法(7/8)例例3-6 q 解 1. 先求A的特征值。由特征方程可求得特征

9、值为由特征方程可求得特征值为 1=2 2= 3=-1 2. 由于矩阵由于矩阵A为友矩阵为友矩阵,故将故将A变换成约旦矩阵的变换矩阵变换成约旦矩阵的变换矩阵P和和 其逆阵其逆阵P-1分别为分别为 136 128 121 9 1 214 112 011 1 PP t tt t tA tAPPA e00 ee0 00e e 100 110 002 2 1 3. 由系统矩阵和矩阵指数函数的变换关系,分别有 状态转移矩阵计算 约旦规范形法约旦规范形法(8/8)-例例3-6 tt tt tt tttt tttt tttt tAAt t t t tt tt tt PP e )3-5(4e e )32- (2

10、e e )3-1- (e e )38- (8ee )64- (4e e )3-5(4ee )62(-2e e )32- (2ee )68(e 9 1 ee 2 2 2 22 22 22 1 状态转移矩阵计算 塞尔维斯特内插法塞尔维斯特内插法(1/1) 3.2.3 塞尔维斯特内插法 q 在讨论塞尔维斯特在讨论塞尔维斯特(Sylvester)内插法计算矩阵指数函数内插法计算矩阵指数函数eAt时时, 需要用到关于矩阵特征多项式的凯莱需要用到关于矩阵特征多项式的凯莱-哈密顿哈密顿(Cayley- Hamilton)定理以及最小多项式的概念。定理以及最小多项式的概念。 因此因此,首先给出凯莱首先给出凯莱

11、-哈密顿定理及最小多项式的概念哈密顿定理及最小多项式的概念,再再 讨论塞尔维斯特内插法。讨论塞尔维斯特内插法。 下面依次介绍下面依次介绍: 凯莱-哈密顿定理 最小多项式 塞尔维斯特内插法计算矩阵指数函数 状态转移矩阵计算 凯莱-哈密顿定理(1/4) 1. 凯莱凯莱-哈密顿定理哈密顿定理 q 凯莱凯莱-哈密顿定理是矩阵方程分析和求解中非常重要的定理哈密顿定理是矩阵方程分析和求解中非常重要的定理, 其表述和证明如下。其表述和证明如下。 q 定理定理3-1(凯莱凯莱-哈密顿定理哈密顿定理) 设设n n矩阵矩阵A的特征多项式为的特征多项式为 f( )=| I-A|= n+a1 n-1+an-1 +an

12、 则矩阵则矩阵A必使由上述特征多项式决定的矩阵多项式函数必使由上述特征多项式决定的矩阵多项式函数 f(A)=An+a1An-1+an-1A+anI=0 上述特征多项式亦称为矩阵上述特征多项式亦称为矩阵A的零化特征多项式。的零化特征多项式。 状态转移矩阵计算 凯莱-哈密顿定理(2/4) q 证明证明 因为因为 I=( I-A)-1( I-A)=adj( I-A)/| I-A|( I-A) 故故 | I-A|I=adj( I-A)( I-A) 由伴随矩阵的定义可知由伴随矩阵的定义可知,伴随矩阵伴随矩阵adj( I-A)可表示为如下可表示为如下 多项式矩阵函数多项式矩阵函数: adj( I-A)=

13、n-1I+ n-2B2+ Bn-1+Bn 其中矩阵其中矩阵B2,B3,Bn为为n n维的常数矩阵。维的常数矩阵。 状态转移矩阵计算 凯莱-哈密顿定理(3/4) 因此由前面两式因此由前面两式,有有 ( n+a1 n-1+an-1 +an)I=( n-1I+ n-2B2+ Bn-1+Bn)( I-A) 整理得整理得 ( n+a1 n-1+an-1 +an)I = nI+(B2-A) n-1+(Bn-Bn-1A) -BnA 状态转移矩阵计算 凯莱-哈密顿定理(4/4) 上式中上式中,令等号两边令等号两边 的同幂次项的系数相等的同幂次项的系数相等,则有则有 a1I-B2+A=0 a2I-B3+AB2=

14、0 an-1I-Bn+ABn-1=0 anI+ABn=0 因此因此,将上述各等式从上至下依次右乘以将上述各等式从上至下依次右乘以An-1,A,I,然后然后 将各等式相加将各等式相加,即得即得 An+a1An-1+an-1A+anI=0 故矩阵故矩阵A满足其本身的零化特征多项式。满足其本身的零化特征多项式。 状态转移矩阵计算 最小多项式最小多项式 (1/3) 2. 最小多项式 q 根据根据凯莱凯莱-哈密尔顿定理哈密尔顿定理,任一任一nn维矩阵维矩阵A满足其自身的特满足其自身的特 征方程征方程,即特征多项式为即特征多项式为A的一个零化多项式。的一个零化多项式。 然而特征多项式不一定是然而特征多项式

15、不一定是A的最小阶次的零化多项式。的最小阶次的零化多项式。 将矩阵将矩阵A满足的最小阶次的首一零化多项式称为最小多满足的最小阶次的首一零化多项式称为最小多 项式项式,也就是说也就是说,定义定义nn维矩阵维矩阵A的最小多项式为满足的最小多项式为满足 (A)=Am+ 1Am-1+ m-1A+ mI=0, m n 的阶次最低的首一多项式的阶次最低的首一多项式 ( )= m+ 1 m-1+ m-1 + m 状态转移矩阵计算 最小多项式最小多项式(2/3) q 最小多项式在矩阵多项式的分析与计算中起着重要作用。最小多项式在矩阵多项式的分析与计算中起着重要作用。 定理定理3-2给出了特征多项式与最小多项式

16、的关系。给出了特征多项式与最小多项式的关系。 q 定理3-2 设首一多项式设首一多项式d( )是是 I-A的伴随矩阵的伴随矩阵adj( I-A)的所有的所有 元素的最高公约式元素的最高公约式,则最小多项式为则最小多项式为 )( )( d AI 状态转移矩阵计算 最小多项式最小多项式(3/3) q 证明 由假设知由假设知,矩阵矩阵adj( I-A)的最高公约式为的最高公约式为d( ),故故 adj( I-A)=d( )B( ), 式中式中,B( )的的n2个元素个元素(为为 的函数的函数)的最高公约式为的最高公约式为1。 由于由于 ( I-A)adj( I-A)=| I-A|I 可得可得 d(

17、)( I-A)B( )=| I-A|I 由上式可知由上式可知,特征多项式特征多项式| I-A|可被整除可被整除d( )。 因此设因此设d( )整除整除| I-A|得到的因式记为得到的因式记为 ( ),故有故有 | I-A|=d( ) ( ), 状态转移矩阵计算 最小多项式最小多项式(4/3) q 由于首一多项式由于首一多项式d( )的最高阶次的系数为的最高阶次的系数为1,所以所以 ( )的最高的最高 阶次的系数也应为阶次的系数也应为1。 因此因此,综合上两式综合上两式,可得可得 ( I-A)B( )= ( )I 因而因而 (A)=0 即即 ( )亦为亦为A的零化多项式。的零化多项式。 设设 (

18、 )为为A的最小多项式的最小多项式,因此零化多项式因此零化多项式 ( )可写为可写为 ( )=g( ) ( )+e( ) 其中其中g( )和和e( )分别是多项式分别是多项式 ( )除以除以 ( )的商和余项的商和余项,且且 e( )的阶次低于的阶次低于 ( )。 状态转移矩阵计算 最小多项式最小多项式(5/3) 由于由于 (A)=0和和 (A)=0,所以必然有所以必然有e(A)=0。 考虑到考虑到 ( )为矩阵为矩阵A的最小多项式的最小多项式,所以不存在比所以不存在比 ( ) 阶次还低的阶次还低的A的零化多项式的零化多项式,故故e( )必为零必为零,即有即有 ( )=g( ) ( ) q 又

19、因为又因为 (A)=0,所以所以 ( )可写为可写为 ( )I=( I-A)H( ) 式中式中,H( )为为 ( )的一个因子矩阵的一个因子矩阵,故故 ( )I=g( ) ( )I=g( )( I-A)H( ) 将上式与将上式与( I-A)B( )= ( )I比较比较,有有 B( )=g( )H( ) 状态转移矩阵计算 最小多项式最小多项式(6/3) q 又因为又因为B( )的的n2个元素的最高公约式为个元素的最高公约式为1,因此因此 g( )=1 于是于是 ( )= ( ) 因此因此,由前面证明的由前面证明的| I-A|=d( ) ( )而证明了最小多项式而证明了最小多项式 ( )为为 )(

20、 )( d AI 状态转移矩阵计算 最小多项式最小多项式(7/3) q 根据上述定理根据上述定理3-2,nn维矩阵维矩阵A的最小多项式可按以下步骤的最小多项式可按以下步骤 求出。求出。 1) 根据伴随矩阵根据伴随矩阵adj( I-A),写出作为写出作为 的因式分解多项式的的因式分解多项式的 adj( I-A)的各元素的各元素; 2) 确定作为伴随矩阵确定作为伴随矩阵adj( I-A)各元素的最高公约式各元素的最高公约式d( )。 选取选取d( )的最高阶次系数为的最高阶次系数为1。 如果不存在公约式如果不存在公约式,则则d( )=1; 3) 最小多项式最小多项式 ( )可由可由| I-A|除以

21、除以d( )得到。得到。 状态转移矩阵计算 塞尔维斯特内插法计算矩阵指数函数塞尔维斯特内插法计算矩阵指数函数(1/4) 3. 塞尔维斯特内插法计算矩阵指数函数 q 基于最小多项式基于最小多项式(或特征多项式或特征多项式),塞尔维斯特内插法可以非常塞尔维斯特内插法可以非常 简洁、快速地计算出矩阵指数函数简洁、快速地计算出矩阵指数函数,其计算思想与过程可描其计算思想与过程可描 述如下。述如下。 q 若若 ( )= m+ 1 m-1+ m-1 + m为矩阵为矩阵A的最小多项式的最小多项式,则由则由 (A)=0有有 Am=- 1Am-1- m-1A- mI 即即Am可用有限项可用有限项Am-1,A,I

22、的线性组合来表示。的线性组合来表示。 状态转移矩阵计算 塞尔维斯特内插法计算矩阵指数函数塞尔维斯特内插法计算矩阵指数函数(2/4) q 将上式两边乘以矩阵将上式两边乘以矩阵A,则有则有 即即Am+1可用有限项可用有限项Am-1,A,I的线性组合来表示。的线性组合来表示。 12 11 12 1111 21 12111 . . ().() mm mm m mmmm m mmm AAAA AAIAA AAI 状态转移矩阵计算 塞尔维斯特内插法计算矩阵指数函数塞尔维斯特内插法计算矩阵指数函数(3/4) 其中i(t)(i=0,1,m-1)为待定的关于时间t的函数。 即,矩阵指数函数eAt亦可以用有限项A

23、m-1,A,I的线性函 数组合表示。 q 依次类推依次类推,则可知则可知,Ai(im)可用有限项可用有限项Am-1,A,I的线性组合的线性组合 来表示。来表示。 因此因此,我们有我们有 1 110 22 )(.)()(. ! . ! 2 n n kk At AtAtIt k tAtA AtIe 关键喔 状态转移矩阵计算 塞尔维斯特内插法计算矩阵指数函数塞尔维斯特内插法计算矩阵指数函数(4/4) q 利用上式去计算矩阵指数函数利用上式去计算矩阵指数函数eAt的关键是如何计算待定函的关键是如何计算待定函 数数 i(t)。 下面分下面分 A的特征值互异 A有重特征值 两种情况来讨论如何计算两种情况来

24、讨论如何计算 i(t)以及以及eAt。 状态转移矩阵计算 (1) A的特征值互异 q 设矩阵设矩阵A的的n个互异特征值为个互异特征值为 1, 2, n,则矩阵则矩阵A的最小多项的最小多项 式式 ( )等于特征多项式等于特征多项式f( )=| I-A|= n+a1 n-1+an-1 +an。 因系统的所有特征值因系统的所有特征值 i使特征多项式使特征多项式f( i)=0,故与前面证故与前面证 明过程类似明过程类似,我们亦有我们亦有 A的特征值互异的特征值互异(1/4) nittt niaaa n ini t nin n i n i , 1)(.)()(e , 1. 1 110 1 1 1 i 其

25、中待定函数i(t)(i=0,1,n-1)与矩阵指数函数eAt的表达式中的 i(t)一致。 状态转移矩阵计算 A的特征值互异的特征值互异(2/4) 因此因此,可得如下待定函数可得如下待定函数 i(t)(i=0,1,n-1)的线性方程组的线性方程组: t t t n n nn n n n t t t e . e e )( . )( )( .1 . .1 .1 2 1 1 1 0 1 1 22 1 11 求解上述方程得函数i(t)后,由式(3-49)可计算得矩阵指数 函数eAt。 状态转移矩阵计算 A的特征值互异的特征值互异(3/4)-例3-7 q 例例3-7 试求如下系统矩阵的矩阵指数函数试求如下

26、系统矩阵的矩阵指数函数 6116 100 010 A ttt ttt ttt t t t t t t 32 32 32 3 2 1 2 1 0 ee2e e3e8e5 e2e6e6 2 1 e e e 931 421 111 )( )( )( q 解 由于矩阵A的3个特征值互异,并分别为-1,-2和-3,因此解方 程组(3-52)可得 t t t n n nn n n n t t t e . e e )( . )( )( .1 . .1 .1 2 1 1 1 0 1 1 22 1 11 状态转移矩阵计算 A的特征值互异的特征值互异(4/4) 则系统的状态转移矩阵为则系统的状态转移矩阵为 ttt

27、 ttt ttt tttttt tttttt tttttt At AtAtIt 32 32 32 3232 3232 3232 2 210 e9e8-e e3-e4e- ee2-e e27e32-e5e9e12-e6 e9-e16e5-e6-e12e6- e3e8-e5e2e6-e6 2 1 )()()(e 状态转移矩阵计算 A有重特征值有重特征值(1/4) (2) A有重特征值 q 由于矩阵由于矩阵A与它的约旦矩阵与它的约旦矩阵 具有相同的最小多项式具有相同的最小多项式 ( ),因因 此此由前面的推导过程可知由前面的推导过程可知,约旦矩阵约旦矩阵 也满足也满足A 设A与 的特征值i的代数重数

28、为mi,则由上式很容易证明 i(t)满足 1 011 ( )( ).( ) Atm m t It At A e i i i 1 011 2 121 1 11 e( )( ).( ) e( )2( ).(1)( ) . (1)! e(1)!( ).( ) () ii i tm imi tm imi mm mt immi i ttt tttmt m tmtt mm A A 求解上述方程求解上述方程,则可求得待定函数则可求得待定函数 i(t)。 状态转移矩阵计算 A有重特征值有重特征值(2/4) q 为清楚说明问题为清楚说明问题,设设A和和 有如下有如下6个特征值个特征值: 1, 1, 1, 2,

29、2, 3。 则相应的矩阵指数函数计算式则相应的矩阵指数函数计算式(3-49)中的待定函数中的待定函数 i(t)(i=0, 1,5)的计算式为的计算式为 t t t t t t t t t t t t t t t 3 2 2 1 1 1 e e e ! 1 1 e e ! 1 1 e ! 2 1 1 1 ! 4 ! 5 43210 1 ! 4 ! 5 43210 ! 2 ! 3 ! 5 ! 2 ! 2 ! 4 3100 )( )( )( )( )( )( 2 1 5 3 4 3 3 3 2 33 5 2 4 2 3 2 2 22 4 2 3 2 2 22 5 1 4 1 3 1 2 11 4 1

30、 3 1 2 11 3 1 2 11 5 4 3 2 1 0 A 状态转移矩阵计算 A有重特征值有重特征值(3/4)例例3-8 q 值得指出的是值得指出的是,上述塞尔维斯特内插法不仅对矩阵上述塞尔维斯特内插法不仅对矩阵A的最小多的最小多 项式成立项式成立,而且对所有矩阵而且对所有矩阵A的零化多项式也成立。的零化多项式也成立。 因此因此,在难以求解最小多项式时在难以求解最小多项式时,上述方法中的最小多项上述方法中的最小多项 式可用矩阵式可用矩阵A的特征多项式代替的特征多项式代替,所得结果一致所得结果一致,仅计算量仅计算量 稍大。稍大。 q 例例3-8 试求如下系统矩阵的矩阵指数函数试求如下系统矩阵的矩阵指数函数 010 230 002 A 状态转移矩阵计算 A有重特征值有重