系统仿真技术课程实验报告

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 on

2、e(n)、ones(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

3、=0.7431 0.3922 0.6555PS：randn()服从正态分布的随机矩阵>> randn(2)ans = -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 =

4、 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.6667 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 -

5、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 = 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 =

6、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 Attributes a 1x3 24 double 12di

7、sp、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为底的对数坐标

8、 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*4.15

9、);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掌握循环、分支语句的编写，学会使用l

10、ook for、help命令。二、程序举例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/

11、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); 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的值:

12、2y = 3.51475去掉多项式或数列开头的零项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;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)<n f(i+2)=f(i)+f(i+1); i=i+1;end>> lookfor ffibnoffib

13、no - 计算斐波那契亚数列的函数文件ffibno - 计算斐波那契亚数列的函数文件>> help ffibno ffibno 计算斐波那契亚数列的函数文件 n可取任意自然数 程序如下三、程序设计题 用一个MATLAB语言编写一个程序：输入一个自然数，判断它是否是素数，如果是，输出“It is one prime”，如果不是，输出“It is not one prime.”。要求通过调用子函数实现。最好能具有如下功能：设计较好的人机对话界面，程序中含有提示性的输入输出语句。能实现循环操作，由操作者输入相关命令来控制是否继续进行素数的判断。如果操作者希望停止这种判断，则可以退出程序。

14、如果所输入的自然数是一个合数，除了给出其不是素数的结论外，还应给出至少一种其因数分解形式。例：输入 6， 因为6不是素数。则程序中除了有“It is not one prime”的结论外，还应有：“6=2*3”的说明。>>a=input('请输入一个整数(end in 0):'); while(a=0)     if a=1           disp('1&

15、#160;is not a prime or a composite number.')    elseif isprime(a)=1         fprintf('%d is a prime.n',a)     elseif isprime(a)=0 

16、0;     f=factor(a);m=length(f);          fprintf('%d is not a prime,%d=%d',a,a,f(1)              for i=2:m  

17、0;          fprintf('*%d',f(i);   end         fprintf('n');      end       a=input('请输入一个整数(end in&#

18、160;0):'); end 实验四 MATLAB的符号计算与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

19、)(1/2)2 求方程的解（包括精确解和一定精度的解）r1=solve('x2-x-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 %建立符号

20、矩阵AB=w,x;y,z %建立数值矩阵Bdet(A) %计算符号矩阵A的行列式det(B) %计算数值矩阵B的行列式A = a, b c, dB = 10 5 -8 11ans = 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求解

21、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 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

22、=linsolve(A,b) %调用linsolve函数求的解XX=Ab %用左除运算求解7syms 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 =1/2/(1+exp(x)(1/2)*exp(x)ans = -2*sin(x)-x*cos(x)ans =-3*cos(x)+x*sin(x)ans

23、 =-b*cos(t)/a/sin(t)三、SIMULINK的使用Simulink部分：逐一熟悉各模块的使用，对下面的随动系统模型进行仿真，实验报告中包含：连好的系统模型及用Scope观测的结果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在控制系统稳定性分析中的应用二、基本命令

24、1. step求取连续系统的单位阶跃响应 2. impulse求取连续系统的单位脉冲响应 3. initial求取连续系统的零输入响应 4. lsim求取连续系统的任意输入响应 5. rlocfind求根轨迹某处的增益值6. bode对数频率特性图也称伯德图 7. margin求取幅值稳定裕度、相角稳定裕度、幅值穿越频率、相角穿越频率及系统稳定性等信息 8. nyquist极坐标频率特性图也称奈奎斯特图 9. Nichols对数幅相特性也称尼克尔斯图 10. cloop绘制单位反馈时闭环系统频率特性曲线三、程序举例1求下面系统的单位阶跃响应%程序如下.num=4 ; den=1 , 1 , 4

25、 ;step(num , den)y , x , t=step(num , den) ;tp=spline(y , t , max(y) %计算峰值时间max(y) %计算峰值tp = 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

26、; 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)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秒的方波的输出响应 %程序

27、如下：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,de

28、n1,-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 = 0.6066 + 4.1988ik = 125.2930p = -7.1346 0.571

29、3 + 4.2182i 0.5713 - 4.2182i -2.0080 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

30、,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. 已知系

31、统模型为求它的幅值裕度和相角裕度 %程序如下： 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 , 1 ; d3=1 , 1.414 , 1; d4=1,1,1; nyquist(n,d1) ; hold on nyquist(n,d2) ; nyquist(n,d3) ; nyquist(

32、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 ; w=0.5 : 0.1 : 3 ;nyquist(n , d , w) ; %人工变量14. 一多环系统，其结构图如下，使用Nyquist频率曲线判断系统的稳定性。 %程序如

33、下：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曲线。 %程序如下：. n=1 ; d=1 , 1 , 0 ; ngrid(new) ; nichols(n , d) ;16. 已知系统的开环传递函数为：当k=2时，分别作nic

34、hols曲线和波特图。 %程序如下： 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=

35、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 ;eig(a) %求特征根rank(ctrb(a,b)ans = 0 -1 -3ans = 3实验六 连续系统数字仿真的基本算法一、实验任务1理解欧拉法和龙格-库塔法的基本思想；2

36、理解数值积分算法的计算精度、速度、稳定性与步长的关系；二、程序举例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

37、(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)+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/

38、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 through 4 0 1.0000 1.0000 1.0000 0.0000 1.9923 1.0000 1.6

39、667 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

40、('请输入步长h='); % 输入步长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*s