现代控制理论控-实验报告.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.1056ik = 278.3081poles = -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.0000den = 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 1c = x1 x2 x3 d = u1 y1 1.5 1 0.5 y1 0Continuous-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 0Continuous-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; %输入状态矩阵Asyms 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; %输入状态矩阵Asyms s t;G=inv(s*eye(size(a)-a) %求预解矩阵eat2=ilaplace(G) %求拉氏反变换 G = (s+3)/(s2+3*s+2), 1/(s2+3*s+2) -2/(s2+3*s+2), s/(s2+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; %输入状态矩阵Asyms 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 %求零输入响应x1G1=inv(s*eye(size(a)-a)*bx2=ilaplace(G1/s) %求零状态响应x2x=x1+x2 %系统的状态响应x y=c*x %系统的输出响应yfor 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/23) 利用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)*x0G1=inv(s*eye(size(g)-g)*hx2=ilaplace(G1/s)x=x1+x2y=c*xfor 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*s2+25*s+4), 25/(25*s2+25*s+4) -4/(25*s2+25*s+4),25*s/(25*s2+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*s2+25*s+4)+25/(25*s2+25*s+4) -4/(25*s2+25*s+4)+25*s/(25*s2+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)-1x = 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)-1y = 253/60*exp(-4/5*t)-322/15*exp(-1/5*t)+73/44） 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.0000den = 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)*x0G1=inv(s*eye(size(a)-a)*b;x2=ilaplace(G1/s)x=x1+x2y=c*xfor 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)-1y = 5/4*exp(-2*t)+3/4 实验二 系统的能空性、能观测性、稳定性分析及实现实验日期 2013年5月21日 成绩：_批阅教师签字 1. 实验目的 加深理解能观测性、能控性、稳定性、最小实现等概念。掌握如何使用MATLAB进行以下分析和实现 系统的能观测性、能控性分析； 系统的稳定性分析； 系统的最小实现。2. 实验环境 MATLAB6.53. 实验内容、源程序代码、实验数据及结果分析 能控性、能观测性及系统实现(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!),endwo=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.7000nc = 0.4800System is controllable!wo = 1.2167 0.0500 0.0833 0.0500 1.2250 0.0500 0.0833 0.0500 0.0917no = 0.1253System 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!),endqo=obsv(a,c),no=rank(qo) %求系统能控性判别矩阵及其秩if n=no,disp(System is observable!), %判断系统的能控性else disp(System is uncontrollable!),end 运行结果: n = 3qc = 0 1 1 0 0 1 1 0 0 1 -11 -12 0 1 -11 -12 60 61nc = 3System is controllable!qo = 1 0 1 0 1 0 -6 -10 -6 0 0 1 36 60 26 -6 -11 -6no = 3System 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.4142Cc = 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 0Ao = -1.0000 -1.3416 -3.8341 Bo = 1.2247 0.0000 -0.4000 -0.7348 -0.5477 0 0.4899 -1.6000 -0.4472Co = 0 0.0000 -2.2361To = 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.aBm=Gm.bCm=Gm.cDm=Gm.d运行结果: 1 state removed.Am = 2.0121 -5.7764 Bm = 0.1212 3.4812 -7.0121 -0.4848Cm = 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!),endQo=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 = 3Qc = 1 -10 73 0 1 -10 0 0 1 nc = 3System is controllable!Qo = 0 1 -1 1 -1 0 -11 -27 -18 no = 3System 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 = 3Qc = 1 -10 73 0 1 -10 0 0 1 nc = 3System 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 = 3Qc = 1 -10 73 0 1 -10 0 0 1 nc = 3System is controllable!Qo = 0 1 1 1 1 0 -9 -27 -18 no = 2System 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!),endQo=obsv(A,C)no=rank(Qo)if n=no,disp(System is observabel!),else disp(System is unobservable!),end运行结果: n = 3Qc = 0 -11.0000 -84.9926 1.0000 1.0000 -8.0000 1.0000 3.0000 7.0000 nc = 3System is controllable!Qo = 1.0000 0 2.0000 6.6660 -8.6667 3.6667 35.7689 -67.4375 -3.5551no = 3System 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.aBm=Gm.bCm=Gm.cDm=Gm.d运行结果: 1 state removed.Am = 3.5391 -12.1540 Bm = 0.0606 5.1323 -12.5391 -0.2425Cm = 0.2500 0.0625 Dm = 0(2) 稳定性 (a) 代数法稳定性判据