控制系统计算机辅助设计：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 4
4 3 2 1
2 3 4 1
3 2
4 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.0000i
2.0000 +
3.0000i 3.0000 + 2.0000i
4.0000 + 1.0000i
4.0000 + 1.0000i 3.0000 + 2.0000i 2.0000 + 3.0000i 1.0000 + 4.0000i
2.0000 +
3.0000i 3.0000 + 2.0000i
4.0000 + 1.0000i 1.0000 + 4.0000i
3.0000 + 2.0000i 2.0000 + 3.0000i
4.0000 + 1.0000i 1.0000 + 4.0000i
>> A(5,6)=5
A =
1 2 3 4 0 0
4 3 2 1 0 0
2 3 4 1 0 0
3 2
4 1 0 0
0 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 57
9 55 54 12 13 51 50 16
17 47 46 20 21 43 42 24
40 26 27 37 36 30 31 33
32 34 35 29 28 38 39 25
41 23 22 44 45 19 18 48
49 15 14 52 53 11 10 56
8 58 59 5 4 62 63 1
>> B=A(2:2:end,:)
B =
9 55 54 12 13 51 50 16
40 26 27 37 36 30 31 33
41 23 22 44 45 19 18 48
8 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+6
f =
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)^5
m =
simplify(100)
【4】
>> i=0:63; s=sum(2.^sym(i))
s =
0615
【5】
>> for i=1:120
if(i==1|i==2) a(i)=1;
else a(i)=a(i-1)+a(i-2);end
if(i==120) a=sym(a); disp(a); end
end
[ 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:1000
for j=2:i
if rem(i,j)==0
if j<i, break;end
if j==i, A(k)=i; k=k+1; break; end
end
end
end
disp(A);
Columns 1 through 13
2 3 5 7 11 13 17 19 23 29 31 37 41 Columns 14 through 26
43 47 53 59 61 67 71 73 79 83 89 97 101 Columns 27 through 39
103 107 109 113 127 131 137 139 149 151 157 163 167 Columns 40 through 52
173 179 181 191 193 197 199 211 223 227 229 233 239 Columns 53 through 65
241 251 257 263 269 271 277 281 283 293 307 311 313 Columns 66 through 78
317 331 337 347 349 353 359 367 373 379 383 389 397 Columns 79 through 91
401 409 419 421 431 433 439 443 449 457 461 463 467 Columns 92 through 104
479 487 491 499 503 509 521 523 541 547 557 563 569 Columns 105 through 117
571 577 587 593 599 601 607 613 617 619 631 641 643 Columns 118 through 130
647 653 659 661 673 677 683 691 701 709 719 727 733 Columns 131 through 143
739 743 751 757 761 769 773 787 797 809 811 821 823 Columns 144 through 156
827 829 839 853 857 859 863 877 881 883 887 907 911 Columns 157 through 168
919 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！');end
if nargout>1,error('出错：输出变量个数过多！');end
if k<=0,error('出错：输入序列应为正整数！');end
if k==1|k==2,y=1;
else y=fib(k-1)+fib(k-2);end
end
【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.6
Sampling 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.7598
x2 4 0 0 0 0 x3 0 2 0 0 0 x4 0 0 2 0 0
x5 0 0 0 0.5 0
b =
u1
x1 2
x2 0
x3 0
x4 0
x5 0
c =
x1 x2 x3 x4 x5
y1 0 0 0.6 0.4273 0.3799
d =
u1
y1 0
Continuous-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.0000
0.1667 + 0.7993i
0.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.4356i
0.2449 + 0.5688i
0.2449 - 0.5688i
分析：由以上信息可知，系统的极点有2个是在s域的右半平面的，因此系统是不稳定的(3) >> num=[1];den=[1 1 -3 -1 2];
G=tf(num,den);
eig(G)
ans =
-2.0000
-1.0000
1.0000
1.0000
分析：由以上信息可知，系统的极点有2个是在s域的右半平面的，因此系统是不稳定的(4) >> num=[3 1];den=[300 600 50 3 1];
G=tf(num,den);
eig(G)
ans =
-1.9152
-0.1414
0.0283 + 0.1073i
0.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.50
x 10-6
P 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 0
Continuous-time state-space model. a =
x1 x2 x3 x4
x1 0 0 0 -0.4
x2 1 0 0 -1.4
x3 0 1 0 -4.3
x4 0 0 1 -4.3
b =
u1
x1 0.4
x2 0.2
x3 0
x4 0
c =
x1 x2 x3 x4
y1 0 0 0 1
d =
u1
y1 0
Continuous-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.0000
0.2000 -6.0000
0.7361 -5.0000
-2.8889 -4.0000
2.2250 -
3.0000
-2.0222 -2.0000
3.0004 -1.0000
1.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.0000
7.3306 -7.0000
-1.2000 -6.0000
-3.6806 -5.0000
11.5556 -4.0000
-6.6750 -3.0000
4.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.0000
0.9559 -7.0000
-0.1761 -6.0000
-0.6111 -5.0000
2.1663 -4.0000
-1.1973 - 0.0010i 0.0000 + 3.0000i
-1.1973 + 0.0010i 0.0000 - 3.0000i
-1.3824 -3.0000
0.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 在积分控制中，控制器的输出与输入误差信号的积分成正比关系。

对一个自动控制系统，如果在进入稳态后存在稳态误差，则称这个控制系统是有稳态误差的或简称有差系统。

为了消除稳态误差，在控制器中必须引入“积分项”。

积分项对误差取决于时间的积分，随着时间的增加，积分项会增大。

这样，即便误差很小，积分项也会随着时间的增加而加大，它推动控制器的输出增大使稳态误差进一步减小，直到等于零。

因此，比例+积分(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 值范围是不存在的
(2) >> num=[1 2 2];den=[1 1 14 8 0];
G=tf(num,den);
rlocus(G),grid
分析：根据根轨迹图可知，使单位负反馈系统稳定的K 值范围为0<K<5.42
(3) >> num=[1];den=[1/2600 1/26 1 0];
G=tf(num,den);
rlocus(G),grid
分析：根据根轨迹图可知，使单位负反馈系统稳定的K 值范围为0<K<101
(4) >> s=tf('s');
G=800*(s+1)/(s^2*(s+10)*(s^2+10*s+50));
rlocus(G),grid
分析：根据根轨迹图可知，使单位负反馈系统稳定的K 值范围为0<K<3.23
00.20.4
0.60.8
11.21.41.6
1.82Step Response
Time (seconds)A m p l i t u d e
>> P=[0,-10,-20,-40];Z=[];
G=zpk(Z,P ,1);
rlocus(G),grid
分析：从右图的根轨迹图可知，存在能够使得闭环系统主导极点大约ξ=0.707阻尼比的K 值，且其位于最右边的那段根轨迹上，下面将利用等阻尼线来寻找符合要求的点
结论：由上图可知，能使得闭环系统主导极点有大约ξ=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);
>> eig(G)
ans =
-15.0000
-3.0000 + 1.0000i
-3.0000 - 1.0000i
由上述开环系统极点分布可知，开环系统是稳定的；再结合Nyquist 图，可以发现Nyquist 图不包围(-1,j0)点，因此闭环系统是稳定的 (2) >> s=tf('s'); G=4*(s/ margin(G)
>> nyquist(G)
>> eig(G) ans =
0 -50.0000 -20.0000
-10.0000 图不包围(-1,j0)闭环系统的阶跃响应 结论：频域分析与阶跃响应所得的结论一致
(3) >> A=[0 2 1;-3 -2 0;1 3 4];B=[4;3;2];C=[1 2 3];D=[0];
G=ss(A,B,C,D);
margin(G)
>> nyquist(G),grid
>> eig(G)
-0.9463 + 1.8633i
-0.9463 - 1.8633i
3.8926
由上述开环系统极点分布可知，因为系统中有一个极点位于s域的右半平面，故该开环系统是不稳定的；再结合Nyquist图，可以发现Nyquist图不包围(-1,j0)点，因此闭环系统是不稳定的
>> step(feedback(G,1))
闭环系统的阶跃响应
结论：频域分析与阶跃响应所得的结论一致
(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),grid
>> eig(G)
ans =
1.0000
0.3680
0.9900
由上述开环系统极点分布可知，因为除了一个点位于单位圆上，其他所有极点均位于单位圆内，故该开环系统是不稳定的；再结合Nyquist图，可以发现Nyquist图不包围(-1,j0)点，因此闭环系统是不稳定的
>> step(feedback(G,1))
闭环系统的阶跃响应
结论：频域分析与阶跃响应所得的结论一致
(5) >> s=tf('s');
G=6*(-s+4)/(s^2*(0.5*s+1)*(0.1*s+1));
margin(G)
>> nyquist(G),grid
>> eig(G)
ans =
-10.0000
-2.0000
由上述开环系统极点分布可知，开环系统是稳定的；再结合Nyquist图，可以发现Nyquist 图顺时针包围(-1,j0)点1圈，因此闭环系统是不稳定的
>> step(feedback(G,1))
闭环系统的阶跃响应
结论：频域分析与阶跃响应所得的结论一致
(6) >> num=[10 -60 110 60];den=[1 17 82 130 100];
G=tf(num,den);
margin(G)
>> nyquist(G);grid
>> eig(G)
ans =
-10.0000
-5.0000
-1.0000 + 1.0000i
-1.0000 - 1.0000i
由上述开环系统极点分布可知，开环系统是稳定的；再结合Nyquist图，可以发现Nyquist 图顺时针包围(-1,j0)点2圈，因此闭环系统是不稳定的
>> step(feedback(G,1))
闭环系统的阶跃响应
结论：频域分析与阶跃响应所得的结论一致
第5章 Simulink在系统仿真中的应用
【1】
1.输入模块组Sources
2.输出池模块组sbf Sinks
3.连续系统模块组Continuous
4.离散系统模块组Discrete
5.非线性模块组Discontinuities
6.数学函数模块组Math Operations
7.查表模组块Look-up Tables
8.用户自定义函数模块组User-defined Functions
9.信号模块组Signal Routing
10.信号属性模块组Signal Attributes
【2】
➢方法一：运用Simulink搭建系统的仿真模型进行分析
原微分方程y(4)(t)+5y(3)(t)+63y’’(t)+4y’(t)+2y(t)=e-3t+e-5t sin(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-5t sin(4t+π/3)
根据初值条件y(0)=1，y’(0)= y’’(0)=1/2，y(3)(t)=0.2即可得到如下的simulink仿真模型
(1)从simulink仿真模型的示波器Scope上观察到的输出波形y(t)
(2)通过MATLAB命令得到的输出波形y(t)
>> [t,x,y]=sim('simulink_Chapter5_2',250);
plot(t,y); grid
➢方法二：运用微分方程数值解进行分析
第一步：首先编写一个MATLAB函数来描述该微分方程，将原微分方程化成“方法一”中的四变量形式后即可编写出如下代码，将其在工作目录下保存为vdp_eq5_2.m文件%——————————————————————————————————————function y=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的工作区Command Window编写如下代码来求出微分方程的数值解，并将数值解[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 %——————————————————————————————————————➢分析结果
从以上仿真来看，无论是采用“方法一”还是“方法二”，均能够得到非常理想的仿真结果曲线。

因此对于该微分方程而言，这两种方法对于数值解的求解均是不错的选择
【3】
➢方法一：运用Simulink搭建系统的仿真模型进行分析
原微分方程y(4)(t)+5t y(3)(t)+6t2y’’(t)+4y’(t)+2e-2t y(t)=e-3t+e-5t sin(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-2t x1(t)-4x2(t)- 6t2x3(t)-5t x4(t)+e-3t+e-5t sin(4t+π/3)
根据初值条件y(0)=1，y’(0)= y’’(0)=1/2，y(3)(t)=0.2即可得到如下的simulink仿真模型
(1)从simulink仿真模型的示波器Scope上观察到的输出波形y(t)
(2)通过MATLAB命令得到的输出波形y(t)
>> [t,x,y]=sim('simulink_Chapter5_3',10);
plot(t,y); grid
➢方法二：运用微分方程数值解进行分析
第一步：首先编写一个MATLAB函数来描述该微分方程，将原微分方程化成“方法一”中的四变量形式后即可编写出如下代码，将其在工作目录下保存为vdp_eq5_3.m文件%——————————————————————————————————————function y=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的工作区Command Window编写如下代码来求出微分方程的数值解，并将数值解[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 %——————————————————————————————————————➢分析结果
从以上仿真来看，无论是采用“方法一”还是“方法二”，均能够得到非常理想的仿真结果曲线。

因此对于该微分方程而言，这两种方法对于数值解的求解均是不错的选择
【6】
(1)从示波器Scope观察所得波形
05101520253035404550-8-6-4-2
2468
2
From: In(1)1)From: In(2)From: In(3)From: In(4)
Step Response (2)通过在MATLAB 工作目录下输入命令得到波形
>> [t,x,y]=sim('simulink_Chapter5_6',50);
plot(t,y(:,1),'-',t,y(:,2),'--*');grid %“实线”的为输出1，“虚线且转折带星号”的为输出2
【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));grid
subplot(222),plot(tt,yy(:,2));grid subplot(223),plot(tt,yy(:,3));grid
subplot(224),plot(tt,yy(:,4));grid ➢ 通过step()函数对上述模型进行仿真
(1)方案一，8路输出，MATLAB 命令如下：
>>g11=tf(1,[4,1]);g12=tf(0.7,[5,1]);g13=tf(0.3,[5 1]);g14=tf(0.2,[5 1]);
g21=tf(0.6,[5,1]);g22=tf(1,[4,1]);g23=tf(0.4,[5 1]);g24=tf(0.35,[5 1]);
g31=tf(0.35,[5,1]);g32=tf(0.4,[5,1]);g33=tf(1,[4,1]);g34=tf(0.6,[5 1]);
g41=tf(0.2,[5,1]);g42=tf(0.3,[5,1]);g43=tf(0.7,[5 1]);g44=tf(1,[4,1]); G1=[g11 g12 g13 g14;g21 g22 g23 g24;g31 g32 g33 g34; g41 g42 g43 g44];
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)); 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) (3)最终的4路输出波形，MATLAB 命令如下
150
200250300350>> [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 仿真的示波器上只能看到两个波形的原因，这是由于两对相同
形状的波形叠加在一起所造成的
最终结论：两种情况下的仿真结果差异不大，在实际工程建模上都是适用的
【11】
x 1’=t sin (x 2e -2.3x 4)
x 2’=x 1
x 3’=sin (x 2e -2.3x 4)
x 4’=x 3
【13】
(1)将直流电机拖动系统的模型在Simulink 平台上搭建出来，提取系统的总模型
(2)绘制系统的阶跃响应、频域响应曲线
➢ 阶跃响应
>>[t,x,y]=sim('simulink_Chapter5_13',0.8);
plot(t,y);grid
➢ 频域响应
>> s=tf('s');
G=(1 【1】
-150
-100
-50
050
M a g n i t u d e (d B )10-1100101102103
-270
-180
-90
0P h a s e (d e g )Bode Diagram
Gm = 13.8 dB (at 7.95 rad/s) , P m = 60.1 deg (at 3.92 rad/s)
Frequency (rad/s) 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);hold on;margin(G*Gc) figure(f4);step(G);hold on;step(G*Gc);
调整前后比较：从Gc(s)的表达式不难发现，这是一个超前校正器。

通过比较校正前后系统
的Bode 图，可以发现校正后的系统的相位裕度反而减小了，此时系统稳定性明显减小了，但此时时间响应速度却反而增大了 原系统的
Bode 图 校正后系统的Bode 图
00.51 1.52 2.5
3 3.5
4 4.5502
46
8
10
12
Step Response
Time (seconds)A m p l i t u d e
分析：通过观察阶跃响应曲线，可以发现校正后的系统较原系统幅值明显增大，瞬态响应速度加快，但是缺点是其超调量变得过大，调整时间增加
建议：设计控制器应当充分考虑系统的稳定性，然后再考虑准确性和快速性
【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;
for gam=20:10:90
Gc=leadlagc(G1,wc,gam,1000,3);
figure(1);margin(Gc*G1);hold on; figure(2),step(feedback(Gc*G1,1),70);hold on;
end
校正后系统
的阶跃响应 原系统的
阶跃响应
0.20.40.6
0.811.21.4Step Response A m p l i t u d e 010203040506070
00.2
0.4
0.6
0.8
11.21.4
1.6
1.8
Step Response
Time (seconds)
A m p l i t u d e
分析：通过上图比较可以发现，随着相位裕度值的增加，闭环系统的超调量是始终在减小的，此时我们需要考虑的一个关键因素就是调整时间了，兼顾这两个因素，只有当这两者达到了最优，最终得到的系统才是理想的。

可以发现，当相位裕度值在50~60这个范围的时候其构成的闭环系统的动态特性达到了最好的效果，和原系统的动态特性相比，这个范围的优势在于其调整时间的增加并不明显，且超调量较原系统已经明显减小。

而再往下，调整时间又开始变长了，因此我们将在这个范围内搜素合适的相位裕度值
%第二阶段筛选合适的期望相位裕度值
>> wc=1;
for gam=50:5:60
Gc=leadlagc(G1,wc,gam,1000,3);step(feedback(Gc*G1,1),40);hold on;
end 60
% 分析：-
1.2
1.4Step Response
➢ 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;
for gam=20:10:50
Gc=leadlagc(G2,wc,gam,1000,3);
figure(1);margin(Gc*G2);hold on; figure(2),step(feedback(Gc*G2,1));hold on;
end
0123
45600.20.4
0.6
0.81
1.21.4
Step Response
Time (seconds)A m p l i t u d e
分析：通过上图比较可以发现，系统此时的振荡次数明显减小，调整时间已经从原来的将近200s 缩短到了现在的6s ，因此这种设计非常有效地减小了调整时间；此时我们只需找到在这个条件下合适的超调量足够小的响应曲线，即能确定出符合要求的期望相位裕度值；可以发现一点，如果我们将超调量那部分减小为负值时，则该闭环系统的振荡就变成仅仅一次了，且调整时间未发生明显变化，因此我们可以选择使调整量变为负值的那个响应曲线所对应的相位裕度值，通过比较，我们最终确定gam=45为期望的相位裕度值
%设计：超前-滞后校正器的最终指标定为剪切频率wc=25，相位裕度值gam=45
>> wc=25;gam=45; Gc=leadlagc(G2,wc,gam,1000,3)
step(feedback(G2*Gc,1));
r 增加方向
Transfer function:
1.289e004 s^2 + 8.328e004 s + 1.276e005
---------------------------------------
s^2 + 168.1 s + 1612
分析：此时闭环系统的超调量已经
消除，振荡已经减少到了一次，且
调整时间也由原来的200s缩短到了
6s，响应速度变快，且动态特性也
达到了最优，因此这个设计是合理
的
【8】
(1)Ga(s)=1/(s+1)^3
○1利用“整定公式”设计控制器
STEP 1:求出受控对象的一阶近似模型
分析：从得出的拟合结果可以看出，对于这个式子，采用次最优降阶方法得出的拟合结果最接近，因此就选择该方法所得的表达式
STEP 2:根据整定公式设计PID控制器
根据由该受控对象得出的较好的FOLPD，可知K=1；T=1.997；L=1.13
结论：所以由此设计的PID控制器模型为Gc(s)=2.1207(1+1/(2.2600s)+0.5650s/(0.0565s+1)) STEP 3:分析给出的受控对象模型在控制器下的阶跃响应曲线
分析：可见最终响应的超调量小，振荡少，调整时间快，达到了最终的设计目的
○2利用PID控制器设计控制器
根据由该受控对象得出的较好的FOLPD，可知K=1；T=1.997；L=1.13，因此可以设计如下PID 控制器
分析：在该PID控制器作用下，可以发现最终响应的超调量小，振荡少，调整时间快，达到了最终的设计目的
(2)Gb(s)=1/(s+1)^5
○1利用“整定公式”设计控制器
STEP 1:求出受控对象的一阶近似模型
分析：从得出的拟合结果可以看出，对于这个式子，采用次最优降阶方法得出的拟合结果最接近，因此就选择该方法所得的表达式
STEP 2:根据整定公式设计PID控制器
根据由该受控对象得出的较好的FOLPD，可知K=1；T=2.624；L=2.59
结论：所以由此设计的PID控制器模型为Gc(s)=1.2158(1+1/(5.1800s)+1.2950s/(0.12950s+1)) STEP 3:分析给出的受控对象模型在控制器下的阶跃响应曲线
分析：可见最终响应的超调量小，振荡较少，调整时间快，达到了最终的设计目的
○2利用PID控制器设计控制器
根据由该受控对象得出的较好的FOLPD，可知K=1；T=2.624；L=2.59，因此可以设计如下PID 控制器
分析：在该PID控制器作用下，可以发现最终响应的超调量小，振荡较少，调整时间快，达
到了最终的设计目的
(3)Gb(s)=(-1.5s+1)/(s+1)^3
○1利用“整定公式”设计控制器
STEP 1:求出受控对象的一阶近似模型
分析：从得出的拟合结果可以看出，对于这个式子，采用次最优降阶方法得出的拟合结果最接近，因此就选择该方法所得的表达式
STEP 2:根据整定公式设计PID控制器
根据由该受控对象得出的较好的FOLPD，可知K=1；T=1.769；L=2.52
结论：所以由此设计的PID控制器模型为Gc(s)=0.8424(1+1/(5.0400s)+1.2600s/(0.12600s+1)) STEP 3:分析给出的受控对象模型在控制器下的阶跃响应曲线
分析：可见最终效果并不是特别明显，阶跃响应的超调量显然已经减小不少，但是其振荡还是过多，调整时间还算比较快的，达到了最终的设计目的，总的来说还是不错的
○2利用PID控制器设计控制器
根据由该受控对象得出的较好的FOLPD，可知K=1；T=1.769；L=2.52，因此可以设计如下PID 控制器
分析：在该PID控制器作用下，可以发现最终效果并不是特别明显，阶跃响应的超调量显然已经减小不少，但是其振荡还是过多，调整时间还算比较快的，达到了最终的设计目的，总的来说还是不错的。