计算机仿真实验报告 MATLAB

注：以下所有程序均在MATLAB7.0下运行通过。

实验一MATLAB语言编程

一、实验目的：

熟悉MATLAB语言及其环境，掌握编程方法。

要求认真听取实验指导老师讲解与演示。

二、具体实验内容、步骤、要求：

1.运行交互式学习软件，学习MATLAB语言

2.在MATLAB的命令窗口下键入如下命令：

INTRO

(注意：intro为一个用MATLAB语言编写的幻灯片程序，主要演示常用的MATLAB语句运行的结果。)

然后，根据显示出来的幻灯片右图按钮进行操作，可按START→NEXT→NEXT按钮一步步运行，观察。

3.自编程序并完成上机编辑，调试，运行，存盘：

（1）、用MATLAB命令完成矩阵的各种运算。

例如：

求出下列运算结果，并上机验证。

（1）A(:,1) %取矩阵A的第一列元素

ans =

11

21

31

41

（2）A(2,:) %取矩阵A的第二行元素

ans =

21 22 23 24

（3）A(1:2,2:3) %取矩阵A第一二行第二三列的元素

ans =

12 13

22 23

（4）A(2:3,2:3) %取矩阵A第二三行第二三列的元素

ans =

22 23

32 33

（5）A(:,1:2) %取矩阵A第一列与第二列元素

ans =

11 12

21 22

31 32

41 42

（6）A(2:3) %取矩阵A第二行与第三行的首列元素

ans =

21 31

（7）A(:) %将矩阵A的所有元素按一列排列

ans =

11

21

31

41

12

22

32

42

13

23

33

43

14

24

34

44

（8）A(:,:) %显示矩阵A

ans =

11 12 13 14

21 22 23 24

31 32 33 34

41 42 43 44

（9）ones(2,2) %建立一个两行两列的全1矩阵

ans =

1 1

1 1

（10）eye(2) %建立一个二维的单位矩阵

ans =

1 0

0 1

（2）、绘制数学函数的图形：

例如：，理解数组运算与矩阵运算的功能。MATLAB程序如下：

t=0:0.1:8; %建立向量t

y=1-2*exp(-t.*sin(t)); %计算向量t的函数向量y

plot(t,y); %利用plot命令绘图

xlabel('t');ylabel('y=1-2*e^(-t*sin(t))'); %注释横坐标与纵坐标

图1.1 所对应的函数图像

4.理解函数文件与命令文件的区别，并自编函数文件并调用。

命令文件可以理解为简单的M文件，命令文件中的变量都是全局变量。函数文件是在命令文件的基础之上多添加了一行函数定义行，其代码组织结构和调用方式与对应的命令文件截然不同。函数文件是以函数声明行“function...”作为开始的，其实质就是用户往MATLAB函数库里边添加了子函数，函数文件中的变量都是局部变量，除非使用了特别声明。函数运行完毕之后，其定义的变量将从工作区间中清除。而命令文件只是将一系列相关的代码结合封装，没有输入参数和输出参数，即不自带参数，也不一定要返回结果。而多数函数文件一般都有输入和输出变量，并见有返回结果。

举例如下：已知非线性方程组

求系统在初始条件为的数值解。

1）定义方程组的m函数equ.m

function xs=equ(xi)

x=xi(1);y=xi(2);z=xi(3);

xs=zeros(3,1);

xs(1)=sin(x)+y^2+log(z)-7;

xs(2)=3*x+2^y-z^3+1;

xs(3)=x+y+z-5;

2）方程组求解

x0=[1 1 1];

xyz=fsolve(@equ,x0)

xyz =

0.5991 2.3959 2.0050

在调用fsolve命令求解方程时调用了自定义的equ函数。

5.学会通过Help，熟悉MATLAB中为用户提供的功能各异的函数与命令。

只需在MATLAB命令窗口中键入help 命令名或函数名，按回车键后即会显

示该命令或函数的相关说明。根据说明即可自学如何使用该命令或函数。

例如键入help eig

实验二数值积分算法练习与函数调用

一、实验目的：

理解数值积分法，熟练掌握MATLAB的函数调用。

二、实验示例介绍：

1.用Euler法求初值问题的数值解：

设方程如下：

取步长h=0.1，名为FZSYZ1.M

上机用如下程序名FZSYZ1.M可求出数值解。

2.在MATLAB中提出了现成的数值积分函数，如ode1、ode23、ode45求解微分方程组，下面介绍ode23它为二阶/三阶的RKF方法在MATLAB的ToolBox。

