MATLAB语言与控制系统仿真-参考答案-第5章
控制系统计算机辅助设计：MATLAB语言与应用(第2版)薛定宇-课后习题答案
第1章控制系统计算机辅助设计概述第2章 MATLAB语言程序设计基础第3章线性控制系统的数学模型第4章线性控制系统的计算机辅助分析第5章 Simulink在系统仿真中的应用第6章控制系统计算机辅助设计第1章控制系统计算机辅助设计概述【1】已阅，略【2】已阅，略【3】已经掌握help命令和Help菜单的使用方法【4】区别：MATLAB语言实现矩阵的运算非常简单迅速，且效率很高，而用其他通用语言则不然，很多通用语言所实现的矩阵运算都是对矩阵维数具有一点限制的，即使限制稍小的，但凡维数过大，就会造成运算上的溢出出错或者运算出错，甚至无法处理数据的负面结果【5】【8】(1)输入激励为正弦信号(2)输入激励为脉冲模拟信号(3)输入激励为时钟信号(4) 输入激励为随机信号(5) 输入激励为阶跃信号δ=0.3δ=0.05δ=0.7结论：随着非线性环节的死区增大，阶跃响应曲线的范围逐渐被压缩，可以想象当死区δ足够大时，将不再会有任何响应产生。
所以可以得到结论，在该非线性系统中，死区的大小可以改变阶跃响应的幅值和超调量。
死区越大，幅值、超调量将越小，而调整时间几乎不受其影响第2章 MATLAB语言程序设计基础【1】>> A=[1 2 3 4;4 3 2 1;2 3 4 1;3 2 4 1]A =1 2 3 44 3 2 12 3 4 13 24 1>> B=[1+4i,2+3i,3+2i,4+i;4+i,3+2i,2+3i,1+4i;2+3i,3+2i,4+i,1+4i;3+2i,2+3i,4+i,1+4i]B =1.0000 + 4.0000i2.0000 +3.0000i 3.0000 + 2.0000i4.0000 + 1.0000i4.0000 + 1.0000i 3.0000 + 2.0000i 2.0000 + 3.0000i 1.0000 + 4.0000i2.0000 +3.0000i 3.0000 + 2.0000i4.0000 + 1.0000i 1.0000 + 4.0000i3.0000 + 2.0000i 2.0000 + 3.0000i4.0000 + 1.0000i 1.0000 + 4.0000i>> A(5,6)=5A =1 2 3 4 0 04 3 2 1 0 02 3 4 1 0 03 24 1 0 00 0 0 0 0 5∴若给出命令A(5,6)=5则矩阵A的第5行6列将会赋值为5，且其余空出部分均补上0作为新的矩阵A，此时其阶数为5×6【2】相应的MATLAB命令：B=A(2:2:end,:)>> A=magic(8)A =64 2 3 61 60 6 7 579 55 54 12 13 51 50 1617 47 46 20 21 43 42 2440 26 27 37 36 30 31 3332 34 35 29 28 38 39 2541 23 22 44 45 19 18 4849 15 14 52 53 11 10 568 58 59 5 4 62 63 1>> B=A(2:2:end,:)B =9 55 54 12 13 51 50 1640 26 27 37 36 30 31 3341 23 22 44 45 19 18 488 58 59 5 4 62 63 1∴从上面的运行结果可以看出，该命令的结果是正确的【3】>> syms x s; f=x^5+3*x^4+4*x^3+2*x^2+3*x+6f =x^5 + 3*x^4 + 4*x^3 + 2*x^2 + 3*x + 6>> [f1,m]=simple(subs(f,x,(s-1)/(s+1)))f1 =19 - (72*s^4 + 120*s^3 + 136*s^2 + 72*s + 16)/(s + 1)^5m =simplify(100)【4】>> i=0:63; s=sum(2.^sym(i))s =0615【5】>> for i=1:120if(i==1|i==2) a(i)=1;else a(i)=a(i-1)+a(i-2);endif(i==120) a=sym(a); disp(a); endend[ 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987, 1597, 2584, 4181, 6765, 10946, 17711, 28657, 46368, 75025, 121393, 196418, 317811, 514229, 832040, 1346269, 2178309, 3524578, 5702887, 9227465, , , , , 5, 1, 6, 7, 3, 70, 03, 73, 76, 49, , 074, 099, 173, 272, 2445, 3717, 6162, 9879, 6041, 55920, 81961, 37881, 19842, 106, 177565, 035288, 212853, 248141, 0460994, , 1170129, 1879264, 8065, , , , 00884757, , 0, 5, 6, 1, 0, 88, , 673, 58, 931, , 120, , 029, 4, 2, 9905, 3072, 2977, 46049, 69026, 15075, 40, 99176, 083277, 082453, 165730, 248183, 7576, 62096, , 4738105, 5814114, 9, 186333, , 284885, 9, 3488322, 9, 0, 0]【6】>>k=1;for i=2:1000for j=2:iif rem(i,j)==0if j<i, break;endif j==i, A(k)=i; k=k+1; break; endendendenddisp(A);Columns 1 through 132 3 5 7 11 13 17 19 23 29 31 37 41 Columns 14 through 2643 47 53 59 61 67 71 73 79 83 89 97 101 Columns 27 through 39103 107 109 113 127 131 137 139 149 151 157 163 167 Columns 40 through 52173 179 181 191 193 197 199 211 223 227 229 233 239 Columns 53 through 65241 251 257 263 269 271 277 281 283 293 307 311 313 Columns 66 through 78317 331 337 347 349 353 359 367 373 379 383 389 397 Columns 79 through 91401 409 419 421 431 433 439 443 449 457 461 463 467 Columns 92 through 104479 487 491 499 503 509 521 523 541 547 557 563 569 Columns 105 through 117571 577 587 593 599 601 607 613 617 619 631 641 643 Columns 118 through 130647 653 659 661 673 677 683 691 701 709 719 727 733 Columns 131 through 143739 743 751 757 761 769 773 787 797 809 811 821 823 Columns 144 through 156827 829 839 853 857 859 863 877 881 883 887 907 911 Columns 157 through 168919 929 937 941 947 953 967 971 977 983 991 997【7】说明：h和D在MATLAB中均应赋值，否则将无法实现相应的分段函数功能syms x; h=input(‘h=’); D=input(‘D=’);y=h.*(x>D)+(h.*x/D).*(abs(x)<=D)-h.*(x<-D)【10】function y=fib(k)if nargin~=1,error('出错：输入变量个数过多，输入变量个数只允许为1！');endif nargout>1,error('出错：输出变量个数过多！');endif k<=0,error('出错：输入序列应为正整数！');endif k==1|k==2,y=1;else y=fib(k-1)+fib(k-2);endend【13】【14】>> t=[-1:0.001:-0.2,-0.1999:0.0001:0.1999,0.2:0.001:1];y=sin(1./t);plot(t,y);grid on;【15】(1) >> t=-2*pi:0.01:2*pi;r=1.0013*t.^2;polar(t,r);axis('square')(2) >> t=-2*pi:0.001:2*pi;r=cos(7*t/2);polar(t,r);axis('square')(3) >> t=-2*pi:0.001:2*pi;r=sin(t)./t;polar(t,r);axis('square')(4) >> t=-2*pi:0.001:2*pi;r=1-cos(7*t).^3;polar(t,r);axis('square')【17】(1)z=xy>> [x,y]=meshgrid(-3:0.01:3,-3:0.01:3);z=x.*y;mesh(x,y,z);>> contour3(x,y,z,50);(1)z=sin(xy)>> [x,y]=meshgrid(-3:0.01:3,-3:0.01:3);z=sin(x.*y);mesh(x,y,z);>> contour3(x,y,z,50);第3章线性控制系统的数学模型【1】(1) >> s=tf('s');G=(s^2+5*s+6)/(((s+1)^2+1)*(s+2)*(s+4))Transfer function:s^2 + 5 s + 6--------------------------------s^4 + 8 s^3 + 22 s^2 + 28 s + 16(2) >> z=tf('z',0.1);H=5*(z-0.2)^2/(z*(z-0.4)*(z-1)*(z-0.9)+0.6)Transfer function:5 z^2 - 2 z + 0.2---------------------------------------z^4 - 2.3 z^3 + 1.66 z^2 - 0.36 z + 0.6Sampling time (seconds): 0.1【2】(1)该方程的数学模型>> num=[6 4 2 2];den=[1 10 32 32];G=tf(num,den)Transfer function:6 s^3 + 4 s^2 + 2 s + 2------------------------s^3 + 10 s^2 + 32 s + 32(2)该模型的零极点模型>> G=zpk(G)Zero/pole/gain:6 (s+0.7839) (s^2 - 0.1172s + 0.4252)-------------------------------------(s+4)^2 (s+2)(3)由微分方程模型可以直接写出系统的传递函数模型【5】(1) >> P=[0;0;-5;-6;-i;i];Z=[-1+i;-1-i];G=zpk(Z,P,8)Zero/pole/gain:8 (s^2 + 2s + 2)-------------------------s^2 (s+5) (s+6) (s^2 + 1)(2) P=[0;0;0;0;0;8.2];Z=[-3.2;-2.6];H=zpk(Z,P,1,'Ts',0.05,'Variable','q')Zero/pole/gain:(q+3.2) (q+2.6)---------------q^5 (q-8.2)Sampling time (seconds): 0.05【8】(1)闭环系统的传递函数模型>> s=tf('s');G=10/(s+1)^3;Gpid=0.48*(1+1/(1.814*s)+0.4353*s/(1+0.4353*s));G1=feedback(Gpid*G,1)Transfer function:7.58 s^2 + 10.8 s + 4.8--------------------------------------------------------------0.7896 s^5 + 4.183 s^4 + 7.811 s^3 + 13.81 s^2 + 12.61 s + 4.8 (2)状态方程的标准型实现>> G1=ss(G1)a =x1 x2 x3 x4 x5 x1 -5.297 -2.473 -2.186 -0.9981 -0.7598x2 4 0 0 0 0 x3 0 2 0 0 0 x4 0 0 2 0 0x5 0 0 0 0.5 0b =u1x1 2x2 0x3 0x4 0x5 0c =x1 x2 x3 x4 x5y1 0 0 0.6 0.4273 0.3799d =u1y1 0Continuous-time state-space model.(3)零极点模型>> G1=zpk(G1)Zero/pole/gain:9.6 (s^2 + 1.424s + 0.6332)--------------------------------------------------------(s+3.591) (s^2 + 1.398s + 0.6254) (s^2 + 0.309s + 2.707)【11】>> Ga=feedback(s/(s^2+2)*1/(s+1),(4*s+2)/(s+1)^2);Gb=feedback(1/s^2,50);G=3*feedback(Gb*Ga,(s^2+2)/(s^3+14))Transfer function:3 s^6 + 6 s^5 + 3 s^4 + 42 s^3 + 84 s^2 + 42 s---------------------------------------------------------------------------s^10 + 3 s^9 + 55 s^8 + 175 s^7 + 300 s^6 + 1323 s^5 + 2656 s^4 + 3715 s^3+ 7732 s^2 + 5602 s + 1400【13】c1=feedback(G5*G4,H3)=G5*G4/(1+G5*G4*H3)c2=feedback(G3,H4*G4)=G3/(1+G3*H4*G4)c3=feedback(c2*G2,H2)=c2*G2/(1+c2*G2*H2)=G3*G2/(1+G3*H4*G4+G3*G2*H1)G=feedback(G6*c1*c3*G1,H1)=G6*c1*c3*G1/(1+ G6*c1*c3*G1*H1)=G6*G5*G4*G3*G2*G1/(1+G3*H4*G4+G3*G2*H1+G5*G4*H3+G5*G4*H3*G3*H4*G4+G5*G4* H3*G3*G2*H1+G6*G5*G4*G3*G2*G1*H1)【14】>> s=tf('s');c1=feedback(0.21/(1+0.15*s),0.212*130/s);c2=feedback(c1*70/(1+0.0067*s)*(1+0.15*s)/(0.051*s),0.1/(1+0.01*s));G=(1/(1+0.01*s))*feedback(130/s*c2*1/(1+0.01*s)*(1+0.17*s)/(0.085*s),0.0044/(1+0.01*s)) Transfer function:0.004873 s^5 + 1.036 s^4 + 61.15 s^3 + 649.7 s^2 + 1911 s---------------------------------------------------------------------------4.357e-014 s^10 + 2.422e-011 s^9 +5.376e-009 s^8 +6.188e-007 s^7+ 4.008e-005 s^6 + 0.001496 s^5 + 0.03256 s^4 + 0.4191 s^3+ 2.859 s^2 + 8.408 s 第4章线性控制系统的计算机辅助分析【1】(1) >> num=[1];den=[3 2 1 2];G=tf(num,den);eig(G)ans =-1.00000.1667 + 0.7993i0.1667 - 0.7993i分析：由以上信息可知，系统的极点有2个是在s域的右半平面的，因此系统是不稳定的(2) >> num=[1];den=[6 3 2 1 1];G=tf(num,den);eig(G)ans =-0.4949 + 0.4356i-0.4949 - 0.4356i0.2449 + 0.5688i0.2449 - 0.5688i分析：由以上信息可知，系统的极点有2个是在s域的右半平面的，因此系统是不稳定的(3) >> num=[1];den=[1 1 -3 -1 2];G=tf(num,den);eig(G)ans =-2.0000-1.00001.00001.0000分析：由以上信息可知，系统的极点有2个是在s域的右半平面的，因此系统是不稳定的(4) >> num=[3 1];den=[300 600 50 3 1];G=tf(num,den);eig(G)ans =-1.9152-0.14140.0283 + 0.1073i0.0283 - 0.1073i分析：由以上信息可知，系统的极点有2个是在s域的右半平面的，因此系统是不稳定的(5) >> s=tf('s');G=0.2*(s+2)/(s*(s+0.5)*(s+0.8)*(s+3)+0.2*(s+2));eig(G)ans =-3.0121-1.0000-0.1440 + 0.3348i-0.1440 - 0.3348i分析：由以上信息可知，系统的所有极点都在s域的左半平面，因此系统是稳定的【2】(1) >> num=[-3 2];den=[1 -0.2 -0.25 0.05];H=tf(num,den,'Ts',0.5);abs(eig(H)')ans =0.5000 0.5000 0.2000分析：由以上信息可知，所有特征根的模均小于1，因此该系统是稳定的(2) >> num=[3 -0.39 -0.09];den=[1 -1.7 1.04 0.268 0.024];H=tf(num,den,'Ts',0.5);abs(eig(H)')ans =1.1939 1.1939 0.1298 0.1298分析：由以上信息可知，由于前两个特征根的模均大于1，因此该系统是不稳定的(3) >> num=[1 3 -0.13];den=[1 1.352 0.4481 0.0153 -0.01109 -0.001043];H=tf(num,den,'Ts',0.5);abs(eig(H)')ans =0.8743 0.1520 0.2723 0.2344 0.1230分析：由以上信息可知，所有特征根的模均小于1，因此该系统是稳定的(4) >> num=[2.12 11.76 15.91];den=[1 -7.368 -20.15 102.4 80.39 -340];H=tf(num,den,'Ts',0.5,'Variable','q');abs((eig(H))')ans =8.2349 3.2115 2.3415 2.3432 2.3432分析：由以上信息可知，所有特征根的模均大于1，因此该系统是不稳定的【3】(1) >> A=[-0.2,0.5,0,0,0;0,-0.5,1.6,0,0;0,0,-14.3,85.8,0;0,0,0,-33.3,100;0,0,0,0,-10];eig(A)ans =-0.2000-0.5000-14.3000-33.3000-10.0000分析：由以上信息可知，该连续线性系统的A矩阵的所有特征根的实部均为负数，因此该系统是稳定的(2)>>F=[17,24.54,1,8,15;23.54,5,7,14,16;4,6,13.75,20,22.5589;10.8689,1.2900,19.099,…-4-3.5-3-2.5-2-1.5-1-0.50x 10-6P ole-Zero Map Real Axis (seconds -1)I m a g i n a r y A x i s (s e c o n d s -1)21.896,3;11,18.0898,25,2.356,9];abs(eig(F)') ans =63.7207 23.5393 12.4366 13.3231 19.7275分析：由以上信息可知，该离散系统的F 矩阵的所有特征根的模均大于1，因此该系统是不稳定的 【4】>> A=[-3 1 2 1;0 -4 -2 -1;1 2 -1 1;-1 -1 1 -2]; B=[1 0;0 2;0 3;1 1];C=[1 2 2 -1;2 1 -1 2];D=[0 0;0 0];G=ss(A,B,C,D); tzero(G)pzmap(G)ans =-3.6124-1.2765结论：∴可以得到该系统的 零点为-3.6124、-1.2765 分析：由以上信息可知，【5】>> s=tf('s');Gc=sscanform(G,'ctrl') Go=sscanform(G,'obsv') a =x1 x2 x3 x4 x1 0 1 0 0 x2 0 0 1 0 x3 0 0 0 1 x4 -0.4 -1.4 -4.3 -4.3 b =u1 x1 0 x2 0 x3 0 x4 1 c =x1 x2 x3 x4 y1 0.4 0.2 0 0 d =u1 y1 0Continuous-time state-space model. a =x1 x2 x3 x4x1 0 0 0 -0.4x2 1 0 0 -1.4x3 0 1 0 -4.3x4 0 0 1 -4.3b =u1x1 0.4x2 0.2x3 0x4 0c =x1 x2 x3 x4y1 0 0 0 1d =u1y1 0Continuous-time state-space model.【9】(1)>> num=[18 514 5982 36380 122664 222088 185760 40320];den=[1 36 546 4536 22449 67284 118124 109584 40320];[R1,P1,K1]=residue(num,[den 0]);[R1,P1]ans =-1.2032 -8.0000-1.0472 -7.00000.2000 -6.00000.7361 -5.0000-2.8889 -4.00002.2250 -3.0000-2.0222 -2.00003.0004 -1.00001.0000 0>> [n,d]=rat(R1);sym([n./d]')ans =[ -379/315, -377/360, 1/5, 53/72, -26/9, 89/40, -91/45, 7561/2520, 1][阶跃响应的解析解]y(t)=(-379/315)*e-8t+(-377/360)*e-7t+(1/5)*e-6t+(53/72)*e-5t+(-26/9)*e-4t+(89/40)*e-3t+ (-90/45)*e-2t+(7561/2520)*e-t+1(2) >> num=[18 514 5982 36380 122664 222088 185760 40320];den=[1 36 546 4536 22449 67284 118124 109584 40320];[R2,P2,K2]=residue(num,den);[R2,P2]ans =9.6254 -8.00007.3306 -7.0000-1.2000 -6.0000-3.6806 -5.000011.5556 -4.0000-6.6750 -3.00004.0444 -2.0000-3.0004 -1.0000>> [n,d]=rat(R2);sym([n./d]')ans =[ 3032/315, 887/121, -6/5, -265/72, 104/9, -267/40, 182/45, -7561/2520][脉冲响应的解析解]y(t)=(3032/315)*e-8t+(887/121)*e-7t+(-6/5)*e-6t+(-265/72)*e-5t+(104/9)*e-4t+(-267/40)*e-3t+(182/45)*e-2t+(-7561/2520)*e-t(3) >> syms t;u=sin(3*t+5);Us=laplace(u)Us =(3*cos(5) + s*sin(5))/(s^2 + 9)>> s=tf('s');Us=(3*cos(5)+s*sin(5))/(s^2+9);num=[18 514 5982 36380 122664 222088 185760 40320];den=[1 36 546 4536 22449 67284 118124 109584 40320];G=tf(num,den); Y=Us*G;num=Y.num{1}; den=Y.den{1};[R3,P3,K3]=residue(num,den); [R3,P3]ans =1.1237 -8.00000.9559 -7.0000-0.1761 -6.0000-0.6111 -5.00002.1663 -4.0000-1.1973 - 0.0010i 0.0000 + 3.0000i-1.1973 + 0.0010i 0.0000 - 3.0000i-1.3824 -3.00000.8614 -2.0000-0.5430 -1.0000>> [n,d]=rat(R3);sym([n./d]')ans =[109/97, 282/295, -59/335, -965/1579, 951/439, - 449/375 + (18*i)/17981, - 449/375 - (18*i)/17981, -1663/1203, 317/368, -82/151]Linear Simulation Results Time (seconds)A m p l i t u d e [正弦信号时域响应的解析解]y(t)=(109/97)*e -8t +(282/295)*e -7t +(-59/335)*e -6t +(-965/1579)*e -5t +(-449/375)*e -4t +(-1663/1203)*e -3t +(317/368)*e -2t +(-82/151)*e -t -2.3947sin(3t)[输出波形]>> num=[18 514 5982 36380 122664 222088 185760 40320];den=[1 36 546 4536 22449 67284 118124 109584 40320]; G=tf(num,den); t=[1:.1:20]';u=sin(3*t+5); lsim(G,u,t);分析：由解析解可知，输出信号的稳态部分是振荡的，并且其幅值与相位始终 在到达稳态的时候保持不变，因此 右图所示的输出波形与解析解所得的结论是一致的【10】(1)因为PI 或PID 控制器均含有Ki/s 节，则当Kp →∞，即|e(t)|一环节后，如果要求|e(t)|→0(2)不稳定系统能用PI 或PID 在积分控制中，控制器的输出与输入误差信号的积分成正比关系。
计算机控制与仿真技术(第二版)-第5章 控制系统的MATLAB 7建模
5.2.2 环节并联连接的化简
环节并联是指多个环节的输入信号相同，所有环
节输出的代数和为其总输出。两环节相并联如图5-2
所示。
G1 ( s ) y1 (t )
+
u (t )
y (t )
G2 ( s)
图 5-2
y 2 (t )
+
并联连接结构
11
parallel()函数命令可以将两个环节并联连接进行 等效化简。它既适用于连续时间系统也适用于离散时
feedback()函数命令格式为： G=feedback(G1,G2,sign)
[a,b,c,d]=feedback(a1,b1,c1,d1,a2,b2,c2,d2,sign)
[num,den]=feedback(num1,den1,num2,den2,sign)
14
5.2.4 闭环系统的化简
通过将所有的输出反馈到输入，从而构成闭
num [c 0 , c1 , , c m ] den [a1 , a 2 , , a a ]
3
在MATLAB中，用函数命令tf()来建立控制系 统的传递函数模型，tf()函数命令常用的调用格式
为：
sys= tf(num, den)
sys= tf(num, den, Ts)
5.2 环节方框图模型的化简
5.2.1 环节串联连接的化简
控制系统的环节串联及其化简就是模块方框图模 型的串联及其化简。可以用MATLAB的函数命令 series()将串联模块化简。
u1 (t )
G1 ( s ) y1 (t )
u 2 (t )
G 2 ( s)
y 2 (t )
图 5-1
MATLAB语言与控制系统仿真参考答案第5章
5、6控制系统的时域响应MATLAB 仿真实训5、6、1实训目的 123 4 5学会利用MATLAB绘制系统的单位阶跃响应曲线 ，掌握读取系统动态性能指标的方法 学会利用MATLAB 绘制系统的单位脉冲响应曲线的方法 掌握利用MATLAB 绘制系统的零输入响应曲线的方法 掌握利用MATLAB 绘制系统的一般输入响应曲线的方法 学会通过仿真曲线读取相关信息 ，并依据有关信息进行系统的时域分析。
5、6、2实训内容 1、编写程序求取下列各系统的单位阶跃响应 ，完成表5-5并记录相关曲线。
G(s) 16 ~2 s 3.2s 16 G 2(s) 16~2 s 2.4s 16 G a (s)16 s 21.6s 16G 4(s)16 -2s s 16解:>> n1=16; >> d 仁[1,3、2,16]; >> sys1=tf( n1,d1); >> step(sys1)-X r倚 25 3 7百'imeQ -4o.a■rp3±±-i_l r v>> n2=16; >> d2=[1,2、4,16]; >> sys2=tf( n2,d2); >> step(sys2)041.4D Q粘EEdE'luaJ 05 1 15 2 2.3 3 55 4<1.5Tim (lie)11-o02:>> n3=16;>> d3=[1,1、6,16];>> sys3=tf( n3,d3);>> step(sys3)lim* (t«c) >> n4=16; >> d4=[1,1,16];>> sys4=tf( n4,d4);>> step(sys4)表5-5序号nc maxt pt s (5%)计算值实验计算值实验计算值实验值1 0、4 4 1、2538 1、250、85690、8632、 18752、12 0、34 1 3723 1、370、 82330、8282、 91672、813 0、24 1、52661、530、80160、8 4 3750 4、940、125 4 1、6731 1、670、 79160、8037、 00007、33>> z1=0、4; w=4;cmax 1=1+exp(-z1*pi/sqrt(1-z1A2)); tp1=pi/(w*sqrt(1-z1A 2)); ts1=3、5/(z1*w); [cmax1,tp1,ts1] ans =>> z2=0、3;w=4;cmax2=1+exp(-z2*pi/sqrt(1-z2A2)); tp2=pi/(w*sqrt(1-z2A2)); ts2=3、5/(z2*w); [cmax2,tp2,ts2]ELIIds1、 2538 0、 85692、 1875ans =1、3723 0、82332、 9167>> z3=0、2;w=4;cmax3=1+exp(-z3*pi/sqrt(1-z3A2));tp3=pi/(w*sqrt(1-z3A2));ts3=3、5/(z3*w);[cmax3,tp3,ts3] ans =1、5266 0、8016 4、3750 >> z4=0、125;w=4;cmax4=1+exp(-z4*pi/sqrt(1-z4A2));tp4=pi/(w*sqrt(1-z4A2));ts4=3、5/(z4*w);[cmax4,tp4,ts4]ans =1、6731 0、7916 7、0000说明：对于二阶欠阻尼系统（0 1）,若系统的闭环传递函数为则系统单位阶跃响应的输出最大值cmaxt pn 1 2调整时间估算值t s3.5n（以5%为误差带）t s 4.4（以2%为误差带）2、已知二阶系统的闭环传递函数如下，编程求取系统的单位阶跃响应并完成表峰值时间曲线。
MATLAB语言基础与应用(第二版)第5章 习题答案
第5章习题与答案5.1用矩阵三角分解方法解方程组123123123214453186920x x x x x x x x x +-=⎧⎪-+=⎨⎪+-=⎩ 解答：>>A=[2 1 -1;4 -1 3;6 9 -1] A =2 1 -1 4 -13 6 9 -1 >>b=[14 18 20]; b =14 18 20 >> [L, U, P]=lu(A) L =1.0000 0 0 0.6667 1.0000 0 0.3333 0.2857 1.0000 U =6.0000 9.0000 -1.0000 0 -7.0000 3.6667 0 0 -1.7143 P =0 0 1 0 1 0 1 0 0 >> y=backsub(L,P*b’) y =20.0000 4.6667 6.0000 >> x=backsub(U,y) x =6.5000 -2.5000 -3.5000 5.2 Cholesky 分解方法解方程组123121332352233127x x x x x x x ++=⎧⎪+=⎨⎪+=⎩ 解答：>> A=[3 2 3;2 2 0;3 0 12] A =3 2 32 2 03 0 12>> b=[5;3;7]b =537>> L=chol(A)L =1.7321 1.1547 1.73210 0.8165 -2.44950 0 1.7321>> y=backsub(L,b)y =-11.6871 15.7986 4.0415>> x=backsub(L',y)x =-6.7475 28.8917 49.93995.3解答:观察数据点图形>> x=0:0.5:2.5x =0 0.5000 1.0000 1.5000 2.0000 2.5000 >> y=[2.0 1.1 0.9 0.6 0.4 0.3]y =2.0000 1.1000 0.9000 0.6000 0.4000 0.3000 >> plot(x,y)图5.1 离散点分布示意图从图5.1观察数据点分布，用二次曲线拟合。
控制系统matlab仿真讲义
系统辨识与自适应控制MATLAB仿真第5章模型参考自适应控制
[ky]
k
(t)
(t)
(5-13)
10
由式（5-13）和上式有
(t) k(t) (t)
(5-14)
式中，k(t) k(t) k。 于是可以看到，当 k(t) k 时， (t) 0
现给出规范化的性能指标函数：
J
(k)
1 2
2 (t)
m2
式中，m 1 2(t) 为规范化信号。
按优化理论，k变化使 J (k) 极小的方向应按负梯度确定：
梯度设计法来叙述。
例5.1.3 设有被控对象式（5-1），仍采用参考模型式（5-2）
和控制器结构式（5-6）。设aˆp (t) 是未知对象参数 ap的估计值，式
（5-6）中的 k(t) 由下式计算：
k(t) aˆp (t) am
(5-16)
为了产生参数估计aˆp (t) ，选择一个稳定的滤波器
系统辨识与自适应控制 MATLAB仿真
第5章模型参考自适应控制
1
模型参考自适应控制是一种不同于自校正控制的另一类自适应 控制形式。根据被控对象结构和控制要求，设计参考模型，使其输 出表达对输入指令的期望响应，然后通过模型输出与被控对象输出 之差来调整控制器参数，使差值趋向于零，也就是使对象输出向模 型输出靠近，最终达到完全一致。根据控制器参数更新方法的不 同，模型参考自适应控制可分为直接自适应控制和间接自适应控制 两种。推演参数自适应规律的方法有两种：梯度法和稳定理论法。 5.1简单自适应控制系统
本节目的：1）给出直接自适应控制和间接自适应控制的概 念；2）自适应控制系统的两种基本设计方法：李亚普诺夫法和2梯
度法。
5.1.1直接自适应控制
直接自适应控制是指控制器参数直接从一个自适应规律中获取
控制系统数字仿真与CAD第五章习题答案
5-1 设控制系统的开环传递函数为2(1)()()(1)(416)K s G s H s s s s s +=-++试画出该系统的根轨迹。
解: 在Matlab 窗口中输入下列命令： num=[1 1]; a=[1 0]; b=[1 -1]; c=[1 4 16]; d=conv(a,b); den=conv(d,c); rlocus(num,den) grid on可得到系统的根轨迹如下图所示：5-2 某反馈控制系统的开环传递函数为()()()()24420KG s H s s s s s =+++ 试绘制其根轨迹。
解：在MATLAB 命令窗口中输入下列命令： num=1;den=conv(conv([1,0],[1,4]),[1,4,20]); rlocus(num,den) grid on 运行结果为：5-3.已知某系统传递函数为2180(1)100()11(1)[()20.31]40200200s W s s s s +=++⨯⨯+ 试绘制其伯德图。
解：分子分母同乘100*200得到280200(100)()(2.5100)(20.3200)200s W s s s s ⨯+=++⨯+在Matlab 窗口中输入下列命令：k=80*200; num=[1 100]; a=[ 100];b=[(1/200) 2* 200]; den=conv(a,b); w=logspace(-1,1,100); [m,p]=bode(k*num,den,w); subplot(2,1,1);semilogx(w,20*log10(m)); grid;xlabel('Frequency(rad/sec)'); ylabel('Gain(dB)'); subplot(2,1,2); semilogx(w,p); grid;xlabel('Frequency(rad/s)'); ylabel('Phase(deg)'); 可绘制该系统的伯德图如下所示。
matlab的控制系统仿真与应用 第五章
[ n a n , n1 a n1 ,,1 a1 ] T 1
a n 2 a1 a n 3 1 0 1 0 0
a n 1 a n2 T b A b A n 1 b a1 1 函数bass_pp( ) 调用格式为:
1
2 n K v1 v2 vn
n v1 v2 vn
1
i K ，于是 vi
vi V，对于给定 i ，可以求出 i i
vi
一般说来
v1
v2 v n 可逆，否则重新选择 1 2 n 。
n V1 1 V2 2 Vn n
ai
(i=1，2,…，n - 1)为特征多项式系数，即
sI a s n a1 s n 1 an 1 s an 0
定义新的状态
ˆ x
x=T
ˆ x
1．变换法
由于系统可控，因此T可逆，系统变换成
1 1 ˆ T aTx ˆ T bu x
其中
0 0 1 T aT 0 an 1 0 0 an 1 0 1 0 a n 2 0 0 0 0 1 T b 1 0 a1 1
若假定 i 与A的特征值有相同的，或 i 中有重根时， i 1,2,, n 则可以对特征值相同的一个或几个加上一定的微小偏量，使之满足上面第一种 情形的条件。然后，再重新进行极点配置。如果效果不够理想，那么还可重新选择 1 2 n 阵来进行配置。 控制系统工具箱中place( )函数是基于鲁棒极点配置的算法，用来求取状态反馈 阵K，使得多输入系统具有指定的闭环极点P，即 p eig( A B * K ) 。 place( )函数调用格式为：
控制系统仿真实验5：答案
实验5：利用MATLAB实现经典控制系统分析Array实验目的：1.练习并掌握求取多个模块串联、并联、反馈后的总的模型方法；2.练习复杂模型的连接方法。
实验内容：1.已知系统的结构图如图1所示，求整个系统的TF模型。
>> num1=[5,1];den1=[1,6,111];sys1=tf(num1,den1);>> num2=[15.6,29.32,1];den2=[12,26,37,102,1];sys2=tf(num2,den2);>> Sys_series=series(sys1,sys2);>> Sys_series1=sys1*sys2;>> tf(Sys_series1)>> tf(Sys_series)78 s^3 + 162.2 s^2 + 34.32 s + 1----------------------------------------------------------------12 s^6 + 98 s^5 + 1525 s^4 + 3210 s^3 + 4720 s^2 + 11328 s + 1112.系统如图2所示，求整个系统的TF模型和ZPK模型。
>> num1=[3,1];den1=[1,3,15];sys1=tf(num1,den1);>> z=[-2 -4];p=[0 -3 -5];k=20;sys2=zpk(z,p,k);>> num2=[0.1,1];den2=[2,1];sys3=tf(num2,den2);>> sys=sys1*sys2*sys3>> tf(sys)Zero/pole/gain:3 (s+0.3333) (s+2) (s+4) (s+10)--------------------------------------s (s+3) (s+5) (s+0.5) (s^2 + 3s + 15)Transfer function:3 s^4 + 49 s^3 + 220 s^2 + 308 s + 80---------------------------------------------------------s^6 + 11.5 s^5 + 59.5 s^4 + 192 s^3 + 307.5 s^2 + 112.5 s3.已知系统的结构图如图3所示，求整个系统的TF模型。
《控制系统计算机辅助设计MATLAB语言与应用第2版》薛定宇_课后习题答案
【17】
（1)z＝ｘy
＞>[x,y]=mesｈgrid(-3:0.01:3，-3：0.０1:３);
z=x.*y；
mesh(x,ｙ，z);
>> contoｕr3(x,ｙ,z,５0)；
(１)ｚ=sin（xｙ)
>> [x,y］=mesｈgrid(-3:０.01:3，-3:０.０1：3);
【２】
相应的ＭATＬAB命令:B=A(２：2:end,：）
>>Ａ=ｍagic(8)
A=
64２361６06７57
95５5412１3 51５01６
1７47４6 20２１43 42２4
40262737３6 30313３
3２３435 29２838 392５
41２32244４5１91８48
4９１514 5253 11 10 5６
【１0】
fｕｎcｔion y＝fｉb(ｋ)
if nargin~=1,ｅrror（＇出错：输入变量个数过多，输入变量个数只允许为１!')；ｅｎdﻭif narｇout＞１,ｅｒrｏr(＇出错：输出变量个数过多！'）;end
if k<=０,error('出错:输入序列应为正整数！＇）;endﻭiｆk==1|k＝=2,y＝１;ﻭｅｌsｅy=ｆib（k-１)＋fib(k-２);ｅnｄﻭeｎｄ
8５8 5954 62 6３1
＞＞Ｂ=Ａ(2:2:ｅnd，:)
B =
9５5 54１21３515０16
402６27３7 3６30３1３3
4１２3２24４4519１8 4８
858 5954 6２631
matlab语言第5章 Simulink仿真设计
信号源模块（Sources）
l Band-Limited White Noise：宽带限幅白噪声； l Chirp Signal：线性调频信号（频率随时间线性变 化的正弦波）；
l Clock：时钟信号； l Constant：常数输入； l Counter Free-Running：自动计数器，发生溢出后， 从0开始重新计数；
分； lPID Controller：PID控制； lPID Controller（2DOF）：二维PID控制； lState-Space：状态空间模型； lTransfer-Fcn：传递函数模型； lTransport Delay：输入信号延迟一个固定时间输出； lVariable Time Delay：输入可变时间信号延迟输出； lVariable Transport Delay：输入信号延迟可变输出； lZero-Pole：零极点模型。
3. Simulink模型的创建和仿真
以图示系统建立Simulink模型
模型建立
（1）启动工具箱 （2）建立Simulink空白模型 （3）根据系统模型选择模块 首先要确定所需模块所在的子模块库名称。例子中用 到的模块有单位阶跃信号﹑符号比较器﹑传递函数模 型和信号输出模块，分别属于信号源模块库﹑数学运 算模块库﹑连续模块库和输出模块库。在模块库浏览 器中打开相应的模块库，选择所需模块。
常用模块库（Commonly Used Blocks）
l Bus Creator：将输入信号合并为总线信号； l Bus Selector：由总线信号选择需要的信号输出； l Constant：常数信号； l Data Type Conversion：数据类型转换模块； l Delay：延迟模块； l Delux：信号分解模块； l Discrete-Time Integrator：离散时间积分器； l Gain：增益模块； l Ground：接地模块； l In1：输入模块； l Integrator：输入信号积分；
MATLAB语言与控制系统仿真-参考答案-第5章
控制系统的时域响应MATLAB 仿真实训1. 学会利用MATLAB 绘制系统的单位阶跃响应曲线，掌握读取系统动态性能指标的方法；2. 学会利用MATLAB 绘制系统的单位脉冲响应曲线的方法；3. 掌握利用MATLAB 绘制系统的零输入响应曲线的方法；4. 掌握利用MATLAB 绘制系统的一般输入响应曲线的方法；5. 学会通过仿真曲线读取相关信息，并依据有关信息进行系统的时域分析。
1.编写程序求取下列各系统的单位阶跃响应，完成表5-5并记录相关曲线。
162.316)(21++=s s s G 164.216)(22++=s s s G 166.116)(23++=s s s G 1616)(24++=s s s G 解：>> n1=16;>> d1=[1,3.2,16]; >> sys1=tf(n1,d1); >> step(sys1)>> n2=16;>> d2=[1,2.4,16]; >> sys2=tf(n2,d2); >> step(sys2)>> n3=16;>> d3=[1,1.6,16]; >> sys3=tf(n3,d3); >> step(sys3) >> n4=16;>> d4=[1,1,16]; >> sys4=tf(n4,d4); >> step(sys4)序号ξnωm axcptst(%5=∆) 计算值实验计算值实验计算值实验值1 42 4 1.37233 4 1.52664 4 1.6731 w=4;cmax1=1+exp(-z1*pi/sqrt(1-z1^2));tp1=pi/(w*sqrt(1-z1^2));ts1=3.5/(z1*w);[cmax1,tp1,ts1]ans =>> z2=0.3;w=4;cmax2=1+exp(-z2*pi/sqrt(1-z2^2));tp2=pi/(w*sqrt(1-z2^2));ts2=3.5/(z2*w);[cmax2,tp2,ts2]ans =>> z3=0.2; w=4;cmax3=1+exp(-z3*pi/sqrt(1-z3^2)); tp3=pi/(w*sqrt(1-z3^2)); ts3=3.5/(z3*w); [cmax3,tp3,ts3]ans =>> z4=0.125; w=4;cmax4=1+exp(-z4*pi/sqrt(1-z4^2)); tp4=pi/(w*sqrt(1-z4^2)); ts4=3.5/(z4*w); [cmax4,tp4,ts4] ans =说明：对于二阶欠阻尼系统（10<<ξ），若系统的闭环传递函数为2222)(nn ns s s Φωξωω++= 则系统单位阶跃响应的输出最大值21max 1ξξπ--+=ec峰值时间21ξωπ-=n p t调整时间估算值ns t ξω5.3= （以5%为误差带） ns t ξω4.4=（以2%为误差带）2.已知二阶系统的闭环传递函数如下，编程求取系统的单位阶跃响应并完成表5-6，记录相关曲线。
自动控制原理的MATLAB仿真与实践第5章 线性系统的频域分析
7
【例5-1】 试绘制惯性环节G(jω)=1/(2s+1)的Nyquist曲线 和Bode图。
解：程序如下：
>>clear
G=tf(1,[2,1]); %建立模型
nyquist (G); %绘制Nyquist图
figure(2); bode (G); %绘制Bode图
4
ngrid；ngrid(‘new’)：绘制尼科尔斯坐标网格即等 20lgM圆和等角曲线组成的网格。‘new’代表清除以前 的图形，与后一个nichols()一起绘制网格。
semilogx(w,20*log10(mag))：绘制半对数坐标下的幅 频特性曲线。
semilogx(w,phase*180/pi)：绘制半对数坐标下的相频 特性曲线。
MATLAB提供了许多用于线性系统频率分析 的函数命令，可用于系统频域的响应曲线、参数分析 和系统设计等。常用的频率特性函数命令格式及其功 能见表5-1。 bode (G)：绘制传递函数的伯德图。其中：G为传递
函数模型，如：tf(), zpk(), ss()。 bode(num,den)：num,den分别为传递函数的分子与
本节分别介绍利用MATLAB进行频域绘图和频 率分析的基本方法。
6
5.2.1 Nyquist曲线和Bode图
MATLAB频率特性包括幅频特性和相频特性。 当用极坐标图描述系统的幅相频特性时，通常称为 奈奎斯特（Nyquist）曲线；用半对数坐标描述系 统的幅频特性和相频特性时，称为伯德（Bode） 图；在对数幅值-相角坐标系上绘制等闭环参数（ M和N）轨迹图，称为尼克尔斯（Nichols）图。
运行结果如图5-2所示。
matlab的控制系统仿真与应用第五章
x Ac x Bc v x(t0 ) x0 , t t0 y Cc x
v G
Ac A B K , Bc B G, Cc C
x ∫
+ A K
+
-
u
B
+
x
C
y
状态反馈极点配置：通过状态反馈矩阵K的选取，使闭环系统的极点，即 ( A B K ) 的特征值 i (i 1,2,, n) 恰好处于所希望的一组给定闭环极点的位置上。 线性定常系统可以用状态反馈任意配置极点的充分必要条件是：该系统必须是完全 能控的。所以，在实现极点的任意配置之前，必须判别受控系统的能控性。
i K ，于是 vi
vi V，对于给定 i ，可以求出 i i
vi
K 1 2 n1 V2 2 Vn n
1 n1
1 疋田算法的具体步骤：首先，适当选择 i R r，从而计算特征向量
5.2.1单输入系统的极点配置
Bass-Gura算法：设受控系统的闭环特征多项式分别为：
(s) det(sI n A) s n a1s n1 an1s an
(s) (s 1 )(s 2 )(s n ) s n 1s n1 n1s n
, n) ，表示闭环系统的极点及其相对应的特征向量。
假定 i 与A阵的特征值相异，且 i j (i j)i, j 1, 2, , n
有 则
(i I n A B K ) vi 0即
x A x B u u K x
(i I n A) vi B K vi
x
∫
控制系统的MATLAB仿真与设计课后答案
第二章1>>x=[15 22 33 94 85 77 60]>>x(6)>>x([1 3 5])>>x(4:end)>>x(find(x>70))2>>T=[1 -2 3 -4 2 -3] ;>>n=length(T);>>TT=T';>>for k=n-1:-1:0>>B(:, n-k)=TT.%>>end>>B>>test=vander(T)3>>A=zeros(2,5);>>A(:)=-4:5>>L=abs(A)>3>>islogical(L)>>X=A(L)4>>A=[4,15,-45,10,6;56,0,17,-45,0]>>find(A>=10&A<=20)5>>p1=conv([1,0,2],conv([1,4],[1,1]));>>p2=[1 0 1 1];>>[q,r]=deconv(p1,p2);>>cq=' 商多项式为' ; cr= ' 余多项式为' ;>>disp([cq,poly2str(q, 's' )]),disp([cr,poly2str(r,'s' )]) 6>>A=[11 12 13;14 15 16;17 18 19];>>PA=poly(A)>>PPA=poly2str(PA, 's' )第三章1>>n=(-10:10)';>>y=abs(n);>>plot(n,y, 'r.' , 'MarkerSize' ,20)>>axis equal>>grid on>>xlabel( 'n' ) 2>>x=0:pi/100:2*pi;>>y=2*exp(-0.5*x).*sin(2*pi*x); >>plot(x,y),grid on;3>>t=0:pi/50:2*pi;>>x=8*cos(t);>>y=4*sqrt(2)*sin(t);>>z=-4*sqrt(2)*sin(t);>>plot3(x,y,z, 'p' );>>title( 'Line in 3-D Space' );>>text(0,0,0, 'origin' );>>xlabel( 'X' ),ylable( 'Y' ),zlable( 'Z' );grid;4>>theta=0:0.01:2*pi;>>rho=sin(2*theta).*cos(2*theta);>>polar(theta,rho, 'k' );5>>[x,y,z]=sphere(20);>>z1=z;>>z1(:,1:4)=NaN;>>c1=ones(size(z1));>>surf(3*x,3*y,3*z1,c1);>>hold on>>z2=z;>>c2=2*ones(size(z2));>>c2(:,1:4)=3*ones(size(c2(:,1:4)));>>surf(1.5*x,1.5*y,1.5*z2,c2);>>colormap([0,1,0;0.5,0,0;1,0,0]);>>grid on>>hold off第四章1>>for m=100:999m1=fix(m/100);m2=rem(fix(m/10),10);m3=rem(m,10);if m==m1*m1*m1+m2*m2*m2+m3*m3*m3disp(m) endend2 M文件：function [s,p]=fcircle(r)s=pi*r*r;p=2*pi*r; 主程序：[s,p]=fcircle(10) 3>>y=0;n=100; for i=1:ny=y+1/i/i;end>>y4 M 文件：function f=factor(n)if n<=1f=1;elsef=factor(n-1)*n;end 主程序：>>s=0;for i=1:5s=s+factor(i);end>>s5>>sum=0;i=1;while sum<2000sum=sum+i;i=i+1;end;>>n=i-26 for循环M文件：function k=jcsum(n)k=0;for i=0:nk=k+2A i;end 主程序：>>jcsum(63)While循环M文件：function k=jcsum1(n)k=0;i=0;while i<=nk=k+2Ai;i=i+1;end 主程序：>>jcsum1(63) 第五章1>>A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4]; >>b=[13,-9,6,0]';>>x=A\b2 M文件：function f=fxyz(u)x=u(1);y=u(2);z=u(3);f=x+y42./x/4+z42./y+2./z;主程序：[U,fmi n]=fmi nsearch( 'fxyz' ,[0.5,0.5,0.5])3>>X=linspace(0,2*pi,50);>>Y=sin(X);>>P=polyfit(X,Y,3)>>AX=linspace(0,2*pi,50);>>Y=sin(X);>>Y1=polyval(P,X)>>plot(X,Y, ':O' ,X,Y1, '-*' )4>>x=0:2.5:10;>>h=[0:30:60]'; >>T=[95,14,0,0,0;88,48,32,12,6;67,64,54, 48,41];>>xi=[0:0.5:10];>>hi=[0:10:60]';>>temps=interp2(x,h,T,xi,hi,'cubic' ); >>mesh(xi,hi,temps);第六章1>>syms x>>y=finverse(1/tan(x))2>>syms x y>>f=1/(1+x A2); g=s in (y);>>fg=compose(f,g)3>>syms x>>g=(exp(x)+x*sin(x))A(1/2);>>dg=diff(g)4>>F=int(int( 'x*exp(-x*y)' , 'x' ), 'y' )5>>syms x>>F=ztrans(x*exp(-x*10))6>>a=[0 1;-2 -3];>>syms s>>inv(s*eye(2)-a);7>>f=solve( 'a*xA2+b*x+c' )8>>f=solve( 'x+y+z=1' , 'x-y+z=2' , '2*x-y-z=1' )9>>y=dsolve( 'D2y+2*Dy+2*y=0' ,'y(0)=1' , 'Dy(0)=0' ) >>ezplot(y),grid on10>>a=maple('simplify(sin(x)A2+cos(x)A2);' )1 1 >>f=maple( 'laplace(exp(-3*t)*sin(t),t,s);' )12>>syms t x>>F=sin(x*t+2*t); >>L=laplace(F) 第七章13第八章1-1>>h=tf([5,0],[1,2,2])1-2>>s = tf( 's' );>>H = [5/(sH+2*s+2)];>>H.inputdelay =21-3>>h=tf([0.5,0],[1,-0.5,0.5],0.1)2>>num=2*[1,0.5];den=[1,0.2,1.01]; >>sy s=tf(num,den)>>[z,p,k]=tf2zp(num,den);>>zpk(z,p,k)>>[A,B,C,D]=tf2ss(num,den);>>ss(A,B,C,D)3 >>num=[1,5]; den=[1,6,5,1];ts=0.1; >>sysc=tf(num,den);>>sysd=c2d(sysc,ts, 'tustin' )4>>r1=1;r2=2;c1=3;c2=4;>>[A,B,C,D]=linmod( 'x84' );>>[num,den]=ss2tf(A,B,C,D);>>sys=tf(num,den)5>>A=[1,1,0;0,1,0;0,0,2];B=[0,0;1,0;0,-2]; >>n=size(A)>>Tc=ctrb(A,B);if n==rank(Tc)disp( ' 系统完全能控' );elsedisp( ' 系统不完全能控' );end第九章1>>num=[2,5,1];den=[1,2,3];>>bode(num,den);grid on ;>>figure;>>nyquist(num,den);2>>num=5*[1,5,6];den=[1,6,10,8];>>step(num,den);grid on ;>>figure;>>impulse(num,den);grid on ;3>>kosi=0.7;wn=6;>>n um=w n^2;de n=[1,2*kosi*w n,wn A2]; >>step(num,den);grid on ;>>figure;>>impulse(num,den);grid on ;4 M 文件：function [rtab,info]=routh(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),a3(i)=rtab(k,i+1)-alpha(k)*rtab(k+1,i+1);endif sum(abs(a3))==0a3=polyder(vec2);info=[info, 'All elements in row ' , ... int2str(k+2) ' are zeros;' ];elseif abs(a3(1))<epsa3(1)=1e-6;info=[info, 'Replaced first element;' ];endrtab=[rtab; a3, zeros(1,nrT-length(a3))];vec1=vec2; vec2=a3;end主程序：>>den=[1,2,8,12,20,16,16];>>[rtab,info]=routh(den)>>a=rtab(:,1)if all(a>0)disp( ' 系统是稳定的' );elsedisp( ' 系统是不稳定的' );end5>>num=7*[1,5];den=conv([1,0,0],conv([1,10],[1,1]));>>[gm,pm,wg,wc]=margin(num,den)第十章1 M 文件：function varargout=rg_lead(ng0,dg0,s1)if nargout==1ngv=polyval(ng0,s1);dgv=polyval(dg0,s1);g=ngv/dgv;thetal=pi-angle(g);zc=real(s1)-imag(s1)/tan(thetal);t=-1/zc;varargout{1}=[t,1];elseif nargout==2ngv=polyval(ng0,s1); dgv=polyval(dg0,s1);g=ngv/dgv; theta=angle(g); phi=angle(s1);if theta>0phi_c=pi-theta;endif theta<0;phi_c=-thetaendtheta_z=(phi+phi_c)/2; theta_p=(phi-phi_c)/2; z_c=real(s1)-imag(s1)/tan(theta_z); p_c=real(s1)-imag(s1)/tan(theta_p);nk=[1 -z_c]; varargout{2}=[1 -p_c]; kc=abs(p_c/z_c);if theta<0kc=-kcendvarargout{1}=kc*nk;elseerror( ' 输出变量数目不正确！' );end主程序：>> ng0=[1]; dg0=10000*[1 0 -1.1772];>>g0=tf(ng0,dg0); %满足开环增益的为校正系统的传递函数>>s=kw2s(0.7,0.5) %期望的闭环主导极点>>ngc=rg_lead(ng0,dg0,s);>>gc=tf(ngc,1)>>g0c=tf(g0*gc);>>rlocus(g0,g0c);>>b1=feedback(g0,1); %未校正系统的闭环传递函数>>b2=feedback(g0c,1); %校正后系统的闭环传递函数>>figure,step(b1, 'r--' ,b2, 'b' ); grid on %绘2 M文件：function [ngc,dgc,k]=rg_lag(ng0,dg0,KK,s1,a)ngv=polyval(ng0,s1);dgv=polyval(dg0,s1);g=dgv/ngv;k=abs(g); %期望主导极点处的根轨迹增益beta=k/KK;[kosi1,wn1]=s2kw(s1);zc=-w n1*si n(a*pi/180)/si n(pi-ata n(sqrt(1-kosi1A2)/kosi1)-(a*pi/180)); %利用正弦定理pc=beta*zc;ngc=beta*[1,-zc];dgc=[1,-pc];主程序：>>KK=20;s1=-2+i*sqrt(6);a=1 >>ng0=[10];dg0=conv([1,0],[1,4]);>>g0=tf(ng0,dg0);>>[ngc,dgc,k]=rg_lag(ng0,dg0,KK,s1,a); >>gc=tf(ngc,dgc)>>g0c=tf(KK*g0*gc);>>b1=feedback(k*g0,1);>>b2=feedback(g0c,1);>>step(b1, 'r--' ,b2, 'b' );grid on3 M文件：function [ngc,dgc,k]=rg_lag(ng0,dg0,KK,s1,a)ngv=polyval(ng0,s1);dgv=polyval(dg0,s1); g=dgv/ngv;k=abs(g); %期望主导极点处的根轨迹增益beta=k/KK;[kosi1,wn1]=s2kw(s1);zc=-w n1*si n(a*pi/180)/si n(pi-ata n(sqrt(1-kosi1A2)/kosi1)-(a*pi/180)); %利用正弦定理pc=beta*zc; ngc=beta*[1,-zc];dgc=[1,-pc];主程序：>>KK=128;s1=-2+i*2*sqrt(3);a=2 >>ng0=[10];dg0=conv([1,0],conv([1,2],[1,8]));>>g0=tf(ng0,dg0);>>[ngc,dgc,k]=rg_lag(ng0,dg0,KK,s1,a);>>gc=tf(ngc,dgc)>>g0c=tf(KK*g0*gc);>>rlocus(g0,g0c);>>b1=feedback(k*g0,1);>>b2=feedback(g0c,1);>>figure,step(b1, 'r--' ,b2, 'b' );grid on4 M文件：function [ngc,dgc]=lead4(ng0,dg0,KK,Pm,w)[mu,pu]=bode(KK*ng0,dg0,w); [gm,pm,wcg,wcp]=margin(mu,pu,w); alf=ceil(Pm-pm+5); phi=(alf)*pi/180;a=(1+sin(phi))/(1-sin(phi)), dbmu=20*log10(mu);mm=-10*log10(a); wgc=spline(dbmu,w,mm), T=1/(wgc*sqrt(a)), ngc=[a*T,1]; dgc=[T,1]; 主程序：>>ng0=[1];dg0=conv([1,0,0],[1,5]);>>g0=tf(ng0,dg0);>>w=logspace(-3,3);>>KK=1;Pm=50;>>[ngc,dgc]=lead4(ng0,dg0,KK,Pm,w);>>gc=tf(ngc,dgc);g0c=tf(KK*g0*gc);>>bode(KK*g0,w);hold on,bode(g0c,w);grid on,hold off>>[gm,pm,wcg,wcp]=margin(g0c) >>Kg=20*log10(gm) >>g1=feedback(g0c,1);>>bode(g1),gridon, >>[mag,phase,w]=bode(g1); >>a=find(mag<=0.707*mag(1)); >>wb=w(a(1)) >>max(mag) >>b=find(mag==max(mag)) >>wr=w(b) 5 M 文件： function [mu,pu]=bode(ng0,dg0,w);[gm,pm,wcg,wcp]=margin(mu,pu,w); 率alf=ceil(Pm-pm+5);度， phi=(alf)*pi/180; a=(1+sin(phi))/(1-sin(phi)); dbmu=20*log10(mu); mm=-10*log10(a); wgc=spline(dbmu,w,mm); T=1/(wgc*sqrt(a)); ngc=[a*T,1]; dgc=[T,1]; 主程序： >>KK=40; Pm=50; >>ng0= KK *[1];dg0=conv([1,0],conv([1,1],[1,4])); >>g0=tf(ng0,dg0) ;>>w=logspace(-2,4); >>[ngc,dgc]=fg_lead_pm(ng0,dg0,Pm,w) >>gc=tf(ngc,dgc),g0c=tf(g0*gc);>>b1=feedback(g0,1); b2=feedback(g0c,1); >>step(b1, 'r--' , b2, 'b' ); grid on >>figure, bode(g0, 'r--' ,g0c, 'b' ,w),grid >>[gm,pm,wcg,wcp]=margin(g0c), Km=20*log10(gm)[ngc,dgc]=fg_lead_pm(ng0,dg0,Pm,w)%计算原系统的对数频率响应数据%求取原系统的相角裕度和剪切频%计算控制器提供的最大超前角%将最大超前角转换为弧度单位灿算a 值%系统的对数幅值%w 处的控制器对数幅值%差值求取wm 认为wm= wc%+算T on,。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
控制系统的时域响应MATLAB 仿真实训实训目的 1. 学会利用MATLAB 绘制系统的单位阶跃响应曲线，掌握读取系统动态性能指标的方法； 2. 学会利用MATLAB 绘制系统的单位脉冲响应曲线的方法； 3. 掌握利用MATLAB 绘制系统的零输入响应曲线的方法； 4. 掌握利用MATLAB 绘制系统的一般输入响应曲线的方法；5.学会通过仿真曲线读取相关信息，并依据有关信息进行系统的时域分析。
实训内容1.编写程序求取下列各系统的单位阶跃响应，完成表5-5并记录相关曲线。
162.316)(21++=s s s G 164.216)(22++=s s s G 166.116)(23++=s s s G 1616)(24++=s s s G 解：>> n1=16;>> d1=[1,,16]; >> sys1=tf(n1,d1); >> step(sys1)>> n2=16;>> d2=[1,,16]; >> sys2=tf(n2,d2); >> step(sys2)>> n3=16;>> d3=[1,,16]; >> sys3=tf(n3,d3); >> step(sys3)>> n4=16;>> d4=[1,1,16]; >> sys4=tf(n4,d4); >> step(sys4)序号ξnωm axcptst(%5=∆)计算值实验计算值实验计算值实验值14243444w=4;cmax1=1+exp(-z1*pi/sqrt(1-z1^2)); tp1=pi/(w*sqrt(1-z1^2));ts1=(z1*w);[cmax1,tp1,ts1]ans =>> z2=;w=4;cmax2=1+exp(-z2*pi/sqrt(1-z2^2)); tp2=pi/(w*sqrt(1-z2^2));ts2=(z2*w);[cmax2,tp2,ts2]ans =>> z3=; w=4;cmax3=1+exp(-z3*pi/sqrt(1-z3^2)); tp3=pi/(w*sqrt(1-z3^2)); ts3=(z3*w); [cmax3,tp3,ts3]ans =>> z4=; w=4;cmax4=1+exp(-z4*pi/sqrt(1-z4^2)); tp4=pi/(w*sqrt(1-z4^2)); ts4=(z4*w); [cmax4,tp4,ts4] ans =说明：对于二阶欠阻尼系统（10<<ξ），若系统的闭环传递函数为2222)(nn ns s s Φωξωω++= 则系统单位阶跃响应的输出最大值21max 1ξξπ--+=ec峰值时间21ξωπ-=n p t调整时间估算值ns t ξω5.3= （以5%为误差带） ns t ξω4.4=（以2%为误差带）2.已知二阶系统的闭环传递函数如下，编程求取系统的单位阶跃响应并完成表5-6，记录相关曲线。
（1）828)(21++=s s s Φ （2）10210)(22++=s s s Φ（3）12212)(23++=s s s Φ（4）16216)(24++=s s s Φ>> n1=8;>> d1=[1,2,8]; >> sys1=tf(n1,d1); >> impulse(sys1)>> n2=10;>> d2=[1,2,10]; >> sys2=tf(n2,d2); >> impulse(sys2)>> n3=12;>> d3=[1,2,12]; >> sys3=tf(n3,d3); >> impulse(sys3)>> n4=16;>> d4=[1,2,16]; >> sys4=tf(n4,d4); >> impulse(sys4)表5-6序号 ξn ωm ax c 实验p t 实验1 2 3 443.已知某单位负反馈系统的开环传递函数为)11.0)(1(5)(++=s s s s G若系统的输入信号分别为（1）t u =1， )20,0(∈t （2）tteu 5.02-=， )20,0(∈t（3）te t u 5.03)2sin(-=，)20,0(∈t（4）te t u 34)2cos(-=，)20,0(∈t编写程序分别求取系统的在给定的输入信号下的响应，记录相应的曲线。
>> n=5;>> d=conv([1,1,0],[,1]); >> [nc,dc]=cloop(n,d) >>sysc=tf(nc,dc); >> t=0::20; >> u1=t;>> u2=t.*exp*t);>> u3=sin(2*t).*exp*t); >> u4=cos(2*t).*exp(-3*t); >>subplot(221);lsim(sysc,u1,t);subplot(222);lsim(sysc,u2,t);subplot(223);lsim(sys c,u3,t);subplot(224);lsim(sysc,u4,t);4．编写程序绘制以下系统的单位阶跃响应曲线和单位脉冲响应曲线。
（1）u x x⎥⎦⎤⎢⎣⎡+⎥⎦⎤⎢⎣⎡--=114001&[]u x y 65.30+=>> a1=[-1,0;0,-4];>> b1=[1;1]; >> c1=[0,]; >> d1=6;>> sys1=ss(a1,b1,c1,d1);>> subplot(1,2,1);step(sys1);subplot(1,2,2);impulse(sys1);（2）u x x ⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡+⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡-=733412219553641251308342& x y ⎥⎦⎤⎢⎣⎡=05132301>> a2=[2,4,3,8;0,3,1,5;2,1,4,6;3,5,-5,9]; >> b2=[1,2;2,1;4,3;3,7]; >> c2=[1,0,3,2;3,1,5,0]; >> d2=0;>> sys2=ss(a2,b2,c2,d2); >> step(sys2)>> impulse(sys2)（3）)53)(4)(5()8.13.6)(2(5)(23++++++=s s s s s s s G >> sysa=zpk([-2],[-5;-4],[5]) Zero/pole/gain: 5 (s+2) ----------- (s+5) (s+4)>> sysb=tf([,],[1,3,0,5]) Transfer function: s +--------------- s^3 + 3 s^2 + 5>> sys3=sysa*sysbZero/pole/gain:(s+2) (s+-------------------------------------------------------(s+5) (s+4) (s+ (s^2 - +>> step(sys3)>> impulse(sys3)5.有三个系统如图5-17（a）、（b）、（c）所示，编程分别绘制各系统的单位阶跃响应和单位脉冲响应。
（a）>> sysa1=tf([5],[1,3,5]);>> sysa2=zpk([-2],[0;-3;-5],[20]);>> sysa3=tf([,1],[2,1]);>> sysa=sysa1*sysa2*sysa3Zero/pole/gain:5 (s+2) (s+10)------------------------------------------------s (s+3) (s+5) (s+ (s^2 + 3s + 5)>> step(sysa)>> impulse(sysa)（b）>> sysb1=tf([22,1],[1,5,23,15]);>> sysb2=zpk([-6],[0;-3;-10],[]);>> sysb3=5*tf([,1],[2,1])*tf([,1],[3,7,1]);>> sysb=sysb1-sysb2+sysb3Zero/pole/gain:(s+ (s+ (s+ (s+ (s+ (s^2 + +--------------------------------------------------------------------------------------------------------------------s (s+3) (s+ (s+ (s+ (s+10) (s+ (s^2 + +>> step(sysb)>> impulse(sysb)(c)>> G1=tf([10],[1,2,0]);>> G2=cloop(tf([1],[,1]),1)Transfer function:1------------s + 1>> G2=feedback(tf([1],[,1]),1,-1) Transfer function:1------------s + 2>> G=G1*G2Transfer function:10---------------------------s^3 + 3 s^2 + 4 s>> H=tf([3],[1,6,0])+tf([1],[,1])Transfer function:s^2 + s + 3-----------------------------s^3 + s^2 + 6 s>> sys=G/(1+G*H)Transfer function:s^6 + 11 s^5 + 82 s^4 + 244 s^3 + 240 s^2------------------------------------------------------------------------------------------------------------------s^9 + s^8 + s^7 + s^6 + 123 s^5 + s^4 + 340 s^3 + 342 s^2 + 120 s>> step(sys)>> impulse(sys)6.设三阶系统闭环传递函数为81066)52()()(232+++++==Φs s s s s s R s C s ）（试绘制其单位阶跃响应在]3.5,2.1[∈t 内的部分响应曲线。