MATLAB实验.docx

1、MATLAB实验实验一一、打开Matlab软件，点击，进入simulink Library Browser，点击文件-新建-Model 选择建立如下simulink模型图。修改单位阶跃模块的参数：step time:0 (阶跃信号起始时间为0时刻)这个值需要设置一下，其余的值都不用改动。intitial value:0（初值为0）；Final value:1（终值为1）；Sample time：0 （采样时间）（这一项我没太弄明白，变了几个值，得到的结果没什么变化。）将求和模块改为+-；将Transfer Fcn模块的参数设为num:10，den:0.02 0.3 1 0；（2）simulat

2、ion菜单star命令开始仿真；（也可使用那个小黑三角的图标工具）（3）双击示波器模块观察波形；（可以使用Autoscale工具，就是那个像望远镜的图标，得到最佳观测效果）二、选择建立如下simulink模型图。（做书上P139页到145页即可）最后生成答案：先画出simulink模型图，然后全选点击右键弹出菜单选Create Subsystem,生成Subsystem模块，然后选中Subsystem功能模块，右键弹出菜单选中Mask Subsystem进入Mask编辑窗口，如下图将模块下Subsystem名字改为disk(PID),点击OK生成PID Controller。实验二1、时域分析

3、（1）根据下面传递函数模型：绘制其单位阶跃响应曲线并在图上读标注出峰值，求出系统的性能指标。答案：程序shiyan3_1_1.m（本程序中用到了stepchar函数，需要按照书后附录D自己建立stepchar.m）function pos,tr,ts,tp=stepchar(g0,delta)y,t=step(g0);figer=plot(t,y),grid onmp,ind=max(y); dimt=length(t); yss=y(dimt);pos=100*(mp-yss)/ysstp=t(ind)for i=1:dimt if y(i)=yss tr=t(i); break; ende

4、nd;Tr=trfor i=1:length(y) if y(i)=(1+delta)*yss ts=t(i); endendTs=ts以下是主程序： G=tf(5 25 30,1 6 10 8);step(G)pos,tr,ts,tp=stepchar(G,0.02)figer = 1.9670e+003pos = 7.5374tp = 2.2086Tr = 1.4356Ts =3.6994（2）已知两个线性定常连续系统的传递函数分别为和，绘制它们的单位脉冲响应曲线。， 答案：程序（shiyan312.m）g1=tf(1 2 4,1 10 5 4);g2=tf(3 2,2 7 2);impu

5、lse(g1,g2) %脉冲响应曲线 运行结果：（3）已知线性定常系统的状态空间模型和初始条件，绘制其零输入响应曲线。 答案： %教材P184例题。A=-0.5572,-0.7814;0.7814,0;B=0;0;C=1.9691,6.4493;D=0;x0=1;0;t=0:0.1:20;initial(A,B,C,D,x0,t),grid on 2、频域分析 设线性定常连续系统的传递函数分别为、和，将它们的Bose图绘制在一张图中。， 答案：G1=tf(1,5 1);G2=tf(0.3,5 1 4);G3=tf(0.6,1 1);bode(G1,r-,G2,b-,G3,k:)运行结果：3、根

