myfunction.c

用C编写的部落寻优的优化算法

    
 total.f[0]=f1;   
total.f[1]=f2;
   break;

     case 38:   // Multiobjective F4
 // x1 in (0,1]
 //   x2 in (0,1]

    f1=pos.p.x[ 0];
    x2=pos.p.x[ 1];
    x=1+9*x2;
    x3=f1/x;
    f2=x*(1-pow(x3,0.25)-pow(x3,4));

 total.f[0]=f1;
total.f[1]=f2;
   break;

         case 39:   // Multiobjective F5
 // x1 in (0,1]
 //   x2 in (0,1]

    f1=pos.p.x[ 0];
    x2=pos.p.x[ 1];
    x=1+9*x2;
    x3=f1/x;
    f2=x*(1-pow(x3,0.5)-x3*sin(5*two_pi*f1));

 total.f[0]=f1;
total.f[1]=f2; 
   break;

   case 40:   // Multiobjective F6 (Deb)
 // x1 in (0,1]
 //   x2 in (0,1]

    f1=pos.p.x[ 0];
    x2=pos.p.x[ 1];  // For 2D
 //x2=0;  For 1D
    x=1+10*x2;
    x3=f1/x;

    f2=x*(1-x3*x3-x3*sin(4*two_pi*f1));

 total.f[0]=f1;
total.f[1]=f2;
   break;
   
case 41:
// Multiobjective   Coello F3
    f1=pos.p.x[ 0];
    x2=pos.p.x[ 1];

    x=11+x2*x2-10*cos(two_pi*x2);
    if (x>=f1)
    {
      f2=x*(1-sqrt(f1/x));
    }
    else
       f2=0;


 total.f[0]=f1;
total.f[1]=f2;
 break;
      
 case 42: // For implicit function
    problem[1].target=0;
    problem[1].eps=0.1;
    problem[1].printlevel=0;
    problem[1].funct=42;
	problem[1].Max_Eval=problem[0].Max_Eval;
    coeff=pos;
   post=PSO(1,problem[1].Max_Eval);
   total.f[0]=post.f.f[0];


 break;

 case 43: // Part 2 for implicit function . Here a spherical one
// coeff is a global variable

    DD= coeff.p.size;
    total.f[0]=0;
  for (d=0;d<DD;d++)
  {
     x=coeff.p.x[d];
     total.f[0]=total.f[0]+x*x;
  }
  total.f[0]=total.f[0]-2+ pos.p.x[ 0];
   break;

case 44:   // Jeannet_Messine
/*
  Ref. ROADEF 2003, p 273
  i ={0,1,2,3,4,5} (see static data a1 and a2)
  x2 in [-15,25]
  x3 in [3,10]
  min -112.5   in position (3,-7.5, 10)
  The authors define 5 different methods, and
their best result is: 3271 evaluations
TRIBES is _far_ better
*/

 i=pos.p.x[ 0];
x2=pos.p.x[ 1];

x3= pos.p.x[ 2];

total.f[0]=20*a1[i]*x2*x2 + 2*a2[i]*x2*x3;
break;

case 45: // MINLP X. Yan
total.f[0]=MINLP(pos,1); // Constraints
total.f[1]=MINLP(pos,0);// Function to minimize
break;

case 46: // Tripod function (Louis Gacogne)
// on [-100, 100], min 0
x1=pos.p.x[0];
x2= pos.p.x[1];

if(x2<0)
{
	total.f[0]=fabs(x1)+fabs(x2+50);	
}
else
{
	if(x1<0)
		total.f[0]=1+fabs(x1+50)+fabs(x2-50);
	else
		total.f[0]=2+fabs(x1-50)+fabs(x2-50);
}
 break;
 
 case 47: // Tripod function (Louis Gacogne)
// on [-100, 100], min 0
// using hierarchical multiobjective optimization
// constraint= in circle radius 50, center ( 0,-50)
n1=1;
n2=0;
x1=pos.p.x[0];
x2= pos.p.x[1];


 r=    x1*x1+(x2+50)*(x2+50);
 r=sqrt(r);

total.f[n1]=constrain_positive(50-r);
   //printf("\n%f, %f, %f,  %f",x1,x2, r,total.f[0]);
