基于MATLAB的系统仿真实验指导书

基于MATLAB的系统仿真实验

实验指导书

说明

MATLAB是MATrix LABoratory的缩写，是一种基于矩阵的数学与工程计算系统，可以用作动态系统的建模与仿真。研究系统的结构和参数的变化对系统性能的影响可利用MATLAB强大的计算和作图功能，因此本实验采用MATLAB仿真研究连续控制系统和离散控制系统的性能分析过程。通过该实验提高学生对控制系统的分析与设计能力，加深对《自动控制原理》课程内容的理解。

该实验教学大纲中只安排了两个学时，本实验的内容比较多，要求同学们参考有关MATLAB应用的书籍利用课后时间完成。同时可以将MATLAB用于习题的验证和课程的学习中。

仿真实验

基于MATLAB 连续系统的仿真

一、 实验目的

在研究系统的结构和参数的变化对系统性能的影响时，采用解析和作图的方法比较麻烦，而且误差也大，用MATLAB 仿真实现则简单方便，精度高。本实验采用MATLAB 实现控制系统的数学描述、控制系统的时域分析及根轨迹和频率特性分析。通过该实验，加深学生对系统阶次，型号，参数与系统性能的关系的理解。 二、实验环境

在计算机Windows 环境下安装好MATLAB6.3以上版本后，双击MATLAB 图标或成“开始”菜单打开MATLAB ，即可进入MATLAB 集成环境。 三、实验原理

三、MATLAB 应用实例

1．拉氏变换和反变换

例 求22)(2

++=t t t f 的拉氏变换 解

键入 syms s t;

ft=t^2+2*t+2; st=laplace(ft,t,s)

运行结果为 st=

2/s^3+2/s^2+2/s

例 求)

