现代控制理论控-实验报告.doc

实验一 MATLAB控制工具箱的应用及线性系统的运动分析

实验日期 2013年5月21日 成绩：______

批阅教师签字

1. 实验目的
⑴ 学习掌握MATLAB控制工具箱中的基本命令
⑵ 掌握线性系统的运动分析方法
2. 实验环境
MATLAB6.5
3. 实验内容、源程序代码、实验数据及结果分析
⑴自选控制对象模型，应用以下命令，并写出结果
1) step, damp, pzmap, rlocus, rlocfind, bode, margin, nyquist
2) tf2ss, ss2tf, tf2zp, zp2ss
3) ss2ss, Jordan, canon, eig
实现如下:
1）① step 单位阶跃响应
当输入为单位阶跃信号时，系统的输出为单位阶跃响应。
其调用格式为: [y, x, t] = step(num, den, t) 或 step(num, den)
例1: 求系统传递函数为 W=的单位阶跃响应
num = [1]; den = [1 0.5 1];
sys = tf(num,den);
step(sys)
grid; %绘制表格
xlabel('t'); ylabel('y');
title('单位阶跃响应');

② damp 返回系统特征值、衰减与频（借用例1）
num = [1]; den = [1 0.5 1];
sys = tf(num,den);
damp(sys)
其输出结果为:
Eigenvalue Damping Freq.(rad/s)
-2.50e-001 + 9.68e-001i 2.50e-001 1.00e+000
-2.50e-001 - 9.68e-001i 2.50e-001 1.00e+000
③ pzmap 绘制系统零极点图
其调用格式为: [p, z] = pzmap(sys) 绘图并返回极点p, 零点z
或者 pzmap(sys) 直接绘图
④ rlocus 绘制系统根轨迹
例2: 给定系统G= , 绘制零极点图
程序如下:
num = [1 2 1 3]; den = [1 0.5 2 1];
sys = tf(num,den);
[p,z] = pzmap(sys)
rlocus(sys);
title('系统闭环根轨迹');

>> p = 0 + 1.4142i >>z = -2.1746
0 - 1.4142i 0.0873 + 1.1713i
-0.5000 0.0873 - 1.1713i
⑤ rlocusfind 计算与根轨迹上极点相应的根轨迹增益
例3: 已知一个单位负反馈系统开环传递函数为: W= 试绘制闭环系统的根轨迹；在根轨迹上任选一点，并计算该点的增益k及其所有极点的位置
num = 1; den = conv([1 0],conv([0.5 1],[4 1]));
sys = tf(num,den);
rlocus(sys);
[k,poles] = rlocfind(sys)
title('系统闭环根轨迹')

>> Select a point in the graphics window (用光标选)
selected_point = -0.7583 - 5.1056i
k = 278.3081
poles = -6.0141 1.8821 + 4.4267i 1.8821 - 4.4267i
⑥ bode 绘制给定系统的伯德图 （以例3为例）
num = 1;
den = conv([1 0],conv([0.5 1],[4 1]));
sys = tf(num,den);
bode(sys);
grid on;
title('Bode Diagram');
xlabel('Frequency(rad/sec)');
ylabel('Gain dB');

⑦ margin 求解稳定裕度 (以例1为例 )
num = [1]; den = [1 0.5 1];
sys = tf(num,den);
[Gm,Pm,Wcg,Wcp] = margin(num,den);
>> Gm = Inf Pm = 41.4091 Wcg = Inf Wcp = 1.3229
返回的变量对 (Gm, Wcg) 为系统的增益裕度的值和与之对应的频率
而 (Pm, Wcp) 为系统的相位裕度的值和与之对应的频率
⑧ nyquist 绘制奈氏曲线（以例1为例）
num = [1]; den = [1 0.5 1];
nyquist(num,den);

