现代控制理论实验09

1、dfsndzdfn杨晓丹|现代控制理论| 2015年5月26日目录,卖脸内彖及目的2二、实瞼方象内家21-水笳卖验22.卖脸数据处理31j选取数据32丿去初始值33丿平滑铠34丿去零值35丿去异常值33糸统一次辨识41j 一阶有勺衡对象辨识42丿无勺術对象辨识53丿壽阶有衡对象辨识74丿最、二采比辨识85丿适推最小二乘廉辨识io4糸统二次辨识121丿劳孝汾132）随机法143丿混沌比巧4丿粒子群弘16三、卖验结果及分析191试验数据处理192. 糸统一次辨识213. 糸统二次辫识231j劳孝般232）随机法233丿混沌比264丿粒子群法284结果分析3°、卖验内彖及目的本次实验测试水

2、箱特性。通过从一稳定状态到另一稳定状态得到水箱的阶跃响 应曲线，通过阶跃响应曲线求岀水箱的传递函数。通过本次试验掌握处理数据的方 法和系统的辨识方法。二、卖验方案内彖1水箱实验如图所示，将水箱开关v4调到合适位置，电动阀调到合适开度，等待一段时 间后上水箱水位稳定。增大电动阀开度,不改变v4开度,等待一段时间后,上水 箱水位稳定，记下电动阀开度和上水箱水位的数据。2. 实验数据处理将水箱实验中得到的数据按以下步骤处理，得到阶跃响应曲线。1) 选取数据测量的数据需要从一个稳态到达另一个稳态。通过观察输入u的曲线和输出y的曲线z选择一段从一个稳态到达另一个稳态的数据。2) 去初始值将输出曲线y的值

3、减去开始时稳态的值,既将曲线向下平移。3 )平滑化输出信号中含有噪声，用五值一平均的方法把噪声去掉。y(t) = y(t-2)+y(t- l)+y(t)+y(t+l)+y(t+2)/5。4 )去零值将曲线开始部分的0值去除”既曲线向左平移至坐标原点。5)去异常值用差分法检测异常值，若满足ly(/ + l)-y1>21)心)-曲-1)1 ,则认为y(t+l) 点为异常点，2为分辨系数，越大越严格。y(t4-l)=2y(t)-y(t-l)e代码如下clc;clearall;closeall;loadwtank;y=y;sa=314;en=842;y1 (1: en-sa + 1)=y(sa:

4、en);lp,m=size(y1);iflp<m lp=m; end dt=2;y2=y1;fori=3:lp-2y2(i) = (y1(i-2)+y1(i-1)+y1(i)+y1 (i + 1)+y1 (i + 2) )/5;endy2=de_abnormal(y2,10);yo=y2 (1);y3=y2-yo;nzero=l;y4=y3(1+nzero:lp);t(l:lp)=0:lp-l;t=t*dt;tl=t(1:lp-nzero);plot(t,yl);holdon;plot(t1,y4);3系统一次辨识系统的阶次和系统中的独立储能元件的个数有关。上水箱的储能元件个数为, 所以

