系统仿真实验报告.docx
《系统仿真实验报告.docx》由会员分享，可在线阅读，更多相关《系统仿真实验报告.docx（44页珍藏版）》请在冰豆网上搜索。
系统仿真实验报告
中南大学系
统仿真技术实验报告
作者：
黄伟
学号：
0909122327
班级：
自动化1205
指导老师：
胡杨
时间：
2014-6-5
一、目的要求
1．了解MATLAB命令窗口和程序文件的调用。
2．熟悉如下MATLAB的基本运算：
①矩阵的产生、数据的输入、相关元素的显示；
②矩阵的加法、乘法、左除、右除；
③特殊矩阵：
单位矩阵、“1”矩阵、“0”矩阵、对角阵、随机矩阵的产生和运算；
④多项式的运算：
多项式求根、多项式之间的乘除。
3熟悉MATLAB基本绘图命令，掌握如下绘图方法：
①坐标系的选择、图形的绘制；
②图形注解（题目、标号、说明、分格线）的加入；
③图形线型、符号、颜色的选取。
4．熟悉MATLAB程序设计的方法和思路；
5．掌握循环、分支语句的编写，学会使用lookfor、help命令。
6．参照已知程序，改动程序中的参数和输入量，验证输出结果。
7．使用lookfor、help命令，验证输出结果
8．掌握MATLAB符号计算的特点和常用基本命令；
9．掌握SIMULINK的使用。
10．掌握MATLAB在控制系统时间响应分析中的应用；
11．掌握MATLAB在系统根轨迹分析中的应用；
12.掌握MATLAB控制系统频率分析中的应用；
13.掌握MATLAB在控制系统稳定性分析中的应用
14．理解欧拉法和龙格-库塔法的基本思想；
15．理解数值积分算法的计算精度、速度、稳定性与步长的关系；
二、仪器设备
计算机1台
Matlab软件1个
三、内容步骤及数据处理
I.MATLAB中矩阵与多项式的基本运算
1．eye（m）
eye（3）
ans=
100
010
001
2．one（n）、ones（m,n）
one
（1）
ans=
1
ones（2,3）
ans=
111
111
3．zeros（m,n）
zeros（3,4）
ans=
0000
0000
0000
4．rand（m,n）
5．diag（v）
p=rand（4,4）;%3,4题%随即数列
>>disp（p）
0.64770.74470.36850.9294
0.45090.18900.62560.7757
0.54700.68680.78020.4868
0.29630.18350.08110.4359
>>diag（p）%提取对角线数据
ans=
0.6477
0.1890
0.7802
0.43596．A\B、A/B、inv（A）*B、B*inv（A）
A=rand（4,4）;
>>B=rand（4,4）;
>>A\B
ans=
-0.17642.15690.04780.2260
0.09261.04880.65841.2567
1.4410-2.3537-0.1777-1.2059
0.1421-1.88930.42390.0172
A/B
ans=
-0.44100.39650.9073-0.5299
-1.56360.94080.54770.9849
3.24340.5827-1.3008-1.3753
1.78290.5311-0.3391-0.8652
>>inv（A）*B
ans=
-0.17642.15690.04780.2260
0.09261.04880.65841.2567
1.4410-2.3537-0.1777-1.2059
0.1421-1.88930.42390.0172
>>B*inv（A）
ans=
-1.2220-0.2555-1.84303.3873
0.45350.76600.8263-0.7193
-0.8123-0.4002-2.58684.1538
-1.92150.1006-2.27683.7547
7．roots（p）
P=[121];%求多项式根
>>roots（P）
ans=
-1
-1
8．Poly
P=[1234;5678;9012;1111];%求特征多项式
>>Poly（P）
ans=
1.0000-9.0000-30.000030.0000-0.00009．conv、deconv
10．A*B与A.*B的区别
11．who与whos的使用
P=[1234;5678;9012;1111];
A=[1111;0000;1234;1222];
>>A*P%向量
ans=
1691215
0000
42182430
31162126
>>A.*P%对应项相乘
ans=
1234
0000
9038
1222
Who%显示参数
Yourvariablesare:
ABPansnp
Whos%显示参数及属性
NameSizeBytesClassAttributes
A4x4128double
B4x4128double
P4x4128double
ans4x4128double
n1x18double
p3x372double
12．disp、size（a）、length（a）的使用
P=[1234;5678];
>>disp（P）
1234
5678
>>size（P）
ans=
24
>>length（P）
ans=
4
II.MATLAB绘图命令
1．plot
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,’r’）%红色
2．loglog
t=[0123456789101112];
x=10.^t;
y=100.^t;1
>>loglog（x,y,’g’）%绿色
3．semilogx
4．Semilogy
x=0:
.1:
10;semilogx（10.^x，x）
x=0:
.1:
10;semilogy（x,10.^x）
5．polar
th=[pi/200:
pi/200:
2*pi]';
r=cos（2*th）;
polar（th,r）,grid
6．title
th=[pi/200:
pi/200:
2*pi]';
r=cos（2*th）;
polar（th,r）;grid;title（'极坐标'）
7．xlabel
x=0:
0.1:
10;semilogx（10.^x,x）;xlabel（'x轴对数变换'）
8．Ylabel
x=0:
.1:
10;semilogy（x,10.^x）;ylabel（'y轴对数变换'）
9．text
10．grid
x=0:
.1:
10;semilogy（x,10.^x）;text（2,100,'（2,100）'）;grid%9,10题
11．bar
x=[456787654];
bar（x）%绘制直方图
12．Stairs
x=0:
0.1:
10;%阶梯
y=x;
stairs（x,y）
13．contour
x=1:
1:
10;%等高线
y=1:
1:
10;
[xx,yy]=meshgrid（y,x）;
z=rand（10,10）;
contour（xx,yy,z,15）
III.MATLAB程序设计
1．计算1~1000之内的斐波那契亚数列
f=[1,1];
i=1;
whilef（i）+f（i+1）<10000
f（i+2）=f（i）+f（i+1）;
i=i+1;
end
>>f,i
f=
Columns1through13
1123581321345589144233
Columns14through20
3776109871597258441816765
i=
19
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=
1.00000.50000.33330.2500
0.50000.33330.25000.2000
0.33330.25000.20000.1667
请比较程序2与程序3的区别
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的值：
'）;
if（x==0）
disp（'分母为零不合法，请重新输入！
'）;
else
y=（3^x）/x+cos（x）;
end
y
请输入x的值：
5
y=
48.8837
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
输入完毕后在MATLAB的命令窗口0输入ex2_4（200），得到运行结果。
ans=
1123581321345589144
IV.MATLAB的符号计算与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*5^（1/2）+1/2
rv=
1.6180339887498948482045868343656
-.61803398874989484820458683436560
rv4=
1.618
-.6180
rv20=
1.6180339887498948482
-.61803398874989484820
3.
>>a=sym（'a'）;b=sym（'b'）;c=sym（'c'）;d=sym（'d'）;
w=10;x=5;y=-8;z=11;
A=[a,b;c,d]
B=[w,x;y,z]
det（A）
det（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.
>>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.
symsa11a12a13a21a22a23a31a32a33b1b2b3;
A=[a11,a12,a13;a21,a22,a23;a31,a32,a33];
b=[b1;b2;b3];
XX=A\b%用左除运算求解
（X=linsolve（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*a21*a12-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*a21*a12-a31*a13*a22+a31*a12*a23+a32*a21*a13-a11*a32*a23）
（-a31*b1*a22-a11*a32*b2+a11*b3*a22-b3*a21*a12+a31*a12*b2+a32*a21*b1）/（a11*a33*a22-a33*a21*a12-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）
7.simulink仿真
V.MATLAB在控制系统分析中的应用
1.求下面系统的单位阶跃响应
%程序如下：
num=[4];den=[1,1,4];
step（num,den）
[y,x,t]=step（num,den）;
tp=spline（y,t,max（y））%计算峰值时间
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）
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函数功能
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）'）
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=
1.6541
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）;
holdon
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];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）;
p=
-10.5969+36.2148i
-10.5969-36.2148i
-0.0562
15.已知系统为：
作该系统的nichols曲线。
%程序如下：
n=[1];d=[1,1,0];
nichols（n,d）;
16.已知系统的开环传递函数为：
当k=2时，分别作nichols曲线和波特图。
%程序如下：
num=1;
den=conv（conv（[10],[11]）,[0.51]）;
subplot（1,2,1）;
nichols（num,den）;grid;%nichols曲线
su
升级会员 