《多步最小二乘法程序msls ⅲ.doc》由会员分享，可在线阅读，更多相关《多步最小二乘法程序msls ⅲ.doc（10页珍藏版）》请在三一文库上搜索。
1、13多步最小二乘法程序MSLS 用多步最小二乘法递推算法估计如下模型的参数： 式中 为高斯白噪声，均值为0，方差为 0.1，输入为M序列信号，。本题采用MSLS方法III 估计，用一个扩大的差分方程作辅助模型。在这个差分方程中，当拟合系统的输入输出数据时，残差是不相关的，然后用最小二乘来辨识这个增广系统，接着在第二级、第三级再估计原始系统和噪声系统参数。定义两个新的多项式和，则有：易知这个增广系统（辅助模型）是5阶的。第一级 先估计上面的辅助模型式，令定义参数向量为代入A、B、C计算可得e1=1.9，e2=1.46，e3=0.539，e4=0.0815，e5=0.0082；f0=0，f1=0.
2、7，f2= 0.8，f3= 1.213，f4= 0.615。因f0=0，可以去掉参数向量中的该项，并相应减少数据矩阵中对应的一列。由辅助模型式可得该参数向量的LS估计为式中 第二级 由多项式的定义式可得其中已由第一级LS估计出来，通过上式又可估计出。将上式展开，然后令两边相同z幂次的系数相等，这样就可得到个关于a和b的线性方程组。用所有的e和f的估计来代替e和f项，这些方程可写成如下向量形式：其中，h为方程中的随机误差向量。即 于是系统参数向量的LS估计可表示为：第三级 估计噪声多项式的系数。由多项式的定义式直接展开可得8个关于c的线性方程组。与第二级相同，令两边相同z幂次的系数相等，可得如下
3、向量形式：其中，为方程误差向量有关量。即 于是系统参数向量的LS估计可表示为M序列作为输入U的起始位置不同，同样也会影响辨识精度。本题中，当n=10时，选取白噪声和M序列见数据文件WhiteNoise.txt和Mserials.txt，当M序列的起始点为37时精度最高。本程序可以方便地设置不同的M序列起始位置观察辨识效果。程序运行结果如下图示：运行后将产生数据文件z_msls.txt、h_msls.txt、sita_msls.txt、c_msls.txt分别存放输出序列、第一级的辅助模型参数辨识结果、条二级系统模型参数辨识结果、第三级噪声模型参数辨识结果。源程序：#include stdio.
4、h#include stdlib.h#include math.h#include brmul.c#include yrinv.cint main() FILE *fp1,*fp2,*fp3,*fp4; static double h511,u651,e651,z651,z16011,y651,y16001,v651,v1651,pp55,ss51; static double u160151,u251601,w51,w115,s51,s151,c21,o12,o121,p55;static double q5151,qu51601,w1p15,pw51,k51,g22,c121,gg22;
5、static double a,b,wpw1,w1s1,k1,err,ogo1,o1c1,o1g12,go21,k222,b1; int i,j,n,m; /*if(fp1=fopen(h.dat,w)=NULL) printf(ERROR); exit(1); if(fp2=fopen(M.dat,r)=NULL) printf(ERROR); exit(1);if(fp3=fopen(wnoise1.dat,r)=NULL) printf(ERROR); exit(1); if(fp4=fopen(Z.dat,W)=NULL) printf(ERROR); exit(1); */fp1=f
6、open(h1.dat,w);fp2=fopen(M.dat,r);fp3=fopen(wnoise1.dat,r);fp4=fopen(msls6.dat,w);for(i=0;i651;i+)fscanf(fp2,%lf,&ui);for(i=0;i651;i+)fscanf(fp3,%lf,&ei);v0=e0;v1=-1.0*v0+e1;for(i=2;i651;i+)vi=-1.0*vi-1-0.41*vi-2+ei;z0=v0;z1=-0.9*z0+0.7*u0+v1;z2=-0.9*z1-0.15*z0+0.7*u1-1.5*u0+v2;/for(i=0;i4;i+)/fprin
7、tf(fp4,%lfn,zi);for(i=3;i651;i+)zi=-0.9*zi-1-0.15*zi-2-0.02*zi-3+0.7*ui-1-1.5*ui-2+vi; /fprintf(fp4,%lfn,zi);for(i=0;i601;i+)z1i0=zi+50-1;/printf(%lfn,z1i0);w1s0=0.0;wpw0=0.0;for(i=0;i5;i+)pii=1.0e+6;for(i=0;i=600;i+)for(j=0;j=50;j+)u1ij=u50-j+i;for(i=0;i=50;i+)for(j=0;j=600;j+)u2ij=u1ji;brmul(u2,u1
8、,51,601,51,q);yrinv(q,51);brmul(q,u2,51,51,601,qu);brmul(qu,z1,51,601,1,h);/printf(%lf ,h00);for(i=0;i51;i+)fprintf(fp1,%lfn,hi0);fclose(fp1);fclose(fp2);fclose(fp3);for(i=0;i651;i+)a=0.0; b=0.0;if(i50)for(j=0;j=i;j+)a=a+hj0*ui-j;yi=a;elsefor(j=0;j=50;j+)b=b+hj0*ui-j;yi=b;w00=-z0;w30=u0;/w40=u0;n=0;
9、/*for(i=0;i4;i+)printf(%lf ,wi0);printf(n);*/for(m=0;m600;m+) for(i=0;i5;i+)si0=s1i0;for(i=0;i5;i+) w10i=wi0;/*for(i=0;i4;i+)printf(%lf ,w10i);*/brmul(w1,p,1,5,5,w1p);/printf(%lfn,w1p00);brmul(w1p,w,1,5,1,wpw);/printf(%lfn,wpw0);k1=1.0/(wpw0+1.0);brmul(p,w,5,5,1,pw);/printf(%lf,k1);for(i=0;i5;i+)ki0
10、=pwi0*k1;/printf(%lfn,k00);brmul(w1,s,1,5,1,w1s);b=zn+1-w1s0;for(i=0;i5;i+)s1i0=si0+ki0*b;/printf(%lf ,s1i0);/printf(n);/printf(n);/*if(m=300)getch();if(m=400)getch();if(m=500)getch();*/brmul(pw,w1p,5,1,5,pp);for(i=0;i5;i+)for(j=0;j5;j+)pij=pij-ppij/(1.0+wpw0);n=n+1;w00=-zn;w10=-zn-1;w20=-zn-2;w30=u
11、n;w40=un-1;/w50=un-2;/for(i=0;i5;i+)/ci=(si0-s1i0)/si0; /di=fabs(ci);/printf(#%lf ,ci0);/printf(&%lf ,di);/*t=d0;if(td1)t=d1;else if(td2)t=d2;else if(td3)t=d3;printf(n%lfn,t);*/for(i=0;i5;i+)printf(%lfn,s1i0); fprintf(fp4,%lf ,s1i0);fprintf(fp4,n);ss00=0.9;ss10=0.15;ss20=0.02;/ss30=0.0;ss30=0.7;ss40
12、=-1.5;err=0.0;for(i=0;i5;i+)err=err+(s1i0-ssi0)*(s1i0-ssi0);printf(误差平方和:%lfn,err);o1c0=0.0;ogo0=0.0;n=0;for(i=0;i2;i+)gii=1.0e+6;v10=z0+0.0*u0;v11=z1+s100*z0-0.0*u1-s130*u0;v12=z2+s100*z1+s110*z0-0.0*u2-s130*u1-s140*u0;/v13=z3+s00*z2+s10*z1+s20*z0-s30*u2-s40*u1;for(i=3;i651;i+)v1i=zi+s100*zi-1+s110
13、*zi-2+s120*zi-3-0.0*ui-s130*ui-1-s140*u0;for(i=0;i651;i+)v1i=v1i;o00=v10;o01=0;for(m=0;m600;m+) for(i=0;i2;i+)ci0=c1i0;for(i=0;i2;i+) o1i0=o0i;/*for(i=0;i4;i+)printf(%lf ,w10i);*/brmul(o1,g,1,2,2,o1g);/printf(%lfn,w1p00);brmul(o1g,o,1,2,1,ogo);/printf(%lfn,wpw0);k1=1.0/(ogo0+1.0);brmul(g,o,2,2,1,go)
14、;/printf(%lf,k1);for(i=0;i2;i+)k2i0=goi0*k1;/printf(%lfn,k00);brmul(o1,c,1,2,1,o1c);b1=vn+1-o1c0;for(i=0;i2;i+)c1i0=ci0+k2i0*b1;/printf(%lf ,c1i0);/printf(n);/printf(n);/*if(m=300)getch();if(m=400)getch();if(m=500)getch();*/brmul(go,o1g,2,1,2,gg);for(i=0;i2;i+)for(j=0;j2;j+)gij=gij-ggij/(1.0+ogo0);n
15、=n+1;o00=-vn;o01=-vn-1;/for(i=0;i5;i+)/ci=(si0-s1i0)/si0; /di=fabs(ci);/printf(#%lf ,ci0);/printf(&%lf ,di);/*t=d0;if(td1)t=d1;else if(td2)t=d2;else if(td3)t=d3;printf(n%lfn,t);*/for(i=0;i2;i+)printf(%lfn,c1i0);fprintf(fp4,%lf ,c1i0);fclose(fp4);err=0.0;err=(c100-1.0)*(c100-1.0)+(c110-0.41)*(c110-0.41);printf(误差平方和为：%lfn,err);return 0;