计算机课件系统辨识与建模2

系统辨识与建模

主讲：黄东

第三讲参数估计的批量法

■最小二乘算法

■参数估计的一般性质

■最小二乘、加权最小二乘估计性质

■噪声方差估计

■广义最小二乘

■偏倚校正算法

■辅助变量法

■多步最小二乘

■相关最小二乘

汇总的观

察误差

《最小二乘算法

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

w(k)为白噪声。假定模型的结构已知(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最小二乘算法

若我们的观测数据可写出N个这样的等式,

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)、…相关。当w(k+l)与w(k)、w(k-l)、

w(k-2)、…不相关时，y(k)与w(k+l)也不相关。

4最小二乘算法

以上结果等同于求使：

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

最小的e,因此称为最小二乘算法。

-=0

de

T

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

①T①-----------正贝U矩阵(对称、正定、可逆)

《加权最小二乘

若对各次观测数据加不同的权，即求使

J=Za(y(k)-(pT(k)0)2最小的0,则得到参数的

加a权最k小二乘估并：

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

A的对角元由cik构成

参数估计的一般性质

无偏性

如果E(3)=0,(6=G"或E(d)=6,则称估计为

无偏的。

2.有效性

如果无偏估计在满足Cov0)=M-i,则称估计为有效

的。其中：”日理―自理称为Fisher

信息矩阵，其逆M-i称为Crammer-Rao下界。在

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

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

4有效性例

对于z=HB+w,若噪声w为零均值、协方差阵

±w=E(WWT)的正态分布噪声，即：W〜N(ORw)，

则输出信号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

有效估计也称为最小方差估计、马尔可夫估计。

&参数估计的一般性质_

3.一致性

若估计©为渐近无偏的代叫E(3)=0),且

L巴Var)=0,则称0为一致估计。

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

最小二乘、加权最小二乘估计性质

估计性质最小二乘

