main.cpp

自己做的小波与神经网络组合方法实现负荷预测程序

#include <iostream.h>
#include <iomanip.h>
#include <stdio.h>
#include <stdlib.h>
#include "math.h"
#include  <afx.h>
#include "time.h"
#include "funshow.h"


////////////////////////////////
///////////////c3///////////////
////////////////////////////////
#define N_c 5//学习样本个数
#define IN_c 7 //输入层神经元数目
#define HN_c 8 //隐层神经元数目
#define ON_c 1 //输出层神经元数目
double P_c[IN_c]; //单个样本输入数据
double T_c[ON_c]; //单个样本教师数据
double W_c[HN_c][IN_c]; //输入层至隐层权值
double V_c[ON_c][HN_c]; //隐层至输出层权值
double X_c[HN_c]; //隐层的输入
double Y_c[ON_c]; //输出层的输入
double H_c[HN_c]; //隐层的输出
double O_c[ON_c]; //输出层的输出
double sita_c[HN_c]; //隐层的阈值
double gama_c[ON_c]; //输出层的阈值
double err_m_c[N_c]; //第m个样本的总误差
double err_O_c[N_c][ON_c];  ///m样本第ON输出的误差
double s2_c[ON_c];   //定义敏感度s2
double s1_c[HN_c];  //定义敏感度数组s1
const  R_c=HN_c*(IN_c+1)+ON_c*(HN_c+1);  //列的常数定义
const  R1_c=R_c*N_c;
const  R2_c=R_c*R_c;
double yacobi_c[N_c][R_c];   //N*ON,IN*HN+HN+HN*ON+ON   行，列
double niu_c,beta_c=10.0;   //定义参数
double det_x_c[R_c];      ////定义误差数组
//定义一个放学习样本的结构
struct {
double input[IN_c];
double teach[ON_c];
}Study_Data_c[N_c];


///////////////////////////
//初始化权、阈值子程序/////
///////////////////////////
initial_c()
{
//隐层权、阈值初始化//
	srand( (unsigned)time( NULL ) );

	for(int i=0;i<HN_c;i++)
	{
		for(int j=0;j<IN_c;j++)
			W_c[i][j]= (double)((rand()/32767.0)*2-1); //初始化输入层到隐层的权值，-1到1之间的随机值	
	}
	for(int ii=0;ii<ON_c;ii++)
	{
		for(int jj=0;jj<HN_c;jj++)
			V_c[ii][jj]= (double)((rand()/32767.0)*2-1); //初始化隐层到输出层的权值，-1到1之间的随机值	
	}
	for(int k=0;k<HN_c;k++)
	{
		sita_c[k] = (double)((rand()/32767.0)*2-1);  //隐层阈值初始化 -1到1之间的随机值	
	}
	for(int kk=0;kk<ON_c;kk++)
	{
		gama_c[kk] = (double)((rand()/32767.0)*2-1); //输出层阈值初始化 ,-1到1之间的随机值	
	}
  return 1;
}//子程序initial()结束

////////////////////////////////
////第m个学习样本输入子程序///
///////////////////////////////
input_P_c(int m)
{
  int i;
  for(i=0;i<IN_c;i++)
  P_c[i]=Study_Data_c[m].input[i];
return 1;
}//子程序input_P(m）结束
/////////////////////////////
////第m个样本教师信号子程序//
/////////////////////////////
input_T_c(int m)
{
 int k;
 for(k=0;k<ON_c;k++)
 T_c[k]=Study_Data_c[m].teach[k];
 return 1;
}//子程序input_T(m)结束

/////////////////////////////////
//隐层各单元输入、输出值子程序///
/////////////////////////////////
H_I_O_c()
{
double sigma;
int i,j;
for (j=0;j<HN_c;j++)
{ 
 sigma=0.0;
 for (i=0;i<IN_c;i++)
  sigma+=W_c[j][i]*P_c[i];//求隐层内积
 X_c[j]=sigma+sita_c[i];//求隐层净输入
 H_c[j]=1.0/(1.0+exp(-X_c[j]));//求隐层输出
}
return 1;
}//子程序H_I_O()结束
///////////////////////////////////
//输出层各单元输入、输出值子程序///
///////////////////////////////////
O_I_O_c()
{
 int k,j;
 double sigma;
  for(k=0;k<ON_c;k++)
  {
	sigma=0.0;
    for (j=0;j<HN_c;j++)
      sigma+=V_c[k][j]*H_c[j];//求输出层内积
    Y_c[k]=sigma+gama_c[k]; //求输出层净输入
    O_c[k]=1.0/(1.0+exp(-Y_c[k]));
  }
return 1;
}//子程序O_I_O()结束

