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