GOclasses/algorithms/MPSO.cpp

   1 /*
   2  *  MPSO.cpp
   3  *  PaGMO
   4  *
   5  *  Created by Dario Izzo on 10/23/08.
   6  *  Copyright 2008 __MyCompanyName__. All rights reserved.
   7  *
   8  */
   9
  10 #include "MPSO.h"
  11 #include "vector"
  12
  13 using namespace std;
  14
  15 void MPSOalgorithm::initMPSO(int generationsInit, int SolDimInit, double omegaInit, double eta1Init, double eta2Init,double vcoeffInit, int nswarmsInit, unsigned long randomSeed){
  16         generations = generationsInit;
  17         SolDim = SolDimInit;
  18         omega = omegaInit;
  19         eta1 = eta1Init;
  20         eta2 = eta2Init;
  21         vcoeff = vcoeffInit;
  22         nswarms = nswarmsInit;
  23         rng = Pk::Random32(randomSeed);
  24 }
  25
  26 Population MPSOalgorithm::evolve(Population deme, GOProblem& problem){
  27
  28     const std::vector<double>& LB = problem.getLB();
  29     const std::vector<double>& UB = problem.getUB();
  30
  31         int NP = deme.size()/nswarms;                                                   //potentially dangerous when the deme size is not divisible by the numberof swamrs
  32         int D = LB.size();
  33
  34         vector<double> dummy(D,0);                                                              //used for initialisation purposes
  35         vector<double> dummy2(NP,0);                                                    //used for initialisation purposes
  36         vector< vector<double> > dummy3(NP,dummy);                              //used for initialisation purposes
  37
  38         //Particle position: X[i][j][k] = i-th swarm, j-th individual, k-th phenotype
  39         vector< vector< vector<double> > > X(nswarms,dummy3);
  40         //Particle velocity: V[i][j][k] = i-th swarm, j-th individual, k-th phenotype
  41         vector< vector< vector<double> > > V(nswarms,dummy3);
  42         //Particle fitness: fit[i][j] = i-th swarm, j-th individual
  43         vector< vector<double> > fit(nswarms,dummy2);
  44         //Global best swarm fitness: gbfit[i] = i-th swarm
  45         vector<double> gbfit(nswarms,0);
  46         //Global best swarm individual: gbX[i][j] = i-th swarm, j-th phenotype
  47         vector< vector<double> > gbX(nswarms,dummy);
  48
  49         //Local best fitness: lbfit[i][j] = i-th swarm, j-th individual
  50         vector< vector <double> > lbfit(nswarms,dummy2);
  51         //Local best chromosome: lbX[i][j][k] =  i-th swarm, j-th individual, k-th phenotype
  52         vector< vector< vector<double> > > lbX(nswarms,dummy3);
  53
  54         double vwidth;                                                                          //Width of the search space
  55     vector<double> MINV(D),MAXV(D);                                             //Maximum and minumum velocity allowed
  56
  57         // Initialise the particles (class Individual) positions, their velocities and their fitness to that of the deme
  58         for ( int i = 0; i<nswarms; i++){
  59                 for ( int j = 0; j<NP; j++ ){
  60                                 X[i][j] =       deme[NP*i+j].getDecisionVector();
  61                                 V[i][j] =       deme[NP*i+j].getVelocity();
  62                                 fit[i][j]       =       deme[NP*i+j].getFitness();
  63                 }
  64         }
  65    // Initialise the minimum and maximum velocity
  66    for ( int i = 0; i<D; i++ ) {
  67                 vwidth = (UB[i]-LB[i])/vcoeff;
  68                 MINV[i] = -1.0*vwidth;
  69                 MAXV[i] = vwidth;
  70         }
  71
  72    for (int i=0;i<nswarms;i++){
  73         // Initialise the global and local bests
  74         gbX[i]=X[i][0];
  75         gbfit[i]=fit[i][0];
  76
  77         lbX[i]=X[i];                    //at the first generation the local best position is the particle position
  78         lbfit[i]=fit[i];                //same for the fitness
  79
  80         for (int j = 1; j<NP; j++){             //the int j = 1 jumps the first member as it is already set as the best
  81                 if (fit[i][j] < gbfit[i]){
  82                         gbfit[i] = fit[i][j];
  83                         gbX[i] = X[i][j];
  84                 }
  85         }
  86    }
  87
  88    // Main PSO loop
  89    for (int iter = 0; iter<generations; iter++){
  90         //loop through the swarms
  91         for (int i = 0;i<nswarms;i++){
  92                 //1 - move the particles and check that velocity and positions are in allowed range
  93                 for (int j = 0; j< NP; j++){
  94                         for (int k = 0; k< D; k++){
  95
  96                                 //new velocity
  97                                 V[i][j][k] = omega * V[i][j][k] + eta1 * Pk::nextDouble(rng) * (lbX[i][j][k] - X[i][j][k]) + eta2 * Pk::nextDouble(rng) * (gbX[i][k] - X[i][j][k]);
  98
  99                                 //check that it is within the allowed velocity range
 100                                 if ( V[i][j][k] > MAXV[k] )
 101                                         V[i][j][k] = MAXV[k];
 102
 103                                 else if ( V[i][j][k] < MINV[k] )
 104                                         V[i][j][k] = MINV[k];
 105
 106                                 //update position
 107                                 X[i][j][k] = X[i][j][k] + V[i][j][k];
 108
 109                                 if (X[i][j][k] < LB[k])
 110                                         X[i][j][k] = Pk::nextDouble(rng) * (UB[k] - LB[k]) + LB[k];
 111                                 else if (X[i][j][k] > UB[k])
 112                                         X[i][j][k] = Pk::nextDouble(rng) * (UB[k] - LB[k]) + LB[k];
 113                         }
 114
 115                         //We evaluate the new individual fitness now as to be able to update immediately the global best
 116                         //in case a better solution is found
 117                         fit[i][j] = problem.objfun(X[i][j]);
 118                         //update local and global best
 119                         if (fit[i][j] < lbfit[i][j]){
 120                                 lbfit[i][j] = fit[i][j];        //local best
 121                                 lbX[i][j] = X[i][j];
 122                                 if(fit[i][j] < gbfit[i]){
 123                                         gbfit[i] = fit[i][j];   //global best
 124                                         gbX[i]  = X[i][j];
 125                                 }
 126                         }
 127                 } //End of loop on the population members
 128         } //End of loop on the swarms
 129
 130         //exchanges two random elements from randomly selected swarms
 131         if (iter % (int)5 == 0){
 132                 int sw1 = (int)(Pk::nextDouble(rng)*nswarms);           //select 1st swarm
 133                 int sw2 = (int)(Pk::nextDouble(rng)*nswarms);
 134                 do{                                                                                     //endless loop if nswarms<2
 135                   sw2 = (int)(Pk::nextDouble(rng)*nswarms);             //selects 2nd swarm different from the first
 136                 } while (sw2 == sw1);
 137
 138                 int in1 = (int)(Pk::nextDouble(rng)*NP);                //select 1st individual
 139                 int in2;
 140                 do{                                                                                     //endless loop if nswarms<2
 141                   in2 = (int)(Pk::nextDouble(rng)*NP);              //selects 2nd individual
 142                 } while (in2 == in1);
 143                 //swap position
 144                 dummy = X[sw1][in1];
 145                 X[sw1][in1]=X[sw2][in2];
 146                 X[sw2][in2] = dummy;
 147
 148                 //swap velocity
 149                 dummy = V[sw1][in1];
 150                 V[sw1][in1]=V[sw2][in2];
 151                 V[sw2][in2] = dummy;
 152
 153                 //swap local bests lbX
 154                 dummy = lbX[sw1][in1];
 155                 lbX[sw1][in1]=lbX[sw2][in2];
 156                 lbX[sw2][in2] = dummy;
 157
 158                 //swap local best lbfit;
 159                 dummy[0] = lbfit[sw1][in1];
 160                 lbfit[sw1][in1]=lbfit[sw2][in2];
 161                 lbfit[sw2][in2] = dummy[0];
 162         }//end of swap
 163    } // end of main PSO loop
 164
 165    //we end by constructing the object Population containing the final results
 166    Population popout;
 167    Individual temp;
 168    for (int i=0; i<nswarms;i++){
 169    for (int j=0; j<NP; j++){
 170         temp.setDecisionVector(lbX[i][j]);
 171         temp.setFitness(lbfit[i][j]);
 172         temp.setVelocity(V[i][j]);
 173         popout.addIndividual(temp);
 174    }
 175    }
 176    return popout;
 177 }
 178