文件夹中的MATLAB/funfun下的M文件，在此介绍其调用方法与应用示例如下：

[t,x]=ode23(‘系统函数名’，)

其中，系统函数名为描述系统状态方程的M函数的名称，该函数名在调用时应

该用引号括起来（文件名字串），为起始与终止时间，为系统初始变量的

值（列向量），为控制解的精度（缺省值为=在ode23中），trace为输

出形式控制变量，非零则程序运行的每步都显示出来。trace为缺省值时——不显示中间结果。t为输出参数返回积分的时间离散值（列向量）。X：输出参量，返回每个时间点值的解的列向量。

注意：系统函数的编写格式为固定的。

例如：若MATLAB求解初值问题的解，其方程如下：

第一步：编写如下程序：

function xdot=fun21(t,x)

xdot=x-t^2

并以fun21.m存盘。

第二步：编写如下程序并以fzsy22.m存盘。

t0=0;tf=3;tol=1e-6;

x0=1;trace=1;

[t,x]=ode23('fun21',t0,tf,tol,trace)

plot(t,x)

第三步：在命令窗口下运行fzsy21即可求出x的解，并画出曲线。

3.实验具体内容、步骤、要求：

（1）运行交互式软件中函数调用，学习程序，

（2）试将（2-2）方程改为用Euler编程求解，试比较用ode23求解结果。MATLAB源程序如下：

%ode23求解

t0=0;tf=3;tol=1e-6;

x0=1;trace=1;

[t,x]=ode23('fun21',t0,tf,tol,trace);

plot(t,x,'r');hold on;

%Euler求解

h1=0.1;

t=[t0:h1:tf];

n=length(t);

u=x0;

uu(1)=u;

for i=2:n

du=u-t(i-1)^2;

u=du*h1+u;

uu(i)=u;

end

plot(t,uu,'b--');hold off;

图2.1 式（2-2）用Euler编程与ode23算法求解结果

（3）试将（2-1）方程改为用ode23算法调用求解，并试比较结果。

图2.2 式（2-1）用Euler编程与ode23算法求解结果

从比较可知：Euler需要小步长，不能自动确定步长；ode23算法自动确定步长，如果误差小于指定的精度，进行下一步，否则缩小步长，因此ode23算法比Euler算法更为准确，但占用的计算时间更长。

（4）利用ode23或ode45求解线性时不变系统微分方程

并绘制出的曲线。式中

MATLAB程序清单如下：

%函数(y(t) ) ?=Ay(t)

function ys=fun22(t,y)

ys=zeros(2,1);

ys(1)=-0.5*y(1)+y(2);

ys(2)=-y(1)-0.5*y(2);

%利用ode45求解其数值解并作图

t0=0;tf=4;tol=1e-6;y0=[0;1];trace=1;

y=zeros(2,1);

[t,y]=ode45('fun22',t0,tf,y0,tol,trace);

plot(t,y(:,1),'r');hold on;

plot(t,y(:,2),'b');hold off；

图2.3 第(4)问利用ode45算法求解所得y(t)曲线

（5）求出与的单位阶跃

响应，并分别求出状态空间模型。

MATLAB程序清单如下：

num1=[2];den1=[1 2 1];

g1=tf(num1,den1);%构造连续传函G1(S)

num2=[1];den2=[2 3 3 1];

g2=tf(num2,den2);%构造连续传函G2(S)

step(g1);hold on;%画出G1(S)的阶跃响应

step(g2);hold off;%画出G2(S)的阶跃响应

[A1,B1,C1,D1]=tf2ss(num1,den1)%得到G1(S)的状态空间描述

[A2,B2,C2,D2]=tf2ss(num2,den2)%得到G2(S)的状态空间描述

从MATLAB窗口中得到的结果可知，对应的状态空间模型为：

对应的状态空间模型为：

图2.4 传递函数G1(S)与G2(S)对应的单位阶跃响应

选做题一：已知系统传递函数为求对应的零极点模

型，绘制系统阶跃响应。

MATLAB程序清单如下：

num=[200 400];den=conv([1 1],[1 10 42]);

g=tf(num,den);%构造连续传递函数

[z,p,k]=tf2zp(num,den)%将普通传函转化为零极点模型

step(g);

在MATLAB命令窗中显示如下结果：

z =

-2

p =

-5.0000 + 4.1231i

-5.0000 - 4.1231i