2)
① ss2tf 将状态空间模型转换成传递函数模型
A=[0 1 0;0 0 1;-4 -3 -2];B=[1;3;-6];
C=[1 0 0]; D=0;
[num,den]=ss2tf(A,B,C,D,1)
>> num = 0 1.0000 5.0000 3.0000
den = 1.0000 2.0000 3.0000 4.0000
② tf2ss 将传递函数转换成状态空间模型
num=[1 5 3];den=[1 2 3 4];
[A,B,C,D]=tf2ss(num,den)
>> A = -2 -3 -4 B= 1 C= 1 5 3 D = 0
1 0 0 0
0 1 0 0
③ tf2zp 将传递函数转换成零极点增益模型
num=[1 2 1 3];den=[1 0.5 2 1];
[z,p,k] = tf2zp(num,den)
>> z = -2.1746 p = 0 + 1.4142i k = 1
0.0873 + 1.1713i 0 - 1.4142i
0.0873 - 1.1713i -0.5000
④ zp2ss 将零极点增益模型转换成状态空间模型
z = [-2.1746 0.0873+1.1713i 0.0873-1.1713i];
p = [0+1.4142i 0-1.4142i -0.5000];
k = 1;
[A,B,C,D] = zp2ss(z,p,k)
>> A = -0.5000 0 0 B = 1
1.6746 0 -1.4142 0
0 1.4142 0 0
C = 1.6746 -0.1746 -0.4387 D = 1
3)
① ss2ss 求状态空间模型的坐标
A=[0 1 0;0 0 1;-5 -3 -2];B=[0;0;1];
C=[1.5 1 0.5];D=0;
T=eye(3,3);
sys=ss(A,B,C,D);
SYS=ss2ss(sys,T)
>> a = x1 x2 x3 b = u1
x1 0 1 0 x1 0
x2 0 0 1 x2 0
x3 -5 -3 -2 x3 1
c = x1 x2 x3 d = u1
y1 1.5 1 0.5 y1 0
Continuous-time model.
② Jordan 约旦标准型 eig 求矩阵的特征值
A=sym(gallery(3)); %构造3x3的随即矩阵
E=eig(A) %求A的特征值
J=jordan(A) %求A的约旦标准型
>> E = [ 1] J = [ 1, 0, 0]
[ 2] [ 0, 2, 0]
[ 3] [ 0, 0, 3]
从运算结果可知A的特征值即为A的约旦标准型对角数值
③ canon 将状态方程化成规范性
Matlab 提供的canon函数只有两种模式即model和companion,model为约旦标准型
A=[0 1 0;0 0 1;-5 -3 -2];B=[0;0;1];
C=[1.5 1 0.5];D=0;
sys=ss(A,B,C,D);
[Gt,P]=canon(sys,'model')
>>a = x1 x2 x3 b = u1
x1 -0.07813 1.645 0 x1 0.5106
x2 -1.645 -0.07813 0 x2 -0.6646
x3 0 0 -1.844 x3 0.6859
c = x1 x2 x3 d = u1
y1 -0.06476 -0.4517 0.3395 y1 0
Continuous-time model.
P = -1.9420 -0.1119 0.5106
-1.6442 -2.1171 -0.6646
1.8602 0.1072 0.6859
⑵ 掌握线性系统的运动分析方法
1) 已知A = ，求e。（用三种方法求解）
① 状态转移矩阵的指数矩阵计算法
a=[0 1;-2 -3]; %输入状态矩阵A
syms t; %定义变量
eat1=expm(a*t) %求e
>> eat1 = [ -exp(-2*t)+2*exp(-t), exp(-t)-exp(-2*t)]
[ -2*exp(-t)+2*exp(-2*t), 2*exp(-2*t)-exp(-t)]
② 状态转移矩阵的拉氏变换计算法
a=[0 1;-2 -3]; %输入状态矩阵A
syms s t;
G=inv(s*eye(size(a))-a) %求预解矩阵
eat2=ilaplace(G) %求拉氏反变换
>> G = [ (s+3)/(s^2+3*s+2), 1/(s^2+3*s+2)]
[ -2/(s^2+3*s+2), s/(s^2+3*s+2)]
eat2 = [ -exp(-2*t)+2*exp(-t), exp(-t)-exp(-2*t)]
[ -2*exp(-t)+2*exp(-2*t), 2*exp(-2*t)-exp(-t)]
③ 非奇异变换法-计算特征值及特征向量矩阵
a=[0 1;-2 -3]; %输入状态矩阵A
syms t; %定义变量
[P,D]=eig(a); %求A的特征值和特征向量(P为特征向量构成的矩阵;D为特征值构成的矩阵)
Q=inv(P); %求变换阵的逆阵
eat3=P*expm(D*t)*Q %求 e
>> eat3 = [ -exp(-2*t)+2*exp(-t), exp(-t)-exp(-2*t)]
[ -2*exp(-t)+2*exp(-2*t), 2*exp(-2*t)-exp(-t)]
2) 利用MATLAB求解书上例2.8题，并画出状态响应和输出响应曲线，求解时域性能指标。（加图标题、坐标轴标注及图标）
a=[-1 0;0 -2];b=[1;1];
c=[1.5 0.5];d=0;
[num,den]=ss2tf(a,b,c,d);
sys=tf(num,den);
figure(1);
[p.z]=pzmap(sys);
rlocus(sys);
grid;
title('系统根轨迹图');
x0=[2;3]; %初始状态
syms s t;
G0=inv(s*eye(size(a))-a)
x1=ilaplace(G0)*x0 %求零输入响应x1
G1=inv(s*eye(size(a))-a)*b
x2=ilaplace(G1/s) %求零状态响应x2
x=x1+x2 %系统的状态响应x
y=c*x %系统的输出响应y
for I=1:61; %计算在各时间点状态xt和yt
tt=0.1*(I-1);
xt(:,I)=subs(x(:),'t',tt);
yt(I)=subs(y,'t',tt);
end;
figure(2);
plot(0:60,[xt;yt]); %绘响应曲线
grid;
xlabel('t(s)');
ylabel('ampitude');
title('系统的状态响应(相交叉-下)和输出响应(叠加-上)')