////////////////////////////////////
/////求输出误差、各样本平方误差/////
////////////////////////////////////  
  Err_yangben_c(int m)
   {    
    double abs_err[ON_c];  //每个样本的绝对误差都是从0开始的
    double sqr_err=0;   //每个样本的平方误差计算都是从0开始的
    for (int k=0;k<ON_c;k++)
	{    
	  err_O_c[m][k]=T_c[k]-O_c[k];
	  abs_err[k]=T_c[k]-O_c[k];	  
      sqr_err+=(abs_err[k])*(abs_err[k]);//求第m个样本下输出层的平方误差
    } 
    err_m_c[m]=sqr_err;   
   return 1;
  }//子程序Err_yangben(m)结束5

  /////////////////////////////////
  //N个样本的全局误差计算子程序////
  /////////////////////////////////
  double Err_Sum_c()
  {
   double total_err=0;
   int m;
   for (m=0;m<N_c;m++) 
   total_err+=err_m_c[m];//每个样本的均方误差加起来就成了全局误差 
   return total_err;
  }//子程序Err_sum()结束

////////////////////////////////////
//////////////求LM敏感度////////////
////////////////////////////////////

Err_O_H_s_c()
{
  int k;
  for (k=0;k<ON_c;k++)
  {  
   s2_c[k]=-O_c[k]*(1.0-O_c[k]);   //敏感度s2   
  }
 return 1;
}//子程序Err_O_H_s(m)结束

Err_H_I_s_c()
{
   int j,k;
   double sigma;

   
   for (j=0;j<ON_c;j++) 
   { 
	 sigma=0.0;     
     for (k=0;k<HN_c;k++) 
	 {  
	   sigma=(H_c[k]*(1.0-H_c[k]))*V_c[j][k];   
	   s1_c[k]=sigma*s2_c[j];             //敏感度s1     
     }
   }

   return 1;
}//子程序Err_H_I_s_c()结束