-1.0000

k =

200

所以对应的零极点模型为。

图2.5 选做题一中的传递函数对应的单位响应曲线图

选做题二：已知一个二阶系统为，分

别取R=0与R=0.5时，系统模型取不同步长的数值解（h=0.5~0.001）。MATLAB程序清单如下：

function ys=cstep(t,y)

R=0;

ys=zeros(2,1);

ys(1)=y(2);

ys(2)=-2*R*y(2)-y(1);

[t,y]=ode45(@cstep,[0:0.5:5],[100 0]);

plot(t,y(:,1),'r--');hold on;

tol=1e-6;trace=1;t0=0;tf=5;

[t1,y1]=ode45('cstep',t0,tf,[100 0],tol,trace);

plot(t1,y1(:,1),'b');hold off

在MATLAB命令窗中显示如下结果：

图2.6 R=0时的仿真图2.7 R=0.5时的仿真从图2.6与图2.7可以看出，步长越短，仿真结果越接近准确值。当步长取0.5时，可以看到与步长取0.001时有明显偏离，而采取自动步长的ode45算法与步长取0.001的仿真

结果一致。

（6）设方程式为用各种数值积分法与不同步长求方程式的数

值解，并比较之。

MATLAB程序清单如下：

%欧拉法求解子函数

function [uu,t]=euler(t,h1,y0)

n=length(t);

u=y0;

uu(1)=u;

for i=2:n

du=-40*u;

u=du*h1+u;

uu(i)=u;

end

%主程序

t0=0;tf=0.5;y0=2;h1=0.01;

t1=[0:h1:1];

[uu1,t1]=euler(t,h1,y0);%步长取0.01，用欧拉法求解

plot(t1,uu1,'r');hold on;

h2=0.001;%步长取0.01，用欧拉法求解

t2=[0:h2:1];

[uu2,t2]=euler(t2,h2,y0);

plot(t2,uu2,'b');hold on;

tol=1e-6;trace=1;u0=1;

[t,u]=ode23('fun26',t0,tf,y0,tol,trace);%用ode23求解，步长自定

plot(t,u,'g--');hold off;grid

gtext('Euler0.01');gtext('Euler0.001');gtext('ode23');

在MATLAB命令窗中显示如下结果：

图2.8 采用不同算法对第（6）题微分方程求解结果示意图从图2.5可以看出，采用欧拉法求解时，步长选取越短，其结果越接近真实值，选用龙阶-库塔法时，由于程序自动变步长直到达到预设精度，因此利用该

方法求得的结果较准确。当欧拉法的步长选为0.001时，其计算结果与采用ode23算法结果一致。如果采用ode45算法计算，其结果与ode23相差不大。

实验三控制工具箱与SIMULINK软件使用

一、实验目的：

熟悉工具箱及其使用，进行系统仿真分析，通过仿真对系统进行校正校验。

二、实验预先要求：

必须先复习教材及上课介绍的有关控制工具箱命令与SIMULINK仿真工具的使用，并对实验题目做好准备。

三、学会调出、运行已有SIMULINK建立的仿真模型。

在打开MATLAB的窗口下，

（1）输入SIMULINK或双击命令窗口工具栏右起第三个按钮（NEW SIMULINK MODEL）将会打开，Library SIMULINK然后指向FILE菜单下拉菜单open调出fzsy31.mdl文件。然后在fzsy31.mdl文件的菜单观察并记录有关设置参数，然后指向Start下拉菜单单击一次观察输出波形。该仿真模型（SIMULINK系统结构框图）如下：

在MATLAB主窗口的Workspace中找到y变量，然后双击，可以看到tout=0：0.1：10对于的y输出，示波器绘图正是根据tout数组与y数组绘制。

图3.1 fzsy31.ml仿真所得波形

从图3.1可以看出，加限幅环节后系统的响应变慢，但超调量变小。实际系

统由于输入信号不能太大，避免冲击，往往需要加限幅环节，但一般加至控制器的输出端，使控制量输出作用于被控对象没有较大的冲击。

四、实验设计题目与要求：

（1）参见题目要求，建立仿真模型，进行仿真实验与分析。然后，回答题目提出问题。

本习题介绍了滞后于超前补偿器在不稳定系统中的补偿作用。

以双积分系统为例

其补偿器具有如下三种：

Ⅰ、比例补偿器：

Ⅱ、超前补偿器：

Ⅲ、滞后补偿器：

初步分析：绘制系统根轨迹图，并回答下列问题。

