系统仿真综合实验报告

(此文档为word格式，下载后您可任意编辑修改！) 合肥工业大学电气与自动化工程

学院

综合实验报告

实验名称: 系统仿真综合实验

姓名:

学号:

专业班级:

实验地点:

指导教师:

成绩:

日期: 2012年7月

实验一MATLAB基本操作

实验目的

1．熟悉MATLAB实验环境，练习MATLAB命令、m文件、Simulink的基本操作。2．利用MATLAB编写程序进行矩阵运算、图形绘制、数据处理等。

3．利用Simulink建立系统的数学模型并仿真求解。

实验原理

MATLAB环境是一种为数值计算、数据分析和图形显示服务的交互式的环境。MATLAB有3种窗口，即：命令窗口（The Command Window）、m-文件编辑窗口（The Edit Window）和图形窗口（The Figure Window），而Simulink另外又有Simulink模型编辑窗口。

1．命令窗口（The Command Window）

当MATLAB启动后，出现的最大的窗口就是命令窗口。用户可以在提示符“>>”后面输入交互的命令，这些命令就立即被执行。

在MATLAB中，一连串命令可以放置在一个文件中，不必把它们直接在命令窗口内输入。在命令窗口中输入该文件名，这一连串命令就被执行了。因为这样的文件都是以“.m”为后缀，所以称为m-文件。

2．m-文件编辑窗口（The Edit Window）

我们可以用m-文件编辑窗口来产生新的m-文件，或者编辑已经存在的m-文件。在MATLAB主界面上选择菜单“FileNewM-file”就打开了一个新的m-文件编辑窗口；选择菜单“FileOpen”就可以打开一个已经存在的m-文件，并且可以在这个窗

口中编辑这个m-文件。 3．图形窗口（The Figure Window ）

图形窗口用来显示MATLAB 程序产生的图形。图形可以是2维的、3维的数据

图形，也可以是照片等。

MATLAB 中矩阵运算、绘图、数据处理等内容参见教材《自动控制系统计算机仿真》的相关章节。

Simulink 是MATLAB 的一个部件，它为MATLAB 用户提供了一种有效的对反馈控制系统进行建模、仿真和分析的方式。

有两种方式启动Simulink:

1．在Command window 中，键入simulink ，回车。 2．单击工具栏上Simulink 图标。

启动Simulink 后，即打开了Simulink 库浏览器（Simulink library browser ）。在该浏览器的窗口中单击“Create a new model （创建新模型）”图标，这样就打开一个尚未命名的模型窗口。把Simulink 库浏览器中的单元拖拽进入这个模型窗口，构造自己需要的模型。对各个单元部件的参数进行设定，可以双击该单元部件的图标，在弹出的对话框中设置参数。

实验内容

1 用MATLAB 可以识别的格式输入下面两个矩阵

14

436782

33554226753421

8

9543i i B i +????+?

?=??+??

??

再求出它们的乘积矩阵C ，并将C 矩阵的右下角2×3子矩阵赋给D 矩阵。赋值完成后，调用相应的命令查看MATLAB 工作空间的占用情况。

答案：

1）程序：

A=[1 2 3 3;2 3 5 7;1 3 5 7;3 2 3 9;1 8 9 4];B=[1+4i 4 3 6 7 8;2 3 3 5 5 4+2i;2 6+7i 5 3 4 2;1 8 9 5 4 3];

C=A*B

D=C(4:5,4:6)

结果：

C =

1.0e+002 *

Columns 1 through 4

0.1400 + 0.0400i 0.5200 + 0.2100i 0.5100 0.4000

0.2500 + 0.0800i 1.0300 + 0.3500i 1.0300 0.7700

0.2400 + 0.0400i 0.9900 + 0.3500i 1.0000 0.7100

0.2200 + 0.1200i 1.0800 + 0.2100i 1.1100 0.8200

0.3900 + 0.0400i 1.1400 + 0.6300i 1.0800 0.9300

Columns 5 through 6

0.4100 0.3100 + 0.0400i

0.7700 0.5900 + 0.0600i

0.7000 0.5100 + 0.0600i

0.7900 0.6500 + 0.0400i

0.9900 0.7000 + 0.1600i

D =

82.0000 79.0000 65.0000 + 4.0000i

93.0000 99.0000 70.0000 +16.0000i

工作空间:

Name Value

A <5x4 double>

B <4x6 double>

C <5x6 double>

D [82. 000 + 0. 000i,79. 000 + 0. 000i,65. 000 + 4. 000i;93. 000 + 0. 000i,99. 000 + 0. 000i,70. 000 + 16. 000i]

