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

《系統(tǒng)辨識》大作業(yè)

學(xué)號：12051124

班級：自動化1班

姓名：張青

信息與控制工程學(xué)院自動化系

2015-07-11

第一題

模仿indcx2,搭建對象，由相關(guān)分析法，獲得脈沖響應(yīng)序列反伏），由后供），參照講

義,獲得系統(tǒng)的脈沖傳遞函數(shù)G（z）和傳遞函數(shù)G"）；應(yīng)用最小二乘辨識,獲得脈沖響應(yīng)序

列以幻；同圖顯示兩種方法的辨識效果圖：應(yīng)用相關(guān)最小二乘法,擬合對象的差分方程模型;

構(gòu)建時變對象，用最小二乘法和帶遺忘因子的最小二乘法，（可以用辨識工具箱）辨識模型

的參數(shù)，比較兩種方法的辨識效果差異；

答：根據(jù)index2搭建結(jié)構(gòu)框圖：

相關(guān)分析法：利用結(jié)構(gòu)框圖得到UY和tout

VariableEditor-UY田口,XWorkspace-*,n?X

屯蒲啕出◎小■11Base▼Selectdata...▼口一Xffl]r：畫㈤Sel...▼?

田UY<201x2double>Name▲Value

123456用UY<201

ffitout<100

15-60■

25-34.2558

35-14.1494

45-2.3471

5-54.7736

6-53.3802

7-5-2.4902-

其中的U就是題目中要求得出的M序列，根據(jù)結(jié)構(gòu)框圖可知序列的周期是

N.=2"-1=24-1=15。

在commandwindow中輸入下列指令，既可以得到脈沖相應(yīng)序列g(shù)（k）

aa=5;NNPP=15;ts=2;

RR=ones(15)+eye(15);

fori=15:-l:l

UU(16-i,:)=UY(16+i:30+i,l)';

end

YY=[UY(31:45,2)];

GG=RR*UU*YY/[aa*aa*(NNPP+l)*ts];

plot(0:2:29,GG)

holdon

stem(0:2:29,GG,'filled')

最小二乘法建模的響應(yīng)序列

由于是二階水箱系統(tǒng)，可以假設(shè)系統(tǒng)的傳遞函數(shù)為G(s)=b/屏,,已知g?),求

1+qs+生

%,々，4,4

已知G(s)的結(jié)構(gòu)，用長除法求得G⑸的s展開式，其系數(shù)等于從g(r)求得的各階矩,

然后求G(s)的參數(shù)。

得到結(jié)果：

al=

-1.1561

a2=

0.4283

b0=

-0.0028

bl=

0.2961

在commandwindow中輸入下列指令得到傳遞函數(shù)：

?num=[0.2961-0.0028];den=[0.4283-1.15611];sys=tf(num,den)

Transferfunction:

0.2961s-0.0028

0.4283s'2-1.156s+1

最小二乘一次算法相關(guān)參數(shù)

%最小二乘法一次完成算法

M=UY(：,1)；

z=UY(:,2);

H=zeros(100,4);

fori=l:100

H(i,l)=-z(i+l);

H(i,2)=-z(i);

H(i,3)=M(i+l);

H(i,4)=M(i);

end

Estimate=inv(H'*H)*H?*(z(3:102))

考結(jié)束

得到相關(guān)系數(shù)為：Estimate=

-0.7866

0.1388

0.5707

0.3115

帶遺忘因子最小二乘法：

考帶遺忘因子最小二乘法程序

M=UY(:,1)；

Z=UY(:,2);

P=1000*eye(5);%

Theta=zeros(5,200);%

Theta(:,1)=[0;0;0;0;0];

K=zeros(4,400);%

K=[10;10;10;10;10];

lamda=0.99;告遺忘因數(shù)

fori=3:201

h=[-z(i-l);-z(i-2);M(i);M(i-2)];

K=P*h*inv(h'*P*h+landa);

Theta(:,i-l)=Theta(:,i-2)+K*(z(i)-h'*Theta(:zi-2));

P=(eye(5)-K*h')*P/lamda;

end

i=l:200;

figure(1)

plot(i,Theta(1,:),i,Theta(2,:),i,Theta(3,:),i,Theta(4,:),i,Theta(5,:)

)

title「帶遺忘因子最小二乘法，)

grid

考結(jié)束

Estimate可由仿真圖得出，可知兩種方法參數(shù)確定十分接近。

兩種方法的辨識效果差異：

采用相關(guān)分析法和最小二乘法辨識出的脈沖響應(yīng)圖形可看出，用最小二乘法辨識出的圖形更

像脈沖響應(yīng)，他在最后無偏差，衰減為零，其可實現(xiàn)無偏差估計。而相關(guān)分析法圖形雖與相

關(guān)分析的相近，但它最后有偏差。帶遺忘因子的遞推最小二乘辨識的參數(shù)平均值可隨著實際

參數(shù)變化而突變。但他對噪聲比較敏感，只是參數(shù)圍繞參數(shù)實際值上下波動，而辨識出參數(shù)

的平均值接近實際值，而且比其他方法更加接近，并且可以防止數(shù)據(jù)飽和。

第二題

設(shè)SISO系統(tǒng)差分方程為(40分)

a

y(k)=~\-0-ci2y(k-2)+bxu{k-\)+b2u(k-2)+

辨識參數(shù)向量為

o=[q&b\%廣，

輸入輸出數(shù)據(jù)詳見數(shù)據(jù)文件uyOl.txt-uy03.txt。式外為噪聲方差各異的白噪聲或有色噪

聲。

試求解：

1)用最小二乘及遞推最小二乘法估計°;

2)用輔助變量法及其遞推算法估計〃；

3)用廣義最小二乘法及其遞推算法估計夕；

4)用夏氏偏差修正法、夏氏改良法及其遞推算法估計夕；

