系统仿真实验报告.docx
《系统仿真实验报告.docx》由会员分享，可在线阅读，更多相关《系统仿真实验报告.docx（33页珍藏版）》请在冰豆网上搜索。
系统仿真实验报告
（此文档为word格式，下载后您可任意编辑修改！
）
系统仿真实验报告
实验一MATLAB中矩阵与多项式的基本运算1
实验二MATLAB绘图命令6
实验三MATLAB程序设计8
实验四MATLAB的符号计算与SIMULINK的使用12
实验五MATLAB在控制系统分析中的应用15
实验六连续系统数字仿真的基本算法26
实验一MATLAB中矩阵与多项式的基本运算
一、实验任务
1．了解MATLAB命令窗口和程序文件的调用。
2．熟悉如下MATLAB的基本运算：
①矩阵的产生、数据的输入、相关元素的显示；
②矩阵的加法、乘法、左除、右除；
③特殊矩阵：
单位矩阵、“1”矩阵、“0”矩阵、对角阵、随机矩阵的产生和运算；
④多项式的运算：
多项式求根、多项式之间的乘除。
二、基本命令训练
1．eye（m）
>>eye（4）
ans=
1000
0100
0010
0001
2．one（n）、ones（m,n）
>>ones（3,4）
ans=
1111
1111
1111
3．zeros（m,n）
>>zeros（2,3）
ans=000
000
4．rand（m,n）
>>rand（3,4）
ans=
0.95010.48600.45650.4447
0.23110.89130.01850.6154
0.60680.76210.82140.7919
5．diag（v）
>>v=[111];
X=diag（v）
X=
100
010
001
>>X=diag（v,1）
X=
0100
0010
0001
0000
>>X=diag（v,2）
X=
00100
00010
00001
00000
00000
6．A\B、A/B、inv（A）*B、B*inv（A）
Matlab提供了两种除法运算：
左除（\）和右除（/）。
一般情况下，x=a\b是方程a*x=b的解，而x=b/a是方程x*a=b的解。
a=[1 2 3;4 2 6;7 4 9];b=[4;1;2];x=a\b
x=
-1.5000
2.0000
0.5000
>>a=[1127;385;436];b=flipud（a）
b=
436
385
1127
>>a/b
ans=
001
010
100
如果a为非奇异矩阵，则a\b和b/a可通过a的逆矩阵与b阵得到：
a\b=inv（a）*b
b/a=b*inv（a）
7．roots（p）
求解多项式x^3-6x^2-72x-27=0的根
>>p=[1-6-72-27];roots（p）
ans=
12.1229
-5.7345
-0.3884
8．Poly
>>r=[123];poly（r）
ans=
1-611-6
创建以之为根的方程为：
x^3-6*x^2+11*x-6=0
>>A=[123;456;780];
P=poly（A）
P=
1.0000-6.0000-72.0000-27.0000
创建特征方程为λ^3-6*λ^2-72*λ-27=0
9．conv、deconv
>>a=[1234];b=[3-5823-7];conv（a,b）
ans=
31736438771-28
>>deconv（a,b）
ans=
0
>>deconv（b,a）
ans=
3-11
>>a=[10-1];b=[1-1];deconv（a,b）
ans=
11
>>conv（a,b）
ans=
1-1-11
10．A*B与A.*B的区别
>>a=[12;34];b=[45;67];c=a*b
c=
1619
3643
>>a.*b
ans=
410
1828
11．who与whos的使用
>>who
Yourvariablesare:
APaansbcr
>>whos
NameSizeBytesClass
A3x372doublearray
P1x432doublearray
a2x232doublearray
ans2x232doublearray
b2x232doublearray
c2x232doublearray
r1x324doublearray
Grandtotalis32elementsusing256bytes
12．disp、size（a）、length（a）的使用
>>d=[12;34];disp（d）
12
34
>>disp（'d'）
d
>>size（d）
ans=
22
>>length（d）
ans=
2
实验二MATLAB绘图命令
一、实验任务
熟悉MATLAB基本绘图命令，掌握如下绘图方法：
1．坐标系的选择、图形的绘制；
2．图形注解（题目、标号、说明、分格线）的加入；
3．图形线型、符号、颜色的选取。
二、基本命令训练
1．plot2．loglog3．semilogx4．semilogy
5．polar6．title7．xlabel8．ylabel
9．text10．grid11．bar12．stairs
13．contour
三、实验举例
1．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）,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）,grid
5．th=[0:
pi/10:
2*pi];
x=exp（j*th）;
plot（real（x）,imag（x）,'r*'）;
grid;
实验三MATLAB程序设计
一、实验任务
1．熟悉MATLAB程序设计的方法和思路；
2．掌握循环、分支语句的编写，学会使用lookfor、help命令。
二、程序举例
1．计算1~1000之内的斐波那契亚数列
f=[1,1];
i=1;
whilef（i）+f（i+1）<1000
f（i+2）=f（i）+f（i+1）;
i=i+1;
end
f,i
f=
Columns1through14
1123581321345589144233377
Columns15through16
610987
i=15
2．m=3;
n=4;
fori=1:
m
forj=1:
n
a（i,j）=1/（i+j-1）;
end
end
formatrat
a
a=
11/21/31/4
1/21/31/41/5
1/31/41/51/6
3．m=3;
n=4;
fori=1:
m
forj=1:
n
a（i,j）=1/（i+j-1）;
end
end
a
a=
11/21/31/4
1/21/31/41/5
1/31/41/51/6
4．x=input（'请输入x的值:
'）;
ifx==10
y=cos（x+1）+sqrt（x*x+1）;
else
y=x*sqrt（x+sqrt（x））;
end
y
x=input（'请输入x的值:
'）
请输入x的值:
10
x=10
>>ifx==10
y=cos（x+1）+sqrt（x*x+1）;
else
y=x*sqrt（x+sqrt（x））;
end
y
y=10.0543
>>x=input（'请输入x的值:
'）
请输入x的值:
8
x=8
>>ifx==10
y=cos（x+1）+sqrt（x*x+1）;
else
y=x*sqrt（x+sqrt（x））;
end
y
y=26.3253
5．去掉多项式或数列开头的零项
p=[0001302009];
fori=1:
length（p）,ifp
（1）==0,p=p（2:
length（p））;
end;
end;
p
p=
1302009
6．建立MATLAB的函数文件，程序代码如下，以文件名ex2_4.m存盘
functionf=ffibno（n）
%ffibno计算斐波那契亚数列的函数文件
%n可取任意自然数
%程序如下
f=[1,1];
i=1;
whilef（i）+f（i+1）f（i+2）=f（i）+f（i+1）;
i=i+1;
end
>>exe2_4（200）
ans=
1123581321345589144
>>lookforffibno
exe2_4.m:
%ffibno计算斐波那契亚数列的函数文件
>>helpexe2_4
ffibno计算斐波那契亚数列的函数文件
n可取任意自然数
三、程序设计题
用一个MATLAB语言编写一个程序：
输入一个自然数，判断它是否是素数，如果是，输出“Itisoneprime”，如果不是，输出“Itisnotoneprime.”。
要求通过调用子函数实现。
最好能具有如下功能：
①设计较好的人机对话界面，程序中含有提示性的输入输出语句。
②能实现循环操作，由操作者输入相关命令来控制是否继续进行素数的判断。
如果操作者希望停止这种判断，则可以退出程序。
③如果所输入的自然数是一个合数，除了给出其不是素数的结论外，还应给出至少一种其因数分解形式。
例：
输入6，因为6不是素数。
则程序中除了有“Itisnotoneprime”的结论外，还应有：
“6=2*3”的说明。
>> close all;
c=input（'是否进行素数运算 1为是 0为否:
'）
while c==1
a=input（'请输入一个自然数:
'）;
if a==1
disp（'1既不是质数也不是合数'）;
else
if factor（a）==a
disp（'It IS one prime'）;
else
disp（'It is not one prime'）;
b=factor（a）;
fprintf（'%3d=',a）
for j=1:
（length（b）-1）
fprintf（'%3d *',b（j））
end
fprintf（'%3d \n *',b（length（b）））
end
end
c=input（'是否进行素数运算 1为是 0为否:
'）
end
是否进行素数运算 1为是 0为否:
1
c =
1
请输入一个自然数:
1
1既不是质数也不是合数
是否进行素数运算 1为是 0为否:
1
c =
1
请输入一个自然数:
6
It is not one prime
6= 2 * 3
实验四MATLAB的符号计算与SIMULINK的使用
一、实验任务
1．掌握MATLAB符号计算的特点和常用基本命令；
2．掌握SIMULINK的使用。
二、程序举例
1．求矩阵对应的行列式和特征根
a=sym（'[a11a12;a21a22]'）;
da=det（a）
ea=eig（a）
da=
a11*a22-a12*a21
ea=
1/2*a11+1/2*a22+1/2*（a11^2-2*a11*a22+a22^2+4*a12*a21）^（1/2）
1/2*a11+1/2*a22-1/2*（a11^2-2*a11*a22+a22^2+4*a12*a21）^（1/2）
2．求方程的解（包括精确解和一定精度的解）
r1=solve（'x^2-x-1'）
rv=vpa（r1）
rv4=vpa（r1,4）
rv20=vpa（r1,20）
r1=
1/2*5^（1/2）+1/2
1/2-1/2*5^（1/2）
rv=
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]%建立符号矩阵A
B=[w,x;y,z]%建立数值矩阵B
det（A）%计算符号矩阵A的行列式
det（B）%计算数值矩阵B的行列式
A=
[a,b]
[c,d]
B=
105
-811
ans=
a*d-b*c
ans=
150
4．symsxy;
s=（-7*x^2-8*y^2）*（-x^2+3*y^2）;
expand（s）%对s展开
collect（s,x）%对s按变量x合并同类项（无同类项）
factor（ans）%对ans分解因式
ans=
7*x^4-13*x^2*y^2-24*y^4
ans=
7*x^4-13*x^2*y^2-24*y^4
ans=
（8*y^2+7*x^2）*（x^2-3*y^2）
5．对方程AX=b求解
A=[34,8,4;3,34,3;3,6,8];
b=[4;6;2];
X=linsolve（A,b）%调用linsolve函数求解
A\b%用另一种方法求解
X=
0.0675
0.1614
0.1037
ans=
0.0675
0.1614
0.1037
6．对方程组求解
a11*x1+a12*x2+a13*x3=b1
a21*x1+a22*x2+a23*x3=b2
a31*x1+a32*x2+a33*x3=b3
symsa11a12a13a21a22a23a31a32a33b1b2b3;
A=[a11,a12,a13;a21,a22,a23;a31,a32,a33];
b=[b1;b2;b3];
X=linsolve（A,b）%调用linsolve函数求的解
XX=A\b%用左除运算求解
XX=（a12*a23*b3-a12*b2*a33+a13*a32*b2-a13*b3*a22+b1*a33*a22-b1*a32*a23）/（a11*a33*a22-a33*a12*a21-a31*a13*a22+a31*a12*a23+a32*a21*a13-a11*a32*a23）
-（a11*a23*b3-a11*b2*a33-a21*a13*b3+b2*a31*a13-a23*a31*b1+a21*b1*a33）/（a11*a33*a22-a33*a12*a21-a31*a13*a22+a31*a12*a23+a32*a21*a13-a11*a32*a23）
（-a31*b1*a22-a11*a32*b2+a11*b3*a22-b3*a12*a21+a31*a12*b2+a32*a21*b1）/（a11*a33*a22-a33*a12*a21-a31*a13*a22+a31*a12*a23+a32*a21*a13-a11*a32*a23）
7．symsabtxyz;
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=
-b*cos（t）/a/sin（t）
三、SIMULINK的使用
1．在命令窗口中输入：
simulink（回车）得到如下simulink模块：
2．双击打开各模块，选择合适子模块构造控制系统，例如：
3．双击各子模块可修改其参数，选择Simulation菜单下的start命令运行仿真，在示波器（Scope）中观察结果。
实验五MATLAB在控制系统分析中的应用
一、实验任务
1．掌握MATLAB在控制系统时间响应分析中的应用；
2．掌握MATLAB在系统根轨迹分析中的应用；
3.掌握MATLAB控制系统频率分析中的应用；
4.掌握MATLAB在控制系统稳定性分析中的应用
二、基本命令
1.step2.impulse3.initial4.lsim5.rlocfind
6.bode7.margin8.nyquist9.Nichols10.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=
1.6062
ans=
1.4441
2.求如下系统的单位阶跃响应
%程序如下：
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（1,2,2）;impulse（a,b,c,d）
5：
系统传递函数为：
输入正弦信号时，观察输出信号的相位差。
%程序如下：
num=[1];den=[1,1];
t=0:
0.01:
10;
u=sin（2*t）;
holdon
plot（t,u,‘r:
’）
lsim（num,den,u,t）
6.有一二阶系统，求出周期为4秒的方波的输出响应
%程序如下：
num=[251];
den=[123];
t=（0:
.1:
10）;
period=4;
u=（rem（t,period）>=period./2）;%看rem函数功能
u=double（u）
lsim（num,den,u,t）
7.已知开环系统传递函数，绘制系统的根轨迹，并分析其稳定性
%程序如下：
num=[12];
den1=[143];
den=conv（den1,den1）;
figure
（1）
rlocus（num,den）
[k,p]=rlocfind（num,den）
figure
（2）
k=55;
num1=k*[12];
den=[143];
den1=conv（den,den）;
[num,den]=cloop（num1,den1,-1）;
impulse（num,den）
title（'impulseresponse（k=55）'）
figure（3）
k=56;
num1=k*[12];
den=[143];
den1=conv（den,den）;
[num,den]=cloop（num1,den1,-1）;
impulse（num,den）
title（'impulseresponse（k=56）'）
Selectapointinthegraphicswindow
selected_point=
-1.5687+0.6708i
k=
2.4232
p=
-3.7489
-2.3284
-0.9613+0.8137i
-0.9613-0.8137i
8.作如下系统的bode图
%程序如下：
n=[1,1];d=[1,4,11,7];
bode（n,d）
9.系统传函如下
求有理传函的频率响应，然后在同一张图上绘出以四阶伯德近似表示的系统频率响应
%程序如下：
num=[1];den=conv（[12],conv（[12],[12]））;
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--'）;
gridon;title（'bodeplot'）;xlabel（'frequency'）;ylabel（'gain'）;
subplot（2,1,2）;semilogx（w,p1,w,p2,'g--'）;gridon;
xlabel（'frequency'）;ylabel（'phase'）;
10.已知系统模型为
求它的幅值裕度和相角裕度
%程序如下：
n=[3.5];d=[1232];
[Gm,Pm,Wcg,Wcp]=margin（n,d）
Gm=
1.1433
Pm=
7.1688
Wcg=
1.7323
Wcp
升级会员 