系统辨识与建模.pptx

《系统辨识与建模.pptx》由会员分享，可在线阅读，更多相关《系统辨识与建模.pptx（72页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、会计学1系统辨识与建模系统辨识与建模第三讲第三讲 参数估计的批量法参数估计的批量法 n n最小二乘算法 n n参数估计的一般性质 n n最小二乘、加权最小二乘估计性质 n n噪声方差估计 n n广义最小二乘 n n偏倚校正算法 n n辅助变量法 n n多步最小二乘 n n相关最小二乘 第1页/共72页最小二乘算法最小二乘算法考虑差分方程：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-

2、1),(k)=-y(k-1),-y(k-n),u(k-,-y(k-n),u(k-),),u(k-m-,u(k-m-+1)T T y(k)=T(k)+w(k)（以下令（以下令=1=1）汇总的观察误差第2页/共72页y(k)=-ay(k)=-a1 1y(k-1)-y(k-1)-a-an ny(k-n)+by(k-n)+b1 1u(k-1)+u(k-1)+b+bmmu(k-m)+w(k)u(k-m)+w(k)y(k+1)=-ay(k+1)=-a1 1y(k)-y(k)-a-an ny(k-n+1)+by(k-n+1)+b1 1u(k)+u(k)+b+bmmu(k-m+1)+w(k+1)u(k-m+1

3、)+w(k+1)y(k+2)=-ay(k+2)=-a1 1y(k+1)-y(k+1)-a-an ny(k-n+2)+by(k-n+2)+b1 1u(k+1)+u(k+1)+b+bmmu(k-m+2)+w(k+2)u(k-m+2)+w(k+2)T T(k)=-y(k-1)(k)=-y(k-1)-y(k-n)u(k-1)-y(k-n)u(k-1)u(k-m)u(k-m)T T(k+1)=-y(k)(k+1)=-y(k)-y(k-n+1)u(k)-y(k-n+1)u(k)u(k-m+1)u(k-m+1)T T(k+2)=-y(k+1)(k+2)=-y(k+1)-y(k-n+2)u(k+1)-y(k-

4、n+2)u(k+1)u(k-m+2)u(k-m+2)T T(k+N)=-y(k+N-1)(k+N)=-y(k+N-1)-y(k-n+N)u(k+N-1)-y(k-n+N)u(k+N-1)u(k-m+N)u(k-m+N)第3页/共72页最小二乘算法最小二乘算法若我们的观测数据可写出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)w(k)为为白白色色时时,条条件件1 1满足满足）条件2:T可逆 第4页/共72页当当当当w(k)w(k)为白色时为白色时为白色时为白