无偏6=(①T①尸①TWN，当：1W(k)为零均值独立随机过

性

程(白噪声)时；或者2E(①丁班)力时，有：

E(^)=0o无偏条件较严格。

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

T生

二(①T①尸①TZw①(①T①)T训工若W是服从正态分

布的独立随机过程，则Zw二。：1，

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

若W是服从正态分布的独立随机过程，则

性27r

L四Cov⑴二例"芍尸二o,其中收敛于一个有

界非奇异矩阵。而Var(3)只是Cov")的对角线上元

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

最小二乘、加权最小二乘估计

性质

估计性质加权最小二乘

无偏0=(①TA①尸①TAWN

性

无偏条件同最小二乘估计。

有效

Cov(3)=(0TA0)-10TAzA①(①TA①尸,当

T生W

A二Zw"时，Cov(3)=(eTgw「l①)T=犷，此时不

要求W是服从正态分布的独立随机过程

一致若W是服从正态分布的独立随机过程，则

性

的Cov(^)=^7(7)-1=0,其中收敛于一个有

界非奇异矩阵。而VarG)只是Cov(3)的对角线上元

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

影响最小二乘计算结果的因素:

■u(k)QTQ是否可逆

■噪声w(k)02大，则的方差大

W0

白色零均值，0是无偏估计

■数据总量NN越大，J的方差越小

fork=l:inl%每一行中的变量循环

Ls.mfori=l:lll%列循环

function[zta,m,tao]=ls(tt)forj=l:m(k)%每行变量中的观测

数据循¥

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

第一列是系统输出数据，其它列jtao=j+tao(k);%构造时考虑时

是对应的输入数据延

ll=size(tt);%得到数据维数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

%指走模型结构%输出堤墓变号

tao(l)=0;end;end;

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

始位置

%给出输入编号

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

%跻次加一，表示蓊疑个数fori=l:lll

tao(i+l)=input(‘输入B⑵的时延')；yy(i)=tt(i+n,l);%构造输出向量

endend;

n=m(l)+max(tao);%算出一个fa=fP*ff;%最小二乘计算

方座最多使用的数据

fay=ff'*yy';

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

inl=r+l;%构造观测数据矩阵什zta=zz*fay%显示参数和结构

kn=0;m,tao

Lsl.m

clearall;%最小二乘法forMISO

loady3;

tt(:,l)=uyr(l,：y；%读入数据，并赋给变量tt。

tt(:,2)=uyr(2,:)1;%tt的格式为：第一列是系统

输出数据，其它列是对应的输入数据

clearuyr;

ls(tt);

《仿真例

1.无噪声模型：数据文件y3.mat

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

辨识结果(给定结构：m=21,tao=02)

zta=

1.5000

0.7000

3.2000

2.白噪声模型：数据文件y5-mat

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

辨识结果(给定结构：m=21,tao=02)

zta=

1.5027

0.7032

3.1935

3.有色噪声模型：数据文件y6.mat032

"0.32w（n左H）

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

辨识结果（给定结构：m=21，tao=02）

zta=

0.9757

0.1782

3.3955

4.白噪声模型b：数据文件y5b.mat

3.2/〃（左）

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

辨识结果(给定结构：m=21tao=02)

zta=

1.1627

0.3750

3.3907

追迹:对角线

元素之和

三噪声方差估计>/—

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)・\储2/(N・dime)

在有色噪声环境下，最小二乘估计是有偏的。下面的一些

算法对最小二乘进行改进。

广乂最小―^乘

11i

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

其中Mk)为白噪声。假定模型的结构已知

(n,m)o

如果噪声模型C(小)已知，显然用C(「)对输入/输出

数据进行滤波，则可得到满足最小二乘估计无偏

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

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

在C(小)未知时，我们可考虑采用迭代估计的方法

去求得。

4广义最小二乘的计算步骤

1令品(9)=1,i=0下标表示迭代次数；。令000000

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

3用最小二乘估计A(J)y(k)=BG13(k)+w(k)中的参数；

4用估计模型河q)、灯厂I)以及各时刻的观测数据，计算

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

5计算圻Ef(k)及如果&小于一定数，则结束

辨识。否则底下一步。

6对于噪声模型C(44k)=w(k),用最小二乘估计出参数,

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

以上算法的每一次循环中都要进行滤波和两次求逆。下面

的算法将在计算工作量上有所改进。

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的格式为：第量变导

一列是系统输出数据，其它列是对应

的输入数据end;end;

%m为各多项式中参数个数，应与tt的列kn=kn+m(k);%算出行

数一致；tao为时逅ll=size(tt);变量的官外位置

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

一个方盒最多使用的数据end;fori=l:lll

lll=ll(l)-n;%算出可列出的方程数yy(i)=tt(i+n,l)；%构造输出

向量

inl=ll(2);%构造观测数据矩阵ff

end;

kn=0;

fa=ff'*ff;%最小二乘计算

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

循环fay=ff'*yy';

fori=l:lll%列循环zz=inv(fa);

for1=1:m(k)%每行变量中zta=zz*fay%显示参数和

曲观测数徭循环结构

jtao=j+tao(k);%构造时考虑时延Mtao

fori=l:r

Gls.mi%给出输入编号

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

%阶次加一，表示参数个数

%广义最小二乘法

clearall;for输入的时延

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

loady6;End

%读入数据赋taomax=max(tao+m);

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

给tt。mc=input('输入噪声模型C(z)的阶次');

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

称式为：第一列是系统输出c=[l];d=[l];%初始化滤波器

数据，其它列是对应的输入lb=l;xci=100000;xci0=1000000;

数据whileabs(xci0-xci)>0.001%滤波

clearuy;plot(tt(:,l))循环

ll=size(tt);%得到数据维数xciO=xci;

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

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

%指虎模型结构[zta,m,tao]=lsO(ttl,mztao,ll);

最小二乘

tao(l)=0;

fork=l:taomax%计算输出估计

y(k)=tt(k,l);%设定初始输出%计算方程误差

endee=0;%计算损失函数值

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;%估计噪声模型

end;lc=size(e);

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

forj=l:m(i)c(2:mc+l)=cc总显示参数和

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

结构

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

elselb=lb+l;

fb(mm+j)=0;end%结束滤波循环

end____________plot(e),ylabel(cerror,),

end;end;xlabel('time')

mm=mm+m(i);

end;

y(k)=fb*zta;

end

广乂最小—^乘例

有色噪声模型：数据文件032G

Ly"6.mat0.32w(左)

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

辨识结果(给定结构：m=21,tao=02)

zta=

1.4950

0.6943

3.2075

c=[1.0000-1.71180.8069]

迭代次数：6

2.白噪声模型b：数据文件y5b.mat

（女）_

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

1+1.5「+0.7/

辨识结果（给定结构：m=21,tao=02）

zta=

1.5016

0.6992

3.1881

c=[1.0000-1.50671.4236-1.02150.5519-0.1689]

迭代次数：7

।W=Ay-Bu

T偏倚校正算法

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

w(k)为白噪声。令彳而=Urw(k),财

AG)y(k)=BG)u(k)+^(k)o分别写成回归模型：