小结：matlab通过确认下标，可以对矩阵进行插入子块、提取子块和重排子块的操作。如果提取子块时，n或m是常数，则返回指定的行列；如果n或m是向量，则返回的是指定矩阵的子块。

2 分别用for和while循环结构编写程序，求出

63

2362630

2122222i i K ===++++++∑

答案： 程序1： s=0;

for k=0:63; s=s+2^k; end

disp('The sum is'),s

结果： The sum is s =

1.8447e+019 程序2：

>> s=0; k=0;

while k<=63; s=s+2^k; k=k+1; end

disp('The sum is'),s 结果：

The sum is s =

1.8447e+019

小结：在for 循环语句，循环体内不能出现对循环控制变量的重新设置，否则会出错；while 循环语句中，在语句内必须有可以修改循环控制变量的命令，否则将陷入死循环，除非循环语句中有控制退出循环的语句。

3 选择合适的步距绘制出下面的图形 （1），其中 （2），其中 答案：

（1）

t=0.1:0.001:10; y=sin(t.^(-1));

plot(t,y)

小结：t=0.1：0.001：10是代表一个步长为0.001的向量，t.^(-1)是代表这个向量的点运算，而不是代表这个向量的向量运算，运算才不会出错。（2）

t=linspace(-pi,pi,60);

y=sin(tan(t))-tan(sin(t));

plot(t,y)

小结：linspace函数能够生成线性分度的向量。

4 对下面给出的各个矩阵求取矩阵的行列式、秩、特征多项式、范数。，

，

答案：

程序：

A=[7.5 3.5 0 0;8 33 4.1 0;0 9 103 -1.5;0 0 3.7 19.3]; B=[5 7 6 5;7 10 8 7;6 8 10 9;5 7 9 10];

C=[1 2 3 4;5 6 7 8;9 10 11 12;13 14 15 16];

D=[3 -3 -2 4;5 -5 1 8;11 8 5 -7;5 -1 -3 -1];

detA=det(A),rankA=rank(A),polyA=poly(A),normA=norm(A) detB=det(B),rankB=rank(B),polyB=poly(B),normB=norm(B) detC=det(C),rankC=rank(C),polyC=poly(C),normC=norm(C) detD=det(D),rankD=rank(D),polyD=poly(D),normD=norm(D) 结果：

detA =

4.3222e+005

rankA =

4

polyA =

1.0e+005 *

0.0000 -0.0016 0.0713 -1.0659 4.3222

normA =

103.7228

detB =

1.0000

rankB =

4

polyB =

1.0000 -35.0000 146.0000 -100.0000 1.0000 normB =

30.2887

detC =

4.7332e-030

rankC =

2

polyC =

1.0000 -34.0000 -80.0000 -0.0000 -0.0000

normC =

38.6227

detD =

595.0000

rankD =

4

polyD =

1.0000 -

2.0000 -32.0000 532.0000 595.0000

normD =

16.6958

小结：det、rank、poly、norm分别求出矩阵的行列式、秩、特征多项式以及范数

5 求解下面的线性代数方程，并验证得出的解真正满足原方程。

(a)

72124

915327

221151

132130

X

-

????

????

-

????

=

????

---

????

????

，(b)

1321390

721264

91532117

2211521

X

????

????

-

????

=

????

-

????

----

????

答案：

(a).

程序：

A=[7 2 1 -2;9 15 3 -2;-2 -2 11 5;1 3 2 13]; B=[4;7;-1;0];

X=inv(A)*B

结果：

X =

0.4979

0.1445

0.0629

-0.0813

(b)

A=[1 3 2 13;7 2 1 -2;7 15 3 -2; -2 -2 11 5]; B=[9 0;6 4;11 7;-2 -1];

X=inv(A)*B

结果

X =

0.9247 0.4695 0.4140 0.2186 -0.1907 0.0791

0.5550 -0.0987

小结：求矩阵A*B=C,已知A,C 求B ，则应用逆矩阵来求解即B=(A*(-1))*C

6 假设有一组实测数据 x 0.1 0.2 0.3 0.4 0.5 0.6 0.

7 0.

8 0.

9 1 y

2.320

1

2.647

2.970

7

3.288

5

3.600

8

3.909

4.214

7

4.519

1

4.823

2

5.127

5

用最小二乘法拟合，求出相应的二次函数。 答案：

程序

x=0.1:0.1:1;

y=[2.3201 2.6470 2.9707 3.2885 3.6008 3.9090 4.2147 4.5191 4.8232 5.1275]; p=polyfit(x,y,2) 结果： p =

-0.1563 3.2828 1.9967

