《控制系统计算机辅助设计：MATLAB语言与应用(第2版)》薛定宇_课后习题答案.doc 68页

'第1章控制系统计算机辅助设计概述第2章MATLAB语言程序设计基础第3章线性控制系统的数学模型第4章线性控制系统的计算机辅助分析第5章Simulink在系统仿真中的应用第6章控制系统计算机辅助设计68 第1章控制系统计算机辅助设计概述【1】http://www.mathworks.cn/已阅，略【2】已阅，略【3】已经掌握help命令和Help菜单的使用方法【4】区别：MATLAB语言实现矩阵的运算非常简单迅速，且效率很高，而用其他通用语言则不然，很多通用语言所实现的矩阵运算都是对矩阵维数具有一点限制的，即使限制稍小的，但凡维数过大，就会造成运算上的溢出出错或者运算出错，甚至无法处理数据的负面结果68 【5】【8】(1)输入激励为正弦信号68 (2)输入激励为脉冲模拟信号(3)输入激励为时钟信号68 (4)输入激励为随机信号(5)输入激励为阶跃信号δ=0.368 δ=0.05δ=0.7结论：随着非线性环节的死区增大，阶跃响应曲线的范围逐渐被压缩，可以想象当死区δ足够大时，将不再会有任何响应产生。所以可以得到结论，在该非线性系统中，死区的大小可以改变阶跃响应的幅值和超调量。死区越大，幅值、超调量将越小，而调整时间几乎不受其影响68 第2章MATLAB语言程序设计基础【1】>>A=[1234;4321;2341;3241]A=1234432123413241>>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.0000i3.0000+2.0000i4.0000+1.0000i4.0000+1.0000i3.0000+2.0000i2.0000+3.0000i1.0000+4.0000i2.0000+3.0000i3.0000+2.0000i4.0000+1.0000i1.0000+4.0000i3.0000+2.0000i2.0000+3.0000i4.0000+1.0000i1.0000+4.0000i>>A(5,6)=5A=123400432100234100324100000005∴若给出命令A(5,6)=5则矩阵A的第5行6列将会赋值为5，且其余空出部分均补上0作为新的矩阵A，此时其阶数为5×6【2】相应的MATLAB命令：B=A(2:2:end,:)>>A=magic(8)A=64236160675795554121351501617474620214342244026273736303133323435292838392541232244451918484915145253111056858595462631>>B=A(2:2:end,:)B=95554121351501640262737363031334123224445191848858595462631∴从上面的运行结果可以看出，该命令的结果是正确的68 【3】>>symsxs;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=18446744073709551615【5】>>fori=1:120if(i==1|i==2)a(i)=1;elsea(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,14930352,24157817,39088169,63245986,102334155,165580141,267914296,433494437,701408733,1134903170,1836311903,2971215073,4807526976,7778742049,12586269025,20365011074,32951280099,53316291173,86267571272,139583862445,225851433717,365435296162,591286729879,956722026041,1548008755920,2504730781961,4052739537881,6557470319842,10610209857723,17167680177565,27777890035288,44945570212853,72723460248141,117669030460994,190392490709135,308061521170129,498454011879264,806515533049393,1304969544928657,2111485077978050,3416454622906707,5527939700884757,8944394323791464,14472334024676221,23416728348467685,37889062373143906,61305790721611591,99194853094755497,160500643816367088,259695496911122585,420196140727489673,679891637638612258,1100087778366101931,1779979416004714189,2880067194370816120,4660046610375530309,7540113804746346429,12200160415121876738,19740274219868223167,31940434634990099905,51680708854858323072,83621143489848422977,135301852344706746049,218922995834555169026,354224848179261915075,573147844013817084101,927372692193078999176,1500520536206896083277,2427893228399975082453,3928413764606871165730,6356306993006846248183,10284720757613717413913,16641027750620563662096,26925748508234281076009,43566776258854844738105,70492524767089125814114,114059301025943970552219,184551825793033096366333,298611126818977066918552,483162952612010163284885,781774079430987230203437,1264937032042997393488322,2046711111473984623691759,3311648143516982017180081,5358359254990966640871840]68 【6】>>k=1;fori=2:1000forj=2:iifrem(i,j)==0ifjD)+(h.*x/D).*(abs(x)<=D)-h.*(x<-D)68 【10】functiony=fib(k) ifnargin~=1,error("出错：输入变量个数过多，输入变量个数只允许为1！");end ifnargout>1,error("出错：输出变量个数过多！");end ifk<=0,error("出错：输入序列应为正整数！");end ifk==1|k==2,y=1; elsey=fib(k-1)+fib(k-2);end end【13】68 【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);gridon;【15】(1)>>t=-2*pi:0.01:2*pi;r=1.0013*t.^2;polar(t,r);axis("square")68 (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")68 【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);68 (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);68 第3章线性控制系统的数学模型【1】(1)>>s=tf("s");G=(s^2+5*s+6)/(((s+1)^2+1)*(s+2)*(s+4))Transferfunction:s^2+5s+6--------------------------------s^4+8s^3+22s^2+28s+16(2)>>z=tf("z",0.1);H=5*(z-0.2)^2/(z*(z-0.4)*(z-1)*(z-0.9)+0.6)Transferfunction:5z^2-2z+0.2---------------------------------------z^4-2.3z^3+1.66z^2-0.36z+0.6Samplingtime(seconds):0.1【2】(1)该方程的数学模型>>num=[6422];den=[1103232];G=tf(num,den)Transferfunction:6s^3+4s^2+2s+2------------------------s^3+10s^2+32s+32(2)该模型的零极点模型>>G=zpk(G)Zero/pole/gain:6(s+0.7839)(s^2-0.1172s+0.4252)-------------------------------------(s+4)^2(s+2)(3)由微分方程模型可以直接写出系统的传递函数模型68 【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)Samplingtime(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)Transferfunction:7.58s^2+10.8s+4.8--------------------------------------------------------------0.7896s^5+4.183s^4+7.811s^3+13.81s^2+12.61s+4.8(2)状态方程的标准型实现>>G1=ss(G1)a=x1x2x3x4x5x1-5.297-2.473-2.186-0.9981-0.7598x240000x302000x400200x50000.50b=u1x12x20x30x40x5068 c=x1x2x3x4x5y1000.60.42730.3799d=u1y10Continuous-timestate-spacemodel.(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))Transferfunction:3s^6+6s^5+3s^4+42s^3+84s^2+42s---------------------------------------------------------------------------s^10+3s^9+55s^8+175s^7+300s^6+1323s^5+2656s^4+3715s^3+7732s^2+5602s+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)68 【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))Transferfunction:0.004873s^5+1.036s^4+61.15s^3+649.7s^2+1911s---------------------------------------------------------------------------4.357e-014s^10+2.422e-011s^9+5.376e-009s^8+6.188e-007s^7+4.008e-005s^6+0.001496s^5+0.03256s^4+0.4191s^3+2.859s^2+8.408s68 第4章线性控制系统的计算机辅助分析【1】(1)>>num=[1];den=[3212];G=tf(num,den);eig(G)ans=-1.00000.1667+0.7993i0.1667-0.7993i分析：由以上信息可知，系统的极点有2个是在s域的右半平面的，因此系统是不稳定的(2)>>num=[1];den=[63211];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=[11-3-12];G=tf(num,den);eig(G)ans=-2.0000-1.00001.00001.0000分析：由以上信息可知，系统的极点有2个是在s域的右半平面的，因此系统是不稳定的(4)>>num=[31];den=[3006005031];G=tf(num,den);eig(G)ans=-1.9152-0.14140.0283+0.1073i0.0283-0.1073i分析：由以上信息可知，系统的极点有2个是在s域的右半平面的，因此系统是不稳定的68 (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=[-32];den=[1-0.2-0.250.05];H=tf(num,den,"Ts",0.5);abs(eig(H)")ans=0.50000.50000.2000分析：由以上信息可知，所有特征根的模均小于1，因此该系统是稳定的(2)>>num=[3-0.39-0.09];den=[1-1.71.040.2680.024];H=tf(num,den,"Ts",0.5);abs(eig(H)")ans=1.19391.19390.12980.1298分析：由以上信息可知，由于前两个特征根的模均大于1，因此该系统是不稳定的(3)>>num=[13-0.13];den=[11.3520.44810.0153-0.01109-0.001043];H=tf(num,den,"Ts",0.5);abs(eig(H)")ans=0.87430.15200.27230.23440.1230分析：由以上信息可知，所有特征根的模均小于1，因此该系统是稳定的(4)>>num=[2.1211.7615.91];den=[1-7.368-20.15102.480.39-340];H=tf(num,den,"Ts",0.5,"Variable","q");abs((eig(H))")ans=8.23493.21152.34152.34322.3432分析：由以上信息可知，所有特征根的模均大于1，因此该系统是不稳定的68 【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,…21.896,3;11,18.0898,25,2.356,9];abs(eig(F)")ans=63.720723.539312.436613.323119.7275分析：由以上信息可知，该离散系统的F矩阵的所有特征根的模均大于1，因此该系统是不稳定的【4】>>A=[-3121;0-4-2-1;12-11;-1-11-2];B=[10;02;03;11];C=[122-1;21-12];D=[00;00];G=ss(A,B,C,D);tzero(G)pzmap(G)ans=-3.6124-1.2765结论：∴可以得到该系统的零点为-3.6124、-1.2765分析：由以上信息可知，系统的特征根的实部均位于s域的左半平面，因此该系统是稳定的68 【5】>>s=tf("s");G=0.2*(s+2)/(s*(s+0.5)*(s+0.8)*(s+3)+0.2*(s+2));Gc=sscanform(G,"ctrl")Go=sscanform(G,"obsv")a=x1x2x3x4x10100x20010x30001x4-0.4-1.4-4.3-4.3b=u1x10x20x30x41c=x1x2x3x4y10.40.200d=u1y10Continuous-timestate-spacemodel.a=x1x2x3x4x1000-0.4x2100-1.4x3010-4.3x4001-4.3b=u1x10.4x20.2x30x40c=x1x2x3x4y10001d=u1y10Continuous-timestate-spacemodel.68 【9】(1)>>num=[1851459823638012266422208818576040320];den=[1365464536224496728411812410958440320];[R1,P1,K1]=residue(num,[den0]);[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.00000>>[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=[1851459823638012266422208818576040320];den=[1365464536224496728411812410958440320];[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-t68 (3)>>symst;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=[1851459823638012266422208818576040320];den=[1365464536224496728411812410958440320];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.0010i0.0000+3.0000i-1.1973+0.0010i0.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][正弦信号时域响应的解析解]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=[1851459823638012266422208818576040320];den=[1365464536224496728411812410958440320];G=tf(num,den);t=[1:.1:20]";u=sin(3*t+5);lsim(G,u,t);分析：由解析解可知，输出信号的稳态部分是振荡的，并且其幅值与相位始终在到达稳态的时候保持不变，因此右图所示的输出波形与解析解所得的结论是一致的68 【10】(1)因为PI或PID控制器均含有Ki/s项，这是一个对误差信号的积分环节，假设去掉这一环节，则当Kp→∞，即|e(t)|很小也会存在较大扰动，这会影响到系统的动态特性；当加入这一环节后，如果要求|e(t)|→0，则控制器输出u(t)会由Ki/s环节得到一个常值，此时系统可以获得较好的动态特性，因此这两个控制器可以消除闭环系统的阶跃响应的稳态误差(2)不稳定系统能用PI或PID控制器消除稳态误差。因为PI或PID控制器均含有积分控制(I)，在积分控制中，控制器的输出与输入误差信号的积分成正比关系。对一个自动控制系统，如果在进入稳态后存在稳态误差，则称这个控制系统是有稳态误差的或简称有差系统。为了消除稳态误差，在控制器中必须引入“积分项”。积分项对误差取决于时间的积分，随着时间的增加，积分项会增大。这样，即便误差很小，积分项也会随着时间的增加而加大，它推动控制器的输出增大使稳态误差进一步减小，直到等于零。因此，比例+积分(PI)控制器，可以使系统在进入稳态后无稳态误差，即不稳定系统能用PI或PID控制器消除稳态误差【13】(1)>>P=[0;-3;-4+4i;-4-4i];Z=[-6;6];G=zpk(Z,P,1);rlocus(G),grid分析：根据根轨迹图可知，可知无论K取何值，均无法保证所有极点均在s域左半平面，因此使单位负反馈系统稳定的K值范围是不存在的68 (2)>>num=[122];den=[111480];G=tf(num,den);rlocus(G),grid分析：根据根轨迹图可知，使单位负反馈系统稳定的K值范围为0>num=[1];den=[1/26001/2610];G=tf(num,den);rlocus(G),grid分析：根据根轨迹图可知，使单位负反馈系统稳定的K值范围为0>s=tf("s");G=800*(s+1)/(s^2*(s+10)*(s^2+10*s+50));rlocus(G),grid分析：根据根轨迹图可知，使单位负反馈系统稳定的K值范围为0>P=[0,-10,-20,-40];Z=[];G=zpk(Z,P,1);rlocus(G),grid分析：从右图的根轨迹图可知，存在能够使得闭环系统主导极点大约ξ=0.707阻尼比的K值，且其位于最右边的那段根轨迹上，下面将利用等阻尼线来寻找符合要求的点68 结论：由上图可知，能使得闭环系统主导极点有大约ξ=0.707阻尼比的K值为K≈22400【20】(1)>>s=tf("s");G=8*(s+1)/(s^2*(s+15)*(s^2+6*s+10));margin(G)>>nyquist(G);68 >>eig(G)ans=00-15.0000-3.0000+1.0000i-3.0000-1.0000i由上述开环系统极点分布可知，开环系统是稳定的；再结合Nyquist图，可以发现Nyquist图不包围(-1,j0)点，因此闭环系统是稳定的>>step(feedback(G,1))闭环系统的阶跃响应结论：频域分析与阶跃响应所得的结论一致68 (2)>>s=tf("s");G=4*(s/3+1)/(s*(0.02*s+1)*(0.05*s+1)*(0.1*s+1));margin(G)>>nyquist(G)68 >>eig(G)ans=0-50.0000-20.0000-10.0000由上述开环系统极点分布可知，开环系统是稳定的；再结合Nyquist图，可以发现Nyquist图不包围(-1,j0)点，因此闭环系统是稳定的>>step(feedback(G,1))闭环系统的阶跃响应结论：频域分析与阶跃响应所得的结论一致68 (3)>>A=[021;-3-20;134];B=[4;3;2];C=[123];D=[0];G=ss(A,B,C,D);margin(G)>>nyquist(G),grid68 >>eig(G)ans=-0.9463+1.8633i-0.9463-1.8633i3.8926由上述开环系统极点分布可知，因为系统中有一个极点位于s域的右半平面，故该开环系统是不稳定的；再结合Nyquist图，可以发现Nyquist图不包围(-1,j0)点，因此闭环系统是不稳定的>>step(feedback(G,1))闭环系统的阶跃响应结论：频域分析与阶跃响应所得的结论一致68 (4)>>P=[0;1;0.368;0.99];Z=[-1.31;-0.054;0.957];G=zpk(Z,P,0.45,"Ts",0.1);margin(G)>>nyquist(G),grid68 >>eig(G)ans=01.00000.36800.9900由上述开环系统极点分布可知，因为除了一个点位于单位圆上，其他所有极点均位于单位圆内，故该开环系统是不稳定的；再结合Nyquist图，可以发现Nyquist图不包围(-1,j0)点，因此闭环系统是不稳定的>>step(feedback(G,1))闭环系统的阶跃响应结论：频域分析与阶跃响应所得的结论一致68 (5)>>s=tf("s");G=6*(-s+4)/(s^2*(0.5*s+1)*(0.1*s+1));margin(G)>>nyquist(G),grid68 >>eig(G)ans=00-10.0000-2.0000由上述开环系统极点分布可知，开环系统是稳定的；再结合Nyquist图，可以发现Nyquist图顺时针包围(-1,j0)点1圈，因此闭环系统是不稳定的>>step(feedback(G,1))闭环系统的阶跃响应结论：频域分析与阶跃响应所得的结论一致68 (6)>>num=[10-6011060];den=[11782130100];G=tf(num,den);margin(G)>>nyquist(G);grid68 >>eig(G)ans=-10.0000-5.0000-1.0000+1.0000i-1.0000-1.0000i由上述开环系统极点分布可知，开环系统是稳定的；再结合Nyquist图，可以发现Nyquist图顺时针包围(-1,j0)点2圈，因此闭环系统是不稳定的>>step(feedback(G,1))闭环系统的阶跃响应结论：频域分析与阶跃响应所得的结论一致68 第5章Simulink在系统仿真中的应用【1】1.输入模块组Sources 2.输出池模块组sbfSinks 3.连续系统模块组Continuous 4.离散系统模块组Discrete 5.非线性模块组Discontinuities 6.数学函数模块组MathOperations 7.查表模组块Look-upTables 8.用户自定义函数模块组User-definedFunctions 9.信号模块组SignalRouting 10.信号属性模块组SignalAttributes【2】Ø方法一：运用Simulink搭建系统的仿真模型进行分析原微分方程y(4)(t)+5y(3)(t)+63y’’(t)+4y’(t)+2y(t)=e-3t+e-5tsin(4t+π/3)令x1(t)=y(t)，x2(t)=y’(t)，x3(t)=y’’(t)，x4(t)=y(3)(t)，则原线性微分方程可化简为x1’(t)=x2(t)x2’(t)=x3(t)x3’(t)=x4(t)x4’(t)=-2x1(t)-4x2(t)-63x3(t)-5x4(t)+e-3t+e-5tsin(4t+π/3)根据初值条件y(0)=1，y’(0)=y’’(0)=1/2，y(3)(t)=0.2即可得到如下的simulink仿真模型68 (1)从simulink仿真模型的示波器Scope上观察到的输出波形y(t)(2)通过MATLAB命令得到的输出波形y(t)>>[t,x,y]=sim("simulink_Chapter5_2",250);plot(t,y);grid68 Ø方法二：运用微分方程数值解进行分析第一步：首先编写一个MATLAB函数来描述该微分方程，将原微分方程化成“方法一”中的四变量形式后即可编写出如下代码，将其在工作目录下保存为vdp_eq5_2.m文件%——————————————————————————————————————functiony=vdp_eq5_2(t,x) y=[x(2);x(3);x(4);-2*x(1)-4*x(2)-63*x(3)-5*x(4)+exp(-3*t)+exp(-5*t)*sin(4*t+pi/3)]; end%——————————————————————————————————————第二步：在MATLAB的工作区CommandWindow编写如下代码来求出微分方程的数值解，并将数值解[t,y]用图形直观地表示出来，即该图像即为输出波形y(t)%——————————————————————————————————————>>x0=[1;1/2;1/2;0.2];[t,y]=ode45("vdp_eq5_2",[0,250],x0);plot(t,y(:,1));grid%——————————————————————————————————————Ø分析结果从以上仿真来看，无论是采用“方法一”还是“方法二”，均能够得到非常理想的仿真结果曲线。因此对于该微分方程而言，这两种方法对于数值解的求解均是不错的选择68 【3】Ø方法一：运用Simulink搭建系统的仿真模型进行分析原微分方程y(4)(t)+5ty(3)(t)+6t2y’’(t)+4y’(t)+2e-2ty(t)=e-3t+e-5tsin(4t+π/3)令x1(t)=y(t)，x2(t)=y’(t)，x3(t)=y’’(t)，x4(t)=y(3)(t)，则原线性微分方程可化简为x1’(t)=x2(t)x2’(t)=x3(t)x3’(t)=x4(t)x4’(t)=-2e-2tx1(t)-4x2(t)-6t2x3(t)-5tx4(t)+e-3t+e-5tsin(4t+π/3)根据初值条件y(0)=1，y’(0)=y’’(0)=1/2，y(3)(t)=0.2即可得到如下的simulink仿真模型68 (1)从simulink仿真模型的示波器Scope上观察到的输出波形y(t)(2)通过MATLAB命令得到的输出波形y(t)>>[t,x,y]=sim("simulink_Chapter5_3",10);plot(t,y);grid68 Ø方法二：运用微分方程数值解进行分析第一步：首先编写一个MATLAB函数来描述该微分方程，将原微分方程化成“方法一”中的四变量形式后即可编写出如下代码，将其在工作目录下保存为vdp_eq5_3.m文件%——————————————————————————————————————functiony=vdp_eq5_3(t,x) y=[x(2);x(3);x(4);-2*exp(-2*t)*x(1)-4*x(2)-6*t^2*x(3)-5*t*x(4)+exp(-3*t)+exp(-5*t)*…sin(4*t+pi/3)]; end%——————————————————————————————————————第二步：在MATLAB的工作区CommandWindow编写如下代码来求出微分方程的数值解，并将数值解[t,y]用图形直观地表示出来，即该图像即为输出波形y(t)%——————————————————————————————————————>>x0=[1;1/2;1/2;0.2];[t,y]=ode45("vdp_eq5_3",[0,10],x0);plot(t,y(:,1));grid%——————————————————————————————————————Ø分析结果从以上仿真来看，无论是采用“方法一”还是“方法二”，均能够得到非常理想的仿真结果曲线。因此对于该微分方程而言，这两种方法对于数值解的求解均是不错的选择68 【6】(1)从示波器Scope观察所得波形(2)通过在MATLAB工作目录下输入命令得到波形>>[t,x,y]=sim("simulink_Chapter5_6",50);plot(t,y(:,1),"-",t,y(:,2),"--*");grid%“实线”的为输出1，“虚线且转折带星号”的为输出268 【7】上述仿真模型中对于输入u采用了Kp=[1,0.7,0.3,0.2;0.6,1,0.4,0.35;0.35,0.4,1,0.6;0.2,0.3,0.7,1]对系统进行补偿，用以减小四路输出的耦合效果(1)从示波器Scope观察所得波形(2)通过在MATLAB工作目录下输入命令得到波形>>u1=1;u2=1;u3=1;u4=1;[tt,x,yy]=sim("simulink_Chapter5_7",40);subplot(221),plot(tt,yy(:,1));gridsubplot(222),plot(tt,yy(:,2));gridsubplot(223),plot(tt,yy(:,3));gridsubplot(224),plot(tt,yy(:,4));grid68 Ø通过step()函数对上述模型进行仿真(1)方案一，8路输出，MATLAB命令如下：>>g11=tf(1,[4,1]);g12=tf(0.7,[5,1]);g13=tf(0.3,[51]);g14=tf(0.2,[51]);g21=tf(0.6,[5,1]);g22=tf(1,[4,1]);g23=tf(0.4,[51]);g24=tf(0.35,[51]);g31=tf(0.35,[5,1]);g32=tf(0.4,[5,1]);g33=tf(1,[4,1]);g34=tf(0.6,[51]);g41=tf(0.2,[5,1]);g42=tf(0.3,[5,1]);g43=tf(0.7,[51]);g44=tf(1,[4,1]);G1=[g11g12g13g14;g21g22g23g24;g31g32g33g34;g41g42g43g44];Kp=[1,0.7,0.3,0.2;0.6,1,0.4,0.35;0.35,0.4,1,0.6;0.2,0.3,0.7,1];G2=ss(G1*Kp);[y1,x1,t1]=step(G2.a,G2.b,G2.c,G2.d,1,40);[y2,x2,t2]=step(G2.a,G2.b,G2.c,G2.d,2,40);[y3,x3,t3]=step(G2.a,G2.b,G2.c,G2.d,3,40);[y4,x4,t4]=step(G2.a,G2.b,G2.c,G2.d,4,40);u1=1;u2=0;u3=0;u4=0;[tt1,x1,yy1]=sim("simulink_Chapter5_7",40);u1=0;u2=1;u3=0;u4=0;[tt2,x2,yy2]=sim("simulink_Chapter5_7",40);u1=0;u2=0;u3=1;u4=0;[tt3,x3,yy3]=sim("simulink_Chapter5_7",40);u1=0;u2=0;u3=0;u4=1;[tt4,x4,yy4]=sim("simulink_Chapter5_7",40);subplot(4,4,1),plot(t1,y1(:,1),":",tt1,yy1(:,1));subplot(4,4,2),plot(t1,y1(:,2),":",tt1,yy1(:,2));subplot(4,4,3),plot(t1,y1(:,3),":",tt1,yy1(:,3));subplot(4,4,4),plot(t1,y1(:,4),":",tt1,yy1(:,4));subplot(4,4,5),plot(t2,y2(:,1),":",tt2,yy2(:,1));subplot(4,4,6),plot(t2,y2(:,2),":",tt2,yy2(:,2));subplot(4,4,7),plot(t2,y2(:,3),":",tt2,yy2(:,3));subplot(4,4,8),plot(t2,y2(:,4),":",tt2,yy2(:,4));subplot(4,4,9),plot(t3,y3(:,1),":",tt3,yy3(:,1));subplot(4,4,10),plot(t3,y3(:,2),":",tt3,yy3(:,2));subplot(4,4,11),plot(t3,y3(:,3),":",tt3,yy3(:,3));subplot(4,4,12),plot(t3,y3(:,4),":",tt3,yy3(:,4));subplot(4,4,13),plot(t4,y4(:,1),":",tt4,yy4(:,1));subplot(4,4,14),plot(t4,y4(:,2),":",tt4,yy4(:,2));68 subplot(4,4,15),plot(t4,y4(:,3),":",tt4,yy4(:,3));subplot(4,4,16),plot(t4,y4(:,4),":",tt4,yy4(:,4));%注：仿真结果中，实线表示精确仿真结果，虚线表示近似结果Ø分析：从图中可以看出，采用step()函数进行仿真的精度还是非常高的，在所得图形中几乎看不出二者的差别(2)方案二，8路输出，MATLAB命令如下>>step(G2)68 (3)最终的4路输出波形，MATLAB命令如下>>[y1,x1,t1]=step(G2.a,G2.b,G2.c,G2.d,1,40);[y2,x2,t2]=step(G2.a,G2.b,G2.c,G2.d,2,40);[y3,x3,t3]=step(G2.a,G2.b,G2.c,G2.d,3,40);[y4,x4,t4]=step(G2.a,G2.b,G2.c,G2.d,4,40);subplot(221);plot(t1,y1(:,1)+y1(:,2)+y1(:,3)+y1(:,4));subplot(222);plot(t2,y2(:,1)+y2(:,2)+y2(:,3)+y2(:,4));subplot(223);plot(t3,y3(:,1)+y3(:,2)+y3(:,3)+y3(:,4));subplot(224);plot(t4,y4(:,1)+y4(:,2)+y4(:,3)+y4(:,4));Ø总结通过对于该4输入4输出多变量系统的Simulink仿真和采用step()函数进行仿真，可以发现得到以下结论两种情况下仿真结果有以下的异同点：(1)异：1、Simulink的仿真结果是最终输出y1和y4均在达到5时趋于平稳，而用step()函数仿真的结果是y1和y4在到达5还不到的位置就趋于平稳了2、Simulink的仿真结果是最终输出y2和y3均在达到4.7附近时趋于平稳，而用step()函数仿真的结果是y2和y3在到4.8附近的位置就趋于平稳，两者之间的差别也很小，但是可以很容易发现这个微小的差异(2)同：两种情况下均满足y1和y4、y2和y3得到的输出波形是对应相同的，这就是为什么在Simulink仿真的示波器上只能看到两个波形的原因，这是由于两对相同形状的波形叠加在一起所造成的最终结论：两种情况下的仿真结果差异不大，在实际工程建模上都是适用的68 【11】x1’=tsin(x2e-2.3x4)x2’=x1x3’=sin(x2e-2.3x4)x4’=x3【13】(1)将直流电机拖动系统的模型在Simulink平台上搭建出来，提取系统的总模型(2)绘制系统的阶跃响应、频域响应曲线Ø阶跃响应>>[t,x,y]=sim("simulink_Chapter5_13",0.8);plot(t,y);grid68 Ø频域响应>>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));figure(1);bode(G);figure(2);nyquist(G);figure(3);rlocus(G);68 68 第6章控制系统计算机辅助设计【1】>>G=zpk([-1.5],[-1.75,-16,-1.5+3i,-1.5-3i],210);Gc=zpk([-1.5],[-14.86],52.5);f1=figure;f2=figure;f3=figure;f4=figure;figure(f1);margin(G);figure(f2);margin(G*Gc);figure(f3);margin(G);holdon;margin(G*Gc)figure(f4);step(G);holdon;step(G*Gc);68 原系统的Bode图校正后系统的Bode图调整前后比较：从Gc(s)的表达式不难发现，这是一个超前校正器。通过比较校正前后系统的Bode图，可以发现校正后的系统的相位裕度反而减小了，此时系统稳定性明显减小了，但此时时间响应速度却反而增大了校正后系统的阶跃响应原系统的阶跃响应分析：通过观察阶跃响应曲线，可以发现校正后的系统较原系统幅值明显增大，瞬态响应速度加快，但是缺点是其超调量变得过大，调整时间增加建议：设计控制器应当充分考虑系统的稳定性，然后再考虑准确性和快速性68 【2】ØG(s)=16/(s(s+1)(s+2)(s+8))<1>原系统的特性>>G1=zpk([],[0;-1;-2;-8],16);subplot(121);margin(G1);subplot(122);step(feedback(G1,1));分析：可以发现原系统所构成的闭环系统超调量过大，这对于系统稳定性具有非常大的影响<2>设计超前-滞后校正器，改进系统动态特性，设所期望的剪切频率为wc=1，通过选择不同的期望相位裕度值来选择合适的相位裕度%第一阶段筛选合适的期望相位裕度值>>wc=1;forgam=20:10:90Gc=leadlagc(G1,wc,gam,1000,3);figure(1);margin(Gc*G1);holdon;figure(2),step(feedback(Gc*G1,1),70);holdon;end68 r增加方向分析：通过上图比较可以发现，随着相位裕度值的增加，闭环系统的超调量是始终在减小的，此时我们需要考虑的一个关键因素就是调整时间了，兼顾这两个因素，只有当这两者达到了最优，最终得到的系统才是理想的。可以发现，当相位裕度值在50~60这个范围的时候其构成的闭环系统的动态特性达到了最好的效果，和原系统的动态特性相比，这个范围的优势在于其调整时间的增加并不明显，且超调量较原系统已经明显减小。而再往下，调整时间又开始变长了，因此我们将在这个范围内搜素合适的相位裕度值%第二阶段筛选合适的期望相位裕度值>>wc=1;forgam=50:5:60Gc=leadlagc(G1,wc,gam,1000,3);step(feedback(Gc*G1,1),40);holdon;end68 分析：相比之下，对于相位裕度值50之后的闭环系统的调整时间增大的并不明显，而60之后的开始有调整时间增大的趋势，55对于50而已其超调量又近一步减小了，因此我们最终选择相位裕度值gam=55%设计：超前-滞后校正器的最终指标定为剪切频率wc=1，相位裕度值gam=55>>wc=1;gam=55;Gc=leadlagc(G1,wc,gam,1000,3)step(feedback(G1*Gc,1));Transferfunction:3.726s^2+1.966s+0.1593----------------------------s^2+2.338s+0.0003726分析：此时闭环系统的超调量明显减小，且调整时间较原闭环系统并不明显，因此此时超前-滞后校正器的设计是合理的68 ØG(s)=2(s+1)/(s(47.5s+1)(0.0625s+2)^2)<1>原系统的特性>>s=tf("s");G2=2*(s+1)/(s*(47.5*s+1)*(0.0625*s+1)^2);subplot(121);margin(G2);subplot(122);step(feedback(G2,1));分析：可以发现原系统所构成的闭环系统超调量略大，振荡频率高，且调整时间过长，这对系统的运作是非常不利的<2>设计超前-滞后校正器，改进系统动态特性，设所期望的剪切频率为wc=25，通过选择不同的期望相位裕度值来选择合适的相位裕度值%筛选合适的期望相位裕度值>>wc=25;forgam=20:10:50Gc=leadlagc(G2,wc,gam,1000,3);figure(1);margin(Gc*G2);holdon;figure(2),step(feedback(Gc*G2,1));holdon;end68 r增加方向分析：通过上图比较可以发现，系统此时的振荡次数明显减小，调整时间已经从原来的将近200s缩短到了现在的6s，因此这种设计非常有效地减小了调整时间；此时我们只需找到在这个条件下合适的超调量足够小的响应曲线，即能确定出符合要求的期望相位裕度值；可以发现一点，如果我们将超调量那部分减小为负值时，则该闭环系统的振荡就变成仅仅一次了，且调整时间未发生明显变化，因此我们可以选择使调整量变为负值的那个响应曲线所对应的相位裕度值，通过比较，我们最终确定gam=45为期望的相位裕度值%设计：超前-滞后校正器的最终指标定为剪切频率wc=25，相位裕度值gam=45>>wc=25;gam=45;Gc=leadlagc(G2,wc,gam,1000,3)step(feedback(G2*Gc,1));Transferfunction:1.289e004s^2+8.328e004s+1.276e005---------------------------------------s^2+168.1s+1612分析：此时闭环系统的超调量已经消除，振荡已经减少到了一次，且调整时间也由原来的200s缩短到了6s，响应速度变快，且动态特性也达到了最优，因此这个设计是合理的68 【8】(1)Ga(s)=1/(s+1)^3利用“整定公式”设计控制器STEP1:求出受控对象的一阶近似模型分析：从得出的拟合结果可以看出，对于这个式子，采用次最优降阶方法得出的拟合结果最接近，因此就选择该方法所得的表达式STEP2:根据整定公式设计PID控制器根据由该受控对象得出的较好的FOLPD，可知K=1；T=1.997；L=1.13结论：所以由此设计的PID控制器模型为Gc(s)=2.1207(1+1/(2.2600s)+0.5650s/(0.0565s+1))STEP3:分析给出的受控对象模型在控制器下的阶跃响应曲线68 分析：可见最终响应的超调量小，振荡少，调整时间快，达到了最终的设计目的利用PID控制器设计控制器根据由该受控对象得出的较好的FOLPD，可知K=1；T=1.997；L=1.13，因此可以设计如下PID控制器68 分析：在该PID控制器作用下，可以发现最终响应的超调量小，振荡少，调整时间快，达到了最终的设计目的(2)Gb(s)=1/(s+1)^5利用“整定公式”设计控制器STEP1:求出受控对象的一阶近似模型68 分析：从得出的拟合结果可以看出，对于这个式子，采用次最优降阶方法得出的拟合结果最接近，因此就选择该方法所得的表达式STEP2:根据整定公式设计PID控制器根据由该受控对象得出的较好的FOLPD，可知K=1；T=2.624；L=2.59结论：所以由此设计的PID控制器模型为Gc(s)=1.2158(1+1/(5.1800s)+1.2950s/(0.12950s+1))STEP3:分析给出的受控对象模型在控制器下的阶跃响应曲线分析：可见最终响应的超调量小，振荡较少，调整时间快，达到了最终的设计目的利用PID控制器设计控制器68 根据由该受控对象得出的较好的FOLPD，可知K=1；T=2.624；L=2.59，因此可以设计如下PID控制器分析：在该PID控制器作用下，可以发现最终响应的超调量小，振荡较少，调整时间快，达到了最终的设计目的(3)Gb(s)=(-1.5s+1)/(s+1)^3利用“整定公式”设计控制器STEP1:求出受控对象的一阶近似模型68 分析：从得出的拟合结果可以看出，对于这个式子，采用次最优降阶方法得出的拟合结果最接近，因此就选择该方法所得的表达式STEP2:根据整定公式设计PID控制器根据由该受控对象得出的较好的FOLPD，可知K=1；T=1.769；L=2.52结论：所以由此设计的PID控制器模型为Gc(s)=0.8424(1+1/(5.0400s)+1.2600s/(0.12600s+1))STEP3:分析给出的受控对象模型在控制器下的阶跃响应曲线68 分析：可见最终效果并不是特别明显，阶跃响应的超调量显然已经减小不少，但是其振荡还是过多，调整时间还算比较快的，达到了最终的设计目的，总的来说还是不错的利用PID控制器设计控制器根据由该受控对象得出的较好的FOLPD，可知K=1；T=1.769；L=2.52，因此可以设计如下PID控制器68 分析：在该PID控制器作用下，可以发现最终效果并不是特别明显，阶跃响应的超调量显然已经减小不少，但是其振荡还是过多，调整时间还算比较快的，达到了最终的设计目的，总的来说还是不错的68'