6、轨迹分析 根据下面负反馈系统的开环传递函数，绘制系统根轨迹，并分析系统稳定的K值范围。 答案：sys=tf(1,conv(1,0,conv(1 1,1 2);rlocus(sys) %P186K,ploes=rlocfind(sys)注:根轨迹图出来后，用鼠标点击竖的虚线和根轨迹的相交点。零极点在左平面时，系统稳定。故系统稳定的K值范围是0-6。selected_point = 0 + 1.4286iK = 6.1227ploes = -3.0111 0.0055 + 1.4260i 0.0055 - 1.4260i实验三 控制系统分析(2)4、系统稳定性分析（1）代数法稳定性判据：（用求分母

7、多项式的根和routh函数两种方法）已知系统的开环传递函数为：试对系统闭环判别其稳定性答案：要用routh函数，首先按照书后附录D建立routh . m如下function rtab,info=routh0(den)info=; vec1=den(1:2:length(den); nrT=length(vec1);vec2=den(2:2:length(den)-1); rtab=vec1; vec2, zeros(1,nrT-length(vec2);for k=1:length(den)-2, alpha(k)=vec1(1)/vec2(1); for i=1:length(vec2),

8、a3(i)=rtab(k,i+1)-alpha(k)*rtab(k+1,i+1); end if sum(abs(a3)=0 a3=polyder(vec2); info=info,All elements in row ,. int2str(k+2) are zeros; elseif abs(a3(1) G=zpk(-2,0,-1,-20,100)Zero/pole/gain: 100 (s+2)-s (s+1) (s+20) Gclose=feedback(G ,1) %单位负反馈，等同于feedback(G ,1,-1)Zero/pole/gain: 100 (s+2)-(s+3.10

9、1) (s+5) (s+12.9) G0=tf(Gclose)Transfer function: 100 s + 200-s3 + 21 s2 + 120 s + 200 den=1 21 120 200; rtab,msg=routh(den) rtab= 1 120 21 0 120 0 0 0msg =All elements in row 4 are zeros;（2）根轨迹法判断系统稳定性：已知一个单位负反馈系统开环传递函数为： 试在系统的闭环根轨迹图上选择一点，求出该点的增益及其系统的闭环极点位置，并判断在该点系统闭环的稳定性。num=1 3;den=conv(conv(1 0,

10、1 5),conv(1 6,1 2 2);G=tf(num,den);rlocus(G)K,p=rlocfind(G)Select a point in the graphics windowselected_point = -7.0912 + 3.6957iK = 1.2158e+003p = -7.1115 + 3.6800i -7.1115 - 3.6800i2.0754 + 3.8890i 2.0754 - 3.8890i -2.9277 再选一次：Select a point in the graphics windowselected_point = -0.2784 + 1.024

11、8iK = 15.9165p = -5.5352 + 0.2986i -5.5352 - 0.2986i -1.3640 -0.2828 + 1.0292i -0.2828-1.0292i（3）Bode图法判断系统稳定性：已知两个单位负反馈系统的开环传递函数分别为：用Bode图法判断系统闭环的稳定性。num1=2.7;den1=1 5 4 0;num2=2.7;den2=1 5 -4 0;g1=tf(num1,den1);g2=tf(num2,den2);margin(g1), grid onGm,Pm,Wcg,Wcp=margin(g1)Gm = 7.4074Pm = 51.7321Wcg

12、= 2.0000Wcp = 0.5783 Gm,Pm,Wcg,Wcp=margin(g2)Gm = InfPm = -58.0504Wcg = NaNWcp =0.53465、系统能控性、能观性分析已知连续系统的传递函数模型，当分别取1，0，1时，判别系统的能控性与能观性当等于-1时num=1 -1;den=1 10 27 18;a,b,c,d=tf2ss(num,den);tc=ctrb(a,b) %求能控矩阵rank(tc) %能控矩阵的秩to=obsv(a,c) %求能观矩阵rank(to) %能观矩阵的秩执行结果：tc = 1 -10 73 0 1 -10 0 0 1ans = 3to

13、 = 0 1 -1 1 -1 0 -11 -27 -18ans = 3能观、能控矩阵都是满秩。故能观又能控（2）当等于0时num=1 0;den=1 10 27 18;a,b,c,d=tf2ss(num,den)tc=ctrb(a,b) rank(tc) to=obsv(a,c) rank(to) 执行结果：tc = 1 -10 73 0 1 -10 0 0 1ans = 3to = 0 1 0 1 0 0 -10 -27 -18ans = 3能观、能控矩阵都是满秩。故能观又能控（3）当等于1时num=1 1;den=1 10 27 18;a,b,c,d=tf2ss(num,den)tc=ct

14、rb(a,b) %求能控矩阵rank(tc) %能控矩阵的秩to=obsv(a,c) %求能观矩阵rank(to) %能观矩阵的秩执行结果：tc =1 -10 73 0 1 -10 0 0 1ans =3to = 0 1 1 1 1 0 -9 -27 -18ans = 2所以能控，不能观实验三 基于BODE图法的控制系统设计1、已知单位负反馈控制系统的开环传递函数为：设计超前校正装置，使校正后系统满足：， 答案：为留有裕量取编写Matlab程序：KK=100;wc=50;Pm=50;ng0= KK *1; dg0=conv(1,0,conv(0.01,1,0.1,1); g0=tf(ng0,d

15、g0) ;t=0:0.001:0.5;w=logspace(-1,3);ngc,dgc= fg_lead_pm_wc(ng0,dg0,Pm,wc,w); %ngc,dgc= fg_lead_pd(ng0,dg0,wc)gc=tf(ngc,dgc),g0c=tf(g0*gc);b1=feedback(g0,1); b2=feedback(g0c,1);step(b1,r-, b2,b,t); grid on Transfer function:0.08578 s + 1-0.002402 s + 1figure, bode(g0,r-,g0c,b,w), grid on, pos,tr,ts,t

16、p=stepchar(b2,0.02) gm,pm,wcg,wcp=margin(g0c), Km=20*log10(gm)figer = 352.0022pos = 24.9651tp = 0.0409Tr = 0.0289Ts = 0.0939pos = 24.9651tr = 0.0289ts = 0.0939tp = 0.0409gm = 5.8934pm = 44.3060wcg = 201.9342wcp = 69.6644Km = 15.4074函数程序：stepchar、fg_lead_pm_wc（带惯性环节的超前控制器）、fg_lead_pd（不带惯）见书后附录，需要同学们先

17、编辑好保存成.m文件。两个函数均可完成实验要求。2、已知单位负反馈控制系统的开环传递函数为：设计超前校正装置，使校正后系统满足：， 答案： 为留有裕量取编写matlab程序：KK=100; Pm=40; wc=5;ng0=KK*1; dg0=conv(1,0,conv(0.1,1,0.01,1);g0=tf(ng0,dg0);t=0:0.001:2; w=logspace(-1,4);ngc,dgc= fg_lag_wc (ng0,dg0,w, wc);gc=tf(ngc,dgc),g0c=tf(g0*gc);b1=feedback(g0,1); b2=feedback(g0c,1);step

18、(b1,r-, b2,b,t); axis(0,5,-1,1.8); grid on Transfer function: 2 s + 1-35.73 s + 1figure, bode(g0,r-,g0c,b,w); grid on, pos,tr,ts,tp=stepchar(b2,0.02)gm,pm,wcg,wcp=margin(g0c)figer = 352.0068pos = 16.6273tp = 0.5803Tr = 0.3869Ts = 2.9014pos = 16.6273tr = 0.3869ts = 2.9014tp = 0.5803gm = 18.6322pm = 5

19、5.0990wcg = 30.7909wcp = 5.0205函数程序：fg_lag_wc见书后附录。stepchar.mfunction pos,tr,ts,tp=stepchar(g0,delta)y,t=step(g0);mp,ind=max(y); dimt=length(t); yss=y(dimt);pos=100*(mp-yss)/yss;tp=t(ind);for i=1:dimt if y(i)=yss tr=t(i); break; endend;for i=1:length(y) if y(i)=(1+delta)*yss ts=t(i); endendfg_lead_p

20、m_wc.mfunction ngc,dgc=fg_lead_pm_wc(ng0,dg0,Pm,wc,w)mu,pu=bode(ng0,dg0,w);ngv=polyval(ng0,j*wc);dgv=polyval(dg0,j*wc);g=ngv/dgv; %求原系统在期望的剪切频率处的频率响应数据G0(jwc) theta=180*angle(g)/pi; %原系统在期望的剪切频率处的相角裕度，单位为度alf=ceil(Pm-(theta+180)+5); %最大超前角phi=(alf)*pi/180;a=(1+sin(phi)/(1-sin(phi);dbmu=20*log10(mu);

21、 mm=-10*log10(a);wgc=spline(dbmu,w,mm); T=1/(wgc*sqrt(a);ngc=a*T,1;dgc=T,1;fg_lead_pd.mfunction ngc,dgc=fg_lead_pd(ng0,dg0,wc)ngv=polyval(ng0,j*wc); dgv=polyval(dg0,j*wc);g=ngv/dgv; mg0=abs(g); t=sqrt(1/mg0)2-1)/(wc2); %幅值相加为零ngc=t,1;dgc=1; %Gc(s)=1+Tsfg_lag_wc.mfunction ngc,dgc=fg_lag_wc(ng0,dg0,w,

22、wc)ngv=polyval(ng0,j*wc);dgv=polyval(dg0,j*wc); g=ngv/dgv;alph=abs(1/g); T=10/(alph*wc);ngc=alph*T,1;dgc=T,1;3、已知控制系统方框图如图3-1所示。图中：试设计反馈校正环节Hc（s）,以满足下列要求：，, 。，。 书上p225页例题KK=50;bp=0.3;ts=2.5;delta=0.05;ng0=20;dg0=conv(1,0,conv(0.9,1,0.007,1);g0=tf(ng0,dg0);w=logspace(-4,3);t=0:0.1:3;mag,phase=bode(KK

23、*g0,w);gm0,pm0,wg0,wc0=margin(mag,phase,w),gm0=20*log10(gm0) mr=0.6+2.5*bpwc=ceil(2+1.5*(mr-1)+2.5*(mr-1)2)*pi/ts),h=(mr+1)/(mr-1)w1=2*wc/(h+1),w2=h*w1 w1=wc/10;w2=9;ng1=1/w1,1;dg1=conv(1/w2,1,conv(1,0,1,0);g1=tf(ng1,dg1);g=polyval(ng1,j*wc)/polyval(dg1,j*wc);K=abs(1/g);g1=tf(K*g1) h=tf(dg1,ng1);Kh=

24、1/K;h=tf(Kh*h) g2=feedback(KK*g0,h);b1=feedback(g0,1);b2=feedback(g2,1);bode(KK*g0,r-,g2,b,h,g,w);grid onstep(b1,r-,b2,b,t);grid onpos,tr,ts,tp=stepchar(b2.delta)实验五 控制系统的PID设计和工程整定方法编写一个Ziegler()函数如下：function Gc,Kp,Ti,Td,H=Ziegler(key,vars)Ti=;Td=;H=;if length(vars)=3 K=vars(1);Tc=vars(2);N=vars(3)

25、; if key=1,Kp=0.5*K; elseif key=2,Kp=0.4*K;Ti=0.85*Tc; elseif key=3,Kp=0.6*K;Ti=0.5*Tc;Td=0.125*Tc; endendswitch key case 1 Gc=Kp; case 2 Gc=tf(Kp*Ti,1,Ti 0); case 3 num=Kp*Ti*Td*(N+1)/N,Kp*(Ti+Td/N),Kp; den=Ti*Td/N,1,0; Gc=tf(num,den);end调用函数Ziegler()，可设计各个控制器。程序如下：G=tf(10,1 10 35 50 24);Kc,pp,wg,w

26、p=margin(G);Tc=2*pi/wg;Gc1,Kp1=Ziegler(1,Kc,Tc,10);disp(P控制器参数Kp1为：)Kp1Gc2,Kp2,Ti2=Ziegler(2,Kc,Tc,10);disp(PI控制器参数Kp2和Ti2为：)Kp2Ti2Gc3,Kp3,Ti3,Td3=Ziegler(3,Kc,Tc,10);disp(PID控制器参数Kp3,Ti3和Td3为：)Kp3Ti3Td3G_c1=feedback(G*Gc1,1);G0=feedback(G,1);G_c2=feedback(G*Gc2,1);G_c3=feedback(G*Gc3,1);step(G0,.,G_c1,r,G_c2,g,G_c3,b)实验六 基于SIMULINK的控制系统综合仿真1、建立模型图：（tu6_1）运行matlab程序如下：t=0:0.1:9.9;for i=1:6 ut=t,i*ones(size(t); tt,xx,yy=sim(tu6_1,10,ut); hold on; plot(tt,yy);grid onend2、建立模型图：（tu6_2）运行matlab程序如下：t=0:20:2000;ut=t,ones(size(t);tt,xx,yy=sim(tu6_2,2000,ut);plot(tt,yy);grid on