MATLAB实验.docx
《MATLAB实验.docx》由会员分享，可在线阅读，更多相关《MATLAB实验.docx（22页珍藏版）》请在冰豆网上搜索。
MATLAB实验
实验一
一、打开Matlab软件，点击
，进入simulinkLibraryBrowser，点击文件---新建----Model
选择建立如下simulink模型图。
修改单位阶跃模块的参数：
steptime:
0（阶跃信号起始时间为0时刻）这个值需要设置一下，其余的值都不用改动。
intitialvalue:
0（初值为0）；Finalvalue:
1（终值为1）；Sampletime：
0（采样时间）（这一项我没太弄明白，变了几个值，得到的结果没什么变化。
）
将求和模块改为+-；
将TransferFcn模块的参数设为num:
[10]，den:
[0.020.310]；
（2）simulation菜单star命令开始仿真；（也可使用那个小黑三角的图标工具）
（3）双击示波器模块观察波形；（可以使用Autoscale工具，就是那个像望远镜的图标，得到最佳观测效果）
二、选择建立如下simulink模型图。
（做书上P139页到145页即可）
最后生成
答案：
先画出simulink模型图，然后全选点击右键弹出菜单选CreateSubsystem,生成Subsystem模块，然后选中Subsystem功能模块，右键弹出菜单选中MaskSubsystem进入Mask编辑窗口，如下图将模块下Subsystem名字改为disk（‘PID’）,点击OK生成PIDController。
实验二
1、时域分析
（1）根据下面传递函数模型：
绘制其单位阶跃响应曲线并在图上读标注出峰值，求出系统的性能指标。
答案：
程序shiyan3_1_1.m（本程序中用到了stepchar函数，需要按照书后附录D自己建立stepchar.m）
function[pos,tr,ts,tp]=stepchar（g0,delta）
[y,t]=step（g0）;
figer=plot（t,y）,gridon
[mp,ind]=max（y）;
dimt=length（t）;
yss=y（dimt）;
pos=100*（mp-yss）/yss
tp=t（ind）
fori=1:
dimt
ify（i）>=yss
tr=t（i）;
break;
end
end;
Tr=tr
fori=1:
length（y）
ify（i）<=（1-delta）*yss|y（i）>=（1+delta）*yss
ts=t（i）;
end
end
Ts=ts
以下是主程序：
>>G=tf（[52530],[16108]）;
step（G）
[pos,tr,ts,tp]=stepchar（G,0.02）
figer=
1.9670e+003
pos=
7.5374
tp=
2.2086
Tr=
1.4356
Ts=
3.6994
（2）已知两个线性定常连续系统的传递函数分别为
和
，绘制它们的单位脉冲响应曲线。
，
答案：
程序（shiyan312.m）
g1=tf（[124],[11054]）;
g2=tf（[32],[272]）;
impulse（g1,g2）%脉冲响应曲线
运行结果：
（3）已知线性定常系统的状态空间模型和初始条件，绘制其零输入响应曲线。
答案：
%教材P184例题。
A=[-0.5572,-0.7814;0.7814,0];B=[0;0];C=[1.9691,6.4493];D=[0];
x0=[1;0];t=0:
0.1:
20;
initial（A,B,C,D,x0,t）,gridon
2、频域分析设线性定常连续系统的传递函数分别为
、
和
，将它们的Bose图绘制在一张图中。
，
，
答案：
G1=tf（[1],[51]）;
G2=tf（[0.3],[514]）;
G3=tf（[0.6],[11]）;
bode（G1,'r-',G2,'b-',G3,'k:
'）
运行结果：
3、根轨迹分析根据下面负反馈系统的开环传递函数，绘制系统根轨迹，并分析系统稳定的K值范围。
答案：
sys=tf（[1],conv（[1,0],conv（[11],[12]）））;
rlocus（sys）%P186
[K,ploes]=rlocfind（sys）
注:
根轨迹图出来后，用鼠标点击竖的虚线和根轨迹的相交点。
零极点在左平面时，系统稳定。
故系统稳定的K值范围是0-6。
selected_point=
0+1.4286i
K=
6.1227
ploes=
-3.0111
0.0055+1.4260i
0.0055-1.4260i
实验三控制系统分析
（2）
4、系统稳定性分析
（1）代数法稳定性判据：
（用求分母多项式的根和routh函数两种方法）
已知系统的开环传递函数为：
试对系统闭环判别其稳定性
答案：
要用routh函数，首先按照书后附录D建立routh.m如下
function[rtab,info]=routh0（den）
info=[];
vec1=den（1:
2:
length（den））;nrT=length（vec1）;
vec2=den（2:
2:
length（den）-1）;
rtab=[vec1;vec2,zeros（1,nrT-length（vec2））];
fork=1:
length（den）-2,
alpha（k）=vec1
（1）/vec2
（1）;
fori=1:
length（vec2）,
a3（i）=rtab（k,i+1）-alpha（k）*rtab（k+1,i+1）;
end
ifsum（abs（a3））==0
a3=polyder（vec2）;
info=[info,'Allelementsinrow',...
int2str（k+2）'arezeros;'];
elseifabs（a3
（1））a3
（1）=1e-6;
info=[info,'Replacedfirstelement;'];
end
rtab=[rtab;a3,zeros（1,nrT-length（a3））];
vec1=vec2;vec2=a3;
end
>>G=zpk（[-2],[0,-1,-20],100）
Zero/pole/gain:
100（s+2）
--------------
s（s+1）（s+20）
>>Gclose=feedback（G,1）%单位负反馈，等同于feedback（G,1,-1）
Zero/pole/gain:
100（s+2）
------------------------
（s+3.101）（s+5）（s+12.9）
>>G0=tf（Gclose）
Transferfunction:
100s+200
--------------------------
s^3+21s^2+120s+200
>>den=[121120200];
>>[rtab,msg]=routh（den）
rtab=1120
210
1200
00
msg=
Allelementsinrow4arezeros;
（2）根轨迹法判断系统稳定性：
已知一个单位负反馈系统开环传递函数为：
试在系统的闭环根轨迹图上选择一点，求出该点的增益及其系统的闭环极点位置，
并判断在该点系统闭环的稳定性。
num=[13];
den=conv（conv（[10],[15]）,conv（[16],[122]））;
G=tf（num,den）;
rlocus（G）
[K,p]=rlocfind（G）
Selectapointinthegraphicswindow
selected_point=
-7.0912+3.6957i
K=
1.2158e+003
p=
-7.1115+3.6800i
-7.1115-3.6800i
2.0754+3.8890i
2.0754-3.8890i
-2.9277
再选一次：
Selectapointinthegraphicswindow
selected_point=
-0.2784+1.0248i
K=
15.9165
p=
-5.5352+0.2986i
-5.5352-0.2986i
-1.3640
-0.2828+1.0292i
-0.2828-1.0292i
（3）Bode图法判断系统稳定性：
已知两个单位负反馈系统的开环传递函数分别为：
用Bode图法判断系统闭环的稳定性。
num1=[2.7];den1=[1540];num2=[2.7];den2=[15-40];
g1=tf（num1,den1）;g2=tf（num2,den2）;
margin（g1）,gridon
[Gm,Pm,Wcg,Wcp]=margin（g1）
Gm=7.4074
Pm=51.7321
Wcg=2.0000
Wcp=0.5783
>>[Gm,Pm,Wcg,Wcp]=margin（g2）
Gm=Inf
Pm=-58.0504
Wcg=NaN
Wcp=0.5346
5、系统能控性、能观性分析
已知连续系统的传递函数模型，
当α分别取－1，0，＋1时，判别系统的能控性与能观性
当α等于-1时
num=[1-1];
den=[1102718];
[a,b,c,d]=tf2ss（num,den）;
tc=ctrb（a,b）%求能控矩阵
rank（tc）%能控矩阵的秩
to=obsv（a,c）%求能观矩阵
rank（to）%能观矩阵的秩
执行结果：
tc=
1-1073
01-10
001
ans=
3
to=
01-1
1-10
-11-27-18
ans=
3
能观、能控矩阵都是满秩。
故能观又能控
（2）当α等于0时
num=[10];
den=[1102718];
[a,b,c,d]=tf2ss（num,den）
tc=ctrb（a,b）
rank（tc）
to=obsv（a,c）
rank（to）
执行结果：
tc=1-1073
01-10
001
ans=3
to=010
100
-10-27-18
ans=3
能观、能控矩阵都是满秩。
故能观又能控
（3）当α等于1时
num=[11];
den=[1102718];
[a,b,c,d]=tf2ss（num,den）
tc=ctrb（a,b）%求能控矩阵
rank（tc）%能控矩阵的秩
to=obsv（a,c）%求能观矩阵
rank（to）%能观矩阵的秩
执行结果：
tc=1-1073
01-10
001
ans=3
to=011
110
-9-27-18
ans=2
所以能控，不能观
实验三基于BODE图法的控制系统设计
1、已知单位负反馈控制系统的开环传递函数为：
设计超前校正装置，使校正后系统满足：
，
，
答案：
为留有裕量取
编写Matlab程序：
KK=100;wc=50;Pm=50;
ng0=KK*[1];dg0=conv（[1,0],conv（[0.01,1],[0.1,1]））;
g0=tf（ng0,dg0）;
t=[0:
0.001:
0.5];w=logspace（-1,3）;
[ngc,dgc]=fg_lead_pm_wc（ng0,dg0,Pm,wc,w）;
%[ngc,dgc]=fg_lead_pd（ng0,dg0,wc）
gc=tf（ngc,dgc）,g0c=tf（g0*gc）;
b1=feedback（g0,1）;b2=feedback（g0c,1）;
step（b1,'r--',b2,'b',t）;gridon
Transferfunction:
0.08578s+1
--------------
0.002402s+1
figure,bode（g0,'r--',g0c,'b',w）,gridon,
[pos,tr,ts,tp]=stepchar（b2,0.02）
[gm,pm,wcg,wcp]=margin（g0c）,Km=20*log10（gm）
figer=352.0022
pos=24.9651
tp=0.0409
Tr=0.0289
Ts=0.0939
pos=24.9651
tr=0.0289
ts=0.0939
tp=0.0409
gm=5.8934
pm=44.3060
wcg=201.9342
wcp=69.6644
Km=15.4074
函数程序：
stepchar、fg_lead_pm_wc（带惯性环节的超前控制器）、fg_lead_pd（不带惯）见书后附录，需要同学们先编辑好保存成.m文件。
两个函数均可完成实验要求。
2、已知单位负反馈控制系统的开环传递函数为：
设计超前校正装置，使校正后系统满足：
，
，
答案：
为留有裕量取
编写matlab程序：
KK=100;Pm=40;wc=5;
ng0=KK*[1];dg0=conv（[1,0],conv（[0.1,1],[0.01,1]））;g0=tf（ng0,dg0）;
t=[0:
0.001:
2];w=logspace（-1,4）;
[ngc,dgc]=fg_lag_wc（ng0,dg0,w,wc）;
gc=tf（ngc,dgc）,g0c=tf（g0*gc）;
b1=feedback（g0,1）;b2=feedback（g0c,1）;
step（b1,'r--',b2,'b',t）;axis（[0,5,-1,1.8]）;gridon
Transferfunction:
2s+1
-----------
35.73s+1
figure,bode（g0,'r--',g0c,'b',w）;gridon,
[pos,tr,ts,tp]=stepchar（b2,0.02）
[gm,pm,wcg,wcp]=margin（g0c）
figer=352.0068
pos=16.6273
tp=0.5803
Tr=0.3869
Ts=2.9014
pos=16.6273
tr=0.3869
ts=2.9014
tp=0.5803
gm=18.6322
pm=55.0990
wcg=30.7909
wcp=5.0205
函数程序：
fg_lag_wc见书后附录。
stepchar.m
function[pos,tr,ts,tp]=stepchar（g0,delta）
[y,t]=step（g0）;
[mp,ind]=max（y）;
dimt=length（t）;
yss=y（dimt）;
pos=100*（mp-yss）/yss;
tp=t（ind）;
fori=1:
dimt
ify（i）>=yss
tr=t（i）;
break;
end
end;
fori=1:
length（y）
ify（i）<=（1-delta）*yss|y（i）>=（1+delta）*yss
ts=t（i）;
end
end
fg_lead_pm_wc.m
function[ngc,dgc]=fg_lead_pm_wc（ng0,dg0,Pm,wc,w）
[mu,pu]=bode（ng0,dg0,w）;
ngv=polyval（ng0,j*wc）;dgv=polyval（dg0,j*wc）;
g=ngv/dgv;%求原系统在期望的剪切频率处的频率响应数据G0（jwc）
theta=180*angle（g）/pi;%原系统在期望的剪切频率处的相角裕度，单位为度
alf=ceil（Pm-（theta+180）+5）;%最大超前角
phi=（alf）*pi/180;
a=（1+sin（phi））/（1-sin（phi））;
dbmu=20*log10（mu）;
mm=-10*log10（a）;
wgc=spline（dbmu,w,mm）;
T=1/（wgc*sqrt（a））;
ngc=[a*T,1];dgc=[T,1];
fg_lead_pd.m
function[ngc,dgc]=fg_lead_pd（ng0,dg0,wc）
ngv=polyval（ng0,j*wc）;dgv=polyval（dg0,j*wc）;
g=ngv/dgv;
mg0=abs（g）;
t=sqrt（（（1/mg0）^2-1）/（wc^2））;%幅值相加为零
ngc=[t,1];dgc=[1];%Gc（s）=1+Ts
fg_lag_wc.m
function[ngc,dgc]=fg_lag_wc（ng0,dg0,w,wc）
ngv=polyval（ng0,j*wc）;dgv=polyval（dg0,j*wc）;
g=ngv/dgv;
alph=abs（1/g）;T=10/（alph*wc）;
ngc=[alph*T,1];dgc=[T,1];
3、已知控制系统方框图如图3-1所示。
图中：
试设计反馈校正环节Hc（s）,以满足下列要求：
，
。
，
。
书上p225页例题
KK=50;bp=0.3;ts=2.5;delta=0.05;
ng0=[20];dg0=conv（[1,0],conv（[0.9,1],[0.007,1]））;
g0=tf（ng0,dg0）;
w=logspace（-4,3）;t=[0:
0.1:
3];
[mag,phase]=bode（KK*g0,w）;
[gm0,pm0,wg0,wc0]=margin（mag,phase,w）,gm0=20*log10（gm0）
mr=0.6+2.5*bp'
wc=ceil（（2+1.5*（mr-1）+2.5*（mr-1）^2）*pi/ts）,h=（mr+1）/（mr-1）
w1=2*wc/（h+1）,w2=h*w1
w1=wc/10;w2=9;
ng1=[1/w1,1];dg1=conv（[1/w2,1],conv（[1,0],[1,0]））;
g1=tf（ng1,dg1）;
g=polyval（ng1,j*wc）/polyval（dg1,j*wc）;K=abs（1/g）;
g1=tf（K*g1）
h=tf（dg1,ng1）;Kh=1/K;h=tf（Kh*h）
g2=feedback（KK*g0,h）;
b1=feedback（g0,1）;b2=feedback（g2,1）;
bode（KK*g0,'r--',g2,'b',h,'g',w）;gridon
step（b1,'r--',b2,'b',t）;gridon
[pos,tr,ts,tp]=stepchar（b2.delta）
实验五控制系统的PID设计和工程整定方法
编写一个Ziegler（）函数如下：
function[Gc,Kp,Ti,Td,H]=Ziegler（key,vars）
Ti=[];Td=[];H=[];
iflength（vars）==3
K=vars
（1）;Tc=vars
（2）;N=vars（3）;
ifkey==1,Kp=0.5*K;
elseifkey==2,Kp=0.4*K;Ti=0.85*Tc;
elseifkey==3,Kp=0.6*K;Ti=0.5*Tc;Td=0.125*Tc;
end
end
switchkey
case1
Gc=Kp;
case2
Gc=tf（Kp*[Ti,1],[Ti0]）;
case3
num=[Kp*Ti*Td*（N+1）/N,Kp*（Ti+Td/N）,Kp];
den=Ti*[Td/N,1,0];
Gc=tf（num,den）;
end
调用函数Ziegler（），可设计各个控制器。
程序如下：
G=tf（10,[110355024]）;
[Kc,pp,wg,wp]=margin（G）;
Tc=2*pi/wg;
[Gc1,Kp1]=Ziegler（1,[Kc,Tc,10]）;
disp（'P控制器参数Kp1为：
'）
Kp1
[Gc2,Kp2,Ti2]=Ziegler（2,[Kc,Tc,10]）;
disp（'PI控制器参数Kp2和Ti2为：
'）
Kp2
Ti2
[Gc3,Kp3,Ti3,Td3]=Ziegler（3,[Kc,Tc,10]）;
disp（'PID控制器参数Kp3,Ti3和Td3为：
'）
Kp3
Ti3
Td3
G_c1=feedback（G*Gc1,1）;
G0=feedback（G,1）;
G_c2=feedback（G*Gc2,1）;
G_c3=feedback（G*Gc3,1）;
step（G0,'.',G_c1,'r',G_c2,'g',G_c3,'b'）
实验六基于SIMULINK的控制系统综合仿真
1、建立模型图：
（tu6_1）
运行matlab程序如下：
t=[0:
0.1:
9.9]';
fori=1:
6
ut=[t,i*ones（size（t））];
[tt,xx,yy]=sim（'tu6_1',10,[],ut）;
holdon;
plot（tt,yy）;gridon
end
2、建立模型图：
（tu6_2）
运行matlab程序如下：
t=[0:
20:
2000]';
ut=[t,ones（size（t））];
[tt,xx,yy]=sim（'tu6_2',2000,[],ut）;
plot（tt,yy）;gridon
升级会员 