系統(tǒng)辨識(shí)試驗(yàn)

1、2. 用普通最小二乘法( OLS)法辨識(shí)對(duì)象數(shù)學(xué)模型 選擇的仿真對(duì)象的數(shù)學(xué)模型如下z(k) 1.5z(k 1) 0.7z(k 2) u(k 1) 0.5u(k 2) v(k) 其中，v(k)是服從正態(tài)分布的白噪聲 N(0,1) 。輸入信號(hào)采用 4階 M序列，幅度為 1。選擇 如下形式的辨識(shí)模型z(k) a1z(k 1) a2z(k 2) b1u(k 1) b2u(k 2) v(k) 設(shè)輸入信號(hào)的取值是從 k =1 到 k =16 的 M序列，則待辨識(shí)參數(shù)?LS 為?LS =(H LH L ) 1 H LzL 。其中，被辨識(shí)參數(shù) ?LS 、觀測(cè)矩陣 z L、 HL的表達(dá)式為a1z(3)z(2)z

2、(1)u(2)u(1)?a2 ,LSz(4) ，zL， HLz(3)z(2)u(3)u(2)LSb1b2z(16)z(15)z(14)u(15)u(14)程序框圖如下所示：圖 2 最小二乘一次完成算法程序框圖參考程序：%ols u=-1,1,-1,1,1,1,1,-1,-1,-1,1,-1,-1,1,1; % 系統(tǒng)辨識(shí)的輸入信號(hào)為一個(gè)周期的 M 序列z=zeros(1,16); % 定義輸出觀測(cè)值的長度for k=3:16z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2); %用理想輸出值作為觀測(cè)值endsubplot(3,1,1) % 畫三行一列圖形窗口

3、中的第一個(gè)圖形stem(u) % 畫出輸入信號(hào) u的經(jīng)線圖形 subplot(3,1,2) % 畫三行一列圖形窗口中的第二個(gè)圖形i=1:1:16; % 橫坐標(biāo)范圍是 1到 16，步長為 1plot(i,z) % 圖形的橫坐標(biāo)是采樣時(shí)刻 i, 縱坐標(biāo)是輸出觀測(cè)值 z, 圖形格式為連續(xù)曲線 subplot(3,1,3) % 畫三行一列圖形窗口中的第三個(gè)圖形stem(z),grid on% 畫出輸出觀測(cè)值 z的經(jīng)線圖形，并顯示坐標(biāo)網(wǎng)格 u,z%顯示輸入信號(hào)和輸出觀測(cè)信號(hào)%L=14%數(shù)據(jù)長度HL=-z(2) -z(1) u(2) u(1);-z(3) -z(2) u(3) u(2);-z(4) -z(

4、3) u(4) u(3);-z(5) -z(4) u(5) u(4);-z(6) -z(5) u(6) u(5);-z(7) -z(6) u(7) u(6);-z(8) -z(7) u(8) u(7);-z(9)-z(8) u(9) u(8);-z(10) -z(9) u(10) u(9);-z(11) -z(10) u(11) u(10);-z(12) -z(11) u(12) u(11);-z(13) -z(12) u(13) u(12);-z(14) -z(13) u(14) u(13);-z(15) -z(14) u(15) u(14) % 給樣本矩陣 HL賦值ZL=z(3);z(4)

5、;z(5);z(6);z(7);z(8);z(9);z(10);z(11);z(12);z(13);z(14);z(15);z(16)% 給樣本矩陣 zL賦值%calculating parameters% 計(jì)算參數(shù)c1=HL'*HL; c2=inv(c1); c3=HL'*ZL; c=c2*c3 %計(jì)算并顯示%DISPLAY PARAMETERSa1=c(1), a2=c(2), b1=c(3), b2=c(4) % 從 中分離出并顯示 a1 、a2、 b1、 b2%End注：由于輸出觀測(cè)值沒有任何噪音成分， 所以辨識(shí)結(jié)果也無任何誤差， 同學(xué)們可以在輸 出觀測(cè)值中添加噪音，觀

