stochastic-programming-3.cpp

很好的随机规划例题

static void Simu(double x[D+1][N+1], double y[D+1][O+1])
{
   int	i, j;
   double xx[N+1], yy[O+1], solution[8];
   float read;
   FILE  *fp1;
   FILE  *fp2;

   if(((fp1=fopen("TRAIN1.dat","r"))==NULL)&&((fp2=fopen("TRAIN2.dat","r"))==NULL)) {
     mark:
     fp1=fopen("TRAIN1.dat","w");
     fp2=fopen("TRAIN2.dat","w");
     for(i=1; i<=D; i++) {
       Input_Output(xx, yy);
       if(i%SIMU_PRINT_NUMBER==0) 
		  printf("\n%d input-output Data have been generated",i);
       for(j=1; j<=N; j++) {
          x[i][j]=xx[j];
          fprintf(fp1,"%3.8f ",xx[j]);
       }
       for(j=1; j<=O; j++) {
          y[i][j]=yy[j];
          fprintf(fp1,"%3.4f ",yy[j]);
       }
	   fprintf(fp1,"\n");

//  Convert chromosome to solution.///////////////////////////////////////////
//  Here, xx[] is used as a intermediate variable.////////////////////////////
      solution[1]=xx[1];					             				    //
      solution[2]=xx[2];								                    //
      solution[3]=xx[3];                                     	    	    //
      solution[4]=sqrt(4-xx[1]);							        	    //
      solution[5]=xx[4];									                //
      solution[6]=3-xx[2]*xx[2];									        //
      solution[7]=sqrt(2-xx[3]*xx[3]-xx[4]*xx[4]);							//
// End of conversion./////////////////////////////////////////////////////////

       for(j=1; j<=7; j++) {
          fprintf(fp2,"%3.8f ",solution[j]);
       }
       for(j=1; j<=O; j++) {
          fprintf(fp2,"%3.4f ",yy[j]);
       }
	   fprintf(fp2,"\n");
     }
     fclose(fp1);
     fclose(fp2);
   }
   else {
     printf("The training data <<TRAIN.dat>> already exists. Renew it? (y/n)");
     if(getche()=='y') goto mark;
     for(i=1; i<=D; i++) {
        for(j=1; j<=N; j++) {
          fscanf(fp1,"%f",&read);
          x[i][j]=read;
        }
        for(j=1; j<=O; j++) {
          fscanf(fp1,"%f",&read);
          y[i][j]=read;
        } 
     }
   }
   fclose(fp1);
}


static void Train(double x[D+1][N+1], double y[D+1][O+1],double standard_x[D+1][N+1], double standard_y[D+1][O+1])
{
   int i, j, k;
   double E1, E2, a, xx[N+1], yy[O+1];
   float read;
   FILE *fp;

   if((fp=fopen("WEIGHT.dat","r"))==NULL) {
        mark:
        BPA(standard_x,standard_y);
        E1=0;
        E2=0;
        for(i=1; i<=D; i++) {
           for(j=1; j<=N; j++) xx[j]=x[i][j];
           NN(xx,yy);
           a=0;
           for(j=1; j<=O; j++) a += (yy[j]-y[i][j])*(yy[j]-y[i][j]);
           E1=E1+a;
           E2=E2+sqrt(a);
        }
        E1=E1/2;
        E2=E2/D;
        printf("\n\nSum-Squared Error = %3.4f",E1);
        printf("\nAverage Error = %3.4f\n\n",E2);
        for(i=1; i<=18; i++){
           k=(int)myu(1,D);
           for(j=1; j<=N; j++) xx[j]=x[k][j];
           NN(xx,yy);
           printf("Sample %2d: y=(",i);
		   for(j=1; j<=O; j++) {
			   if(j<O) printf("%3.4f,",y[k][j]);
			   else printf("%3.4f",y[k][j]);
		   }
		   printf(") Dif.=(");
           for(j=1; j<=O; j++) {
			   if(j<O) printf("%3.4f,",yy[j]-y[k][j]);
			   else printf("%3.4f",yy[j]-y[k][j]);
		   }
           printf(")\n");
        }
        printf("\nAre you satisfactory?(y/n)");
        if(getche()=='n') goto mark;
		printf("\n");

        fp=fopen("WEIGHT.dat","w");
        for(i=1; i<=H; i++) {
           for(j=0; j<=N; j++) fprintf(fp,"%3.4f ",WH[i][j]);
           fprintf(fp,"\n");
        }
        for(i=1; i<=O; i++) {
          for(j=0; j<=H; j++) fprintf(fp,"%3.4f ",WO[i][j]);
          fprintf(fp,"\n");
        }
        fclose(fp);

   }
   else {
     printf("\n\nThe weight data <<WEIGHT.dat>> already exists. Renew it? (y/n)");
     if(getche()=='y') goto mark;
     for(i=1; i<=H; i++)
        for(j=0; j<=N; j++) {
 	     fscanf(fp, "%f", &read);
        WH[i][j]=read;
     }
	  for(i=1; i<=O; i++)
		  for(j=0; j<=H; j++) {
			  fscanf(fp,"%f", &read);
           WO[i][j]=read;
        }
     fclose(fp);
   }
}

