合肥工业大学电气与自动化工程学院系统仿真综合实验报告.doc

《合肥工业大学电气与自动化工程学院系统仿真综合实验报告.doc》由会员分享，可在线阅读，更多相关《合肥工业大学电气与自动化工程学院系统仿真综合实验报告.doc（70页珍藏版）》请在三一文库上搜索。

1、合肥工业大学电气与自动化工程学院综合实验报告实 验 名 称:系统仿真综合实验姓 名:学 号:专 业 班 级:实 验 地 点:指 导 教 师:成 绩:日 期:2012年7月实验一 MATLAB基本操作实验目的1熟悉MATLAB实验环境，练习MATLAB命令、m文件、Simulink的基本操作。2利用MATLAB编写程序进行矩阵运算、图形绘制、数据处理等。3利用Simulink建立系统的数学模型并仿真求解。实验原理MATLAB环境是一种为数值计算、数据分析和图形显示服务的交互式的环境。MATLAB有3种窗口，即：命令窗口（The Command Window）、m-文件编辑窗口（The Edit

2、Window）和图形窗口（The Figure Window），而Simulink另外又有Simulink模型编辑窗口。1命令窗口（The Command Window）当MATLAB启动后，出现的最大的窗口就是命令窗口。用户可以在提示符“”后面输入交互的命令，这些命令就立即被执行。在MATLAB中，一连串命令可以放置在一个文件中，不必把它们直接在命令窗口内输入。在命令窗口中输入该文件名，这一连串命令就被执行了。因为这样的文件都是以“.m”为后缀，所以称为m-文件。2m-文件编辑窗口（The Edit Window）我们可以用m-文件编辑窗口来产生新的m-文件，或者编辑已经存在的m-文件。在M

3、ATLAB主界面上选择菜单“File/New/M-file”就打开了一个新的m-文件编辑窗口；选择菜单“File/Open”就可以打开一个已经存在的m-文件，并且可以在这个窗口中编辑这个m-文件。3图形窗口（The Figure Window）图形窗口用来显示MATLAB程序产生的图形。图形可以是2维的、3维的数据图形，也可以是照片等。 MATLAB中矩阵运算、绘图、数据处理等内容参见教材自动控制系统计算机仿真的相关章节。Simulink是MATLAB的一个部件，它为MATLAB用户提供了一种有效的对反馈控制系统进行建模、仿真和分析的方式。有两种方式启动Simulink:1在Command w

4、indow中，键入simulink，回车。2单击工具栏上Simulink图标。启动Simulink后，即打开了Simulink库浏览器（Simulink library browser）。在该浏览器的窗口中单击“Create a new model（创建新模型）”图标，这样就打开一个尚未命名的模型窗口。把Simulink库浏览器中的单元拖拽进入这个模型窗口，构造自己需要的模型。对各个单元部件的参数进行设定，可以双击该单元部件的图标，在弹出的对话框中设置参数。 实验内容1 用MATLAB可以识别的格式输入下面两个矩阵 再求出它们的乘积矩阵C，并将C矩阵的右下角23子矩阵赋给D矩阵。赋值完成后，调

5、用相应的命令查看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*BD=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.0400

6、i 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.1600iD = 82.0000 79.0000 65.0000 + 4.

7、0000i 93.0000 99.0000 70.0000 +16.0000i工作空间:Name Value A B C D 82.000000000000000 + 0.000000000000000i,79.000000000000000 + 0.000000000000000i,65.000000000000000 + 4.000000000000000i;93.000000000000000 + 0.000000000000000i,99.000000000000000 + 0.000000000000000i,70.000000000000000 + 16.0000000000000

8、00i小结：matlab通过确认下标，可以对矩阵进行插入子块、提取子块和重排子块的操作。如果提取子块时，n或m是常数，则返回指定的行列；如果n或m是向量，则返回的是指定矩阵的子块。2 分别用for和while循环结构编写程序，求出答案：程序1：s=0;for k=0:63; s=s+2k;enddisp(The sum is),s结果：The sum iss = 1.8447e+019程序2： s=0;k=0;while k=63; s=s+2k; k=k+1;enddisp(The sum is),s结果：The sum iss = 1.8447e+019小结：在for循环语句，循环体内不能

9、出现对循环控制变量的重新设置，否则会出错；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函数能够

10、生成线性分度的向量。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),p