Y=3+卅CC+W,组合起来有Y=[<D,C]+W/

「八i-iL。」

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

其最小二乘解为：7=利用分块矩

阵求逆公式有：

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

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

须要注意，c的求取仍然是一个迭代过程。

、偏倚校正算法的计算步骤

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

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

2计算j=丫-①©,并依据尸Co+W构造。，计算

D=QTMQO

-1T

3计算d=DQMY,并计算A6»=「及^=0LS-N。

4若参数已收敛，则结束辨识，否则转2。

以上算法的一次循环中没有滤波，且只有一次求逆。如果

将第3步中°的计算改为：。=（。％）-％上，则还可省去D

的计算。（这一改进由夏天长首先给出。）本法可能会

出现收敛慢的情况，可用对犷求均值来解决

mc=input('输入噪声模型C(z)的阶次')；

n=m(l)+max(tao);%胃由一个方连

Gls2.m最多使用的数据

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

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

loady6;inl=r+l;%构造观测数据矩阵什

%读入数据,kn=0;

fork=l:inl%每一行中的变量循环

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

第一列桌系统输出数据，其它列是:fori=l:lll%列循环

fori=l:m(k)%每行变量中的观测数

据循环

clearuyr;plot(tt(:,l))

%构造时考虑时延

ll=size(tt);%得到数据维数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%输出差及变号

%指走模型结构

tao(l)=0;end;end;

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

置

i%给出输入编号

end;

m(i+l)=inputC输入B(z)的阶次)+1;

%跻次加一，表示蓑藏个数fori=l:lll

%构造输出向量

tao(i+l)=input(‘输入B(z)的时延')；yy(i)=tt(i+n,l);

endend:

%最小二乘计算

end,end

fay=ff*yy';

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

zz=inv(fa);

fd=fc'*fc;%最小二乘计算

zaa=zz*ff';

fdcy=fc'*cy;

zta=zz*fay%显示参数和结构

zzc=inv(fd);

M,tao

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

ztab=zta-zta;显示参颠和结构_____________

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;%计算输出误差lb=lb+l;

e(mc+l:lll+mc)=yy'-ff*zta;end%结束滤波循环

%计算损失函数值

ee=O;

fork=l:lll+mc

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

end

xci=ee/(lll+mc)______________

%估计噪声模型

fori=l:mc

forj=l:lll

《偏倚校正算法例

工.有色噪声模型：数据文件y6.mat032G

"0.32)

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

辨识结果(给定结构：m=21,tao=02)

zta=

1.4879

0.6831

3.1976

c=[1.0000-1.70300.7983]

迭代次数：20

辅助变量法

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

望值为零得到参数的无偏估计。受此启发，若在Y=o)e+的合项

乘上中丁，使其满足以下两个条件：1UT的喇望值为零；2.w

可逆，则也可得到参数的无偏估计。下面讨论辅助变量的选取：

设模型为A()y(k)=B()u(k)+(k),若u(k)与(k)不相关：

一1-1

a选取辅助模型D()2(k)=F(《)u(k),用z(k)3u(k)构造肌

-1-1

b若系统的纯时延r已知，则可用u(k-T)、u(k)构造中；

c用u(k)、u(k)构造中

d(k)=D()w(k),^D()的阶次n已知,则可用y(k-n)、u(k)构造中；

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

()11也)计算出输出估计28),再用Z(k)、u(k)构造史；

--1

n=m(l)+max(tao);%算出一个方

Ivjs程盘多使用的效据

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

clearall;%辅助变量最小二乘法forinl=r+l;%构造观测数据矩阵"

MISOkn=0;

loady6;fork=l:inl%每一行中的变量循环

%读入数据,

m：牖盘黯fori=l:lll%列循环

fori=l:m(k)%每行变量中的观测数

我■匕费滁%输出数据%,tt其的它格列式是为:据循环

对应的输入数据jtao=j+tao(k);%构造时考虑时延

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

ll=size(tt);%得到数据维数jtao+l,k);end

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

r=H(2)-l;%输出z亚量变号z

m(l)=inputC输入A(z)的阶次')；

%指走模型结构end;end;

kn=kn+m(k);%算出行变量的启

tao(l)=0;始位置