6、察 ols 的辨識(shí)效果。同時(shí)，可以嘗試增加輸入信號(hào)的數(shù)量，看辨識(shí) 結(jié)果有何變化。實(shí)驗(yàn)二 基于 RLS法的系統(tǒng)辨識(shí)數(shù)字仿真實(shí)驗(yàn)一、實(shí)驗(yàn)?zāi)康?、深入理解系統(tǒng)辨識(shí)中相關(guān)分析法及最小二乘法的相關(guān)內(nèi)容2、學(xué)會(huì)用 Matlab 或 C 語言等進(jìn)行系統(tǒng)辨識(shí)的仿真研究二、實(shí)驗(yàn)設(shè)備裝有相應(yīng)軟件的計(jì)算機(jī)。三、實(shí)驗(yàn)原理1. 考慮如下圖所示的仿真對(duì)象：N(z1)v(k)e(k)1+G(z )u(k)y(k)圖 1 遞推最小二乘算法辨識(shí)z(k)圖中， v(k)是服從N(0,1) 分布的不相關(guān)隨機(jī)噪聲。且111)G(z 1) B(z ) ，N (z 1)D(z )A(z 1) C (z 1)A(z 1) 1 1.5a1z

7、 1 0.7z 2 C(z 1)B(z 1) 1.0z 1 0.5z 2D(z 1) 1選擇上圖所示的辨識(shí)模型。仿真對(duì)象選擇如下的模型結(jié)構(gòu)：z(k) a1z(k 1) a2 z(k 2) b1u(k 1) b2u(k 2) v(k)(2)其中， v(k )是服從正態(tài)分布的白噪聲 N( 0,1) 。輸入信號(hào)采用 4 位移位寄存器產(chǎn)生的 M 序列，幅度為 0.03 。按式 (3)z(k) 1.5z(k 1) 0.7z(k 2) u(k 1) 0.5u(k 2) v(k) (3)構(gòu)造 h ( k) ；加權(quán)陣取單位陣 L I ；利用如下公式計(jì)算 K(k)、?(k)和 P( k) ，計(jì)算各次參數(shù)辨識(shí)的相

8、對(duì)誤差，精度滿足要求后停機(jī)。遞推最小二乘法的推導(dǎo)公式如下：(k 1) (k) k(k 1) y(k 1) xT(k 1) (k)K(k 1) P(k)x(k 1)1 x (k 1)P(k)x(k 1) (4)P(k 1) P(k) k(k 1)x (k 1)P(k)2. 階的辨識(shí)前面所討論的系統(tǒng)辨識(shí)方法， 都是假定模型的階次是已知的， 因此僅僅要求估計(jì)差分方 程的系數(shù)。 但實(shí)際上， 系統(tǒng)的階次是很難被準(zhǔn)確知道的。 因?yàn)閷?duì)階次的理解程度是直接與一 個(gè)線性差分方程的準(zhǔn)確結(jié)構(gòu)有關(guān)的， 所以有關(guān)階次的確定也可以稱為系統(tǒng)結(jié)構(gòu)的確定。 經(jīng)驗(yàn) 指出， 一個(gè)模型的階次不準(zhǔn)， 就可能在控制系統(tǒng)設(shè)計(jì)時(shí)發(fā)生嚴(yán)重問題

9、。 故在辨識(shí)過程中， 模 型的階次是否合適是必須加以檢驗(yàn)的。 一般階的方法中， 常用的有這么幾種： 零極點(diǎn)相消法、 目標(biāo)函數(shù)法和 F 檢驗(yàn)法。下面只介紹其中的目標(biāo)函數(shù)法。當(dāng)我們用不同階的模型給系統(tǒng)的輸入輸出觀測(cè)數(shù)據(jù)進(jìn)行最小二乘擬合時(shí)，會(huì)得到 不同的估計(jì)誤差：N2TJe2(k) ET Ek1因此利用 J 極小化確定階是很自然的。實(shí)驗(yàn)表明，假設(shè)模型具有大于 1 而小于 Nmax 的 階 n，當(dāng) n 取 1,2, 時(shí) , 若隨著 n 的增加 , 在 n? ( 階的估計(jì)量 )-1 時(shí)， J 最后一次出現(xiàn)陡峭 的下降，往后 J就近似地保持不變或者只有微小的下降 (見下圖), 則取n n?。也就是說， 模

