閉環(huán)系統(tǒng)可辨識性實驗研究

1、研1206 劉新菊 模式識別與智能系統(tǒng) 2012020176實驗三 閉環(huán)系統(tǒng)可辨識性實驗研究實驗目的：通過實驗掌握閉環(huán)系統(tǒng)可辨識性條件及其應(yīng)用。實驗內(nèi)容：1. 模擬如下閉環(huán)系統(tǒng)前向通道方程其中v(k)是服從正態(tài)分布N（0，1）的高斯白噪聲。反饋通道方程1） 二階反饋且存在反饋通道噪聲其中w(k)是服從正態(tài)分布N（0,1）的高斯白噪聲。2） 二階反饋但不加反饋噪聲3） 一階反饋并加反饋噪聲4） 一階反饋不加反饋噪聲5） 比例反饋加反饋噪聲6） 比例反饋不加反饋噪聲2. 在以上各種不同的閉環(huán)狀態(tài)下，利用最小二乘法直接辨識方法估計方法估計前向通道的模型參數(shù)。3. 根據(jù)實驗結(jié)果分析閉環(huán)系統(tǒng)的可辨識性條

2、件。實驗步驟：（1） 熟悉閉環(huán)系統(tǒng)的可辨識性概念及條件。（2） 設(shè)計實驗方案，編制試驗程序，上機調(diào)劑，記錄結(jié)果。（3） 分析實驗結(jié)果，完成書面實驗報告，要求同時提供程序框圖及其說明。實驗三 閉環(huán)系統(tǒng)可辨識性實驗報告一、熟悉閉環(huán)系統(tǒng)的可辨識性概念及條件系統(tǒng)辨識中的閉環(huán)問題作為一類特殊的辨識問題，近年來越來越受到人們的關(guān)注。Soderstrom 指出，由于輸出信號的干擾噪聲通過反饋環(huán)節(jié)與輸入信號相關(guān)，直接采用頻譜分析法，辨識結(jié)果將是對象傳遞函數(shù)與反饋傳遞函數(shù)倒數(shù)間的一個加權(quán)平均值；Ljung 的理論分析表明，如果采用預報誤差法進行開環(huán)辨識時，只要對象模型集包含真實對象的動態(tài)特性

3、，即使噪聲模型不足以描述噪聲的真實動態(tài)特性，仍可獲得對象的一致無偏估計，而將此方法直接用于閉環(huán)辨識時，只要噪聲模型不能精確描述真實噪聲，即使對象模型集包含了真實對象動態(tài)特性，得到的將是對象參數(shù)的有偏估計；Gustavsson 則舉出了反饋環(huán)節(jié)是比例調(diào)節(jié)器時，直接采用預報誤差法將導致對象不可辨識性的實例。  對上述問題的解決，存在三種經(jīng)典的解決方法：直接法、間接法、聯(lián)合輸入輸出法。 (1)直接辨識法。直接法是利用閉環(huán)條件下得到的數(shù)據(jù)，直接把對象當作在開環(huán)環(huán)境中進行辨識。但要保證系統(tǒng)的可辨識性及對象參數(shù)估計的一致無偏性，采用直接法時必須滿足條件：存在足夠的外

4、部激勵信號；控制器的階數(shù)足夠高；控制器在不同的模式間切換；控制器時變、非線性。 (2)間接辨識法。間接法的思想是先根據(jù)閉環(huán)系統(tǒng)的輸入輸出信號用開環(huán)辨識方法得到閉環(huán)系統(tǒng)的傳遞函數(shù)，再由已知的控制器傳遞函數(shù)推導出對象的辨識模型。 (3)聯(lián)合輸入輸出法。將閉環(huán)系統(tǒng)的輸入輸出統(tǒng)一看作為白噪聲驅(qū)動系統(tǒng)產(chǎn)生的輸出信號。先辨識出整個系統(tǒng)的傳遞函數(shù)矩陣，然后利用傳遞函數(shù)矩陣中的分塊矩陣與系統(tǒng)各環(huán)節(jié)傳遞函數(shù)的關(guān)系，求出各環(huán)節(jié)的傳遞函數(shù)。 經(jīng)典閉環(huán)辨識法在許多方面并不盡如人意：直接法雖然簡潔,但必須滿足較為苛刻的條件；間接法與聯(lián)合輸入輸出法實驗前要知道反饋控制器的傳遞函數(shù)，并且需要至

