1、系統辨識作業方法一 : 遞推最小二乘法NUM=xlsread( 'shuju' , 'B4:B257');u=NUM;NUM1=xlsread( 'shuju' , 'C4:C257');z=NUM1;%直接給出被辨識參數的初始值，取一個充分小的實N=length(u);c0=0.001,0.001,0.001,0.001'向量p0=10A7*eye(4,4);%初始犬態P0也采用直接取方式，取一個充分大的實數單位矩陣E=0.00000005;%相對誤差E 參考值c=c0,zeros(4,253);%被辨識參數矩陣的初始值

2、及大小e=zeros(4,254);%相對誤差的初始值及大小%開始遞推運算for k=3:N;h1=-z(k-1),-z(k-2),u(k-1),u(k-2)'%求h1x=h1'*p0*h1+1;x1=inv(x);k1=p0*h1*x1;%求k1c1=c0+k1*(z(k)-h1'*c0);%求ce1=(c1-c0)./c0;%求參數的相對變化e(:,k)=e1;%把當前相對變化的列向量加入誤差矩陣的最后一列c0=c1;%新獲得的參數作為下一次遞推的舊參數c(:,k)=c1;%把當前所辨識參數的c1 列向量加入辨識參數矩陣的最后一列p1=p0-k1*h1'*p

3、0;%求p1 值p0=p1;%把當前值給下次用if norm(e1)<=Ebreak ;%若參數收斂滿足要求，終止計算endend%分離參數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(1);%畫第1 個圖形i=1:254;%橫坐標從1 到 254plot(i,a1, 'k' ,i,a2,'b' ,i,b1,'r',i,b2,'g' ) %畫出a1 ， a2 ， b1 ， b2 的各次辨識結

4、果legend( 'a1' , 'a2' , 'b1' , 'b2');xlabel( 'k');ylabel('辨識參數)；%B注縱軸變量title('最小二乘各次遞推參數估計值)圖形標題a1=c(1,254)a2=c(2,254)b1=c(3,254)b2=c(4,254)figure(2);%®第 2 個圖形i=1:254;%橫坐標從1到254plot(i,ea1, 'k' ,i,ea2,'b' ,i,eb1,'r',i,eb2,

5、9;g' );%出 al , a2 , bl , b2 的各次辨識結果的收斂情況legend( 'ea1', 'ea2', 'eb1', 'eb2');xlabel( 'k' ) %B注橫軸變量ylabel('參數誤差)%標注縱軸變量title('參數的誤差U斂情況')%圖形標題ea1=e(1,254)ea2=e(2,254)eb1=e(3,254)eb2=e(4,254)最小二乘各次遞推參數估計值6000400Qo o0 2o o o2 一一o o o Ho ooo6參數的誤差收斂

6、情況o20o52-最小二乘遞推辨識結果：al = -0.6418 a2 =-0.3451 bl = 0.0015 b2 =-0.0014ea1 =-1.4483e-004ea2 =3.5589e-004eb1 = -0.0020eb2 = -0.0021二、似然法辨識v=randn(254,1);%產生正態分布隨機數V=0；%計算噪聲方差for i=1:254V=V+v(i)*v(i);endV1=V/254;N=xlsread( 'shuju.xls' ) ;%>數據作為輸入輸出矩陣A=N(1 : 254,1,2,3);x=A(：,2);y=A(:,3);u=x'

7、;z=y'%導到輸入輸出數據o1=0.001*ones(6,1);p0=eye(6,6);%M 初值zf(1)=01zf(2)=0.1;vf(2)=01vf(1)=01uf(2)=01uf(1)=01% 迭代計算參數值和誤差值for k=3:254h=-z(k-1);-z(k-2);u(k-1);u(k-2);v(k-1);v(k-2);hf=h;K=p0*hf*inv(hf'*p0*hf+1);p=eye(6,6)-K*hf'*p0;v(k)=z(k)-h'*o1;o=o1+K*v(k);p0=p；o1=o;a1(k)=o(1);a2(k)=o(2);b1(k

8、)=o(3);b2(k)=o(4);d1(k)=o(5);d2(k)=o(6);e1(k)=abs(a1(k)+1.2);e2(k)=abs(a2(k)-0.6);e3(k)=abs(b1(k)-1.0);e4(k)=abs(b2(k)-0.5);e5(k)=abs(d1(k)+1.0);e6(k)=abs(d2(k)-0.2);zf(k)=z(k)-d1(k)*zf(k-1)-d2(k)*zf(k-2);uf(k)=u(k)-d1(k)*uf(k-1)-d2(k)*uf(k-2);vf(k)=v(k)-d1(k)*vf(k-1)-d2(k)*vf(k-2);hf=-zf(k-1);-zf(k

9、-2);uf(k-1);uf(k-2);vf(k-1);vf(k-2);end%若權人'可以在巴字運行也觀測到參數就圖figure(1)k=1:254;plot(k,a1, 'k:',k,a2,'b' ,k,b1,'r',k,b2,'m:' ,k,d1,'g' ,k,d2,'k');xlabel( 'k')ylabel('參數辨識')legend( 'a1=-1.2,', 'a2=0.6', 'b1=1.0',

