Matlab控制系统计算机辅助设计

实 验 目 录

实验一：Matlab环境熟悉与基本运算（设计型） 实验二：Matlab语言程序设计(设计型） 实验三：控制系统模型的建立（设计型） 实验四：Simulink仿真入门（验证型） 实验五：控制系统时域仿真分析（设计型） 实验六：Simulink环境下时域仿真 实验七：控制系统根轨迹仿真分析

实验八：控制系统频域仿真分析（设计型）

实验名称：Matlab环境熟悉与基本运算（设计型）

1、矩阵运算 （1）矩阵的乘法 A=[1 2;3 4]; B=[5 5;7 8]; y=A^2*B

（2）矩阵除法

A=[1 2 3;4 5 6;7 8 9]; B=[1 0 0;0 2 0;0 0 3]; y1=A\\B

警告: 矩阵接近奇异值，或者缩放错误。结果可能不准确。RCOND = 1.541976e-18。

y2=A/B

（3）矩阵的转置及共轭转置

A=[5+i,2-i,1;6*i,4,9-i]; y1=A.'

y2=A'

y =

105 115 229 251

y1 =

1.0e+16 *

-0.4504 1.8014 -1.3511 0.9007 -3.6029 2.7022 -0.4504 1.8014 -1.3511 y2 =

1.0000 1.0000 1.0000 4.0000 2.5000 2.0000 7.0000 4.0000 3.0000 y1 =

5.0000 + 1.0000i 0.0000 + 6.0000i 2.0000 - 1.0000i 4.0000 + 0.0000i 1.0000 + 0.0000i 9.0000 - 1.0000i

y2 =

5.0000 - 1.0000i 0.0000 - 6.0000i 2.0000 + 1.0000i 4.0000 + 0.0000i 1.0000 + 0.0000i 9.0000 + 1.0000i

（4）使用冒号选出指定元素

A=[1 2 3;4 5 6;7 8 9]; y1=A(1:2,3)

y2=A(2:3,:)

（5）复数矩阵的生成： syms a1 a2 a3 a4 b1 b2 b3 b4 A=[a1+b1*j a2+b2*j;a3*exp(b3*j) a4*exp(b4*j)]

a1=3;a2=-2;a3=9;a4=23;b1=5;b2=3;b3=6;b4=33;

A=[a1+b1*j a2+b2*j;a3*exp(b3*j) a4*exp(b4*j)]

2、多项式 p =

1 0 2 4

roots(p)%用roots函数求多项式的根

①A=[1.2 3 5 0.9;5 1.7 5 6;3 9 0 1;1 2 3 4]; P=poly(A); poly2sym(P)

②polyval(P,20)

③polyvalm(P,A)

y1 =

3 6

y2 =

4 5 6 7 8 9

A =

[ a1 + b1*1i, a2 + b2*1i] [ a3*exp(b3*1i), a4*exp(b4*1i)] A =

3.0000 + 5.0000i -2.0000 + 3.0000i 8.15 - 2.5147i -0.3054 +22.9980i

ans =

0.58 + 1.7445i 0.58 - 1.7445i -1.1795 + 0.0000i ans =

x^4 - (69*x^3)/10 - (3863*x^2)/50 - (8613*x)/100 + 12091/20

ans =

7.2778e+04 ans =

1.0e-11 *

-0.4093 -0.4849 -0.3876 -0.4519 -0.5002 -0.8072 -0.6004 -0.5684 -0.4704 -0.6196 -0.80 -0.5230 -0.3297 -0.4455 -0.3595 -0.3638

3、基本绘图命令 （1）t=0:0.01:2*pi; y=cos(t); plot(t,y)

运行结果如右图所示：

（2）

t=0:0.01:2*pi; y1=cos(t-0.25); y2=sin(t-0.5); plot(t,y1,t,y2)

运行结果如右图所示：

4、基本绘图控制 t=0:0.01:4*pi; x1=10*sin(t); plot(t,x1,'-.b*')

axis([-20 20 -20 20]) grid on

title('正弦曲线x1=10*sin(t)图') xlabel('t=0:0.01:4*pi') ylabel('x1')

text(pi/2,10,'最大值') 运行结果如右图所示：

实验名称：Matlab语言程序设计（设计型）

1、编写命令文件：计算 1+2+…+n<5000 时的最大 n 值； 解：

建立M文件，文件名称记为two1，编

辑程序内容如下所示： 保存，并在命令窗口中输入two1，运行s=0; 结果如下所示： n=1; two1 while(s<5000), s=s+n; n = n=n+1; end, 101

n

2、编写函数文件： 解：

①用for循环结构编写程序：

建立函数文件，文件名称记为findsum1，编辑程序内容如下所示：

function [ s1 ] =findsum1( ) %findsum1 用for循环编写的程序 % 求2的0到30次幂的和,s表示和，i表示执行次数 s1=0; for i=0:30, s1=s1+2^i; End end

保存，并在命令窗口中输入

[s1]=findsum1( )，运行结果如下所示： [s1]=findsum1( ) s1 =

2.1475e+09

②用while循环结构编写程序：

建立函数文件，文件名记为findsum2，编辑程序内容如下所示：

function [ s2 ] =findsum2( ) %findsum2 用while循环编写的程序 % 求2的0到30次幂的和,s表示和，i表示执行次数 s2=0; i=0;

while(i<=30),

s2=s2+2^i; i=i+1; end end

保存，并在命令窗口中输入[ s2 ] =findsum2( )，运行结果如下所示： [ s2 ] =findsum2( ) s2 =

3、解：

建立M文件，文件名记为two3，编辑程序内容如下：

x=input('please input a character:\\n','s'); if x=='y'||x=='Y', disp('x=1'); else if x=='n'||x=='N',

2.1475e+09

disp('x=0'); else

disp('its erro'); end end

保存，并在命令窗口中输入two2，运行结果如下所示： its erro two2

please input a character: Y x=1

two2

please input a character: n x=0 two2

please input a character: S

实验名称：控制系统模型的建立和分析（设计型）

1、在MATLAB环境中输入下面的系统模型 解：

建立M文件，文件名称记为thr1，编辑程序内容如下所示：

num=[1,4,3,2]; %系

s34s23s2G(s)2s(s1)(s4)24统模型的分子系数

den=conv([1,0,0],conv([1,1],[1,8,20]));%分母系数 G1=tf(num,den);%得到开环系统模型 G2=feedback(G1,1);%得到闭环系统模型

[p1,z1]=pzmap(G1)%不加分号，在命令窗口可显示出开环系统的零极点 [p2,z2]=pzmap(G2)% 不加分号，在命令窗口可显示出闭环系统的零极点

在命令窗口输入thr1，则得出以下结果： thr1

p1 =

0 0 -4.0000 + 2.0000i -4.0000 - 2.0000i -1.0000

z1 =

-3.2695 -0.3652 + 0.6916i -0.3652 - 0.6916i

2、假设

s23s5;串联连G1(s)4s4s33s22s1

且

p2 =

-3.9174 + 2.1007i -3.9174 - 2.1007i -1.1449 -0.0102 + 0.2972i -0.0102 - 0.2972i

z2 =

-3.2695 -0.3652 + 0.6916i -0.3652 - 0.6916i

系接

s1G1(s)2;而成，已知

s3s4统由两个模块和

解：建立M文件，文件名称为thr2，编辑程序内容如下所示：

G3=tf([1,1],[1,3,4]);%实验三的第二题，题中的G1（s）模型 G4=tf([1,3,5],[1,4,3,2,1]);% 题中的G2（s）模型

G5=series(G3,G4)%用串联函数，将G3和G4模型串联，得到G5模型，在命令窗口显示出模型 Gss=ss(G5)%用ss得到G5模型的状态方程模型

在命令窗口中输入thr2，得到以下运行结果： Transfer function:

s^3 + 4 s^2 + 8 s + 5

------------------------------------------------- s^6 + 7 s^5 + 19 s^4 + 27 s^3 + 19 s^2 + 11 s + 4 a =

x1 x2 x3 x4 x5 x6

x1 -7 -2.375 -1.688 -0.5938 -0.3438 -0.25 x2 8 0 0 0 0 0 x3 0 2 0 0 0 0 x4 0 0 2 0 0 0 x5 0 0 0 1 0 0 x6 0 0 0 0 0.5 0 b =

u1 x1 1 x2 0 x3 0 x4 0 x5 0 x6 0 c =

x1 x2 x3 x4 x5 x6 y1 0 0 0.0625 0.125 0.25 0.3125 d =

u1 y1 0

Continuous-time model. 3、假设系统的对象模型为 G(s)10

(s1)3,并定义一个PID控制器10.4353sGPID(s)0.4811.814s0.04353

解：

建立M文件，文件名称为thr3，编辑程序内容如下所示：

s=tf('s');%拉氏算子

G6=10/((s+1)^3);%实验三第3题的对象模型

Gpid=0.48*(1+1/(1.814*s)+0.4353*s/0.04353);%PID控制器 G7=series(G6,Gpid);%控制器与对象模型进行串联得到的开环系统模型 G8=feedback(G7,1)%闭环系统的传递函数模型

[z3,p3,k3]=zpkdata(G7,'v')%得到开环系统的零极点及增益 [z4,p4,k4]=zpkdata(G8,'v') %得到闭环系统的零极点及增益

在命令窗口中输入thr3，得到以下运行结果 Transfer function:

3.79 s^2 + 0.379 s + 0.20

------------------------------------------------------- 0.076 s^4 + 0.2369 s^3 + 4.027 s^2 + 0.458 s + 0.20 z3 = z4 = -0.0500 + 0.2294i -0.0500 + 0.2294i -0.0500 - 0.2294i -0.0500 - 0.2294i p3 = p4 = 0 -1.4442 + 6.9670i -1.0000 -1.4442 - 6.9670i -1.0000 + 0.0000i -0.0558 + 0.2217i -1.0000 - 0.0000i -0.0558 - 0.2217i k3 = k4 = 48 48

4、解：三种典型模型建模函数：

①传递函数模型：G=tf(num,den)； ②零极点增益模型：G=zpk(z,p,k)；

各种模型之间的转换语句： [A,B,C,D]=tf2ss(num,den); [num,den]=ss2tf(A,B,C,D,iu); [z,p,k]=tf2zp(num,den);

③状态空间模型：G=ss(A,B,C,D)；

[num,den]=zp2tf(z,p,k); [A,B,C,D]=zp2ss(z,p,k); [z,p,k]=ss2zp(A,B,C,D,i)

实验名称：Simulink仿真环境入门（验证型）

1、搭建下图所示控制系统模型，自己设计输入正弦曲线参数，观察示波器及输出数据。 解：画出的仿真图如下，Sine Wave:O(t)=3*sin(2*t)，增益为5。

出现以下波形图：左图为使用plot函数画出的figure图，右图为示波器显示的图。

实验名称：控制系统时域仿真（设计型）

1、时域分析 （1）解：

建立M文件，写入以下程序： G1=tf([5,25,30],[1,6,10,8]); t=[0:0.0005:25]; [y1,t]=step(G1,t); [ymax,tp]=max(y1); max_overshoot=ymax-1 [y2,t]=impulse(G1,t); plot(t,y1,t,y2)

figure图如右图所示：

（2）解：

建立M文件，写入以下程序：

wn=[2:2:12]; z=0.7;

t=0:0.0005:25; hold on

for i=1:length(wn)

G2=tf(wn(i)^2,[1,2*z*wn(i),wn(i)^2]); step(G2,t) end

grid on,hold off

legend('wn=2','wn=4','wn=6','wn=8','wn=10','wn=12') figure图如右图所示： （3）解：建立M文件，写入以下程序： wn=6; zetas=[0.2 0.4 0.6 0.8 1.0 1.5 2.0]; t=0:0.0005:25; hold on for i=1:length(zetas)

G3=tf(wn^2,[1,2*zetas(i)*wn,wn^2]); step(G3,t) end

grid on,hold off

legend('zetas=0.2','zetas=0.4','zetas=0.6','zetas=0.8','zetas=1.0','zetas=1.5','zetas=2.0')

figure图如右上图所示：

2、系统稳定性分析 （1）解：

①零极点分布图法程序：

G4=zpk(-2,[0;-1;-20],100) G5=feedback(G4,1) pzmap(G5)

figure图：

由图可知该系统的极点均位于左半平面，所以该系统稳定。

②程序判定法程序：

G4=zpk(-2,[0;-1;-20],100) G5=feedback(G4,1) p=pole(G5); for i=1:length(p) if real(p(i))>0 disp('系统不稳定') break end end

disp('系统稳定') pzmap(G5)

figure图：

实验名称：Simulink环境下时域仿真

1、解：仿真图：

阶跃响应： 方波响应： 2、解：仿真图：

阶跃响应： 正弦信号输入响应：

斜坡输入响应：

3、 解：仿真图：

建立M文件，编写以下程序：

G1=tf([1 1],[2 1]); G2=tf([5],[2 3 1]); H1=tf([1],[2 1]); G3=series(G1,G2); G4=feedback(G3,H1) isstable(G4)

得到以下实验结果： g3

Transfer function:

8 s^4 + 20 s^3 + 18 s^2 + 12 s + 6

ans =

10 s^2 + 15 s + 5 ----------------------------------

1

实验名称：控制系统根轨迹分析 结论：编写程序判断该系统的稳定性结果是该系统是稳定的。

在Simulink环境下验证编程结果： 给出阶跃信号，画出以下仿真图进行仿真：

示波器显示结果：

从上图可以看出该系统是是稳定的，所以编程结果是正确的。

1、

（1）解：编写程序如下所示：

G=zpk(-0.5,[0;-1;-2;-5],1)

（2）解：编写程序如下所示：

G=zpk([],[0;-2],1) rlocus(G) [k,p]=rlocfind(G)

Figure图：系统的根轨迹图：

rlocus(G)

[k,p]=rlocfind(G)

Figure图：系统的根轨迹图：

2、根轨迹分析

（1）解：编写程序如下所示：

G=zpk(-3,[0;-5;-6;-1+j;-1-j],1) rlocus(G)

[k,p]=rlocfind(G)

（2）解：编写以下程序：

G=zpk([],[0;-1;-2],1) rlocus(G) k=0:0.01:98; [r,k]=rlocus(G,k); [m,n]=size(r); for i=1:n

if real(r(2,i))>0 break end end kg=k(i-1)

Figure图：系统的根轨迹图：

由根轨迹图可以看出点击的这一点的增益k=35.2；闭环极点的位置

p=-0.00587+1.35i，因为p<0，所以该点系统闭环是稳定的。

程序的执行结果：

Zero/pole/gain: 1

------------- s (s+1) (s+2) kg =

6

由上述的执行结果可以看出该系统稳定的k值范围是当k<6时，系统是稳定的。 Figure图如右图所示：

实验名称：控制系统频域仿真分析（设计型）

1、（1）解：编写以下程序：

zeta=0.7;

wn=[2,4,6,8,10,12]; for i=1:length(wn) num=[wn(i)^2]; den=[1,2*zeta*wn(i),wn(i)^2]; G1=tf(num,den) bode(G1) end legend('wn=2','wn=4','wn=6','wn=8','wn=10','wn=12') 得到右图所示运行结果： （2）解：编写以下程序： wn=6; zeta=[0.2,0.4,0.6,0.8,1.0,1.5,2.0]; hold on for i=1:length(zeta) num=[wn^2]; den=[1,2*zeta(i)*wn,wn^2]; G2=tf(num,den) bode(G2) end hold off

legend('zeta=0.2','zeta=0.4','zeta=0.6','zeta=0.8','zeta=1.0','zeta=1.5','zeta=2.0')

得到右图所示运行结果： （3）

（a）编写以下程序：

num=[10];

den=conv([1,0,0],conv([5,1],[1,5])); G1=tf(num,den) nyquist(G1) bode(G1)

得到右图所示运行结果： （b）编写以下程序：

num=conv(8,[1,1]);

den=conv([1,0,0],conv([1,15],[1,6,10]));

G2=tf(num,den) nyquist(G2) bode(G2)

得到右图所示运行结果： （c）编写以下程序：

num=conv(4,[1/3,1]);

den=conv([1,0],conv([0.02,1],conv([0.05,1],[0.1,1]))); G3=tf(num,den) nyquist(G3) bode(G3)

得到右图所示运行结果：

2、（1）解：编写以下程序：

G1=tf(2.7,[1,5,4,0]) G2=tf(2.7,[1,5,-4,0]) figure(1) margin(G1) figure(2) margin(G2)

得到以下运行结果：

由图可以分析出G1的系统是稳定的，G2系统不稳定。

Bode DiagramGm = 17.4 dB (at 2 rad/sec) , Pm = 51.7 deg (at 0.578 rad/sec)50Magnitude (dB)Phase (deg)0-50-100-150-90-135-180-225-27010-210-1100101102Frequency (rad/sec)Bode DiagramGm = Inf , Pm = -58.1 deg (at 0.535 rad/sec)50Magnitude (dB)Phase (deg)0-50-100-150-210-240-27010-210-1100101102Frequency (rad/sec)