5、少三個可檢測的信號，此外，閉環(huán)系統(tǒng)與開環(huán)對象模型間的參數(shù)轉(zhuǎn)換導致了繁重的計算工作量。1、二階反饋且存在反饋通道噪聲%*%閉環(huán)系統(tǒng)%最小二乘辨識的遞推算法%Z(k)=-1.4*Z(k-1)-0.45*Z(k-2)+u(k-1)+0.7*u(k-2)+v(k)%u(k)=0.33*z(k)+0.033*z(k-1)-0.4*z(k)+w(k) %*clearclc%*產(chǎn)生均值為0，方差為1 的正態(tài)分布噪聲*v=random('Normal',0,1,1,402);w=random('Normal',0,1,1,402);%*產(chǎn)生觀測序列z*z=zeros(402,1

6、);u=zeros(402,1);z(1)=-1'z(2)=0;u(1)=-1'u(2)=0;for i=3:402 z(i)=-1.4*z(i-1)-0.45*z(i-2)+u(i-1)+0.7*u(i-2)+v(i);u(i)=0.33*z(i)+0.033*z(i-1)-0.4*z(i-2)+w(i);end%*遞推求解*P=100*eye(4); %估計方差Pstore=zeros(4,401);Pstore(:,1)=P(1,1),P(2,2),P(3,3),P(4,4);Theta=zeros(4,401); %參數(shù)的估計值，存放中間過程估值Theta(:,1)=3

