實驗2---相關分析法辨識脈沖響應

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

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

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

4、 %初始化Np=26-1; %M序列長度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)實際輸出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、%補償量c%計算脈沖響應估計值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%計算脈沖響應估計誤差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('互相關函數(shù)','脈沖響應估計值','脈沖響應理論值')end產(chǎn)生白噪聲的函數(shù)：function sig=whitenoise(sigma,len) %白噪聲產(chǎn)生函數(shù)，sigma為均方差，len為白噪聲序列數(shù)據(jù)長度a=65539;M=2147483647;b=100;x(1)=12345;r(1)=x(1)/M; % 第一部分為產(chǎn)生0-1的均勻分布隨機數(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的隨機序列 endendsig = (sig - 12*0.5)*sigma; end 7實驗結果及分析另r1, 白噪聲均方差sigma=0.5,運行命令：sigma=response(1)運行結果如下：得sigma = 0.0416另r3,白噪聲均方差sigma=0.5,運行命令：sigma=response(3,0.5)運行結果如下：得sigma = 0.0373圖像如下：比較發(fā)現(xiàn)r=1和r=3時產(chǎn)生的曲線基本相似,但是脈沖響應誤差在r=3時更小。另r1, 白噪聲均方差sigma=1,運行命令：sigma= response(1,1)運行結果如下：得sigma = 0.0594比較第一個與第三個結果,發(fā)現(xiàn)脈沖響應的估計誤差是隨著輸入白噪聲標準差的增大而增大的，白噪聲標準差越小，對系統(tǒng)的輸出干擾越小.8結論在本次系統(tǒng)辨