使用最小二乘法法進行系統(tǒng)辨識的兩種方法

1、遞推最小二乘法辨識與仿真現(xiàn)在有如下的辨識仿真對象：w(k)e(k)+y(k)u(k)圖1遞推最小二乘算法辨識G(z-1)N(z-)1)圖中，v(k)是服從N(0,1)分布的不相關(guān)隨機噪聲。且G(z1)=B(7-1)，N(z1)=D(z-1)A(z1)C(z1)A(z1)=11.5a1z1+0.7z2=C(z1)B(z1)=1.0z1+0.5z2D(z1)=1選擇上圖所示的辨識模型。仿真對象選擇如下的模型結(jié)構(gòu)：y(k)+ay(k1)+ay(k2)=bu(k1)+bu(k2)+w(k)1212可得系統(tǒng)模型為：y(k)1.5y(k1)+0.7y(k2)=u(k1)+0.5u(k2)+w(k)遞推最小

2、二乘法的推導(dǎo)公式如下：6(k)=6(k1)+K(k)z(k)he(k)6(k1)P(k)=P-1(k1)+h(k)he(k)-1=IK(k)ht(k)P(k1)K(k)=P(k1)h(k)he(k)P(k1)h(k)+1-1相關(guān)程序如下：%exp053%遞推最小二乘法程序%clear%清理工作間變量L=55;%M序列的周期yl=l;y2=l;y3=l;y4=0;%四個移位寄存器的輸出初始值fori=l:L;%開始循環(huán)，長度為Lxl=xor(y3,y4);%第一個移位積存器的輸入是第3個與第4個移位積存器的輸出的“或”x2=y1;%第二個移位積存器的輸入是第3個移位積存器的輸出x3=y2;%第三

3、個移位積存器的輸入是第2個移位積存器的輸出x4=y3;%第四個移位積存器的輸入是第3個移位積存器的輸出y(i)=y4;%取出第四個移位積存器幅值為0和1的輸出信號，ify(i)0.5,u(i)=-0.03;%如果M序列的值為1時,辨識的輸入信號取“-0.03”elseu(i)=0.03;%當(dāng)M序列的值為0時,辨識的輸入信號取“0.03”end%小循環(huán)結(jié)束yl=xl;y2=x2;y3=x3;y4=x4;%為下一次的輸入信號做準備end%大循環(huán)結(jié)束，產(chǎn)生輸入信號uw=normrnd(0,sqrt(0.1),1,55);%加入白噪聲figure(1);%第1個圖形，偽隨機序列stem(u),grid

4、on%以徑的形式顯示出輸入信號并給圖形加上網(wǎng)格z(2)=0;z(1)=0;%取z的前兩個初始值為零fork=3:55;%循環(huán)變量從3到55z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2)+w(k);%給出理想的辨識輸出采樣信號endc0=0.0010.0010.0010.001,0.001;%直接給出被辨識參數(shù)的初始值,即一個充分小的實向量p0=106*eye(5,5);%直接給出初始狀態(tài)P0,即一個充分大的實數(shù)單位矩陣E=0.000000005;%相對誤差E=0.000000005c=c0,zeros(5,54);%被辨識參數(shù)矩陣的初始值及大小e=ze

5、ros(5,55);%相對誤差的初始值及大小fork=3:55;%開始求Kh1=-z(k-1),-z(k-2),u(k-1),u(k-2),w(k);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ù)ce1=c1-c0;%求參數(shù)當(dāng)前值與上一次的值的差值e2=e1./c0;%求參數(shù)的相對變化e(:,k)=e2;%把當(dāng)前相對變化的列向量加入誤差矩陣的最后一列c0=c1;%新獲得的參數(shù)作為下一次遞推的舊參數(shù)c(:,k)=c1;%把辨識參數(shù)c列向量加入辨識參數(shù)矩陣的最后一列p1=p0-k1