////////////////////////////////////
/////////////求雅可比矩阵///////////
////////////////////////////////////
qiuyacobi_c(int m)         //求雅可比矩阵各元素
{
  int i,j,jj;
    jj=0;
    
    for(i=0;i<HN_c;i++)                        
	 {
      for(j=0;j<IN_c;j++)
		{ 
		  yacobi_c[m][jj]=s1_c[i]*P_c[j];
	      jj=jj+1;
		}
      yacobi_c[m][jj]=s1_c[i];			
	  jj=jj+1;
     }

    for(j=0;j<ON_c;j++)
	 {
	  for(i=0;i<HN_c;i++)
	    {
		 yacobi_c[m][jj]=s2_c[j]*H_c[i];
         jj=jj+1;
	    }
      yacobi_c[m][jj]=s2_c[j];                    
	  jj=jj+1;
     } 
  return 1;
 } 

     
        //////////////////////////////////////////  
        ///////////////////////矩阵求逆///////////
        //////////////////////////////////////////  
 qiu_detx_c()
 {  
  int i,j,k;
  double tran[R_c][N_c];    //转置矩阵   
  double tran1[R1_c],yacobi1[R1_c];
  double jz1[R2_c],jz_ji[R_c][R_c];  
  double I[R_c][R_c]; //单位矩阵  
  double jz_he[R_c][R_c];
  double a[R_c][R_c];       //逆矩阵
  double jz2[R2_c];
  double mulmatrix[R_c][N_c];   ////a[R][R]*tran[R][N]
  double a1[R2_c];
  double jz[R1_c];
  double v_x[N_c][1],det_x1[R_c][1];   
  double jz3[R1_c],jz4[N_c];
  double det_xx[R_c];	
  int m;
   
                ///////////////////求矩阵转置tran[R][N]
  for(i=0;i<N_c;i++)
  {
    for(j=0;j<R_c;j++) 
    {
	 tran[j][i]=yacobi_c[i][j];
	}
  }
  
                 ////////////实矩阵相乘tran[R][N]*yacobi[N][R]
  
  
         ////////////二维化一维，用于矩阵相乘
  k=0;
  for(i=0;i<R_c;i++)     
  {
	for(j=0;j<N_c;j++)
	{	
	  tran1[k]=tran[i][j];
	  k=k+1;
    }
  }
  k=0;
  for(i=0;i<N_c;i++)
  {
	for(j=0;j<R_c;j++)
	{	
	  yacobi1[k]=yacobi_c[i][j];
	  k=k+1;
    }
  }  
    
   
                                  ///////矩阵相乘
  trmul(tran1,yacobi1,R_c,N_c,R_c,jz1);

   
                   ////////////一维化二维，还原
   
  for(i=0;i<R_c;i++)
  {
   for(j=0;j<R_c;j++)
     jz_ji[i][j]=jz1[i*R_c+j];
  }

             
                 ////////////求正定矩阵I[R][R]
  for(i=0;i<R_c;i++)    
  {
   for(j=0;j<R_c;j++)
     {
	   if(i==j)
	      I[i][j]=1;
       else
		   I[i][j]=0;  
      }
  }

  
  for(i=0;i<R_c;i++)
   for(j=0;j<R_c;j++)
    I[i][j]=niu_c*I[i][j];

                        ///////实矩阵相加
  for(i=0;i<R_c;i++)
   for(j=0;j<R_c;j++)
    jz_he[i][j]=jz_ji[i][j]+I[i][j];

                         //////////求逆
  
  for (i=0;i<R_c;i++)
    for (j=0;j<R_c;j++)
      a[i][j]=jz_he[i][j];



         ////////////二维化一维，用于对称正定矩阵求逆
  
  k=0;
  for(i=0;i<R_c;i++)     
  {
	for(j=0;j<R_c;j++)
	{	
	  jz2[k]=a[i][j];
	  k=k+1;
    }
  }
 
    
                     //////////////矩阵求逆
  i=ssgj(jz2,R_c);

   
                   ////////////一维化二维，还原
   
  for(i=0;i<R_c;i++)
  {
   for(j=0;j<R_c;j++)
     a[i][j]=jz2[i*R_c+j];
  }
               
           /////////////矩阵乘积//////////////    
            ////////////二维化一维，用于矩阵相乘
  k=0;
  for(i=0;i<R_c;i++)     
  {
	for(j=0;j<R_c;j++)
	{	
	  a1[k]=a[i][j];
	  k=k+1;
    }
  }
  

   
                                  ///////矩阵相乘
  trmul(a1,tran1,R_c,R_c,N_c,jz);

   
                   ////////////一维化二维，还原
   
  for(i=0;i<R_c;i++)
  {
   for(j=0;j<N_c;j++)
     mulmatrix[i][j]=jz[i*N_c+j];
  }
  
                                           //求det_x            (3)
  
  for(m=0;m<N_c;m++)
   {
	  for(i=0;i<1;i++)               //ON即输出为1 ,求v(x)
	  {   
		  v_x[m][i]=err_O_c[m][i];           /////////////误差T[k]-O[k];
	  }
  }



           ////////////二维化一维，用于矩阵相乘
  k=0;
  for(i=0;i<R_c;i++)     
  {
	for(j=0;j<N_c;j++)
	{	
	  jz3[k]=mulmatrix[i][j];
	  k=k+1;
    }
  }
  k=0;
  for(i=0;i<N_c;i++)
  {
	for(j=0;j<1;j++)
	{	
	  jz4[k]=v_x[i][j];             
	  k=k+1;
    }
  }  
    
                                  ///////矩阵相乘
  
  trmul(jz3,jz4,R_c,N_c,1,det_xx);

   
                   ////////////一维化二维，还原
   
  for(i=0;i<R_c;i++)
  {
   for(j=0;j<1;j++)
     det_x1[i][j]=det_xx[i+j];
  }
	    
  for(i=0;i<R_c;i++)
	for(j=0;j<1;j++)
      det_x_c[i]=-det_x1[i][j];    ////得到 det_x
 
  return 1; 
 }




