⭐ 欢迎来到虫虫下载站! | 📦 资源下载 📁 资源专辑 ℹ️ 关于我们
⭐ 虫虫下载站

📄 stochastic-programming-3.cpp

📁 很好的随机规划例题
💻 CPP
📖 第 1 页 / 共 2 页
字号:
// Stochastic Programming, Dependent-Chance Programming
// Written by Microsoft Visual C++
// Copyright by UTLab @ Tsinghua University
// http://orsc.edu.cn/UTLab

#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <conio.h>
#include "UTLab.h"

#define N 4     // number of input neurons (decision variables)
#define H 10   // number of hidden neurons
#define O 3     // number of output neurons
#define D 3000  // number of training data
#define M 3     // number of objectives
#define TYPE -1 // 1=max; -1=min
#define GEN 1000
#define POP_SIZE 30
#define P_MUTATION 0.2
#define P_CROSSOVER 0.3
#define SIMU_PRINT_NUMBER 100
#define GA_PRINT_NUMBER 100

static double BOUND[N+O+1][2], WH[H+1][N+1], WO[O+1][H+1];
static double OBJECTIVE[POP_SIZE+1][M+1], q[POP_SIZE+1], CHROMOSOME[POP_SIZE+1][N+1];
static void   Simu(double x[D+1][N+1], double y[D+1][O+1]);
static void   Input_Output(double x[N+1], double y[O+1]);
static void   Train(double original_x[D+1][N+1], double original_y[D+1][O+1],
			  double standard_x[D+1][N+1], double standard_y[D+1][O+1]);
static void   Standarlize(double original_x[D+1][N+1], double original_y[D+1][O+1],
			  double standard_x[D+1][N+1], double standard_y[D+1][O+1]);
static void   BPA(double standard_x[D+1][N+1], double standard_y[D+1][O+1]);
static void   Init_NN(void);
static void   NN(double x[N+1], double y[O+1]);
static void   Init_GA(void);
static void   evaluation(int gen);
static void   selection(void);
static void   crossover(void);
static void   mutation(void);
static void   Objectives(void);
static int    constraint_check(double x);

static void Input_Output(double x[N+1], double y[O+1])
{
   int i, Num[4];
   double sum[4], solution[8], xi[4];

   for(  ;   ;  ){
	   x[1]=myu(0, 4);
	   x[2]=myu(0, sqrt(3.));
	   x[3]=myu(0, sqrt(2.));
	   x[4]=myu(0, sqrt(2.));
       sum[1]=x[1]+x[2]+x[3];
	   sum[2]=sqrt(4-x[1])+x[4];
       sum[3]=3-x[2]*x[2]+sqrt(2-x[3]*x[3]-x[4]*x[4]);
	   if(sum[1]>3.4) continue;
	   if(sum[2]>2.65) continue;
	   if(sum[3]>2.) continue;
	   if(3-x[2]*x[2]<0) continue;
	   if(2-x[3]*x[3]-x[4]*x[4]<0) continue;
	   break;
   }
	   
/// Convert chromosome, x[], to solution.  ///////////////////////////
   solution[1]=x[1];            			             			//
   solution[2]=x[2];								                //
   solution[3]=x[3];                                     	    	//
   solution[4]=sqrt(4-x[1]);							    	    //
   solution[5]=x[4];									            //
   solution[6]=3-x[2]*x[2];									        //
   solution[7]=sqrt(2-x[3]*x[3]-x[4]*x[4]);							//
/// End of conversion.////////////////////////////////////////////////   

   sum[1]=solution[1]+solution[2]+solution[3];
   sum[2]=solution[4]+solution[5];
   sum[3]=solution[6]+solution[7];       
   Num[1]=0; Num[2]=0; Num[3]=0;     
   for(i=1; i<=6000; i++) {
       xi[1]=myu(3,5);
       xi[2]=myn(3.5,1);
       xi[3]=myexp(9);
	   if((sum[1]<=xi[1])&&(sum[2]<=xi[2])) Num[1]++;
       if((sum[1]<=xi[1])&&(sum[3]<=xi[3])) Num[2]++;
       if((sum[1]<=xi[1])&&(sum[2]<=xi[2])&&(sum[3]<=xi[3])) Num[3]++;
   }
   for(i=1;i<=O;i++) y[i] = Num[i]/6000.;
}