5、色时,条件条件条件条件1 1满足满足满足满足T WN=-y(k-1)-y(k)-y(k+1)-y(k-1)-y(k)-y(k+1)-y(k+N-2)w(k)-y(k+N-2)w(k)-y(k-2)-y(k-1)-y(k)-y(k-2)-y(k-1)-y(k)-y(k+N-1)w(k+1)-y(k+N-1)w(k+1).-y(k-n)-y(k-n+1)-y(k-n+2)-y(k-n)-y(k-n+1)-y(k-n+2)-y(k-n+N-1).-y(k-n+N-1).u(k-1)u(k)u(k+1)u(k-1)u(k)u(k+1)u(k+N-2).u(k+N-2).u(k-m)u(k-m+1)u(

6、k-m+2)u(k-m)u(k-m+1)u(k-m+2)u(k-m+N-1)w(k+N-1)u(k-m+N-1)w(k+N-1)y(k)y(k)与与w(k)w(k)、w(k-1)w(k-1)、w(k-2)w(k-2)、相关。相关。相关。相关。当当当当w(k+1)w(k+1)与与w(k)w(k)、w(k-1)w(k-1)、w(k-w(k-2)2)、不不相关时，相关时，相关时，相关时，y(k)y(k)与与w(k+1)w(k+1)也不也不相关相关相关相关。第5页/共72页最小二乘算法最小二乘算法以上结果等同于求使：J=(y(k)-T(k)2最小的，因此称为最小二乘算法。=0=0T YN=T-正则方程

7、T-正则矩阵(对称、正定、可逆)第6页/共72页加权最小二乘加权最小二乘 若对各次观测数据加不同的权，即求使J=k(y(k)-T(k)2最小的,则得到参数的加权最小二乘估计：LSLS=(=(T T)-1-1T TYYN N 的对角元由k构成第7页/共72页参数估计的一般性质参数估计的一般性质 1无偏性 如果E()=0,(=-)或E()=，则称估计为无偏的。2有效性如果无偏估计 满足Cov()=M-1，则称估计为有效的。其中：称为Fisher信息矩阵，其逆M-1称为Crammer-Rao 下界。在一般情况下，有Crammer-Rao不等式：Cov()=E(-)(-)TM-1 第8页/共72页有效

8、性例有效性例对于z=H+w，若噪声w为零均值、协方差阵w=E(WWT)的正态分布噪声，即：WN(0,w)，则输出信号ZN(E(H),w)，即，因此：于是，M=E(HTw-1H)，其Crammer-Rao不等式为：Cov()E(HTw-1H)-1。有效估计也称为最小方差估计、马尔可夫估计。第9页/共72页参数估计的一般性质参数估计的一般性质 3一致性若估计 为渐近无偏的(E()=0），且 Var()=0，则称 为一致估计。Var()=E(-)T(-)第10页/共72页最小二乘、加权最小二乘最小二乘、加权最小二乘估计性质估计性质估计性质估计性质 最小二乘最小二乘 无偏性 =(T)-1T WN，当:

9、1 w(k)为零均值独立随机过程（白噪声）时；或者2 2 E(T WN)=0时，有：E()=0。无偏条件较严格。有效性 Cov()=E(T)-1TWNWTN(T)-1=(T)-1Tw(T)-1M-1,若W是服从正态分布的独立随机过程，则w=I，因而Cov()=I(T)-1=M-1。一致性 若W是服从正态分布的独立随机过程，则 Cov()=0,其中 收敛于一个有界非奇异矩阵。而Var()只是Cov()的对角线上元的和，因而也有：Var()=0。第11页/共72页最小二乘、加权最小二乘估计性最小二乘、加权最小二乘估计性质质估计性质估计性质 加权最小二乘加权最小二乘 无偏性 =(T)-1TWN无偏条

10、件同最小二乘估计。有效性 Cov()=(T)-1Tw(T)-1，当=w-1时，Cov()=(Tw-1)-1=M-1，此时不要求W是服从正态分布的独立随机过程 一致性 若W是服从正态分布的独立随机过程，则 Cov()=0,其中 收敛于一个有界非奇异矩阵。而Var()只是Cov()的对角线上元的和，因而也有：Var()=0。第12页/共72页影响最小二乘计算结果的因素影响最小二乘计算结果的因素:n nu(k)u(k)T T T T是否可逆是否可逆是否可逆是否可逆 n n噪声噪声噪声噪声w(k)w(k)大大大大,则则则则 的方差大的方差大的方差大的方差大 白色零均值白色零均值白色零均值白色零均值,是

11、无偏估计是无偏估计是无偏估计是无偏估计n n数据总量数据总量数据总量数据总量N NN N越大越大越大越大,的的的的方差越小方差越小方差越小方差越小第13页/共72页最小二乘算法的最小二乘算法的MATLAB程序程序 读入数据读入数据读入结构读入结构构造矩阵构造矩阵和和Y Y计算计算T T和和T TY Y计算计算Ls第14页/共72页function zta,m,tao=ls(tt)function zta,m,tao=ls(tt)%最小二乘法最小二乘法for MISO,ttfor MISO,tt的格式为：的格式为：第一列是系统输出数据，其它第一列是系统输出数据，其它列是对应的输入数据列是对应的输

12、入数据ll=size(tt);ll=size(tt);%得到数据维数得到数据维数r=ll(2)-1;r=ll(2)-1;m(1)=input(m(1)=input(输入输入 A(z)A(z)的阶次的阶次););%指定模型结构指定模型结构tao(1)=0;tao(1)=0;for i=1:rfor i=1:r i i%给出输入编号给出输入编号m(i+1)=input(m(i+1)=input(输入输入 B(z)B(z)的阶次的阶次)+1;)+1;%阶次加一，表示参数个数阶次加一，表示参数个数tao(i+1)=input(tao(i+1)=input(输入输入B(z)B(z)的时延的时延););e

13、ndendn=m(1)+max(tao);n=m(1)+max(tao);%算出一个方程算出一个方程最多使用的数据最多使用的数据lll=ll(1)-n;lll=ll(1)-n;%算出可列出的方程数算出可列出的方程数in1=r+1;in1=r+1;%构造观测数据矩阵构造观测数据矩阵ff ffkn=0;kn=0;for k=1:in1 for k=1:in1%每一行中的变量循环每一行中的变量循环for i=1:lll for i=1:lll%列循环列循环for j=1:m(k)for j=1:m(k)%每行变量中的观测数据每行变量中的观测数据循环循环jtao=j+tao(k);jtao=j+tao

14、(k);%构造时考虑时延构造时考虑时延if k1 ff(i,j+kn)=tt(i+n-jtao+1,k);endif k1 ff(i,j+kn)=tt(i+n-jtao+1,k);endif k=1,ff(i,j)=-tt(i+n-jtao,k);end if k=1,ff(i,j)=-tt(i+n-jtao,k);end%输输出变量变号出变量变号end;end;end;end;kn=kn+m(k);kn=kn+m(k);%算出行变量的启始位算出行变量的启始位置置end;end;for i=1:lllfor i=1:lllyy(i)=tt(i+n,1);yy(i)=tt(i+n,1);%构造输

15、出向量构造输出向量end;end;fa=ff*ff;fa=ff*ff;%最小二乘计算最小二乘计算fay=ff*yy;fay=ff*yy;zz=inv(fa);zz=inv(fa);zta=zz*fay zta=zz*fay%显示参数和结构显示参数和结构mm，taotaoLs.m第15页/共72页Ls1.mclear all;clear all;%最小二乘法最小二乘法for MISOfor MISOload y3;load y3;tt(:,1)=uyr(1,:);tt(:,1)=uyr(1,:);%读入数据，并赋给变量读入数据，并赋给变量tt tt。tt(:,2)=uyr(2,:);tt(:,2

16、)=uyr(2,:);%tt%tt的格式为：第一列是系统输出的格式为：第一列是系统输出数据，其它列是对应的输入数据数据，其它列是对应的输入数据clear uyr;clear uyr;plot(tt(:,1)plot(tt(:,1)ls(tt);ls(tt);第16页/共72页仿真例仿真例1.1.无噪声模型：数据文件无噪声模型：数据文件 y3.maty3.mat(1+1.5q(1+1.5q-1-1+0.7q+0.7q-2-2)y(k)=q)y(k)=q-2-23.2u(k)3.2u(k)辨识结果辨识结果辨识结果辨识结果（给定结构：（给定结构：（给定结构：（给定结构：m=2 1m=2 1，tao=

17、0 2tao=0 2）zta=zta=1.5000 1.5000 0.7000 0.7000 3.2000 3.2000第17页/共72页2.2.白噪声模型：数据文件白噪声模型：数据文件白噪声模型：数据文件白噪声模型：数据文件 y5.maty5.mat(1+1.5q(1+1.5q-1-1+0.7q+0.7q-2-2)y(k)=q)y(k)=q-2-23.2u(k)+3.2u(k)+w(k)w(k)辨识结果辨识结果辨识结果辨识结果（给定结构：（给定结构：（给定结构：（给定结构：m=2 1m=2 1，tao=0 2tao=0 2）zta=zta=1.5027 1.5027 0.7032 0.703

18、2 3.1935 3.1935第18页/共72页3.3.有色噪声模型：数据文件有色噪声模型：数据文件有色噪声模型：数据文件有色噪声模型：数据文件 y6.maty6.mat(1+1.5q(1+1.5q-1-1+0.7q+0.7q-2-2)y(k)=q)y(k)=q-2-23.2u(k)+3.2u(k)+辨识结果辨识结果辨识结果辨识结果（给定结构：（给定结构：（给定结构：（给定结构：m=2 1m=2 1，tao=0 2tao=0 2）zta=zta=0.97570.9757 0.1782 0.1782 3.3955 3.3955第19页/共72页4.4.白噪声模型白噪声模型白噪声模型白噪声模型b

19、b：数据文件：数据文件：数据文件：数据文件 y5y5b b.mat.maty(k)=+y(k)=+w(k)w(k)辨识结果辨识结果辨识结果辨识结果（给定结构：（给定结构：（给定结构：（给定结构：m=2 1m=2 1，tao=0 2tao=0 2）zta=zta=1.16271.1627 0.3750 0.3750 3.3907 3.3907第20页/共72页噪声方差估计噪声方差估计e=Ye=YNN-=+W-=+WNN-=+W-=+WNN-(-(T T)-1-1T TY YNN=W=WNN-(-(T T)-1-1T T W WNN=(I-(=(I-(T T)-1-1T T)W)WNN=B W=B

20、 WNN，因为因为因为因为B BT T=B,B=B,B2 2=B=B所以，所以，所以，所以，E(eE(eT Te)=E(We)=E(WNN T TB WB WNN)=E(W =E(WNN T T W WNN-W-WNN T T(T T)-1-1T T W WNN)=N-(Tr(=N-(Tr(T T)-1-1T T)=(N-dim)=(N-dim)=(y(k)-)=(y(k)-)2 2/(N-dim)/(N-dim)在有色噪声环境下，最小二乘估计是有偏的。下面的一些在有色噪声环境下，最小二乘估计是有偏的。下面的一些在有色噪声环境下，最小二乘估计是有偏的。下面的一些在有色噪声环境下，最小二乘估计是

21、有偏的。下面的一些算法对最小二乘进行改进。算法对最小二乘进行改进。算法对最小二乘进行改进。算法对最小二乘进行改进。追迹:对角线元素之和第21页/共72页广义最小二乘广义最小二乘考虑差分方程：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()未知时，我们可考虑采用迭代估计的方法去求得。第22页/共72页广义最小二乘的计算步骤广义最小二乘的计算步骤 11令

22、令令令C C0 0()=1()=1，i=0i=0下标表示迭代次数；下标表示迭代次数；下标表示迭代次数；下标表示迭代次数；0 0=1000000=100000022计算计算计算计算(k)=C(k)=Ci i()y(k)()y(k)，(k)=C(k)=Ci i()u(k)()u(k)；i=i+1;i=i+1;3 3用最小二乘估计用最小二乘估计用最小二乘估计用最小二乘估计A()(k)=B()(k)+w(k)A()(k)=B()(k)+w(k)中的参数；中的参数；中的参数；中的参数；44用用用用估估估估计计计计模模模模型型型型 、以以以以及及及及各各各各时时时时刻刻刻刻的的的的观观观观测测测测数数数数

23、据据据据，计计计计算算算算出残差出残差出残差出残差:(k)=()y(k)-()u(k):(k)=()y(k)-()u(k)5 5 计计计计算算算算 i i=22(k)(k)及及及及=i-1i-1-i i，如如如如果果果果 小小小小于于于于一一一一定定定定数数数数，则则则则结结结结束束束束辨识。否则转下一步。辨识。否则转下一步。辨识。否则转下一步。辨识。否则转下一步。6 6 对对对对于于于于噪噪噪噪声声声声模模模模型型型型C(C()(k)=(k)=w(k)w(k)，用用用用最最最最小小小小二二二二乘乘乘乘估估估估计计计计出出出出参参参参数数数数，得到更新的得到更新的得到更新的得到更新的C Ci

24、i()()后，返回后，返回后，返回后，返回2 2。以上算法的每一次循环中都要进行滤波和两次求逆。下面以上算法的每一次循环中都要进行滤波和两次求逆。下面以上算法的每一次循环中都要进行滤波和两次求逆。下面以上算法的每一次循环中都要进行滤波和两次求逆。下面的算法将在计算工作量上有所改进。的算法将在计算工作量上有所改进。的算法将在计算工作量上有所改进。的算法将在计算工作量上有所改进。第23页/共72页function zta,m,tao=ls0(tt,m,tao,ll)function zta,m,tao=ls0(tt,m,tao,ll)%最小二乘法最小二乘法for MISO,ttfor MISO,t

25、t的格式为：第一的格式为：第一列是系统输出数据，其它列是对应列是系统输出数据，其它列是对应的输入数据的输入数据%m%m为各多项式中参数个数，应与为各多项式中参数个数，应与tt tt的列的列数一致；数一致；taotao为时延；为时延；ll=size(tt);ll=size(tt);n=max(m)+max(tao);n=max(m)+max(tao);%算出一个算出一个方程最多使用的数据方程最多使用的数据lll=ll(1)-n;lll=ll(1)-n;%算出可列出的方程数算出可列出的方程数in1=ll(2);in1=ll(2);%构造观测数据矩阵构造观测数据矩阵ff ffkn=0;kn=0;fo

26、r k=1:in1 for k=1:in1%每一行中的变量循每一行中的变量循环环for i=1:lll for i=1:lll%列循环列循环for j=1:m(k)for j=1:m(k)%每行变量中的每行变量中的观测数据循环观测数据循环jtao=j+tao(k);jtao=j+tao(k);%构造时考虑时延构造时考虑时延if k1if k1,ff(i,j+kn)=tt(i+n-ff(i,j+kn)=tt(i+n-jtao+1,k);endjtao+1,k);endif k=1,ff(i,j)=-tt(i+n-if k=1,ff(i,j)=-tt(i+n-jtao,k);end jtao,k)

27、;end%输出变量输出变量变号变号end;end;end;end;kn=kn+m(k);kn=kn+m(k);%算出行变算出行变量的启始位置量的启始位置end;for i=1:lllend;for i=1:lllyy(i)=tt(i+n,1);yy(i)=tt(i+n,1);%构造输出向构造输出向量量end;end;fa=ff*ff;fa=ff*ff;%最小二乘计算最小二乘计算fay=ff*yy;fay=ff*yy;zz=inv(fa);zz=inv(fa);zta=zz*fay zta=zz*fay%显示参数和结显示参数和结构构M taoM taoLs0.m第24页/共72页clear al

28、l;clear all;%广义最小二乘法广义最小二乘法for for MISOMISOload y6;load y6;tt(:,1)=uyr(1,:);tt(:,1)=uyr(1,:);%读入数据赋给读入数据赋给tt tt。tt(:,2)=uyr(2,:);tt(:,2)=uyr(2,:);%tt%tt的格式的格式为：第一列是系统输出数据，为：第一列是系统输出数据，其它列是对应的输入数据其它列是对应的输入数据clear uy;plot(tt(:,1)clear uy;plot(tt(:,1)ll=size(tt);ll=size(tt);%得到数据维数得到数据维数in1=ll(2);r=in1

29、-1;in1=ll(2);r=in1-1;m(1)=input(m(1)=input(输入输入 A(z)A(z)的阶次的阶次););%指定模型结构指定模型结构tao(1)=0;tao(1)=0;for i=1:r for i=1:r i i%给出输入编号给出输入编号m(i+1)=input(m(i+1)=input(输入输入 B(z)B(z)的阶次的阶次)+1;)+1;%阶次加一，表示参数个数阶次加一，表示参数个数tao(i+1)=input(tao(i+1)=input(输入输入B(z)B(z)的时延的时延););EndEndtaomax=max(tao+m);taomax=max(tao+

30、m);mc=input(mc=input(输入噪声模型输入噪声模型C(z)C(z)的阶次的阶次););c=1;d=1;c=1;d=1;%初始化滤波器初始化滤波器lb=1;xci=100000;xci0=1000000;lb=1;xci=100000;xci0=1000000;while abs(xci0-xci)0.001 while abs(xci0-xci)0.001%滤波循环滤波循环xci0=xci;xci0=xci;tt1=filter(c,d,tt);tt1=filter(c,d,tt);%输入输出滤波输入输出滤波.zta,m,tao=ls0(tt1,m,tao,ll);%zta,m

31、,tao=ls0(tt1,m,tao,ll);%主最小主最小二乘二乘Gls.m第25页/共72页for k=1:taomax for k=1:taomax%计算输出估计计算输出估计 y(k)=tt(k,1);y(k)=tt(k,1);%设定初始输出设定初始输出 endendfor k=1+taomax:ll(1)for k=1+taomax:ll(1)mm=0;mm=0;for i=1:in1for i=1:in1 if i=1 if i=1 for j=1:m(i)for j=1:m(i)fb(mm+j)=-tt(k-j,1);fb(mm+j)=-tt(k-j,1);end;end;else

32、 else for j=1:m(i)for j=1:m(i)if k-tao(i)-j+10 if k-tao(i)-j+10 fb(mm+j)=tt(k-tao(i)-j+1,i);fb(mm+j)=tt(k-tao(i)-j+1,i);else else fb(mm+j)=0;fb(mm+j)=0;end end end;end;end;end;mm=mm+m(i);mm=mm+m(i);end;end;y(k)=fb*zta;y(k)=fb*zta;%用用LSLS估计参数求输出估计估计参数求输出估计endende=tt(:,1)-ye=tt(:,1)-y;%计算方程误差计算方程误差ee=

33、0;ee=0;%计算损失函数值计算损失函数值for k=1:ll(1)for k=1:ll(1)ee=ee+e(k)*e(k);ee=ee+e(k)*e(k);endendxci=ee/ll(1)xci=ee/ll(1)taoc=0;taoc=0;%估计噪声模型估计噪声模型lc=size(e);lc=size(e);cc,mc,taoc=ls0(e,mc,taoc,lc);cc,mc,taoc=ls0(e,mc,taoc,lc);c(2:mc+1)=cc c(2:mc+1)=cc%显示参数和结显示参数和结构构lb=lb+1;lb=lb+1;end end%结束滤波循环结束滤波循环plot(e)

34、plot(e)，ylabel(ylabel(errorerror)，xlabel(time)xlabel(time)第26页/共72页广义最小二乘例广义最小二乘例1.1.有色噪声模型：数据文件有色噪声模型：数据文件有色噪声模型：数据文件有色噪声模型：数据文件 y6.maty6.mat(1+1.5q(1+1.5q-1-1+0.7q+0.7q-2-2)y(k)=q)y(k)=q-2-23.2u(k)+3.2u(k)+辨识结果辨识结果辨识结果辨识结果（给定结构：（给定结构：（给定结构：（给定结构：m=2 1m=2 1，tao=0 2tao=0 2）zta=zta=1.49501.4950 0.694

35、3 0.6943 3.2075 3.2075c=1.0000 -1.7118 0.8069c=1.0000 -1.7118 0.8069迭代次数：迭代次数：迭代次数：迭代次数：6 6第27页/共72页2.2.白噪声模型白噪声模型白噪声模型白噪声模型b b：数据文件：数据文件：数据文件：数据文件 y5y5b b.mat.maty(k)=+y(k)=+w(k)w(k)辨识结果辨识结果辨识结果辨识结果（给定结构：（给定结构：（给定结构：（给定结构：m=2 1m=2 1，tao=0 2tao=0 2）zta=zta=1.50161.5016 0.6992 0.6992 3.1881 3.1881c=1

36、.0000 -1.5067 1.4236 -1.0215 0.5519 -0.1689c=1.0000 -1.5067 1.4236 -1.0215 0.5519 -0.1689迭代次数：迭代次数：迭代次数：迭代次数：7 7第28页/共72页偏倚校正算法偏倚校正算法仍仍仍仍考考考考虑虑虑虑差差差差分分分分方方方方程程程程：A(A()y(k)=B()y(k)=B()u(k)+)u(k)+w(k)w(k)，其其其其中中中中w(k)w(k)为白噪声。令为白噪声。令为白噪声。令为白噪声。令(k)=w(k)(k)=w(k)，则，则，则，则A()y(k)=B()u(k)+(k)A()y(k)=B()u(k

37、)+(k)。分别写成回归模型：。分别写成回归模型：。分别写成回归模型：。分别写成回归模型：Y=+,=C+W,Y=+,=C+W,组合起来有组合起来有组合起来有组合起来有Y=,+WY=,+W，其其其其最小二乘解为：最小二乘解为：最小二乘解为：最小二乘解为：=，利用分块矩，利用分块矩，利用分块矩，利用分块矩阵求逆公式有：阵求逆公式有：阵求逆公式有：阵求逆公式有：=(=(T T)-1-1T TY-(Y-(T T)-1-1T T =D =D-1 T-1 TMY M=I-(MY M=I-(T T)-1-1T T D=D=T TMM须要注意，须要注意，须要注意，须要注意，的求取仍然是一个迭代过程。的求取仍然

38、是一个迭代过程。的求取仍然是一个迭代过程。的求取仍然是一个迭代过程。=Ay-Bu第29页/共72页偏倚校正算法的计算步骤偏倚校正算法的计算步骤 1 1 令令C(C()=1)=1，用用最最小小二二乘乘法法求求 =LSLS=(=(T T)-1-1T TY Y，并并保保留留、=(=(T T)-1-1T T以及以及M=I-M=I-。2 2 计算计算 =Y-=Y-，并依据，并依据 =C +W=C +W构造构造 ，计算，计算D=D=T TM M 。3 3 计算计算 =D=D-1 T-1 TMYMY，并计算，并计算 =及及 =LSLS-。4 4 若参数已收敛，则结束辨识，否则转若参数已收敛，则结束辨识，否则

39、转2 2。以以上上算算法法的的一一次次循循环环中中没没有有滤滤波波，且且只只有有一一次次求求逆逆。如如果果将将第第3 3步步中中 的的计计算算改改为为：=(=(T T )-1-1 T T ，则则还还可可省省去去D D的的计计算算。（这这一一改改进进由由夏夏天天长长首首先先给给出出。）本本本本法法法法可可可可能能能能会出现收敛慢的情况，可用对会出现收敛慢的情况，可用对会出现收敛慢的情况，可用对会出现收敛慢的情况，可用对 求均值来解决求均值来解决求均值来解决求均值来解决第30页/共72页Gls2.mclear all;clear all;%偏倚校正法偏倚校正法for MISOfor MISOloa

40、d y6;load y6;tt(:,1)=uyr(1,:);tt(:,1)=uyr(1,:);%读入数据，并读入数据，并赋给变量赋给变量tt tt。tt(:,2)=uyr(2,:);tt(:,2)=uyr(2,:);%tt%tt的格式为：第的格式为：第一列是系统输出数据，其它列是一列是系统输出数据，其它列是对应的输入数据对应的输入数据clear uyr;plot(tt(:,1)clear uyr;plot(tt(:,1)ll=size(tt);ll=size(tt);%得到数据维数得到数据维数r=ll(2)-1;r=ll(2)-1;m(1)=input(m(1)=input(输入输入 A(z)

41、A(z)的阶次的阶次););%指定模型结构指定模型结构tao(1)=0;tao(1)=0;for i=1:rfor i=1:ri i%给出输入编号给出输入编号m(i+1)=input(m(i+1)=input(输入输入 B(z)B(z)的阶次的阶次)+1;)+1;%阶次加一，表示参数个数阶次加一，表示参数个数tao(i+1)=input(tao(i+1)=input(输入输入B(z)B(z)的时延的时延););endendmc=input(mc=input(输入噪声模型输入噪声模型C(z)C(z)的阶次的阶次););n=m(1)+max(tao);n=m(1)+max(tao);%算出一个方程

42、最算出一个方程最多使用的数据多使用的数据lll=ll(1)-n;lll=ll(1)-n;%算出可列出的方程数算出可列出的方程数lb=1;xci=100000;xci0=1000000;c=1;lb=1;xci=100000;xci0=1000000;c=1;in1=r+1;in1=r+1;%构造观测数据矩阵构造观测数据矩阵ff ffkn=0;kn=0;for k=1:in1 for k=1:in1%每一行中的变量循环每一行中的变量循环for i=1:lll for i=1:lll%列循环列循环for j=1:m(k)for j=1:m(k)%每行变量中的观测数据每行变量中的观测数据循环循环jt

43、ao=j+tao(k);jtao=j+tao(k);%构造时考虑时延构造时考虑时延if k1 ff(i,j+kn)=tt(i+n-jtao+1,k);endif k1 ff(i,j+kn)=tt(i+n-jtao+1,k);endif k=1,ff(i,j)=-tt(i+n-jtao,k);end if k=1,ff(i,j)=-tt(i+n-jtao,k);end%输输出变量变号出变量变号end;end;end;end;kn=kn+m(k);kn=kn+m(k);%算出行变量的启始位置算出行变量的启始位置end;end;for i=1:lll for i=1:lllyy(i)=tt(i+n,

44、1);yy(i)=tt(i+n,1);%构造输出向量构造输出向量end;end;第31页/共72页fa=ff*ff;fa=ff*ff;%最小二乘计算最小二乘计算fay=ff*yy;fay=ff*yy;zz=inv(fa);zz=inv(fa);zaa=zz*ff;zaa=zz*ff;zta=zz*fay zta=zz*fay%显示参数和结构显示参数和结构M,taoM,taoztab=zta-zta;ztab=zta-zta;ztaa=zta;ztaa=zta;while abs(xci0-xci)0.05 while abs(xci0-xci)0.05 xci0=xci;xci0=xci;e(

45、1:mc)=0;e(1:mc)=0;%计算输出误差计算输出误差e(mc+1:lll+mc)=yy-ff*zta;e(mc+1:lll+mc)=yy-ff*zta;%计算损失函数值计算损失函数值ee=0;ee=0;for k=1:lll+mcfor k=1:lll+mc ee=ee+e(k)*e(k);ee=ee+e(k)*e(k);endendxci=ee/(lll+mc)xci=ee/(lll+mc)%估计噪声模型估计噪声模型for i=1:mcfor i=1:mcfor j=1:lllfor j=1:lllfc(j,i)=-e(j+mc-i);fc(j,i)=-e(j+mc-i);end,

46、endend,endcy=e(mc+1:lll+mc);cy=e(mc+1:lll+mc);fd=fc*fc;fd=fc*fc;%最小二乘计算最小二乘计算fdcy=fc*cy;fdcy=fc*cy;zzc=inv(fd);zzc=inv(fd);c(2:mc+1)=zzc*fdcy c(2:mc+1)=zzc*fdcy%显示显示参数和结构参数和结构ztab=ztab+zaa*fc*c(2:mc+1);ztab=ztab+zaa*fc*c(2:mc+1);zta=ztaa-ztab/lbzta=ztaa-ztab/lb%zta=ztaa-zaa*fc*c(2:mc+1)%zta=ztaa-zaa

47、*fc*c(2:mc+1)lb=lb+1;lb=lb+1;end end%结束滤波循环结束滤波循环第32页/共72页偏倚校正算法例偏倚校正算法例1.1.有色噪声模型：数据文件有色噪声模型：数据文件有色噪声模型：数据文件有色噪声模型：数据文件 y6.maty6.mat(1+1.5q(1+1.5q-1-1+0.7q+0.7q-2-2)y(k)=q)y(k)=q-2-23.2u(k)+3.2u(k)+辨识结果辨识结果辨识结果辨识结果（给定结构：（给定结构：（给定结构：（给定结构：m=2 1m=2 1，tao=0 2tao=0 2）zta=zta=1.48791.4879 0.6831 0.6831

48、3.1976 3.1976c=1.0000 -1.7030 0.7983c=1.0000 -1.7030 0.7983迭代次数：迭代次数：迭代次数：迭代次数：2020第33页/共72页辅助变量法辅助变量法分分析析最最小小二二乘乘法法中中，在在Y=+WY=+W的的各各项项乘乘上上T T,然然后后利利用用T T WW的的期期望望值值为为零零得得到到参参数数的的无无偏偏估估计计。受受此此启启发发，若若在在Y=+Y=+的的各各项项乘乘上上T T，使使其其满满足足以以下下两两个个条条件件：1.1.T T 的的期期望望值值为为零零；2.2.T T可逆，则也可得到参数的无偏估计。下面讨论辅助变量的选取：可逆

49、，则也可得到参数的无偏估计。下面讨论辅助变量的选取：设模型为设模型为A()y(k)=B()u(k)+(k)A()y(k)=B()u(k)+(k)，若，若u(k)u(k)与与 (k)(k)不相关：不相关：a a 选取辅助模型选取辅助模型D()z(k)=F()u(k)D()z(k)=F()u(k)，用，用z(k)z(k)、u(k)u(k)构造构造；b b 若系统的纯时延若系统的纯时延 已知，则可用已知，则可用u(k-)u(k-)、u(k)u(k)构造构造；c c 用用u(k)u(k)、u(k)u(k)构造构造d (k)=D()w(k),d (k)=D()w(k),若若D()D()的阶次的阶次n n

50、已知已知,则可用则可用y(k-n)y(k-n)、u(k)u(k)构造构造；e e 先求出最小二乘解先求出最小二乘解LSLS=(=(T T)-1-1T TY Y，然后依据，然后依据 ()z(k)=()z(k)=()u(k)()u(k)计算出输出估计计算出输出估计z(k)z(k)，再用，再用z(k)z(k)、u(k)u(k)构造构造；第34页/共72页clear all;clear all;%辅助变量最小二乘法辅助变量最小二乘法for for MISOMISOload y6;load y6;tt(:,1)=uyr(1,:);tt(:,1)=uyr(1,:);%读入数据，并赋读入数据，并赋给变量给变