10、'b2=0.5', 'd1=-1.0', 'd2=0.2');%S標炷title( '遞推極大使然RML');figure(2)k=1:254;plot(k,e1, 'k' ,k,e2,'b' ,k,e3,'r',k,e4,'m' ,k,e5,'g' ,k,e6,'k');xlabel（ 'k'）;ylabel（'誤差過度過程）;title（ 誤差曲線）誤差曲線2.5 I1112 -o20遞推極大使然RML05a

11、1=-1.2, a2=0.6b1-1.0b2=O 石 d1=-1.0d2=0.2-0.5-1.5-2L050100I150I200I250300三、增廣法辨識z=xlsread('shuju', 'B4:B257');u=xlsread('shuju', 'C4:C257');v=randn(1,254);c0=0.001 0.001 0.001 0.001 0.001 0.001 0.001'%直接給出被辨識參數的初始值 , 即一個充分小的實向量p0=10A6*eye(7,7);%直接給出初始狀態P0 ,即一個充分大的實

12、數單位矩陣E=0.00000000005;%相對誤差E=0.000000005c=c0,zeros(7,253);%被辨識參數矩陣的初始值及大小e=zeros(7,254);%相對誤差的初始值及大小for k=3:254;%開始求Kh1=-z(k-1),-z(k-2),u(k-1),u(k-2),v(k),v(k-1),v(k-2)'%為求K(k)作準備x=h1'*p0*h1+1;k1=p0*h1*(inv(x);%Kd1=z(k)-h1'*c0; c1=c0+k1*d1;%辨識參數ce1=c1-c0; e2=e1./c0;%求參數誤差的相對變化e(:,k)=e2;c0

13、=c1;%給下一次用c(:,k)=c1;%把遞推出的辨識參數c 的列向量加入辨識參數矩陣p1=p0-k1*k1'*h1'*p0*h1+1;%find p(k)p0=p1;%給下次用if e<=Ebreak ; %若收斂情況滿足要求，終止計算end %判斷結束end %循環結束%分離變量a1=c(1,:); a2=c(2,:); b1=c(3,:); b2=c(4,:);%分離出a1 、 a2 、 b1 、 b2d1=c(5,:); d2=c(6,:); d3=c(7,:);%分離出d1 、 d2 、 d3ea1=e(1,:); ea2=e(2,:); eb1=e(3,:)

14、; eb2=e(4,:);%分離出a1 、 a2 、 b1 、 b2的收斂情況ed1=e(5,:); ed2=e(6,:); ed3=e(7,:);%分離出d1 、 d2 、 d3 的收斂情況figure(1);%畫第二個圖形i=1:254;plot(i,a1, 'r',i,a2,'r:',i,b1,'b' ,i,b2,'b:',i,d1,'g' ,i,d2,'g:',i,d3, 'g+' ) %畫出各個被辨識參數title( 'Parameter Identificatio

15、nwith Recursive Least Squares Method' ) %標題figure(2);i=1:254;%畫出第三個圖形plot(i,ea1, 'r',i,ea2,'r:',i,eb1,'b' ,i,eb2,'b:',i,ed1,'g' ,i,ed2,'g:' ,i,ed2, 'r+') %畫出各個參數收斂情況title( 'Identification Precision')旃題figure(3); subplot(4,1,1);%0出第四個圖形，第一個子圖i=1:254; plot(i,zs(i),'r'