多步最小二乘法程序msls

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 "

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)