6、*k1*h1*p0*h1+1;%求出p(k)的值p0=p1;%給下次用ife2=Ebreak;%若參數(shù)收斂滿足要求，終止計算end%小循環(huán)結(jié)束end%大循環(huán)結(jié)束c;%顯示被辨識參數(shù)e;%顯示辨識結(jié)果的收斂情況%分離參數(shù)a1=c(1,:);a2=c(2,:);b1=c(3,:);b2=c(4,:);d1=c(5,:);ea1=e(1,:);ea2=e(2,:);eb1=e(3,:);eb2=e(4,:);figure(2);%第2個圖形i=l:55;%橫坐標從1到55plot(i,a1,r,i,a2,:,i,b1,g,i,b2,:,i,b1,k)%畫出a1，a2，b1，b2的各次辨識結(jié)果titl

7、e(系統(tǒng)辨識結(jié)果)圖形標題A(z-1)=11.5az-1+0.7z-2=C(z-1)1B(z-1)=1.0z-1+0.5z-2D(z-1)=1一z-1+0.2z-2那么系統(tǒng)模型機構(gòu)為：y(k)一1.5y(k一1)+0.7y(k一2)=u(k一1)+0.5u(k一2)+w(k)一w(k一1)+0.2w(k一2)相關(guān)程序如下：%exp054.m%遞推最小二乘法編程%clear%清理工作間變量L=55;%M序列的周期yl=l;y2=l;y3=l;y4=0;%四個移位寄存器的輸出初始值fori=1:L;%開始循環(huán)，長度為Lxl=xor(y3,y4);%第一個移位積存器的輸入是第3個與第4個移位積存器的

8、輸出的“或”x2=y1;%第二個移位積存器的輸入是第3個移位積存器的輸出x3=y2;%第三個移位積存器的輸入是第2個移位積存器的輸出x4=y3;%第四個移位積存器的輸入是第3個移位積存器的輸出y(i)=y4;%取出第四個移位積存器幅值為0和T的輸出信號，ify(i)0.5,u(i)=-0.03;%如果M序列的值為1時,辨識的輸入信號取“-0.03”elseu(i)=0.03;%當(dāng)M序列的值為0時,辨識的輸入信號取“0.03”end%小循環(huán)結(jié)束yl=xl;y2=x2;y3=x3;y4=x4;%為下一次的輸入信號做準備end%大循環(huán)結(jié)束，產(chǎn)生輸入信號uw=normrnd(0,sqrt(0.1),1

9、,55);figure(1);%第1個圖形，偽隨機序列stem(u),gridon%以徑的形式顯示出輸入信號并給圖形加上網(wǎng)格z(2)=0;z(1)=0;%取z的前兩個初始值為零fork=3:55;%循環(huán)變量從3到55z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2)+w(k)-w(k-1)+0.2*w(k-2);%給出理想的辨識輸出采樣信號endc0=0.0010.0010.0010.001,0.001;%直接給出被辨識參數(shù)的初始值，即一個充分小的實向量p0=106*eye(5,5);%直接給出初始狀態(tài)P0,即一個充分大的實數(shù)單位矩陣E=0.0000000

10、05;%相對誤差E=0.000000005c=c0,zeros(5,54);%被辨識參數(shù)矩陣的初始值及大小e=zeros(5,55);%相對誤差的初始值及大小fork=3:55;%開始求Kh1=-z(k-1),-z(k-2),u(k-1),u(k-2),w(k),w(k-1),w(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ù)ce1=c1-c0;%求參數(shù)當(dāng)前值與上一次的值的差值e2=e1./c0;%求參數(shù)的相對變化e(:,k)=e2;%把當(dāng)前相對變化的列向量加入誤差矩陣的最后一列c0=c1;%新獲得的參數(shù)作為下一次遞推的舊參數(shù)c(:,k)=c1;%把辨識參數(shù)c列向量加入辨識參數(shù)矩陣的最后一列p1=p0-k1*k1*h1*p0*h1+1;%求出p(k)的值pO=pl;%給下次用fe2=Ebreak;%若參數(shù)收斂滿足要求，終止計算end%小循環(huán)結(jié)束end%大循環(huán)結(jié)束c;%顯示被辨識參數(shù)e;%顯示辨識結(jié)果的收斂情況%分離參數(shù)a1=c(1