系统仿真技术实验报告

1、系统仿真技术实验报告学 院：信息科学与工程学院专 业：电气工程及其自动化班 级：XXXXX 姓 名：XXXXX 学 号：XXXXXXXXXXX 目 录一、 实验一（MATLAB中矩阵与多项式的基本运算）.41.1 实验题目.41.2 实验扩展.131.3 实验心得.20二、 实验二（MATLAB绘图命令）.212.1 实验题目.212.2 实验扩展.262.3 实验心得.29三、 实验三（MATLAB程序设计）.303.1 实验题目.303.2 实验扩展.373.3 实验心得.38四、 实验四（MATLAB的符号计算与SIMULINK的使用）.394.1 实验题目.394.2 实验扩展.484

2、.3 实验心得.49五、实验五（MATLAB在控制系统分析中的应用）.505.1 实验题目.505.2 实验心得.71六、 实验六（连续系统数字仿真的基本算法）.726.1实验题目.72七、总结.78 实验一：MATLAB中矩阵与多项式的基本运算1.1 实验题目1.1.1 实验任务： 1了解MATLAB命令窗口和程序文件的调用。2熟悉如下MATLAB的基本运算： 矩阵的产生、数据的输入、相关元素的显示； 矩阵的加法、乘法、左除、右除； 特殊矩阵：单位矩阵、“1”矩阵、“0”矩阵、对角阵、随机矩阵的产生和运算； 多项式的运算：多项式求根、多项式之间的乘除。1.1.2 基本命令训练：1eye(m)

3、2one(n)、ones(m,n) 3zeros(m,n)实验结果：4rand(m,n) 5diag(v)分析：rand(m,n):产生m*n的随机阵diag（v）：以v的元素作为对角阵的元素产生对角阵实验结果：6AB 、A/B、 inv(A)*B 、B*inv(A)分析：X=AB是AX=B的解，且AB=A(-1)BX=A/B是XB=A的解，且A/B=AB(-1)Inv(A)*B:A的逆矩阵与左乘BB*inv(A)：B左乘A的逆矩阵实验结果：7roots(p)8Poly分析：roots(p):多项式p的根poly(C):用C的元素为根构造多项式实验结果：9conv 、deconv分析：conv

4、:进行多项式的乘法deconv:进行多项式的除法实验结果：10A*B 与 A.*B的区别分析：A*B：进行矩阵运算A.*B：进行数组运算实验结果：11who与whos的使用分析：who:列出当前工作空间中所有变量名字whos:列出当前工作空间中所有变量，以及它们的名字、尺寸、所占字节数、属性等信息。实验结果：12disp、size(a)、length(a)的使用分析：disp:屏幕输出函数size(a):求矩阵的大小(获取数组的行数和列数)length(a):数组长度(即行数或列数中较大者)实验结果：12 实验扩展1.2.1 向量及其运算1）、向量产生实验结果：2）、向量基本运算分析：dot:

5、维数相同的两个向量的点积，结果是一个标量cross:两个向量的叉积就是一个过两个相交向量的交点且垂直于两个向量所在平面的向量，结果是一个矢量实验结果：1.2.2 矩阵及其运算1）、矩阵产生通过冒号可以截取指定矩阵中的部分2）、矩阵基本运算det(A):返回矩阵A的行列式的值eig(A):返回矩阵A的特征值inv(A):求A的逆矩阵fliplr(A):矩阵A关于中轴的反转矩阵(左右反转)flipud(A):上下反转rot90(A):矩阵A逆时针旋转90度diag(A):抽取主对角线元素构成向量rem(A,x):A矩阵除以模数x后的余数1.2.3 数组及其运算实验结果：1.2.4 多项式及其运算分

6、析：polyval(P,X):计算多项式P在X点的值（若X是向量或矩阵，则采用数组运算（点运算）polyvalm(P,X):以方阵X为自变量，计算多项式的值，采用矩阵运算1.3、 实验心得MATLAB中数值运算的部分较为繁多，但和之前学过的C语言等也有相同的部分。但MATLAB的语言更为精简，比较方便。MATLAB是一款工程性质十分强的软件，只有大量的应用实践，才能够对其中的函数运用自如，才能理解其中的内涵。实验二 MATLAB绘图命令2.1 实验题目21.1实验任务熟悉MATLAB基本绘图命令，掌握如下绘图方法：1坐标系的选择、图形的绘制；2图形注解（题目、标号、说明、分格线）的加入；3图形