if(x2<0)
{                           
	total.f[n2]=fabs(x1)+fabs(x2+50);
}
else
{
	if(x1<0)
		total.f[n2]=1+fabs(x1+50)+fabs(x2-50);
	else
		total.f[n2]=2+fabs(x1-50)+fabs(x2-50);
}
break;

case 48: // DeJong f4
	//      Global min f(x)=0, x(i)=0

	total.f[0] = 0;
	for( d=0;d<DD;d++)
	{
		total.f[0] = total.f[0] + pos.p.x[ d]*pos.p.x[ d]*pos.p.x[ d]*pos.p.x[ d];
	}
	break;

 case 49: // Pressure vessel 
	 	  // Ref New Optim. Tech. in Eng. p 638
/* D=4
  1.1 <= x1 <= 12.5        granularity 0.0625
  0.6 <= x2 <= 12.5         granularity 0.0625
  0 .0 <= x3 <= 240
  0.0 <= x4 <= 240
  constraints
  g1:= 0.0193*x3-x1 <=0
  g2 := 0.00954*x3-x2<=0
  g3:= 750*1728-pi*x3*x3*(x4-(4/3)*x3)  <=0
*/           
	x1=pos.p.x[0];
	x2= pos.p.x[1];
	x3=pos.p.x[2];
	x4= pos.p.x[3];
 
   f=0.6224*x1*x3*x4 + 1.7781*x2*x3*x3 + 3.1611*x1*x1*x4 + 19.84*x1*x1*x3;
 
  // Constraints
  y=0.0193*x3-x1; if ( y>0) {c= 1+pow(10,10)*y;  f=f*c*c; }
  y=  0.00954*x3-x2; if (y>0) {c=1+y; f=f*c*c;  }
  y = 750*1728-pi*x3*x3*(x4+(4.0/3)*x3);  if (y>0) {c=1+y; f=f*c*c;  }
  
  total.f[0]=f;
 break;  
 
 case 50: // Coil compression spring 
	 	  // Ref New Optim. Tech. in Eng. p 644

	x1=pos.p.x[0];

	x2= pos.p.x[1];
	x3= pos.p.x[2];

	f=pi*pi*x2*x3*x3*(x1+2)*0.25;
	
	// Constraints
	Cf=1+0.75*x3/(x2-x3) + 0.615*x3/x2;
	K=0.125*G*pow(x3,4)/(x1*x2*x2*x2);
	sp=Fp/K;
	lf=Fmax/K + 1.05*(x1+2)*x3;

	y=8*Cf*Fmax*x2/(pi*x3*x3*x3) -S;
	if (y>0) {c=1+y; f=f*c*c*c;}

	y=lf-lmax;
	if (y>0) {c=1+y; f=f*c*c*c;}

	y=sp-spm;
	if (y>0) {c=1+y; f=f*c*c*c;}

	//y=sp+(Fmax-Fp)/K + 1.05*(x1+2)*x3 - lf;
	y=sp-Fp/K;

	if (y>0) {c=1+pow(10,10)*y; f=f*c*c*c;}

	y=sw- (Fmax-Fp)/K;
	if (y>0) {c=1+pow(10,10)*y; f=f*c*c*c;}

  total.f[0]=f;

	 break;


 case 51: // Gear train
	  // Ref New Optim. Tech. in Eng. p 634

	x1=pos.p.x[0];
	x2= pos.p.x[1];
	x3=pos.p.x[2];
	x4= pos.p.x[3];

	f=1/6.931 - x1*x2/(x3*x4); f=f*f;
	total.f[0]=f;
	break;

 case 52: // Pressure vessel as multiobjective (4 functions)
             
	x1=pos.p.x[0];
	x2= pos.p.x[1];
	x3=pos.p.x[2];
	x4= pos.p.x[3];
  
	 total.f[0]=0.6224*x1*x3*x4 + 1.7781*x2*x3*x3 + 3.1611*x1*x1*x4 + 19.84*x1*x1*x3;
 
  // Constraints

	y=0.0193*x3-x1; 
	total.f[1]=y+fabs(y);



	y=  0.00954*x3-x2;
	total.f[2]=y+fabs(y);

	y = 750*1728-pi*x3*x3*(x4+(4.0/3)*x3); 
	total.f[3]=y+fabs(y);  

 break;

 case 53: // Coil compression spring as multiobjective (6 functions)

	x1=pos.p.x[0];
	x2= pos.p.x[1];
	x3= pos.p.x[2];

	 total.f[0]=pi*pi*x2*x3*x3*(x1+2)*0.25;
	
	// Constraints
	Cf=1+0.75*x3/(x2-x3) + 0.615*x3/x2;
	K=0.125*G*pow(x3,4)/(x1*x2*x2*x2);
	sp=Fp/K;
	lf=Fmax/K + 1.05*(x1+2)*x3;

	y=8*Cf*Fmax*x2/(pi*x3*x3*x3) -S;
	total.f[1]=y+fabs(y);

	y=lf-lmax;
	total.f[2]=y+fabs(y);

	y=sp-spm;
	total.f[3]=y+fabs(y);

	y=sp-Fp/K;
	total.f[4]=y+fabs(y);

	y=sw- (Fmax-Fp)/K;
	total.f[5]=y+fabs(y);

	 break;

 case 54: // Pressure vessel by homeomorphism