①在第Ⅰ种情况下，系统能够稳定吗？增加的作用是什么？

绘制根轨迹：

g1=tf([1],[1 0 0]);rlocus(g1)

图3.2 比例补偿器对应系统的根轨迹图

从图3.2可以看出，串联比例补偿器后，不论开环增益如何变化，闭环特征根总在虚轴上，系统总处于临界稳定。改变的值意义不大，因为系统本身不稳定。

②在第Ⅰ种情况下，系统能够稳定吗？增加的作用是什么？

g2=tf([1 1],[1 5 0 0]);

rlocus(g2)

图3.3 超前补偿器对应系统的根轨迹图

从图3.3可以看出，加超前补偿器后，系统的所有闭环特征根总位于s平面的右半平面，系统稳定，实现了补偿的目的，增大能加快系统的响应速度。

③在第Ⅰ种情况下，系统能够稳定吗？增加的作用是什么？

g3=tf([1 5],[1 1 0 0]);

rlocus(g3)

图3.3 滞后补偿器对应系统的根轨迹图

从图3.3可以看出，串联滞后补偿器后，系统的闭环特征根中存在位于s平面右半平面的根，系统发散，改变的值意义不大，因为系统本身不稳定。

上机实验：

④设={0.1，0.5，1}，对上述每一种情况进行仿真，求得其阶跃响应。

要求将这三种情况的闭环阶跃响应绘制在一张图中。分析该图，证明在初步分析过程中所得的结论。在每种情况中，仿真时间可为{10，50，5}。

Ⅰ、比例补偿器：

图3.4 串联比例补偿器SIMULINK仿真模型

图3.5 串联比例控制器后系统的阶跃响应仿真波形

从上述的波形可以看出，单纯用比例控制器无法使系统稳定，系统总处于临界稳定状态，即等幅振荡，这与根轨迹在虚轴上所得结论一致。增大比例系数仅增大等幅振荡的频率。

Ⅱ、超前补偿器：

图3.6 串联超前补偿器SIMULINK仿真模型

图3.7 串联超前补偿器后系统的阶跃响应仿真波形

从上图可以看出，加大的值能加快系统的响应速度。由于该系统为Ⅱ型系

统，所以能无差跟踪阶跃输入。系统型别判断稳态误差仅在系统稳定的前提下适应，后面的滞后补偿器对应的系统由于本身不稳定，故不能利用此结论。

Ⅲ、滞后补偿器：

SIMULINK模型与串联超前补偿器一样，只需将补偿环节的相关参数加以修改。

从图3.8可以看出，不论值如何变化，系统总不稳定，越大，系统发散

得越厉害，与根轨迹分析结果一致。

图3.8 串联滞后补偿器后系统的阶跃响应仿真波形

（2）系统对象为

PI补偿器

a 若由与串联组成的单位反馈系统，称为系统A。

b 若由放置在的反馈通道称为系统B。

c 若由与串联组成的单位称为系统C。绘制上述每一种情况的阶跃响

应与斜坡响应。求每一种情况下的系统阶跃响应误差与斜坡响应误差，并与理论分析结果进行比较。

a MATLAB程序如下：

g=tf([4],[1 2 0]);%构造G1(S)

k1=tf([1 1],[1 0]);%构造K1(S)

gk1=series(g,k1);%将G1与K1串联起来

sys1=feedback(gk1,1);%构造单位负反馈传递函数

t=0:0.1:10;

[y1,t]=step(sys1);%返回sys1的单位阶跃响应的时间矩阵与对应输出矩阵

figure(1)

plot(t,y1);%绘制sys1单位阶跃响应曲线

u=t;

figure(2)

[y2,t2]=lsim(sys1,u,t,0);%返回sys的斜坡响应的时间矩阵与对应输出矩阵

plot(t2,y2);%绘制sys1斜坡响应曲线

i=length(t);

y1(i)-1%求单位阶跃响应的稳态误差

y2(length(t2))-t2(length(t2))%求斜坡响应的稳态误差

图3.9 系统A单位阶跃响应曲线图3.10 系统A斜坡响应曲线在MATLAB窗口中得到的单位阶跃响应稳态误差是0.0026，斜坡响应的稳态误差是-0.0024，仿真的时间仅为16s，如果在更长的时间内，系统的稳态误差将为零。从理论上分析，该系统为Ⅱ型系统，能无差跟踪单位阶跃信号与斜坡信号，所以理论分析结果与实际仿真结果一致。