10、型階次的確定可以直接依次計(jì)算階次n 1，2，時(shí)的最小二乘估計(jì) n?以及相應(yīng)的損失函數(shù) J ，然后選擇當(dāng) J 下降不明顯時(shí)的階次作為合適的模型階次n ，這種方法也叫確定階的估計(jì)準(zhǔn)則方法，有很廣的應(yīng)用。四、實(shí)驗(yàn)內(nèi)容1. 用遞推最小二乘法（ RLS）法辨識(shí)對(duì)象數(shù)學(xué)模型。2. 對(duì)象階的辨識(shí)。五、實(shí)驗(yàn)要求1. 熟悉系統(tǒng)辨識(shí)中的相關(guān)內(nèi)容。2. 掌握 Matlab 或 C 語言等進(jìn)行系統(tǒng)辨識(shí)仿真研究的一般步驟。3. 實(shí)驗(yàn)前基本應(yīng)完成相關(guān)的編程任務(wù)，實(shí)驗(yàn)時(shí)調(diào)試相應(yīng)程序。4. 修改相應(yīng)參數(shù)與隨機(jī)噪聲幅度，觀察并分析結(jié)果。5. 軟件包人機(jī)界面的開發(fā)與設(shè)計(jì)。 （選做）六、實(shí)驗(yàn)步驟1. 首先要熟悉一下 MATLAB的

11、運(yùn)行環(huán)境 :1）File->New->M-File 打開 M文件編輯窗口2）輸入自己編寫的程序3）點(diǎn)擊 run 按鈕，如果程序出錯(cuò)則調(diào)試程序，如果運(yùn)行正常的話則觀察程序的 運(yùn)行結(jié)果2. 用遞推最小二乘法（ RLS）法辨識(shí)對(duì)象數(shù)學(xué)模型在這個(gè)實(shí)驗(yàn)中，我們采用以下模型進(jìn)行仿真：y=1.5*yk-1-0.7*yk-2+0*uk+1.0*uk-10.5*uk-2+ek(5)圖2 遞推最小二乘法算法辨識(shí)的 Malab程序流程圖其中 uk 是幅值為 1 的 PRBS信號(hào)輸入， ek 是白噪聲，即（ 0， 1）的正態(tài)分布序列，它的方差時(shí)可以調(diào)整的。這個(gè)系統(tǒng)的采樣值y（ k）和 u（k） 作為已知數(shù)

12、據(jù)，采用實(shí)驗(yàn)一的最 小二乘法估計(jì)這個(gè)系統(tǒng)的參數(shù)。面是遞推最小二乘法的流程圖：相關(guān)程序如下：%RLSclear% 清理工作間變量L=15;% M 序列的周期y1=1;y2=1;y3=1;y4=0;% 四個(gè)移位寄存器的輸出初始值for i=1:L;% 開始循環(huán)，長度為 Lx1=xor(y3,y4);% 第一個(gè)移位積存器的輸入是第 3 個(gè)與第 4 個(gè)移位積存器的輸出的 “或”x2=y1;%第二個(gè)移位積存器的輸入是第3 個(gè)移位積存器的輸出x3=y2;%第三個(gè)移位積存器的輸入是第2 個(gè)移位積存器的輸出x4=y3;%第四個(gè)移位積存器的輸入是第3 個(gè)移位積存器的輸出y(i)=y4;% 取出第四個(gè)移位積存器幅

13、值為"0" 和 "1" 的輸出信號(hào)，if y(i)>0.5,u(i)=-0.03;% 如果 M序列的值為 "1" 時(shí),辨識(shí)的輸入信號(hào)取“ -0.03 ” else u(i)=0.03;% 當(dāng) M序列的值為 "0" 時(shí), 辨識(shí)的輸入信號(hào)取“ 0.03 ”end%小循環(huán)結(jié)束y1=x1;y2=x2;y3=x3;y4=x4;% 為下一次的輸入信號(hào)做準(zhǔn)備 end%大循環(huán)結(jié)束，產(chǎn)生輸入信號(hào) ufigure(1);% 第 1 個(gè)圖形stem(u),grid on% 以徑的形式顯示出輸入信號(hào)并給圖形加上網(wǎng)格z(2)=0;z(

14、1)=0;% 取 z 的前兩個(gè)初始值為零for k=3:15;% 循環(huán)變量從 3 到 15z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2);% 給出理想的辨識(shí)輸出采樣信號(hào)end%RLS遞推最小二乘辨識(shí)c0=0.001 0.001 0.001 0.001'% 直接給出被辨識(shí)參數(shù)的初始值 , 即一個(gè)充分小的實(shí)向 量p0=106*eye(4,4);% 直接給出初始狀態(tài) P0，即一個(gè)充分大的實(shí)數(shù)單位矩陣E=0.000000005;% 相對(duì)誤差 E=0.000000005c=c0,zeros(4,14);% 被辨識(shí)參數(shù)矩陣的初始值及大小e=zeros(4,

