實(shí)驗(yàn)2---相關(guān)分析法辨識(shí)脈沖響應(yīng)

1、實(shí)驗(yàn)2 相關(guān)分析法辨識(shí)脈沖響應(yīng)實(shí)驗(yàn)報(bào)告哈爾濱工業(yè)大學(xué)航天學(xué)院控制科學(xué)與工程系專業(yè)： 自動(dòng)化 班級(jí)： 1004102 姓名： 籍洋 日期： 2013 年 10 月 10 日1實(shí)驗(yàn)題目： 相關(guān)分析法辨識(shí)系統(tǒng)脈沖響應(yīng) 2實(shí)驗(yàn)?zāi)康?通過(guò)仿真實(shí)驗(yàn)掌握利用相關(guān)分析法辨識(shí)脈沖響應(yīng)的原理和方法。3實(shí)驗(yàn)主要原理一個(gè)單入單出線性定常系統(tǒng)的動(dòng)態(tài)特性可用它的脈沖響應(yīng)函數(shù)g()來(lái)描述。這樣，只要記錄x(t)、y(t)的值，并計(jì)算它們的互相關(guān)函數(shù)，即可求得脈沖響應(yīng)函數(shù)g()。而在系統(tǒng)有正常輸入的情形下，辨識(shí)脈沖響應(yīng)的原理圖如下圖所示。4實(shí)驗(yàn)對(duì)象或參數(shù)相關(guān)分析法v(k)u(k)z(k)下圖為本實(shí)驗(yàn)的原理框圖。系統(tǒng)的傳遞函數(shù)

2、為,其中；分別為系統(tǒng)的輸入和輸出變量；為測(cè)量白噪聲，服從正態(tài)分布，均值為零，方差為，記作；為系統(tǒng)的脈沖響應(yīng)理論值，為系統(tǒng)脈沖響應(yīng)估計(jì)值，為系統(tǒng)脈沖響應(yīng)估計(jì)誤差。系統(tǒng)的輸入采用M序列（采用實(shí)驗(yàn)1中的M序列即可），輸出受到白噪聲的污染。根據(jù)過(guò)程的輸入和輸出數(shù)據(jù)，利用相關(guān)分析法計(jì)算出系統(tǒng)的脈沖響應(yīng)值，并與系統(tǒng)的脈沖響應(yīng)理論值比較，得到系統(tǒng)脈沖響應(yīng)估計(jì)誤差值，當(dāng)時(shí)，應(yīng)該有。u(k)x(k)y(k)u(k)y(k)1、模擬過(guò)程傳遞函數(shù)，獲得過(guò)程的輸入和輸出數(shù)據(jù)（采樣時(shí)間取1秒）。(1) 慣性環(huán)節(jié)其中，T為慣性環(huán)節(jié)的時(shí)間常數(shù)，K為慣性環(huán)節(jié)的靜態(tài)放大倍數(shù)。若采樣時(shí)間記作，則慣性環(huán)節(jié)的輸出可寫(xiě)成：(2) 傳遞

3、函數(shù)仿真（串聯(lián)） 令，則的表達(dá)框圖為：2、互相關(guān)函數(shù)的計(jì)算其中，r為周期數(shù)，表示計(jì)算互相關(guān)函數(shù)所用的數(shù)據(jù)是從第二個(gè)周期開(kāi)始的，目的是等過(guò)程仿真數(shù)據(jù)進(jìn)入平穩(wěn)狀態(tài)。（可分別令r =1、3，對(duì)比仿真結(jié)果）3、c的補(bǔ)償補(bǔ)償量c應(yīng)取，不能取。因?yàn)槭侵芷诤瘮?shù)，則有，故不能取。4、計(jì)算脈沖響應(yīng)估計(jì)值 脈沖響應(yīng)估計(jì)值 脈沖響應(yīng)估計(jì)誤差 u(k)y(k)5程序框圖結(jié)束 繪圖計(jì)算估計(jì)誤差求理論脈沖響應(yīng)g0求估計(jì)脈沖響應(yīng)G求RmzZ=Y+V求出系統(tǒng)在M序列作用下的輸出y產(chǎn)生M序列和白噪聲v初始化參數(shù)T0,T1,T2,K 開(kāi)始6程序代碼function sigma= response(r)x=0,1,0,1,1,0;

