控制系统设计与仿真实验报告.docx
《控制系统设计与仿真实验报告.docx》由会员分享,可在线阅读,更多相关《控制系统设计与仿真实验报告.docx(32页珍藏版)》请在冰豆网上搜索。
![控制系统设计与仿真实验报告.docx](https://file1.bdocx.com/fileroot1/2023-2/3/516f2dd0-526d-4a1e-bc21-db833cc2ea15/516f2dd0-526d-4a1e-bc21-db833cc2ea151.gif)
控制系统设计与仿真实验报告
控制系统设计与仿真实验报告
学院电气与自动化工程学院
专业自动化
年级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)矫正后系统阶跃响应曲线