ywls.c

遗忘因子法递推算法用于估计数学模型的参数,精确度较高!

//////////////////////////////////////////////
// 用遗忘因子法递推法估计参数：				//	
// 数据读入：m序列和白噪声分别由m.txt文件   //
// 和whitenoise.txt文件中读入，             //
// 程序执行结果存入"递推结果.txt"中。       //
//////////////////////////////////////////////
#include <stdio.h>
#include <math.h>
#include "brmul.c"
#define L 1000
#define N 900
#define a 1000
void main()
{ 
	int i,j,m;
	double s=0.99;
	double d1[2]={0.8,0.5},d2[2]={0.6,0.3};
	double u[L],e[L],y[L];
	double B[2][N],k[2][1],p[2][2];
	double w[2][1],wt[1][2],kwt[2][2];
	double wtb,b[2],pw[2][1],wtp[1][2],wtpw,pwwtp[2][2];
	FILE *fp1,*fp2,*fp3;
	fp1=fopen("m.txt","r");
	fp2=fopen("whitenoise.txt","r");
	fp3=fopen("递推结果.txt","w");
	for(i=0;i<L;i++)
	{ 
		fscanf(fp1,"%lf",&u[i]);//读入输入
		fscanf(fp2,"%lf",&e[i]);//读入白噪声
	}
	fclose(fp1);
	fclose(fp2);
	y[0]=e[0];
	for(i=1;i<L;i++)
	{  
		if (i<300)
			y[i]=-0.8*y[i-1]+0.5*u[i-1]+e[i];//求出输出
		else 
			y[i]=-0.6*y[i-1]+0.3*u[i-1]+e[i];
	}
	for(i=0;i<2;i++)
		for(j=0;j<2;j++)
		{ 
			if(i==j)
				p[i][j]=pow(a,2);//给P赋初值
			else 
				p[i][j]=0;
		}
	for(i=0;i<2;i++)//给b赋初值
		b[i]=0;
	for(m=0;m<N;m++)//循环计算
	{
		wtpw=0;
   		w[0][0]=-y[1+m];
		wt[0][0]=w[0][0];
		w[1][0]=u[1+m];
		wt[0][1]=w[1][0];
		brmul(wt,p,1,2,2,wtp);
		for(i=0;i<2;i++)
			wtpw=wtpw+wtp[0][i]*w[i][0];
		brmul(p,w,2,2,1,pw);
		for(i=0;i<2;i++)
			k[i][0]=pw[i][0]/(s+wtpw);
		brmul(k,wt,2,1,2,kwt);
		brmul(pw,wtp,2,1,2,pwwtp);
		for(i=0;i<2;i++)
			for(j=0;j<2;j++)
			{  
				p[i][j]=p[i][j]-pwwtp[i][j]/(s+wtpw);
			   	p[i][j]=p[i][j]/s;
			}
		wtb=0;
		for(i=0;i<2;i++)
			wtb=wtb+wt[0][i]*b[i];
		for(i=0;i<2;i++)
		b[i]=k[i][0]*(y[m+2]-wtb)+b[i];
		for(i=0;i<2;i++)
			B[i][m]=b[i];
	}
	for(i=0;i<2;i++)//最后结果
		printf("%lf   ",b[i]);
	printf("\n");
	printf("误差:\n");
	for(i=0;i<2;i++)
		printf("%lf   ",b[i]-d2[i]);
	for(i=0;i<N;i++)
	{ 
		for(j=0;j<2;j++)
			fprintf(fp3,"%lf    ",B[j][i]);
		fprintf(fp3,"\n");  
	}
	fclose(fp3);
}