fori=l:rend;

i%给出输入编号fori=l:lll

m(i+l)=inputC输入B(z)的阶次)+1;

%跻次加一，表示塞疑个数yy(i)=tt(i+n,l);%构造输出向量

tao(i+l)=inputC输入B(z)的时延')；end;

end

X

%f-1

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

匕

%最小二乘计算o+le

变,2

量I%以11%430)

为¥

fay=ff*yy';

zz=inv(fa);

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

zta=zz*fay%显示参数和结构%以岁助假型输出为*甫助变量

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);%构造输出向量

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

end;

kn=u；fa=fh'*ff:%辅助变量最小

fork=l:inl%每一行中的变量循环二乘计算

fori=l:lll%列循环fay=fh'*yy';

fori=l:m(k)%每行变量中的观测数

据循环zz=inv(fa);

zta=zz*fay%显示参数和结

jtao=j+tao(k);%构造时考虑时延构

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(k)为辅助变

量

《辅助变量法例

有色噪声模型：数据文件032

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

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

辨识结果(给定结构：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(k)为辅助变量

zta=[2.7082-1.85113.9692]

4.以u(k-tao)为辅助变量(方程病态，溢出)

5.以y(k・tao)为辅助变量

zta=[4.30792.91313.8725]

2-白噪声模型b：数据文件y5b.mat

-2

八、3.2qu(k)

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

1+1.5q+U.7q

辨识结果(给定结构：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解同次第方程组由**式，可得关系：

AU)F(「)=B(「)D(「),两边分别展开,并按J的同

次嘉相等规则，可列出\+找+加个方程：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的无偏估计，即A()、B()O

此法也可用来求C()。

b传递函数等价降阶

由**式，可有：F()/D()=B()/A(),此说明两个

传递函数是等价的。对F()/D()施加激励信号

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

小二乘法处理{z(k),u(k)},选择合适的阶

次，可得A()、B()的无偏估计。

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))'];

式为：第一列是系统输出数据,

其它列是对应的输入数据tf=filter(b,a,tt(:,2));

clearuyr;

ll=size(tt);%得到数据维数clearztataofffafayzz;

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

[zta,m,tao]=IslO(tt);%辨识

高阶模型s一第二步，计算低阶模型'

ls(tt);

多步最小二乘例一传递函数等价

4降阶

有色噪声模型：数据文件

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

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

辨识结果(给定结构：m=21，tao=02)

zta=

1.4938

0.6938

3.2194

相关最小—^乘

设模型为A()y(k)=B()u(k)+E(k),若u(k)与E(k)

不相关，则用u(k・j)乘模型中的各项并求期望,

得:用最小二乘法可得

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

A。、B()的无偏估计。

注意:使用相关最小二乘时，扰动信号序列的周

期不要太长，以保证由相关函数组成的模型在

求解时不发生病态。另外，计算相关函数时

不要以信号序列周期的整数倍来计算。

Cov_ls

