計(jì)算機(jī)課件系統(tǒng)辨識(shí)與建模2

系統(tǒng)辨識(shí)與建模

主講：黃東

第三講參數(shù)估計(jì)的批量法

■最小二乘算法

■參數(shù)估計(jì)的一般性質(zhì)

■最小二乘、加權(quán)最小二乘估計(jì)性質(zhì)

■噪聲方差估計(jì)

■廣義最小二乘

■偏倚校正算法

■輔助變量法

■多步最小二乘

■相關(guān)最小二乘

匯總的觀

察誤差

《最小二乘算法

考慮差分方程：AQ-)y(k)=BG)u(k)+w<k),其中

w(k)為白噪聲。假定模型的結(jié)構(gòu)已知(n,m,T),

將其寫成線性回歸模型：

y(k)=-a1y(k-l)-...-any(k-n)+b1u(k-l)+

...+bmu(k-m)+w(k)

B—[外,???，七力]??.，bmF

cp(k)=[-y(k-l),...,-y(k-n),u(k-l),...,u(k-m)]T

y(k)=(pT(k)0+w(k)

>--------

y(k)=-aiy(k-l)-…-any(k-n)+biU(k-l)+...+bmu(k-m)+w(k)

y(k+l)=-a1y(k)-...-any(k-n+l)+b1u(k)+...+bmu(k-m+l)+w(k+l)

y(k+2)=-a1y(k+l)-...-any(k-n+2)+b1u(k+l)+...+bmu(k-m+2)+w(k+2)

(pT(k)=[-y(k-l)…-y(k-n)u(k-l)...u(k-m)]

(pT(k+l)=[-y(k)...-y(k-n+l)u(k)...u(k-m+l)]

(pT(k+2)=[-y(k+l)...-y(k-n+2)u(k+l)...u(k-m+2)]

(pT(k+3)=[-y(k+2)…-y(k-n+3)u(k+2)...u(k-m+3)]

4最小二乘算法

若我們的觀測(cè)數(shù)據(jù)可寫出N個(gè)這樣的等式,

yN=oe+wN，

OT=[cp(k),……,(p(N+k-l)]

①TYN=(pT00_1.(pTWN

T1T

Q=((D(D)-(DYN-(①T(P尸①TWN

T

eLS=(<P0)iCiJTYN

T

條件1：E(0WN)=0(w(k)白色)

條件2:eT①可逆

①TWN=

-y(k-l)-y(k)-y(k+l)…-y(k+N-2)w(k)

-y(k-2)-y(k-l)-y(k)...-y(k+N-l)w(k+l)

-y(k-n)-y(k-n+l)-y(k-n+2)…-y(k-n+N-l)

u(k-l)u(k)u(k+l)...u(k+N-2)

u(k-m)u(k-m+l)u(k-m+2)...u(k-m+N-l)w(k+N-l)

y(k)與w(k)、w(k-l)、w(k-2)、…相關(guān)。當(dāng)w(k+l)與w(k)、w(k-l)、

w(k-2)、…不相關(guān)時(shí)，y(k)與w(k+l)也不相關(guān)。

4最小二乘算法

以上結(jié)果等同于求使：

J=Z(y(k)-cpT(k)0)2

最小的e,因此稱為最小二乘算法。

-=0

de

T

①丁YN=0cpe--正則方程

①T①-----------正貝U矩陣(對(duì)稱、正定、可逆)

《加權(quán)最小二乘

若對(duì)各次觀測(cè)數(shù)據(jù)加不同的權(quán)，即求使

J=Za(y(k)-(pT(k)0)2最小的0,則得到參數(shù)的

加a權(quán)最k小二乘估并：

0LSa=(SA(l))T<DTAYN

A的對(duì)角元由cik構(gòu)成

參數(shù)估計(jì)的一般性質(zhì)

無(wú)偏性

如果E(3)=0,(6=G"或E(d)=6,則稱估計(jì)為

無(wú)偏的。

2.有效性

如果無(wú)偏估計(jì)在滿足Cov0)=M-i,則稱估計(jì)為有效

的。其中：”日理―自理稱為Fisher

信息矩陣，其逆M-i稱為Crammer-Rao下界。在

一般情況下，有Crammer-Rao不等式:

Cov(3)=E{(d—。)3—e)丁}》MT

4有效性例

對(duì)于z=HB+w,若噪聲w為零均值、協(xié)方差陣

±w=E(WWT)的正態(tài)分布噪聲，即：W?N(ORw)，

則輸出信號(hào)Z?N(E(H0)Rw)，即，

2-11T-1

p(z|8)=(2?")exp(--[Z-E(H3)][Z-

010gp(z?e)_i

----------------------={￡(/T)￡：[z-￡(幺咽}

因此:"于是,T

M=E(HZW

明),其Crammer-Rao不等式為：

T

Cove)>E(HZwiH)-i0

有效估計(jì)也稱為最小方差估計(jì)、馬爾可夫估計(jì)。

&參數(shù)估計(jì)的一般性質(zhì)_

3.一致性

若估計(jì)?為漸近無(wú)偏的代叫E(3)=0),且

L巴Var)=0,則稱0為一致估計(jì)。

Var(<9)=E{-e)T--e)}

最小二乘、加權(quán)最小二乘估計(jì)性質(zhì)

估計(jì)性質(zhì)最小二乘

無(wú)偏6=(①T①尸①TWN，當(dāng)：1W(k)為零均值獨(dú)立隨機(jī)過(guò)

性

程(白噪聲)時(shí)；或者2E(①丁班)力時(shí)，有：

E(^)=0o無(wú)偏條件較嚴(yán)格。

有效CovC)=E{(①T①)1CpTW&N①(①T①)T}

T生

二(①T①尸①TZw①(①T①)T訓(xùn)工若W是服從正態(tài)分

布的獨(dú)立隨機(jī)過(guò)程，則Zw二。：1，

因而Cov(3)=囁1(4)3)-』一1。

若W是服從正態(tài)分布的獨(dú)立隨機(jī)過(guò)程，則

性27r

L四Cov⑴二例"芍尸二o,其中收斂于一個(gè)有

界非奇異矩陣。而Var(3)只是Cov")的對(duì)角線上元

的和，因而也有：“Var(3)=0o

最小二乘、加權(quán)最小二乘估計(jì)

性質(zhì)

估計(jì)性質(zhì)加權(quán)最小二乘

無(wú)偏0=(①TA①尸①TAWN

性

無(wú)偏條件同最小二乘估計(jì)。

有效

Cov(3)=(0TA0)-10TAzA①(①TA①尸,當(dāng)

T生W

A二Zw"時(shí)，Cov(3)=(eTgw「l①)T=獷，此時(shí)不

要求W是服從正態(tài)分布的獨(dú)立隨機(jī)過(guò)程

一致若W是服從正態(tài)分布的獨(dú)立隨機(jī)過(guò)程，則

性

的Cov(^)=^7(7)-1=0,其中收斂于一個(gè)有

界非奇異矩陣。而VarG)只是Cov(3)的對(duì)角線上元

的和，因而也有：々"Var(^)=0o

影響最小二乘計(jì)算結(jié)果的因素:

■u(k)QTQ是否可逆

■噪聲w(k)02大，則的方差大

W0

白色零均值，0是無(wú)偏估計(jì)

■數(shù)據(jù)總量NN越大，J的方差越小

fork=l:inl%每一行中的變量循環(huán)

Ls.mfori=l:lll%列循環(huán)

function[zta,m,tao]=ls(tt)forj=l:m(k)%每行變量中的觀測(cè)

數(shù)據(jù)循￥

%最小二乘法forMISO,tt的格式為:

第一列是系統(tǒng)輸出數(shù)據(jù)，其它列jtao=j+tao(k);%構(gòu)造時(shí)考慮時(shí)

是對(duì)應(yīng)的輸入數(shù)據(jù)延

ll=size(tt);%得到數(shù)據(jù)維數(shù)ifk>lff(i,j+kn)=tt(i+n-

日1(2)-1;jtao+l,k)；end

m(l)=inpu4輸入A(z)的階次》ifk==l,ff(i,j)=-tt(i+n-jtao,k);end

%指走模型結(jié)構(gòu)%輸出堤墓變號(hào)

tao(l)=0;end;end;

fori=l:rkn=kn+m(k);%算出行變量的啟

始位置

%給出輸入編號(hào)

m(i+l)=input(‘輸入B(z)的階次)+1;end;

%躋次加一，表示蓊疑個(gè)數(shù)fori=l:lll

tao(i+l)=input(‘輸入B⑵的時(shí)延')；yy(i)=tt(i+n,l);%構(gòu)造輸出向量

endend;

n=m(l)+max(tao);%算出一個(gè)fa=fP*ff;%最小二乘計(jì)算

方座最多使用的數(shù)據(jù)

fay=ff'*yy';

lll=ll(l)-n;%算出可列出的方程數(shù)zz=mv(ta);

inl=r+l;%構(gòu)造觀測(cè)數(shù)據(jù)矩陣什zta=zz*fay%顯示參數(shù)和結(jié)構(gòu)

kn=0;m,tao

Lsl.m

clearall;%最小二乘法forMISO

loady3;

tt(:,l)=uyr(l,：y；%讀入數(shù)據(jù)，并賦給變量tt。

tt(:,2)=uyr(2,:)1;%tt的格式為：第一列是系統(tǒng)

輸出數(shù)據(jù)，其它列是對(duì)應(yīng)的輸入數(shù)據(jù)

clearuyr;

ls(tt);

《仿真例

1.無(wú)噪聲模型：數(shù)據(jù)文件y3.mat

(l+1.5q1+0,7q-2)y(k)=q23,2u(k)

辨識(shí)結(jié)果(給定結(jié)構(gòu)：m=21,tao=02)

zta=

1.5000

0.7000

3.2000

2.白噪聲模型：數(shù)據(jù)文件y5-mat

(l+1.5q1+0.7q2)y(k)=q23.2u(k)+w(k)

辨識(shí)結(jié)果(給定結(jié)構(gòu)：m=21,tao=02)

zta=

1.5027

0.7032

3.1935

3.有色噪聲模型：數(shù)據(jù)文件y6.mat032

"0.32w（n左H）

（l+LSqT+O7qOMIQnq七.ZMIO+ITrTT^iT7

辨識(shí)結(jié)果（給定結(jié)構(gòu)：m=21，tao=02）

zta=

0.9757

0.1782

3.3955

4.白噪聲模型b：數(shù)據(jù)文件y5b.mat

3.2/〃（左）

y（k）i5「l+w（k）

辨識(shí)結(jié)果(給定結(jié)構(gòu)：m=21tao=02)

zta=

1.1627

0.3750

3.3907

追跡:對(duì)角線

元素之和

三噪聲方差估計(jì)>/—

e=y；-y=(i>o+wN-<i>

T1T

=WN-0(00)0WN//

T1T2

=(I-0(00)0)WN=BWN/41為BT=B,B=B

所以，E(eTe)=E(WjBWy

=E(WNTWN*NT0(0T0)-I0TWJ

=b：N?b：(Tr6(<pT?尸(pT)=^N-dim0)

w：=E(y(k)?\儲(chǔ)2/(N?dime)

在有色噪聲環(huán)境下，最小二乘估計(jì)是有偏的。下面的一些

算法對(duì)最小二乘進(jìn)行改進(jìn)。

廣乂最小―^乘

11i

考慮差分方程：AG)y(k)=BG)u(k)+-77w(k),

其中Mk)為白噪聲。假定模型的結(jié)構(gòu)已知

(n,m)o

如果噪聲模型C(小)已知，顯然用C(「)對(duì)輸入/輸出

數(shù)據(jù)進(jìn)行濾波，則可得到滿足最小二乘估計(jì)無(wú)偏

條件的模型：AG)Rk)：B(qy?+w(k),其中:

(k)久()y(蛇，(k)u=C()u(此。

在C(小)未知時(shí)，我們可考慮采用迭代估計(jì)的方法

去求得。

4廣義最小二乘的計(jì)算步驟

1令品(9)=1,i=0下標(biāo)表示迭代次數(shù)；。令000000

2計(jì)算，k)=C()y(k),G(k)=C6)u(k)；i=i+1;

3用最小二乘估計(jì)A(J)y(k)=BG13(k)+w(k)中的參數(shù)；

4用估計(jì)模型河q)、燈廠I)以及各時(shí)刻的觀測(cè)數(shù)據(jù)，計(jì)算

出殘差，左)：/(k)=2G)y(k)-JG)u(k)

5計(jì)算圻Ef(k)及如果&小于一定數(shù)，則結(jié)束

辨識(shí)。否則底下一步。

6對(duì)于噪聲模型C(44k)=w(k),用最小二乘估計(jì)出參數(shù),

得到更新的C/)后丁返回2。

以上算法的每一次循環(huán)中都要進(jìn)行濾波和兩次求逆。下面

的算法將在計(jì)算工作量上有所改進(jìn)。

LsO.mifk>l,ff(i,j+kn)=tt(i+n-

jtao+lk);end

functionz

[ztam,tao]=lsO(tt,mtaoll)ifk==l,ff(i,j)=-tt(i+n-

zzzjtao,松;end%輸出變

%最小二乘法forMISO,tt的格式為：第量變導(dǎo)

一列是系統(tǒng)輸出數(shù)據(jù)，其它列是對(duì)應(yīng)

的輸入數(shù)據(jù)end;end;

%m為各多項(xiàng)式中參數(shù)個(gè)數(shù)，應(yīng)與tt的列kn=kn+m(k);%算出行

數(shù)一致；tao為時(shí)逅ll=size(tt);變量的官外位置

n=max(m)+max(tao);%算出

一個(gè)方盒最多使用的數(shù)據(jù)end;fori=l:lll

lll=ll(l)-n;%算出可列出的方程數(shù)yy(i)=tt(i+n,l)；%構(gòu)造輸出

向量

inl=ll(2);%構(gòu)造觀測(cè)數(shù)據(jù)矩陣ff

end;

kn=0;

fa=ff'*ff;%最小二乘計(jì)算

fork=l:inl%每一行中的變量

循環(huán)fay=ff'*yy';

fori=l:lll%列循環(huán)zz=inv(fa);

for1=1:m(k)%每行變量中zta=zz*fay%顯示參數(shù)和

曲觀測(cè)數(shù)徭循環(huán)結(jié)構(gòu)

jtao=j+tao(k);%構(gòu)造時(shí)考慮時(shí)延Mtao

fori=l:r

Gls.mi%給出輸入編號(hào)

m(i+l)=input(‘輸入B(z)的階次')+1;

%階次加一，表示參數(shù)個(gè)數(shù)

%廣義最小二乘法

clearall;for輸入的時(shí)延

MISOtao(i+l)=inputCB(z)')；

loady6;End

%讀入數(shù)據(jù)賦taomax=max(tao+m);

tt(:zl)=uyr(l,:)';

給tt。mc=input('輸入噪聲模型C(z)的階次');

tt(-2)=uyr(2,:)';%tt的

稱式為：第一列是系統(tǒng)輸出c=[l];d=[l];%初始化濾波器

數(shù)據(jù)，其它列是對(duì)應(yīng)的輸入lb=l;xci=100000;xci0=1000000;

數(shù)據(jù)whileabs(xci0-xci)>0.001%濾波

clearuy;plot(tt(:,l))循環(huán)

ll=size(tt);%得到數(shù)據(jù)維數(shù)xciO=xci;

inl=ll(2);r=inl-l;ttl=filter(c,d,tt);%輸入輸出濾波.

m(l)=input。輸入A⑵的階次');%主

%指虎模型結(jié)構(gòu)[zta,m,tao]=lsO(ttl,mztao,ll);

最小二乘

tao(l)=0;

fork=l:taomax%計(jì)算輸出估計(jì)

y(k)=tt(k,l);%設(shè)定初始輸出%計(jì)算方程誤差

endee=0;%計(jì)算損失函數(shù)值

fork=l+taomax:ll(l)fork=l:ll(l)

mm=0;ee=ee+e(k)*e(k);

f,i=l:inl

ifi==lend

forj=l:m(i)xci=ee/ll(l)

taoc=0;%估計(jì)噪聲模型

end;lc=size(e);

else[cc,mc,taoc]=lsO(e,mc,taoc,lc);

forj=l:m(i)c(2:mc+l)=cc總顯示參數(shù)和

ifk-tao(i)-j+l>0

結(jié)構(gòu)

fb(mm+j)=tt(k-tao(i)-j+l,i);

elselb=lb+l;

fb(mm+j)=0;end%結(jié)束濾波循環(huán)

end____________plot(e),ylabel(cerror,),

end;end;xlabel('time')

mm=mm+m(i);

end;

y(k)=fb*zta;

end

廣乂最小—^乘例

有色噪聲模型：數(shù)據(jù)文件032G

Ly"6.mat0.32w(左)

(1+L5qT+0.7q-2)y(k)=q-23.2u(k)+rr^^I7

辨識(shí)結(jié)果(給定結(jié)構(gòu)：m=21,tao=02)

zta=

1.4950

0.6943

3.2075

c=[1.0000-1.71180.8069]

迭代次數(shù)：6

2.白噪聲模型b：數(shù)據(jù)文件y5b.mat

（女）_

y（k）=―3.2+w（k）

1+1.5「+0.7/

辨識(shí)結(jié)果（給定結(jié)構(gòu)：m=21,tao=02）

zta=

1.5016

0.6992

3.1881

c=[1.0000-1.50671.4236-1.02150.5519-0.1689]

迭代次數(shù)：7

?W=Ay-Bu

T偏倚校正算法

仍疆差分方程：AG，洪疥才)u(k)t-(k),其中

w(k)為白噪聲。令彳而=Urw(k),財(cái)

AG)y(k)=BG)u(k)+^(k)o分別寫成回歸模型：

Y=3+卅CC+W,組合起來(lái)有Y=[<D,C]+W/

「八i-iL?！?/p>

e「①’①.①1「①’丫1

其最小二乘解為：7=利用分塊矩

陣求逆公式有：

^=(0T0)-10TY-(0T0)-10Toc

C=DTQTMYM=I-0(0T0)10T。D=JMQ

須要注意，c的求取仍然是一個(gè)迭代過(guò)程。

、偏倚校正算法的計(jì)算步驟

1）=1,用最小二乘法求/=%=（①T①尸①,丫，并保

留①、「=（①丁①尸①T以及M=I-①'

2計(jì)算j=丫-①?,并依據(jù)尸Co+W構(gòu)造。，計(jì)算

D=QTMQO

-1T

3計(jì)算d=DQMY,并計(jì)算A6?=「及^=0LS-N。

4若參數(shù)已收斂，則結(jié)束辨識(shí)，否則轉(zhuǎn)2。

以上算法的一次循環(huán)中沒(méi)有濾波，且只有一次求逆。如果

將第3步中°的計(jì)算改為：。=（。％）-％上，則還可省去D

的計(jì)算。（這一改進(jìn)由夏天長(zhǎng)首先給出。）本法可能會(huì)

出現(xiàn)收斂慢的情況，可用對(duì)獷求均值來(lái)解決

mc=input('輸入噪聲模型C(z)的階次')；

n=m(l)+max(tao);%胃由一個(gè)方連

Gls2.m最多使用的數(shù)據(jù)

lll=ll(l)-n;%算出可列出的方程數(shù)

clearall;%偏倚校正法forMISOlb=l;xci=100000;xci0=1000000;c=[l];

loady6;inl=r+l;%構(gòu)造觀測(cè)數(shù)據(jù)矩陣什

%讀入數(shù)據(jù),kn=0;

fork=l:inl%每一行中的變量循環(huán)

tt(:,2)=uyr(2,:)'；％tt的格式為

第一列桌系統(tǒng)輸出數(shù)據(jù)，其它列是:fori=l:lll%列循環(huán)

fori=l:m(k)%每行變量中的觀測(cè)數(shù)

據(jù)循環(huán)

clearuyr;plot(tt(:,l))

%構(gòu)造時(shí)考慮時(shí)延

ll=size(tt);%得到數(shù)據(jù)維數(shù)jtao=j+tao(k);

r=H(2)-l;ifk>lff(i,j+kn)=tt(i+n-jtao+l,k);end

ifk==l,ff(i,j)=-tt(i+n-jtao,k);end

m(l)=inpu4輸入A(z)的階次I%輸出差及變號(hào)

%指走模型結(jié)構(gòu)

tao(l)=0;end;end;

fori=l:rkn=kn+m(k);%算出行變量的啟始位

置

i%給出輸入編號(hào)

end;

m(i+l)=inputC輸入B(z)的階次)+1;

%躋次加一，表示蓑藏個(gè)數(shù)fori=l:lll

%構(gòu)造輸出向量

tao(i+l)=input(‘輸入B(z)的時(shí)延')；yy(i)=tt(i+n,l);

endend:

%最小二乘計(jì)算

end,end

fay=ff*yy';

cy=e(mc+l:lll+mc)';

zz=inv(fa);

fd=fc'*fc;%最小二乘計(jì)算

zaa=zz*ff';

fdcy=fc'*cy;

zta=zz*fay%顯示參數(shù)和結(jié)構(gòu)

zzc=inv(fd);

M,tao

c(2:mc+1)=[zzc*fdcy]'%

ztab=zta-zta;顯示參顛和結(jié)構(gòu)_____________

ztaa=zta;ztab=ztab+zaa*fc*c(2:mc+l)';

whileabs(xci0-xci)>0,05zta=ztaa-ztab/lb

xciO=xci;%zta=ztaa-zaa*fc*c(2:mc+1)'

e(l:mc)=O;%計(jì)算輸出誤差lb=lb+l;

e(mc+l:lll+mc)=yy'-ff*zta;end%結(jié)束濾波循環(huán)

%計(jì)算損失函數(shù)值

ee=O;

fork=l:lll+mc

ee=ee+e(k)*e(k);

end

xci=ee/(lll+mc)______________

%估計(jì)噪聲模型

fori=l:mc

forj=l:lll

《偏倚校正算法例

工.有色噪聲模型：數(shù)據(jù)文件y6.mat032G

"0.32)

(1+L5qY+0.7q-2)y(k)=q23.2u(k)+TTT77rm7

辨識(shí)結(jié)果(給定結(jié)構(gòu)：m=21,tao=02)

zta=

1.4879

0.6831

3.1976

c=[1.0000-1.70300.7983]

迭代次數(shù)：20

輔助變量法

分二乘法中，在Y=o)e+w的各項(xiàng)乘上①丁,然后利用①Tw的期

望值為零得到參數(shù)的無(wú)偏估計(jì)。受此啟發(fā)，若在Y=o)e+的合項(xiàng)

乘上中丁，使其滿足以下兩個(gè)條件：1UT的喇望值為零；2.w

可逆，則也可得到參數(shù)的無(wú)偏估計(jì)。下面討論輔助變量的選?。?/p>

設(shè)模型為A()y(k)=B()u(k)+(k),若u(k)與(k)不相關(guān)：

一1-1

a選取輔助模型D()2(k)=F(《)u(k),用z(k)3u(k)構(gòu)造肌

-1-1

b若系統(tǒng)的純時(shí)延r已知，則可用u(k-T)、u(k)構(gòu)造中；

c用u(k)、u(k)構(gòu)造中

d(k)=D()w(k),^D()的階次n已知,則可用y(k-n)、u(k)構(gòu)造中；

先求出最小二乘解e；s二(①'①)T①，然后依據(jù)()z(k)二

()11也)計(jì)算出輸出估計(jì)28),再用Z(k)、u(k)構(gòu)造史；

--1

n=m(l)+max(tao);%算出一個(gè)方

Ivjs程盤多使用的效據(jù)

lll=H(l)-n;%算出可列出的方程數(shù)

clearall;%輔助變量最小二乘法forinl=r+l;%構(gòu)造觀測(cè)數(shù)據(jù)矩陣"

MISOkn=0;

loady6;fork=l:inl%每一行中的變量循環(huán)

%讀入數(shù)據(jù),

m：牖盤黯fori=l:lll%列循環(huán)

fori=l:m(k)%每行變量中的觀測(cè)數(shù)

我■匕費(fèi)滁%輸出數(shù)據(jù)%,tt其的它格列式是為:據(jù)循環(huán)

對(duì)應(yīng)的輸入數(shù)據(jù)jtao=j+tao(k);%構(gòu)造時(shí)考慮時(shí)延

clearuyr;ifk>lff(ij+kn)=tt(i+n-

ll=size(tt);%得到數(shù)據(jù)維數(shù)jtao+l,k);end

ifk==lff(i,j)=-tt(i+n-jtaok);end

r=H(2)-l;%輸出z亞量變號(hào)z

m(l)=inputC輸入A(z)的階次')；

%指走模型結(jié)構(gòu)end;end;

kn=kn+m(k);%算出行變量的啟

tao(l)=0;始位置

fori=l:rend;

i%給出輸入編號(hào)fori=l:lll

m(i+l)=inputC輸入B(z)的階次)+1;

%躋次加一，表示塞疑個(gè)數(shù)yy(i)=tt(i+n,l);%構(gòu)造輸出向量

tao(i+l)=inputC輸入B(z)的時(shí)延')；end;

end

X

%f-1

k-=21,fh(i7tt(i+n-jtao-

匕

%最小二乘計(jì)算o+le

變,2

量I%以11%430)

為￥

fay=ff*yy';

zz=inv(fa);

ifk==l,fh(izj)=-tf(i+n-j);end

zta=zz*fay%顯示參數(shù)和結(jié)構(gòu)%以歲助假型輸出為*甫助變量

end;end;

%建立輔助變量

kn=kn+m(k);%算出行變量

%a=[lzta(l:m(l))'];的啟始7立置

%b=|~0zta(m⑴；⑴

end;

a=[l1.70.72];fori=l:lll

b=[001]____________________yy(i)=tt(i+n,l);%構(gòu)造輸出向量

tf=，ilter(6,a,tt(:,2));

end;

kn=u；fa=fh'*ff:%輔助變量最小

fork=l:inl%每一行中的變量循環(huán)二乘計(jì)算

fori=l:lll%列循環(huán)fay=fh'*yy';

fori=l:m(k)%每行變量中的觀測(cè)數(shù)

據(jù)循環(huán)zz=inv(fa);

zta=zz*fay%顯示參數(shù)和結(jié)

jtao=j+tao(k);%構(gòu)造時(shí)考慮時(shí)延構(gòu)

ifk>lfh(ij+kn)=tt(i+n-jtao+l,k);endM,tao

%ifk==l,fh(i,j)=-tt(i+n-

jtao+l,2);end%以u(píng)(k)為輔助變

量

《輔助變量法例

有色噪聲模型：數(shù)據(jù)文件032

y"6.mat0.32wd(H左)

(1+L5qY+0.7q-2)y(k)=q-23.2u(k)+TrT77rm7

辨識(shí)結(jié)果(給定結(jié)構(gòu)：m=21,tao=02)

L使用輔助模型：a=[11.70.72];b=[001];

zta=

1.4956

0.6962

3.2234

----------------------------

2.以ds為輔助模型：b=[003.3955];a=[l

0.97570.1782];

zta=[1.49080.68953.2235]

3.以u(píng)(k)為輔助變量

zta=[2.7082-1.85113.9692]

4.以u(píng)(k-tao)為輔助變量(方程病態(tài)，溢出)

5.以y(k?tao)為輔助變量

zta=[4.30792.91313.8725]

2-白噪聲模型b：數(shù)據(jù)文件y5b.mat

-2

八、3.2qu(k)

y(k)=--------------+w(k)

1+1.5q+U.7q

辨識(shí)結(jié)果(給定結(jié)構(gòu)：m=21,tao=02)

1.使用輔助模型：a=[l1,70.72];b=[001];

zta=[1.47400.68453.3692]

2.以?LS為輔助模型：b=[003.3907];

a=[l1.16270.3750];

zta=[1.44560.63893.3777]

多步最小一^乘

考晨分方程：A()y(k)=B()u(k)+”w(k),其中w(k)為

白噪聲。模型可改寫為C()A()y(k『C()B()u(k)+w(k)或

D()y(k)=F()u(k)+w(k),其中:

D()=C()A(),F()=C()B()---------**。

此模型可用最小二乘解出D()、F()o這是第一步。

第二步可有兩種不同方法：

a解同次第方程組由**式，可得關(guān)系：

AU)F(「)=B(「)D(「),兩邊分別展開,并按J的同

次嘉相等規(guī)則，可列出\+找+加個(gè)方程：F=H0,

=T

其中S[ap...ana,bp...bnb],F=[fp...fnb+nc,0,??.0]

000................0100....................00

*00............0-d110...................00

H=-f2-ft0..............0-d2-d11....................00

...............................................j.............................................-41

n+n+n

.......................~f2t...............................................~d2."dlabc行

...........至與...................

-f

1nb+nc............................「…................

n"-f,nb+nc",,"""■""""""""9-rJlna+nc.......................

°°-fpib+nc............."'1,'°-dna+nc....................

°0°................-fpib+nd0...............-dna+nc

心列黑列

用最小二乘可求得e的無(wú)偏估計(jì)，即A()、B()O

此法也可用來(lái)求C()。

b傳遞函數(shù)等價(jià)降階

由**式，可有：F()/D()=B()/A(),此說(shuō)明兩個(gè)

傳遞函數(shù)是等價(jià)的。對(duì)F()/D()施加激勵(lì)信號(hào)

{u(k)},可得輸出z(可二F()/D()u(k)。用最

小二乘法處理{z(k),u(k)},選擇合適的階

次，可得A()、B()的無(wú)偏估計(jì)。

Msjs

clearall;%最小二乘法forMISO

%-----------------------------------

loady6;

a=[lzta(l:m(l))'];%輸出

tt(:.l)=uyr(l:Y;%讀入

藪缸并賦叁(變量tt。仿真

tt(:,2)=uyr(2,:)';%tt的格b=[zta(m(l)+l:m(l)+m(2))'];

式為：第一列是系統(tǒng)輸出數(shù)據(jù),

其它列是對(duì)應(yīng)的輸入數(shù)據(jù)tf=filter(b,a,tt(:,2));

clearuyr;

ll=size(tt);%得到數(shù)據(jù)維數(shù)clearztataofffafayzz;

r=H(2)-l;%-----------------------------------

[zta,m,tao]=IslO(tt);%辨識(shí)

高階模型s一第二步，計(jì)算低階模型'

ls(tt);

多步最小二乘例一傳遞函數(shù)等價(jià)

4降階

有色噪聲模型：數(shù)據(jù)文件

y"6.mat00?.32w..(.左)

(1+L5qT+0.7q-2)y(k)=q-23.2u(k)+rr^77^

辨識(shí)結(jié)果(給定結(jié)構(gòu)：m=21，tao=02)

zta=

1.4938

0.6938

3.2194

相關(guān)最小—^乘

設(shè)模型為A()y(k)=B()u(k)+E(k),若u(k)與E(k)

不相關(guān)，則用u(k?j)乘模型中的各項(xiàng)并求期望,

得:用最小二乘法可得

A()RuyG)=B()Ru0),

A。、B()的無(wú)偏估計(jì)。

注意:使用相關(guān)最小二乘時(shí)，擾動(dòng)信號(hào)序列的周

期不要太長(zhǎng)，以保證由相關(guān)函數(shù)組成的模型在

求解時(shí)不發(fā)生病態(tài)。另外，計(jì)算相關(guān)函數(shù)時(shí)

不要以信號(hào)序列周期的整數(shù)倍來(lái)計(jì)算。

Cov_ls

tao(i+l)=input(‘輸入B(z)的時(shí)延')；

clearall;%相關(guān)最小二乘法forend

MISOn=m(D+max(tao);%算出一個(gè)方

程霰多使用的數(shù)據(jù)

loady6;

tt(:,l)=uyr(lj),;%讀入數(shù)據(jù),mz=2A7-l;

并賦給變量tt。lll=mz-n;%算出可列出的方程數(shù)

tt(:,2)=uyr(2,:),;%tt的格式inl=r+l;

為：第一列是系統(tǒng)輸出數(shù)據(jù)，其%計(jì)算相關(guān)函數(shù)

它列是對(duì)應(yīng)的輸入數(shù)據(jù)

fork=l:inl

clearuyr;forj=l:mz

%得到數(shù)據(jù)維數(shù)

ll=size(tt);rt(j,k)=O;

r=H(2)-l;

fori=j:j+ll(l)-mz

m(l)=inputC輸入A(z)的階次')；

0/0

tao(l)=0;%指定模型結(jié)構(gòu)代嗨矗疆隰甥擷

fori=l:r2人7-1,一個(gè)M序列周期

i%給出輸入編號(hào)end

m(i+l)=inputC輸入B(z)的階次)+1;rt(j,k)=rt(j,k)/(ll(l)-mz);

%躋次加一，表示套藏個(gè)數(shù)

end

end

kn=O;%構(gòu)造觀測(cè)數(shù)據(jù)矩陣ff%最小二乘求解

fork=l:inl%每一行中的變量循環(huán)

fori=l:lll%列循環(huán)

forj=l:m(k)%每行變量中的觀測(cè)數(shù)據(jù)循環(huán)

jtao=j+tao(k);%構(gòu)造時(shí)考慮時(shí)延

ifk>lff(i,j+kn)=rt(i+n-jtao+l,k);end

ifk==l,ff(i,j)=-rt(i+n-jtao,k)；end%輸出變量變號(hào)

end;end;

kn=kn+m(k);%算出行變量的啟始位置

end;

fori=l:lll

yy(i)=rt(i+n,l);%構(gòu)造輸出向量

end;________________________________________

fa=ff'*ff;%最小二乘計(jì)算

fay=ff'*yy';

zz=inv(fa);

zta=zz*fay%顯示參數(shù)和結(jié)構(gòu)

M,tao

相關(guān)最小二乘例

有色噪聲模型：數(shù)據(jù)文件

y"6.mat00?.32vp..(.左)

(l+1.5q1+O.7q-2)y(k)=q-23,2u(k)+7Tr7777^

辨識(shí)結(jié)果(給定結(jié)構(gòu)：m=21,tao=02)

zta=

1-4707

0.6711

3.2319

第四講辨識(shí)原理

-隨機(jī)逼近法

-模型參考自適應(yīng)辨識(shí)方法

-極大似然法

-預(yù)報(bào)誤差估計(jì)法

■Bayse估計(jì)

1極大驗(yàn)后參數(shù)估計(jì)法

2條件期望參數(shù)估計(jì)法

隨機(jī)逼近法

設(shè)/誤差e(k)是參數(shù)估計(jì)值6的函數(shù),參數(shù)辨識(shí)

問(wèn)題可通過(guò)極小化e(k)的方差來(lái)實(shí)現(xiàn)。即求參數(shù)

e使下列準(zhǔn)則函數(shù)最?。簀(e)=i/2E{e2(k)}。j(e)

的負(fù)梯度為：付號(hào)中(k)}o期果可求解

卜用=0,則可求嵋寥數(shù)的估計(jì)。但當(dāng)e(k)的分

布未知時(shí)，實(shí)際上是不可求解的。

在計(jì)算數(shù)學(xué)中，求二次函數(shù)的極小值常采用迭代法。首先

給出參數(shù)的一個(gè)估計(jì)值，以二次函數(shù)在該參數(shù)估計(jì)值處

的負(fù)梯度為修正方向，選取適當(dāng)?shù)牟介L(zhǎng)后，修正參數(shù)估

計(jì)值，直到收斂。

y隨機(jī)逼近法

仿此，我們有：e(k+i)=e(k)+p(k)［*［。如

果在求［-"時(shí)不求期望，則得到一個(gè)隨機(jī)的迭

代算法，稱之為隨機(jī)逼近法。

考慮線性回歸模型：y(k)=(pT(k)e+e(k),其中

e(k)是零均值隨機(jī)噪聲。

J(e)=1/2E{[y(k)-q)T(k)e]2},

「dJ1=E{(p(k)[y(k)-(pT(k)0]}□

de

A矗假定e(k)是各態(tài)遍歷的,則梯度為零可用

E{(p(k)[y(k)-(pT(k)9]]=(p(k)[y(k)-

L卜=i

(pT(k)ero來(lái)代替，由此得到了最小二乘法。

B應(yīng)用隨機(jī)逼近法,可得：

T

0(k+l)=0(k)+p(k)(p(k)[y(k)-(p(k)0(k)]o

p(k)的選取應(yīng)保證迭代收斂，可選取滿足如下條

件的p(k)：p(k)〉0且他。/)=0；

2P(左)=00；Z夕2(左)<8,例如p(k)=l/k；

例

C最準(zhǔn)則函數(shù)的二階導(dǎo)數(shù)(即海賽矩陣)之逆來(lái)

參與選擇修正方向，則稱為牛頓法:

0(k+l)=0(k)+p(k)R(ky1q)(k)[y(k)V(k)0(k)]

其中：

2

R(k)=—=E[(p(k)(pT(k)]

d6

=R(k-l)+p(k)[(p(k)(pT(k)-R(k-l)]

牛頓法的優(yōu)點(diǎn)是收斂速度快。

、模型參考自適應(yīng)辨識(shí)方法

考慮線性回歸模型：y(k)=(pT(k)9+e(k),其中e(k)

是零均值隨機(jī)噪聲。以輸出估計(jì)誤差為反饋信

號(hào)，以PI調(diào)節(jié)器的方式來(lái)修正參數(shù)：

e=eI+ep,eI(k)=eI(k-i)+p(p(k)E(k),

ep(k)=Q(p(k)e(k)

其中P為對(duì)稱正定矩陣，Q滿足P/2+Q>0

E(k)為廣義誤差，最簡(jiǎn)單的取法為

T

￡(k)=y(k)~(p(k)0(k~l)=￡0(k)

$極大似然法

輸出z是一個(gè)隨機(jī)變量，它的概率密度p（z|e）取決于參數(shù)

0O當(dāng)獲得觀測(cè)序列ZL={Z（1）,Z（2）〃..,Z（L）}T時(shí)，由該觀

測(cè)序列組成的聯(lián)合概率密度p（zje）應(yīng)當(dāng)取得最大值。

（當(dāng)一個(gè)隨機(jī)事件發(fā)生了，我們有理由相信，外部條

件一定處于使這隨機(jī)事件發(fā)生的概率最大時(shí)的狀態(tài)。）

那么，e的極大似然估計(jì)就是使p（Zje）|0ML=max的參

數(shù)估計(jì)值。

由于p（zje）中4已知，因而它只是參數(shù)e的函數(shù)，故

稱它為e的似然函數(shù)。有時(shí)也記作L（z」e）o

$極大似然法

極大似然原理可用下列等價(jià)的表示方式：

，P(Z"=0,「Slog

d3

°MLL」”

「”(z「aiog￡(z"。)]

-----二06二U

〔股兀L。L

求解極大似然估計(jì)的下一步是要給出p(zje)的

具體描述。

q獨(dú)立觀測(cè)情況

設(shè)z(l)，z(2)〃..，z(L)是一組在獨(dú)立觀測(cè)條件下獲得

的隨機(jī)序列，即各觀測(cè)值是互相獨(dú)立的，則

p(zje)可簡(jiǎn)化為：

P(zLie)=p(z(i)ie)p(z(2)ie)...p(z(L)ie)

=np(z(k)|B),其對(duì)數(shù)似然函數(shù)為：

k=\

L

L(Z|0)=Sinp(z(k)|0)

Lk—1o

《獨(dú)立觀測(cè)情況

設(shè)罪)?N(m,(j2),gp.

12

P(z(k)ie)=(2叩)2exp[-(C],負(fù)對(duì)數(shù)似然函數(shù)

、r2b

1L

-L(Z.|0)=[z(Q-"+Llno+(L/2)ln2n

」2b4=]L

當(dāng)。已知時(shí)，準(zhǔn)則函數(shù)就是：J(e)=[Z(A)—〃

/k=1

當(dāng)。未知時(shí)，可先由min0(e))求方及七所，再由

min(-L(Zje))求心

21￡

dJL八2m]2

(=---J.——Z[z(左)-

——―L(Z/I8))=2一=0'。minT

dacr3crLLk=i

4非獨(dú)立觀測(cè)

若z(1)，z(2),…,z(L)是非獨(dú)立觀測(cè)條件下獲得的隨

機(jī)序列，即觀測(cè)值z(mì)(k)是在已有觀測(cè)

z(l),z(2)〃..,z(k)的基礎(chǔ)上得到的，則p(Zje)應(yīng)按

條件概率的乘法規(guī)則寫成：

P(zLie)=p(z(L)|zL.1,e)p(zL.1ie),以此類推，有：

P(zLie)=iP(z(k)|zk.lze)o

n

k=\

I.非獨(dú)立觀測(cè)

設(shè)晶?N(rrikQk),并取mk為z(k)的條件均值:

￡(k)=E[z(k)|Zk4],即：

*八2

p(z(k)|zk_pe)=罟「],負(fù)對(duì)數(shù)為：

￡人2A

<(ZL|e)=工產(chǎn)>.+p]+(L/2)ln2n。

%二1左L

當(dāng)外三。已知時(shí)，就是：J(e)=工卬左)-￡")]2

當(dāng)Ok三。未知時(shí)，可先由min(J(e力求”及Jmin，

再由min(-L(Zf|0))求w。

Z

八221人2

o=-Jmin=一][2(左)一2(左)]

LL

例

考慮模型：z吟+千其中w?N(0Q2I),則

C2/CB

u

—Dz^N八(—DA爐1),

cCB2

[—z-u]

11DDA

(Iexp[--]，

P|0)=2

yj17i(y22b

_L1cCBCB

ln--------[—zu]r[—z-u]

J=lnL(0)=~~2g9L

2。~DDADDA

dJdJdJdJdJ

——=0

-----=0---------二0-----=0-------------二0

dAdBdCdBda2

u

(N

例

直彘解以上非線性方程組是困難的,如果已有參

數(shù)的某個(gè)估計(jì)，并用其構(gòu)成如下預(yù)濾波器：

人人人

*C*C*B*

z=~~TZ,u=，y=—u,

則上前前兩個(gè)方程卻可寫成：》--獷:心。，

這等價(jià)于輔助變量法。對(duì)于噪聲模

型，我們已有：

定義：--j,則川=以償=0。12，再次構(gòu)造預(yù)

濾波福v=v/D,w=w/D,則E［。-（0-1）皿*=°,

EC-9T）=0這也等價(jià)于輔助變量法。

4預(yù)報(bào)誤差估計(jì)法

對(duì)于模型結(jié)構(gòu)已知的系統(tǒng)，如果我們獲得參數(shù)的

某個(gè)