小结：polyfit 是最小二乘拟合的函数，本题说明用最小二乘拟合求得的二次函数为y=-0.1563*x^2+3.2828*x+1.9967

7 考虑线性微分方程

(4)

(3)

353345sin(4/3)t t y

y

y y y u

u e e t π???

--++++==++输入信号

(1)(2)(3)(0)1,(0)(0)1/2,0.2,y y y y ====方程初值

(1) 试用Simulink 搭建起系统的仿真模型，并绘制出仿真结果曲线。

(2) 将给定的微分方程转换成状态方程，并建立S 函数，再利用Simulink 进行仿真。 答案： (1)

参数设置：

Integrator ：Initial condition 分别设置成0.2，0.5，0.5，1 Gain0,1,2, 3：gain 分别设置成3，3，4，5

Sine Wave ： Frequency(radsec)设置成 4 ；phase(rad)设置成 pi3

Transfer Fcn ：Numerator coefficients设置成[1] ；Denominator coefficients 设置成[1 5]

Transfer Fcn1：Numerator coefficients设置成[1] ；Numerator coefficients设置成[1 3]

Add ：List of signs 设置成++----

仿真模型：

（2）

建立相应的S函数：

function[sys,x0,str,ts]= reds7_2 (t,x,u,flag,A,B,C,D) switch flag

case 0

[sys,x0,str,ts]=mdlInitializeSizes(A,D);

case 1

sys=mdlDerivatives(t,x,u,A,B);

case 3

sys=mdlOutputs(t,x,u,C,D);

case {2,4,9}

sys=[];

otherwise

error(['Unhandled flag=',num2str(flag)]); end

function[sys,x0,str,ts]=mdlInitializeSizes(A,D) sizes=simsizes;

sizes.NumContStates=size(A,1);

sizes.NumDiscStates=0;

sizes.NumOutputs=size(A,1)+size(D,1);

sizes.NumInputs=size(D,2);

sizes.DirFeedthrough=1;

sizes.NumSampleTimes=1;

sys=simsizes(sizes);

x0=[1;0.5;0.5;0.2];

str=[];

ts=[-1 0];

function sys=mdlDerivatives(t,x,u,A,B)

sys=A*x+B*u;

function sys=mdlOutputs(t,x,u,C,D)

sys=[C*x+D*u;x];

由此建立仿真模型

Scope

reds7_2S-Function

u y fcn Embedded MATLAB Function

Clock

仿真曲线：

小结：可以根据微分方程做出状态图并赋初值，得到系统的仿真模型以及仿真结果曲线；通过构建S 函数，可以将状态方程封装在S 函数，以简化运算。

8 建立下图所示非线性系统的Simulink 模型，并观察在单位阶跃信号输入下系统的输出曲线和误差曲线。

答案：参数设置：

Relay: Switch on point设置成0.5；Switch on point设置成0.5；Output when on 设置成2.5；Output when off 设置成0

Relay1: Switch on point 设置成-0.5；Switch on point设置成-0.5 ；Output when on 设置成2.5；Output when off 设置成0

Transport Delay: Time delay设置成0.4

Transfer Fcn: Numerator coefficients设置成[30];Denominator coefficients设置成[1 0.5 1]

Transfer Fcn1: Numerator coefficients设置成[30];Denominator coefficients设置成[1 6.5 1]

仿真模型：

仿真曲线

小结：1.用两个继电器并联得到具有死去特性的三位置继电特性

2.从图中可以看出输出比较稳定

实验二经典控制理论

实验目的

以MATLAB及Simulink为工具，对控制系统进行时域、频域及根轨迹分析。

实验原理

1、时域分析法是根据系统的微分方程（或传递函数），利用拉普拉斯变换直接解出

动态方

程，并依据过程曲线及表达式分析系统的性能。时域响应指标如图所示。

延迟时间t d，指响应曲线第一次达到其终值一半所需要的时间。

上升时间t r，指响应曲线从终值10%上升到终值90%所需要的时间；对于有振荡的系统，也可定义为响应从零第一次上升到终值所需要的时间。上升时间是系统响应速度的一种度量。

峰值时间t p，指响应超过终值达到第一个峰值所需要的时间。

调节时间t s，指响应达到并保持在终值±5%（或±2%）内所需要的时间。

超调量σ%，指响应的最大偏离量h(t p)与终值h(∞)之差的百分比，即：

稳态误差，描述系统稳态性能的一种性能指标。

2、频域分析法通常从频率特性出发对系统进行研究。在工程分析和设计中，通常把

频率特