tao(i+l)=input(‘输入B(z)的时延')；

clearall;%相关最小二乘法forend

MISOn=m(D+max(tao);%算出一个方

程霰多使用的数据

loady6;

tt(:,l)=uyr(lj),;%读入数据,mz=2A7-l;

并赋给变量tt。lll=mz-n;%算出可列出的方程数

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

为：第一列是系统输出数据，其%计算相关函数

它列是对应的输入数据

fork=l:inl

clearuyr;forj=l:mz

%得到数据维数

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;%指定模型结构代嗨矗疆隰甥撷

fori=l:r2人7-1,一个M序列周期

i%给出输入编号end

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

%跻次加一，表示套藏个数

end

end

kn=O;%构造观测数据矩阵ff%最小二乘求解

fork=l:inl%每一行中的变量循环

fori=l:lll%列循环

forj=l:m(k)%每行变量中的观测数据循环

jtao=j+tao(k);%构造时考虑时延

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

ifk==l,ff(i,j)=-rt(i+n-jtao,k)；end%输出变量变号

end;end;

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

end;

fori=l:lll

yy(i)=rt(i+n,l);%构造输出向量

end;________________________________________

fa=ff'*ff;%最小二乘计算

fay=ff'*yy';

zz=inv(fa);

zta=zz*fay%显示参数和结构

M,tao

相关最小二乘例

有色噪声模型：数据文件

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

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

辨识结果(给定结构：m=21,tao=02)

zta=

1-4707

0.6711

3.2319

第四讲辨识原理

-随机逼近法

-模型参考自适应辨识方法

-极大似然法

-预报误差估计法

■Bayse估计

1极大验后参数估计法

2条件期望参数估计法

随机逼近法

设/误差e(k)是参数估计值6的函数,参数辨识

问题可通过极小化e(k)的方差来实现。即求参数

e使下列准则函数最小：j(e)=i/2E{e2(k)}。j(e)

的负梯度为：付号中(k)}o期果可求解

卜用=0,则可求嵋寥数的估计。但当e(k)的分

布未知时，实际上是不可求解的。

在计算数学中，求二次函数的极小值常采用迭代法。首先

给出参数的一个估计值，以二次函数在该参数估计值处

的负梯度为修正方向，选取适当的步长后，修正参数估

计值，直到收敛。

y随机逼近法

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

果在求［-"时不求期望，则得到一个随机的迭

代算法，称之为随机逼近法。

考虑线性回归模型：y(k)=(pT(k)e+e(k),其中

e(k)是零均值随机噪声。

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

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

de

A矗假定e(k)是各态遍历的,则梯度为零可用

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

L卜=i

(pT(k)ero来代替，由此得到了最小二乘法。

B应用随机逼近法,可得：

T

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

p(k)的选取应保证迭代收敛，可选取满足如下条

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

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

例

C最准则函数的二阶导数(即海赛矩阵)之逆来

参与选择修正方向，则称为牛顿法:

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(k)=(pT(k)9+e(k),其中e(k)

是零均值随机噪声。以输出估计误差为反馈信

号，以PI调节器的方式来修正参数：

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

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

其中P为对称正定矩阵，Q满足P/2+Q>0

E(k)为广义误差，最简单的取法为

T

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

$极大似然法

输出z是一个随机变量，它的概率密度p（z|e）取决于参数

0O当获得观测序列ZL={Z（1）,Z（2）〃..,Z（L）}T时，由该观

测序列组成的联合概率密度p（zje）应当取得最大值。

（当一个随机事件发生了，我们有理由相信，外部条

件一定处于使这随机事件发生的概率最大时的状态。）

那么，e的极大似然估计就是使p（Zje）|0ML=max的参

数估计值。

由于p（zje）中4已知，因而它只是参数e的函数，故

称它为e的似然函数。有时也记作L（z」e）o

$极大似然法

极大似然原理可用下列等价的表示方式：

，P(Z"=0,「Slog

d3

°MLL」”

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

-----二06二U

〔股兀L。L

求解极大似然估计的下一步是要给出p(zje)的

具体描述。

q独立观测情况

设z(l)，z(2)〃..，z(L)是一组在独立观测条件下获得

的随机序列，即各观测值是互相独立的，则

p(zje)可简化为：

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

=np(z(k)|B),其对数似然函数为：

k=\

L

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

Lk—1o

《独立观测情况

设罪)〜N(m,(j2),gp.

12

P(z(k)ie)=(2叩)2exp[-(C],负对数似然函数

、r2b

1L

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

」2b4=]L

当。已知时，准则函数就是：J(e)=[Z(A)—〃

/k=1

当。未知时，可先由min0(e))求方及七所，再由

min(-L(Zje))求心

21£

dJL八2m]2

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

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

dacr3crLLk=i

4非独立观测

若z(1)，z(2),…,z(L)是非独立观测条件下获得的随

机序列，即观测值z(k)是在已有观测

z(l),z(2)〃..,z(k)的基础上得到的，则p(Zje)应按

条件概率的乘法规则写成：

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

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

n

k=\

I.非独立观测

设晶〜N(rrikQk),并取mk为z(k)的条件均值:

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

*八2

p(z(k)|zk_pe)=罟「],负对数为：

£人2A

<(ZL|e)=工产>.+p]+(L/2)ln2n。

%二1左L

当外三。已知时，就是：J(e)=工卬左)-£")]2

当Ok三。未知时，可先由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

例

直彘解以上非线性方程组是困难的,如果已有参

数的某个估计，并用其构成如下预滤波器：

人人人

*C*C*B*

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

则上前前两个方程却可写成：》--犷:心。，

这等价于辅助变量法。对于噪声模

型，我们已有：

定义：--j,则川=以偿=0。12，再次构造预

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

EC-9T）=0这也等价于辅助变量法。

4预报误差估计法

对于模型结构已知的系统，如果我们获得参数的

某个