// WARNING. TO COMPLETE
// May be not possible
	alpha=0.0193;
	beta=0.00954;
	gamma=750*1728;
	delta=240;

	y1=pos.p.x[0];
	y2=pos.p.x[1];
	y3=pos.p.x[2];
	y4=pos.p.x[3];

	x4=y4+delta;

	a_1=0.75*x4;
	b=0.75*(y3-gamma)/pi;

	z=solve_polyn_3(b,0,a_1,1);

	y=z.z[0]; choice=1;

root_choice:

	x3=y-a_1/3;
	x1=alpha*x3-y1;
	x2=beta*x3-y2;
printf("\nx1,x2,x3 %f %f %f",x1,x2,x3);


	if (x1<0 || x2<0)// || x3<0)
	{
		if(z.status>1)
		{
			if (choice==1) {choice=2; y=z.z[1]; goto root_choice;}
			if(choice==2) {choice=3;y=z.z[2]; goto root_choice;}
		}
		printf("\n error my_function case 54. No possible root");
		scanf("%i",&d);
	}


//	x=x3*x3*x3+a_1*x3*x3+b;
//	printf("\n x %f =? 0",x); // Should be equal to 0



   f=0.6224*x1*x3*x4 + 1.7781*x2*x3*x3 + 3.1611*x1*x1*x4 + 19.84*x1*x1*x3;

//	y=gamma-pi*x3*x3*(x4+(4.0/3)*x3); // Should be equal to y3;
//	printf("\nx3 %.12f, x4 %f,  %f =? %f",x3,x4,y,y3);

	printf("\n y %f %f %f %f",y1,y2,y3,y4);
	printf("\n x %f %f %f %f",x1,x2,x3,x4);
	printf("\n f %f",f);

	total.f[0]=f;
	 break;


case 55: // Neural Network Training. 2 Bit Parity (XOR)
	// Dim=9
	// min 0.0
      total.f[0]= ANNXor(pos.p);
break;
      
case 56: // Neural Network Training. 4 Bit Parity
	// Dim=25
	// A good result is 0.11
      total.f[0]= ANNParity4(pos.p);
   break;
case 57: // Neural Network Training. Three Color Cube
	// Dim=46
	// A good result is 0.2
	total.f[0]=ANNCOLORCUBE(pos.p);
	break;
case 58: // Neural Network Training. Diabetes in Pima Indians
	// Dim=64
	// A good result is 0.16
	total.f[0]=ANNPIMA(pos.p);
	break;

case 59: // Neural Network Training. Sin Times Sin
	// Dim=26
	// A good result is 0.23
	total.f[0]= ANNSINSIMP(pos.p);
	break;
case 60: // Neural Network Training. Rise Time Servomechanism
	// Dim=28
	// A good result is 0.45
	total.f[0]= ANNSERVO(pos.p);
	break;