b MATLAB程序如下：

与a中程序大致相同，只需把feedback函数改为B=feedback(g,k1);

图3.11 系统B单位阶跃响应曲线图3.12 系统B斜坡输入响应曲线在MATLAB窗口中得到的稳态误差分别为-0.093与-4.0586，如果增大仿真时间，系统B对于单位阶跃输入的问题误差将为0，对斜坡输入的稳态误差为无穷大，也就是说系统B不能跟踪斜坡输入。

C MATLAB程序如下：

与a中程序大致相同，只需把k2=tf([1 1],[1 0]);函数改为k2=tf([1 0.5],[1 0]);

图3.13 系统C单位阶跃响应曲线图3.14 系统C斜坡输入响应曲线在MATLAB窗口中得到的稳态误差分别为0.0015与-0.0035，如果增大仿真时间，系统C对于单位阶跃输入与斜坡输入的稳态误差均为0。但比较图3.9与图3.13，图3.10与图3.14可知，后者的响应速度明显比前者快，这是因为系统C 相当于把PI调节器的积分系数改小，从1改为0.5，减小积分时间常数，加快了系统动态响应过程，因此后者响应速度加快。

（3）完成以下题目：（本习题是PID控制系统的设计问题），考虑如下一个不稳定系统。

设计要求：闭合稳定，且单位阶跃响应稳态误差为0.

初步分析：

①利用根轨迹分析可知，使用串联PID补偿器可以达到上述两个方面的

要求。注意在考虑系统稳态误差之前，系统必须是稳定的。设系统为单位负反馈系统结构。PID补偿传递函数为

图3.15 串联PID补偿器后系统的根轨迹图

从图3.15可以看出，当开环增益大于2.73时，系统的所有闭环特征根在s 平面右半平面，系统稳定，所以串联PID补偿器能使不稳定系统G(S)稳定。

上机实验：

②使用试探方法选择PID参数，满足上述特性要求，然后求其系统的阶跃响应，会绘制其响应图形。注意到PID补偿引入了两个零点。这些零点可以为复数或实数，请分别对两种情况进行实验。

MATLAB原程序如下：

kp=50;kd=30;ki=100;

g=tf([1],[1 -2 1]);

pid=tf([kd kp ki],[1 0]);

sys1=series(g,pid);

sys=feedback(sys1,1);%求闭环传递函数

t=0:0.1:10;

[y,t1]=step(sys,t);

plot(t1,y);grid

[num,den]=tfdata(sys);

roots(num{1})%求所引入的另个闭环零点

要使补偿后闭环零点为实数只需满足。

当=50，=100，=30时，闭合零点为-0.833+1.6245i与-0.833-1.6245i，所得的仿真波形如下：

图3.16 PID补偿器器零点为复数仿真图图3.17 PID补偿器零点为实数仿真图

当=30，=10，=10时，闭合零点为-0.3820与-2.6180，所得的仿真波形如图3.17所示。

实验四数字控制系统仿真与综合应用

一、实验目的：

掌握用于离散系统的MATLAB命令及数字控制系统仿真方法。应用仿真对系统的设计进行分析比较。

二、实验内容：（分三个部分）

1.实验题目与要求：

若，若使用零阶保持器将G(S)转换为G(Z)。采样周期T=0.1S，

用仿真方法得出G(S)与G(Z)的阶跃响应图并进行比较。

将连续系统的传递函数用下面零阶保持器法转换成离散系统传递函数。

num=10;den=[1 7 10];ts=0.1;

[n_zoh,d_zoh]=c2dm(num,den,ts)

上述程序可以给出如下结果：

n_zoh=

0 0.0398 0.0315

d_zoh=

1.000 -1.4253 0.4966

下面求取G(S)和G(Z)的阶跃响应，首先将G(S)离散化，然后用step命令求取连续系统的阶跃响应，用dstep命令求取离散系统的阶跃响应，并用hold命令将其画在同一坐标内。

连续系统和用零阶保持器发离散后数字系统的阶跃响应是相同的。因此这种数字化方法总是产生与连续系统相同的阶跃响应。

2.仿真研究与系统分析

参见《计算机控制系统》教材P160例5-1。被控对象为

①采用零阶保持器，将W1(S)转换为W1(Z)，串接一个计算机调节模型D(Z)组成一个典型的数字反馈系统（见图4.1）。当采用最少拍设计原则设计出来的。（对于单位阶跃输入下）