5)用增廣矩陣法估L。；

6)用極大似然法估1一夕；

7)分析噪聲夕外特性：

答：1.基本最小二乘法：

%首先將TXT中的數(shù)據(jù)裝載在Matlab中：

clear;clc;

Wbs=0;

whileWbs<l

disp。從下列文件中選擇所需加載文件：uyl.txtuy2.txtuy3.txt')

filename=input('輸入所需打開文件名:\'3);

[Wbs,Message]=fopen(filename廳);％打開文件;

end

Data-fscanf(Wbs/%g%g',[2inf]);%生成數(shù)據(jù)矩陣

Data=Data';

%最小二乘法，取b0=0

n=2;%根據(jù)題目，本系統(tǒng)為2階系統(tǒng)

L=length(Data);%辨識所需數(shù)據(jù)的總長度

N=L-n;%構(gòu)造測量矩陣的數(shù)據(jù)長度

FIA=zeros(N,2*n);%構(gòu)造測量矩陣

fori=1:N%測量矩陣賦值

forI=l:n*2

if(l>n)

FIA(i,l)=Data(n+2+i-l,l);

else

FIA(iJ)=-Data(n+i-lz2);

end

end

end

丫=Data(n+l:n+N,2);%輸出數(shù)據(jù)矩陣

thita=inv(FIA'*FIA)*FIA'*Y;%計算參數(shù)矩陣

dispC使用最小二乘算法所得結(jié)果如下：')

%辨識參數(shù)數(shù)據(jù)輸出如下：

Isal=thita(l)Jsa2=thita(2|Jsbl=thita(3),lsb2=thita(4)

實驗結(jié)果：

Uyl.txtuy2.txtuy3.txt

Isal=Isal=Isal=

-1.4555-1.4579-1.4693

lsa2=lsa2=lsa2=

0.85190.85130.8639

Isbl=Isbl=Isbl=

0.42500.73440.7451

lsb2=lsb2=lsb2=

0.56190.57730.5543

2.遞推最小二乘法：

%首先將TXT中的數(shù)據(jù)裝載在Matlab中：

clear;clc;

Wbs=0;

whileWbs<l

disp(,從下列文件中選擇所需加載文件:uyl.txtuy2.txtuy3.txt')

filename=input('輸入所需打開文件名:'/s');

[Wbs,Message]=fopen(filename,T');%打開文件；

end

Data=fscanf(Wbs/%g%g',[2inf]);%生成數(shù)據(jù)矩陣

Data=Data';

%遞推最小二乘法

n=2;%根據(jù)題目已知該系統(tǒng)階數(shù)為2

L=length(Data);

%為辨識參數(shù)向量賦初值，這里均取為0.001

fora=l:l:n*2

thitaO(a,l)=0.001;

end

%直接給出矩陣Pn的初始狀態(tài)P0

A

P0=109*eye(n*2zn*2);

thita=[thitaO,zeros(n*2,L)]:%需辨識參數(shù)矩陣的初值及大小

fori=n+l:l:L%這里i從第n+1個數(shù)開始

form=l:n*2%輸入矩陣賦值

if(m<=n)

FIAl(m,l)=-Data(i-m,2);

else

FIAl(m,l)=Data(i-m+2,l);

end

end

X=FIA1'*PO*FIA1+1;%引入中間變量X為K遞推公式中的求逆項,1維

K1=PO*FIA1/X;

DI=Data(i,2)-FIAl'*thitaO;

Pl=PO-K1*K1'*X;%求出P(K)矩陣

PO=Pl;%作為下次迭代初值

thital=thitaO+Kl*Dl%估值補償公式，thital為求被辨識的參數(shù)向量

thitaO=thital;%所得的參數(shù)向量作為下次迭代初值

end

dispC遞推最小二乘法所得結(jié)果如下：()

Dtal=thital(l),Dta2=thital⑵,Dtbl=thitaM3),Dtb2=thitaM4)%顯示被辨識參數(shù)

實驗結(jié)果：

Uyl.txtuy2.txtuy3.txt

Dtal=Dtal=Dtal=

-1.4555-1.4579-1.4693

Dta2=Dta2=Dta2=

0.85190.85130.8639

Dtbl=Dtbl=Dtbl=

0.42500.73440.7451

Dtb2=Dtb2=Dtb2=

0.56190.57730.5543

3.輔助變量法：

%首先將TXT中的數(shù)據(jù)裝載在Matlab中：

clear;clc;

Wbs=0;

whileWbs<l

dispC從下列文件中選擇所需加載文件：uyl.txtuy2.txtuy3.txt')

filename=input('輸入所需打開文件名

[Wbs,Message]=fopen(filename-r');%打開文件；

end

Data=fscanf(Wbs/%g%g',[2inf]);%生成數(shù)據(jù)矩陣

Data=Data';

%使用遞推輔助變量法對2階系統(tǒng)進行辨識，取bO=O

n=2;L=length(Data);N=L-n;

%分別取出輸入、輸出數(shù)據(jù)u,y

U=Data(l:Lzl);Y=Data(l:L,2);

%首先使用最小二乘法，粗略得到初始的被辨識參數(shù)向量

fzOL=[-Y(2:L-l),-Y(l:L-2),U(2:L-l),U(l:L-2)];

Yl=Data(3:L,2);%構(gòu)造輸出向量，也用于輔助變量法的輸出向量

thita=inv(fzOL'*fzOL)*fzOL'*Yl;%得到初始參數(shù)向量

%獲得初值，以下使用輔助變量法

fori=l:400

Yf=fzOL*thita;%輔助模型輸出

%構(gòu)造輔助變量矩陣，改變新的矩陣中的輸出量，輸入量不變

Zfz=fzOL;

Zfz(2:N,l)=-Yf(l:(N-l),l);

Zfz(3:N,2)=-Yf(l:(N-2),l);

thita=inv(Zfz*fzOL)*2fz，*Yl;%求取被估計參數(shù)的輔助變量估值

end

dispC使用輔助變量法，得到如下辨識結(jié)果:')

Fzal=thita(l),Fza2=thita(2),Fzbl=thita(3)/Fzb2=thita(4)

實驗結(jié)果：

Uyl.txtuy2.txtuy3.txt

Fzal=Fzal=Fzal=

-1.4755-1.4479-1.4652

Fza2=Fza2=Fza2=

0.86930.83870.8620

Fzbl=Fzbl=Fzbl=

0.43040.73020.7441

Fzb2=Fzb2=Fzb2=

0.55890.58210.5554

4.遞推輔助變量法：

%首先將TXT中的數(shù)據(jù)裝載在Matlab中:

cleacclc,

Wbs=0;

whileWbs<l

disp('從下列文件中選擇所需加載文件：uyl.txtuy2.txtuy3.txt')

filename=input(,輸入所需打開文件名:\,3)；

[Wbs,Message]=fopen(filename,T);%打開文件；

end

Data=fscanf(Wbs/%g%g',[2inf]);%生成數(shù)據(jù)矩陣

Data=Data';

%使用遞推輔助變最法對2階系統(tǒng)進行辨識，取b0=0

n=2;L=length(Data);

%分別取出輸入、輸出數(shù)據(jù)u,y

U=Data(l:L,l);Y=Data(l:Lz2);

N=498;

n0=100;

%取前100個數(shù)據(jù)利用最小二乘辨識

FIA=[-Y(2:nO-l),-Y(l:nO-2),U(2:nO-l),U(l:nO-2)J;

Yl=Data⑶nO,2);%構(gòu)造輸出向量，也用于輔助變量法的輸出向量

thita=inv(FIA'*FIA)*FIA'*Yl;%得到初始參數(shù)向量

Z=[-Y(2:nO-l)-Y(l:n0-2)U(2:nO-l)U(l:n0-2)];

thita=inv(Z'*FIA)*Z**Yl;

P0=inv(7'*FIA);

%后面的數(shù)據(jù)計算，使用遞推輔助變量法

forn=nO+l:N

Yll=[Y(l);Y(2);Z*thita];

Zn=[-Yll(n-l)-Yll(n-2)U(n-l)U(n-2)];

Z=[Z;Zn];

kesy=[-Y(n-l);-Y(n-2);U(n-l);U(n-2)];

K=P0*Zn7(l+kesy'*P0*Zn');

Pl=PO-K*kesy'*PO;

P0=Pl;

thita=thita+K*(Y(n)-kesy'*thita);

end

dispC遞推輔助變量法辨識結(jié)果：')；

Dfal=thita(l)/Dfa2=thita(2)..Dfbl=thita(3)/Dfb2=thita(4)

實驗結(jié)果：

Uyl.txtuy2.txtuy3.txt

Dfal=

Dfal=Dfal=

-1.4671-1.4657-1.4617

Dfa2=Dfa2=Dfa2=

0.86240.86250.8559

Dfbl=Dfbl=Dfbl=

0.42910.74660.7400

Dfb2=Dfb2=Dfb2=

0.55520.56850.5538

5.廣義最小二乘法：

%首先將TXT中的數(shù)據(jù)裝載在Matlab中：

clear;clc;

Wbs=0;

whileWbs<l

disp(,從下列文件中選擇所需加載文件:uyl.txtuy2.txtuy3.txt')

filename=inputs輸入所需打開文件名:、,「s');

[Wbs,Message]=fopen(filename「r');%打開文件;

end

Data=fscanf(Wbs；%g%g',[2inf]);%生成數(shù)據(jù)矩陣

Data=Data';

%使用廣義最小二乘法，辨識2階系統(tǒng)的參數(shù)

n=2;

L=length(Data);

N=L-n;

form=l:l:300

%分別取出輸入、輸出數(shù)據(jù)u,y

U=Data(l:L,l);Y=Data(l:L,2);

%首先使用最小二乘法，求系統(tǒng)參數(shù)向量初值

測量矩陣

gyOL=[-Y(2:L-1),-Y(1:L-2),U(2:L-1),U(1:L-2)];%LS

Zgyl-Data(3:L,2);%輸出矩降y

thita=inv(gyOL，*gyOL)*gyOL'*Zgyl;%參數(shù)向量初值

Gyal=thita(l);Gya2=thita(2);Gybl=thita(3);Gyb2=thita(4);

%得到初值后，以卜使用廣義最小二乘法進行辨識，由殘差代替噪聲誤差

11)=丫⑴;％k=l時,殘差的值

E(2)=Y(2)+Gyal*Y(l)-Gybl*U(l);%k=2時，，殘差的值

fori=3:N+n%殘差公式如下：

E(i)=Y(i)+Gyal*Y(i-l)+Gya2*Y(i-2)-Gybl*U(i-l)-Gyb2*U(i-2);

end

El=E(n+l:n+N)';

omiga(l:N,l)=-E(n:n+N-l)';omiga(l:N,2)=-E(n-l:n+N-2)';

%fl表示由殘差代替噪聲項，而辨識得到的濾波器的估計參數(shù)

fl=inv(omiga'*omiga)*omiga'*El;

%得到fl后，計算經(jīng)過濾波的u,y數(shù)據(jù)；其中，k=l,2時：

Data(lzl)=U(l);

Data(l/2)=Y(1);

Data(2/2)=Y(2)+fl(l)*Y(l);

Data(2,l)=U(2)+fl(l)*U(l);

%k>=3時：

fork=3:N+n

Data(k,l)=U(k)+fl(l)*U(k-l)+fl(2)*U(k-2);

Data化,2)=Y(k)+fl(l)*Y(k-l)+fl(2)*Y(k-2);

end

end

dispC由廣義最小二乘法得到的辨識結(jié)果如下：')

Gyal=thita(l),Gya2=thita(2|,Gybl=thita(3),Gyb2=thita(4)

實驗結(jié)果:

實驗結(jié)果：

Uyl.txtuy2.txtuy3.txt

Gyal=Gyal=

Gyal=

-1.4630

“c,一1?4669

-1.4694

Gya2=Gya2=

Grya2=J

0.86470，85650.8607

Gybl=Gybl=Gybl=

0.42880.74150.7425

Gyb2=Gyb2=Gyb2=

0.55490.56880.5553

6.遞推廣義最小二乘法：

%首先將TXT中的數(shù)據(jù)裝載在Matlab中：

clear;clc;

Wbs=0;

whileWbs<l

disp('從下列文件中選擇所需加載文件：uyl.txtuy2.txtuy3.txt')

filename=input('輸入所需打開文件名：s');

[Wbs,Message]=fopen(filename,*r*);%打開文件:

end

Data=fscanf(Wbs,*%g%g',[2inf]);%生成數(shù)據(jù)矩陣

Data=Data';

ul=Data(:,1);

yl二Data(:,2);

n=2;

N=500-n;

n0=50;

先用前50組數(shù)據(jù)利用Ism估計出初值

%構(gòu)造矩陣phi

phi(:,l)=-yl(n:n+n0-l);

phi(:,2)=-y1(n-l:n+n0~2);

phi(:,3)=ul(n:n+nO-1);

phi(:,4)=ul(n-l:n+nO-2);

y(:,l)=yl(n+1:n+n0);

seta=(inv(phi，*phi))*phi**y;

setaO=seta;

f0=[00]';

plO=inv(phi'*phi);

p20=10-6*eye(2,2);

fori=nO:N-2

y(1:i,l)=yl(n+1:n+i);

u(l:i,l)=ul(n+l:n+i);

%計算殘差

e=y-phi*setaO;

%f估計

omega=[-e(n+i-2)-e(i+lT)]';

K2=p20*omega*inv(1^omega，*p20*omega);

fl=fO+K2*(e(n+i-2)-omega,*f0);

p2=p20-K2*omega**p20;

p20=p2;

fO=f1;

先濾波

yy(1:i,l)=yl(n:n+i-l);

yy(1:i,2)=yl(n-l:nH-2);

uu(l:i,l)=ul(n:(n+i-1));

uu(l:i,2)=ul((n~l):(n+i-2));

yg=y+yy*fl；%yg(k)=y(k)+f(l)*y(k-1)+f(2)*y(k-2)

ug=u+uu*fl;%ug(k)=u(k)+f(l)*u(k-l)+f(2)*u(k-2)

%seta估計

ksai=[-yg(n+i-2)-yg(n+i-3)ug(n+i-2)ug(n+i-3)],;

Kl=plO*ksai*inv(l+ksai,*plO*ksai);%K(N+1)值

pl=plO-Kl*ksai'*plD;%P(N+1)值

plO=pl;

setal=setaO+Kl*(yl(n+i+l)-ksai,*setaO);%seta(N+l)值

setaO=setal;

%構(gòu)造新phi陣

phi(1:i+1,l)=-yl(n:n+i-l+l);

phi(1:i+1,2)=-yl(n-1:n+i-2+l);

phi(1:i+1,3)=ul(n:n+i-l+l);

phi(1:i+1,4)=ul(n-1:n+i-2+l);

end

disp('遞推廣義最小二乘法辨識結(jié)果：')；

al=setaO(l),a2=seta0(2),bl=setaO(3),b2=seta0(4)

實驗結(jié)果：

Uyl.txtuy2.txtuy3.txt

al=a-=al=

-0.7976-0.7872-0.8052

a2=a2=a2=

-0.1941-0.2300-0.1857

bl=bl=bl=

-0.00730.15770.1552

b2=b2=b2=

0.36830.50900.4912

7.夏氏偏差修正法：

%首先將TXT中的數(shù)據(jù)裝載在Matlab中：

clear;clc;

Wbs=0;

whileWbs<l

disp(,從下列文件中選擇所需加載文件：uyl.txtuy2.txtuy3.txt')

filename=input('輸入所需打開文件名's');

[Wbs,Message]=fopen(filename,T);%打開文件;

end

Data=fscanf(Wbs；%g%g',[2inf]);%生成數(shù)據(jù)矩陣

Data=Data';

%使用夏氏修正法，對2階系統(tǒng)進行參數(shù)辨識

n=2;L=length(Data);N=L-n;

U=Data(:,l)；

Y=Data(:z2);

%首先使用最小二乘法，求系統(tǒng)參數(shù)向量初值

glOL=(-Y(2:L-l)/-Y(l:L-2),U(2:L-l),U(l:L-2)];

Zgll=Data(3:L,2)；

Sgll=glOL'*glOL;Sgl2=inv(Sgll);Sgl3=glOL,*Zgll;

Xsthita=Sgl2*Sgl3;%計算參數(shù)矩陣

%得到初俏后，以下使用夏氏修正法進行參數(shù)辨識；

thitabO=0;%設(shè)偏差項的偏差初值為0

Fa=Sgl2*glOL';

M=eye(N)-glOL*Sgl2*glOL;

F=eye(N);

if(F>=10A-6*eye(N))

El=Zgll-glOL*Xsthita;%計算殘差E

omiga(2:N,l)=-El(l:N-l);

omiga(3:N,2)=-El(l:N-2);

D=omiga'*M*omiga;

Fx=inv(D)*omiga'*M*Zgll;

thitab=Fa*omiga*Fx;

Xsthita=Xsthita-thitab;

F=thitab-thitabO;

thitabO=thitab;

end

dispC夏氏修正法,)

Xsal=Xsthita(l),Xsa2=Xsthita(2),Xsbl=Xsthita(3),Xsb2=Xsthita(4)

實驗結(jié)果:

Uyl.txtuy2.txtuy3.txt

Xsal=Xsal=Xsal=

-1.4789-1.4680-1.4651

Xsa2=Xsa2=Xsa2=

0.86590.85480.8603

Xsbl=Xsbl=Xsbl=

0.42590.73920.7423

Xsb2=Xsb2=Xsb2=

0.55530.57120.5562

8.夏氏改良法：

%首先將TXT中的數(shù)據(jù)裝載在Matlab中：

clear;clc;

Wbs=0;

whileWbs<l

disp('從下列文件中選擇所需加載文件：uyl.txtuy2.txtuy3.txt')

filename=input(輸入所需打開文件名N,'s');

[Wbs,Message]=fopen(filename,T);%打開文件；

end

Data=fscanf(Wbs；%g%g',[2inf]);%生成數(shù)據(jù)矩陣

Data=Data';

%最小二乘法，取b0=0

n=2;%根據(jù)題目,本系統(tǒng)為2階系統(tǒng)

L=length(Data);%辨識所需數(shù)據(jù)的總長度

N=L-n;%構(gòu)造測量矩陣的數(shù)據(jù)長度

FIA=zeros(N,2*n);%構(gòu)造測量矩陣

fori=1:N%測量矩陣賦值

forI=l:n*2

if(l>n)

FIA(iJ)=Data(n+2+i-l,l);

else

FIA(iJ)=-Data(n+i-l,2);

end

end

end

Y=Data(n+l:n+N,2);%輸出數(shù)據(jù)矩陣

thita=inv(FIA'*FIA)*FW*Y;%計算參數(shù)矩陣

thitaLS=thita;%最小二乘估值

m=inv(FIA,*FIA)*FIA,;

fori=1:1:1000%構(gòu)造Omega

Err=Y-FIA*thitaLS;%計算殘差E

omiga(2:N,l)=-Err(l:N-l);

omiga(3:N,2)=-Err(l:N-2);

F=inv(omiga'*omiga)xomiga'*Err;

thita=thitaLS-m*omiga*F;

end

dispC使用夏氏改良法辨識結(jié)果為：,)

Xgal=thita(l),Xga2=thita|2),Xgbl=thita(3),Xgb2=thita(4)

實驗結(jié)果：

Uyl.txtuy2.txtuy3.txt

Xgal=Xgal=Xgal=

-1.4779-1.4674-1.4651

Xga2=Xga2=Xga2=

0.86540.85440.8603

Xgbl=Xgbl=Xgbl=

0.42580.73900.7424

Xgb2=Xgb2=Xgb2=

0.55570.57160.5562

9.夏氏遞推法：

%首先將TXT中的數(shù)據(jù)裝載在Matlab中：

clear;clc;

Wbs=0;

whileWbs<l

dispC從下列文件中選擇所需加載文件：uyl.txtuy2.txtuy3.txt')

filename=input('輸入所需打開文件名：s');

[Wbs,Message]=fopen(filename,'r');%打開文件;

end

Data=fscanf(Wbs,*%g,[2inf]);%生成數(shù)據(jù)矩陣

Data=Data,;

%這里選用450組數(shù)據(jù)進行遞推迭代

n=2;N=450-n;

U=Data(:,1);

Y=Data(:,2);

%初始賦值如下：

beta=[000000]*;P0=10"10*eye(6,6);

El=0;E2=0;

%循環(huán)迭代

fori=3:N

El

Y(n+i)-beta(l)*Y(n+i-l)-beta(2)*Y(n+i-2)-beta(3)*Y(n+i-l)-beta(4)*U(n+i-2);

E2

Y(n+i-1)-beta(l)*Y(n+i-2)-beta(2)*Y(n+i-3)-beta(3)*U(n+i-2)-beta(4)*U(n+i-3);

FIA=[-Y(n+i)-Y(n+i-l)U(n+i)U(n+i-l)-El-E2]';

K=P0*FIA*inv(l+FIA，*P0*FIA);

Pl=P0-K*FIA'*P0;

P0=Pl;

betal=beta+K*(Y(n+i+1)-FIA**beta);

beta=betal;

end

dispC取450組數(shù)據(jù)進行遞推夏氏辨識的結(jié)果為:‘)；

Dxal=beta(l),Dxa2=beta(2),Dxbl=beta(3),Dxb2=beta(4)

實驗結(jié)果:

Uyl.txtuy2.txtuy3.txt

Dxal=Dxal=Dxal=

-1.4067-1.5289-1.4919

Dxa2=Dxa2=Dxa2=

0.80490.94510.8810

Dxbl=Dxbl=Dxbl=

0.43000.74930.7385

Dxb2=Dxb2=Dxb2=

0.56500.54320.5410

10.增廣矩陣法：

%首先將TXT中的數(shù)據(jù)裝載在Matlab中：

clear;clc;

Wbs=0;

whileWbs<l

dispC從下列文件中選擇所需加載文件：uyl.txtuy2.txtuy3.txt')

filename=input('輸入所需打開文件名：s');

[Wbs,Message]=fopen(filename,'r');%打開文件;

end

Data=fscanf(Wbs,*%g,[2inf]);%生成數(shù)據(jù)矩陣

Data=Data,;

%使用增廣矩陣法進行2階系統(tǒng)的參數(shù)辨識

n=2;L=length(Data);

U=Data(l:L,1)；Y=Data(l:L,2);

fori=1:3*n%首先給辨識參數(shù)向量賦初值，這里取thitaO為0

thitaO(i,1)=0;

end

P0=1010*eye(3*n,3*n)，然后設(shè)定初始狀態(tài)P0為充分大的對角陣

ERRO=0.00000000001;%收斂情況指標(biāo)，相對誤差ERRO

thita=[thitaO,zeros(3*n,L)];%設(shè)定參數(shù)向量的初值及大小

zs=zeros(L,1);為構(gòu)造初始噪聲矩陣

hn=[0,0,0,0,0,0]';

forj=n+1:1:L

zs(j-l)=Y(j-l)-hn，*thitaO;

hn=[-Y(j-1),-Y(j-2),U(j-1),U(j-2),zs(j-1),zs(j-2)]*;

K=PO*hn*inv(hn'*PO*hn+l)計算修正矩陣K

thitaO=thitaO+K*[Y(j)-hn，*thitaO];

P0=P0-K*hn'*P0;%計算狀態(tài)矩陣P(N)

thita(:,j)=thitaD;

E=abs((thita(:,j-1)-thitaO)./thitaO);

ifECERRObreak;%限據(jù)設(shè)定的誤差限，判斷何時跳出

end

end

Zgmal=thita0(l),Zgma2=thita0(2),Zgmbl=thita0(3),Zgmb2=thita0(4),

Zgmcl=thita0(5),Zgmc2=thita0(6)

Uyl.txtuy2.txtuy3.txt

Zgmal=Zgmal=Zgmal=

-1.4683-1.4629-1.4573

Zgma2=Zgma2=Zgma2=

0.86260.85200.8543

Zgmbl=Zgmbl=Zgmbl=

0.42460.73960.7427

Zgmb2=Zgmb2=Zgmb2=

0.55400.57510.5563

Zgmcl=Zgmcl=Zgmcl=

-0.92050.02820.5183

Zgmc2=Zgmc2=Zgmc2=

0.0436-0.26890.2806

11.極大似然法：

%首先將TXT中的數(shù)據(jù)裝載在Matlab中：

clear;clc;

Wbs=0;

whileWbs<l

disp。從下列文件中選擇所需加載文件：uyl.txtuy2.txtuy3.txt')

filename=input(輸入所需打開文件名N,'s');

[Wbs,Message]=fopen(filename,T);%打開文件；

end

Data=fscanf(Wbs；%g%g\[2inf]);%生成數(shù)據(jù)矩陣

Data=Data';

n=2;L=length(Data);N=L-n;

U=Data(:,l);Y=Data(:,2);

mthita=zeros(3*n,l);%給被辨識參數(shù)矩陣賦初始值

%用最小二乘法求初始的被辨識參數(shù)矩陣

OL=(-Y(2:L-l),-Y(l:L-2)/U(2:L-l)/U(l:L-2)];

Zl=Data(3:L,2);

thita=inv(OL'*OL)*OL'*Zl;

mthita(l:2*n,l)=thita;

%得到出之后，使用牛頓―拉布森極大似然法編程

F=7Pros(l,1);J==1;

Eal=zeros(L,l);Ea2=zeros(L,l)；

Ebi=zeros(L,l);Eb2=zeros(L,l);

Eel=zeros(L,l);Ec2=zeros*L,l)；

%給J的梯度和Hassian矩陣賦初始值

Ethita=zeros(3*n,l);

fd=zeros(3*n,l)；

Hassian=zeros(3*n,3*n);

forp=1:500

Jal=mthita(l);Ja2=mthita(2);Jbl=mthita⑶;Jb2=mthita(4);Jcl=mthita(5);Jc2=

mthita(6);

fori=n+l;L

y2(i)=-Jal*Y(i-l)-Ja2*Y(i-2)+Jbl*U(i-l)+Jb2*U(i-2)+Jcl*E(i-l)+Jc2*E(i-2);

E(i)=Y(i)-y2(i);

J=0.5*(E(i)A2);

sgma=(l/N)*(E(i)A2);

end

forj=n+l:l:L

Eal(j)=Y(j-l)-[JclJc2]*[Eal(j-l)/Eal(j-2)],;Ea2(j)=Y(j-2)-[Jcl/Jc2]*[Ea2(j-l),Ea2(j-2)]';

Ebl(j)=?U(j-lHJcl,Jc2]*[Ebl(j-l),Ebl(j?2)了;Eb2(j)=

-U(j-2)-[JclJc2]*[Eb2(j-l),Eb2(j-2)]';

Ecl(j)=-E(j-l)-[Jcl,Jc21*[Ecl(j-l),Ecl(j-2)],;Ec2(j)=-E(j-2)-[Jcl,Jc2]*[Ec2(j-l),Ec2(j-2)]';

Ethita=lEal(j),Ea2(j),Ebl(j),Eb2(j),Ecl(j),Ec2(j)];

fd=fd+E(j)*Ethita,;%迭代計算J的梯度矩陣

Hassian=Hassian+Ethita'*Ethita;%迭代計算Hassian矩陣

end

%用牛頓-拉卜森法計算被辨識的參數(shù)的估值

mthita=mthita-inv(Hassian)*fd;

JFIA=(sgma-sgmaO)/sgmaO;

sgmaO=sgma;

%跳出迭代的判斷

ifJFIA<0.000000001break

end

end

dispC極大似然法(牛頓-拉卜森方法)，得到的結(jié)果為：1)

Jal=mthita(l)Ja2=mthita(2)Jbl=mthita⑶,Jb2=mthita(4)Jcl=mthita(5)Jc2=mthita(6i

實驗結(jié)果:

Uyl.txtuy2.txtuy3.txt

Jal=Jal=Jal=

-1.4789-1.4680-1.4651

Ja2=Ja2=Ja2=

0.86590.85480.8603

Jbl=Jbl=Jbl=

0.42590.73920.7423

Jb2=Jb2=Jb2=

0.55530.57120.5562

Jcl=Jcl=Jcl=

-0.9669-0.04110.5206

Jc2=Jc2=Jc2=

-0.4682-0.3384-0.0694

12.分析噪聲式幻特性：

RLS法：當(dāng)噪聲比較小時，辨識精度較高；當(dāng)噪聲比較大時，收斂點可能

不唯一，參數(shù)估計值往往是有偏的；但是運用此方法時，數(shù)據(jù)要充分多，否則辨識精度明顯

下降；另外噪聲模型階次不宜過高，也會影響精度，所以誤差的原因主要是數(shù)據(jù)量的大小以

及噪聲模型階次。

IV法：辨識效果較好，能適應(yīng)較大的范圍的噪聲特性。對初值P(0)敏感，

可用LS法起步，否則沒有可靠的收斂性；對數(shù)據(jù)u(k;在z(k)的直流分量敏感，對

z(k)的直流分量不敏感。所以誤差的主要原因是初值P(0)以及數(shù)據(jù)計算過程產(chǎn)生的誤差。

ELS法：ELS是極大似然法的一種近似形式，當(dāng)數(shù)據(jù)長度較小時，辨識精

度可能優(yōu)于極大似然法，但數(shù)據(jù)長度較大時，精度低于極大似然法。

第三題

以上題的結(jié)果為例，進行：(10分)

1.分析比較各種方法估”的精度；

2.分析其計算量；

3.分析噪聲方差的影響；

4.比較白噪聲和有色噪聲對辨識的影響。

答：

1..分析比較各種方法估L的精度；

最小二乘估計雖然不能滿足量測方程中的每一個方程，使每個方程都有偏差，但它使所

有方程偏差的平方和達到最小，兼顧了所有方程的近似程度，使整體誤差達到最小。軾助變

量法、最小二乘法、極大似然法對輸入輸出模型有相同的精度.

增廣最小二乘法所得結(jié)果中b2誤差較大；

輔助變量法辨識精度高于最小二乘法

夏氏法辨識精度較高；克眼有偏估計。

2.分析其計算量；

最小二乘估計計算量大、存儲大；

增廣最小二乘遞推算法，擴充了最小二乘法的參數(shù)向量和數(shù)據(jù)向量的維數(shù),把噪聲模型

的辨識同時考慮進去，所以計算最較大；

遞推廣義LS中的計算分為二部分：

第一部分：按遞推LS法，隨N增大不斷計算參數(shù)所以計算量較大，

第二部分：在遞推過程中，由于是變化的故而需要不斷計算；

輔助變量法：(1)計算與基本最小二乘估計同樣簡單;

⑵辨識精度高于最小LS估計法；

⑶是一種無偏估計方法；

(4)參數(shù)估計時需構(gòu)造輔助變量矩陣。

夏氏法：(1)無偏的估計方法；

(2)計算量較廣義LS的要小許多；

(3)不需進行I/O數(shù)據(jù)的反復(fù)過濾，計算效率高：

(4)其遞推算法可推廣至MIM0系統(tǒng)，而廣義LS法則不行；

(5)是一種循環(huán)迭代方法，收斂速度較LS要慢；

(6)估計結(jié)果較好，可分為夏式修正法和夏式改良法兩種

極大似然;計算量較小。

3.分析噪聲方差的影響；

噪聲越大，對最小二乘法辨識效果影響較大，對于增廣遞推最小二乘法，抗噪聲能力較

強。

通江編程計算，獲得在噪聲方差比較小的情況下，各種方法所獲得的估值比較理想，但隨著

噪聲方差的增大，估值的偏差隨之增大，橫向比較看來復(fù)式法與廣義最小二乘法能夠

更好地還原參數(shù)值，當(dāng)觀測值足夠多時，各種方法都能很好地反映參數(shù)真實值。

4.比較白噪聲和有色噪聲對辨識的影響“

由于所用的輸出觀測值有有色噪聲成分和白噪聲成分，所以辨識結(jié)果有誤差，且用最小

二乘法計算得到的參數(shù)估計誤差較大。但是有色噪聲誤差明顯大于白噪聲。

第四題.

系統(tǒng)模型階次的辨識：(ID分)

用三種方法確定系統(tǒng)的階次并辨識；

比較所用三種方法的優(yōu)劣及有效性；

答：

1.按殘差方差定階法(含估計值誤差方差最小法和F檢驗法)和殘差白色定階法:

%首先將TXT中的數(shù)據(jù)裝載在Matlab中：

clear;clc;

Wbs=0;

whileWbs<1

disp('從下列文件中選擇所需加載文件:uyl.txtuy2.txtuy3.txt,)

filename=input('輸入所需打開文件名：s');

[Wbs,Message]=fopen(fi1ename,'r');%打開文件;

end

Data=fscanf(Wbs/%g,[2inf]);%生成數(shù)據(jù)矩陣

Data=Data>;

%以估計誤差方差最小為標(biāo)準(zhǔn)的按殘差定階

L=length(Data);%輸入輸出數(shù)據(jù)長度

Jn=zeros(1,10);

t=zeros(1,10)：

rm=zeros(10,10);

%由階次n=l開始計算

forn=l:1:10;

N=L-n;

FIA=zeros(N,2*n)G構(gòu)造測量矩陣

du=zeros(n,1);

dy=zeros(n+1,1);

rl-0;r0=0;

先使用LS法進行參數(shù)辨識

fori=1:N嬌則量矩陣賦值

for1=1:n*2

if(K=n)

FIA(i,1)=-Data(n+i-l,2);

elseif(n+i-l+2>0)

FIA(i,1)=DataCn+i-1+2,1);

end

end

end

Y=Data(n+l:n+N,2);與輸出數(shù)據(jù)矩陣

t