>> G0 = [ 1/(s+1), 0] x1 = [ 2*exp(-t)]
[ 0, 1/(s+2)] [ 3*exp(-2*t)]
G1 = [ 1/(s+1)] x2 = [ 1-exp(-t)]
[ 1/(s+2)] [ 1/2-1/2*exp(-2*t)]
x = [ exp(-t)+1] y = 3/2*exp(-t)+7/4+5/4*exp(-2*t)
[ 5/2*exp(-2*t)+1/2]
3) 利用MATLAB求解书上例2.12题，并画出状态响应和输出响应曲线。（加图标题、坐标轴标注及图标）
g=[0 1;-0.16 -1];h=[1;1];
G=ss(g);
x0=[1;-1];
syms s t;
G0=inv(s*eye(size(g))-g)
x1=ilaplace(G0)*x0
G1=inv(s*eye(size(g))-g)*h
x2=ilaplace(G1/s)
x=x1+x2
y=c*x
for I=1:61;
tt=0.1*(I-1);
xt(:,I)=subs(x(:),'t',tt);
yt(I)=subs(y,'t',tt);
end;
plot(0:60,[xt;yt]);
xlabel('t(s)');
ylabel('ampitude');
title('系统的状态响应和输出响应')

>> G0 = [ 25*(s+1)/(25*s^2+25*s+4), 25/(25*s^2+25*s+4)]
[ -4/(25*s^2+25*s+4),25*s/(25*s^2+25*s+4)]
x1 = [ 4/3*exp(-4/5*t)-1/3*exp(-1/5*t)]
[ -16/15*exp(-4/5*t)+1/15*exp(-1/5*t)]
G1 = [ 25*(s+1)/(25*s^2+25*s+4)+25/(25*s^2+25*s+4)]
[ -4/(25*s^2+25*s+4)+25*s/(25*s^2+25*s+4)]
x2 = [ 5/2*exp(-4/5*t)-15*exp(-1/5*t)+25/2]
[ -2*exp(-4/5*t)+3*exp(-1/5*t)-1]
x = [ 23/6*exp(-4/5*t)-46/3*exp(-1/5*t)+25/2]
[ -46/15*exp(-4/5*t)+46/15*exp(-1/5*t)-1]
y = 253/60*exp(-4/5*t)-322/15*exp(-1/5*t)+73/4
4）
① P36 1.4(2) 已知系统的状态空间表达式为: =+
Y= 用MATLAB编程求系统的传递函数矩阵
A=[2 1 4;0 2 0;0 0 1];B=[1 0;3 4;2 1];C=[3 5 1];D=[0 0];
[num,den]=ss2tf(A,B,C,D,1)
>> num = 0 20.0000 -29.0000 -13.0000
den = 1 -5 8 -4
② P36 1.5(3) 已知系统传递函数阵: G(s)=
求系统的状态空间模型
num=[1 4 2 2;3 1 1 0];den=[1 2 3 2];
[A,B,C,D]=tf2ss(num,den)
>> A = -2 -3 -2 B = 1 C = 2 -1 0 D = 1
1 0 0 0 -5 -8 -6 3
0 1 0 0
③ P56 2.3(3) 已知系统方成如下:=x + u y=x 求输入和初值为以下值时的状态响应和输出响应:u(t)=1(t), x(0)=
a=[0 1;-6 -5];b=[1;0];c=[1.5 0.5];d=0;
G=ss(a,b,c,d);
x0=[1;1];
syms s t;
G0=inv(s*eye(size(a))-a);
x1=ilaplace(G0)*x0
G1=inv(s*eye(size(a))-a)*b;
x2=ilaplace(G1/s)
x=x1+x2
y=c*x
for I=1:61;
tt=0.1*(I-1);
xt(:,I)=subs(x(:),'t',tt);
yt(I)=subs(y,'t',tt);
end;
plot(0:60,[xt;yt]);
xlabel('t(s)');
ylabel('ampitude');
title('系统的状态响应和输出响应');

