多步最小二乘法程序msls ⅲ.docVIP

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.7，f2= – 0.8，f3= – 1.213，f4= – 0.615。 因f0=0，可以去掉参数向量中的该项，并相应减少数据矩阵中对应的一列。 由辅助模型式可得该参数向量的LS估计为 式中 第二级 由多项式的定义式可得 其中已由第一级LS估计出来，通过上式又可估计出。 将上式展开，然后令两边相同z幂次的系数相等，这样就可得到７个关于a和b的线性方程组。用所有的e和f的估计来代替e和f项，这些方程可写成如下向量形式： 其中，?为方程中的随机误差向量。 即 于是系统参数向量的LS估计可表示为： 第三级 估计噪声多项式的系数。 由多项式的定义式直接展开可得8个关于c的线性方程组。与第二级相同，令两边相同z幂次的系数相等，可得如下向量形式： 其中，ξ为方程误差向量有关量。 即 于是系统参数向量的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.h #include stdlib.h #include math.h #include brmul.c #include yrinv.c int main() { FILE *fp1,*fp2,*fp3,*fp4; static double h[51][1],u[651],e[651],z[651],z1[601][1],y[651],y1[600][1],v[651],v1[651],pp[5][5],ss[5][1]; static double u1[601][51],u2[51][601],w[5][1],w1[1][5],s[5][1],s1[5][1],c[2][1],o[1][2],o1[2][1],p[5][5]; static double q[51][51],qu[51][601],w1p[1][5],pw[5][1],k[5][1],g[2][2],c1[2][1],gg[2][2]; static double a,b,wpw[1],w1s[1],k1,err,ogo[1],o1c[1],o1g[1][2],go[2][1],k2[2][2],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=fopen(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,u[i]); for(i=0;i651;i++) fscanf(fp3,%lf,e[i]); v[0]=e[0]; v[1]=-1.0*v[0]+e[1]; for(i=2;i651