case 61: // Moving Peaks 
	// Warning. Only for mono-objective non recursive optimization
	// (rank 0 in .f and level=0)
	for (d=0;d<DD;d++) gen[d]=pos.p.x[d];
	total.f[0]=eval_movpeaks (gen);
  //printf("\n %f %f",global_max,  total.f[0]);


	if(recent_change==1 && recurs==0)
	{
//printf("\n %f", best_result.f.f[0]);		
		recurs=1; // myfunction will be called again in reinit_swarm
				// so, it is will be then necessary to skip this part
		n_change=n_change+1;
	// Recompute the memorized f values

	reinit_swarm(level,0); // Re-initialise the swarm but keeping the best result
//	reinit_swarm(level,1); // Re-initialise the swarm without keeping the best result

	offline_error+=total_error(best_result.f);
	offline_error_cont+=total_error(best_result.f);
	recent_change=0;
	recurs=0;

	// The current particule has been evaluated twice
	eval_f_tot[level]=eval_f_tot[level]-1; 
	}	 
	break;

case 62: // Goldberg抯 order 3 deceptive problem
	// Be sure dimension=3*k
	// Interval inside ]-0.5 1.5[
	// and granularity=1, so that the only possible x values
	// are 0 and 1
	// The min value is 0
f=0;
/* // Exact binary form
for (d=0;d<DD-2;d=d+3)

{
	bit_sum=0;
	for (i=d;i<d+3;i++) bit_sum+=pos.p.x[i];
	if (bit_sum==0) {f+=0.9; continue;}
	if (bit_sum==1) {f+=0.6; continue;}
	if (bit_sum==2) {f+=0.3; continue;}
	if (bit_sum==3) {f+=1; continue;}
}
*/
/* // Continuous form
for (d=0;d<DD-2;d=d+3)
{
	s=0;
	for (i=d;i<d+3;i++) s+=pos.p.x[i];
	s=MIN(s,3); s=MAX(s,0);
	if(s<=2) f+=0.9-0.3*s;
	else f+=0.3+0.7*(s-2);
}
*/

// Continuous integer form
for (d=0;d<DD-2;d=d+3)
{
	s=0;
	for (i=d;i<d+3;i++) s+=floor(pos.p.x[i]+0.5);
	s=MIN(s,3); s=MAX(s,0);
	if(s<=2) f+=0.9-0.3*s;
	else f+=0.3+0.7*(s-2);
}

total.f[0]=DD/3-f;
break;

case 63: // M黨lenbein's order 5 (MC variant), seen as a 1D function
        // Binary optim. Min 0 on 1057 = 100010001000
D=15; // Length of the bit string
f=0;
bit_x=bits(pos.p.x[0],D);
bit.size=5;
for (d=0;d<D-3;d=d+5)
{
	for (i=d;i<d+5;i++) bit.x[i-d]=bit_x.x[i];
   num=number(0,bit);

   if (num==0) {f+=3.0; continue;}
   if (num==1) {f+=4.0; continue;}
   if (num==3) {f+=2.0; continue;}
   if (num==7) {f+=1.0; continue;}
   if (num==31) {f+=3.5; continue;}   
}

total.f[0]=4*D/5-f;
	break;

 case 64: // Schaffer'f6  (2D)
 	x1=pos.p.x[0];
	x2= pos.p.x[1];
  x1=x1-1; x2=x2-1;    // Optional shift
  r = x1*x1 + x2*x2;
  y=  1 + 0.001*r ;
  y1 = y*y;
  y=sin(sqrt(r));
  y2 = y*y - 0.5;
    total.f[0] = 0.5 + y2/y1;
 break;

//-----------------------------------------------------------
case 99: //  Test functions
// Test
	x1=pos.p.x[0];
	x2= pos.p.x[1];
	f1=2*x1*x2*(1-x2)+ x2*x2;
	f2=2*(1-x1)*x2*(1-x2)+(1-x2)*(1-x2);
	total.f[0]=1-f1;
	break;
// polynomial fitting problem 
		// Cf. Differential Evolution Homepage, C code
		// on [-300 200]^9

	f=0; y=-1;
   dx = 2/dx;
   for (i=0;i<=M;i++)
   {
      py = pos.p.x[0];
      for (j=1;j<DD;j++)
      {
		py = y*py + pos.p.x[j];
      }
      if (py<-1 || py>1) f+=(1-py)*(1-py);
      y+=dx;
   }
   py = pos.p.x[0];
   for (j=1;j<DD;j++) py=1.2*py+pos.p.x[j];
   py = py-72.661;
   if (py<0) f+=py*py;
   py = pos.p.x[0];
   for (j=1;j<DD;j++) py=-1.2*py+pos.p.x[j];
   py =py-72.661;
   if (py<0) f+=py*py;
		total.f[0]=f;
   break;
