多步最小二乘法程序msls ⅲ

1、13多步最小二乘法程序MSLS 用多步最小二乘法遞推算法估計如下模型的參數(shù)： 式中 為高斯白噪聲，均值為0，方差為 0.1，輸入為M序列信號，。本題采用MSLS方法III 估計，用一個擴大的差分方程作輔助模型。在這個差分方程中，當(dāng)擬合系統(tǒng)的輸入輸出數(shù)據(jù)時，殘差是不相關(guān)的，然后用最小二乘來辨識這個增廣系統(tǒng)，接著在第二級、第三級再估計原始系統(tǒng)和噪聲系統(tǒng)參數(shù)。定義兩個新的多項式和，則有：易知這個增廣系統(tǒng)（輔助模型）是5階的。第一級 先估計上面的輔助模型式，令定義參數(shù)向量為代入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，可以去掉參數(shù)向量中的該項，并相應(yīng)減少數(shù)據(jù)矩陣中對應(yīng)的一列。由輔助模型式可得該參數(shù)向量的LS估計為式中 第二級 由多項式的定義式可得其中已由第一級LS估計出來，通過上式又可估計出。將上式展開，然后令兩邊相同z冪次的系數(shù)相等，這樣就可得到個關(guān)于a和b的線性方程組。用所有的e和f的估計來代替e和f項，這些方程可寫成如下向量形式：其中，h為方程中的隨機誤差向量。即 于是系統(tǒng)參數(shù)向量的LS估計可表示為：第三級 估計噪聲多項式的系數(shù)。由多項式的定義式直接展開可得8個關(guān)于c的線性方程組。與第二級相同，令兩邊相同z冪次的系數(shù)相等，可得如下

3、向量形式：其中，為方程誤差向量有關(guān)量。即 于是系統(tǒng)參數(shù)向量的LS估計可表示為M序列作為輸入U的起始位置不同，同樣也會影響辨識精度。本題中，當(dāng)n=10時，選取白噪聲和M序列見數(shù)據(jù)文件WhiteNoise.txt和Mserials.txt，當(dāng)M序列的起始點為37時精度最高。本程序可以方便地設(shè)置不同的M序列起始位置觀察辨識效果。程序運行結(jié)果如下圖示：運行后將產(chǎn)生數(shù)據(jù)文件z_msls.txt、h_msls.txt、sita_msls.txt、c_msls.txt分別存放輸出序列、第一級的輔助模型參數(shù)辨識結(jié)果、條二級系統(tǒng)模型參數(shù)辨識結(jié)果、第三級噪聲模型參數(shù)辨識結(jié)果。源程序：#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=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