5、系统阶次应为一阶。水箱的出口流量变化大致为°(严牛=如如，dt a芈丄h所以水箱应为一阶有自衡对象，传递函数应为g(s) = £。dt a ats + l1)一阶有自衡对象辨识1 k 阶有自衡对象的终值y = lim$*_* =k ,时间t为0.632ys时对应的值。代码如下clc;clearall;closeall;loadwtank jy;y1=y;lp,m=size(yl);iflp<mlp=m;endvl, v2, v3, v4, v5, v6f v7 =value (yl, dt); ys=v6;t632=ys*0.632;fori=l:lperr(i)=(

6、yl(i)-ys)*(yl(i)-ys);ende,p=min(err);t=p*dt;k=ys;2)无自衡对象辨识无自衡对象可近似看为g($)二仍：)”,沿着渐近线做斜线，k为该斜线的斜 率，nt二ota=丄*(丝)八2-丄,b点为y周上x0的位置。2 乃 0b6代码如下clc;clearall;closeall;loadsou3;y1=y;lp,m=size(yl);if m>lplp=m;endsum=0;fori=lp:-l:fix(lp*0.95) sum=sum+(yl(i)-yl(i-1)/dt; endk=sum/(lp-fix(lp*0.95); oh=yl(lp)-k

7、*dt*lp; ota=-oh/k;t=ota;ta=fix(ota/dt);ob=y1(ta+1);y=o；yi=o；r=l;n=l;fori=l:lpy=a*y+k*b*r;yl=yl+y*dt;y2(i)=yl;t(i)=i*dt;enderr=zeros(lp,1);fori=l:lperr(i)=4bs(y2(i)-y (i);endaerr=max(err);ifaerr/k>05n=round(1/2/pi*(-oh/ob)*(-oh/ob)-1/6);if n<ln=l;endt=ota/n;aa=exp(-dt/t);bb=l-aa;x=zeros (n + 1,

8、1);fori=l:lpx(1)=a*x(1)+k*b*r;if n>lx(2:n) =aa*x(2:n)+bb*x(l:n-l);endx(n+1)=x(n+1)+x(n)*dt;y2(i)=x(n + l);endend ns=num2str(n);ts=num2str(t);ks=num2str(k);textl=k=2.3,t=16,n=3 ;text2=1k=1,ks,1,t=!,ts,1,n=1,ns; plot(t,y,t,y2)legend (text1,text2,4);3）高阶有衡对象辨识对于z型曲线，可以估计出较为近似的詁而等容模型的k、t、n值。若 终值为ys ,

9、 t04为0.4ys时对应的时间点，t08为0.8ys时对应的时间点,则有7；）4 4-» i/1.07574 cut =_, k=ys , vn = +0.5。2.16s抵讥代码如下clc;clearall;closeall;loadwtank jy;y1=y;lp, m=size(yl);iflp<mlp=m;endvlr v2f v3f v4,v5z v6f v7=value (ylr dt);ys=v6;k=ys;t04=ys*0 4;t08=ys*0.8;e04=zeros(lp,1);e08=zeros(lp,1);fori=l:lpe04(i) = yl(i)-t

10、04*yl (i)-t04;e08 (i) = yl (i)-t08*yl(i)-t08;ende,pos4=min(e04);e,pos8=min(e08);pt04=(pos4-l)*dt;pt08=(pos8-l)*dt;en=zeros(10,1);fori=l:10en(i)=abs(1.075*pt04/(pt08-pt04)+0.5-sqrt(i); endez n=min (en);t=(pt04+pt08)/2.16/n;x=zeros(nz1);r=l;fori=l:lpx(1)=a*x(1)+k*b*r;if n>lfor j=2:nx(j)=a*x(j)+b*x(

11、j-l);endendt(i)=i*dt;y2(i)=x(n);endu(1:lp)=u2-ul;un=0;du=u2-ul;k=k/du;q=qv(y/u,dt,krt/nzun)4) 最小二乘法辨识若y随%,%， xht变化而变化/则有y = %+/?內+ b2xx +»£(<m)>y =/ b =.x =严丿5丿111x11x2nnn)/ 、*；o可写成y = xb + e。采用最小二乘法进行参数估计,则儿二+勺心+亿心+乞«儿=4+勺®+»%+勺 乞=儿一儿若使ek最小,令戶=儿-)+勺弘+化心)2 ,则有ap加ar&quo

12、t;2工以-儿=0 %驾=-2工取 -儿 = °8bn z 带入(1)式得公|jt-1t-1mahahahth工心】b<»+d；ja +【zx皿1念+乏3久=2xi儿 i-l上二 1上二ik=ak=aimah，+xxia=£上二 1jk=lmmam几诫+心心仏日工仏和介+、i-ijt二 1上二i/ a r 111、3、%令丫 =a,b =,x =兀21x2i a x枷丿则式可写成(xtx)b = xty ,贝b = (xtxyxxty可将控制对象写成差分方程的形式亍吋伙-d =亍s心-八则x为y和u/=0 ；=1组成的矢巨阵 y伙)=-ay(k-1)陽丁伙-

13、n) + bu(k -1)bnu(k -n)。代码如下clc;clearall;closeall;loadwtankjy;lpr m=size(y);iflp<mlp=m;endn=l;n=l;u(1:lp)=u2-ul;x=zeros (nf1); t(1:lp)=0:lp-1; t=t *dt;h=zeros(lp-n,2*n);z=zeros(lp-nz1); fori=n+l:lp for j=l:nh(i, j)=-y(i-j);h(i, j + n)=u(i-j);endz(i)=y(i);endxita = inv(hf*h)*hf *z;y1(1:n)=y(1:n);fo