运行结果如下:
>> x1 = [ -3*exp(-3*t)+4*exp(-2*t)]
[ -8*exp(-2*t)+9*exp(-3*t)]
x2 = [ 2/3*exp(-3*t)-3/2*exp(-2*t)+5/6]
[ -1+3*exp(-2*t)-2*exp(-3*t)]
x = [ -7/3*exp(-3*t)+5/2*exp(-2*t)+5/6]
[ -5*exp(-2*t)+7*exp(-3*t)-1]
y = 5/4*exp(-2*t)+3/4

实验二 系统的能空性、能观测性、稳定性分析及实现

实验日期 2013年5月21日 成绩：______

批阅教师签字

1. 实验目的
加深理解能观测性、能控性、稳定性、最小实现等概念。掌握如何使用MATLAB进行以下分析和实现
⑴ 系统的能观测性、能控性分析；
⑵ 系统的稳定性分析；
⑶ 系统的最小实现。
2. 实验环境
MATLAB6.5
3. 实验内容、源程序代码、实验数据及结果分析
⑴ 能控性、能观测性及系统实现
(a) 了解以下命令的功能，自选对象模型，进行运算，并写出结果。
gram, ctrb, obsv, lyap, ctrbf, obsvf, mineral
实现如下:
① gram
%利用gram矩阵判断系统的能控性和能观测性
a=[0 1 0;0 0 1;-6 -11 -6];b=[0 1;1 0;0 1];
c=[1 0 1;0 1 0];d=0;
G=ss(a,b,c,d);
wc=gram(G,'c'),nc=det(wc) %求能控性Gram矩阵及其行列式的值
if nc~=0,disp('System is controllable!'),%判断系统的能控性
else disp('System is uncontrollable!'),
end
wo=gram(G,'o'),no=det(wo) %求能控性Gram矩阵及其行列式的值
if no~=0,disp('System is observable!'), %判断系统的能控性
else disp('System is uncontrollable!'),
end
运行结果:
>> wc = 1.7000 -0.5000 -0.7000
-0.5000 0.7000 -0.5000
-0.7000 -0.5000 1.7000
nc = 0.4800
System is controllable!
wo = 1.2167 0.0500 0.0833
0.0500 1.2250 0.0500
0.0833 0.0500 0.0917
no = 0.1253
System is observable!
② ctrb obsv
%利用ctrb和obsv矩阵判断系统的能控性和能观测性
a=[0 1 0;0 0 1;-6 -11 -6];b=[0 1;1 0;0 1];
c=[1 0 1;0 1 0];d=0;
n=length(a) %求系统阶次
qc=ctrb(a,b),nc=rank(qc) %求能控性判别矩阵及其秩
if n==nc,disp('System is controllable!'), %判断系统的能控性
else disp('System is uncontrollable!'),
end
qo=obsv(a,c),no=rank(qo) %求系统能控性判别矩阵及其秩
if n==no,disp('System is observable!'), %判断系统的能控性
else disp('System is uncontrollable!'),
end
运行结果:
>> n = 3
qc = 0 1 1 0 0 1
1 0 0 1 -11 -12
0 1 -11 -12 60 61
nc = 3
System is controllable!
qo = 1 0 1
0 1 0
-6 -10 -6
0 0 1
36 60 26
-6 -11 -6
no = 3
System is observable!
③ lyap
X=lyap(A,C) 求解满足李雅普诺夫方程的对称矩阵X
其中A,C为给定矩阵，且C为对称矩阵。
A=[0 1;-0.5 -1];
C=eye(2,2);
P=lyap(A,C)
>>P = 2.2 2.2
2.2 2.2
④ ctrbf obsvf
%能控性结构分解和能观测性结构分解
a=[0 0 -1;1 0 -3;0 1 -3];
b=[1;1;0];c=[0 1 -2];d=0;
[Ac,Bc,Cc,Tc,Kc]=ctrbf(a,b,c)
[Ao,Bo,Co,To,Ko]=obsvf(a,b,c)
运行结果:
>> Ac = -1.0000 0.0000 -0.0000 Bc = 0
-2.1213 -2.5000 0.8660 0
-1.2247 -2.5981 0.5000 -1.4142
Cc = 1.7321 1.2247 -0.7071
Tc = -0.5774 0.5774 -0.5774 Kc = 1 1 0
0.4082 -0.4082 -0.8165
-0.7071 -0.7071 0
Ao = -1.0000 -1.3416 -3.8341 Bo = 1.2247
0.0000 -0.4000 -0.7348 -0.5477
0 0.4899 -1.6000 -0.4472
Co = 0 0.0000 -2.2361
To = 0.4082 0.8165 0.4082 Ko = 1 1 0
-0.9129 0.3651 0.1826
0 -0.4472 0.8944
⑤ mineral 最小实现
num=[2 2];den=[1 6 11 6];
G=tf(num,den);Gs=ss(G);
Gm=minreal(Gs);
Am=Gm.a
Bm=Gm.b
Cm=Gm.c
Dm=Gm.d
运行结果:
>> 1 state removed.
Am = 2.0121 -5.7764 Bm = 0.1212
3.4812 -7.0121 -0.4848
Cm = 0.5000 0.1250 Dm = 0
(b) 已知连续系统的传递函数模型G(s)=，当a分别取-1、0、1时，判别系统的能控性与能观测性
程序如下:
① a = -1
num=[1 -1];den=[1 10 27 18];
[A,B,C,D]=tf2ss(num,den)
n=length(A)
Qc=ctrb(A,B)
nc=rank(Qc)
if n==nc,disp('System is controllable!'),
else disp('System is uncontrollable!'),
end
Qo=obsv(A,C)
no=rank(Qo)
if n==no,disp('System is observable!'),
else disp('System is unobservable!'),
end
运行结果:
>> A = -10 -27 -18 B = 1 C = 0 1 -1 D = 0
1 0 0 0
0 1 0 0 n = 3
Qc = 1 -10 73
0 1 -10
0 0 1 nc = 3
System is controllable!
Qo = 0 1 -1
1 -1 0
-11 -27 -18 no = 3
System is observable!
② a = 0
>> A = -10 -27 -18 B = 1 C = 0 1 0 D = 0
1 0 0 0
0 1 0 0 n = 3
Qc = 1 -10 73
0 1 -10
0 0 1 nc = 3
System is controllable!
Qo = 0 1 0
1 0 0
-10 -27 -18 no = 3
System is observable!
③ a = 1
>> A = -10 -27 -18 B = 1 C = 0 1 -1 D = 0
1 0 0 0
0 1 0 0 n = 3
Qc = 1 -10 73
0 1 -10
0 0 1 nc = 3
System is controllable!
Qo = 0 1 1
1 1 0
-9 -27 -18 no = 2
System is unobservable!
(c) 已知系统阵为A=, B=, C=, 判别系统的能控性与能观测性
程序如下:
A=[6.666 -10.6667 -0.3333;1 0 1;0 1 2];
B=[0;1;1];C=[1 0 2];
n=length(A)
Qc=ctrb(A,B)
nc=rank(Qc)
if n==nc,disp('System is controllable!'),
else disp('System is uncontrollable!'),
end
Qo=obsv(A,C)
no=rank(Qo)
if n==no,disp('System is observabel!'),
else disp('System is unobservable!'),
end
运行结果:
>> n = 3
Qc = 0 -11.0000 -84.9926
1.0000 1.0000 -8.0000
1.0000 3.0000 7.0000
nc = 3
System is controllable!
Qo = 1.0000 0 2.0000
6.6660 -8.6667 3.6667
35.7689 -67.4375 -3.5551
no = 3
System is observabel!
(d) 求系统G(s)=的最小实现
程序如下:
num=[1 1];den=[1 10 27 18];
G=tf(num,den);
Gs=ss(G);
Gm=minreal(Gs);
Am=Gm.a
Bm=Gm.b
Cm=Gm.c
Dm=Gm.d
运行结果:
>> 1 state removed.
Am = 3.5391 -12.1540 Bm = 0.0606
5.1323 -12.5391 -0.2425
Cm = 0.2500 0.0625 Dm = 0
(2) 稳定性
(a) 代数法稳定性判据