// Stochastic
	total.f[0] = 0;
	for( d=0;d<DD;d++)
	{
		xd=pos.p.x[d];
		total.f[0]=	total.f[0]+(d+1)*pow(xd,4)+alea_normal(0,1);
	}
	break;	

// Step
	total.f[0] = 0;
	for( d=0;d<DD;d++)
	{
		total.f[0]=	total.f[0]+(int)pos.p.x[d];
	}
	break;

// Beatle (Luc Heylen)
 // Min 0 on (1,1)
 // search space [-3 5]^2 => succes rate 76%
	x1=pos.p.x[0];
	x2= pos.p.x[1];
f = 100*(x2-x1*x1)*(x2-x1*x1)+(1-x1)*(1-x1);
f2=100*(-x2+2-x1*x1)*(-x2+2-x1*x1)+(1-x1)*(1-x1);
//f=log(0.2+f+10*(f*(5.1-mod(sqrt(f+f2),5))))+1.61;
y1=sqrt(f+f2);
y2=y1-5*((int)y1/5);
y3= 0.2+f+10*(f*(5.1-y2));
//f=log(y3)-log(0.2); // Modified function (just one solution)
 f=log(y3)+1.61; // Original function
  total.f[0]=f;
break;
 // Ripple (Luc Heylen)
// Min 0 on (1,1)
// search space [-3 5]^2  => success rate 64%
	x1=pos.p.x[0];
	x2= pos.p.x[1];

y1=1-exp(-3*(x1-1)*(x1-1)-0.3*(x2-1)*(x2-1));
y2 = 100*(x2-x1*x1)*(x2-x1*x1)+(1-x1)*(1-x1);
f=(12+sqrt(y2)/2-10*cos(20*y1)-abs((cos(3*(x1-1))+cos(10*(x2-1)))))/10;
total.f[0]=f;
break;
 

// Dynamic Parabola (Sphere)
  sigma=0.01;
  pr=0.5;
  if (alea(0,1)<pr) r=0; else r=1;
	total.f[0] = 0;
	for( d=0;d<DD;d++)
	{
     //xd= pos.p.x[ d]*(1+r*alea_normal(0,sigma));
     xd= pos.p.x[ d] + r*alea_normal(0,sigma); 
    total.f[0] = total.f[0] +xd*xd ;
	}
  DYN=1;
	break;
// Leon 2D
// on [-1.2 1.2]^2 min 0 on (1 1)
	x1=pos.p.x[0];
	x2= pos.p.x[1];
	f1=x2-x1*x1*x1;
	 total.f[0]=100*f1*f1 +(1-x1)*(1-x1);
	 break;

// Schaffer 2D
// on [-100 100]^2  min 0 on (0 0)
	x1=pos.p.x[0];
	x2= pos.p.x[1];
	f1=1+0.001*(x1*x1+x2*x2);
	f1=f1*f1;
	 total.f[0]=0.5+(sin(sqrt(x1*x1+x2*x2))-0.5)/f1;
	 break;

// Bukin6 2D
// on [-15 -5]x[-3 3] min 0 on (-10 1)
	x1=pos.p.x[0];
	x2= pos.p.x[1];
	 total.f[0]=100*sqrt(abs(x2-0.01*x1*x1)) +0.01*abs(x1+10);
	 break;

// Schwefel 2D
// on [-500 500]^2, min=-837.9658 on (420.9687, 420.9687)
	x1=pos.p.x[0];
	x2= pos.p.x[1];
	 total.f[0]=-x1*sin(sqrt(abs(x1)))-x2*sin(sqrt(abs(x2)));
break;


  // somme i*x_i
  total.f[0] = 0;
	for( d=0;d<DD;d++)
	{
		total.f[0] = total.f[0] + (d+1)*pos.p.x[ d];
	 }

   break;

	 // Livre sur l'OEP. Chapitre Contraintes
	// hom閛morphisme disque/4=carr