中南大学自动化胡杨系统仿真实验报告完整版.doc

中南大学系统仿真实验报告 指导老师： 实验者： 学号： 专业班级： 完成时间： 50实验一 MATLAB中矩阵与多项式的基本运算基本命令训练：1 eye(m) 取n=3，程序如下： eye(3)ans = 1 0 0 0 1 0 0 0 1结论：eye(n)用于产生nn维的单位矩阵，在这里n取3，故产生33维单位矩阵。2 one(n)、ones(m,n) 对ones(n) 取n=5，对ones(m,n)取m=3，n=5，程序如下： ones(5)ans = 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ones(3,5)ans = 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1结论：ones(n)用于产生nn维的全1矩阵，在这里n取5，故产生5行5列全1矩阵。ones(m,n)用于产生mn维的全1矩阵，在这里m取3，n取5，故产生3行5列的全1矩阵。3 zeros(m,n)取m=3，n=2，程序如下： zeros(3,2)ans = 0 0 0 0 0 0结论：zeros(m,n)用于产生mn维全0矩阵，在这里m取3，n取2，故产生3行2列全0矩阵。4rand(m,n)取m=3，n=4，程序如下： rand(3,4)ans = 0.9501 0.4860 0.4565 0.4447 0.2311 0.8913 0.0185 0.61540.6068 0.7621 0.8214 0.7919结论：rand(m,n)用于产生mn维平均分布的随机矩阵，在这里m取3，n取4，故产生了3行4列的随机矩阵5diag(v)先创建33的魔方矩阵v，在进行diag(v)运算，程序如下： v=magic(3)diag(v)v = 8 1 6 3 5 7 4 9 2ans = 8 5 2结论：diag(v)用于得到矩阵v的对角元素6AB 、A/B、 inv(A)*B 、B*inv(A)先创建A、B两个矩阵，在进行运算，程序如下： A=1,2;3,4; B=5,6;7,8; a=ABb=A/Bc=inv(A)*Bd=B*inv(A)a = -3 -4 4 5b = 3.0000 -2.0000 2.0000 -1.0000c = -3.0000 -4.0000 4.0000 5.0000d = -1.0000 2.0000 -2.0000 3.0000结论： / 表示矩阵右除， 表示矩阵左除，inv(A)表示求A的逆矩阵，由实验结果可知，矩阵左除与右除结果不一样，矩阵左乘与右乘结果也不一样，AB是求AX=B的解，A/B是求XB=A的解。所以编程求解的时候要注意区分他们的区别。 7、roots(p) syms x; a=3*x.3+2*x+5; p=3,0,2,5 roots(p)p = 3 0 2 5ans = 0.5000 + 1.1902i 0.5000 - 1.1902i -1.0000 结论：roots(p)函数用于求多项式的根，以向量形式输入多项式的系数，对应降幂排列，然后调用函数，即可求得对应多项式的根。 8、 poly A=1,2;3,4; poly(A)ans =1.0000 -5.0000 -2.0000 结论：poly(A)用于求矩阵A的特征多项式的系数 9conv 、deconv A=1,2; B=3,4; a=conv(A,B)b=deconv(A,B)a = 3 10 8b = 0.3333结论：使用conv函数对多项式进行乘法运算，其使用格式为c=conv(a,b),其中a和b为两个多项式的系数向量，c为相乘所生成的多项式的系数向量。使用deconv(a,b)完成除法运算。10 A*B 与 A.*B的区别 A=1,2; B=5,6; a=A*B A=1,2; B=5,6; b=A.*B a = 17 b = 5 12 结论：A.*B称为“点乘”、“位乘“，即为两个行列数相同的矩阵，对应位置一一相乘，得到的结果依位置对应到结果矩阵中，而A*B为矩阵乘法，要求前者A的列数与后者B行数对应。11who与whos的使用 A=1,2;3,4; who whos Your variables are: A Name Size Bytes Class Attributes A 2x2 32 double 结论：who给出变量的名称清单；而whos给出所有变量的详细信息。12 disp、size(a)、length(a)的使用 a=helloworld; disp(a) a=1,2,3,4; B=size(a) C=length(a) helloworld B = 1 4 C = 4 结论：disp函数的作用是直接将内容输出在Matlab命令窗口中，size(a)表示矩阵每个维度的长度，length(a)表示矩阵a的最大的长度。实验二 MATLAB绘图命令基本命令训练1plot 2loglog 3semilogx 4semilogy5polar 6title 7xlabel 8ylabel9text 10grid 11bar 12stairs13contour1t=0:pi/360:2*pi*22/3;x=93*cos(t)+36*cos(t*4.15);y=93*sin(t)+36*sin(t*4.15);plot(y,x),grid; 实验结果为：t=0:pi/360:2*pi*22/3;x=93*cos(t)+36*cos(t*4.15);y=93*sin(t)+36*sin(t*4.15);plot(y,x) 实验结果为： 实验结论：plot()用于绘制二维曲线，grid用于切换有无网格的状态。2 t=0:0.05:100;x=t;y=2*t;z=sin(2*t);plot3(x,y,z,b:) 实验结果为： 实验结论：plot3(x,y,z)用于绘制三维曲线，b 表示设置曲线的颜色为蓝色，：表示曲线线型为点线，格式为plot3(函数参数，函数参数，曲线参数设置)3 t=0:pi/20:2*pi;y=sin(x);stairs(x,y) 实验结果为： 实验结论：stairs(x,y)表示绘制出的二维曲线为阶梯图。4 th=pi/200:pi/200:2*pi;r=cos(2*th);polar(th,r),grid 实验结果为： 实验结论：polar()用于绘制二维曲线的极坐标图。5 th=0:pi/10:2*pi;x=exp(j*th);plot(real(x),imag(x),r*);grid;实验结果为： 实验结论：r表示设置曲线颜色为红色，*表示曲线的数据点形为星号。6、 x=0:1000; y=0:1000; loglog(x,y);title(Loglog);grid on;实验结果为： 实验结论：loglog（）用于绘制横纵轴均为对数刻度的图形，title（）用于为图形添加标题，本例为添加Loglog作为标题。7、 x=0:1000; y=0:1000; semilogx(x,y);title(Loglog);grid on;实验结果为：将semilogx换成semilogy程序如下： x=0:1000; y=0:1000; semilogy(x,y);title(Loglog);grid on;实验结果为： 实验结论：semilogx（）用于绘制半对数图，其中x轴坐标为对数，若换成semilogy则表示y轴坐标为对数。8、 x=0:1000; y=0:1000; plot(x,y); x=0:1000; y=0:1000; plot(x,y);grid on;xlabel(fontsize20itxrm);ylabel(fontsize20y);text(500,500,中点)实验结果为：实验结论：xlabel和ylabel分别表示给x轴和y轴添加标注，text(x,y,string)用于给图形坐标（x,y）处书写注释，本程序给x轴和y轴分别标注x,y,，在（500,500）坐标处注释“中点”。9、 t=0:pi/100:2*pi; alpha=3; y=sin(alpha*t); bar(t,y);grid on;实验结果为：实验结论：bar(x,y)用于绘制二维条形图。10、 x=-8:0.5:8; y=-8:0.5:8; xx,yy=meshgrid(x,y); c=sqrt(xx.2+yy.2)+eps; z=sin(c)./c; contour(xx,yy,z)实验结果为： 实验结论：contour(x,y,z)用于绘制等高线。补充实验：多窗口绘制图形subplot() subplot(2,2,1);t=0:pi/200:2*pi;y=sin(t);plot(t,y);subplot(2,2,2);t=0:pi/200:2*pi;y=cos(t);plot(t,y);subplot(2,2,4);t=0:pi/200:2*pi;y=t;plot(t,y);实验结果为：实验结论：本实验测试subplot()函数，由实验结果可知，subplot（）函数中某一个未编写并不会影响整个函数的运行，只是未编写的那个部分不显示，其他的照常显示，比如编写了subplot(2,2,1),subplot(2,2,2),subplot(2,2,4)，但是未编写subplot(2,2,3),那么结果只显示subplot(2,2,1),subplot(2,2,2),subplot(2,2,4)中的结果，而且顺序按原位置，而subplot(2,2,3)的不会显示。实验三 MATLAB程序设计 1计算11000之内的斐波那契亚数列 f=1,1; i=1; while f(i)+f(i+1) f,i f = Columns 1 through 10 1 1 2 3 5 8 13 21 34 55 Columns 11 through 16 89 144 233 377 610 987 i = 15 2. m=3; n=4; for i=1:m for j=1:n a(i,j)=1/(i+j-1); end end format rat a a = 1 1/2 1/3 1/4 1/2 1/3 1/4 1/5 1/3 1/4 1/5 1/6 3. m=3; n=4; for i=1:m for j=1:n a(i,j)=1/(i+j-1); end end a a = 1 0.5 0.33333 0.25 0.5 0.33333 0.25 0.2 0.33333 0.25 0.2 0.16667 实验2用了format rat,结果为分数表示，实验3没有则用小数表示。4、 x=input(请输入x的值:); if x15 y=x*sqrt(x+sqrt(x); else y=x; end y 请输入x的值:10 y = 10.054 x=input(请输入x的值:); if x15 y=x*sqrt(x+sqrt(x); else y=x; end y 请输入x的值:11 y = 115、去掉多项式或数列开头的零项 . 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; p p = 1 3 0 2 0 0 96、建立MATLAB的函数文件，程序代码如下，以文件名ex2_4.m存盘点击File-New-Function建立文件，文件名为ex2_4.m,结果如下：在MATLAB的命令窗口输入ex2_4(200)，得到运行结果： ex2_4(200)ans = 1 1 2 3 5 8 13 21 34 55 89 144在MATLAB的命令窗口输入lookfor ffibno，得到结果： lookfor ffibnoex2_4 - ffibno 计算斐波那契亚数列的函数文件在MATLAB的命令窗口输入help ex2_4，得到结果： help ex2_4 ffibno 计算斐波那契亚数列的函数文件 n可取任意自然数 程序如下3、 程序设计题 function sushu n=input(请输入一个数n:); if(n=1) fprintf(1既不是素数也不是合数n); disp(是否继续？);disp(1.是； 2.否 );b=input( );if b=1 sushu;else disp(谢谢使用！); break;endend if(n=2) fprintf(2是素数n); disp(是否继续？);disp(1.是； 2.否 );b=input( );if b=1 sushu;else disp(谢谢使用！); break;endendif(n=3) fprintf(3是素数n); disp(是否继续？);disp(1.是； 2.否 );b=input( );if b=1 sushu;else disp(谢谢使用！); break;endendif(n3) for i=2:(n-1) if mod(n,i)=0 a=n/i; t=0; fprintf(%d=%d*%dn,n,a,i); else t=1; end end if (t=1) clear t; fprintf(%d不是素数n,n); else fprintf(%d是素数n,n); clear t; enddisp(是否继续？);disp(1.是； 2.否 );b=input( );if b=1 sushu;else disp(谢谢使用！);endend运行结果为： sushu请输入一个数n:11既不是素数也不是合数是否继续？1.是； 2.否 1请输入一个数n:22是素数是否继续？1.是； 2.否 1请输入一个数n:3838=19*238=2*1938不是素数是否继续？1.是； 2.否 1请输入一个数n:2323不是素数是否继续？1.是； 2.否 2谢谢使用！实验四 MATLAB的符号计算与SIMULINK的使用程序举例1 求矩阵对应的行列式和特征根 a=sym(1 0 0 1); da=det(a) ea=eig(a) da = 1 ea = 1实验结论：det()函数用于计算矩阵对于的行列式的值，eig()函数用于计算矩阵的特征值和特征向量2. 求方程的解（包括精确解和一定精度的解） r1=solve(x2-x-1) rv=vpa(r1) rv4=vpa(r1,4) rv20=vpa(r1,20) r1 = 1/2 - 5(1/2)/2 5(1/2)/2 + 1/2 rv = -0.61803398874989484820458683436564 1.6180339887498948482045868343656 rv4 = -0.618 1.618 rv20 = -0.6180339887498948482 1.6180339887498948482实验结论：vpa（s,n）称为变精度算法函数，表示将s表示为n位有效数的符号对象3、 a=sym(a);b=sym(b);c=sym(c);d=sym(d); %定义4个符号变量w=10;x=5;y=-8;z=11; %定义4个数值变量A=a,b;c,d %建立符号矩阵AB=w,x;y,z %建立数值矩阵Bdet(A) %计算符号矩阵A的行列式det(B) %计算数值矩阵B的行列式 A = a, b c, d B = 10 5 -8 11 ans = a*d - b*c ans = 150 4、 syms x y; s=(-7*x2-8*y2)*(-x2+3*y2); expand(s) %对s展开 collect(s,x) %对s按变量x合并同类项(无同类项) factor(ans) % 对ans分解因式 ans = 7*x4 - 13*x2*y2 - 24*y4 ans = 7*x4 + (-13*y2)*x2 - 24*y4 ans = (7*x2 + 8*y2)*(x2 - 3*y2)实验结论：expand函数用于多项式的展开运算，collect函数用于符号表达式的展开运算和合并同类项，factor用于对函数进行因式分解。5、对方程 AX=b求解 A=34,8,4;3,34,3;3,6,8;b=4;6;2;X=linsolve(A,b) %调用linsolve函数求解Ab %用另一种方法求解X = 0.067482 0.16137 0.10367ans = 0.067482 0.16137 0.10367实验结论：运用linsove(A,b)与Ab结果一样，都用于对AX=b进行求解 6、对方程组求解 a11*x1+a12*x2+a13*x3=b1 a21*x1+a22*x2+a23*x3=b2 a31*x1+a32*x2+a33*x3=b3 A=a11,a12,a13;a21,a22,a23;a31,a32,a33; b=b1;b2;b3; XX=Ab XX = (a12*a23*b3-a12*b2*a33+a13*a32*b2-a13*a22*b3+b1*a22*a33-b1*a32*a23)/(a11*a22*a33-a11*a32*a23-a21*a12*a33+a32*a21*a13-a22*a31*a13+a31*a12*a23) -(a11*a23*b3-a11*b2*a33-a21*a13*b3-a23*a31*b1+b2*a31*a13+a21*b1*a33)/(a11*a22*a33-a11*a32*a23-a21*a12*a33+a32*a21*a13-a22*a31*a13+a31*a12*a23) (a32*a21*b1-a11*a32*b2+a11*a22*b3-a22*a31*b1-a21*a12*b3+a31*a12*b2)/(a11*a22*a33-a11*a32*a23-a21*a12*a33+a32*a21*a13-a22*a31*a13+a31*a12*a23) 7、 syms a b t x y z; f=sqrt(1+exp(x);diff(f) %未指定求导变量和阶数，按缺省规则处理f=x*cos(x);diff(f,x,2) %求f对x的二阶导数diff(f,x,3) %求f对x的三阶导数f1=a*cos(t);f2=b*sin(t);diff(f2)/diff(f1) %按参数方程求导公式求y对x的导数 ans = exp(x)/(2*(exp(x) + 1)(1/2) ans = - 2*sin(x) - x*cos(x) ans = x*sin(x) - 3*cos(x) ans = -(b*cos(t)/(a*sin(t)3、 SIMULINK的使用选择合适子模块构造控制系统如下：选择Simulation菜单下的start命令运行仿真，在示波器（Scope）中观察结果：R(s)为阶跃输入，C(s)为输出 实验结果为：实验五 MATLAB在控制系统分析中的应用基本命令 1. step 2. impulse 3. initial 4. lsim 5. rlocfind 6. bode 7. margin 8. nyquist 9. Nichols 10. cloop1、求下面系统的单位阶跃响应 num=4 ; den=1 , 1 , 4 ;step(num , den)y , x , t=step(num , den) ;tp=spline(y , t , max(y) max(y)tp = 1.6062ans =1.4441实验结论：step(num,den)用于绘制系统阶跃响应曲线，spline(y,t,max(y)函数是由y,t的值计算max(y)对应的函数值t，在本例中即峰值时间，而max(y)用于计算峰值。2、求如下系统的单位阶跃响应 a=0,1;-6,-5;b=0;1;c=1,0;d=0;y,x=step(a,b,c,d);plot(y)程序结果为：实验结论：step(a,b,c,d)用于绘制系统阶跃响应曲线3、求下面系统的单位脉冲响应%程序如下： num=4 ; den=1 , 1 ,4 ;impulse(num,den)实验结果：实验结论：impulse(num,den)函数用于求取系统单位脉冲响应，其用法基本同step函数。4、已知二阶系统的状态方程为：求系统的零输入响应和脉冲响应。 %程序如下： a=0 , 1 ; -10 , -2 ; b=0 ; 1 ; c=1 , 0 ; d=0 ; x0=1 ,0; subplot(1 , 2 , 1) ; initial(a , b , c ,d,x0) subplot(1 , 2 , 2) ; impulse(a , b , c , d)实验结果为：实验结论：initial(a , b , c ,d,x0)用于求解零输入响应系统，x0为初始条件，impulse(a , b , c , d)用于求单位脉冲响应。5、系统传递函数为：输入正弦信号时，观察输出信号的相位差。 %程序如下： num=1 ; den=1 , 1 ; t=0 : 0.01 : 10 ; u=sin(2*t) ; hold on plot(t , u , r) lsim(num , den , u , t)程序结果为：实验结论：lsim(num , den , u , t)用于求系统对任意输入u的响应。6、有一二阶系统，求出周期为4秒的方波的输出响应%程序如下：num=2 5 1;den=1 2 3;t=(0:.1:10);period=4;u=(rem(t,period)=period./2);%看rem函数功能lsim(num,den,u,t); 实验结果为：7. 已知开环系统传递函数，绘制系统的根轨迹，并分析其稳定性 %程序如下： num=1 2;den1=1 4 3;den=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(impulse response(k=56)实验结果：Select a point in the graphics windowselected_point = -3.3886 + 2.3602ik = 23.5611p = -5.0868 -0.4355 + 2.2831i -0.4355 - 2.2831i -2.0423 实验结论：den=conv(A,B)用于多项式A,B以系数行向量表示，进行相乘，rlocus(num,den)用于绘制指定系统的根轨迹。分析得系统是不稳定的。函数rlocfind用于计算给定一组根的根轨迹增益8、作如下系统的bode图 %程序如下： n=1 , 1 ; d=1 , 4 , 11 , 7 ; bode(n , d) 实验结果为：实验结论：bode(n,d)用于绘制Bode图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);m2,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)实验结果为：Gm = 1.1433Pm = 7.1688Wcg = 1.7323Wcp = 1.6541 实验结论：Gm,Pm,Wcg,Wcp=margin(n,d)给出系统相对稳定参数，返回参数分别为幅值裕度、相角裕度、幅值穿越频率、相角穿越频率。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) ; 实验结果为： 实验结论：nyquist(sys,w)函数用于绘制系统Nyquist图，由用户指定选取频率范围。12、已知系统的开环传递函数为 绘制系统的Nyqusit图，并讨论系统的稳定性 %程序如下： G=tf(1000,conv(1,3,2,1,5);nyquist(G);axis(square)实验结果为：Warning: This plot type does not support this option for the AXIS command.InD:MATLAB6p5p1toolboxcontrolctrlguisctrluisaxesgroupaddbypass.m at line 114 In D:MATLAB6p5p1toolboxmatlabgraph2dprivatemwbypass.p at line 24 In D:MATLAB6p5p1toolboxmatlabgraph2daxis.m at line 7613. 分别由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) ; %人工变量实验结果为：实验结论：nyquist(n , d , w)用于求指定范围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,den);p figure(1) nyquist(num,den)figure(2)num2,den2=cloop(num,den);impulse(num2,den2); 实验结果为：p = -10.597 + 36.215i -10.597 - 36.215i -0.056187 实验结论：num1,den1=zp2tf(z1,p1,k1)用于将系统函数的零极点转化为系统函数一般形式的系数，cloop()函数用于将系统通过正负反馈连接成闭环系统，分析得系统稳定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; % 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)实验结果为：ans = -2.5214 -0.23931 + 0.85787i -0.23931 - 0.85787ians = -3.3089 0.