多步最小二乘法程序msls

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

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

4、stdio.h"#include "stdlib.h"#include "math.h"#include "brmul.c"#include "yrinv.c"int 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,p5

5、5;static double q5151,qu51601,w1p15,pw51,k51,g22,c121,gg22; 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&q

6、uot;); 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=fopen("h1.dat","w");fp2=fopen("m.dat","r");fp3=fopen(&quo

7、t;wnoise1.dat","r");fp4=fopen("msls6.dat","w");for(i=0;i<651;i+)fscanf(fp2,"%lf",&ui);for(i=0;i<651;i+)fscanf(fp3,"%lf",&ei);v0=e0;v1=-1.0*v0+e1;for(i=2;i<651;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.

8、15*z0+0.7*u1-1.5*u0+v2;/for(i=0;i<4;i+)/fprintf(fp4,"%lfn",zi);for(i=3;i<651;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;i<601;i+)z1i0=zi+50-1;/printf("%lfn",z1i0);w1s0=0.0;wpw0=0.0;for(i=0;i<5;i+)pii=1.0e+6;for(

9、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,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;i<51;i+)fprintf(fp1,"%lfn",hi0);fclose(fp1);fclose(fp2);fclose

10、(fp3);for(i=0;i<651;i+)a=0.0; b=0.0;if(i<50)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;/*for(i=0;i<4;i+)printf("%lf ",wi0);printf("n");*/for(m=0;m<600;m+) for(i=0;i<5;i+)si0=s1i0;for(i=0;i<5;i+) w10i

11、=wi0;/*for(i=0;i<4;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;i<5;i+)ki0=pwi0*k1;/printf("%lfn",k00);brmul(w1,

12、s,1,5,1,w1s);b=zn+1-w1s0;for(i=0;i<5;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;i<5;i+)for(j=0;j<5;j+)pij=pij-ppij/(1.0+wpw0);n=n+1;w00=-zn;w10=-zn-1;w20

13、=-zn-2;w30=un;w40=un-1;/w50=un-2;/for(i=0;i<5;i+)/ci=(si0-s1i0)/si0; /di=fabs(ci);/printf("#%lf ",ci0);/printf("&&%lf ",di);/*t=d0;if(t<d1)t=d1;else if(t<d2)t=d2;else if(t<d3)t=d3;printf("n%lfn",t);*/for(i=0;i<5;i+)printf("%lfn",s1i0); f

14、printf(fp4,"%lf ",s1i0);fprintf(fp4,"n");ss00=0.9;ss10=0.15;ss20=0.02;/ss30=0.0;ss30=0.7;ss40=-1.5;err=0.0;for(i=0;i<5;i+)err=err+(s1i0-ssi0)*(s1i0-ssi0);printf("誤差平方和:%lfn",err);o1c0=0.0;ogo0=0.0;n=0;for(i=0;i<2;i+)gii=1.0e+6;v10=z0+0.0*u0;v11=z1+s100*z0-0.0*u1-s

15、130*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;i<651;i+)v1i=zi+s100*zi-1+s110*zi-2+s120*zi-3-0.0*ui-s130*ui-1-s140*u0;for(i=0;i<651;i+)v1i=v1i;o00=v10;o01=0;for(m=0;m<600;m+) for(i=0;i<2;i+)ci0=c1i0;for(i=0;i<2;i+) o1i0=o0i;/*fo

16、r(i=0;i<4;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);/printf("%lf",k1);for(i=0;i<2;i+)k2i0=goi0*k1;/printf("%lfn",k00);brmul(o1,c,1,2,1,

17、o1c);b1=vn+1-o1c0;for(i=0;i<2;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;i<2;i+)for(j=0;j<2;j+)gij=gij-ggij/(1.0+ogo0);n=n+1;o00=-vn;o01=-vn-1;/for(i=0;i<5;i+)/ci=(si0-s1i0)/si0; /di=fabs(ci);/printf("#%lf ",ci0);/printf("&&%lf ",di);/*t=d0;if(t<d1)