14、ri=n+l:lpfor j=l:nhl(j)=-yl (i-j);hl(j+n)=u(i-j);endy1(i)=hl*xitm;endq=0;fori=l:lpq=q+(yl(i)-y(i) )*(yl(i)-y(i);endq=q/lp;fori=l:nnum(1,i)=xita(n+i);den(1,i)=xita (i);end den= 1 den;sys=tf(num,dena dt);sysc=d2c(sys,1zoh1);numcz dene=tfdata(sysc, 1v1);kt=numc(n+l)/dene(n+l);tt= (dene(1)/dene(n + l)a(

15、1/n);ks=num2str(kt);ts=num2str(tt);ns=num2str(n);qs=num2str(q);textl=k=',ks, ',t= *,ts, 'z n= *,nsz 1,q=',qs;plot (t,y,tzyl);title (textl)5) 递推最小二乘法辨识递推最小二乘法的步骤为(1 )设定遗忘因子久和假定的nk吕星 (m二煎擦匚丈丈(6) m(i + 4)0ss(8)<u(i+y)丈一+2酬七(卜)1+(y)02+y)0la+rx7+(y)0h(i+s0m44 (9)(t+>l)ntc+>img<

16、;( s )q+(i+*)wqa+r)fa+ry+uhu+)4«±(寸)s+e 丄 m+：+q)=+(i+eymqmina+亠富*叵ihz 旺也屋薜菸津(m )(n)、02h(0)£ oh(o)0w強(z)pu 二7-h)亍(n+rq 二riahu) v n:e aoj dt:- + n2ojm*0i oihd 二n*z)aqhi (i sobzheq-hxcn-luepsajtrpn 二dt:o丄dt:t)a 二nlznu (qt:i) npu vemdt uivde 二 a) z-hs 丄 uiqu f arue-umpeot 二hesoto 二tepetof

