现代控制理论实验06.docx
- 文档编号:6261770
- 上传时间:2023-01-04
- 格式:DOCX
- 页数:19
- 大小:423.84KB
现代控制理论实验06.docx
《现代控制理论实验06.docx》由会员分享,可在线阅读,更多相关《现代控制理论实验06.docx(19页珍藏版)》请在冰豆网上搜索。
现代控制理论实验06
杨晓丹 | 现代控制理论 | 2015年4月28日
实验6系统仿真方法
一、实验内容及目的
实验用不同的方法对系统进行仿真,通过比较花费时间和仿真精度来评估各种仿真方法。
通过本次实验学会用单环节离散、整体离散、梯形法、四阶龙格库塔法、高阶龙格库塔法仿真,并了解各种仿真方法的精度。
二、实验方案内容
kp1=0.32,ki1=0.0015,kp2=20,ki2=0.0008,k1=0.93,T1=73.3,k2=2.086,T2=96.1。
将各个环节转换为下图的形式
1.系统的状态空间方程
根据上图写出系统的状态空间方程
a=[0000000-ki1;
ki200-ki2000-kp1*ki2;
kp2*k1/T1k1/T1-1/T1-kp2*k1/T1000-kp1*kp2*k1/T1;
001/T1-1/T10000;
000k2/T2-1/T2000;
00001/T2-1/T200;
000001/T2-1/T20;
0000001/T2-1/T2];
b=[ki1;kp1*ki2;kp1*kp2*k1/T1;0;0;0;0;0;];
c=[00000001];
d=[00000000];
根据
,
进行仿真
2.单环节离散化
由于
,
,
,
,推出
。
泰勒展开后得
,取第一阶
令aa1=exp(-dt/T1),bb1=1-aa1,aa2=exp(-dt/T2),bb2=1-aa2。
对于各个环节有
e1=r-x(8);
x
(1)=x
(1)+ki1*dt*e1;
u1=x
(1)+kp1*e1;
e2=u1-x(4);
x
(2)=x
(2)+e2*ki2*dt;
u2=x
(2)+e2*kp2;
x(3)=aa1*x(3)+k1*bb1*u2;
x(4)=aa1*x(4)+bb1*x(3);
x(5)=aa2*x(5)+bb2*k2*x(4);
x(6)=aa2*x(6)+bb2*x(5);
x(7)=aa2*x(7)+bb2*x(6);
x(8)=aa2*x(8)+bb2*x(7);
3.欧拉法
,
。
泰勒展开后得
,取第一阶
。
4.二阶欧拉法
二阶欧拉法先算出y(k)的导数E1,再用y(k)的导数算出y(k+1)的导数E2,并对两个导数求平均值算出y(k+1/2)的导数E=(E1+E2)/2,相当于取到泰勒展开的第二阶。
5.四阶龙格库塔法
龙格库塔法先算出y(k)的导数E1,再用y(k)的导数E1算出y(k+1/2)的导数E2,再用E2再次算出y(k+1/2)的导数E3,用E3算出y(k+1)的导数E4。
令E=(E1+2*E2+2*E3+E4)/6,相当于取到泰勒展开的第四阶。
6.自己的仿真方法
,
,
离散后如下
e1=r-x(8);
x
(1)=x
(1)+ki1*dt*e1;
u1=x
(1)+kp1*e1;
e2=u1-x(4);
x
(2)=x
(2)+e2*ki2*dt;
u2=x
(2)+e2*kp2;
x(3)=x(3)+(k1*u2-x(3))/T1*dt;
x(4)=x(4)+(x(3)-x(4))/T1*dt;
x(5)=x(5)+(k2*x(4)-x(5))/T2*dt;
x(6)=x(6)+(x(5)-x(6))/T2*dt;
x(7)=x(7)+(x(6)-x(7))/T2*dt;
x(8)=x(8)+(x(7)-x(8))/T2*dt;
7.改变仿真时间,观察各种仿真方法的精度
将仿真时间dt增大和减小,观察仿真精度和仿真花费的时间的变化。
代码如下:
clc;
clearall;
closeall;
st=3000;
dt=10;
dtb=0.001;
lp=fix(st/dt);
lpb=fix(st/dtb);
lpk=fix(lpb/lp);
r=1;
kp1=0.32;
ki1=0.0015;
kp2=20;
ki2=0.0008;
k1=0.93;
T1=73.3;
k2=2.086;
T2=96.1;
a=[0000000-ki1;
ki200-ki2000-kp1*ki2;
kp2*k1/T1k1/T1-1/T1-kp2*k1/T1000-kp1*kp2*k1/T1;
001/T1-1/T10000;
000k2/T2-1/T2000;
00001/T2-1/T200;
000001/T2-1/T20;
0000001/T2-1/T2];
b=[ki1;kp1*ki2;kp1*kp2*k1/T1;0;0;0;0;0;];
c=[00000001];
d=zeros(1,8);
x=zeros(8,1);
t=zeros(lp,1);
y=zeros(lp,1);
%tb=zeros(lpb,1);
%faiM=zeros(8);
%nn=20;
%forii=1:
20
%faiM=faiM+dtb^ii/factorial(ii)*a^(ii-1);
%end
%fai=eye(8)+faiM*a;
%faiM=faiM*b;
%
%fori=1:
lpb
%tb(i)=i*dtb;
%x=fai*x+faiM*r;
%y(i)=c*x;
%end
%savestandard1tby
loadstandard1
x=zeros(8,1);
y1=zeros(lp,1);
aa1=exp(-dt/T1);
bb1=1-aa1;
aa2=exp(-dt/T2);
bb2=1-aa2;
t1=clock;
fori=1:
lp
t(i)=i*dt;
e1=r-x(8);
x
(1)=x
(1)+ki1*dt*e1;
u1=x
(1)+kp1*e1;
e2=u1-x(4);
x
(2)=x
(2)+e2*ki2*dt;
u2=x
(2)+e2*kp2;
x(3)=aa1*x(3)+k1*bb1*u2;
x(4)=aa1*x(4)+bb1*x(3);
x(5)=aa2*x(5)+bb2*k2*x(4);
x(6)=aa2*x(6)+bb2*x(5);
x(7)=aa2*x(7)+bb2*x(6);
x(8)=aa2*x(8)+bb2*x(7);
y1(i)=x(8);
end
time1=etime(clock,t1);
x=zeros(8,1);
y2=zeros(lp,1);
t2=clock;
fori=1:
lp
E=a*x+b*r;
x=x+dt*E;
y2(i)=c*x;
end
time2=etime(clock,t2);
x=zeros(8,1);
y3=zeros(lp,1);
t3=clock;
fori=1:
lp
x0=x;
E1=a*x+b*r;
x=x0+dt*E1;
E2=a*x+b*r;
E=(E1+E2)/2;
x=x0+dt*E;
y3(i)=c*x;
end
time3=etime(clock,t3);
x=zeros(8,1);
y4=zeros(lp,1);
t4=clock;
fori=1:
lp
x0=x;
E1=a*x+b*r;
x=x0+dt/2*E1;
E2=a*x+b*r;
x=x0+dt/2*E2;
E3=a*x+b*r;
x=x0+dt*E3;
E4=a*x+b*r;
E=(E1+2*E2+2*E3+E4)/6;
x=x0+dt*E;
y4(i)=c*x;
end
time4=etime(clock,t4);
x=zeros(8,1);
y5=zeros(lp,1);
t5=clock;
fori=1:
lp
e1=r-x(8);
x
(1)=x
(1)+ki1*dt*e1;
u1=x
(1)+kp1*e1;
e2=u1-x(4);
x
(2)=x
(2)+e2*ki2*dt;
u2=x
(2)+e2*kp2;
x(3)=x(3)+(k1*u2-x(3))/T1*dt;
x(4)=x(4)+(x(3)-x(4))/T1*dt;
x(5)=x(5)+(k2*x(4)-x(5))/T2*dt;
x(6)=x(6)+(x(5)-x(6))/T2*dt;
x(7)=x(7)+(x(6)-x(7))/T2*dt;
x(8)=x(8)+(x(7)-x(8))/T2*dt;
y5(i)=x(8);
end
time5=etime(clock,t5);
sum1=0;
sum2=0;
sum3=0;
sum4=0;
sum5=0;
fori=1:
lp
sum1=sum1+(y(i*lpk)-y1(i))^2;
sum2=sum2+(y(i*lpk)-y2(i))^2;
sum3=sum3+(y(i*lpk)-y3(i))^2;
sum4=sum4+(y(i*lpk)-y4(i))^2;
sum5=sum5+(y(i*lpk)-y5(i))^2;
end
sum1=sum1/lp;
sum2=sum2/lp;
sum3=sum3/lp;
sum4=sum4/lp;
sum5=sum5/lp;
sum1s=num2str(sum1);
sum2s=num2str(sum2);
sum3s=num2str(sum3);
sum4s=num2str(sum4);
sum5s=num2str(sum5);
time1s=num2str(time1);
time2s=num2str(time2);
time3s=num2str(time3);
time4s=num2str(time4);
time5s=num2str(time5);
text0=['standard'];
text1=['t=',time1s,'sum=',sum1s];
text2=['t=',time2s,'sum=',sum2s];
text3=['t=',time3s,'sum=',sum3s];
text4=['t=',time4s,'sum=',sum4s];
text5=['t=',time5s,'sum=',sum5s];
dts=num2str(dt);
testt=['dt=',dts];
plot(tb,y,t,y1,t,y2,t,y3,t,y4,t,y5);
legend(text0,text1,text2,text3,text4,text5,4);
title(testt);
三、实验结果及分析
1.改变dt仿真
图例中从上至下分别为基准线,单环节离散法,欧拉法,梯型法,龙格库塔法,自己的仿真方法。
2.结果分析
可以看出,随着仿真步距dt的增大,仿真花费的时间减小,仿真精度降低。
去掉dt较大时记下的仿真花费时间t较小的几项,则t*dt的值基本不变,既仿真步距和仿真时间成反比。
由于仿真的误差和仿真的阶次有关,对于一阶的单环节离散法、欧拉法以及自己的仿真方法,s/dt为一常数,s为标准差。
对于二阶的梯形法,s/dt^2为一常数,对于四阶的龙格库塔法,s/dt^4为一常数(由于累加存在误差,舍去dt<1之前的值)。
若仿真方法为N阶,则仿真误差的标准差s和dt^n成正比。
仿真时间相同时,仿真精度从大到小排列:
自己的方法、单环节离散法、欧拉法、梯形法、龙格库塔法。
可以推出,仿真精度相同时,仿真花费的时间从多到少排列:
自己的方法、单环节离散法、欧拉法、梯形法、龙格库塔法。
高阶的仿真方法仿真的效率更高。
四、实验中遇到的问题
1.作为比较基准的二十阶仿真y的仿真步距dt也随其他仿真dt的变化而变化,应更改代码,将二十阶仿真y的仿真步距固定。
实验时使用了修改后的代码。
2.dt=10时,整体离散法发散,但精度较低的单环节离散法却没有发散。
3.自己实验时,仿真精度只有10^-23,并且dt<1时受dt影响较小。
可能是MATLAB计算精度不足,但老师课上仿真时,仿真精度达到了10^-28,减小dt却达不到这样的精度。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 现代 控制 理论 实验 06