性画成一些曲线，从频率特性曲线出发进行研究。这些曲线包括幅频特性和相频特性曲线，幅相频率特性曲线，对数频率特性曲线以及对数幅相曲线等，其中以幅相频率特性曲线，对数频率特性曲线应用最广。对于最小相位系统，幅频特性和相频特性之间存在着唯一的对于关系，故根据对数幅频特性，可以唯一地确定相应的相频特性和传递函数。根据系统的开环频率特性去判断闭环系统的性能，并能较方便地分析系统参量对系统性能的影响，从而指出改善系统性能的途径。

3、根轨迹是求解闭环系统特征根的图解方法。由于控制系统的动态性能是由系统闭

环零极点共同决定，控制系统的稳定性由闭环系统极点唯一确定，利用根轨迹确定闭环系统的零极点在s平面的位置，分析控制系统的动态性能。

实验内容

*1．控制系统数学模型的转换

《自动控制系统计算机仿真》教材第4章中的所有例题

【例4-1】已知系统的传递函数如下，利用MATLAB建立其相应的传递函数系统模型

答案

程序：

clc,clear

num=5*[2 0 3];

den=conv(conv(conv([1 0 0],[3 1]),conv([1 2],[1 2])),[5 0 3 8]); printsys(num,den,'s')

G=tf(num,den)

结果：

小结：G=tf(num,den)是用来求传递函数的系统模型的函数

【例4-2】已知系统传递函数如下

应用MATLAB函数将其转换为状态方程形式的模型答案：

程序：

clc

clear

num=[12 24 0 20];

den=[2 4 6 2 2];

sys=tf(num,den)

s=ss(sys)

结果：

Transfer function:

12 s^3 + 24 s^2 + 20

2 s^4 + 4 s^

3 + 6 s^2 + 2 s + 2

a =

x1 x2 x3 x4

x1 -2 -1.5 -0.5 -0.5

x2 2 0 0 0

x3 0 1 0 0

x4 0 0 1 0

b =

u1

x1 4

x2 0

x3 0

x4 0

c =

x1 x2 x3 x4

y1 1.5 1.5 0 1.25

u1

y1 0

Continuous-time model.

小结：tf(num,den)函数用来求传递函数的状态模型

【例4-3】某线性定常系统的状态空间表达式如下，求该系统的传递函数。

答案：

clc

clear

A=[0 1 1;0 0 1;-10 -17 -8];

B=[0;0;1];C=[5 6 1];D=0;

SS=ss(A,B,C,D)

G1=tf(SS)

结果：

Transfer function:

s^2 + 11 s + 5

s^3 + 8 s^2 + 27 s + 10

小结：ss(A,B,C,D)函数是用来求状态方程的传递函数

【例4-4】对于例4-3题中的线性定常系统，将其转换为zpk形式

答案：

程序：

clc

clear

A=[0 1 1;0 0 1;-10 -17 -8];

B=[0;0;1];C=[5 6 1];D=0;

[z,p,k]=ss2zp(A,B,C,D)

G1=zpk(z,p,k)

z =

-10.5249

-0.4751

p =

-0.4199

-3.7901 + 3.0745i

-3.7901 - 3.0745i

k =

1.0000

Zeropolegain:

(s+10.52) (s+0.4751)

(s+0.4199) (s^2 + 7.58s + 23.82)

小结：zpk(z,p,k)是用来求状态方程的zpk形式的函数【例4-5】某系统的传递函数如下，求它的部分分式形式。

答案：

程序：

clc

clear

num=[20 10];den=[1 15 74 120];

[R,P,H]=residue(num,den)

结果

R =

-55.0000

90.0000

-35.0000

P =

-6.0000

-5.0000

-4.0000

H =

[]

小结：[R,P,H]=residue(num,den)函数是用来求传递函数的部分分式函数【例4-6】某连续系统的状态空间表达式如下，采用零阶保持器将其离散化，设采样周期为0.1秒，求离散化的系统方程。

答案：

程序：

A=[0 1 0;0 0 1;-6 -11 -6];

B=[1 0;2 -1;0 2];

C=[1 -1 0;2 1 -1];

D=0;

T=0.1;

G=ss(A,B,C,D);

Gn=c2d(G,T)

结果

a =

x1 x2 x3

x1 0.9991 0.0984 0.004097

x2 -0.02458 0.9541 0.07382

x3 -0.4429 -0.8366 0.5112

b =

u1 u2

x1 0.1099 -0.004672

x2 0.1959 -0.0902

x3 -0.1164 0.1936

c =

x1 x2 x3

y1 1 -1 0

y2 2 1 -1

d =

u1 u2

y1 0 0