控制系统设计与仿真实验报告.docx
《控制系统设计与仿真实验报告.docx》由会员分享，可在线阅读，更多相关《控制系统设计与仿真实验报告.docx（32页珍藏版）》请在冰豆网上搜索。
控制系统设计与仿真实验报告
控制系统设计与仿真实验报告
学院电气与自动化工程学院
专业自动化
年级2010级
班级3班
姓名
学号
指导教师邓斌
2013年5月2日
一、第一次上机实验
2、采用四阶龙格库塔法求如下二阶系统的在幅值为1脉宽为1刺激下响应的数值解。
，
1）程序代码
T=10;e=0.5;w=10;
dt=0.05;N=T/dt;¨
y0=[0;0];yy=y0;
A=[01;-w^2-2*e*w];B=[00;01];C=[10];U0=[0;w^2];U1=[0;0];
fori=1:
20;
K1=A*yy+B*U0;
K2=A*（yy+dt/2*K1）+B*U0;
K3=A*（yy+dt/2*K2）+B*U0;
K4=A*（yy+dt*K3）+B*U0;
yy=yy+dt/6*（K1+2*K2+2*K3+K4）;
y0=[y0yy];
end
fori=21:
N;
K1=A*yy+B*U1;
K2=A*（yy+dt/2*K1）+B*U1;
K3=A*（yy+dt/2*K2）+B*U1;
K4=A*（yy+dt*K3）+B*U1;
yy=yy+dt/6*（K1+2*K2+2*K3+K4）;
y0=[y0yy];
end
figure
plot（[0:
dt:
T],C*y0,'k.'）;
xlabel（'时间t'）
ylabel（'输出y'）
title（'输出响应曲线'）
gridon
2）输出曲线
如右图1所示。
图1
3、采用四阶龙格库塔法求高阶系统阶单位跃响应曲线的数值解。
，
，
1）程序代码
T=50;dt=0.05;N=50/0.05;
w=10;e=0.5;t=5;
a=[010;001;-w*w/t-（t*w*w+2*e*w）/t（-2*t*e*w+1）/t];
b=[000;000;001];u=[0;0;w^2/t];c=[100];
yy=[0;0;0];y1=yy;
fori=1:
N;
k1=a*y1+b*u;
k2=a*（y1+dt*k1/2）+b*u;
k3=a*（y1+dt/2*k2）+b*u;
k4=a*（y1+dt*k3）+b*u;
y1=y1+dt*（k1+2*k2+2*k3+k4）/6;
yy=[yyy1];
end
figure
plot（[0:
dt:
T],c*yy,'b'）;
xlabel（'时间t'）
ylabel（'输出y'）
title（‘阶跃输入响应曲线'）
gridon
2）输出曲线（图2-2）
图2-2
4、自学OED45指令用法，并求解题2中二阶系统的单位阶跃响应。
1）程序代码
t=[0T];
x0=[00];
[ty]=ode45（'second_order',t,x0,[],e,w）;
figure;
plot（t,y（:
1）,'b.'）;
xlabel（'时间t'）
ylabel（'输出y'）
title（'ode45求取系统阶跃响应曲线'）
grid;
second_order：
2）输出曲线
二、第二次上机实验
1、试用simulink方法解微分方程，并封装模块，输出为
。
得到各状态变量的时间序列，以及相平面上的吸引子。
参数入口为
的值以及
的初值。
（其中
，以及初值分别为
）提示：
模块输入是输出量的微分。
1）simulink
2）输出曲线
2、用simulink搭建PI控制器的控制回路，被控对象传递函数：
，分别分析
（1）、比例系数由小到大以及积分时间由小到大对阶跃响应曲线的影响。
（2）、控制器输出有饱和以及反馈有时滞情况下，阶跃响应曲线的变化。
（3）、主控制回路传递函数为：
，副回路为：
，主回路采用PI控制器，副回路采用P控制器，分析控制系统对主回路以及副回路的阶跃扰动的抑制。
注：
PI控制器表达式为
，串级控制如图所示。
1）simulink
（1）
（2）
（3）
2）输出曲线
（1）设ti=5，改变kc
ti=5;forkc=2:
7
[t,x,y]=sim（'exp2_2',20）;
plot（t,y）;
holdon
end
grid
设kc=5，改变ti
kc=5;
forti=2:
7
[t,x,y]=sim（'exp2_2',20）;
plot（t,y）;
holdon
end
grid
（2）控制器输出有饱和以及反馈有时滞情况下，阶跃响应曲线
（3）串级控制下输出曲线
3、编写S函数模块，实现两路正弦信号的叠加，正弦信号相位差为60度。
1）程序代码及simulink模块
function[sys,x0,str,ts]=sfunction_exp（t,x,u,flag）
switchflag,
case0
[sys,x0,str,ts]=mdlInitializeSizes（）;
case1,
sys=mdlDerivatives（t,x,u）;
case2,
sys=mdlUpdate（t,x,u）;
case3,
sys=mdlOutputs（t,x,u）;
case4,
sys=mdlGetTimeOfNextVarHit（t,x,u）;
case9,
sys=mdlTerminate（t,x,u）;
otherwise
error（['Unhandledflag=',num2str（flag）]）;
end;
function[sys,x0,str,ts]=mdlInitializeSizes（）
sizes=simsizes;
sizes.NumContStates=0;%
sizes.NumDiscStates=0;
sizes.NumOutputs=1;
sizes.NumInputs=1;
sizes.DirFeedthrough=1;
sizes.NumSampleTimes=1;
sys=simsizes（sizes）;
x0=[];str=[];ts=[00];
functionsys=mdlDerivatives（t,x,u）
sys=[];
functionsys=mdlUpdate（t,x,u）
sys=[];
functionsys=mdlOutputs（t,x,u）
b=sin（t-pi/3）;
a=u+b;
sys=a;
functionsys=mdlGetTimeOfNextVarHit（t,x,u）
sampleTime=1;%Example,setthenexthittobeonesecondlater.
sys=t+sampleTime;
functionsys=mdlTerminate（t,x,u）
sys=[]；
2）输出曲线
5、利用触发子系统构建以零阶保持器，实现对正弦信号的采样，并比较不同采用周期下的采样波形。
1）simulink模块
2）输出曲线
（1）T=1s
（2）T=0.1s
（3）T=0.01s
6、若被控对象传递函数为
控制器为
，试用simulink搭建一单位反馈控制系统，分析采用周期T对系统单位阶跃响应的影响。
1）simulink
2）输出曲线
3）程序代码
fori=[0.20.40.60.8]
T=i;z1=exp（-0.1*T）;p1=exp（-T）;K=（1-p1）/（1-z1）;
[t,x,y]=sim（'exp2_6',20）;
switchi
case0.2,
plot（t,y,'r'）;holdon
case0.4,
plot（t,y,'g'）;holdon
case0.6,
plot（t,y,'b'）;holdon
case0.8,
plot（t,y,'k'）;holdon
otherwise
[];
end
end
grid
7、设一单位反馈控制系统，控制器采用PI控制，Kp=200,Ki=10,控制器饱和非线性宽度为2，受控对象为时变模型，由微分方程给出，如下：
求系统单位阶跃响应，并分析不同Kp取值对响应曲线的影响。
1）simulink模块
2）输出曲线
三、第三次上机实验
2、试用至少三种方法，判断一下系统的稳定性:
:
1）程序代码
2）输出结果
3、试产生一周期为5秒，时长为30秒，最大值为1，最小值为0的三角波；得到如下一阶系统在三角波输入下的时间响应曲线。
1）程序代码
clear
num=1;den=[21];
y=[];
fort=0:
0.01:
30%三角波信号
yy=mod（t,5）;
yy=yy/5;
y=[yyy];
end
sys=tf（num,den）;
t=0:
0.01:
30;
figure
（2）
lsim（sys,y,t,'r'）%求取输出曲线
grid
2）输出曲线
5、对如下二阶系统做时域分析，得到阻尼比在0~1之间变化的时候，阶跃响应的上升时间，调节时间，峰值时间，超调量以及衰减比（第一个峰值与稳态值之差与第二个峰值与稳态值之差的比）其中
。
1）程序代码
w=5;num=[w^2];T=0.05;tt=10;N=tt/T;
fore=0:
0.2:
1
den=[12*e*ww*w];
sys=tf（num,den）;
t=0:
0.05:
tt;
title（'单位阶跃响应曲线（e=0-1）'）
step（sys,t）
y=step（sys,t）;
holdon
[maxy,imax]=max（y）;
yss=y（N）;
e
if（abs（yss-1）<0.001）
if（maxy>1）
disp（'峰值时间'）
tp=0.05*imax%峰值时间
ttt=find（y>1）;
tr=0.05*ttt
（1）%上升时间
overshoot=（maxy-1）/1*100%超调量
forn=1:
N
if（y（n）>1.02*yss）
ts=n*0.05;
elseif（y（n）<0.98*yss）
ts=n*0.05;
end
end
end
ts%调节时间
n_max=0;
forn=2:
N
maxy2=y（n）;
if（y（n）>y（n+1）&y（n）>y（n-1））
n_max=n_max+1;
end
if（n_max==2）
break;
end
end
shuaijianbi=（maxy-1）/（maxy2-1）%衰减比
yss%稳态值
else
disp（'´ËʱϵͳÎÞ³¬µ÷'）
fori=1:
1:
num
if（y（i）<0.1*yss）
t1=t（i）;
elseif（y（i）>=yss）%避免运算时间过长，到达最大值就返回
break;
end
end
end
fori=1:
1:
num
if（y（i）<0.9*yss）
t2=t（i）;
elseif（y（i）>=yss）%避免运算时间过长，到达最大值就返回break;
end
end
end
tr=t2-t1%上升时间
forn=1:
N
if（y（n）>1.02*yss）
ts=n*0.05;
elseif（y（n）<0.98*yss）
ts=n*0.05;
end
end
end
ts%%调节时间
yss
end
else
disp（'´此系统不稳定'）
end
end
grid
2）程序输出
e=0此时系统不稳定
e=0.2000有超调
tp=0.7000
tr=0.4500
overshoot=52.6122
ts=3.9500
shuaijianbi=3.6285
yss=1.0000
e=0.4000有超调
tp=0.7500
tr=0.5000、
overshoot=25.3177
ts=1.7000
shuaijianbi=15.4901
yss=1.0000
e=0.6000有超调
tp=0.8500
tr=0.6500
overshoot=9.4535
ts=1.2000
shuaijianbi=111.0836
yss=1.0000
e=0.8000有超调
tp=1.1000
tr=0.9000
overshoot=1.5163
ts=0.8000
shuaijianbi=4.3518e+003
yss=1.0000
e=1此时系统无超调
tr=0.6500
ts=1.2000
yss=1
6、已知开环传递函数如下，1）试用根轨迹方法得到其临界稳定增益。
2）若k=10,试用伯德图方法，判断其稳定性。
1）程序代码
2）输出程序
Transferfunction:
10
-----------------------------
0.2s^3+2.3s^2+3.1s+1
Gm=3.4650；Pm=27.6659；Wcg=3.9370；Wcp=2.0701
7、已知系统开环传递函数如下
试设计一超前校正环节，使得超调量为20%，调节时间为1s。
系统单位斜坡稳态响应误差为10%。
并作出校正前后后的系统单位阶跃响应时域曲线加以比较。
1）程序代码
clear;ts=1;%%系统调节时间设为1s
d1=40;d2=conv（conv（[10],[12]）,[110]）;
sys1=tf（d1,d2）;sys=feedback（sys1,1）;
figure
（1）
step（sys）
grid
hold
kc=5;%%矫正环节增益设为5
sys2=tf（d1*kc,d2）;sysn=feedback（sys2,1）;
step（sysn）%%矫正前阶跃响应曲线
sigma=0.2;%%由超调量为20%得到
zeta=（（log（1/sigma））^2/（（pi）^2+（log（1/sigma））^2））^（1/2）;%超调量计算式
wn=4.5/（zeta*ts）;%%调节时间计算式（误差带=2%）
p=[12*wn*zetawn^2];
r=roots（p）;
s_1=r
（1）;%%理想主导极点
nk1=2;
dk1=conv（conv（[10],[0.51]）,[0.11]）;
ngv=polyval（nk1,s_1）;dgv=polyval（dk1,s_1）;g=ngv/dgv;zetag=angle（g）;
mg=abs（g）;ms=abs（s_1）;zetas=angle（s_1）;
tz=（sin（zetas）-kc*mg*sin（zetag-zetas））/（kc*mg*ms*sin（zetag））;
tp=-（kc*mg*sin（zetas）+sin（zetag+zetas））/（ms*sin（zetag））;
nk=[tz,1];dk=[tp,1];Gc=tf（nk,dk）
sysf=feedback（sys1*kc*Gc,1）;
[ty]=step（sysf）;
figure
（2）
plot（y,t,'k'）
title（'矫正后系统阶跃响应曲线'）
xlabel（'t/s'）;ylabel（'y'）;
grid
2）输出结果及曲线
0.4927s+1
Gc=--------------
0.009068s+1
（A）较之前及加入比例调节后阶跃响应曲线
（B）矫正后系统阶跃响应曲线
升级会员 