static void Objectives(void)
{
	double x[N+1], y[O+1];
	int i,j;

	for(i=1; i<=POP_SIZE; i++) {
		for(j=1;j<=N;j++) x[j] = CHROMOSOME[i][j];
		NN(x,y);
        y[1] = 0.92-y[1];
        y[2] = 0.85-y[2];
        y[3] = 0.85-y[3];
		for(j=1;j<=O;j++) if(y[j]<0) y[j]=0;
		for(j=1;j<=O;j++) OBJECTIVE[i][j] = y[j];

	}
	for(i=1; i<=POP_SIZE; i++)
	  OBJECTIVE[i][0]= 10000*OBJECTIVE[i][1]+1000*OBJECTIVE[i][2]
	                    +OBJECTIVE[i][3];
}


static int constraint_check(double x[])
{
	double sum[4];
	if( x[1]<0||x[1]>4) return 0;
	if( x[2]<0||x[2]>sqrt(3.)) return 0;
	if( x[3]<0||x[3]>sqrt(2.)) return 0;
	if( x[4]<0||x[4]>sqrt(2.)) return 0;
	if( x[3]*x[3]+x[4]*x[4]>2) return 0;
	if( x[1]+x[2]+x[3]>2) return 0;
    sum[1]=x[1]+x[2]+x[3];
	sum[2]=sqrt(4-x[1])+x[4];
	sum[3]=3-x[2]*x[2]+sqrt(2-x[3]*x[3]-x[4]*x[4]);
	if(sum[1]>3.40) return 0;
	if(sum[2]>2.65) return 0;;
	if(sum[3]>2.) return 0;;
    return 1;
}


static void Init_GA(void)
{
  double x[N+1];
  int i,j;

  for(i=1; i<=POP_SIZE; i++) {
	  mark:
	    x[1]=myu(0, 4);
	    x[2]=myu(0, sqrt(3.));
	    x[3]=myu(0, sqrt(2.));
	    x[4]=myu(0, sqrt(2.));
	  if(constraint_check(x) == 0) goto mark;
	  for(j=1; j<=N; j++) CHROMOSOME[i][j] = x[j];
  }
}



int main( void )
{
   int i, j;
   double a, x[D+1][N+1], y[D+1][O+1], standard_x[D+1][N+1], standard_y[D+1][O+1];
   double xx[8], solution[8];
   FILE *fp;

   srand(100);
   Simu(x,y);
   Standarlize(x,y,standard_x,standard_y);
   Train(x,y,standard_x,standard_y);
   q[0]=0.05; a = 0.05;
   for(i=1;i<=POP_SIZE;i++) {a=a*0.95; q[i]=q[i-1]+a;}

   fp=fopen("RESULT.dat","w");
   Init_GA();
   evaluation(0);
   for(i=1; i<=GEN; i++) {
	  selection();
	  crossover();
	  mutation();
	  evaluation(i);
	  if(i%GA_PRINT_NUMBER != 0) continue;

//  Convert chromosome to solution.///////////////////////////////////////////
//  Here, xx[] is used as an intermediate variable.////////////////////////////
	  for(j=1; j<=N; j++) xx[j]=CHROMOSOME[0][j];							//
      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(fp,"%3.4f ",solution[j]);
	  for(j=1; j<=M; j++) fprintf(fp,"%3.4f ",OBJECTIVE[0][j]);
      fprintf(fp,"%3.4f\n",OBJECTIVE[0][0]);
	  printf("\nGeneration NO.%d\n", i);
	  printf("x=(");
	  for(j=1; j<=7; j++) { 
		  if(j<7) printf("%3.4f,",solution[j]);
		  else printf("%3.4f",solution[j]);
	  }
	  	  if(M==1) printf(")\nf=%3.6f\n", OBJECTIVE[0][1]);
	  else {
	      printf(")\nf=(");
	      for(j=1; j<=M; j++) {
		     if(j<M) printf("%3.4f,", OBJECTIVE[0][j]);
		     else printf("%3.4f", OBJECTIVE[0][j]);
		  }
          printf(")  Aggregating Value=%3.4f\n",OBJECTIVE[0][0]);
	  }
   }
   printf("\nThe final result has benn written in RESULT.dat.\n\n");
   fclose(fp);
   return 1;
}