4、 %初始化Np=26-1; %M序列長(zhǎng)度a=1; %振幅t=1;for i=1:Np*(r+1) y(i)=x(6); temp=xor(x(5),x(6); for j=5:-1:1 x(j+1)=x(j); end x(1)=temp; endfor i=1:Np*(r+1) if(y(i)=0) u(i)=a; else u(i)=-a; endend K=120;T1=8.3;T2=6.2;T0=1;K1=K/T1/T2;x(1)=0;y(1)=0;for k=2:Np*(r+1) x(k)=exp(-T0/T1)*x(k-1)+T1*K1*(1-exp(-T0/T1)*u(k-1)+

5、T1*K1*(T1*(exp(-T0/T1)-1)+T0)*(u(k)-u(k-1)/T0; y(k)=exp(-T0/T2)*y(k-1)+T2*(1-exp(-T0/T2)*x(k-1)+T1*(T2*(exp(-T0/T2)-1)+T0)*(x(k)-x(k-1)/T0; %未經(jīng)白噪聲污染的輸出endv=whitenoise(1,length(y); %產(chǎn)生白噪聲z=y+v; %系統(tǒng)實(shí)際輸出for k=1:Np sum=0; for i=Np+1:(r+1)*Np sum=u(i-k)*z(i)+sum; end Rmz(k)=1/(r*Np)*sum;endc=-Rmz(Np-1);

6、%補(bǔ)償量c%計(jì)算脈沖響應(yīng)估計(jì)值for k=1:Np G(k)=Np/(Np+1)*a2*t)*(Rmz(k)+c); g0(k)=K/(T1-T2)*(exp(-k*t/T1)-exp(-k*t/T2);end%計(jì)算脈沖響應(yīng)估計(jì)誤差SUM1=0;SUM2=0;for k=1:Np e(k)=g0(k)-G(k); SUM1=e(k)2+SUM1; SUM2=g0(k)2+SUM2;endsigma=sqrt(SUM1/SUM2);step=0:Np-1;plot(step,Rmz(63),Rmz(1:62);hold on;plot(step,G(63),G(1:62),'r'

7、);plot(step,g0(63),g0(1:62),'g');legend('互相關(guān)函數(shù)','脈沖響應(yīng)估計(jì)值','脈沖響應(yīng)理論值')end產(chǎn)生白噪聲的函數(shù)：function sig=whitenoise(sigma,len) %白噪聲產(chǎn)生函數(shù)，sigma為均方差，len為白噪聲序列數(shù)據(jù)長(zhǎng)度a=65539;M=2147483647;b=100;x(1)=12345;r(1)=x(1)/M; % 第一部分為產(chǎn)生0-1的均勻分布隨機(jī)數(shù)for i=1:12*len x(i+1)=mod(a*x(i)+b,M); r(i+1)=x(i+

8、1)/M;end %-% n = 12;for i = 1 : len sig(i)=0; for j = 1 : n sig(i) = sig(i) + r(n * (i - 1) + j);%第二部分產(chǎn)生正態(tài)分布，方差為sigma的隨機(jī)序列 endendsig = (sig - 12*0.5)*sigma; end 7實(shí)驗(yàn)結(jié)果及分析另r1, 白噪聲均方差sigma=0.5,運(yùn)行命令：sigma=response(1)運(yùn)行結(jié)果如下：得sigma = 0.0416另r3,白噪聲均方差sigma=0.5,運(yùn)行命令：sigma=response(3,0.5)運(yùn)行結(jié)果如下：得sigma = 0.0373圖像如下：比較發(fā)現(xiàn)r=1和r=3時(shí)產(chǎn)生的曲線基本相似,但是脈沖響應(yīng)誤差在r=3時(shí)更小。另r1, 白噪聲均方差sigma=1,運(yùn)行命令：sigma= response(1,1)運(yùn)行結(jié)果如下：得sigma = 0.0594比較第一個(gè)與第三個(gè)結(jié)果,發(fā)現(xiàn)脈沖響應(yīng)的估計(jì)誤差是隨著輸入白噪聲標(biāo)準(zhǔn)差的增大而增大的，白噪聲標(biāo)準(zhǔn)差越小，對(duì)系統(tǒng)的輸出干擾越小.8結(jié)論在本次系統(tǒng)辨