系統(tǒng)辨識(shí)大作業(yè)2

班級(jí)：自動(dòng)化1班張青義，獲得系統(tǒng)的脈沖傳遞函數(shù)G(2)和傳遞函數(shù)G(S);應(yīng)用最小二乘辨識(shí)，獲得脈沖響應(yīng)序列8(k);同圖顯示兩種方法的辨識(shí)效果圖；應(yīng)用相關(guān)最小二乘法，擬合對(duì)象的差分方程模型；構(gòu)建時(shí)變對(duì)象，用最小二乘法和帶遺忘因子的最小二乘法，(可以用辨識(shí)工具箱)辨識(shí)模型的參數(shù)，比較兩種方法的辨識(shí)效果差異；答：根據(jù)index2搭建結(jié)構(gòu)框圖：aulKKNOTaaa相關(guān)分析法：利用結(jié)構(gòu)框圖得到UY和toutX×1234567X6由于是二階水箱系統(tǒng)，可以假設(shè)系統(tǒng)的傳遞函數(shù)已知G(s)的結(jié)構(gòu)，用長(zhǎng)除法求得G(s)的s展開式，其系數(shù)等于從g(r)求得的各階矩，然后求G(s)的參數(shù)。得到結(jié)果：在commandwindow中輸入下列指令得到傳遞函數(shù)：最小二乘一次算法相關(guān)參數(shù)Estimate=inv(H'*H)*H'*得到相關(guān)系數(shù)為：Estimate=h=[-z(i-1);-z(i-2);M(i);M(i-1);M(iTheta(:,i-1)=Theta(:,i-2)+K*(z(i)-h'plot(i,Theta(1,:),i,Theta(2,:),i,Theta(3,:),i,Theta(4,:),i,The)第二題設(shè)SISO系統(tǒng)差分方程為(40分)y(k)=-a,y(k-1)-a?Y(k-2)+bu(k-D)+b?μ(k-2)+聲。1)用最小二乘及遞推最小二乘法估計(jì)0;2)用輔助變量法及其遞推算法估計(jì)0;答：1.基本最小二乘法：n=2;%根據(jù)題目，本系統(tǒng)為2階系統(tǒng)Y=Data(n+1:n+N,2);%輸出數(shù)據(jù)矩陣disp('使用最小二乘算法所得結(jié)果如下：)%辨識(shí)參數(shù)數(shù)據(jù)輸出如下：lsa1=thita(1),Isa2=thita(22.遞推最小二乘法：%首先將TXT中的數(shù)據(jù)裝載在Matlab中：whileWbs<1disp('從下列文件中選擇所需加載文件：uy1.txtuy2.txtuy3.txt)Data=fscanf(Wbs,%g%g',[2inf]);%生成數(shù)據(jù)矩陣%遞推最小二乘法n=2;%根據(jù)題目已知該系統(tǒng)階數(shù)為2%為辨識(shí)參數(shù)向量賦初值，這里均取為0.001%直接給出矩陣Pn的初始狀態(tài)POthita=[thita0,zeros(n*2,L)];%需辨識(shí)參數(shù)矩陣的初值及大小fori=n+1:1:L%這里i從第n+1個(gè)數(shù)開始form=1:n*2%輸入矩陣賦值X=FIA1*PO*FIA1+1;%引入中間變量X為K遞推公式中的求逆項(xiàng)，1維PO=P1;%作為下次迭代初值thita1=thitaO+K1*D1;%估值補(bǔ)償公式，thita1為求被辨識(shí)的參數(shù)向量thitaO=thita1;%所得的參數(shù)向量作為下次迭代初值disp('遞推最小二乘法所得結(jié)果如下：)Dta1=thita1(1),Dta2=thita1(2),Dtb1=thita1(3),Dtb2=thita1(4)%顯示被辨識(shí)參數(shù)%首先將TXT中的數(shù)據(jù)裝載在Matlab中：disp('從下列文件中選擇所需加載文件：uy1.txtuy2.txtuy3.txt)Data=fscanf(Wbs,'%g%g’[2inf]);%生成數(shù)據(jù)矩陣%使用遞推輔助變量法對(duì)2階系統(tǒng)進(jìn)行辨識(shí)，取bO=0%分別取出輸入、輸出數(shù)據(jù)u,y%首先使用最小二乘法，粗略得到初始的被辨識(shí)參數(shù)向量Y1=Data(3:L,2);%構(gòu)造輸出向量，也用于輔助變量法的輸出向量%獲得初值，以下使用輔助變量法Yf=fzOL*thita;%輔助模型輸出%構(gòu)造輔助變量矩陣，改變新的矩陣中的輸出量，輸入量不變thita=inv(Zfz'*fzOL)*Zfz'*Y1;%求取被估計(jì)參數(shù)的輔助變量估值disp('使用輔助變量法，得到如下辨識(shí)結(jié)果：)4.遞推輔助變量法：%首先將TXT中的數(shù)據(jù)裝載在Matlab中：disp('從下列文件中選擇所需加載文件：uy1.txtuy2.txtuy3.txt)Data=fscanf(Wbs,%g%g',[2inf]);%生成數(shù)據(jù)矩陣%使用遞推輔助變量法對(duì)2階系統(tǒng)進(jìn)行辨識(shí)，取bO=0%分別取出輸入、輸出數(shù)據(jù)u,y%取前100個(gè)數(shù)據(jù)利用最小二乘辨識(shí)Y1=Data(3:n0,2);%構(gòu)造輸出向量，也用于輔助變量法的輸出向量Z=[-Y(2:n0-1)-Y(1:n0-2)U(2:n0-1)U(1%后面的數(shù)據(jù)計(jì)算，使用遞推輔助變量法Y11=[Y(1);Y(2);Z*thikesy=[-Y(n-1);-Y(n-2disp(遞推輔助變量法辨識(shí)結(jié)果：);%首先將TXT中的數(shù)據(jù)裝載在Matlab中：disp('從下列文件中選擇所需加載文件：uy1.txtuy2.txtuy3.txt)[Wbs,Message]=fopen(filename,'r');%Data=fscanf(Wbs,%g%g',[2inf]);%生成數(shù)據(jù)矩陣%使用廣義最小二乘法，辨識(shí)2階系統(tǒng)的參數(shù)%分別取出輸入、輸出數(shù)據(jù)u,y%首先使用最小二乘法，求系統(tǒng)參數(shù)向量初值Zgy1=Data(3:L,2);%輸出矩陣y%得到初值后，以下使用廣義最小二乘法進(jìn)行辨識(shí)，由殘差代替噪聲誤差E(1)=Y(1);%k=1時(shí)，殘差的值%fl表示由殘差代替噪聲項(xiàng)，而辨識(shí)得到的濾波器的估計(jì)參數(shù)fl=inv(omiga'*omiga)*omiga%得到fl后，計(jì)算經(jīng)過(guò)濾波的u,y數(shù)據(jù)；其中，k=1,2時(shí)：Data(k,1)=U(k)+fl(1)*Udisp('由廣義最小二乘法得到的辨識(shí)結(jié)果如下：)實(shí)驗(yàn)結(jié)果：Data=fscanf(WbsN=500-n:%用前50組數(shù)據(jù)利用1sm估計(jì)出初值phi(:,1)=-y1(n:n+n0-1phi(:,2)=-y1(n-1:n+nO-2phi(:,3)=ul(n:n+nO-1phi(:,4)=ul(n-1:n+n0-2%計(jì)算殘差K2=p20*omega*inv(1+omega'*p20*f1=f0+K2*(e(n+i-2)-p2=p20-K2*omega'*p2p2O=p2;%濾波yy(1:i,2)=y1(n-1:n+iuu(1:i,1)=ul(n:(n+i-1uu(1:i,2)=ul((n-1):(n+i-2%yg(k)=y(k)+f(1)*y(k-1)+f(%ug(k)=u(k)+f(1)*u(k-1)+f(ksai=[-yg(n+i-2)-yg(n+i-3)ug(n+i-2)ug(n+iK1=p10*ksai*inv(1+ksai’*p10*ksai);%K(N+1)值phi(1:i+1,3)=ul(n:n+i-1+phi(1:i+1,4)=ul(n-1:n+i-2disp('遞推廣義最小二乘法辨識(shí)結(jié)果：’);al=seta0(1),a2=seta0(2),bl=seta0實(shí)驗(yàn)結(jié)果：7.夏氏偏差修正法：%首先將TXT中的數(shù)據(jù)裝載在Matlab中：[Wbs,Message]=fopen(filename,'r');%打開文件；Data=fscanf(Wbs,%g%g,[2inf]);%生成數(shù)據(jù)矩陣%使用夏氏修正法，對(duì)2階系統(tǒng)進(jìn)行參數(shù)辨識(shí)%首先使用最小二乘法，求系統(tǒng)參數(shù)向量初值Xsthita=Sgl2*Sgl3;%計(jì)算參數(shù)矩陣%得到初值后，以下使用夏氏修正法進(jìn)行參數(shù)辨識(shí)：thitab0=0;%設(shè)偏差項(xiàng)的偏差初值為0E1=Zgl1-glOL*Xsthita;%計(jì)算殘差Edisp('夏氏修正法)Xsa1=Xsthita(1),Xsa2=Xsthita(2),Xsb1=Xsthita(3),Xsb2=Xsthita(4)實(shí)驗(yàn)結(jié)果：5Xsb2=8.夏氏改良法：%首先將TXT中的數(shù)據(jù)裝載在Matlab中：disp('從下列文件中選擇所需加載文件：uy1.txtuy2.txtuy3.txt)[Wbs,Message]=fopen(filename,'r);%打開文件；Data=fscanf(Wbs,%g%g',[2inf]);%生成數(shù)據(jù)矩陣%最小二乘法，取bO=0n=2;%根據(jù)題目，本系統(tǒng)為2階系統(tǒng)L=length(Data);%辨識(shí)所需數(shù)據(jù)的總長(zhǎng)度N=L-n;%構(gòu)造測(cè)量矩陣的數(shù)據(jù)長(zhǎng)度FIA=zeros(N,2*n);%構(gòu)造測(cè)量矩陣fori=1:N%測(cè)量矩陣賦值Y=Data(n+1:n+N,2);%輸出數(shù)據(jù)矩陣thitaLS=thita;%最小二乘估值disp('使用夏氏改良法辨識(shí)結(jié)果為：)Xga1=thita(1),Xga2=thita(2),Xgb1=thit%首先將TXT中的數(shù)據(jù)裝載在Matlab中：disp('從下列文件中選擇所需加載文件：uy1.txtuy2.txtuy3.txt')Data=fscanf(Wbs,’%g%g',[2inf]);%生成數(shù)據(jù)矩陣%這里選用450組數(shù)據(jù)進(jìn)行遞推迭代n=2;N=450-n;%初始賦值如下：=[000000]’;PO=10^10*ey%循環(huán)迭代FIA=[-Y(n+i)-Y(n+i-1)U(n+i)U(n+K=P0*FIA*inv(1+FIA'*POEQ\*jc3\*hps20\o\al(\s\up9(P0),be)EQ\*jc3\*hps20\o\al(\s\up9(=),l)EQ\*jc3\*hps20\o\al(\s\up9(1:),b)disp('取450組數(shù)據(jù)進(jìn)行遞推夏氏辨識(shí)的結(jié)果為：’);Dxal=beta(1),Dxa2=beta(2),Dxb1=beta(3),Dxb2Dxb2=%首先將TXT中的數(shù)據(jù)裝載在Matlab中：disp('從下列文件中選擇所需加載文件：uy1.txtuy2.txtuy3.txt')Data=fscanf(Wbs,’%g%g',[2inf]);%生成數(shù)據(jù)矩陣%使用增廣矩陣法進(jìn)行2階系統(tǒng)的參數(shù)辨識(shí)fori=1:3*n%首先給辨識(shí)參數(shù)向量賦初值，這里取thita0為0P0=10^10*eye(3*n,3*n);%然后設(shè)定初始狀態(tài)P0為充分大的對(duì)角陣thita=[thita0,zeros(3*n,L)];%設(shè)定參數(shù)向量的初值及大小zs=zeros(L,1);%構(gòu)造初始噪聲矩陣forj=n+1:1:Lhn=[-Y(j-1),-Y(j-2),U(j-1),U(j-2),zs(j-1),zs(j-K=P0*hn*inv(hn'*PO*hn+1);%計(jì)算修正矩陣Kthita0=thitaO+K*[Y(j)-hn'*tP0=PO-K*hn'*P0;%計(jì)算狀態(tài)矩陣P(N)E=abs((thita(:,j-1)-thita0)ifE<=ERRObreak;%根據(jù)設(shè)定的誤差限，判斷何時(shí)跳出EQ\*jc3\*hps24\o\al(\s\up10(end),Zgm)Uy1.txtuy2.txtuy3.txt11.極大似然法%首先將TXT中的數(shù)據(jù)裝載在Matlab中：[Wbs,Message]=fopen(filename,'r');%打開文件；Data=fscanf(Wbs,%g%g,[2inf]);%生成數(shù)據(jù)矩陣mthita=zeros(3*n,1);%給被辨識(shí)參數(shù)矩陣賦初始值%用最小二乘法求初始的被辨識(shí)參數(shù)矩陣OL=[-Y(2:L-1),-Y(1:L-2),U(2:L-thita=inv(OL'*OL)*OL'*Z%得到出之后，使用牛頓—拉布森極大似然法編程%給J的梯度和Hassian矩陣賦初始值Ja1=mthita(1);Ja2=mthita(2);Jby2(i)=-Ja1*Y(i-1)-Ja2*Y(i-2)+Jb1*U(i-1)+Jb2*U(i-2)+Jc1*E(i-1)+Jc2*forj=n+1:1:lEa1(j)=Y(j-1)-[Jc1,Jc2]*[Ea1(j-1),Ea1(j-2)]';Ea2(j)=Y(j-2)-[Jc1,Jc2]Eb1(j)-U(j-1)--U(j-2)-[Jc1,Jc2]*[Eb2(j-1),EEc1(j)=-E(j-1)-[Jc1,Jc2]*[Ec1(j-1),Ec1(j-2)]';Ec2(j)=-E(j-2)-[Jc1,Jc2]Ethita=[Ea1(j),Ea2(j),fd=fd+E(j)*Ethita';%迭代計(jì)算」的梯度矩陣Hassian=Hassian+Ethita'*Ethita;%迭代計(jì)算Hassian矩陣%用牛頓-拉卜森法計(jì)算被辨識(shí)的參數(shù)的估值%跳出迭代的判斷disp('極大似然法(牛頓-拉卜森方法),得到的結(jié)果為：)可用LS法起步，否則沒(méi)有可靠的收斂性；對(duì)數(shù)據(jù)u(k)在z(以上題的結(jié)果為例，進(jìn)行：(10分)第二部分：在遞推過(guò)程中，由于是變化的故而需要不斷計(jì)算；(2)辨識(shí)精度高于最小LS估計(jì)法；(4)參數(shù)估計(jì)時(shí)需構(gòu)造輔助變量矩陣。第四題.系統(tǒng)模型階次的辨識(shí)：(10分)1.按殘差方差定階法(含估計(jì)值誤差方差最小法和F檢驗(yàn)法)和殘差白色定階法：disp(從下列文件中選擇所需加載文件：uy1.txtuy2.txtuy3.txt')