static void evaluation(int gen)
{
  double a;
  int   i, j, k, label;

  Objectives();

  if(gen == 0){
	 for(k=0; k<=M; k++) OBJECTIVE[0][k] = OBJECTIVE[1][k];
	 for(j = 1; j <= N; j++) CHROMOSOME[0][j] = CHROMOSOME[1][j];
  }

  for(i=0; i<POP_SIZE; i++){
	  label = 0;  a = OBJECTIVE[i][0];
	  for(j=i+1; j<=POP_SIZE; j++)
		 if((TYPE*a)<(TYPE*OBJECTIVE[j][0])) {
			 a = OBJECTIVE[j][0];
			 label = j;
		 }
	  if(label != 0) {
		 for(k=0; k<=M;k++) {
			 a = OBJECTIVE[i][k];
			 OBJECTIVE[i][k] = OBJECTIVE[label][k];
			 OBJECTIVE[label][k] = a;
		 }
		 for(j = 1; j <= N; j++) {
			 a = CHROMOSOME[i][j];
			 CHROMOSOME[i][j] = CHROMOSOME[label][j];
			 CHROMOSOME[label][j] = a;
		 }
	  }
  }
}

static void selection()
{
  double r, temp[POP_SIZE+1][N+1];
  int   i, j, k;

  for(i=1; i<=POP_SIZE; i++) {
	  r = myu(0, q[POP_SIZE]);
	  for(j=0; j<=POP_SIZE; j++) {
		  if(r <= q[j]) {
			  for(k=1; k<=N; k++) temp[i][k] = CHROMOSOME[j][k];
			  break;
		  }
	  }
  }

  for(i=1; i<=POP_SIZE; i++)
	 for(k=1; k<=N; k++)
		 CHROMOSOME[i][k] = temp[i][k];
}

static void crossover()
{
  int   i, j, jj, k;
  double r, x[N+1], y[N+1];

  for(i = 1; i <= POP_SIZE/2; i++) {
	 if(myu(0,1)>P_CROSSOVER) continue;
	 j = (int)myu(1,POP_SIZE);
	 jj = (int)myu(1,POP_SIZE);
	 r = myu(0,1);
	 for(k = 1; k <= N; k++) {
		 x[k] = r*CHROMOSOME[j][k]+(1-r)*CHROMOSOME[jj][k];
		 y[k] = r*CHROMOSOME[jj][k]+(1-r)*CHROMOSOME[j][k];
	 }
	 if(constraint_check(x) == 1)
		 for(k=1; k<=N; k++) CHROMOSOME[j][k] = x[k];
     if(constraint_check(y) == 1)
		 for(k=1; k<=N; k++) CHROMOSOME[jj][k] = y[k];
  }
}



static void mutation(void)
{
  int i, j;
  double y[N+1], infty, direction[N+1];
  double precision = 0.0001;
  double INFTY = 10;

  for(i=1; i<=POP_SIZE; i++) {
	  if(myu(0,1)>P_MUTATION) continue;
	  for(j=1; j<=N; j++) direction[j] = myu(-1,1);
	  infty = myu(0,INFTY);
	  while(infty>precision) {
		  for(j=1; j<=N; j++) y[j] = CHROMOSOME[i][j]+infty*direction[j];
		  if(constraint_check(y) == 1) {
			 for(j=1; j<=N; j++) CHROMOSOME[i][j] = y[j];
			 break;
		  }
		  infty = myu(0,infty);
	  }
  }
}

static void Standarlize(double x[D+1][N+1],double y[D+1][O+1],double xx[D+1][N+1],double yy[D+1][O+1])
{  
	int i,k;
	for(i=1; i<=N; i++) {
		BOUND[i][0]=x[1][i];
		BOUND[i][1]=x[1][i];
	}
	for( i=1; i<=O; i++ ) {
		BOUND[N+i][0]=y[1][i];
		BOUND[N+i][1]=y[1][i];
	}
	for( k=2; k<=D; k++ ) {
		for( i=1; i<=N; i++ )
		{   if(BOUND[i][0]>x[k][i])	
				BOUND[i][0]=x[k][i];
			if(BOUND[i][1]<x[k][i])	
				BOUND[i][1]=x[k][i];
		}
		for( i=1; i<=O; i++ )
		{   if(BOUND[N+i][0]>y[k][i])	
				BOUND[N+i][0]=y[k][i];
			if(BOUND[N+i][1]<y[k][i])	
				BOUND[N+i][1]=y[k][i];
		}
	}
	for( i=1; i<=N+O; i++ )
		BOUND[i][1]=BOUND[i][1] - BOUND[i][0];
	for( k=1; k<=D; k++ )
	{	for( i=1; i<=N; i++ )
			xx[k][i]=(x[k][i]-BOUND[i][0])/BOUND[i][1];
		for( i=1; i<=O; i++ )
			yy[k][i]=(y[k][i]-BOUND[N+i][0])/BOUND[N+i][1];
	}
}

static void NN( double x[], double y[] )
{	
	double h[H+1];
	int i,j;
	for( i=1; i<=H; i++ )
	{	h[i] = WH[i][0];
		for( j=1; j<=N; j++ )
			h[i] += WH[i][j] * (x[j]-BOUND[j][0])/BOUND[j][1];
		h[i] = tanh( h[i] );
	}
	for( i=1; i<=O; i++ )
	{	y[i] = WO[i][0];
		for( j=1; j<=H; j++)
			y[i] += WO[i][j] * h[j];
		y[i] = y[i]*BOUND[N+i][1] + BOUND[N+i][0];
	}
}

static void Init_NN(void)
{
   int i, j;
   for(i=1; i<=H; i++)
      for(j=0; j<=N; j++)
         WH[i][j]= myu(-0.3,0.3);
   for(i=1; i<=O; i++)
      for(j=0; j<=H; j++)
         WO[i][j]= myu(-0.3,0.3);
}