/* Standard PSO 2007 Contact for remarks, suggestions etc.: MauriceClerc@WriteMe.com Last update 2010-01-04 Fixed wrong fitness evaluation for G3 (function code 10) 2009-12-29 Random number generator KISS (option). For more reproducible results 2009-07-12 The initialisation space may be smaller than the search space (for tests) 2009-06-03 Fixed a small mistake about the first best position 2009-05-05 Step function 2009-04-19 Schaffer f6, 2D Goldstein-Price 2009-03-31 A small network optimisation 2009-03-12 Schwefel 2.2, Neumaier 3, G3 (constrained) 2008-10-02 Two Schwefel functions 2008-08-12 For information: save the best position over all runs 2007-12-10 Warning about rotational invariance (valid here only on 2D) 2007-11-22 stop criterion (option): distance to solution < epsilon and log_progress evaluation 2007-11-21 Ackley function -------------------------------- Contributors The works and comments of the following persons have been taken into account while designing this standard. Sometimes this is for including a feature, and sometimes for leaving out one. Auger, Anne Blackwell, Tim Bratton, Dan Clerc, Maurice Croussette, Sylvain Dattasharma, Abhi Eberhart, Russel Hansen, Nikolaus Keko, Hrvoje Kennedy, James Krohling, Renato Langdon, William Li, Wentao Liu, Hongbo Miranda, Vladimiro Poli, Riccardo Serra, Pablo Stickel, Manfred -------------------------------- MotivationQuite often, researchers claim to compare their version of PSO with the "standard one", but the "standard one" itself seems to vary!Thus, it is important to define a real standard that would stay unchanged for at least one year.This PSO version does not intend to be the best one on the market(in particular, there is no adaptation of the swarm size nor of thecoefficients). This is simply very near to the original version (1995),with just a few improvements based on some recent works. --------------------------------- Metaphorsswarm: A team of communicating people (particles)At each time step Each particle chooses a few informants at random, selects the best one from this set, and takes into account the information given by the chosen particle. If it finds no particle better than itself, then the "reasoning" is: "I am the best, so I just take my current velocity and my previous best position into account" ----------------------------------- Parameters/Optionsclamping := true/false => whether to use clamping positions or notrandOrder:= true/false => whether to avoid the bias due to the loop on particles "for s = 1 to swarm_size ..." or notrotation := true/false => whether the algorithm is sensitive to a rotation of the landscape or not You may also modify the following ones, although suggested valuesare either hard coded or automatically computed:S := swarm sizeK := maximum number of particles _informed_ by a given onew := first cognitive/confidence coefficientc := second cognitive/confidence coefficient ----------------------------------- EquationsFor each particle and each dimensionEquation 1: v(t+1) = w*v(t) + R(c)*(p(t)-x(t)) + R(c)*(g(t)-x(t))Equation 2: x(t+1) = x(t) + v(t+1)wherev(t) := velocity at time tx(t) := position at time tp(t) := best previous position of the particleg(t) := best position amongst the best previous positions of the informants of the particleR(c) := a number coming from a random distribution, which depends on cIn this standard, the distribution is uniform on [0,c]Note 1:When the particle has no informant better than itself,it implies p(t) = g(t)Therefore, Equation 1 gets modified to:v(t+1) = w*v(t) + R(c)*(p(t)-x(t))Note 2:When the "non sensitivity to rotation" option is activated(p(t)-x(t)) (and (g(t)-x(t))) are replaced by rotated vectors, so that the final DNPP (Distribution of the Next Possible Positions)is not dependent on the system of co-ordinates. ----------------------------------- Information links topology A lot of work has been done about this topic. The main result is this: There is no "best" topology. Hence the random approach used here. ----------------------------------- InitialisationInitial positions are chosen at random inside the search space (which is supposed to be a hyperparallelepiped, and often evena hypercube), according to a uniform distribution.This is not the best way, but the one used in the original PSO.Each initial velocity is simply defined as the half-difference of tworandom positions. It is simple, and needs no additional parameter.However, again, it is not the best approach. The resulting distributionis not even uniform, as is the case for any method that uses auniform distribution independently for each component.The mathematically correct approach needs to use a uniformdistribution inside a hypersphere. It is not very difficult,and was indeed used in some PSO versions. However, it is quitedifferent from the original one. Moreover, it may be meaningless for some heterogeneous problems,when each dimension has a different "interpretation".------------------------------------ From SPSO-06 to SPSO-07The main differences are:1. option "non sensitivity to rotation of the landscape" Note: although theoretically interesting, this option is quite computer time consuming, and the improvement in result may only be marginal. 2. option "random permutation of the particles before each iteration" Note: same remark. Time consuming, no clear improvement3. option "clamping position or not" Note: in a few rare cases, not clamping positions may induce an infinite run, if the stop criterion is the maximum number of evaluations 4. probability p of a particular particle being an informant of another particle. In SPSO-06 it was implicit (by building the random infonetwork) Here, the default value is directly computed as a function of (S,K), so that the infonetwork is exactly the same as in SPSO-06. However, now it can be "manipulated" ( i.e. any value can be assigned) 5. The search space can be quantised (however this algorithm is _not_ for combinatorial problems)Also, the code is far more modular. It means it is slower, but easierto translate into another language, and easier to modify. ----------------------------------- Use Define the problem (you may add your own one in problemDef() and perf()) Choose your options Run and enjoy! */ #include "stdio.h" #include "math.h" #include #include #define D_max 85 // Max number of dimensions of the search space #define R_max 500 // Max number of runs #define S_max 910 // Max swarm size #define zero 0 // 1.0e-30 // To avoid numerical instabilities #define ulong unsigned long // To generate pseudo-random numbers with KISS#define RAND_MAX_KISS ((unsigned long) 4294967295) // Structuresstruct quantum { double q[D_max]; int size; }; struct SS { int D; double max[D_max]; double maxInit[D_max]; double min[D_max]; double minInit[D_max]; struct quantum q; // Quantisation step size. 0 => continuous problem }; struct param { double c; // Confidence coefficient int clamping; // Position clamping or not int K; // Max number of particles informed by a given one double p; // Probability threshold for random topology // (is actually computed as p(S,K) ) int randOrder; // Random choice of particles or not int rand; // 0 => use KISS. Any other value: use the standard C RNG int initLink; // How to re-init links int rotation; // Sensitive to rotation or not int S; // Swarm size int stop; // Flag for stop criterion double w; // Confidence coefficient }; struct position { double f; int improved; int size; double x[D_max];}; struct velocity { int size; double v[D_max]; };struct problem { double epsilon; // Admissible error int evalMax; // Maximum number of fitness evaluations int function; // Function code double objective; // Objective value // Solution position (if known, just for tests) struct position solution; struct SS SS; // Search space }; struct swarm { int best; // rank of the best particle struct position P[S_max]; // Previous best positions found by each particle int S; // Swarm size struct velocity V[S_max]; // Velocities struct position X[S_max]; // Positions }; struct result { double nEval; // Number of evaluations struct swarm SW; // Final swarm double error; // Numerical result of the run};struct matrix // Useful for "non rotation sensitive" option{ int size; double v[D_max][D_max]; }; // Sub-programs ulong rand_kiss(); // For the pseudo-random number generator KISSvoid seed_rand_kiss(ulong seed); double alea (double a, double b);int alea_integer (int a, int b);double alea_normal (double mean, double stdev);double distanceL(struct position x1, struct position x2, double L);struct velocity aleaVector(int D,double coeff);struct matrix matrixProduct(struct matrix M1,struct matrix M2);struct matrix matrixRotation(struct velocity V);struct velocity matrixVectProduct(struct matrix M,struct velocity V);double normL (struct velocity v,double L);double perf (struct position x, int function,struct SS SS, double objective); // Fitness evaluationstruct position quantis (struct position x, struct SS SS);struct problem problemDef(int functionCode);struct result PSO ( struct param param, struct problem problem);int sign (double x); // Global variables long double E; // exp(1). Useful for some test functions long double pi; // Useful for some test functions int randOption;struct matrix reflex1;long double sqrtD; // For Network problem int bcsNb; int btsNb; // File(s); FILE * f_run;FILE * f_synth; // ================================================= int main () { struct position bestBest; // Best position over all runs int d; // Current dimension double error; // Current error double errorMean; // Average error double errorMin; // Best result over all runs double errorMeanBest[R_max]; double evalMean; // Mean number of evaluations int functionCode; int i,j; int nFailure; // Number of unsuccessful runs double logProgressMean; struct param param; struct problem pb; int run, runMax; struct result result; double successRate; double variance; f_run = fopen ("f_run.txt", "w"); f_synth = fopen ("f_synth.txt", "w"); E = exp ((long double) 1); pi = acos ((long double) -1); // ----------------------------------------------- PROBLEM functionCode = 10; /* (see problemDef( ) for precise definitions) 0 Parabola (Sphere) 1 Griewank 2 Rosenbrock (Banana) 3 Rastrigin 4 Tripod (dimension 2) 5 Ackley 6 Schwefel 7 Schwefel 1.2 8 Schwefel 2.21 9 Neumaier 3 10 G3 11 Network optimisation (Warning: see problemDef() and also perf() for problem elements (number of BTS and BSC) 12 Schwefel 13 2D Goldstein-Price 14 Schaffer f6 15 Step 100 Shifted Parabola 99 Test */ runMax = 1; // Numbers of runs if (runMax > R_max) runMax = R_max; // ----------------------------------------------------- // PARAMETERS // * means "suggested value" param.clamping =1; // 0 => no clamping AND no evaluation. WARNING: the program // may NEVER stop (in particular with option move 20 (jumps)) 1 // *1 => classical. Set to bounds, and velocity to zero param.initLink = 0; // 0 => re-init links after each unsuccessful iteration // 1 => re-init links after each successful iteration param.rand=0; // 0 => Use KISS as random number generator. // Any other value => use the standard C one param.randOrder=1; // 0 => at each iteration, particles are modified // always according to the same order 0..S-1 //*1 => at each iteration, particles numbers are // randomly permutated param.rotation=0; // WARNING. Experimental code, completely valid only for dimension 2 // 0 => sensitive to rotation of the system of coordinates // 1 => non sensitive (except side effects), // by using a rotated hypercube for the probability distribution // WARNING. Quite time consuming! param.stop = 0; // Stop criterion // 0 => error < pb.epsilon // 1 => eval >= pb.evalMax // 2 => ||x-solution|| < pb.epsilon // ------------------------------------------------------- // Some information printf ("\n Function %i ", functionCode); printf("\n (clamping, randOrder, rotation, stop_criterion) = (%i, %i, %i, %i)", param.clamping, param.randOrder, param.rotation, param.stop); if(param.rand==0) printf("\n WARNING, I am using the RNG KISS"); // =========================================================== // RUNs // Initialize some objects pb=problemDef(functionCode); // You may "manipulate" S, p, w and c // but here are the suggested values param.S = (int) (10 + 2 * sqrt(pb.SS.D)); // Swarm size if (param.S > S_max) param.S = S_max; printf("\n Swarm size %i", param.S); param.K=3; param.p=1-pow(1-(double)1/(param.S),param.K); // (to simulate the global best PSO, set param.p=1) //param.p=1; // According to Clerc's Stagnation Analysis param.w = 1 / (2 * log ((double) 2)); // 0.721 param.c = 0.5 + log ((double) 2); // 1.193 // According to Poli's Sampling Distribution of PSOs analysis //param.w = ??; // in [0,1[ //param.c = // smaller than 12*(param.w*param.w-1)/(5*param.w -7); printf("\n c = %f, w = %f",param.c, param.w); randOption=param.rand; // Global variable //--------------- sqrtD=sqrt((long double) pb.SS.D); // Define just once the first reflexion matrix if(param.rotation>0) { reflex1.size=pb.SS.D; for (i=0;i pb.epsilon) // Failure { nFailure = nFailure + 1; } // Memorize the best (useful if more than one run) if(run==0) bestBest=result.SW.P[result.SW.best]; else if(error 1) { printf ("\n Best min value = %f", errorMin); printf ("\nPosition of the optimum: "); for (d=0;dpb.evalMax) R.SW.S=pb.evalMax; R.nEval = R.SW.S; // Find the best R.SW.best = 0; switch (param.stop) { default: errorPrev =R.SW.P[R.SW.best].f; // "distance" to the wanted f value (objective) break; case 2: errorPrev=distanceL(R.SW.P[R.SW.best],pb.solution,2); // Distance to the wanted solution break; } for (s = 1; s < R.SW.S; s++) { switch (param.stop) { default: zz=R.SW.P[s].f; if (zz < errorPrev) { R.SW.best = s; errorPrev=zz; } break; case 2: zz=distanceL(R.SW.P[R.SW.best],pb.solution,2); if (zz= R.SW.S) s1 = s; // Find the best informant g = s1; for (m = s1; m < R.SW.S; m++) { if (LINKS[m][s] == 1 && R.SW.P[m].f < R.SW.P[g].f) g = m; } //.. compute the new velocity, and move // Exploration tendency for (d = 0; d < pb.SS.D; d++) { R.SW.V[s].v[d]=param.w *R.SW.V[s].v[d]; } // Prepare Exploitation tendency p-x for (d = 0; d < pb.SS.D; d++) { PX.v[d]= R.SW.P[s].x[d] - R.SW.X[s].x[d]; } PX.size=pb.SS.D; if(g!=s) { for (d = 0; d < pb.SS.D; d++) // g-x { GX.v[d]= R.SW.P[g].x[d] - R.SW.X[s].x[d]; } GX.size=pb.SS.D; } // Option "non sentivity to rotation" if (param.rotation>0) { normPX=normL(PX,2); if (g!=s) normGX=normL(GX,2); if(normPX>0) { RotatePX=matrixRotation(PX); } if(g!= s && normGX>0) { RotateGX=matrixRotation(GX); } } // Exploitation tendencies switch (param.rotation) { default: for (d = 0; d < pb.SS.D; d++) { R.SW.V[s].v[d]=R.SW.V[s].v[d] + + alea(0, param.c)*PX.v[d]; } if (g!=s) { for (d = 0; d < pb.SS.D; d++) { R.SW.V[s].v[d]=R.SW.V[s].v[d] + alea(0,param.c) * GX.v[d]; } } break; case 1: // First exploitation tendency if(normPX>0) { zz=param.c*normPX/sqrtD; aleaV=aleaVector(pb.SS.D, zz); expt1=matrixVectProduct(RotatePX,aleaV); for (d = 0; d < pb.SS.D; d++) { R.SW.V[s].v[d]=R.SW.V[s].v[d]+expt1.v[d]; } } // Second exploitation tendency if(g!=s && normGX>0) { zz=param.c*normGX/sqrtD; aleaV=aleaVector(pb.SS.D, zz); expt2=matrixVectProduct(RotateGX,aleaV); for (d = 0; d < pb.SS.D; d++) { R.SW.V[s].v[d]=R.SW.V[s].v[d]+expt2.v[d]; } } break; } // Update the position for (d = 0; d < pb.SS.D; d++) { R.SW.X[s].x[d] = R.SW.X[s].x[d] + R.SW.V[s].v[d]; } if (R.nEval >= pb.evalMax) { //error= fabs(error - pb.objective); goto end; } // -------------------------- noEval = 1; // Quantisation R.SW.X[s] = quantis (R.SW.X[s], pb.SS); switch (param.clamping) { case 0: // No clamping AND no evaluation outside = 0; for (d = 0; d < pb.SS.D; d++) { if (R.SW.X[s].x[d] < pb.SS.min[d] || R.SW.X[s].x[d] > pb.SS.max[d]) outside++; } if (outside == 0) // If inside, the position is evaluated { R.SW.X[s].f = perf (R.SW.X[s], pb.function, pb.SS,pb.objective); R.nEval = R.nEval + 1; } break; case 1: // Set to the bounds, and v to zero for (d = 0; d < pb.SS.D; d++) { if (R.SW.X[s].x[d] < pb.SS.min[d]) { R.SW.X[s].x[d] = pb.SS.min[d]; R.SW.V[s].v[d] = 0; } if (R.SW.X[s].x[d] > pb.SS.max[d]) { R.SW.X[s].x[d] = pb.SS.max[d]; R.SW.V[s].v[d] = 0; } } R.SW.X[s].f =perf(R.SW.X[s],pb.function, pb.SS,pb.objective); R.nEval = R.nEval + 1; break; } // ... update the best previous position if (R.SW.X[s].f < R.SW.P[s].f) // Improvement { R.SW.P[s] = R.SW.X[s]; // ... update the best of the bests if (R.SW.P[s].f < R.SW.P[R.SW.best].f) { R.SW.best = s; } } } // End of "for (s0=0 ... " // Check if finished switch (param.stop) { default: error = R.SW.P[R.SW.best].f; break; case 2: error=distanceL(R.SW.P[R.SW.best],pb.solution,2); break; } //error= fabs(error - pb.epsilon); if (error < errorPrev) // Improvement { initLinks = 0; } else // No improvement { initLinks = 1; // Information links will be reinitialized } if(param.initLink==1) initLinks=1-initLinks; errorPrev = error; end: switch (param.stop) { case 0: case 2: if (error > pb.epsilon && R.nEval < pb.evalMax) { noStop = 0; // Won't stop } else { noStop = 1; // Will stop } break; case 1: if (R.nEval < pb.evalMax) noStop = 0; // Won't stop else noStop = 1; // Will stop break; } } // End of "while nostop ... // printf( "\n and the winner is ... %i", R.SW.best ); // fprintf( f_stag, "\nEND" ); R.error = error; return R; } // =========================================================== double alea (double a, double b) { // random number (uniform distribution) in [a b] // randOption is a global parameter double r; if(randOption==0) r=a+(double)rand_kiss()*(b-a)/RAND_MAX_KISS; else r=a + (double) rand () * (b - a)/RAND_MAX; // printf("\nalea 929 r %f",r); return r; } // =========================================================== int alea_integer (int a, int b) { // Integer random number in [a b] int ir; double r; r = alea (0, 1); ir = (int) (a + r * (b + 1 - a)); if (ir > b) ir = b; return ir; } // =========================================================== double alea_normal (double mean, double std_dev) { /* Use the polar form of the Box-Muller transformation to obtain a pseudo random number from a Gaussian distribution */ double x1, x2, w, y1; // double y2; do { x1 = 2.0 * alea (0, 1) - 1.0; x2 = 2.0 * alea (0, 1) - 1.0; w = x1 * x1 + x2 * x2; } while (w >= 1.0); w = sqrt (-2.0 * log (w) / w); y1 = x1 * w; // y2 = x2 * w; if(alea(0,1)<0.5) y1=-y1; y1 = y1 * std_dev + mean; return y1; } // =============================================================struct velocity aleaVector(int D,double coeff){ struct velocity V; int d; int i; int K=2; // 1 => uniform distribution in a hypercube // 2 => "triangle" distribution double rnd; V.size=D; for (d=0;d Euclidean int d; double n; n = 0; for (d = 0; d < x1.size; d++) n = n + pow (fabs(x1.x[d] - x2.x[d]), L); n = pow (n, 1/L); return n; } //============================================================struct matrix matrixProduct(struct matrix M1,struct matrix M2){ // Two square matrices of same size struct matrix Product; int D; int i,j,k; double sum; D=M1.size; for (i=0;i V where V'=(1,1,...1)*normV/sqrt(D) (i.e. norm(V') = norm(V) ) */ struct velocity B={0}; int i,j,d, D; double normB,normV,normV2; //struct matrix reflex1; // Global variable struct matrix reflex2; struct matrix rotateV; double temp; D=V.size; normV=normL(V,2); normV2=normV*normV; reflex2.size=D; // Reflection relatively to the vector V'=(1,1, ...1)/sqrt(D) // norm(V')=1 // Has been computed just once (global matrix reflex1) //Define the "bisectrix" B of (V',V) as an unit vector B.size=D; temp=normV/sqrtD; for (d=0;d0) { for (d=0;d rotation rotateV=matrixProduct(reflex2,reflex1); return rotateV;}//==========================================================struct velocity matrixVectProduct(struct matrix M,struct velocity V){ struct velocity Vp; int d,j; int Dim; double sum; Dim=V.size; for (d=0;d zero) // Note that qd can't be < 0 { //qd = qd * (SS.max[d] - SS.min[d]) / 2; quantx.x[d] = qd * floor (0.5 + x.x[d] / qd); } } return quantx; } //================================================== KISS/* A good pseudo-random numbers generator The idea is to use simple, fast, individually promising generators to get a composite that will be fast, easy to code have a very long period and pass all the tests put to it. The three components of KISS are x(n)=a*x(n-1)+1 mod 2^32 y(n)=y(n-1)(I+L^13)(I+R^17)(I+L^5), z(n)=2*z(n-1)+z(n-2) +carry mod 2^32 The y's are a shift register sequence on 32bit binary vectors period 2^32-1; The z's are a simple multiply-with-carry sequence with period 2^63+2^32-1. The period of KISS is thus 2^32*(2^32-1)*(2^63+2^32-1) > 2^127 */static ulong kiss_x = 1;static ulong kiss_y = 2;static ulong kiss_z = 4;static ulong kiss_w = 8;static ulong kiss_carry = 0;static ulong kiss_k;static ulong kiss_m;void seed_rand_kiss(ulong seed) { kiss_x = seed | 1; kiss_y = seed | 2; kiss_z = seed | 4; kiss_w = seed | 8; kiss_carry = 0;}ulong rand_kiss() { kiss_x = kiss_x * 69069 + 1; kiss_y ^= kiss_y << 13; kiss_y ^= kiss_y >> 17; kiss_y ^= kiss_y << 5; kiss_k = (kiss_z >> 2) + (kiss_w >> 3) + (kiss_carry >> 2); kiss_m = kiss_w + kiss_w + kiss_z + kiss_carry; kiss_z = kiss_w; kiss_w = kiss_m; kiss_carry = kiss_k >> 30; return kiss_x + kiss_y + kiss_w;} // =========================================================== double perf (struct position x, int function, struct SS SS, double objective) { // Evaluate the fitness value for the particle of rank s int d; double DD; int k; int n; double f, p, xd, x1, x2; double s11, s12, s21, s22; double sum1,sum2; double t0, tt, t1; struct position xs; // Shifted Parabola/Sphere (CEC 2005 benchmark) static double offset_0[30] = { -3.9311900e+001, 5.8899900e+001, -4.6322400e+001, -7.4651500e+001, -1.6799700e+001, -8.0544100e+001, -1.0593500e+001, 2.4969400e+001, 8.9838400e+001, 9.1119000e+000, -1.0744300e+001, -2.7855800e+001, -1.2580600e+001, 7.5930000e+000, 7.4812700e+001, 6.8495900e+001, -5.3429300e+001, 7.8854400e+001, -6.8595700e+001, 6.3743200e+001, 3.1347000e+001, -3.7501600e+001, 3.3892900e+001, -8.8804500e+001, -7.8771900e+001, -6.6494400e+001, 4.4197200e+001, 1.8383600e+001, 2.6521200e+001, 8.4472300e+001 }; /* static float bts[5][2]= { {6.8, 9.0}, {8.3, 7.9}, {6.6, 5.6}, {10, 5.4}, {8, 3} }; */ static float bts [19][2]= { {6, 9}, {8, 7}, {6, 5}, {10, 5}, {8, 3} , {12, 2}, {4, 7}, {7, 3}, {1, 6}, {8, 2}, {13, 12}, {15, 7}, {15, 11}, {16, 6}, {16, 8}, {18, 9}, {3, 7}, {18, 2}, {20, 17} }; float btsPenalty= 100; double z1,z2; xs = x; switch (function) { case 100: for (d = 0; d < xs.size; d++) { xs.x[d]=xs.x[d]-offset_0[d]; } case 0: // Parabola (Sphere) f = 0; for (d = 0; d < xs.size; d++) { xd = xs.x[d]; f = f + xd * xd; } break; case 1: // Griewank f = 0; p = 1; for (d = 0; d < xs.size; d++) { xd = xs.x[d]; f = f + xd * xd; p = p * cos (xd / sqrt ((double) (d + 1))); } f = f / 4000 - p + 1; break; case 2: // Rosenbrock f = 0; t0 = xs.x[0] + 1; // Solution on (0,...0) when // offset=0 for (d = 1; d < xs.size; d++) { t1 = xs.x[d] + 1; tt = 1 - t0; f += tt * tt; tt = t1 - t0 * t0; f += 100 * tt * tt; t0 = t1; } break; case 3: // Rastrigin k = 10; f = 0; for (d = 0; d < xs.size; d++) { xd = xs.x[d]; f =f+ xd * xd - k * cos (2 * pi * xd); } f =f+ xs.size * k; break; case 4: // 2D Tripod function // Note that there is a big discontinuity right on the solution // point. x1 = xs.x[0] ; x2 = xs.x[1]; s11 = (1.0 - sign (x1)) / 2; s12 = (1.0 + sign (x1)) / 2; s21 = (1.0 - sign (x2)) / 2; s22 = (1.0 + sign (x2)) / 2; //f = s21 * (fabs (x1) - x2); // Solution on (0,0) f = s21 * (fabs (x1) +fabs(x2+50)); // Solution on (0,-50) f = f + s22 * (s11 * (1 + fabs (x1 + 50) + fabs (x2 - 50)) + s12 * (2 + fabs (x1 - 50) + fabs (x2 - 50))); break; case 5: // Ackley sum1=0; sum2=0; DD=x.size; pi=acos(-1); for (d=0;d1+zero) f=f+btsPenalty; } // Distances for(d=0;d