static void BPA(double x[D+1][N+1],double y[D+1][O+1])
{	
	
	double E0 = 0.04;
    double EE0 = 0.007;
	double TOLERANCE=0.0001*((N+1+O)*H+O);
	int    i,j,k,step=0;
	double ALPHA=0.05;
	double BETA=4./3;
	double ETA=0.1;
	double MU=1.;
	double DETA_WH[H+1][N+1],DETA_WO[O+1][H+1];
	double E_WO[O+1],E_WH[H+1],h[H+1],NN_y[O+1];
	double e[O+1],lambda,E,EE,TD,temp;
	
	Init_NN();
    for( i=1; i<=H; i++ )
		for( j=0; j<=N; j++ )
			DETA_WH[i][j] = 0.;
	for( i=1; i<=O; i++ )
		for( j=0; j<=H; j++ )
			DETA_WO[i][j] = 0.;

	lambda=1.;
	do
	{	step++; 
	    E=0.;
		TD=0.;
		EE=0.;
		for( k=1; k<=D; k++)
		{	/*calculate the output of the hidden neurons*/
			h[0]=1.;
			for( i=1; i<=H; i++ )
			{	h[i] = WH[i][0];
				for( j=1; j<=N; j++ )
					h[i] += WH[i][j] * x[k][j];
				h[i] = tanh(h[i]);
			}

			/*calculate the output of the outlayer neurons*/
			for( i=1; i<=O; i++ )
			{	NN_y[i] = WO[i][0];
				for( j=1; j<=H; j++)
					NN_y[i] += WO[i][j] * h[j];
			}

			/*adjust the weights of WO:*/
			for( i=1; i<=O; i++ )
			{	e[i] = y[k][i] - NN_y[i];
				E_WO[i] = lambda*e[i] + (1-lambda)*tanh(BETA*e[i]);
				for( j=0; j<=H; j++ )
				{	DETA_WO[i][j] = ETA*DETA_WO[i][j] + ALPHA*E_WO[i]*h[j];
					TD += fabs(DETA_WO[i][j]);
					WO[i][j] += DETA_WO[i][j];
				}
			}

			/*adjust the weights of WH:*/
			for( j=1; j<=H; j++ )
			{	for( i=1,temp=0.; i<=O; i++ )
					temp += E_WO[i] * WO[i][j];
				E_WH[j] = ( 1-h[j]*h[j])*temp;
				DETA_WH[j][0] = ETA*DETA_WH[j][0] +ALPHA*E_WH[j];
				TD += fabs(DETA_WH[j][0]);
				WH[j][0] += DETA_WH[j][0];
				for( i=1; i<=N; i++ )
				{   DETA_WH[j][i]= ETA*DETA_WH[j][i] +ALPHA*E_WH[j]*x[k][i];
					TD += fabs(DETA_WH[j][i]);
					WH[j][i] += DETA_WH[j][i];
				}
			}

			/*calculate the total errors*/
			for( i=1; i<=H; i++ )
			{	h[i] = WH[i][0];
				for( j=1; j<=N; j++ )
					h[i] += WH[i][j]*x[k][j];
				h[i] = tanh(h[i]);
			}
			for( i=1; i<=O; i++ )
			{	NN_y[i] = WO[i][0];
				for( j=1; j<=H; j++)
					NN_y[i] += WO[i][j] * h[j];
			}
			for( i=1; i<=O; i++ )
			{	e[i] = y[k][i] - NN_y[i];
				EE += fabs(e[i])/D;
				E += 0.5*e[i]*e[i];
			}
		}
		lambda = exp( -MU/(E*E));
		TD=TD/D;
		if(step%100==0) printf("\nBPA_%d  (Std)Sum-squared-E=%f  Average-E=%f  Gradient=%f;",step,E,EE,TD);
	} while ( step<2000 && TD>TOLERANCE && E>E0 && EE>EE0);
	printf("\nBPA_%d  (Std)Sum-squared-E=%f  Average-E=%f  Gradient=%f;",step,E,EE,TD);
}

⌨️ 快捷键说明

复制代码 Ctrl + C
搜索代码 Ctrl + F
全屏模式 F11
切换主题 Ctrl + Shift + D
显示快捷键 ?
增大字号 Ctrl + =
减小字号 Ctrl + -