15、15);% 相對(duì)誤差的初始值及大小for k=3:15; % 開始求 Kh1=-z(k-1),-z(k-2),u(k-1),u(k-2)' x=h1'*p0*h1+1; x1=inv(x); %開始求 K(k)k1=p0*h1*x1;% 求出 K 的值d1=z(k)-h1'*c0; c1=c0+k1*d1;% 求被辨識(shí)參數(shù) ce1=c1-c0;% 求參數(shù)當(dāng)前值與上一次的值的差值e2=e1./c0;% 求參數(shù)的相對(duì)變化e(:,k)=e2; % 把當(dāng)前相對(duì)變化的列向量加入誤差矩陣的最后一列c0=c1;% 新獲得的參數(shù)作為下一次遞推的舊參數(shù)c(:,k)=c1;% 把辨識(shí)參數(shù)

16、c 列向量加入辨識(shí)參數(shù)矩陣的最后一列 p1=p0-k1*k1'*h1'*p0*h1+1;% 求出 p(k) 的值p0=p1;% 給下次用if e2<=E break;% 若參數(shù)收斂滿足要求，終止計(jì)算end% 小循環(huán)結(jié)束end%大循環(huán)結(jié)束c%顯示被辨識(shí)參數(shù)e%顯示辨識(shí)結(jié)果的收斂情況%分離參數(shù)a1=c(1,:); a2=c(2,:); b1=c(3,:); b2=c(4,:); ea1=e(1,:); ea2=e(2,:); eb1=e(3,:); eb2=e(4,:);figure(2);% 第 2 個(gè)圖形i=1:15;% 橫坐標(biāo)從 1 到 15plot(i,a1,'

17、;r',i,a2,':',i,b1,'g',i,b2,':')%畫出 a1，a2，b1，b2 的各次辨識(shí)結(jié)果title('Parameter Identification with Recursive Least Squares Method')%圖形標(biāo)題figure(3); % 第 3 個(gè)圖形i=1:15; % 橫坐標(biāo)從 1 到 15plot(i,ea1,'r',i,ea2,'g',i,eb1,'b',i,eb2,'r:') % 畫出 a1，a2，b1，b2

18、 的各次辨 識(shí)結(jié)果的收斂情況title('Identification Precision') %圖形標(biāo)題注：同樣這個(gè)程序使用的輸出信號(hào)也沒有噪音， 所以辨識(shí)的結(jié)果沒有誤差， 請(qǐng)同學(xué)們?cè)谳敵鲂盘?hào)中加入噪音，再使用 RLS對(duì)其辨識(shí)，觀察辨識(shí)結(jié)果，進(jìn)行分析。2. 同樣采用這個(gè)模型采用 MATLAB或者 C 語言進(jìn)行階的辨識(shí)：y=1.5*yk-1-0.7*yk-2+0*uk+1.0*uk-1 0.5*uk-2+ek其中 uk 是幅值為 1 的 PRBS信號(hào)輸入， ek 是白噪聲，即（ 0， 1）的正態(tài)分布序列， 它的方差時(shí)可以調(diào)整的。這個(gè)系統(tǒng)的采樣值y（ k）和 u（k） 作為已知數(shù)據(jù)，采用實(shí)驗(yàn)一的最小二乘法估計(jì)這個(gè)系統(tǒng)的參數(shù)。令模型的階次分別為n1， 2，3，其 J 值如下表所示。從表中看出，由 n? 1到n? 2 ，J發(fā)生顯著的變化。而由 n? 2到n? 3，J 沒有多大的變化， 故可確定這系統(tǒng)的階次 n? 2辨識(shí)的結(jié)果如下，具體編程由同學(xué)們來實(shí)現(xiàn)。噪聲水平損失函數(shù) Jn1n2n30.0265.8630.0000.1248.1470.9870.9830.5337.8482