《系统仿真技术课程实验报告(共44页).docx》由会员分享，可在线阅读，更多相关《系统仿真技术课程实验报告(共44页).docx（44页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、精选优质文档-倾情为你奉上 中 南 大 学课题名称： 系统仿真技术课程实验报告 学 院： 信息科学与工程学院 目录实验一 MATLAB中矩阵与多项式的基本运算一、实验任务1了解MATLAB命令窗口和程序文件的调用。2熟悉如下MATLAB的基本运算： 矩阵的产生、数据的输入、相关元素的显示； 矩阵的加法、乘法、左除、右除； 特殊矩阵：单位矩阵、“1”矩阵、“0”矩阵、对角阵、随机矩阵的产生和运算； 多项式的运算：多项式求根、多项式之间的乘除。二、基本命令训练1 eye(m)建立单位矩阵 eye(4)ans = 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 12 one(n)、one
2、s(m,n)建立全一矩阵one(n)运行未成功 ones(2)ans = 1 1 1 1 ones(2,4)ans = 1 1 1 1 1 1 1 13 zeros(m,n)全零矩阵 zeros(3)ans = 0 0 0 0 0 0 0 0 0 zeros(2,7)ans = 0 0 0 0 0 0 0 0 0 0 0 0 0 04 rand(m,n)服从均匀分布的随机矩阵 rand(2)ans = 0.8491 0.6787 0.9340 0.7577 rand(1,3)ans =0.7431 0.3922 0.6555PS：randn()服从正态分布的随机矩阵 randn(2)ans =
3、 -0.4326 0.1253 -1.6656 0.28775diag(v)对角矩阵 diag(0)ans = 0 diag(2,2)ans = 0 0 2 0 0 0 0 0 0 diag(4,-3)ans = 0 0 0 0 0 0 0 0 0 0 0 04 0 0 06.AB 、A/B、 inv(A)*B 、B*inv(A)A/B=A*inv(B),AB= inv(A)*B a=1 4;2 5a = 1 4 2 5 b=3 7;1 3b = 3 7 1 3 a/bans = -0.5000 2.5000 0.5000 0.5000 abans = -3.6667 -7.6667 1.66
4、67 3.6667 inv(a)*bans = -3.6667 -7.6667 1.6667 3.6667 b*inv(a)ans = -0.3333 1.66670.3333 0.33337.roots(p)求多项式p的解 a=1 2 1a = 1 2 1 roots(a)ans = -1 -18poly(a)求以a为根的多项式 a=1 0 1a = 1 0 1 poly(a)ans = 1 -2 1 09conv 、deconvconv(a,b)将多项式a,b相乘deconv（a,b）将多项式a,b相除 a=1 0 1a = 1 0 1 b=1 1b = 1 1 conv(a,b)ans
5、= 1 1 1 1 a=1 0 -1a = 1 0 -1 b=1 -1b = 1 -1 deconv(a,b)ans = 1 1 c=1 1 1c = 1 1 1 deconv(c,b)ans = 1 210A*B 与 A.*B的区别A*B按照矩阵乘法相乘A.*B对应位置相乘 a=1 2;3 4a = 1 2 3 4 b(1:2,1:2)=1b = 1 1 1 1 c=a*bc = 3 3 7 7 c1=a.*bc1 = 1 2 3 411who与whos的使用后不加变量就是查询所有变量who内存变量列表whos内存变量详细信息 whos a Name Size Bytes Class Att
6、ributes a 1x3 24 double 12disp、size(a)、length(a)的使用disp(a)显示变量asize(a)显示a的行，列数length(a)显示a的列数 A(1:2,1:3)=1A = 1 1 1 1 1 1 size(A)ans =2 3 length(A)ans = 3实验二 MATLAB绘图命令一、实验任务熟悉MATLAB基本绘图命令，掌握如下绘图方法：1坐标系的选择、图形的绘制；2图形注解（题目、标号、说明、分格线）的加入；3图形线型、符号、颜色的选取。二、基本命令训练1plot绘制平面曲线 2loglog绘制全对数坐标图形，X、Y轴都是取以10为底的
7、对数坐标 3semilogx绘制半对数坐标，其中X轴取以10为底的对数坐标，Y轴为线性坐标 4semilogy绘制半对数坐标，其中Y轴取以10为底的对数坐标，X轴为线性坐标 5polar绘制相角为theta，半径为radius的极坐标图形。 6title添加图片标示 7xlabel添加x坐标轴标示 8ylabel添加y坐标轴标示9text在坐标处书写注解 10grid绘制网格 11bar条形图 12stairs阶梯图13contour绘制等高线三、实验举例1t=0:pi/360:2*pi*22/3;x=93*cos(t)+36*cos(t*4.15);y=93*sin(t)+36*sin(t*
8、4.15);plot(y,x),grid;2 t=0:0.05:100;x=t;y=2*t;z=sin(2*t);plot3(x,y,z,b:)3 t=0:pi/20:2*pi;y=sin(x);stairs(x,y)4 th=pi/200:pi/200:2*pi;r=cos(2*th);polar(th,r),grid5 th=0:pi/10:2*pi;x=exp(j*th);plot(real(x),imag(x),r*);grid;实验三 MATLAB程序设计一、实验任务1熟悉MATLAB程序设计的方法和思路；2掌握循环、分支语句的编写，学会使用look for、help命令。二、程序举
9、例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,if = Columns 1 through 14 1 1 2 3 5 8 13 21 34 55 89 144 233 377 Columns 15 through 16 610 987i = 152 m=3;n=4;for i=1:m for j=1:n a(i,j)=1/(i+j-1); endendformat rataa = 1 1/2 1/3 1/4 1/2 1/3 1/4 1/5 1/3 1/4 1/5 1/6 3
10、m=3;n=4;for i=1:m for j=1:n a(i,j)=1/(i+j-1); endendaa = 1.0000 0.0000 0.3333 0.0000 0.0000 0.3333 0.0000 0.0000 0.3333 0.0000 0.0000 0.6667程序二结果为分数，程序三结果为小数4 x=input(请输入x的值:); if x=10 y=cos(x+1)+sqrt(x*x+1); else y=x*sqrt(x+sqrt(x); end y请输入x的值:10y = 10.8941请输入x的值:2y = 3.51475去掉多项式或数列开头的零项p=0 0 0 1
11、 3 0 2 0 0 9;for i=1:length(p),if p(1)=0,p=p(2:length(p); end;end;pp = 1 3 0 2 0 0 96. 建立MATLAB的函数文件，程序代码如下，以文件名ex2_4.m存盘function f=ffibno(n)%ffibno 计算斐波那契亚数列的函数文件%n可取任意自然数%程序如下f=1,1;i=1;while f(i)+f(i+1) lookfor ffibnoffibno - 计算斐波那契亚数列的函数文件ffibno - 计算斐波那契亚数列的函数文件 help ffibno ffibno 计算斐波那契亚数列的函数文件
12、n可取任意自然数 程序如下三、程序设计题 用一个MATLAB语言编写一个程序：输入一个自然数，判断它是否是素数，如果是，输出“It is one prime”，如果不是，输出“It is not one prime.”。要求通过调用子函数实现。最好能具有如下功能：设计较好的人机对话界面，程序中含有提示性的输入输出语句。能实现循环操作，由操作者输入相关命令来控制是否继续进行素数的判断。如果操作者希望停止这种判断，则可以退出程序。如果所输入的自然数是一个合数，除了给出其不是素数的结论外，还应给出至少一种其因数分解形式。例：输入 6， 因为6不是素数。则程序中除了有“It is not one pr
13、ime”的结论外，还应有：“6=2*3”的说明。a=input(请输入一个整数(endin0):);while(a=0)ifa=1disp(1isnotaprimeoracompositenumber.)elseifisprime(a)=1fprintf(%disaprime.n,a)elseifisprime(a)=0f=factor(a);m=length(f);fprintf(%disnotaprime,%d=%d,a,a,f(1)fori=2:mfprintf(*%d,f(i);endfprintf(n);enda=input(请输入一个整数(endin0):);end 实验四 MAT
14、LAB的符号计算与SIMULINK的使用一、实验任务1掌握MATLAB符号计算的特点和常用基本命令；2掌握SIMULINK的使用。二、程序举例1求矩阵对应的行列式和特征根a=sym(a11 a12;a21 a22);da=det(a)ea=eig(a)da = a11*a22-a12*a21 ea = 1/2*a11+1/2*a22+1/2*(a112-2*a11*a22+a222+4*a12*a21)(1/2) 1/2*a11+1/2*a22-1/2*(a112-2*a11*a22+a222+4*a12*a21)(1/2)2 求方程的解（包括精确解和一定精度的解）r1=solve(x2-x-
15、1)rv=vpa(r1)rv4=vpa(r1,4)rv20=vpa(r1,20)r1 = 1/2*5(1/2)+1/2 -1/2*5(1/2)+1/2rv = 1. -. rv4 = 1.618 -.6180 rv20 = 1. -.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, dB = 10 5 -8 11ans =
16、 a*d-b*cans = 1504 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*x2*y2-24*y4ans =(8*y2+7*x2)*(x2-3*y2)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.1687 0.4034
17、0.2592ans = 0.1687 0.4034 0.25926 对方程组求解a11*x1+a12*x2+a13*x3=b1a21*x1+a22*x2+a23*x3=b2a31*x1+a32*x2+a33*x3=b3syms a11 a12 a13 a21 a22 a23 a31 a32 a33 b1 b2 b3;A=a11,a12,a13;a21,a22,a23;a31,a32,a33;b=b1;b2;b3;X=linsolve(A,b) %调用linsolve函数求的解XX=Ab %用左除运算求解7syms a b t x y z;f=sqrt(1+exp(x);diff(f) %未指定
18、求导变量和阶数，按缺省规则处理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 =1/2/(1+exp(x)(1/2)*exp(x)ans = -2*sin(x)-x*cos(x)ans =-3*cos(x)+x*sin(x)ans =-b*cos(t)/a/sin(t)三、SIMULINK的使用Simulink部分：逐一熟悉各模块的使用，对下面的随动系统模型进行仿真，实验报告中包含：连好的系统模型及用Scop
19、e观测的结果G1(s)G2(s)R(s)C(s) 其中：R(s)为阶跃输入，C(s)为输出 plot(tout,simout.signals.values)实验五 MATLAB在控制系统分析中的应用一、实验任务1掌握MATLAB在控制系统时间响应分析中的应用；2掌握MATLAB在系统根轨迹分析中的应用； 3. 掌握MATLAB控制系统频率分析中的应用； 4. 掌握MATLAB在控制系统稳定性分析中的应用二、基本命令1. step求取连续系统的单位阶跃响应 2. impulse求取连续系统的单位脉冲响应 3. initial求取连续系统的零输入响应 4. lsim求取连续系统的任意输入响应 5.
20、 rlocfind求根轨迹某处的增益值6. bode对数频率特性图也称伯德图 7. margin求取幅值稳定裕度、相角稳定裕度、幅值穿越频率、相角穿越频率及系统稳定性等信息 8. nyquist极坐标频率特性图也称奈奎斯特图 9. Nichols对数幅相特性也称尼克尔斯图 10. cloop绘制单位反馈时闭环系统频率特性曲线三、程序举例1求下面系统的单位阶跃响应%程序如下.num=4 ; den=1 , 1 , 4 ;step(num , den)y , x , t=step(num , den) ;tp=spline(y , t , max(y) %计算峰值时间max(y) %计算峰值tp
21、= 1.5708ans =1.44192. 求如下系统的单位阶跃响应 %程序如下a=0,1;-6,-5;b=0;1;c=1,0;d=0;y,x=step(a,b,c,d);plot(y)3. 求下面系统的单位脉冲响应：%程序如下 num=4 ; den=1 , 1 ,4 ;impulse(num,den)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(
22、1 , 2 , 2) ; 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)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.
23、 已知开环系统传递函数，绘制系统的根轨迹，并分析其稳定性%程序如下：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,de
24、n);num,den=cloop(num1,den1,-1);impulse(num,den)title(impulse response(k=56)Select a point in the graphics windowselected_point = 0.6066 + 4.1988ik = 125.2930p = -7.1346 0.5713 + 4.2182i 0.5713 - 4.2182i -2.0080 8. 作如下系统的bode图%程序如下：n=1 , 1 ; d=1 , 4 , 11 , 7 ; bode(n , d)9. 系统传函如下求有理传函的频率响应，然后在同一张图上绘
25、出以四阶伯德近似表示的系统频率响应 %程序如下：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)
26、;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.654111. 二阶系统为：令wn=1,分别作出=2 , 1 , 0.707 , 0.5时的nyquist曲线。 %程序如下： n=1 ; d1=1 , 4 , 1 ; d2=1 , 2
27、 , 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)13. 分别由w的自动变量和人工变量作下列系统的nyquist曲线： %程序如下：n=1 ; d=1 , 1 ,0 ;nyquist(n ,d) ; %自动变量n=1 ; d=1 , 1 ,0
28、; 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,den);p figure(1)nyquist(num,den)figure(2)num2,den2=cloop(num,den);impulse(num2,den2);15已知系统为：作该系统的nichols曲
29、线。 %程序如下：. 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时闭环系统的稳定性。%程序如下： d
30、1=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.2393 + 0.8579i -0.2393 - 0.8579ians = -3.3089 0.1545 + 1.7316i 0.1545 - 1.7316i18. 系统的状态方程为： 试确定系统的稳定性。 %程序如下： a=-4,-3,0 ; 1,0,0 ; 0,1,0 ; b=1;0;0 ; c=0,1,2 ; d=0 ;
31、eig(a) %求特征根rank(ctrb(a,b)ans = 0 -1 -3ans = 3实验六 连续系统数字仿真的基本算法一、实验任务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)+
32、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% 利用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)+
33、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*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)程序运行结果如下： 时间 解析解 欧拉法 RK2法 RK4法 Columns 1 th
34、rough 4 0 1.0000 1.0000 1.0000 0.0000 1.9923 1.0000 1.6667 0.0000 1.9874 1.3333 1.6793 0.0000 1.9133 1.1068 1.7113 0.0000 1.9710 1.0625 1.1403 1.0000 1.8877 1.8238 1.5690 Column 5 1.0000 1.5307 1.2607 1.0262 1.7527 1.11932. 考虑如下二阶系统： 在上的数字仿真解（已知：，），并将不同步长下的仿真结果与解析解进行精度比较。 说明：已知该微分方程的解析解分别为： 采用RK4法进行计算，选择状态变量： 则有如下状态空间模型及初值条件 采用RK4法进行计算。程序如下：clearh=input(请输入步长h=); % 输入步长M=round(10/h); % 置总计算步数t(1)=0; % 置自变量初值y_0(1)=100; y_05(1)=100;