7、线型、符号、颜色的选取。2.1.2 基本命令训练1plot 2loglog 3semilogx 4semilogy5polar 6title 7xlabel 8ylabel9text 10grid 11bar 12stairs13contour实验结果：2.2 实验扩展 2.2.1 图形的修饰 1. 曲线颜色: b(blue),c(cyan),g(green),k(black),m(magenta),r(red),w(white),y(yellow)2.曲线线型：-（实线，默认）， ：（点线）， -.(点画线), -(虚线) 3.数据点形：.(实点), +(十字形), 。(圆圈) , *(星号

8、), x(叉号), s(正方形), d(棱形), h(六角形), p(五角形), (右三角), (左三角) 实验举例：2.2.2 三维图形的绘制仿真举例：t=0:pi/30:10*pi; %设定时间范围plot3(cos(t),sin(t),t,-b,Linewidth,4)%绘制三维图形，并作修饰grid on %加网格axis square %坐标为方形figure(2) %新建图形窗口2stem3(cos(t),sin(t),t,-.k) %绘制三维火柴杆图实验结果：2.3 实验心得用户应用MATLAB，只需指定绘图方式，并提供充足的绘图数据，可以得到所需的图形，也可对图形加以修饰。二维

9、图形的绘制基本函数为plot。较为复杂的三维图形分为三维曲线，三维网线图和三维曲面图。对符号函数或是无法求出函数参数之间显示关系的隐函数，可以使用符号函数绘制图形的函数。实验三 MATLAB程序设计3.1 实验题目3.1.1实验任务1熟悉MATLAB程序设计的方法和思路；2掌握循环、分支语句的编写，学会使用look for、help命令。31.2 程序训练1计算11000之内的斐波那契亚数列 f=1,1;i=1;while f(i)+f(i+1)1000 f(i+2)=f(i)+f(i+1); i=i+1;endf,i实验结果：2 m=3;n=4;for i=1:m for j=1:n a(i

10、,j)=1/(i+j-1); endendformat rata实验结果：3 m=3;n=4;for i=1:m for j=1:n a(i,j)=1/(i+j-1);endenda实验结果：程序2和程序3的区别：format rat 小数分数表示。4 x=input(请输入x的值:); if x=10 y=cos(x+1)+sqrt(x*x+1); else y=x*sqrt(x+sqrt(x); end y实验结果：5去掉多项式或数列开头的零项p=0 0 0 1 3 0 2 0 0 9;for i=1:length(p),if p(1)=0,p=p(2:length(p); end;end

11、;p实验结果：6 建立MATLAB的函数文件，程序代码如下，以文件名ex2_4.m存盘function f=ffibno(n)%ffibno 计算斐波那契亚数列的函数文件%n可取任意自然数%程序如下f=1,1;i=1;while f(i)+f(i+1)n f(i+2)=f(i)+f(i+1); i=i+1;end输入完毕后在MATLAB的命令窗口输入ex2_4(200)，得到运行结果。在MATLAB的命令窗口输入lookfor ffibno，得到结果：ex2_4.m: %ffibno 计算斐波那契亚数列的函数文件在MATLAB的命令窗口输入help ex2_4，得到结果： ffibno 计算斐

12、波那契亚数列的函数文件 n可取任意自然数 程序如下实验结果：3.1.3 程序设计用一个MATLAB语言编写一个程序：输入一个自然数，判断它是否是素数，如果是，输出“It is one prime”，如果不是，输出“It is not one prime.”。要求通过调用子函数实现。最好能具有如下功能：设计较好的人机对话界面，程序中含有提示性的输入输出语句。能实现循环操作，由操作者输入相关命令来控制是否继续进行素数的判断。如果操作者希望停止这种判断，则可以退出程序。如果所输入的自然数是一个合数，除了给出其不是素数的结论外，还应给出至少一种其因数分解形式。例：输入 6， 因为6不是素数。则程序中除

13、了有“It is not one prime”的结论外，还应有：“6=2*3”的说明。实验结果：3.2 实验扩展3.2.1 选择结构举例：g=input(please input grade:,s); %由用户交互输入成绩等级switch(g) %判断用户输入caseA,adisp(85100);caseB,bdisp(7084);caseC,cdisp(6069);caseD,ddisp(=period./2);%看rem函数功能lsim(num,den,u,t);实验结果： 7. 已知开环系统传递函数，绘制系统的根轨迹，并分析其稳定性 %程序如下： num=1 2;den1=1 4 3;d

14、en=conv(den1,den1);figure(1)rlocus(num,den)k,p= rlocfind(num,den) figure(2)k=55;num1=k*1 2;den=1 4 3;den1=conv(den,den);num,den=cloop(num1,den1,-1);impulse(num,den)title(impulse response (k=55) )figure(3)k=56;num1=k*1 2;den=1 4 3;den1=conv(den,den);num,den=cloop(num1,den1,-1);impulse(num,den)title(i

15、mpulse response(k=56)实验结果： 8. 作如下系统的bode图 %程序如下： n=1 , 1 ; d=1 , 4 , 11 , 7 ; bode(n , d)实验结果：9. 系统传函如下 求有理传函的频率响应，然后在同一张图上绘出以四阶伯德近似表示的系统频率响应 %程序如下： num=1;den=conv(1 2,conv(1 2,1 2); w=logspace(-1,2); t=0.5;m1,p1=bode(num,den,2);p1=p1-t*w*180/pi;n2,d2=pade(t,4);numt=conv(n2,num);dent=(conv(den,d2);m

16、2,p2=bode(numt,dent,w);subplot(2,1,1);semilogx(w,20*log10(m1),w,20*log10(m2),g-);grid on ; title(bode plot);xlabel(frequency);ylabel(gain);subplot(2,1,2);semilogx(w,p1,w,p2,g-);grid on;xlabel(frequency);ylabel(phase);实验结果：10. 已知系统模型为 求它的幅值裕度和相角裕度 %程序如下： n=3.5; d=1 2 3 2; Gm,Pm,Wcg,Wcp=margin(n,d)实验结

17、果：11. 二阶系统为：令wn=1,分别作出=2 , 1 , 0.707 , 0.5时的nyquist曲线。 %程序如下：n=1 ; d1=1 , 4 , 1 ; d2=1 , 2 , 1 ; d3=1 , 1.414 , 1; d4=1,1,1; nyquist(n,d1) ; hold on nyquist(n,d2) ; nyquist(n,d3) ; nyquist(n,d4) ;实验结果：12. 已知系统的开环传递函数为 绘制系统的Nyqusit图，并讨论系统的稳定性 %程序如下： G=tf(1000,conv(1,3,2,1,5);nyquist(G);axis(square)实验

18、结果：13. 分别由w的自动变量和人工变量作下列系统的nyquist曲线： %程序如下： n=1 ; d=1 , 1 ,0 ;nyquist(n ,d) ; %自动变量n=1 ; d=1 , 1 ,0 ; w=0.5 : 0.1 : 3 ;nyquist(n , d , w) ; %人工变量实验结果：14. 一多环系统，其结构图如下，使用Nyquist频率曲线判断系统的稳定性。 %程序如下： k1=16.7/0.0125;z1=0;p1=-1.25 -4 -16;num1,den1=zp2tf(z1,p1,k1);num,den=cloop(num1,den1);z,p,k=tf2zp(num

19、,den);p figure(1)nyquist(num,den)figure(2)num2,den2=cloop(num,den);impulse(num2,den2);实验结果：15. 已知系统为：作该系统的nichols曲线。 %程序如下： n=1 ; d=1 , 1 , 0 ; ngrid(new) ; nichols(n , d) ;实验结果：16. 已知系统的开环传递函数为：当k=2时，分别作nichols曲线和波特图。 %程序如下：num=1;den=conv(conv(1 0,1 1),0.5 1);subplot(1,2,1);nichols(num,den);grid; %

20、 nichols曲线subplot(1,2,2);g=tf(num,den);bode(feedback(g,1,-1);grid; %波特图实验结果：17. 系统的开环传递函数为： 分别确定k=2和k=10时闭环系统的稳定性。%程序如下：d1=1 , 3 , 2 , 0 ; n1=2 ;nc1 , dc1=cloop(n1 , d1 ,-1) ;roots(dc1)d2=d1 ; n2=10 ;nc2 , dc2=cloop(n2 , d2,-1) ; roots(dc2)实验结果：18. 系统的状态方程为： 试确定系统的稳定性。 %程序如下： a=-4,-3,0 ; 1,0,0 ; 0,1

21、,0 ; b=1;0;0 ; c=0,1,2 ; d=0 ;eig(a) %求特征根rank(ctrb(a,b)实验结果：5.2 实验心得通过本次试验，熟悉了解了MATLAB中SIMULINK的应用，特别是在自动控制系统建模中，它非常的简便易操作，相比输入语句的建模来说，它形象而快速，不仅便于理解和分析，而且大大提高了仿真速度。同时，也巩固了MATLAB的画图功能，同样为之后对控制系统的仿真打下基础，对于相关传递函数图形的绘制有很大的帮助。实验六 连续系统数字仿真的基本算法6.1 实验题目6.1.1 实验任务1理解欧拉法和龙格-库塔法的基本思想；2理解数值积分算法的计算精度、速度、稳定性与步长

22、的关系。6.1.2 程序练习1. 取h=0.2，试分别用欧拉法、RK2法和RK4法求解微分方程的数值解，并比较计算精度。 注：解析解： %程序如下cleart(1)=0; % 置自变量初值y(1)=1; y_euler(1)=1; y_rk2(1)=1; y_rk4(1)=1; % 置解析解和数值解的初值h=0.2; % 步长% 求解析解for k=1:5 t(k+1)=t(k)+h; y(k+1)=sqrt(1+2*t(k+1);end% 利用欧拉法求解for k=1:5 y_euler(k+1)=y_euler(k)+h*(y_euler(k)-2*t(k)/y_euler(k);end%

23、 利用RK2法求解for k=1:5 k1=y_rk2(k)-2*t(k)/y_rk2(k); k2=(y_rk2(k)+h*k1)-2*(t(k)+h)/(y_rk2(k)+h*k1); y_rk2(k+1)=y_rk2(k)+h*(k1+k2)/2;end% 利用RK4法求解for k=1:5 k1=y_rk4(k)-2*t(k)/y_rk4(k); k2=(y_rk4(k)+h*k1/2)-2*(t(k)+h/2)/(y_rk4(k)+h*k1/2); k3=(y_rk4(k)+h*k2/2)-2*(t(k)+h/2)/(y_rk4(k)+h*k2/2); k4=(y_rk4(k)+h*

24、k3)-2*(t(k)+h)/(y_rk4(k)+h*k3); y_rk4(k+1)=y_rk4(k)+h*(k1+2*k2+2*k3+k4)/6;end% 输出结果disp( 时间 解析解 欧拉法 RK2法 RK4法)yt=t, y, y_euler, y_rk2, y_rk4;disp(yt)实验结果：2. 考虑如下二阶系统： 在上的数字仿真解（已知：，），并将不同步长下的仿真结果与解析解进行精度比较。 说明：已知该微分方程的解析解分别为： 采用RK4法进行计算，选择状态变量： 则有如下状态空间模型及初值条件 采用RK4法进行计算。程序如下：clearh=input(请输入步长h=); %

25、 输入步长M=round(10/h); % 置总计算步数t(1)=0; % 置自变量初值y_0(1)=100; y_05(1)=100; % 置解析解的初始值(y_0和y_05分别对应于为R=0和R=0.5)x1(1)=100; x2(1)=0; % 置状态向量初值y_rk4_0(1)=x1(1); y_rk4_05(1)=x1(1); % 置数值解的初值 % 求解析解for k=1:M t(k+1)=t(k)+h; y_0(k+1)=100*cos(t(k+1); y_05(k+1)=100*exp(-t(k+1)/2).*cos(sqrt(3)/2*t(k+1)+100*sqrt(3)/3

26、*exp(-t(k+1)/2).*sin(sqrt(3)/2*t(k+1);end% 利用RK4法求解% R=0for k=1:M k11=x2(k); k12=-x1(k); k21=x2(k)+h*k12/2; k22=-(x1(k)+h*k11/2); k31=x2(k)+h*k22/2; k32=-(x1(k)+h*k21/2); k41=x2(k)+h*k32; k42=-(x1(k)+h*k31); x1(k+1)=x1(k)+h*(k11+2*k21+2*k31+k41)/6; x2(k+1)=x2(k)+h*(k12+2*k22+2*k32+k42)/6; y_rk4_0(k+

27、1)=x1(k+1);end% R=0.5for k=1:M k11=x2(k); k12=-x1(k)-x2(k); k21=x2(k)+h*k12/2; k22=-(x1(k)+h*k11/2)-(x2(k)+h*k12/2); k31=x2(k)+h*k22/2; k32=-(x1(k)+h*k21/2)-(x2(k)+h*k22/2); k41=x2(k)+h*k32; k42=-(x1(k)+h*k31)-(x2(k)+h*k32); x1(k+1)=x1(k)+h*(k11+2*k21+2*k31+k41)/6; x2(k+1)=x2(k)+h*(k12+2*k22+2*k32+k42)/6; y_rk4_05(k+1)=x1(k+1);end% 求出误差最