////////////////////////////////////得到det_x,下一步修正权、阀值  
   ////////////////////////////////////
   ///////////////修正权阀值///////////
   ///////////////////////////////////
 JZ_quanfazhi_c()
 {   
   int i,j,jj;
   jj=0;
    for(i=0;i<HN_c;i++)                       
	 {
      for(j=0;j<IN_c;j++)
		{ 
		  W_c[i][j]=W_c[i][j]+det_x_c[jj];    ///W[HN][IN]; 输入层至隐层权值
	      jj=jj+1;
		}
      sita_c[i]=sita_c[i]+det_x_c[jj];		//sita[HN]; 隐层的阈值
	  jj=jj+1;
     }
    for(j=0;j<ON_c;j++)
	{
	 for(i=0;i<HN_c;i++)
	  {
		V_c[j][i]=V_c[j][i]+det_x_c[jj];     //V[ON][HN]; 隐层至输出层权值
        jj=jj+1;
	  }
      gama_c[j]=gama_c[j]+det_x_c[jj];       //gama[ON]; 输出层的阈值               
	  jj=jj+1;
    } 
  return 1;
 }


////////////////////////////
////////c3神经网络预测//////
////////////////////////////
void netc3()
{
FILE *fp,*in;
double Pre_error;  //预定误差
int study,m,k,i,j;
double sum_err,sum_err1;   //定义误差平方和
niu_c=0.01;
if ((fp=fopen("c3bp.txt","w+"))==NULL)
{
printf("不能创建c3bp.txt文件!\n");
exit(1);
}
if ((in=fopen("c3in.txt","a+"))==NULL)
{
printf("不能打开c3in.txt文件!\n");
exit(1);
}

fscanf(in,"%lf",&Pre_error);
study=0;
for (m=0;m<N_c;m++)
{ 
   for (i=0;i<IN_c;i++)
	   fscanf(in,"%lf",&Study_Data_c[m].input[i]);
}
 
for (m=0;m<N_c;m++)
{
	for (k=0;k<ON_c;k++)
	{
	  fscanf(in,"%lf",&Study_Data_c[m].teach[k]);     
	} 
}

initial_c(); //隐层、输出层权、阈值初始化 
	
 
 for (m=0;m<N_c;m++) 
 {
  input_P_c(m); //输入第m个学习样本 
  input_T_c(m);//输入第m个样本的教师信号 
  H_I_O_c(); //第m个学习样本隐层各单元输入、输出值 
  O_I_O_c(); //第m个学习样本输出层各单元输入、输出值   
  Err_yangben_c(m);    ///第m个样本的误差平方和                
  Err_O_H_s_c();    ///求敏感度s2
  Err_H_I_s_c();     ///求敏感度s1
  qiuyacobi_c(m);     ////求雅可比矩阵                    (2)
 }
 sum_err1=Err_Sum_c();       //////所有样本误差平方和         (1)	     
    
 qiu_detx_c();              //求det_x            (3)   
  
 JZ_quanfazhi_c();  ///修正权、阀值          (4)    
      
  
   /////////////////////////////////////////////
   ////////////求相应的网络输出及相应误差/////// (4)
   /////////////////////////////////////////////
 
 for (m=0;m<N_c;m++) 
 {  
  H_I_O_c(); //第m个学习样本隐层各单元输入、输出值 
  O_I_O_c(); //第m个学习样本输出层各单元输入、输出值 
  Err_yangben_c(m);   ///第m个样本的误差平方和
 }
 sum_err=Err_Sum_c();    //修正权、阀值后所有样本误差平方和

   
 while (sum_err>Pre_error)
 {   
     if(sum_err<sum_err1)
	   {         
	     niu_c=niu_c/beta_c;              ///参数修正
		 
         for (m=0;m<N_c;m++) 
		  {  			             
             H_I_O_c(); //第m个学习样本隐层各单元输入、输出值 
             O_I_O_c(); //第m个学习样本输出层各单元输入、输出值 
             Err_yangben_c(m);   ///第m个样本的误差平方和
             Err_O_H_s_c();    ///求敏感度s2
             Err_H_I_s_c();     ///求敏感度s1            
             qiuyacobi_c(m);     ////求雅可比矩阵                    (2)
		  }
          sum_err1=Err_Sum_c();          //(1)
                                                       //求det_x            (3)
          qiu_detx_c();              //求det_x            (3) 
          
                
         JZ_quanfazhi_c();  ///修正权、阀值              (4)
	     for (m=0;m<N_c;m++) 
		  {            
            H_I_O_c(); //第m个学习样本隐层各单元输入、输出值 
            O_I_O_c(); //第m个学习样本输出层各单元输入、输出值