GOclasses/algorithms/ASA.cpp

   1 /*
   2  *  ASA.cpp
   3  *  SeGMO, a Sequential Global Multiobjective Optimiser
   4  *
   5  *  Created by Dario Izzo on 5/16/08.
   6  *  Copyright 2008 ¿dvanced Concepts Team (European Space Agency). All rights reserved.
   7  *
   8  */
   9
  10 #include "ASA.h"
  11 #include "PkRandom.h"
  12 #include <iostream>
  13
  14 Population ASAalgorithm::evolve(Individual x0, GOProblem& problem) {
  15
  16     const std::vector<double>& LB = problem.getLB();
  17     const std::vector<double>& UB = problem.getUB();
  18
  19         std::vector<double> xNEW = x0.getDecisionVector();
  20         std::vector<double> xOLD = xNEW;
  21         double fNEW;
  22         fNEW = x0.getFitness();
  23         double fOLD;
  24         fOLD = fNEW;
  25         SolDim = xNEW.size();
  26         std::vector<double> Step(SolDim,StartStep);
  27         std::vector<int> acp(SolDim,0) ;
  28         double ratio=0;
  29         double currentT = T0;
  30         double prob=0;
  31         double r=0;
  32
  33         //Main SA loops
  34         for ( int jter=0; jter < niterOuter; jter++){
  35                 for ( int mter = 0; mter < niterTemp; mter++){
  36                         for ( int kter = 0 ; kter < niterRange; kter++){
  37                                 int nter =  (int)(drng()*SolDim);
  38                                 for ( int numb = 0; numb < SolDim ;numb++){
  39                                         nter=(nter+1) % SolDim;
  40                                         //We modify the current point actsol by mutating its nter component within
  41                                         //a Step that we will later adapt
  42                                         r = 2.0*drng()-1.0; //random number in [-1,1]
  43                                         xNEW[nter] = xOLD[nter] + r * Step[nter] * ( UB[nter] - LB[nter] );
  44
  45                                         // If new solution produced is infeasible ignore it
  46                                         if ((xNEW[nter] > UB[nter]) || (xNEW[nter] < LB[nter])){
  47                                                 xNEW[nter]=xOLD[nter];
  48                                                 continue;
  49                                         }
  50                                         //And we valuate the objective function for the new point
  51                                         fNEW = problem.objfun(xNEW);
  52
  53                                         // We decide wether to accept or discard the point
  54                                         if (fNEW < fOLD)  //accept
  55                                         {
  56                                                 xOLD[nter] = xNEW[nter];
  57                                                 fOLD = fNEW;
  58                                                 acp[nter]++;    //Increase the number of accepted values
  59                                         }
  60                                         else  //test it with Boltzmann to decide the acceptance
  61                                         {
  62
  63                                                 prob = exp ( (fOLD - fNEW )/ currentT );
  64
  65                                                 // we compare prob with a random probability.
  66                                                 if (prob > drng())
  67                                                 {
  68                                                         xOLD[nter] = xNEW[nter];
  69                                                         fOLD = fNEW;
  70                                                         acp[nter]++;    //Increase the number of accepted values
  71                                                 }
  72                                                 else{
  73                                                         xNEW[nter] = xOLD[nter];
  74                                                 }
  75                                         } // end if
  76                                 } // end for(nter = 0; ...
  77                         } // end for(kter = 0; ...
  78
  79
  80                         // adjust the step (adaptively)
  81                         for (int iter=0; iter < SolDim; iter++)
  82                         {
  83                                 ratio = (double)acp[iter]/(double)niterRange;
  84                                 acp[iter]=0;  //reset the counter
  85                                 if (ratio > .6) //too many acceptances, increase the step by a factor 3 maximum
  86                                 {
  87                                         Step[iter] = Step [iter] * (1 + 2 *(ratio - .6)/.4);
  88                                 }
  89                                 else
  90                                 {
  91                                         if (ratio < .4) //too few acceptance, decrease the step by a factor 3 maximum
  92                                         {
  93                                                 Step [iter]= Step [iter] / (1 + 2 * ((.4 - ratio)/.4));
  94                                         };
  95                                 };
  96                                 //And if it becomes too large, reset it to its initial value
  97                                 if ( Step[iter] > StartStep ){
  98                                                 Step [iter] = StartStep;
  99                                 };
 100                         }
 101                 }
 102                 // Cooling schedule
 103                 currentT *= Tcoeff;
 104         }
 105
 106         Population newpop;
 107         if (fOLD < x0.getFitness()){
 108         x0.setDecisionVector( xOLD );
 109         x0.setFitness( fOLD );
 110         }
 111         newpop.addIndividual( x0 );
 112         return newpop;
 113         }
 114
 115         void ASAalgorithm::initASA(int niterTotInit, int niterTempInit, int niterRangeInit, int SolDimInit, double T0Init, double TcoeffInit, double StartStepInit, uint32_t randomSeed){
 116                 niterTot=niterTotInit;
 117                 niterTemp=niterTempInit;
 118                 niterRange=niterRangeInit;
 119                 SolDim=SolDimInit;
 120                 T0=T0Init;
 121                 Tcoeff=TcoeffInit;
 122                 StartStep=StartStepInit;
 123                 niterOuter = niterTot / (niterTemp * niterRange * SolDim);
 124                 drng.seed(randomSeed);
 125         }
 126
 127         void ASAalgorithm::initASA(int niterTotInit,int SolDimInit, double Ts, double Tf, uint32_t randomSeed){
 128                 niterTot=niterTotInit;
 129                 niterTemp=1;
 130                 niterRange=20;
 131                 SolDim=SolDimInit;
 132                 niterOuter = niterTot / (niterTemp * niterRange * SolDim);
 133                 T0=Ts;
 134                 Tcoeff=pow(Tf/Ts,1.0/(double)(niterOuter));
 135                 StartStep=1;
 136                 drng.seed(randomSeed);
 137         }