7、;3;3;3;% K=zeros(4,400); %增益矩陣K=10;10;10;10;for i=3:402 h=-z(i-1);-z(i-2);u(i-1);u(i-2); K=P*h*inv(h'*P*h+1); Theta(:,i-1)=Theta(:,i-2)+K*(z(i)-h'*Theta(:,i-2); P=(eye(4)-K*h')*P; Pstore(:,i-1)=P(1,1),P(2,2),P(3,3),P(4,4);end%*disp('參數(shù)a1、a2、b1、b2 估計結(jié)果：')Theta(:,401)i=1:401;figure

8、(1)plot(i,Theta(1,:),i,Theta(2,:),i,Theta(3,:),i,Theta(4,:)title('待估參數(shù)過渡過程')gridfigure(2)plot(i,Pstore(1,:),i,Pstore(2,:),i,Pstore(3,:),i,Pstore(4,:)title('估計方差變化過程')Grid參數(shù)a1、a2、b1、b2 估計結(jié)果：ans = 1.3842 0.4367 1.0082 0.6828二階反饋且存在反饋通道噪聲模型參數(shù)辨識過程2、 二階反饋但不加反饋噪聲%*%閉環(huán)系統(tǒng)%最小二乘辨識的遞推算法%Z(k)=-1

9、.4*Z(k-1)-0.45*Z(k-2)+u(k-1)+0.7*u(k-2)+v(k)%u(k)=0.33*z(k)+0.033*z(k-1)-0.4*z(k) %*clearclc%*產(chǎn)生均值為0，方差為1 的正態(tài)分布噪聲*v=random('Normal',0,1,1,2002);%w=random('Normal',0,1,1,402);%*產(chǎn)生觀測序列z*z=zeros(2002,1);u=zeros(2002,1);z(1)=-1'z(2)=0;u(1)=-1'u(2)=0;for i=3:2002 z(i)=-1.4*z(i-1)-

10、0.45*z(i-2)+u(i-1)+0.7*u(i-2)+v(i); u(i)=0.33*z(i)+0.033*z(i-1)-0.4*z(i-2);end%*遞推求解*P=100*eye(4); %估計方差Pstore=zeros(4,2001);Pstore(:,1)=P(1,1),P(2,2),P(3,3),P(4,4);Theta=zeros(4,2001); %參數(shù)的估計值，存放中間過程估值Theta(:,1)=3;3;3;3;% K=zeros(4,400); %增益矩陣K=10;10;10;10;for i=3:2002 h=-z(i-1);-z(i-2);u(i-1);u(i-

11、2); K=P*h*inv(h'*P*h+1); Theta(:,i-1)=Theta(:,i-2)+K*(z(i)-h'*Theta(:,i-2); P=(eye(4)-K*h')*P; Pstore(:,i-1)=P(1,1),P(2,2),P(3,3),P(4,4);end%*disp('參數(shù)a1、a2、b1、b2 估計結(jié)果：')Theta(:,2001)i=1:2001;figure(1)plot(i,Theta(1,:),'r',i,Theta(2,:),'r*',i,Theta(3,:),'b'

12、,i,Theta(4,:),'b*')legend('a1','a2','b1','b2')title('待估參數(shù)過渡過程')gridfigure(2)plot(i,Pstore(1,:),i,Pstore(2,:),i,Pstore(3,:),i,Pstore(4,:)title('估計方差變化過程')grid參數(shù)a1、a2、b1、b2 估計結(jié)果：ans = 1.3812 0.4188 0.9284 0.7231二階反饋但不加反饋噪聲模型辨識過程3、 一階反饋并加反饋噪聲%*%閉環(huán)

13、系統(tǒng)%最小二乘辨識的遞推算法%Z(k)=-1.4*Z(k-1)-0.45*Z(k-2)+u(k-1)+0.7*u(k-2)+v(k)%u(k)=z(k)+0.2*z(k-1)+w(k) %*clearclc%*產(chǎn)生均值為0，方差為1 的正態(tài)分布噪聲*v=random('Normal',0,1,1,2002);w=random('Normal',0,1,1,2002);%*產(chǎn)生觀測序列z*z=zeros(2002,1);u=zeros(2002,1);z(1)=-1'z(2)=0;u(1)=-1'u(2)=0;for i=3:2002 z(i)=-

14、1.4*z(i-1)-0.45*z(i-2)+u(i-1)+0.7*u(i-2)+v(i); u(i)=z(i)+0.2*z(i-1)+w(i);end%*遞推求解*P=100*eye(4); %估計方差Pstore=zeros(4,2001);Pstore(:,1)=P(1,1),P(2,2),P(3,3),P(4,4);Theta=zeros(4,2001); %參數(shù)的估計值，存放中間過程估值Theta(:,1)=3;3;3;3;% K=zeros(4,400); %增益矩陣K=10;10;10;10;for i=3:2002 h=-z(i-1);-z(i-2);u(i-1);u(i-2)

15、; K=P*h*inv(h'*P*h+1); Theta(:,i-1)=Theta(:,i-2)+K*(z(i)-h'*Theta(:,i-2); P=(eye(4)-K*h')*P; Pstore(:,i-1)=P(1,1),P(2,2),P(3,3),P(4,4);end%*disp('參數(shù)a1、a2、b1、b2 估計結(jié)果：')Theta(:,2001)i=1:2001;figure(1)plot(i,Theta(1,:),'r',i,Theta(2,:),'r*',i,Theta(3,:),'b',i

16、,Theta(4,:),'b*')legend('a1','a2','b1','b2')title('待估參數(shù)過渡過程')gridfigure(2)plot(i,Pstore(1,:),i,Pstore(2,:),i,Pstore(3,:),i,Pstore(4,:)title('估計方差變化過程')grid參數(shù)a1、a2、b1、b2 估計結(jié)果：ans = 1.4047 0.4392 1.01270.6569一階反饋并加反饋噪聲模型參數(shù)辨識過程4、 一階反饋不加反饋噪聲%*%閉環(huán)系統(tǒng)

17、%最小二乘辨識的遞推算法%Z(k)=-1.4*Z(k-1)-0.45*Z(k-2)+u(k-1)+0.7*u(k-2)+v(k)%u(k)=z(k)+0.2*z(k-1)%*clearclc%*產(chǎn)生均值為0，方差為1 的正態(tài)分布噪聲*v=random('Normal',0,1,1,2002);%w=random('Normal',0,1,1,2002);%*產(chǎn)生觀測序列z*z=zeros(2002,1);u=zeros(2002,1);z(1)=-1'z(2)=0;u(1)=-1'u(2)=0;for i=3:2002 z(i)=-1.4*z(i

18、-1)-0.45*z(i-2)+u(i-1)+0.7*u(i-2)+v(i); u(i)=z(i)+0.2*z(i-1);end%*遞推求解*P=100*eye(4); %估計方差Pstore=zeros(4,2001);Pstore(:,1)=P(1,1),P(2,2),P(3,3),P(4,4);Theta=zeros(4,2001); %參數(shù)的估計值，存放中間過程估值Theta(:,1)=3;3;3;3;% K=zeros(4,400); %增益矩陣K=10;10;10;10;for i=3:2002 h=-z(i-1);-z(i-2);u(i-1);u(i-2); K=P*h*inv(

19、h'*P*h+1); Theta(:,i-1)=Theta(:,i-2)+K*(z(i)-h'*Theta(:,i-2); P=(eye(4)-K*h')*P; Pstore(:,i-1)=P(1,1),P(2,2),P(3,3),P(4,4);end%*%*disp('參數(shù)a1、a2、b1、b2 估計結(jié)果：')Theta(:,2001)i=1:2001;figure(1)plot(i,Theta(1,:),'r',i,Theta(2,:),'r*',i,Theta(3,:),'b',i,Theta(4,:

20、),'b*')legend('a1','a2','b1','b2')title('待估參數(shù)過渡過程')gridfigure(2)plot(i,Pstore(1,:),i,Pstore(2,:),i,Pstore(3,:),i,Pstore(4,:)title('估計方差變化過程')grid參數(shù)a1、a2、b1、b2 估計結(jié)果：ans = 8.4624 1.7891 8.05780.6345一階反饋并加反饋噪聲模型參數(shù)辨識過程5、 比例反饋加反饋噪聲%*%閉環(huán)系統(tǒng)%最小二乘辨識的遞推算

21、法%Z(k)=-1.4*Z(k-1)-0.45*Z(k-2)+u(k-1)+0.7*u(k-2)+v(k)%u(k)=1.2*z(k)+w(k)%*clearclc%=產(chǎn)生均值為0，方差為1 的正態(tài)分布噪聲=v=random('Normal',0,1,1,2002);w=random('Normal',0,1,1,2002);%*產(chǎn)生觀測序列z*z=zeros(2002,1);u=zeros(2002,1);z(1)=-1'z(2)=0;u(1)=-1'u(2)=0;for i=3:2002 z(i)=-1.4*z(i-1)-0.45*z(i-2

22、)+u(i-1)+0.7*u(i-2)+v(i); u(i)=1.2*z(i)+w(i);end%*遞推求解*P=100*eye(4); %估計方差Pstore=zeros(4,2001);Pstore(:,1)=P(1,1),P(2,2),P(3,3),P(4,4);Theta=zeros(4,2001); %參數(shù)的估計值，存放中間過程估值Theta(:,1)=3;3;3;3;% K=zeros(4,400); %增益矩陣K=10;10;10;10;for i=3:2002 h=-z(i-1);-z(i-2);u(i-1);u(i-2); K=P*h*inv(h'*P*h+1); T

23、heta(:,i-1)=Theta(:,i-2)+K*(z(i)-h'*Theta(:,i-2); P=(eye(4)-K*h')*P; Pstore(:,i-1)=P(1,1),P(2,2),P(3,3),P(4,4);end%*%*disp('參數(shù)a1、a2、b1、b2 估計結(jié)果：')Theta(:,2001)i=1:2001;figure(1)plot(i,Theta(1,:),'r',i,Theta(2,:),'r*',i,Theta(3,:),'b',i,Theta(4,:),'b*')l

24、egend('a1','a2','b1','b2')title('待估參數(shù)過渡過程')gridfigure(2)plot(i,Pstore(1,:),i,Pstore(2,:),i,Pstore(3,:),i,Pstore(4,:)title('估計方差變化過程')grid參數(shù)a1、a2、b1、b2 估計結(jié)果：ans = 1.4508 0.4730 1.0414 0.7135比例反饋加反饋噪聲模型參數(shù)辨識過程6、 比例反饋不加反饋噪聲%*%閉環(huán)系統(tǒng)%最小二乘辨識的遞推算法%Z(k)=-1.4*Z(k

25、-1)-0.45*Z(k-2)+u(k-1)+0.7*u(k-2)+v(k)%u(k)=1.2*z(k)%*clearclc%*產(chǎn)生均值為0，方差為1 的正態(tài)分布噪聲*v=random('Normal',0,1,1,2002);w=random('Normal',0,1,1,2002);%*產(chǎn)生觀測序列z*z=zeros(2002,1);u=zeros(2002,1);z(1)=-1'z(2)=0;u(1)=-1'u(2)=0;for i=3:2002 z(i)=-1.4*z(i-1)-0.45*z(i-2)+u(i-1)+0.7*u(i-2)+

26、v(i); u(i)=1.2*z(i);end%*遞推求解*P=100*eye(4); %估計方差Pstore=zeros(4,2001);Pstore(:,1)=P(1,1),P(2,2),P(3,3),P(4,4);Theta=zeros(4,2001); %參數(shù)的估計值，存放中間過程估值Theta(:,1)=3;3;3;3;% K=zeros(4,400); %增益矩陣K=10;10;10;10;for i=3:2002 h=-z(i-1);-z(i-2);u(i-1);u(i-2); K=P*h*inv(h'*P*h+1); Theta(:,i-1)=Theta(:,i-2)+K*(z(i)-h'*Theta(:,i-2); P=(eye(4)-K*h')*P; Pstore(:,i-1)=P(1,1),P(2,2),P(3,3),P(4,4);end%*disp('參數(shù)a1、a2、b1、b2 估計結(jié)果：')Theta(:,2001)i=1:2001;figure(1)plot(i,Theta(1,:),'r',i,Theta(2,:),'r*',i,Theta(3,:),'b',i,Theta(4,:),'b*')legend