现代控制理论 ——03 状态方程的解
证明
证明
eAt=I+At+2!1A2t2+⋯
=⎣⎢⎢⎢⎢⎡10⋮001⋮0⋯⋯⋱⋯00⋮1⎦⎥⎥⎥⎥⎤+⎣⎢⎢⎢⎢⎡λ10⋮00λ2⋮0⋯⋯⋱⋯00⋮λn⎦⎥⎥⎥⎥⎤t+2!1⎣⎢⎢⎢⎢⎡λ120⋮00λ22⋮0⋯⋯⋱⋯00⋮λn2⎦⎥⎥⎥⎥⎤t2+⋯
=⎣⎢⎢⎢⎢⎡∑n=0+∞n!1λ1ntn0⋮00∑n=0+∞n!1λ2ntn⋮0⋯⋯⋱⋯00⋮∑n=0+∞n!1λnntn⎦⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎡eλ1t0⋮00eλ2t⋮0⋯⋯⋱⋯00⋮eλnt⎦⎥⎥⎥⎥⎤。
- 若A 为m×m 的若尔当块,即A=⎣⎢⎢⎢⎢⎢⎢⎡λ0⋮001λ⋮0001λ⋯⋯⋯⋯⋱λ000⋮1λ⎦⎥⎥⎥⎥⎥⎥⎤m×m,则eAt=eλt⎣⎢⎢⎢⎢⎢⎢⎢⎡10⋮00t1⋮002!t2t⋱⋯⋯⋯⋯⋱10(m−1)!tm−1(m−2)!tm−2⋮t1⎦⎥⎥⎥⎥⎥⎥⎥⎤m×m。
若尔当块
形如⎣⎢⎢⎢⎢⎢⎢⎡λ0⋮001λ⋮0001λ⋯⋯⋯⋯⋱λ000⋮1λ⎦⎥⎥⎥⎥⎥⎥⎤m×m 为m 阶若尔当矩阵,1 阶若尔当矩阵为λ。
- 若A 为一个有多个若尔当块的若尔当矩阵(即若当标准型),即A=⎣⎢⎢⎢⎢⎡A10⋮00A2⋮0⋯⋯⋱⋯00⋮An⎦⎥⎥⎥⎥⎤,则eAt=⎣⎢⎢⎢⎢⎡eA1t0⋮00eA2t⋮0⋯⋯⋱⋯00⋮eAnt⎦⎥⎥⎥⎥⎤。
# 矩阵指数的计算
- 定义计算:eAt=∑n=0+∞n!1An(t)n。该方法适用于计算机运算。
例题
已知A=[0−110],求eAt。
由定义,
eAt=I+At+2!1+⋯=[1001]+[0−tt0]+2!1[−t200−t2]+⋯=[1−2!t2+⋯−t+3!t3−⋯t−3!t3+⋯1−2!t2+⋯]=[cost−sintsintcost].
- 拉氏变换法:利用拉氏变换在频域中求解齐次状态方程的解。
设线性时不变齐次状态方程为 x˙=Ax(t),x(0)=x0,t≥t0。
作拉氏变换有 sX(s)−x(0)=AX(s),即 (sI−A)X(s)=x(0),那么
X(s)=(sI−A)−1x(0)
取拉氏逆变换有 x(0)=L−1[(sI−A)−1x(0)]=L−1[(sI−A)−1]x(0),因此
eAt=L−1[(sI−A)−1](2)
例题
计算矩阵 A=[0−21−3] 的矩阵指数。
由拉氏变换法,(sI−A)=[s2−1s+3],则(sI−A)−1=[(s+1)(s+2)s+3(s+1)(s+2)−2(s+1)(s+2)1(s+1)(s+2)s],
则 eAt=L−1[(s+1)(s+2)s+3(s+1)(s+2)−2(s+1)(s+2)1(s+1)(s+2)s]=[2e−t−2e−t+2e−2te−t−e−2t−e−t+2e−2t]。
- 将矩阵化为对角标准型或若尔当标准型。
若A=⎣⎢⎢⎢⎢⎡λ10⋮00λ2⋮0⋯⋯⋱⋯00⋮λn⎦⎥⎥⎥⎥⎤ 为对角矩阵,则eAt 也为对角矩阵(性质 8),即eAt=⎣⎢⎢⎢⎢⎡eλ1t0⋮00eλ2t⋮0⋯⋯⋱⋯00⋮eλnt⎦⎥⎥⎥⎥⎤。
(1)当矩阵A 的 n 个特征值 λ1,λ2…λn 均两两互异时,则可确定变换阵 P 及其逆矩阵 P−1 ,使得矩阵A 对角化:A=P⎣⎢⎢⎢⎢⎡λ10⋮00λ2⋮0⋯⋯⋱⋯00⋮λn⎦⎥⎥⎥⎥⎤P−1,则有
eAt=P⎣⎢⎢⎢⎢⎡eλ1t0⋮00eλ2t⋮0⋯⋯⋱⋯00⋮eλnt⎦⎥⎥⎥⎥⎤P−1(3)
解题步骤
- 求解系统矩阵A 的特征值 λ1,λ2…λn 。(特征值两两互异)
- 求解特征值对应的特征向量p1,p2…pn,构造变换阵 P 并求解其逆矩阵 P−1 。
- 求解矩阵指数 eAt=P⎣⎢⎢⎢⎢⎡eλ1t0⋮00eλ2t⋮0⋯⋯⋱⋯00⋮eλnt⎦⎥⎥⎥⎥⎤P−1。
例题
试用化为对角标准型法求解矩阵A=[0−21−3] 的矩阵指数 eAt。
求解特征值∣λI−A∣=∣∣∣∣∣λ2−1λ+3∣∣∣∣∣=(λ+1)(λ+2),得到特征值为λ1=−1,λ2=−2。继而求解特征向量p1=[1−1],p2=[1−2]。
故变换矩阵 P=[1−11−2],求逆有 P−1=[2−11−1]。
则矩阵指数为 eAt=P[e−t00e−2t]P−1=[2e−t−e−2t−2e−t+2e−2te−t−e−2t−e−t+2e−2t]。
试用化为对角标准型法求解矩阵A=⎣⎢⎡0−6−61−11−11−165⎦⎥⎤ 的矩阵指数 eAt。
求解特征值∣λI−A∣=∣∣∣∣∣∣∣λ66−1λ+11111−6λ−5∣∣∣∣∣∣∣=(λ+1)(λ+2)(λ+3),得到特征值为λ1=−1,λ2=−2,λ3=−3。继而求解特征向量p1=⎣⎢⎡101⎦⎥⎤,p2=⎣⎢⎡124⎦⎥⎤,p3=⎣⎢⎡169⎦⎥⎤。
故变换矩阵 P=⎣⎢⎡101124169⎦⎥⎤,求逆有 P−1=⎣⎢⎡3−3125−423−23−1⎦⎥⎤。
则矩阵指数为
eAt=P⎣⎢⎡e−t000e−2t000e−3t⎦⎥⎤P−1=⎣⎢⎡101124169⎦⎥⎤⎣⎢⎡e−t000e−2t000e−3t⎦⎥⎤⎣⎢⎡3−3125−423−23−1⎦⎥⎤=⎣⎢⎡3e−t−3e−2t+e−3t−6e−t+6e−3t3e−t−12e−2t+9e−3t25e−t−4e−2t+23e−3t−8e−2t+9e−3t25e−t−16e−2t+227e−3t−2e−t+3e−2t−e−3t6e−2t−6e−3t−2e−t+12e−2t−9e−3t⎦⎥⎤。
(2)当 n×n 矩阵A 有n 重特征根时,存在线性非奇异变换 P 及其逆矩阵 P−1 ,将矩阵 A 转化为若尔当标准型:A=P⎣⎢⎢⎢⎢⎡λ0⋮01λ⋮0⋯⋯⋱⋯001λ⎦⎥⎥⎥⎥⎤n×nP−1,则有
eAt=Peλt⎣⎢⎢⎢⎢⎢⎢⎢⎡10⋮00t1⋮002!t2t⋮00⋯⋯⋱⋯⋯(n−1)!tn−1(n−2)!tn−2⋮t1⎦⎥⎥⎥⎥⎥⎥⎥⎤n×nP−1(4)
拓展到一般情况,矩阵A 同时存在重特征根和单特征根时,以有三重根λ1、两重根λ2 和单根λ3 的矩阵A 为例,若存在变换阵 P 及其逆矩阵 P−1 ,将矩阵 A 转化为若尔当标准型:A=P⎣⎢⎢⎢⎢⎢⎢⎢⎡λ101λ11λ1λ21λ20λ1⎦⎥⎥⎥⎥⎥⎥⎥⎤P−1,则有
eAt=P⎣⎢⎢⎢⎢⎢⎢⎢⎡eλ1t00000teλ1teλ1t000021t2eλ1tteλ1teλ1t000000eλ2t00000teλ2teλ2t000000eλ3t⎦⎥⎥⎥⎥⎥⎥⎥⎤P−1(5)
例题
试求矩阵A=⎣⎢⎡013602−524⎦⎥⎤ 的矩阵指数。
求解特征值∣λI−A∣=∣∣∣∣∣∣∣λ−1−3−6λ−25−2λ−4∣∣∣∣∣∣∣=(λ−1)2(λ−2),得到特征值为λ1=λ2=1,λ3=2。继而求解特征向量和广义特征向量p1=⎣⎢⎡1−73−75⎦⎥⎤,p2=⎣⎢⎡1−4922−4946⎦⎥⎤,p3=⎣⎢⎡2−1−2⎦⎥⎤。
故变换矩阵 P=⎣⎢⎡1−73−751−4922−49462−1−2⎦⎥⎤,求逆有 P−1=⎣⎢⎡27−4−628−115−71⎦⎥⎤。
则矩阵指数为
eAt=P⎣⎢⎡e−t00te−tet000e2t⎦⎥⎤P−1=⎣⎢⎡9et+7tet−8e2t−4et−3tet+4e2t−8et−5tet+8e2t22et+28tet+−22e2t−10et−12tet+11e2t−22et−20tet−22e2t−2et−7tet+2e2tet+3tet−e2t3et+5tet−2e2t⎦⎥⎤。
- 化矩阵指数为矩阵A 的有限项。
该方法将矩阵指数表示为eAt=a0(t)I+a1(t)A+⋯+an−1An−1。
当特征值两两互异时,
⎣⎢⎢⎢⎢⎡a0(t)a1(t)⋮an−1(t)⎦⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎡11⋮1λ1λ2⋮λn⋯⋯⋱⋯λ1n−1λ2n−1⋮λnn−1⎦⎥⎥⎥⎥⎤−1⎣⎢⎢⎢⎢⎡eλ1teλ2t⋮eλnt⎦⎥⎥⎥⎥⎤(6)
当存在重特征值时(以三重根λ1 和二重根λ2,其余根为单根为例),
⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡a0(t)a1(t)a2(t)a3(t)a4(t)a5(t)⋮an−1(t)⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡001011⋮101λ11λ2λ3⋮λn12λ1λ122λ2λ22λ32⋮λn23λ13λ12λ133λ22λ23λ33⋮λn3⋯⋯⋯⋯⋯⋯⋱⋯2!(n−1)(n−2)λ1n−31!(n−1)λ1n−2λ1n−11!(n−1)λ2n−2λ2n−1λ3n−1⋮λnn−1⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤−1⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡2!1t2eλ1t1!1teλ1teλ1t1!1teλ2teλ2teλ3t⋮eλn−3t⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤(7)
证明:Cayley-Hamilton定理
解题步骤
求解系统矩阵A 的特征值 λ1,λ2…λn 。
求解有限项,根据特征值的互异性分情况分析:
- 当特征值两两互异时,直接根据⎣⎢⎢⎢⎢⎡a0(t)a1(t)⋮an−1(t)⎦⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎡11⋮1λ1λ2⋮λn⋯⋯⋱⋯λ1n−1λ2n−1⋮λnn−1⎦⎥⎥⎥⎥⎤−1⎣⎢⎢⎢⎢⎡eλ1teλ2t⋮eλnt⎦⎥⎥⎥⎥⎤ 求解有限项。
- 当特征值存在重根时,对于单根部分列写方程:
eλit=a0(t)+a1(t)λi+⋯+an−1(t)λin−1
而对于k 重根部分在列写方程eλit=a0(t)+a1(t)λi+⋯+ak−1(t)λik−1 外还需要补充方程:
⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧teλit=a1(t)+2a2(t)λi+⋯+(k−1)ak−1(t)λik−2t2eλit=2a2(t)+6a3(t)λi+⋯+(k−1)(k−2)ak−1(t)λik−3⋮tk−1eλit=(k−1)!ak−1(t)
联立n 条方程求解有限项
- 代入求解矩阵指数:
eAt=a0(t)I+a1(t)A+⋯+an−1An−1
例题
试求矩阵A=⎣⎢⎡002103010⎦⎥⎤ 的矩阵指数,利用化为有限项法求解。
求解特征值∣λI−A∣=∣∣∣∣∣∣∣λ0−2−1λ−30−1λ∣∣∣∣∣∣∣=(λ+1)2(λ−2),得到特征值为λ1,2=−1,λ3=2。
对于单根 λ3=2,有 e2t=a0(t)+2a1(t)+4a2(t),
对于二重根λ1,2=−1,有 e−t=a0(t)−a1(t)+a2(t),还需要补充方程:
te−t=a1(t)−2a2(t)
联立三组方程解得:
⎩⎪⎨⎪⎧a0(t)=91(e2t+8e−t+6te−t)a1(t)=91(2e2t−2e−t+3te−t)a3(t)=91(e2t−e−t−3te−t)
eAt=a0(t)I+a1(t)A+⋯+an−1An−1=91⎣⎢⎡e2t+(8+6t)e−t2e2t−(2+6t)e−t4e2t+(6−4t)e−te2t−(2−3t)e−t4e2t+(5−3t)e−t8e2t+(3−8t)e−te2t−(1+3t)e−t2e2t−(2−3t)e−t4e2t+(5−3t)e−t⎦⎥⎤
# 线性时不变系统非齐次状态方程的解
动态系统在控制的作用下的运动称为受控运动。线性时不变系统非齐次状态方程的解即为线性时不变系统的受控运动。考虑系统 x˙(t)=Ax(t)+Bu(t),x(0),t≥0,其动态响应形式为:
x(t)=eA(t−t0)x(t0)+∫t0teA(t−τ)Bu(τ)dτ,t≥0(8)
可理解为由两部分组成:一部分是由初始状态引起的系统自由运动,即零输入响应;另外一部分是由控制输入所产生的受控运动,即零状态响应。
推导过程
对于系统x˙(t)=Ax(t)+Bu(t),x(0),t≥0,左乘e−At 后求导可得:
dtd[e−Atx(t)]=e−At[x˙(t)−Ax(t)]=e−AtBu(t)
两边积分得:
∫0t{dtd[e−Atx(t)]}dτ=∫0te−AtBu(t)dτ
e−Atx(t)−x(0)I=∫0te−AtBu(t)dτ
x(t)=eA(t−t0)x(t0)+∫t0teA(t−τ)Bu(τ)dτ,t≥0
# 线性时不变系统的状态转移矩阵
在线性时不变系统解 x(t)=eA(t−t0)x(t0)+∫t0teA(t−τ)Bu(τ)dτ,t≥0 中,定义状态转移矩阵Φ(t,t0)=eA(t−t0)。
注
- 线性时不变系统的状态转移矩阵可记为:Φ(t,t0)=Φ(t−t0)。
- x(t) 是由初始值引起的零输入解和控制产生的零状态解的叠加。
- 解的结构显示了从x(t0) 到x(t) 的一种变换关系。
线性连续系统的状态转移矩阵
- 定义
对于线性连续系统的状态方程:x˙(t)=A(t)x(t)+B(t)u(t),x(t0)=x0,A(t)∈Rn×n,那么称满足以下矩阵方程的解Φ(t,t0) 为系统的状态转移矩阵。
Φ˙(t,t0)=A(t)Φ(t,t0),Φ(t0,t0)=I,t≥t0(9)
- 状态转移矩阵的性质
- dtdΦ(t,t0)=A(t)Φ(t,t0),Φ(t0,t0)=I
- Φ(t2,t1)Φ(t1,t0)=Φ(t2,t0)
- Φ(mt)=Φ(t+t+⋯+t)=[Φ(t)]m
# 线性时变系统状态方程的解 *
# 线性时变系统齐次状态方程的解
# 线性时变系统的状态转移矩阵
# 线性时变系统非齐次状态方程的解
# 线性连续系统的时间离散化
线性连续系统的时间离散化问题本质上就是在一定的采样方式和保持方式下,由系统的连续时间状态空间描述来得到对应的离散时间状态空间描述,并建立两者的系数矩阵间的关系式。
# 近似离散化
考虑以下线性时变系统:x˙(t)=A(t)x(t)+B(t)u(t),当采样周期T 较小且精度要求不高时,可将其离散化为:
x˙(kT)≈T1[x((k+1)T)−x(kT)](10)
令t=kT,有
T1[x((k+1)T)−x(kT)]=A(kT)x(kT)+B(kT)u(kT)
x[(k+1)T]=[I+TA(kT)]x(kT)+TB(kT)u(kT)=G(kT)x(kT)+H(kT)u(kT)
其中,G(kT)=I+TA(kT),H(kT)=TB(kT)。
注:一般而言,当采样周期为系统最小时间系数的101 左右,近似度已经足够。
例题
系统的状态方程为x˙(t)=A(t)x(t)+B(t)u(t),其中A(t)=[005(1−e−5t)5(e−5t−1)],B(t)=[505e−5t5(1−e−5t)]。试求采样周期为T=0.2s 时的离散状态方程。
直接代入公式有:
G(kT)=I+TA(kT)=[101−e−ke−k]
H(kT)=TB(kT)=[10e−k1−e−k]
那么,离散状态方程为:x[(k+1)T]=G(kT)x(kT)+H(kT)u(kT)
将状态方程x˙=[0−21−3]x+[01]u 近似离散化,T=0.2s。
由题:G=I+TA=[1001]+0.2[0−21−3]=[1−0.40.20.4],H=0.2[01]=[00.2]。
故离散状态方程为:
$x[0.2(k+1)]=\begin{bmatrix}1&0.2\\-0.4&0.4\end{bmatrix}x(0.2k)+\begin{bmatrix}0\\0.2\end{bmatrix}u(0.2k)$
# 线性时不变系统状态方程的离散化
在线性时不变系统中,x˙(t)=A(x)+B(u),其时间离散化状态方程为:
x[(k+1)T]=Gx(kT)+Hu(kT)(11)
其中G=eAT,H=(∫0TeATdt)B。假设条件为:(1) 等采样周期T;(2)u(t)≡u(kT),kT≤t≤(k+1)T。
推导证明
对于线性时不变系统 x˙(t)=A(x)+B(u),其状态方程的解为:
x(t)=eA(t−t0)x(t0)+∫eA(t−τ)Bu(τ)dτ(12)
假设:(1) 等采样周期T;(2)x(k)=[x(t)]t=kT,u(k)=[u(t)]_
那么令 t=(k+1)T,t0=kT,有:
x[(k+1)T]=eATx(kT)+∫kT(k+1)TeA[(k+1)T−τ]Bu(τ)dτ=eATx(kT)+∫kT(k+1)TeA[(k+1)T−τ]Bdτ⋅u(kT)
令 t=(k+1)T−τ,dτ=−dt,有:
x[(k+1)T]=eATx(kT)+∫0τeA(t)Bdt⋅u(kT)=eATx(kT)+∫0τeA(t)dt⋅Bu(kT)
令G=eAT,H=(∫0TeATdt)B,有线性时不变系统的离散状态方程为:
x[(k+1)T]=Gx(kT)+Hu(kT)
解题步骤
- 求解矩阵指数,方法见矩阵指数的计算。
- 求解系数矩阵:G=eAT,H=(∫0TeATdt)B。
- 列写时间离散化状态方程:x[(k+1)T]=Gx(kT)+Hu(kT)
例题
将状态方程x˙=[001−2]x+[01]u 离散化,T=0.1s。
利用拉氏变换法求解矩阵指数函数。取拉氏变换有:
[sI−A]−1=[s0−1s+2]−1=[s10s(s+2)1s+21]
取拉氏逆变换得到矩阵指数函数:
eAt=L−1[sI−A]−1=[100.5(1−e−2T)e−2T]
进而求解系数矩阵:
G=eAT=[100.5(1−e−2T)e−2T]=[100.0910.819]
H=(∫0TeATdt)B=[∫0T[100.5(1−e−2T)e−2T]dt][01]=[T00.5T+0.25e−2T−0.25−0.5e−2T+0.5][01]=[0.0050.091]
故时间离散化状态方程为:
x[0.1(k+1)]=[100.0910.819]x(0.1k)+[0.0050.091]u(0.1k)
# 线性离散系统状态方程的解
离散系统的差分方程形状态方程有两种解法:递推法和 z 变换法。其中递推法在时变系统和时不变系统中都适用,而 z 变换法只适用于时不变系统。
# 递推法
- 在线性时变系统中,x(k+1)=G(k)x(k)+H(k)u(k),有:
⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧x(1)=G(0)x(0)+H(0)u(0)x(2)=G(1)x(1)+H(1)u(1)x(3)=G(2)x(2)+H(2)u(2)⋮
给定初始条件x(0) 和输入序列u(0),u(1),⋯ 后即可求解x(k)。
- 在线性时不变系统中,x(k+1)=Gx(k)+Hu(k),其中G,H 均为常数矩阵,因此:
x(k)=Gkx(0)+i=0∑k−1Gk−1−iHu(i)(13)
上式称为线性时不变离散系统的状态转移方程,其中Φ(k)=Gk 称为线性时不变离散系统的状态转移矩阵。
状态转移矩阵的性质:
- Φ(k+1)=GΦk,Φ(0)=I
- Φ(k2−k0)=Φ(k2−k1)Φ(k1−k0)
- Φ−1(k)=Φ(−k)
# z 变换法
考虑时不变离散系统:x(k+1)=Gx(k)+Hu(k),取 z 变换有:
zx(z)−zx(0)=Gx(z)+Hu(z)
z(z)=(zI−G)−1zx(0)+(zI−G)−1Hu(z)(14)
取 z 逆变换有:
x(k)=z−1[(zI−G)−1z]x(0)+z−1[(zI−G)−1Hu(z)](15)
对比公式(13)和公式(15),由解的唯一性可知,
z−1[(zI−G)−1z]=Gk(16)
z−1[(zI−G)−1Hu(z)]=i=0∑k−1Gk−1−iHu(i)(17)
例题
考虑离散系统:x(k+1)=Gx(k)+Hu(k),其中G=[0−0.161−1],H=[11],初始条件为x(0)=[1−1],试求当u(k)=1 时状态方程的解。
用 z 变换法求解,先计算(zI−G)−1,有
(zI−G)−1=[z0.16−1z+1]−1=(z+0.2)(z+0.8)1[z+1−0.161z]=[34×z+0.21−31×z+0.81−30.8×z+0.21+30.8×z+0.8135×z+0.21−35×z+0.81−31×z+0.21+34×z+0.81]
由于u(k)=1,则u(z)=z−1z,故zx(0)+Hu(z)=[z−z]+[z−1zz−1z]=[z−1z2z−1−z2+2z]。
那么代入公式(15)有:
x(z)=(zI−G)−1[zx(0)+Hu(z)]=[−617×z+0.2z+922×z+0.8z+1825×z−1z63.4×z+0.2z−917.6×z+0.8z+187×z−1z]
求 z 逆变换有:
x(k)=[−617(−0.2)k+922(−0.8)k+182563.4(−0.2)k−917.6(−0.2)k+187]