实验八MATLAB状态空间分析
实验八 线性系统的状态空间分析
§8.1 用MATLAB 分析状态空间模型
1、状态空间模型的输入
线性定常系统状态空间模型
x Ax Bu y Cx Du
=+=+ 将各系数矩阵按常规矩阵形式描述。
[][][]11
121120
10
1;;;n n n nn n n A a a a a a a B b b b C c c c D d ====
在MATLAB 里，用函数SS()来建立状态空间模型
(,,,)sys ss A B C D =
例8.1 已知某系统微分方程
22d d 375d d y y y u t t
++= 求该系统的状态空间模型。
解：将上述微分方程写成状态空间形式
0173A ??=??--??，01B ??=????
[]50C =，0D =
调用MATLAB 函数SS()，执行如下程序
% MATLAB Program example 6.1.m
A=[0 1;-7 -3];
B=[0;1];
C=[5 0];
D=0;
sys=ss(A,B,C,D)
运行后得到如下结果
a =
x1 x2
x1 0 1
x2 -7 -3
b =
u1
x1 0
x2 1
c =
x1 x2
y1 5 0
d =
u1
y1 0
Continuous-time model.
2、状态空间模型与传递函数模型转换
状态空间模型用sys 表示，传递函数模型用G 表示。
G=tf(sys)
sys=ss(G)
状态空间表达式向传递函数形式的转换
G=tf(sys)
Or [num,den]=ss2tf(A,B,C,D) 多项式模型参数
[num,den]=ss2tf(A,B,C,D,iu)
[z,p,k]=ss2zp(A,B,C,D,iu) 零、极点模型参数
iu 用于指定变换所需的输入量，iu 默认为单输入情况。
传递函数向状态空间表达式形式的转换
sys=ss(G)
or [A,B,C,D]=tf2ss(num,den)
[A,B,C,D]=zp2ss(z,p,k)
例 8.2
11122211220.560.050.03 1.140.2500.1101001x x u x x u y x y x -??????????=+??????????-????????????????=???????
????? 试用矩阵组[a ，b ，c ，d]表示系统，并求出传递函数。
% MATLAB Program example 6.2.m
a=[-0.56 0.05;-0.25 0];
b=[0.03 1.14;0.11 0];
c=[1 0;0 1];
d=zeros(2,2);
sys=ss(a,b,c,d)
G1=tf(sys)
G2=zpk(sys)
运行后得到如下结果
a =
x1 x2
x1 -0.56 0.05
x2 -0.25 0
b =
u1 u2
x1 0.03 1.14
x2 0.11 0
c =
x1 x2
y1 1 0
y2 0 1
d =
u1 u2
y1 0 0
y2 0 0
Continuous-time model.
Transfer function from input 1 to output...
0.03 s + 0.0055
#1: ---------------------
s^2 + 0.56 s + 0.0125
0.11 s + 0.0541
#2: ---------------------
s^2 + 0.56 s + 0.0125
Transfer function from input 2 to output...
1.14 s
#1: ---------------------
s^2 + 0.56 s + 0.0125
-0.285
#2: ---------------------
s^2 + 0.56 s + 0.0125
Zero/pole/gain from input 1 to output...
0.03 (s+0.1833)
#1: ----------------------
(s+0.5367) (s+0.02329)
0.11 (s+0.4918)
#2: ----------------------
(s+0.5367) (s+0.02329)
Zero/pole/gain from input 2 to output...
1.14 s
#1: ----------------------
(s+0.5367) (s+0.02329)
-0.285
#2: ----------------------
(s+0.5367) (s+0.02329)
例8.3 考虑下面给定的单变量系统传递函数
3243272424()10355024
s s s G s s s s s +++=++++ 由下面的MATLAB 语句直接获得状态空间模型。
>> num=[1 7 24 24];
>> den=[1 10 35 50 24];
>> G=tf(num,den);
>> sys=ss(G)
运行后得到如下结果：
a =
x1 x2 x3 x4
x1 -10 -4.375 -3.125 -1.5
x2 8 0 0 0
x3 0 2 0 0
x4 0 0 1 0
b =
u1
x1 2
x2 0
x3 0
x4 0
c =
x1 x2 x3 x4
y1 0.5 0.4375 0.75 0.75
d =
u1
y1 0
Continuous-time model.
3. 线性系统的非奇异变换与标准型状态空间表达式
syst=ss2ss(sys,T)
sys, syst 分别为变换前、后系统的状态空间模型，T 为非奇异变换阵。
[At,Bt,Ct,Dt]=ss2ss(A,B,C,D,T)
(A,B,C,D)、（At,Bt,Ct,Dt ）分别为变换前、后系统的状态空间模型的系数矩阵。
§8.2 利用MATLAB 求解系统的状态方程
线性定常连续系统状态方程
x Ax Bu =+，0(0)x x =，0t ≥
状态响应
00()()()()d t
x t t x t Bu φφτττ=+-?， 0t ≥ 式中状态转移矩阵()At t e φ=，则有
()0()(0)()d t
At A t x t e x e Bu τττ-=+?， 0t ≥ 1． 用MATLAB 中expm(A)函数计算状态转移矩阵At e
例8.4 022130x x u -????=+????-????，1(0)1x ??=????
，0u = ①求当0.2t =时，状态转移矩阵即0.2At
t e
=；
>> A=[0 -2;1 -3];
>> dt=0.2;
>> phi=expm(A*dt)
得到如下结果
phi =
0.9671 -0.2968
0.1484 0.5219
②计算0.2t =时系统的状态响应
110.2220.2(0)0.96710.29680.6703(0)(0)0.14840.52190.6703At t t x x e x x x ==-????????=?==????????????
????
2． 用step()，impulse() 求阶跃输入，脉冲输入响应
例8.5 连续二阶系统
[]111222120.75240.7268110.72680022.87768.9463x x u x x u x y x ---??????????=+??????????????????????=????
求系统的单位阶跃响应
% MATLAB Program of example 4.5.m
A=[-0.7524 -0.7268;0.7268 0];
B=[1 -1;0 2];
C=[2.8776 8.9463];
D=0;
step(A,B,C,D);
figure(1)
grid on ;
title('单位阶跃响应')
xlabel('时间')
ylabel('振幅')
运行结果
3． 用initial()函数，求系统的零输入响应
[y,t,x]=initial(sys,x 0)
6.5例中，当输入0u =时，状态初值[](0)0.2
0.2x =
A=[-0.7524 -0.7268;0.7268 0];
B=[1 -1;0 2];
C=[2.8776 8.9463];
D=0;
t=[0:0.01:15];u=0;
sys=ss(A,B,C,D);
x0=[0.2 0.2];
[y,t,x]=initial(sys,x0,t)
plot(t,x) 运行结果
§8.3 系统的可控性与可观性分析
1. 线性定常系统的可控性分析
x Ax Bu y Cx Du
=+=+ 可控性矩阵
21[,,,,]n c u B AB A B A B -=，
系统完全可控 rank c u n =。
在MATLAB 中，可用(,)ctrb A B 函数求可控性矩阵c u
ctrb(A,B)
例 8.6 120011101000111x x u ????????=+????????????
， 判断系统的可控性。 ℅MATLAB program of example 6.6.m
A=[1 2 0;1 1 0;0 0 1];
B=[0 1;1 0;1 1];
n=3;
CAM=ctrb(A,B);
rcam=rank(CAM);
if rcam==n
disp('system is controlled') elseif rcam disp('system is not controlled') end 执行结果 system is controlled 例 8.7 2220 ,010,1 2612 x Ax bu A b -- ???? ????=+=-= ???? ???? - ???? 将该系统状态方程转换为可控标准型。 变换矩阵 1 11 1 1 1 ,[0,,0,1] c n P P A P P u P A - - ?? ?? ?? == ?? ?? ?? ℅MATLAB Program of example 6.7.m A=[-2 2 -2;0 -1 0;2 -6 1]; b=[0;1;2]; s=ctrb(A,b); if det(s)~=0 s1=inv(s); end P=[s1(3,:);s1(3,:)*A;s1(3,:)*A*A]; PT=inv(P); A1=P*A*PT%(Ac=PAP^) b1=P*b%(bc=P*b) 运行结果 A1 = 0.0000 1.0000 -0.0000 -0.0000 0 1.0000 -2.0000 -3.0000 -2.0000 b1 = 1.0000 这样可得可控标准型矩阵 110100001,02321c c A A b b ????????====????????---???? 2. 线性定常系统的可观性分析 x Ax Bu y Cx Du =+??=+? 可观性矩阵01n C CA U CA -??????=?????? 系统可观 0rankU n = 在MATLAB 中，可用函数obsv(A,C)确定可观性矩阵。 例8.8 23112211x x u -????=+????-????，2112y x ??=??-?? 确定可观性。 % MATLAB Program of example4.8.m A=[-2 3;3 -2]; B=[1 1;1 1]; C=[2 1;1 -2]; D=0; n=2; ob=obsv(A,C); roam=rank(ob); if roam==n disp('system is observable') elseif roam~=n disp('system is no observable') end 运行结果 system is observable §8.4 用MATLAB 实现极点配置 1. 调用place 函数进行极点配置 k=place(A,B,P) A,B 为系统系数矩阵，P 为配置极点，k 为反馈增益矩阵。 例8.9 给定状态方程 010000010100010001101x x u ???? ???? -????=+????????-???? ，[]1234y x = 将极点配置在1,2,1s j =---±，确定反馈增益矩阵k 。 % MATLAB Program of example 4.9.m A=[0 1 0 0;0 0 -1 0;0 0 0 1;0 0 11 0]; B=[0;1;0;-1]; eig(A)'; P=[-1;-2;-1+sqrt(-1);-1-sqrt(-1)]; k=place(A,B,P) eig(A-B*k)' 运行结果如下： k = -0.4000 -1.0000 -21.4000 -6.0000 ans = -2.0000 -1.0000 - 1.0000i -1.0000 + 1.0000i -1.0000 2. 调用Ackerann 公式计算状态反馈矩阵k A=[0 1 0 0;0 0 -1 0;0 0 0 1;0 0 11 0]; b=[0;1;0;-1]; eig(A)' P=[-1;-2;-1+sqrt(-1);-1-sqrt(-1)]; k = ACKER(A,b,P) eig(A-b*k)' 运行结果 k = -0.4000 -1.0000 -21.4000 -6.0000 §8.5 用MATLAB 设计状态观测器 例6.10 已知系统状态方程 010000010100010001101x x u ??? ? ?? ?? -????=+????????-???? ， []1000y x = （1） 判别可观性； （2） 若系统可观，设计全维状态观测器，使闭环极点为2,3,2,2j j ---+--。 %example4.10 %输入系统状态方程 a=[0 1 0 0;0 0 -1 0;0 0 0 1;0 0 11 0]; b=[0;1;0;-1]; c=[1 0 0 0]; n=4; %计算可观性矩阵 ob=obsv(a,c); roam=rank(ob); %判断可观性 if roam==n disp('system is observable') elseif roam~=n disp('system is no observable') end %求解反馈增益矩阵 a=[0 1 0 0;0 0 -1 0;0 0 0 1;0 0 11 0]; b=[0;1;0;-1]; c=[1 0 0 0]; p1=[-2; -3; -2+sqrt(-1); -2-sqrt(-1)] a1=a'; b1=c'; c1=b'; k=acker(a1,b1,p1) %求解系统矩阵 h=(k)' ahc=a-h*c 运行结果 system is observable h = 9 42 -148 -492 ahc = -9 1 0 0 -42 0 -1 0 148 0 0 1 492 0 11 0 （注：专业文档是经验性极强的领域，无法思考和涵盖全面，素材和资料部分来自网络，供参考。可复制、编制，期待你的好评与关注）