《系统辨识与建模》由会员分享，可在线阅读，更多相关《系统辨识与建模（72页珍藏版）》请在金锄头文库上搜索。
1、系统辨识与建模系统辨识与建模 1第三讲 参数估计的批量法 n最小二乘算法 n参数估计的一般性质 n最小二乘、加权最小二乘估计性质 n噪声方差估计 n广义最小二乘 n偏倚校正算法 n辅助变量法 n多步最小二乘 n相关最小二乘 2最小二乘算法考虑差分方程：A( )y(k)=B( )u(k)+ w(k)，其中w(k)为白噪声。假定模型的结构已知（n,m,），将其写成线性回归模型： y(k)=-a1y(k-1)-any(k-n)+b1u(k-)+ +bmu(k-m-+1)+ w(k) =a1, , an,b1 , bmT (k)=-y(k-1),-y(k-n),u(k-),u(k-m-+1)T y(k
2、)=T(k)+w(k) （以下令=1）汇总的观察误差3y(k)=-a1y(k-1)-any(k-n)+b1u(k-1)+bmu(k-m)+ w(k)y(k+1)=-a1y(k)-any(k-n+1)+b1u(k)+bmu(k-m+1)+w(k+1)y(k+2)=-a1y(k+1)-any(k-n+2)+b1u(k+1)+bmu(k-m+2)+w(k+2) T(k)=- y(k-1)-y(k-n) u(k-1)u(k-m)T(k+1)=- y(k)-y(k-n+1) u(k)u(k-m+1)T(k+2)=- y(k+1)-y(k-n+2) u(k+1)u(k-m+2) T(k+N)=- y(k+
3、N-1)-y(k-n+N) u(k+N-1)u(k-m+N)4最小二乘算法若我们的观测数据可写出N个这样的等式， YN=+WN ， T=(k)，(N+k-1) T YN=T+T WN =(T)-1T YN - (T)-1T WN LS=(T)-1TYN 条件1：E(T WN)=0 （当w(k)为白色时,条件1满足）条件2: T可逆 5当w(k)为白色时,条件1满足T WN= - y(k-1) -y(k) -y(k+1) -y(k+N-2) w(k) - y(k-2) -y(k-1) -y(k) -y(k+N-1) w(k+1) . - y(k-n) -y(k-n+1) -y(k-n+2)-y(
4、k-n+N-1) . u(k-1) u(k) u(k+1) u(k+N-2) . . u(k-m) u(k-m+1) u(k-m+2) u(k-m+N-1) w(k+N-1)y(k)与w(k)、w(k-1)、w(k-2)、相关。相关。当当w(k+1)与w(k)、w(k-1)、w(k-2)、不相关时，相关时， y(k)与w(k+1)也不相关相关。6最小二乘算法以上结果等同于求使： J=(y(k)-T(k)2最小的，因此称为最小二乘算法。 =0T YN=T-正则方程T-正则矩阵(对称、正定、可逆)7加权最小二乘加权最小二乘 若 对 各 次 观 测 数 据 加 不 同 的 权 ， 即 求 使J=k(
5、y(k)-T(k)2最小的,则得到参数的加权最小二乘估计：LSLS=(=(T T)-1-1T TYYN N 的对角元由k构成8参数估计的一般性质 1无偏性 如果E( )=0,( = - ) 或E( )= ，则称估计为无偏的。 2有效性如果无偏估计 满足Cov( )=M-1，则称估计为有效的。其中： 称为Fisher信息矩阵，其逆M-1称为Crammer-Rao 下界。在一般情况下，有Crammer-Rao不等式：Cov( )=E( - )( - )TM-1 9有效性例对于z=H+w，若噪声w为零均值、协方差阵w=E(WWT)的正态分布噪声，即：WN(0,w)，则输出信号ZN(E(H),w)，即
6、，因此： 于是，M=E(HTw-1H)， 其 Crammer-Rao不 等 式 为 ： Cov( )E(HTw-1H)-1。有效估计也称为最小方差估计、马尔可夫估计。 10参数估计的一般性质 3一致性若估计 为渐近无偏的( E( )=0），且 Var( )=0，则称 为一致估计。 Var( )=E( - )T( - ) 11最小二乘、加权最小二乘估计性质12最小二乘、加权最小二乘估计性质13影响最小二乘计算结果的因素:nu(k) T T是否可逆是否可逆 n噪声噪声w(k) 大大,则则 的方差大的方差大 白色零均值白色零均值, 是无偏估计是无偏估计n数据总量数据总量N N越大越大, 的方差越小的
7、方差越小14最小二乘算法的MATLAB程序 读入数据读入数据读入结构读入结构构造矩阵构造矩阵和和Y Y计算计算T T和和T TY Y计算计算Ls15function zta,m,tao=ls(tt)%最小二乘法for MISO, tt的格式为：第一列是系统输出数据，其它列是对应的输入数据ll=size(tt); %得到数据维数r=ll(2)-1; m(1)=input(输入 A(z)的阶次); %指定模型结构tao(1)=0;for i=1:r i %给出输入编号m(i+1)=input(输入 B(z)的阶次)+1; %阶次加一，表示参数个数tao(i+1)=input(输入B(z)的时延);
8、endn=m(1)+max(tao); %算出一个方程最多使用的数据lll=ll(1)-n; %算出可列出的方程数in1=r+1; %构造观测数据矩阵ffkn=0;for k=1:in1 %每一行中的变量循环for i=1:lll %列循环for j=1:m(k) %每行变量中的观测数据循环jtao=j+tao(k); %构造时考虑时延if k1 ff(i,j+kn)=tt(i+n-jtao+1,k);endif k=1, ff(i,j)=-tt(i+n-jtao,k);end %输出变量变号end;end;kn=kn+m(k); %算出行变量的启始位置end;for i=1:lllyy(i)
9、=tt(i+n,1); %构造输出向量end;fa=ff*ff; %最小二乘计算fay=ff*yy;zz=inv(fa);zta=zz*fay %显示参数和结构m，taoLs.m16Ls1.mclear all; %最小二乘法for MISOload y3;tt(:,1)=uyr(1,:); %读入数据，并赋给变量tt。tt(:,2)=uyr(2,:); %tt的格式为：第一列是系统输出数据，其它列是对应的输入数据clear uyr;plot(tt(:,1) ls(tt);17仿真例1. 无噪声模型：数据文件 y3.mat(1+1.5q-1+0.7q-2)y(k)=q-23.2u(k)辨识结果
10、辨识结果（给定结构：（给定结构：m =2 1，tao = 0 2）zta = 1.5000 0.7000 3.2000182. 白噪声模型：数据文件白噪声模型：数据文件 y5.mat(1+1.5q-1+0.7q-2)y(k)=q-23.2u(k)+w(k)辨识结果辨识结果（给定结构：（给定结构：m =2 1，tao = 0 2）zta = 1.5027 0.7032 3.1935193. 有色噪声模型：数据文件有色噪声模型：数据文件 y6.mat(1+1.5q-1+0.7q-2)y(k)=q-23.2u(k)+辨识结果辨识结果（给定结构：（给定结构：m =2 1，tao =0 2）zta =
11、0.9757 0.1782 3.3955204. 白噪声模型白噪声模型b：数据文件：数据文件 y5b.maty(k)= + w(k)辨识结果辨识结果（给定结构：（给定结构：m =2 1，tao =0 2）zta = 1.1627 0.3750 3.390721噪声方差估计e=YN- =+WN- =+WN-(T)-1TYN=WN-(T)-1T WN =(I-(T)-1T) WN=B WN， 因为因为BT=B ,B2=B所以，所以，E(eTe)=E(WN TB WN) = E(WN T WN- WN T(T)-1T WN) = N- (Tr(T)-1T)= (N-dim) =(y(k)- )2/(
12、N-dim)在有色噪声环境下，最小二乘估计是有偏的。下面的一些在有色噪声环境下，最小二乘估计是有偏的。下面的一些算法对最小二乘进行改进。算法对最小二乘进行改进。 追迹:对角线元素之和22广义最小二乘考虑差分方程：A( )y(k)=B( )u(k)+ w(k)，其中w(k)为白噪声。假定模型的结构已知（n,m）。如果噪声模型C( )已知，显然用C( )对输入/输出数据进行滤波，则可得到满足最小二乘估计无偏条件的模型：A( ) (k)= B( ) (k)+ w(k)，其中： (k) =C( )y(k)， (k) =C( )u(k)。在C( )未知时，我们可考虑采用迭代估计的方法去求得。 23广义最
13、小二乘的计算步骤 1令令C0()=1，i=0下标表示迭代次数；下标表示迭代次数；0=10000002计算计算(k)=Ci()y(k)，(k)=Ci()u(k)；i=i+1;3用最小二乘估计用最小二乘估计A()(k)=B()(k)+w(k)中的参数；中的参数；4用用估估计计模模型型、以以及及各各时时刻刻的的观观测测数数据据，计计算算出残差出残差:(k)=()y(k)-()u(k)5计计算算i=2(k)及及=i-1-i，如如果果小小于于一一定定数数，则则结结束束辨识。否则转下一步。辨识。否则转下一步。6对对于于噪噪声声模模型型C()(k)=w(k)，用用最最小小二二乘乘估估计计出出参参数数，得到更
14、新的得到更新的Ci()后，返回后，返回2。以上算法的每一次循环中都要进行滤波和两次求逆。下面以上算法的每一次循环中都要进行滤波和两次求逆。下面的算法将在计算工作量上有所改进。的算法将在计算工作量上有所改进。24function zta,m,tao=ls0(tt,m,tao,ll)%最小二乘法for MISO,tt的格式为：第一列是系统输出数据，其它列是对应的输入数据%m为各多项式中参数个数，应与tt的列数一致；tao为时延；ll=size(tt);n=max(m)+max(tao); %算出一个方程最多使用的数据lll=ll(1)-n; %算出可列出的方程数in1=ll(2); %构造观测数据
15、矩阵ffkn=0;for k=1:in1 %每一行中的变量循环for i=1:lll %列循环for j=1:m(k) %每行变量中的观测数据循环jtao=j+tao(k); %构造时考虑时延if k1, ff(i,j+kn)=tt(i+n-jtao+1,k);endif k=1, ff(i,j)=-tt(i+n-jtao,k);end %输出变量变号end;end;kn=kn+m(k); %算出行变量的启始位置end; for i=1:lllyy(i)=tt(i+n,1); %构造输出向量end;fa=ff*ff; %最小二乘计算fay=ff*yy;zz=inv(fa);zta=zz*fay
16、 %显示参数和结构M taoLs0.m25clear all; %广义最小二乘法for MISOload y6;tt(:,1)=uyr(1,:); %读入数据赋给tt。tt(:,2)=uyr(2,:); %tt的格式为：第一列是系统输出数据，其它列是对应的输入数据clear uy; plot(tt(:,1)ll=size(tt); %得到数据维数in1=ll(2); r=in1-1;m(1)=input(输入 A(z)的阶次); %指定模型结构tao(1)=0; for i=1:r i %给出输入编号m(i+1)=input(输入 B(z)的阶次)+1; %阶次加一，表示参数个数tao(i+1
17、)=input(输入B(z)的时延);Endtaomax=max(tao+m);mc=input(输入噪声模型C(z)的阶次); c=1;d=1; %初始化滤波器lb=1;xci=100000;xci0=1000000;while abs(xci0-xci)0.001 %滤波循环xci0=xci;tt1=filter(c,d,tt); %输入输出滤波.zta,m,tao=ls0(tt1,m,tao,ll); %主最小二乘Gls.m26for k=1:taomax %计算输出估计 y(k)=tt(k,1); %设定初始输出 endfor k=1+taomax:ll(1)mm=0;for i=1:
18、in1 if i=1 for j=1:m(i) fb(mm+j)=-tt(k-j,1); end; else for j=1:m(i) if k-tao(i)-j+10 fb(mm+j)=tt(k-tao(i)-j+1,i); else fb(mm+j)=0; end end;end; mm=mm+m(i);end;y(k)=fb*zta; %用LS估计参数求输出估计ende=tt(:,1)-y; %计算方程误差ee=0; %计算损失函数值for k=1:ll(1) ee=ee+e(k)*e(k);endxci=ee/ll(1)taoc=0; %估计噪声模型lc=size(e);cc,mc,t
19、aoc=ls0(e,mc,taoc,lc);c(2:mc+1)=cc %显示参数和结构lb=lb+1;end %结束滤波循环plot(e)，ylabel(error)，xlabel(time)27广义最小二乘例1. 有色噪声模型：数据文件有色噪声模型：数据文件 y6.mat(1+1.5q-1+0.7q-2)y(k)=q-23.2u(k)+辨识结果辨识结果（给定结构：（给定结构：m =2 1，tao =0 2）zta = 1.4950 0.6943 3.2075c = 1.0000 -1.7118 0.8069迭代次数：迭代次数： 6282. 白噪声模型白噪声模型b：数据文件：数据文件 y5b.
20、maty(k)= + w(k)辨识结果辨识结果（给定结构：（给定结构：m =2 1，tao =0 2）zta = 1.5016 0.6992 3.1881c = 1.0000 -1.5067 1.4236 -1.0215 0.5519 -0.1689迭代次数：迭代次数： 729偏倚校正算法仍仍考考虑虑差差分分方方程程：A( )y(k)=B( )u(k)+ w(k)，其其中中w(k)为白噪声。令为白噪声。令(k)= w(k)，则，则A( )y(k)= B( )u(k)+(k)。分别写成回归模型：。分别写成回归模型：Y=+,=C+W,组合起来有组合起来有Y=, +W，其其最小二乘解为：最小二乘解为
21、： = ，利用分块矩，利用分块矩阵求逆公式有：阵求逆公式有： =(T)-1TY-(T)-1T =D-1 TMY M=I-(T)-1T D= TM须要注意，须要注意， 的求取仍然是一个迭代过程。的求取仍然是一个迭代过程。 =Ay-Bu30偏倚校正算法的计算步骤 1 令C( )=1，用最小二乘法求=LS=(T)-1TY，并保留、=(T)-1T以及M=I-。2 计算 =Y- ，并依据 =C +W构造 ，计算D= TM 。3 计算 =D-1 TMY，并计算 = 及 =LS- 。4 若参数已收敛，则结束辨识，否则转2。以上算法的一次循环中没有滤波，且只有一次求逆。如果将第3步中 的计算改为： =( T
22、)-1 T ，则还可省去D的计算。（这一改进由夏天长首先给出。）本本法法可可能能会会出现收敛慢的情况，可用对出现收敛慢的情况，可用对 求均值来解决求均值来解决31Gls2.mclear all; %偏倚校正法for MISOload y6;tt(:,1)=uyr(1,:); %读入数据，并赋给变量tt。tt(:,2)=uyr(2,:); %tt的格式为：第一列是系统输出数据，其它列是对应的输入数据clear uyr;plot(tt(:,1)ll=size(tt); %得到数据维数r=ll(2)-1;m(1)=input(输入 A(z)的阶次); %指定模型结构tao(1)=0;for i=1:
23、ri %给出输入编号m(i+1)=input(输入 B(z)的阶次)+1; %阶次加一，表示参数个数tao(i+1)=input(输入B(z)的时延);endmc=input(输入噪声模型C(z)的阶次); n=m(1)+max(tao); %算出一个方程最多使用的数据lll=ll(1)-n; %算出可列出的方程数lb=1;xci=100000;xci0=1000000;c=1;in1=r+1; %构造观测数据矩阵ffkn=0;for k=1:in1 %每一行中的变量循环for i=1:lll %列循环for j=1:m(k) %每行变量中的观测数据循环jtao=j+tao(k); %构造时考
24、虑时延if k1 ff(i,j+kn)=tt(i+n-jtao+1,k);endif k=1, ff(i,j)=-tt(i+n-jtao,k);end %输出变量变号end;end;kn=kn+m(k); %算出行变量的启始位置end; for i=1:lllyy(i)=tt(i+n,1); %构造输出向量end;32fa=ff*ff; %最小二乘计算fay=ff*yy;zz=inv(fa);zaa=zz*ff;zta=zz*fay %显示参数和结构M,taoztab=zta-zta;ztaa=zta;while abs(xci0-xci)0.05 xci0=xci;e(1:mc)=0; %计
25、算输出误差e(mc+1:lll+mc)=yy-ff*zta;%计算损失函数值ee=0;for k=1:lll+mc ee=ee+e(k)*e(k);endxci=ee/(lll+mc)%估计噪声模型for i=1:mcfor j=1:lllfc(j,i)=-e(j+mc-i);end,endcy=e(mc+1:lll+mc);fd=fc*fc; %最小二乘计算fdcy=fc*cy;zzc=inv(fd);c(2:mc+1)=zzc*fdcy %显示参数和结构ztab=ztab+zaa*fc*c(2:mc+1);zta=ztaa-ztab/lb%zta=ztaa-zaa*fc*c(2:mc+1)
26、lb=lb+1;end %结束滤波循环33偏倚校正算法例1. 有色噪声模型：数据文件有色噪声模型：数据文件 y6.mat(1+1.5q-1+0.7q-2)y(k)=q-23.2u(k)+辨识结果辨识结果（给定结构：（给定结构：m =2 1，tao =0 2）zta = 1.4879 0.6831 3.1976c =1.0000 -1.7030 0.7983迭代次数：迭代次数： 2034辅助变量法分析最小二乘法中，在Y=+W的各项乘上T,然后利用T W的期望值为零得到参数的无偏估计。受此启发，若在Y=+ 的各项乘上T，使其满足以下两个条件：1.T 的期望值为零；2. T可逆，则也可得到参数的无偏
27、估计。下面讨论辅助变量的选取：设模型为A( )y(k)= B( )u(k)+ (k)，若u(k)与 (k)不相关：a 选取辅助模型D( )z(k)= F( )u(k)，用z(k)、u(k)构造；b 若系统的纯时延已知，则可用u(k-)、u(k)构造；c 用u(k)、u(k)构造d (k)=D( )w(k),若D( )的阶次n已知,则可用y(k-n)、u(k)构造；e 先求出最小二乘解LS=(T)-1TY ，然后依据 ( )z(k)= ( )u(k)计算出输出估计z(k)，再用z(k)、u(k)构造； 35clear all; %辅助变量最小二乘法for MISOload y6;tt(:,1)=
28、uyr(1,:); %读入数据，并赋给变量tt。tt(:,2)=uyr(2,:); %tt的格式为：第一列是系统输出数据，其它列是对应的输入数据clear uyr;ll=size(tt); %得到数据维数r=ll(2)-1;m(1)=input(输入 A(z)的阶次); %指定模型结构tao(1)=0;for i=1:ri %给出输入编号m(i+1)=input(输入 B(z)的阶次)+1; %阶次加一，表示参数个数tao(i+1)=input(输入B(z)的时延);endn=m(1)+max(tao); %算出一个方程最多使用的数据lll=ll(1)-n; %算出可列出的方程数in1=r+1
29、; %构造观测数据矩阵ffkn=0;for k=1:in1 %每一行中的变量循环for i=1:lll %列循环for j=1:m(k) %每行变量中的观测数据循环jtao=j+tao(k); %构造时考虑时延if k1 ff(i,j+kn)=tt(i+n-jtao+1,k);endif k=1, ff(i,j)=-tt(i+n-jtao,k);end %输出变量变号end;end;kn=kn+m(k); %算出行变量的启始位置end;for i=1:lllyy(i)=tt(i+n,1); %构造输出向量end;Iv_ls36fa=ff*ff; %最小二乘计算fay=ff*yy;zz=inv(
30、fa);zta=zz*fay %显示参数和结构%建立辅助变量%a=1 zta(1:m(1);%b=0 zta(m(1)+1:m(1)+m(2);a=1 1.7 0.72;b=0 0 1;tf=filter(b,a,tt(:,2);kn=0;for k=1:in1 %每一行中的变量循环for i=1:lll %列循环for j=1:m(k) %每行变量中的观测数据循环jtao=j+tao(k); %构造时考虑时延if k1 fh(i,j+kn)=tt(i+n-jtao+1,k);end%if k=1, fh(i,j)=-tt(i+n-jtao+1,2);end %以u(k)为辅助变量%if k=
31、1, fh(i,j)=-tt(i+n-jtao-tao(2)+1,2);end %以u(k-tao)为辅助变量if k=1, fh(i,j)=-tf(i+n-j);end %以辅助模型输出为辅助变量end;end;kn=kn+m(k); %算出行变量的启始位置end; for i=1:lllyy(i)=tt(i+n,1); %构造输出向量end;fa=fh*ff; %辅助变量最小二乘计算fay=fh*yy;zz=inv(fa);zta=zz*fay %显示参数和结构M，tao37辅助变量法例有色噪声模型：数据文件有色噪声模型：数据文件 y6.mat(1+1.5q-1+0.7q-2)y(k)=q
32、-23.2u(k)+辨识结果辨识结果（给定结构：（给定结构：m =2 1，tao =0 2）1. 使用辅助模型：使用辅助模型：a=1 1.7 0.72;b=0 0 1;zta = 1.4956 0.6962 3.2234382. 以以LS为辅助模型：为辅助模型： b=0 0 3.3955; a=1 0.9757 0.1782;zta =1.4908 0.6895 3.22353.以以u(k)为辅助变量为辅助变量X zta = 2.7082 -1.8511 3.96924.以以u(k-tao)为辅助变量（方程病态，溢出）为辅助变量（方程病态，溢出）X5.以以y(k-tao)为辅助变量为辅助变量X
33、 zta =4.3079 2.9131 3.8725392. 白噪声模型白噪声模型b：数据文件：数据文件 y5b.maty(k)= + w(k)辨识结果辨识结果（给定结构：（给定结构：m =2 1，tao =0 2）1. 使用辅助模型：使用辅助模型：a=1 1.7 0.72;b=0 0 1;zta =1.4740 0.6845 3.3692 2. 以以LS为辅助模型：为辅助模型： b=0 0 3.3907; a=1 1.1627 0.3750;zta =1.4456 0.6389 3.377740多步最小二乘考虑差分方程：A()y(k)= B()u(k)+ w(k)，其中w(k)为白噪声。模型
34、可改写为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()。这是第一步。第二步可有两种不同方法：a a 解解同同次次幂幂方方程程组组 由*式，可得关系：A( )F( )=B( )D( )，两边分别展开，并按 的同次幂相等规则，可列出na+nb+nc个方程：F=H，其中=a1,.ana,b1,.bnbT，F=f1,.fnb+nc,0,.0 41 0 0 0. . 0 1 0 0. . 0 0 -f1 0 0. . 0 -d1 1 0. .0 0H= -f2 f
35、1 0. . .0 -d2 -d1 1. .0 0 . . . . .-d1 1 . -f2 -f1 . . .-d2 d1 na+nb+nc 行 . . -f3 -f2. . . . -fnb+nc. . 0 -fnb+nc. . -dna+nc. 0 0 -fnb+nc. . 0 -dna+nc. . 00 0. . -fnb+nc 0. .-dna+nc n na a列列 n nb b列列 42用最小二乘可求得的无偏估计，即A()、B()。此法也可用来求C()。b 传递函数等价降阶传递函数等价降阶 由*式，可有： F()/D()=B()/A()，此说明两个传递函数是等价的。对F()/D(
36、)施加激励信号u(k)，可得输出z(k)= F()/D()u(k)。用最小二乘法处理 z(k)，u(k)，选择合适的阶次，可得A()、B()的无偏估计。 43clear all; %最小二乘法for MISOload y6;tt(:,1)=uyr(1,:); %读入数据，并赋给变量tt。tt(:,2)=uyr(2,:); %tt的格式为：第一列是系统输出数据，其它列是对应的输入数据clear uyr;ll=size(tt); %得到数据维数r=ll(2)-1;zta,m,tao=ls10(tt); %辨识高阶模型%-a=1 zta(1:m(1); %输出仿真b=zta(m(1)+1:m(1)+
37、m(2);tf=filter(b,a,tt(:,2);tt(:,1)=tf;clear zta tao ff fa fay zz;%-s=第二步，计算低阶模型ls(tt);Ms_ls44多步最小二乘例-传递函数等价降阶 有色噪声模型：数据文件有色噪声模型：数据文件 y6.mat(1+1.5q-1+0.7q-2)y(k)=q-23.2u(k)+辨识结果辨识结果（给定结构：（给定结构：m =2 1，tao =0 2）zta = 1.4938 0.6938 3.219445相关最小二乘设模型为A()y(k)=B()u(k)+E(k)，若u(k)与E(k)不相关，则用u(k-j)乘模型中的各项并求期望
38、，得：A()Ruy(j)= B()Ru(j)，用最小二乘法可得A()、B()的无偏估计。注意:使用相关最小二乘时,扰动信号序列的周期不要太长,以保证由相关函数组成的模型在求解时不发生病态。另外，计算相关函数时不要以信号序列周期的整数倍来计算。46clear all; %相关最小二乘法for MISOload y6;tt(:,1)=uyr(1,:); %读入数据，并赋给变量tt。tt(:,2)=uyr(2,:); %tt的格式为：第一列是系统输出数据，其它列是对应的输入数据clear uyr;ll=size(tt); %得到数据维数r=ll(2)-1;m(1)=input(输入 A(z)的阶次)
39、;tao(1)=0; %指定模型结构for i=1:ri %给出输入编号m(i+1)=input(输入 B(z)的阶次)+1; %阶次加一，表示参数个数tao(i+1)=input(输入B(z)的时延);endn=m(1)+max(tao); %算出一个方程最多使用的数据mz=27-1;lll=mz-n; %算出可列出的方程数in1=r+1;%计算相关函数for k=1:in1for j=1:mzrt(j,k)=0;for i=j:j+ll(1)-mzrt(j,k)=rt(j,k)+tt(i,k)*tt(i-j+1,2); %所有变量对U1求相关，长度为27-1，一个M序列周期endrt(j,
40、k)=rt(j,k)/(ll(1 )-mz); endendCov_ls47kn=0; %构造观测数据矩阵ff %最小二乘求解for k=1:in1 %每一行中的变量循环for i=1:lll %列循环for j=1:m(k) %每行变量中的观测数据循环jtao=j+tao(k); %构造时考虑时延if k1 ff(i,j+kn)=rt(i+n-jtao+1,k);endif k=1, ff(i,j)=-rt(i+n-jtao,k);end %输出变量变号end;end;kn=kn+m(k); %算出行变量的启始位置end;for i=1:lllyy(i)=rt(i+n,1); %构造输出向量
41、end;fa=ff*ff; %最小二乘计算fay=ff*yy;zz=inv(fa);zta=zz*fay %显示参数和结构M，tao48相关最小二乘例有色噪声模型：数据文件有色噪声模型：数据文件 y6.mat(1+1.5q-1+0.7q-2)y(k)=q-23.2u(k)+辨识结果辨识结果（给定结构：（给定结构：m =2 1，tao =0 2）zta = 1.4707 0.6711 3.231949第四讲 辨识原理 n随机逼近法n模型参考自适应辨识方法 n极大似然法n预报误差估计法nBayse估计1 极大验后参数估计法 2 条件期望参数估计法 50随机逼近法设广义误差e(k)是参数估计值的函数
42、，参数辨识问题可通过极小化e(k)的方差来实现。即求参数使下列准则函数最小：J()=1/2Ee2(k)。J()的负梯度为： =E-e(k) 。如果可求解 =0，则可求得参数的估计。但当e(k)的分布未知时，实际上是不可求解的。 在计算数学中，求二次函数的极小值常采用迭迭代代法法。首先给出参数的一个估计值，以二次函数在该参数估计值处的负梯度为修正方向，选取适当的步长后，修正参数估计值，直到收敛。51随机逼近法仿此，我们有：(k+1)= (k)+(k) 。如果在求如果在求 时不求期望，则得到一个随机的迭时不求期望，则得到一个随机的迭代算法，称之为随机逼近法。代算法，称之为随机逼近法。考虑线性回归模
43、型：y(k)=T(k)+e(k)，其中e(k)是零均值随机噪声。J()=1/2E y(k)-T(k)2, =E(k) y(k)-T(k)。 52例A 如果假定e(k)是各态遍历的，则梯度为零可用E(k) y(k)-T(k)= (k) y(k)-T(k)=0来代替，由此得到了最小二乘法最小二乘法。B 应用随机逼近法随机逼近法，可得：(k+1)=(k)+(k)(k)y(k)-T(k)(k)。(k)的选取应保证迭代收敛，可选取满足如下条件的(k)： (k)0且 ； ； ，例如(k)=1/k； 53例C 若以准则函数的二阶导数（即海赛矩阵）之逆来参与选择修正方向，则称为牛顿法牛顿法：(k+1)=(k)
44、+(k)R(k)-1(k)y(k)-T(k)(k)其中：R(k)= =E(k)T(k)=R(k-1)+(k)(k)T(k)-R(k-1)牛顿法的优点是收敛速度快。54模型参考自适应辨识方法 考虑线性回归模型：y(k)=T(k)+e(k)，其中e(k)是零均值随机噪声。以输出估计误差为反馈信号，以PI调节器的方式来修正参数：=I+P， I(k)=I(k-1)+P(k)(k)，P(k)=Q(k)(k)其中P为对称正定矩阵，Q满足P/2+Q0 (k)为广义误差，最简单的取法为(k)= y(k)-T(k)(k-1)= 0(k)55极大似然法 输出z是一个随机变量，它的概率密度p(z|)取决于参数。当获
45、得观测序列ZL=z(1),z(2),.,z(L)T时，由该观测序列组成的联合概率密度p(ZL|)应当取得最大值。（当一个随机事件发生了，我们有理由相信，外部条件一定处于使这随机事件发生的概率最大时的状态。）那么，的极大似然估计就是使p(ZL|)|ML=max的参数估计值。由于p(ZL|)中ZL已知，因而它只是参数的函数，故称它为的似然函数。有时也记作L(ZL|)。 56极大似然法极大似然原理可用下列等价的表示方式： ， ，求解极大似然估计的下一步是要给出p(ZL|)的具体描述。 57独立观测情况设z(1),z(2),.,z(L)是一组在独立观测条件下获得的随机序列，即各观测值是互相独立的，则p
46、(ZL|)可简化为：p(ZL|)= p(z(1)|)p(z(2)|).p(z(L)|) = p(z(k)|)，其对数似然函数为：L(ZL|)= ln p(z(k)|)。58独立观测情况设z(k)N(m,2)，即：p(z(k)|)= ，负对数似然函数为：- L(ZL|)= +Lln+(L/2)ln2当已知时，准则函数就是：J()=当未知时，可先由min(J()求 及Jmin，再由min(-L(ZL|)求 。 ， 59非独立观测若z(1),z(2),.,z(L)是非独立观测条件下获得的随机序 列 ， 即 观 测 值 z(k)是 在 已 有 观 测z(1),z(2),.,z(k)的基础上得到的，则p
47、(ZL|)应按条件概率的乘法规则写成：p(ZL|)= p(z(L)|ZL-1,)p(ZL-1 |)，以此类推，有：p(ZL|)= p(z(k)|Zk-1,)。60非独立观测 设z(k)N(mk,k)，并取mk为z(k)的条件均值： (k)= Ez(k)| Zk-1，即：p(z(k)|Zk-1,)= ，负对数为：-L(ZL|)= +(L/2)ln2。当k已知时，就是：J()=当k未知时，可先由min(J()求 及Jmin，再由min(- L(ZL|)求 。 61例 考虑模型：z= + ，其中wN(0,2I),则 N( ,2I)，p( |)= ，J=lnL()=62例 i=1,2,.,n 63例直
48、接求解以上非线性方程组是困难的，如果已有参数的某个估计，并用其构成如下预滤波器： ， ， ，则上面前两个方程组可写成： ， ，这等价于辅助变量法。对于噪声模型，我们已有：定义： ，则w=Cv/D=Cv-(D-1)w，再次构造预滤波器：v*=v/D，w*=w/D,则 ， 这也等价于辅助变量法。 6465预报误差估计法 对于模型结构已知的系统，如果我们获得参数的某个估计，则可对输出进行预报，那么用预报误差的大小来衡量参数估计的优劣也是合理的。定义输出估计误差：e(k)=y(k)- (k)， 将收敛于e(k)的协方差阵。通常采用的预报误差准则有：J1()=Tr(D(), J2()=log det(D
49、()。可以证明，当e(k)N(0,e)时，若J1()中的取为=Le-1，则预报误差估计等价于极大似然估计； 66例 考虑模型：z= + 。利用已经获得的参数估计和噪声w的到k-1时刻为止的估计，我们有输出估计 ,并且可得噪声估计 =由此可见，以 为输出估计误差时，预报误差法与极大似然法的准则函数是等价的。67Bayse估计 参数是与观测值关联的随机变量。在在给给定定的的观观测测数数据据下下，使使得得验验后后的的条条件件概概率率密密度度p(|ZL)最最大大的的参参数数最最有有可可能能是真实参数。是真实参数。根据全概率公式：p(|ZL)= p(z(L)|ZL-1,)p(|ZL-1)/p(z(L)|
50、ZL-1) = p(|Z0)其中p(|ZL-1)为参数的验前概率密度函数，p(z(L)|ZL-1, )为数据的条件概率密度函数，p(z(L)|ZL-1)是与参数无关的项。根据验后条件概率密度p(|ZL)，有两种方法可求得参数估计，一是极大验后参数估计法，另一是条件期望参数估计法，它们统称为Bayes方法。 681 极大验后参数估计法 参数的极大验后估计必定满足： =0= +若令上式右边第一项为零，所得的参数估计即为极大似然估计。可见，极大验后参数估计比极大似然估计多考虑了参数的验前概率知识，因而一般情况下极大验后参数估计优于极大似然估计。但是参数的验前概率密度不易获得，因此，极大似然估计比极大
51、验后参数估计用得更普遍。 692 条件期望参数估计法 直接以参数的条件期望作为参数估计值是合乎逻辑的。 (L)=E|ZL。可以证明，条件期望参数估计等价于极小化参数估计误差的方差： minE- (L)T- (L)|ZL因此也称为最小方差估计。 70例 考虑线性回归模型：y(k)=T(k)+e(k)，其中e(k)是零均值随机噪声，e(k)N(0, )。设参数在数据Z0条件下的验前概率密度为：p(|Z0)= , 其中N=dim。那么在数据ZL-1条件下概率密度为：p(|ZL-1)=由噪声e(k)N(0, )可知，p(z(L)|ZL-1,)= ，于是有： 71log p(|ZL)=const- -=const-若令则log p(|ZL)= const-显然极大验后参数估计应该是 ，而从正态分布的形式来看，它也是参数的条件期望估计。 72