2)(34(6

)(2

++++=

s s s s s F 的拉氏反变换 解

键入 syms s t;

Fs=(s+6)/(s^2+4*s+3)/(s+2); ft=ilaplace(Fs,s,t)

运行结果为 ft=

3/2*exp(-3*t)+5/2*exp(-t)-4*exp(-2*t)

2。求根运算

例 求多项式 43)(23++=s s s p 的根，再由根建多项式。 解

键入 p=[1 3 0 4]; r=root(p) 运行结果为 r=

-3.3553

0.1777+1.0773i 1.7777-1.0773i

键入 p=poly(r)

运行结果为 p=

1.0000 3.0000 0.0000 40000

例 实现多项式相乘：)4)(123(2+++s s s ，并求5-=s 时的值。 解

键入 p=[3 2 1]; q=[1 4]; n=conv(p,q) 运行结果为 n=

3 1

4 9 4 键入 vlaue=polyval(n,-5) 运行结果为 value=

-66

3．微分方程求解

例 解下列微分方程：

1)(2)(3)(32

2=++t y dt t dy dt

t y d 初始条件 0)0()0('==y y 解

键入 y=dsolve('3*D2y+3*Dy+2*y=1','y(0)=0,Dy(0)=0') 运行结果为 y=

1/2-1/2*exp(-1/2*t)*cos(1/6*15^(1/2)*t)-1/10*15^(1/2)*exp(-1/2*t)*sin(1/6*15^(1/2)*t)

4．传递函数

例 求如图所示系统的传递函数。

解

键入numg=[1];deng=[500 0 0]; numh=[1 1];denh=[1 2];

[num,den]=feedback(numg,deng,numh,denh,-1); printsys(num,den) 运行结果为 num/den=

s +2

------------------------------

500s^3 + 1000 s^2 + s+1

其中， 由函数 printsys(num,den) 打印出传递函数。

5。系统输出响应及性能分析

例 二阶系统闭环传递函数的标准形式为 2

2

22)()

()(n

n n s s s R s C s ωξωω++==Φ 若n ω确定，系统的瞬态响应和ξ的取值有关。下面用MATLAB 分析在不同的ξ值时，系统的单位阶跃响应。所用的MA TLAB 程序如下：

%n ω=1

t=0:0.1:12;num=[1];

zeta1=0;den1=[1 2*zeta1 1]; zeta3=0.3;den3=[1 2*zeta3 1]; zeta5=0.5;den5=[1 2*zeta5 1]; zeta7=0.7;den7=[1 2*zeta7 1]; zeta9=1.0;den9=[1 2*zeta9 1]; [y1,x,t]=step(num,den1,t); [y3,x,t]=step(num,den3,t); [y5,x,t]=step(num,den5,t); [y7,x,t]=step(num,den7,t); [y9,x,t]=step(num,den9,t); plot(t,y1,t,y3,t,y5,t,y7,t,y9,) grid on; 运行结果见图。

例 已知A

A

K s s K s 55.345)(2++=

Φ，分别计算K A ＝200时，系统的性能指标t p ，t s ，%σ。

解 MA TLAB 程序如下：

t=0:0.01:2; num=[1000];

den=[1 34.5 1000]; [y,x,t]=step(num,den,t); plot(t,y); ％求超调量 maxy=max(y); yss=y(length(t));

pos=100*(maxy-yss)/yss; ％求峰值时间 for i=1:1:201

if y(i)==maxy,n=i;end end

tp=(n-1)*0.01; ％求调节时间 for i=1:1:201

if (y(i)<1.05 & y(i)>0.95),m=i;end break; end

ts=(m-1)*0.01;

计算结果为

t p =0.12s ，t s =0.17s ，%σ＝12.93%

例 系统闭环特征方程为 0242)(2

3

=+++=s s s s q 用MA TLAB 判系统稳定性。 解

键入 den=[1 1 2 24]; roots(den) 运行结果为 ans=

-3.0000

1. 000+

2.6458i 1.0000-2.6458i

可见，系统有两个右根，故系统不稳定。 例 求出下列系统11)(1+=

s s G ，)1(1)(2+=s s s G ，)

1(14)(23++=s s s s G 单位阶跃 响应及稳态误差。 解

MATLAB 程序如下：

t=0:0.1:20;

[num1,den1]=cloop([1],[1 1]); [num2,den2]=cloop([1],[1 1 0]); [num3,den3]=cloop([4 1],[1 1 0 0]); y1=step(num1,den1,t); y2=step(num2,den2,t); y3=step(num3,den3,t); subplot(311); plot(t,y1);

subplot(312);plot(t,y2); subplot(313);plot(t,y3); er1=y1(length(t))-1; er2=y2(length(t))-1; er3=y3(length(t))-1;

运行结果见图。

在命令窗口可得

er1 = -0.5000 % 0型系统的稳态误差 er2 = 2.4294e-005 % Ⅰ型系统的稳态误差 er3 = -3.3525e-004 % Ⅱ型系统的稳态误差

例 用SIMULINK 对第七节中单闭环无静差调速系统进行仿真。 解 在SIMULINK 下建立系统仿真结构图如图3-46所示。

仿真的结果见图3-47。

6。绘制系统的根轨迹图

例 设一单位反馈系统如图4-31所示，试绘制该系统的根轨迹。

解

键入num=[1];

den=[1 3 2 0]; rlocus(num,den) 运行结果见图4-32。

例 已知系统前向通道传递函数为5)(+=

s K s G c ，)

8(1

)(0++=s s s s G ，反馈通道传递函数为2

1

)(+=

s s H ，试绘出系统闭环根轨迹图。

解

键入 Gc=tf(1,[1 5]);

Go=tf([1 1],[1 8 0]); H=tf(1,[1 2]); rlocus(Gc*Go*H);

v=[-10 10; –10 13];axis(v); grid on

运行结果见图4-33。

图4-33 例4-18系统根轨迹

7．绘制系统的频率特性图

例5 振荡环节传递函数的标准形式为

2

2

22)()

()(n

n n s s s R s C s ωζωω++==Φ 当n ω确定时，系统的频率特性与ζ的取值有关。用MATLAB 绘出1=n ω时，在不同ζ取值下，系统的Bode 图和Nyquist 图。 解 1）作Bode 图程序： num=[1];

zeta1=0.1;den1=[1 2*zeta1 1]; zeta3=0.3;den3=[1 2*zeta3 1];

zeta5=0.5;den5=[1 2*zeta5 1]; zeta7=0.7;den7=[1 2*zeta7 1]; zeta9=1.0;den9=[1 2*zeta9 1];

[mag1,phase1,w1]=bode(num,den1); [mag3,phase3,w3]=bode(num,den3); [mag5,phase5,w5]=bode(num,den5); [mag7,phase7,w7]=bode(num,den7); [mag9,phase9,w9]=bode(num,den9); subplot(211);

semilogx(w1,20*log10(mag1),w3,20*log10(mag3),…

w5,20*log10(mag5),w7,20*log10(mag7), w9,20*log10(mag9)); subplot(212);

semilogx(w1,phase1,w3,phase3,w5,phase5,w7,phase7,w9,phase9); 运行结果见图。

2）作Nyquist 图程序： num=[1];

zeta1=0.4;den1=[1 2*zeta1 1]; zeta6=0.6;den6=[1 2*zeta6 1]; zeta8=0.8;den8=[1 2*zeta8 1]; [re1,im1]=nyquist(num,den1); [re2,im2]=nyquist(num,den6);

[re3,im3]=nyquist(num,den8);

plot(re1,im1,re2,im2,re3,im3); 运行结果见图。

例 已知系统的开环传递函数)

12()

10()(++=s s s s G ，求开环对数频率特性曲线。

解 MATLAB 程序如下：

z=[-10];p=[0,-1/2];k=1/2; h=zpk(z,p,k); bode(h); grid

运行结果见图。

例 开环传递函数为)

100167.0)(10368.0)(1162.0(33

.2)()(+++=

s s s s H s G

作出开环伯德图，并求系统的稳定裕量。

解 命令窗口输入

h1=tf([2.33],[0.162 1]); h2=tf([1],[0.0368 1]); h3=tf([1],[0.00167 1]); h=h1*h2*h3;

[num,den]=tfdata(h);

[mag,phase,w]=bode(num,den); subplot(211);

semilogx(w,20*log10(mag));grid; subplot(212);

semilogx(w,phase);grid;

[gm,pm,wcg,wcp]=margin(mag,phase,w) 运行结果见图。

同时，在MATLAB 命令窗口中可以得到系统的稳定裕量：

gm = 53.9692 pm = 93.6601 wcg = 141.9068 wcp = 11.6305

若要生成带有裕量标记的伯德图，则可将程序中最后一条改为

margin(mag,phase,w); 此时，所得图形见图。

8．系统性能的校正

例 试用MA TLAB 对教材中例6-1中系统设计超前校正装置。 解 如例6-1中所述，校正前系统的开环传递函数为 )

15.0(20

)(+=

s s s G

编写MATLAB 程序如下： num=20;den=[0.5 1 0];

[gm,pm,wcg,wcp]=margin(num,den); %求校正前系统的相位裕量pm

dpm=50-pm+5; %根据性能指标要求确定m ? phi=dpm*pi/180;

a=(1+sin(phi))/(1-sin(phi)); %根据式（6-6）求α mm=–10*log10(a); %计算-10lg α幅值 [mu,pu,w]=bode(num,den);

mu_db=20*log10(mu); ％在未校正系统的幅频特性上 wc=spline(mu_db,w,mm); %找到幅值为mm 处的频率 T=1/(wc*sqrt(a)); ％根据式（6-11）求T p=a*T;

nk=[p,1];dk=[T,1]; %根据式（6-8）求校正装置 gc=tf(nk,dk);

printsys(nk,dk,'s');

可从命令窗口得到校正装置的传递函数为

0.2268 s + 1 ------------ 0.0563 s + 1

再输入下面的命令

h=tf(num,den);h1=tf(nk,dk); g=h*h1;

[gm1,pm1,wcg1,wcp1]=margin(g); 可以得到校正后系统的相位裕量

pm1 = 49.7706

满足设计要求。

例 试用MA TLAB 对教材中例6-2中的系统设计滞后校正装置。 解 如例6-2中所述，校正前系统的开环传递函数为

)

15.0)(1(5

)(++=

s s s s G

编写MATLAB 程序如下： num=5;

den=[0.5 1.5 1 0];

[gm,pm,wcg,wcp]=margin(num,den); %求校正前的相位裕量 dpm=-180+40+12; %求?

[mu,pu,w]=bode(num,den);

wc=spline(pu,w,dpm); %在校正前曲线上找到与?对应的频率wc mu_db=20*log10(mu); %将幅值转化成以dB 为单位 m_wc=spline(w,mu_db,wc); %求未校正系统在wc 处的幅值 beta=10^(-m_wc/20); %求β值

w2=0.2*wc; %求滞后校正网络的转折频率2ω T=1/(beta*w2); %求T

nk=[beta*T 1];dk=[T 1]; %求校正网络的分子和分母系数

h1=tf(nk,dk); printsys(nk,dk,'s');

可得校正装置的传递函数为

10.76 s + 1 ----------- 102.3 s + 1

再输入下面的命令 h=tf(num,den); g=h*h1;

[gm1,pm1,wcg1,wcp1]=margin(g); 可得校正后的相角裕量为

pm1 = 41.5141

四、实验内容

1. 设典型二阶系统的闭环传递函数为2

222)(n

n n

s s s ωξωω++=Φ， 改变系统的参数 n ωξ,，观察并记录输出响应曲线，并分析n ωξ,的变化对系统性能的影响和性能

指标的变化。

2. 设系统的闭环传递函数为42)()(2+++=

Φs s a s K s ，5

.023)

()(2

3++++=Φs s s a s K s 改变系统的参数a K ,，观察并记录输出响应曲线，并求系统的时域性能指标，分析参数的变化对系统性能的影响和性能指标的变化。

3. 绘制下列系统的根轨迹图，分析开环零极点对系统性能的影响。 （1）25.33)1()(2+++=

s s s K s G r (2) )

3)(1()(++=

s s s K s G r (3) )3)(2()

5()(+++=

s s s s K s G r (4) )

52)(5()2()(2

++++=s s s s s K s G r 4. 绘制下列系统的伯德图和乃奎斯特图，并求系统的频域指标相位裕量和幅值裕量， 根据乃奎斯特图分析系统的稳定性，根据仿真结果分析系统的结构和参数与系统性能 的关系。 (1) 1

2.001.01)(2

++=

s s s G (2) )17.0)(15.0(1

)(++=s s s s G

(3) )15.0)(1()12.0(10)(+++=

s s s s s G (4) )

12.001.0)(15.0(10

)(2

+++=s s s s s G 5. 设系统的开环传递函数)

18.0()(+=

s s K

s G ，采用超前校正，使校正后系统速度误

差系数100=v K ，相位裕量 45≥γ，用MATLAB 设计校正装置的参数。

第二部分 基于MATLAB 离散系统的仿真

一、 实验目的

本实验采用MATLAB 实现对离散控制系统的性能分析。通过该实验，加深学生对系统阶次，型号，参数与系统性能的关系的理解。 二、实验环境

在计算机Windows 环境下安装好MATLAB6.3以上版本后，双击MATLAB 图标或成“开始”菜单打开MATLAB ，即可进入MATLAB 集成环境。 三、MATLAB 应用实例 1．z 变换和z 反变换

例 求单位斜坡函数t t f =)(的z 变换 解 命令窗口中键入

syms t T; ztrans(t*T)

运行结果为

ans =

T*z/(z-1)^2

即有

2

)

1(-z Tz

求表达式F 的z 反变换可用 f=iztrans(F)

2. 连续环节的离散化 例：已知：)(,1,1

)(2z G s T s

s s G 求=+=

h=tf([1 ],[1 1 0]); hd=c2d(h,1) 运行结果为： Transfer function:

0.3679 z + 0.2642 ---------------------- z^2 - 1.368 z + 0.3679

Sampling time: 1 3.脉冲传递函数

例：368.0)(,1)(21

-=-=z z

z G z z z G ，试求两个环节串联的开环脉冲传递函数

用MA TLAB 编写程序如下： h1=tf([1 0],[1 -1],[]);

h2=tf([1 0],[1 -0.368],[]); h3=series(h1,h2) 运算结果： Transfer function: z^2 --------------------- z^2 - 1.368 z + 0.368

例：368.0)(,1)(21

-=-=z z

z G z z z G ，试求两个环节并联的开环脉冲传递函数。

用MA TLAB 编写程序如下： h1=tf([1 0],[1 -1],[]); h2=tf([1 0],[1 -0.368],[]); h3=parallel(h1,h2) 运行结果如下： Transfer function: 2 z^2 - 1.368 z --------------------- z^2 - 1.368 z + 0.368 例 设)(,1z G s T 求=。

键入MA TLAB 命令

sys=tf([1],[1 1 0]); c2d(sys,1) 运行结果为

Transfer function: 0.3679 z + 0.2642 --------------------------

z^2 - 1.368 z + 0.3679

4.系统的输出响应

例 如图所示的采样控制系统，设T=1s ，单位阶跃输入，求系统的动态响应。

求得系统的闭环脉冲传递函数为

368

.0736.0632.0)()(2

+-=z z z

z R z C 令输入形式为 )51,1(ones u =

MATLAB 程序如下：

num=[0.632,0];den=[1,-0.736,0.368]; u=ones(1,51);k=0:50; y=filter(num,den,u); plot(k,y),grid;

xlabel('k');ylabel('y(k)');

运行结果见图：

例 系统结构如图所示，求系统的单位阶跃响应。

令输入为)51,1(ones u =

MATLAB 程序如下：

g=tf([1],[1 1 0]); %对象传递函数

d=c2d(g,1); %用零阶保持器离散化 cd=d/(1+d); %求闭环传递函数

cd1=minreal(cd); %传递函数约去公因子 [num,den]=tfdata(cd1,'v'); %求得分子、分母系数 u=ones(1,51);k=0:50;

y=filter(num,den,u); plot(k,y),grid;

xlabel('k');ylabel('y(k)');

运行结果见图。

例 系统脉冲传递函数为

6322

.02644

.03678.0)()(2+-+=z z z z R z C 求离散单位阶跃响应。 MATLAB 命令窗口键入：

num=[0.3678 0.2644];

den=[1 –1 0.6322]; dstep(num,den) 结果见图。

例 设某单位反馈离散系统的开环传递函数为 368

.0368.1264

.0368.0)(2

+-+=z z z z G 运用SIMULINK 进行仿真分析。

在SIMULINK 环境下构建该系统的结构图如图所示：

仿真结果如图所示：

5．稳定性分析

例 设控制系统的特征方程为 0391*******)(23=-+-=z z z z D 运用MATLAB 判断系统的稳定性。 编写MATLAB 程序如下

p=[45 -117 119 -39]; %特征方程系数多项式

r=roots(p); %求特征方程根r x=[-1:0.01:1]'; y=sqrt(1-x.^2);

plot(x,y,x,-y); %绘制单位圆

hold,plot(r,'xr'); %以红色的“×”表示特征方程根 程序运行结果见图：

由图可见，特征方程有两个根在z 平面单位圆外，故系统为不稳定。 四、实验内容

1．根据下面给出的系统结构图，用MATLAB 的语言实现连续系统的离散化。求出系统的开环脉冲传递函数和闭环脉冲传递函数。

（1） T=1

（2）T=0.1

（3）

T=0.5

（4）

T=0.2

(5)T=1

2．根据给出的系统结构图，用MA TLAB 语言或SIMULINK 绘制出系统的单位阶跃响应曲线， 并分析零阶保持器对系统性能的影响。 （1）

（2）

（3）

（4） T=0.1

T=0.5T=0.2s

(5)

（6）

3．根据给出的系统结构图，用MATLAB 的语言实现系统的稳定性判断并分析零阶保持器

对系统稳定性的影响。

，

(1)

(2)

(3) （4）

T=1s

T=1s

T=1s

T=1s

T=0.5T=0.2s