17、 otok=p*h1*inv(h*p*h +lamda); xita=xita+k*(y(i)-h*xita);p=l/lamda*(i-k*h)*p;endy1(1:n)=y(1:n);fori=n+l:lpfor j=l:nhl(j)=-yl (i-j)；hl(j+n)=u(i-j);endy1(i)=hl*xita;endq=0;fori=l:lpq=q+(yl(i)-y(i)*(yl(i)-y(i);endq=q/lp;fori=l:nnum(1,i)=xita(n+i);den(1,i)=xita (i);endden=1 den;sys=tf(num,den,dt); sysc=d

18、2c(sys, f zoh1);numca dene=tfdata(sysca 1v1 );kt=numc(n+l)/dene(n+l);tt=(dene(1)/dene(n+l)a(1/n);ks=num2str(kt);ts=num2str(tt);ns=num2str(n);qs=num2str(q);text 1=1k=1f ksa ,t=!,tsz 1z n=1,ns, 1,q=1z qs;plot(tz y,t, yl)；title (textl)4.系统二次辨识参考一次辨识时得到的参数k、t、n ,确立二次辨识的辨识范o评价指标q定为预测曲线丫人和原曲线y的方差。评价程序代码如下

19、function q=qv(y,u,dt,k,t,n, un)lp,m=size(y);if m>lplp=m;endq=0;a=exp(-dt/t);b=l-a;x=zeros(nf1);y=o；fori=l:lpx (1) =a*x (1) +k*bt(i);if n>lx(2:n)=a*x(2:n)+b*x(l:n-l);endif un=ly=y+x (n)*dt;y1 (i)=y;elsey1 (i)=x(n);endq=q+(y(i)-yl(i)*(y(i)-yl(i);endq=q/lp;1 )穷举法把k和t可能的取值逐个列出，取q最小时的k和t值。代码如下clc;c

20、learall;closeall;loadwtankjy;lp, m=size(y);iflp<mlp=m;endkmax=18;kmin=10;tmax=370;tmin=220;step=100;qb=10e40;u(1:lp)=u2-ul;khi=;thi=;fori=kmin:(kmax-kmin)/step:kmax for j=tmin:(tmax-tmin)/step:tmax for n=1:2q=qv (y, u, dt, i, j, nz 0);khi=khiz i;thi=thi,j;if qb>qqb=q;kb=i;tb=j;nb=n;endendenden

21、dks=num2str(kb);ts=num2str(tb);ns=num2str(nb);qbs=num2str(qb);text1=1k=1,ks,1,t=1,ts,1,n=1,ns,1,q=1,qbs; plot (khi,thi, f+ f);holdon;plot (kb,tb, 1r+1z 1 linewidth 1z 5);title(text 1)2 ）随机法随机法和穷举法类似,每次产生一个随机的值,q最小时的k和t值。代码如下clc;clearall; closeall; loadwtank jy;lp, m=size (y); iflp<m lp=m;endkmax=

22、18;kmin=10;tmax=37 0;tmin=220; step=10000;qb=10e4 0;u(1:lp)=u2-ul;khi=;thi=;fori=l:stepfor n=l:1k=rand()*(kmax-kmin)+kmin; t=rand()*(tmax-tmin)+tmin; q=qv(y,uzdt,kzt,n,o);khi=khi,k；thi=thi,t;if qb>qqb=q;kb=k;tb=t;nb=n;endendendks=num2str(kb);ts=num2stt(tb);ns=num2str(nb);qbs=num2str(qb);text1=1k=

23、1,ks, 1,t=1,ts, 1,n= f,ns, 1z q=1,qbs; plot (khi,thi, 1 + 1);holdon;plot (kb,tb, 1r+ 1 , 1 linewidth 1z 5);title (text1)3 )混沌法混沌法与随机法类似，通过函数x=4*x*(l-x)来产生无规律的数据,小时的k和t值。代码如下clc;clearall; closeall; loadwtankjy;ta=clock;lp, m=size (y); iflp<m lp=m;endkmax=18;kmin=10;tmax=370;tmin=220;step=10000;qb=

24、10e40;u(1:lp)=u2-ul;khi=;thi=j;templ=rand();temp2=rand();fori=l:stepfor n=l:ltempl=40*templ*(1-templ);temp2=40*temp2*(l-temp2); k=templ*(kmax-kmin)+kmin; t=temp2*(tmax-tmin)+tmin; q=qv(yzu/dt/kzt,nzo)；khi=khi,k;thi=thi,t;if qb>qqb=q;kb=k;tb=t;nb=n;endendendtimel=etime(clock,ta);ks=num2str(kb);ts=

25、num2str(tb);ns=num2str(nb);qbs=num2str(qb);text 1=1k=1z ks, 1,t= f z ts, 'z n=1z ns, 1,q=1,qbs; plot (khi,thi, 1+ 1)；holdon;plot (kb,tb, 1r+1, 'linewidth 1z 5);title (textl)4）粒子群法先在搜索空间内随机放置n个点.产生随机的移动速度来移动点。每次移动 的移动速度由历史上最佳点、当前最佳点、上一次的移动速度决定。以这样的搜索 方式找到合适的k、t、n值。代码如下 clc;clearall;closeall;l

26、oadwtankjy;lp,m=size(y);iflp<mlp=m;endu(1:lp)=u2-ul;kmax=18;kmin=10;tmax=370;tmin=220;num=50;maxv(1)=(kmax-kmin)/10;maxv(2)=(tmax-tmin)/10;cl=0.7;c2=l4;c3=1.4;fori=l:numpso(i,1)=rand()*(kmax-kmin)+kmin;pso(i, 2)=rand()*(tmax-tmin)+tmin;psov(i,1)=rand()* (kmax-kmin)+kmin;psov (i,2)=rand()*(tmax-tm

27、in)+kmin;endbestpso(1:num)=10e40;solpso=10e40;q (1:num)=10e40;bestq(1:num)=10e40;khi=;thi=j;for n=l:lfor ab=l:200fori=l:numq(i)=qv(yzuzdt,pso(i, 1) zpso(iz2) zn, 0);khi=khi pso(iz1);thi=thi pso(i,2);if q(i)<bestq(i)bestq(i)=q (i);bestpso(ir1)=pso (i,1);bestpsofi,2)=pso(i,2);endendbestp=q (1);fori

28、=2:numif q(i)<bestpbestp=q(i);bestppso(l)=pso(i,1);bestppso(2)=pso(i,2);endendifbestp<solpsosolpso=bestp;solk=bestppso(1);sojlt= bestppso (2);soln=n;endfori=l:numpsov(i,1)=cl*psov(iz1)+c2*rand()*(bestpso(i,1)-pso (i,1)+c3*rand()*(bestppso(1)-pso(i, 1);psov(i,2)=cl*psov(ir 2)+c2*rand()*(bestpso

29、(i,2)-pso(iz 2)+c3*rand()*(bestppso(2)-pso(i, 2);if psov(i,1)>maxv(1)psov(i,l)=maxv(l);endif psov(iz 2)>maxv(2)psov(iz2)=maxv(2);endif psov(izl)<-maxv(l)psov(iz1)=-maxv(1);endif psov(i,2)<maxv(2)psov(iz2)=-maxv(2);endpso(i)=pso(i)+psov(i);if pso (i,1)>kmaxpso(i,1)=kmax;endif pso (i,2)

30、>tmaxpso(i,2)=tmax;endif pso(if1)<-kmaxpso (i,1)=-kmax;endif ps0(i,2)<-tmaxpso(i,2)=-tmax;endendendendks=num2str(solk);ts=num2stt(solt);ns=num2str(soln);qbs=num2str(solpso);text 1=k= j ks, 1r t=!,ts, 1z n=1,ns, 1,q=',qbs; plot (khi,thi, 1 + 1);holdon;plot (solk,sol1r+1z 'linewidth 1

31、z 5);titie (text 1)三.卖验结果及分析l试验数据处理原始数据截取其中一段的数据020040060080010001200平滑化及去除异常点将曲线移动至坐标原点2.系统一次辨识最小二乘法递推最小二乘法辨识方法ktnq阶14.033829410.0364无自衡两点法14.0338295.370410.0368最小二乘法14.6864322.761110.0088递推最小二乘 法14.6864322.761110.0088由于曲线后部较平,无自衡算法计算得k二0,无法得到结果，从响应曲线上也 可以看出系统是有自衡对象。3. 系统二次辨识1）穷举法100*100k=15.04,t=3

32、44.5,n=1,0=0.0063309380 i11111左s至a=aa=芸一耋#蓋芸至8 o o o o o6 4 2 0 8 3 3 3 3 22 ）随机法100001003803403203002802602401213151617181114380360k=14.92361=360.6094 .n=1.0=0.02481340+ +320300280260240220w4-+ + +11+1213+m-14+卄+# +# +1516+i1718380360340+320300280260240220101003803603403203002802601112+午+ +13+ +q14

33、15161718k=15.1938.t=356.3616ln=1.0=0.0076607+ + i_i-13+ 卄+ + +丰+ +3 )混沌法10000k=15.0417 .t=343.119(n=1.0=0.0063798380 i11111r-+k=14.8463.t=335.4139.n=1.0=0.0070801380+360'+340-：+320+300280 -260240+220 10141516+h=±17100k=15.1441 .t=352,1044.n=1 .q=0.00695380+360 卜+ +340320300280260240卜+100+ +

34、44-220w1112380360+4-340320300280260+ +l +i + + i 十131415+i + +古161718k=14.6906.7=329.6293.0=1.0=0.009731+-a4-+4l±l.1516+ + +-j_4±-174）粒子群法10000k=14.998.t=341.7974,0=1.0=0.0062971380360340320300. p.28026024022010+100380360340320300280260240+.+11ii i.11:- +i1 >1 1十j h+ + + +电+禪ihi+ + +x+

35、+12+ + +卄 +卄131415161718k=15.1415 j=350.8671 .n=1.0=0.0068035+ + +4-+ + + + + + + +0羊為+ +14100k=15.0082 7=353.5694.0=1.0=0.010763803601卜 +1 1+1 1+1+ :+ + 340- + + + + + +320+ +卄卄 + +-300+才+ + 尸+甘丄+丄+ +280-+ + + 260+ +% + + +240- .+ + + + + + + +22011j1101112131415161718100k=15.232,t=356.0526.n=1.0=0.0075566380i111 1360. + +-+340 # + +t+ + +320 + + +300-+ 暫 k+ + 姑 + -+村 + +4-280-+260+ + + + + + + + + + + + + 240-+ +-+ +2201 11 +