11、olyB=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+005rankA = 4polyA = 1.0e+005 * 0.0000 -0.0016 0.0713 -1.0659 4.3222normA = 103.7228detB = 1.0000rankB = 4polyB = 1.0000 -35.0000 146.0000 -100.0000 1.

12、0000normB = 30.2887detC = 4.7332e-030rankC = 2polyC = 1.0000 -34.0000 -80.0000 -0.0000 -0.0000normC = 38.6227detD = 595.0000rankD = 4polyD = 1.0000 -2.0000 -32.0000 532.0000 595.0000normD = 16.6958小结：det、rank、poly、norm分别求出矩阵的行列式、秩、特征多项式以及范数5 求解下面的线性代数方程，并验证得出的解真正满足原方程。(a)， (b)答案：(a).程序：A=7 2 1 -2;9

13、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.07910.5550 -0.0987小结：求矩阵A*B=C,已知A,C求B，则应用逆矩阵来求解即B=(A*(-1)*C6 假设有一组实测数据x0.10.20.30.40.50.60.70.8

14、0.91y2.32012.64702.97073.28853.60083.90904.21474.51914.82325.1275用最小二乘法拟合，求出相应的二次函数。答案：程序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*x2+3.2828*x+1.99677 考虑线性微分方程 (1)

15、试用Simulink搭建起系统的仿真模型，并绘制出仿真结果曲线。(2) 将给定的微分方程转换成状态方程，并建立S函数，再利用Simulink进行仿真。答案：(1)参数设置：Integrator1,2,3,4：Initial condition分别设置成0.2，0.5，0.5，1Gain0,1,2, 3：gain 分别设置成3，3，4，5Sine Wave： Frequency(rad/sec)设置成 4 ；phase(rad)设置成 pi/3Transfer Fcn ：Numerator coefficients设置成 1 ； Denominator coefficients 设置成1 5Tr

16、ansfer Fcn1：Numerator coefficients设置成1 ； Numerator coefficients设置成1 3 Add ： List of signs 设置成 +-仿真模型：（2）建立相应的S函数：functionsys,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

17、4,9 sys=; otherwise error(Unhandled flag=,num2str(flag);endfunctionsys,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)

18、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;由此建立仿真模型仿真曲线：小结：可以根据微分方程做出状态图并赋初值，得到系统的仿真模型以及仿真结果曲线；通过构建S函数，可以将状态方程封装在S函数，以简化运算。8 建立下图所示非线性系统的Simulink模型，并观察在单位阶跃信号输入下系统的输出曲线和误差曲线。答案：参数设置：Relay: Switch on point设置成0.5；S

19、witch on point设置成0.5；Output when on 设置成2.5；Output when off 设置成0Relay1: Switch on point 设置成-0.5；Switch on point设置成-0.5 ；Output when on 设置成2.5； Output when off 设置成0Transport Delay: Time delay设置成0.4Transfer Fcn: Numerator coefficients设置成30;Denominator coefficients设置成1 0.5 1Transfer Fcn1: Numerator coef

20、ficients设置成30;Denominator coefficients设置成1 6.5 1仿真模型：仿真曲线小结：1.用两个继电器并联得到具有死去特性的三位置继电特性 2.从图中可以看出输出比较稳定实验二 经典控制理论实验目的以MATLAB及Simulink为工具，对控制系统进行时域、频域及根轨迹分析。实验原理1、 时域分析法是根据系统的微分方程（或传递函数），利用拉普拉斯变换直接解出动态方程，并依据过程曲线及表达式分析系统的性能。时域响应指标如图所示。延迟时间td，指响应曲线第一次达到其终值一半所需要的时间。上升时间tr，指响应曲线从终值10%上升到终值90%所需要的时间；对于有振荡的

21、系统，也可定义为响应从零第一次上升到终值所需要的时间。上升时间是系统响应速度的一种度量。峰值时间tp，指响应超过终值达到第一个峰值所需要的时间。调节时间ts，指响应达到并保持在终值5%（或2%）内所需要的时间。超调量%，指响应的最大偏离量h(tp)与终值h()之差的百分比，即：稳态误差，描述系统稳态性能的一种性能指标。2、 频域分析法通常从频率特性出发对系统进行研究。在工程分析和设计中，通常把频率特性画成一些曲线，从频率特性曲线出发进行研究。这些曲线包括幅频特性和相频特性曲线，幅相频率特性曲线，对数频率特性曲线以及对数幅相曲线等，其中以幅相频率特性曲线，对数频率特性曲线应用最广。对于最小相位系

22、统，幅频特性和相频特性之间存在着唯一的对于关系，故根据对数幅频特性，可以唯一地确定相应的相频特性和传递函数。根据系统的开环频率特性去判断闭环系统的性能，并能较方便地分析系统参量对系统性能的影响，从而指出改善系统性能的途径。3、 根轨迹是求解闭环系统特征根的图解方法。由于控制系统的动态性能是由系统闭环零极点共同决定，控制系统的稳定性由闭环系统极点唯一确定，利用根轨迹确定闭环系统的零极点在s平面的位置，分析控制系统的动态性能。实验内容*1控制系统数学模型的转换自动控制系统计算机仿真教材第4章中的所有例题【例4-1】已知系统的传递函数如下，利用MATLAB建立其相应的传递函数系统模型答案程序：clc

23、clearnum=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函数将其转换为状态方程形式的模型答案：程序：clcclearnum=12 24 0 20;den=2 4 6 2 2;sys=tf(num,den)s=ss(sys) 结果：Transfer function: 12 s3 + 24 s2 + 20-2 s4 + 4 s3

24、 + 6 s2 + 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 d = u1 y1 0 Continuous-time model.小结：tf(num,den)函数用来求传递函数的状态模型【例4-3】某线性定常系统的状态空间表达式如下，求该系统的传递函数。答案：clcclearA=0 1 1;0 0 1;-10 -17 -8;B=0;0;1;C=5 6 1;D=0

25、SS=ss(A,B,C,D)G1=tf(SS)结果：Transfer function: s2 + 11 s + 5-s3 + 8 s2 + 27 s + 10小结：ss(A,B,C,D)函数是用来求状态方程的传递函数【例4-4】对于例4-3题中的线性定常系统，将其转换为zpk形式 答案：程序：clcclearA=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.4751p = -0.4199 -3.7901 + 3.0745i -3.7901 - 3.

26、0745ik = 1.0000 Zero/pole/gain: (s+10.52) (s+0.4751)-(s+0.4199) (s2 + 7.58s + 23.82)小结：zpk(z,p,k)是用来求状态方程的zpk形式的函数【例4-5】 某系统的传递函数如下，求它的部分分式形式。答案：程序：clcclearnum=20 10;den=1 15 74 120;R,P,H=residue(num,den)结果R = -55.0000 90.0000 -35.0000P = -6.0000 -5.0000 -4.0000H = 小结：R,P,H=residue(num,den)函数是用来求传递函

27、数的部分分式函数【例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.0

28、902 x3 -0.1164 0.1936 c = x1 x2 x3 y1 1 -1 0 y2 2 1 -1 d = u1 u2 y1 0 0 y2 0 0 Sampling time: 0.1Discrete-time model.小结：Gn=c2d(G,T)采用零阶保持器将系统模型离散化，求出采样周期为T离散化的系统方程的函数【例4-7】 已知系统如图所示，利用MATLAB求出系统的状态空间表达式。其中，sys1: sys2: 答案：程序：A1=-9 17;-1 3;B1=0 -1;-1 0;C1=-3 2;-13 18;D1=-1 0;-1 0;num=2;den=1 2;G1=ss(A

29、1,B1,C1,D1);G2=tf(num,den);sys=feedback(G1,G2,2,2,-1)结果a = x1 x2 x3 x1 -9 17 1 x2 -1 3 0 x3 -26 36 -2 b = u1 u2 x1 0 -1 x2 -1 0 x3 -2 0 c = x1 x2 x3 y1 -3 2 0 y2 -13 18 0 d = u1 u2 y1 -1 0 y2 -1 0 Continuous-time model.小结：sys=feedback(G1,G2,2,2,-1)是用来求反馈系统的函数【例4-8】 已知系统的状态空间表达式如下，求线性变换，将其变换成能控标准形。答案

30、程序：A=1 2 -1;0 2 1;1 -3 2;B=0;1;1;C=1 0 1;Qc=ctrb(A,B);A1=poly(A)if rank(Qc)=3 disp(The system is controllable ) Q=Qc*A1(3) A1(2) 1;A1(2) 1 0;1 0 0,P=inv(Q) Qb=P*A*Q,Bb=P*B,Cb=C*Qelse disp(The system is uncontrollable)end结果A1 = 1.0000 -5.0000 12.0000 -11.0000The system is controllable Q = 3.0000 1.0

31、000 0 2.0000 -2.0000 1.0000 7.0000 -6.0000 1.0000P = 0.2353 -0.0588 0.0588 0.2941 0.1765 -0.1765 0.1176 1.4706 -0.4706Qb = -0.0000 1.0000 0.0000 -0.0000 0 1.0000 11.0000 -12.0000 5.0000Bb = 0 0 1Cb = 10.0000 -5.0000 1.0000小结：Qc=ctrb(A,B)是用来求能控矩阵的秩，从如判断其能控性，进而求出其线性变换，并将其变换成能控标准形【例4-9】 已知系统传递函数如下，用不同状

32、态空间实现函数进行转换。答案：num=4 5 1;den=1 6 11 6;G=tf(num,den);A,B,C,D=tf2ss(num,den);G1=ss(A,B,C,D)G2=ss(G)G3=canon(G,modal)G4=canon(G,companion)结果：a = x1 x2 x3 x1 -6 -11 -6 x2 1 0 0 x3 0 1 0 b = u1 x1 1 x2 0 x3 0 c = x1 x2 x3 y1 4 5 1 d = u1 y1 0 Continuous-time model. a = x1 x2 x3 x1 -6 -2.75 -1.5 x2 4 0 0

33、x3 0 1 0 b = u1 x1 2 x2 0 x3 0 c = x1 x2 x3 y1 2 0.625 0.125 d = u1 y1 0 Continuous-time model. a = x1 x2 x3 x1 -3 0 0 x2 0 -2 0 x3 0 0 -1 b = u1 x1 -15.52 x2 -19.6 x3 5.745 c = x1 x2 x3 y1 -0.7086 0.3572 8.604e-016 d = u1 y1 0 Continuous-time model. a = x1 x2 x3 x1 0 0 -6 x2 1 0 -11 x3 0 1 -6 b = u

34、1 x1 1 x2 0 x3 0 c = x1 x2 x3 y1 4 -19 71 d = u1 y1 0 Continuous-time model. 小结：系统传递函数可以用不同状态空间实现函数进行转换2已知二阶系统闭环传递函数为：（1） 编写程序求解系统的阶跃响应；计算系统的闭环根、阻尼比、无阻尼振荡频率；（2）编程计算并在同一幅图上绘出和的阶跃响应曲线；（3）编程计算并在同一幅图上绘出和的阶跃响应曲线。（4）绘出以下系统的阶跃响应曲线，并与原系统响应曲线比较（将5条曲线画在同一幅图上） ；记录各系统的峰值和峰值时间，计算超调量。答案：（1）、（2）、（3）：程序：num1=10;den

35、1=1 2 10;num2=10;den2=1 2*sqrt(10) 10;num3=10;den3=1 4*sqrt(10) 10;num4=10;den4=4 4 10;num5=40;den5=1 4 40;G1=tf(num1,den1);G2=tf(num2,den2);G3=tf(num3,den3);G4=tf(num4,den4);G5=tf(num5,den5);wn,kai,p=damp(G1);wn=wn,ksai=kai,p=pstep(G1,-);hold ongtext(原系统)pausestep(G2,:);hold ongtext(ksai1=1)pausest

36、ep(G3,*);hold ongtext(ksai2=2)pausestep(G4,.);hold ongtext(wn1=wn/2)pausestep(G5,-.);gtext(wn2=2wn)结果：wn = 3.1623 3.1623ksai = 0.3162 0.3162p = -1.0000 - 3.0000i -1.0000 + 3.0000i小结：1）wn,kai,p=damp(G1)函数是用来求系统的闭环根、阻尼比、无阻尼振荡频率的2）阻尼比在（0,0.7）范围内，阻尼比越大，调节时间越短，上升时间越长 在大于0.7范围内，阻尼比越大，调节时间越长，上升时间越长3）无阻尼振荡频

37、率越大峰值时间越大（4）程序:num1=2 10;den1=1 2 10;num2=1 0.5 10;den2=1 2 10;num3=1 0.5 0;den3=1 2 10;num4=1 0;den4=1 2 10;num5=10;den5=1 2 10;G1=tf(num1,den1);G2=tf(num2,den2);G3=tf(num3,den3);G4=tf(num4,den4);G5=tf(num5,den5);step(G1,-);hold ongtext(G1)pausestep(G2,-);hold ongtext(G2)pausestep(G3,:);hold ongtext(G3)pausestep(G4,+);hold ongtext(G4)pausestep(G5,m-);gtext(原系统)结果：系统峰值峰值时间超调量G(s)1.35