Particle Swarm Optimizer

PSO.h:

// PSO.h
//    Particle Swarm Optimization implementation
//    topology: fips = fully informed particles swarm
//
// original publication:
//    Kennedy, J. and Eberhart, R. (1995)
//    “Particle Swarm Optimization”
//    Proc. of the 1995 IEEE International Conference on Neural Networks
//    pp. 1942-1948, IEEE Press
//
// excellent introduction (for the different velocity update rules see this document):
//    Poli, R.; Kennedy, J. & Blackwell, T.
//    Particle Swarm Optimisation: an overview
//    Swarm Intelligence Journal, 2007, 1, 33-57
//
// further references:
//    e.g. http://en.wikipedia.org/wiki/Particle_swarm_optimization
//
// note
//    this PSO implementations tries to maximize the value given by your custom objective function f(ParticleLocation)
//    if you need to minimize, design your objective function correspondingly!
//    e.g. f(ParticleLocation) --> -f(ParticleLocation)
//
// details
//
// contact persons for this code (add your name here if you augment this code)
//   Jürgen Brauer

#pragma once

#undef max

#include <vector>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <limits.h>


#include "PSO_ObjectiveFunction.h"


using namespace std;



class PSO
{
public:

    

	PSO(int Dims, PSO_ObjectiveFunction* f)
	{		
		NrDims = Dims;
		yourObjectiveFunction = f;

		MinParamVal = new float[NrDims];
		MaxParamVal = new float[NrDims];
		ParamRange  = new float[NrDims];
		ConvergenceThreshold = new float[NrDims];

		OptStep = 0;

		GlobalBestLoc    = new float[NrDims];
		GlobalBestLocVal = -std::numeric_limits<float>::max();
		//GlobalBestLocVal = -1000000.0f;

		total_velVecLen = 0.0f;
		ratio_of_particles_already_converged = 0.0f;

		// default particle velocity update formula to use
		PSOType = ConstrictionCoefficient;
		//PSOType = InertiaWeight;
		//PSOType = RandomWalker;

		// default parameters for inertia weight velocity update formula
		GammaStart = 0.9f;
		GammaStop = 0.4f;
		Gamma = GammaStart;
		GammaReductionRate = 0.001f;
		Phi1  = 1.49618f;
		Phi2  = 1.49618f;

		// default parameters for constriction coefficient velocity update formula
		Xi = 0.72984f;
		//Phi1 = 2.05f;
		//Phi2 = 2.05f;
    Phi1 = 1.75f;
    Phi2 = 1.75f;

    myOutputFileStream = nullptr;
	}



	void setPSOParams_InertiaWeight(float IntertiaStartWeight,
																	float IntertiaStopWeight,
																	float InertiaWeightReductionRate,
																	float AccelerationPersonalBestDir,
																	float AccelerationGlobalBestDir)
	{
		GammaStart = IntertiaStartWeight;
		GammaStop  = IntertiaStopWeight;
		GammaReductionRate = InertiaWeightReductionRate;
		Gamma = GammaStart;

		Phi1 = AccelerationPersonalBestDir;
		Phi2 = AccelerationGlobalBestDir;

	} // setPSOParams_InertiaWeight



	void setPSOParams_ConstrictionCoefficient(float ConstrictionCoefficient,
																						float AccelerationPersonalBestDir,
																						float AccelerationGlobalBestDir)
	{
		Xi = ConstrictionCoefficient;

		Phi1 = AccelerationPersonalBestDir;
		Phi2 = AccelerationGlobalBestDir;

	} // setPSOParams_ConstrictionCoefficient




	void setParamRange(int ParamNr,  float MinVal, float MaxVal, float ConvergedThreshold)
	{
		MinParamVal[ParamNr] = MinVal;
		MaxParamVal[ParamNr] = MaxVal;
		ParamRange[ParamNr]  = MaxVal - MinVal;
		ConvergenceThreshold[ParamNr] = ConvergedThreshold;
	}


	void reset_but_keep_particle_locations(float probToKeepLoc)
	{
		// reset optimization step nr
		OptStep = 0;

		// reset global best location value
		GlobalBestLocVal = -std::numeric_limits<float>::max();

		for (int i=0; i<Particles.size(); i++)
		{			
			// get the i-th particle
			Particle* p = Particles[i];

			// reset personal best location value
			p->PersonalBestLocVal = -std::numeric_limits<float>::max() + 1;

			// reset its velocity vector
			for (int d=0; d<NrDims; d++)
			{
				float ParamRange = MaxParamVal[d]/5 - MinParamVal[d]/5;

				// initialize particle's velocity in that dimension
				float rndVal2 = (float)rand() / (float)RAND_MAX;
				p->vel[d] = -ParamRange + rndVal2 * (2*ParamRange); // get random value in [-range,range]

			} // for (all dimensions of search space)

			// keep the location for this particle or recompute?
			float rnd = (float)rand() / (float)RAND_MAX;
			//if (rnd<probToKeepLoc)
				//continue;
			
		} // for (all particles living in the swarm's population)

	} // reset_but_keep_particle_locations
	


	void generateStartPop(int NrParticles)
	{
		// evaluation function specified?
		if (yourObjectiveFunction == nullptr)
			return;

		// delete particles, if they already exist
		Particles.clear();

		// reset global best location value
		GlobalBestLocVal = -std::numeric_limits<float>::max();
		
		// reset optimization step nr
		OptStep = 0;

		// 1. generate <NrParticles> new particles
		for (int i=0; i<NrParticles; i++)
		{
			// 1.1 create one new particle ...
			Particle* p = new Particle;

			// 1.2 ... and add it to the list of all particles
			Particles.push_back( p );

			currParticleNr = i;

			// 1.3 prepare current + personal best particle location vectors
			p->Loc                = new float[NrDims];
			p->PersonalBestLoc    = new float[NrDims];
			//p->PersonalBestLocVal = -std::numeric_limits<float>::max();
			p->PersonalBestLocVal = -std::numeric_limits<float>::max() + 1;
			p->vel                = new float[NrDims];

			// 1.4 generate random start location for that particle
			for (int d=0; d<NrDims; d++)
			{
				float ParamRange = MaxParamVal[d] - MinParamVal[d];

				// initialize particle's start location
				float rndVal1 = (float)rand() / (float)RAND_MAX;
				float RndVal1InThatRange = rndVal1 * ParamRange;
				float startval_dth_arg = MinParamVal[d] + RndVal1InThatRange;
				p->Loc[d] = startval_dth_arg;

				// initialize particle's best location to start location                    
				p->PersonalBestLoc[d] = p->Loc[d];

				// initialize particle's velocity in that dimension
				float rndVal2 = (float)rand() / (float)RAND_MAX;
				p->vel[d] = -ParamRange + rndVal2 * (2*ParamRange); // get random value in [-range,range]

			} // for (all dimensions of search space)


      // 1.5 compute particle's velocity vector length
      float vecLen = 0.0;
      for (int d=0; d<NrDims; d++)
			{
				vecLen += p->vel[d] * p->vel[d];
			}
      p->velVecLen = sqrt( vecLen );


			// 1.6 evaluate objective function at particle's current location			
			p->PersonalBestLocVal = yourObjectiveFunction->evaluate( p->PersonalBestLoc );

			// 1.7 found a better global best?
			if (p->PersonalBestLocVal > GlobalBestLocVal)
			{
				// yes, found better maximum! -> update global best
				GlobalBestLocVal = p->PersonalBestLocVal;
				for (int d=0; d<NrDims; d++)
					GlobalBestLoc[d] = p->PersonalBestLoc[d];
			}


		} // for (all particles to create)

	} // generateStartPop




	void freeResources()
	{
		for (int i=0; i<Particles.size(); i++)
		{
			// get that particle
			Particle* p = Particles.at(i);

			// free arrays
			delete[] p->Loc;
			delete[] p->PersonalBestLoc;
			delete[] p->vel;
		}

		Particles.clear();

		delete[] MinParamVal;
		delete[] MaxParamVal;
		delete[] GlobalBestLoc;

	} // freeResources



	void optimizeOneStep()
	{
		// evaluation function specified?
		if (yourObjectiveFunction == nullptr)
			return;


		// keep track of nr of optimization steps done
		OptStep++;



		// 1. update particle locations
    //#pragma omp parallel for
		for (int i=0; i<Particles.size(); i++)
		{
			// 1.1 get that particle
			Particle* p = Particles.at(i);

			currParticleNr = i;

			// 1.2 for each dimension of search space:
			// update velocity + location
			for (int d=0; d<NrDims; d++)
			{
				float velOld = p->vel[d];
				float vecToPersonalBest = p->PersonalBestLoc[d] - p->Loc[d];
				float vecToGlobalBest   = GlobalBestLoc[d]      - p->Loc[d];

				float Rnd1 = (float)rand() / (float)RAND_MAX;
				float Rnd2 = (float)rand() / (float)RAND_MAX;

				float ParamRange = MaxParamVal[d] - MinParamVal[d];

				//
				// PSO velocity update formula with inertia weight control for velocity:
				float velNew = 0.0f;
				if (PSOType == InertiaWeight)
					velNew = Gamma * velOld +
					Phi1 * Rnd1 * vecToPersonalBest +
					Phi2 * Rnd2 * vecToGlobalBest;
				else
				//
				// PSO velocity update formula with constriction coefficient for velocity:
				if (PSOType == ConstrictionCoefficient)
					velNew = Xi * (velOld + Phi1 * Rnd1 * vecToPersonalBest + Phi2 * Rnd2 * vecToGlobalBest);
				else
				//
				// Particle performs a Random Walk
				if (PSOType == RandomWalker)
					velNew = (-0.5f + Rnd1) * (ParamRange * 1.0f);

				// make sure, the particle's velocity is not larger than the dimension				
				if (velNew > ParamRange)
					velNew = ParamRange;
				if (velNew < -ParamRange)
					velNew = -ParamRange;

				// update particle's velocity
				p->vel[d] = velNew;




				// update particle's location
				p->Loc[d] = p->Loc[d] + velNew;

				// make sure, d-th coordinate stays in [ MinParamVal[d], MaxParamVal[d] ]
				if (p->Loc[d] < MinParamVal[d])
					p->Loc[d] = MinParamVal[d];
				if (p->Loc[d] > MaxParamVal[d])
					p->Loc[d] = MaxParamVal[d];

			} // for (all dimensions)


      // 1.3 compute particle's velocity vector length
      float vecLen = 0.0;
      for (int d=0; d<NrDims; d++)
			{
				vecLen += p->vel[d] * p->vel[d];
      }
      p->velVecLen = sqrt( vecLen );

            


			// 1.4 now that the particle is at a new location:
			//     evaluate objective function there & check for a better personal/global best

			// get objective function value at current particle's location
            

      // here it is where we normally spent a lot of time!
			// how much time depends on your evaluation / objective function ...
      /////////////////////////////////////////
			float val = yourObjectiveFunction->evaluate( p->Loc );
      /////////////////////////////////////////




			// found better personal best?
			if (val > p->PersonalBestLocVal)
			{
				// yes! store this better personal best value + location
				p->PersonalBestLocVal = val;
				for (int d=0; d<NrDims; d++)
					p->PersonalBestLoc[d] = p->Loc[d];

				// found even a better global best?
        // note: using OpenMP only one particle can change the global best entry! no simultaneous changes!
        //#pragma omp critical
        {   
				    if (val > GlobalBestLocVal)
				    {
					    GlobalBestLocVal = val;
					    for (int d=0; d<NrDims; d++)
						    GlobalBestLoc[d] = p->Loc[d];

				    } // if (better global best found)

         } // pragma omp critical


			} // if (better personal best found)


		} // for (all particles)



		// 2. reduce inertia weight
		Gamma -= GammaReductionRate;
		if (Gamma < GammaStop)
			Gamma = GammaStop;



		// 3. record best values found for each optimization step
		if (OptStep < MAX_NR_STEPS_TO_RECORD_OPT_PROGRESS)
			OptimizationProgress[ OptStep ] = GlobalBestLocVal;

	} // optimizeOneStep




  void writeParticleVelocitiesToFile(string fname)
  {
      // 1. output file stream already created?
      if (myOutputFileStream == nullptr)
      {
          // no! so create it
          myOutputFileStream = new ofstream(fname, ios::out);
      }

      *myOutputFileStream << setw(5) << setfill('0') << OptStep << "\t";

      // 2. for each particle compute its velocity vector length and write it to the output file
      for (int i=0; i<Particles.size(); i++)
			{
				// 2.1 get that particle
				Particle* p = Particles.at(i);

				// 2.2 compute velocity vector length
				float vecLen = 0.0;
				for (int d=0; d<NrDims; d++)
				{
					vecLen += p->vel[d] * p->vel[d];
				}
				vecLen = sqrt( vecLen );

				// 2.3 write vector length to file
				*myOutputFileStream << setw(10) << setfill(' ') << setprecision(4) << vecLen << "\t";

			} // for (all particles)

      // 3. put a new line after writing all N particle velocities
      *myOutputFileStream << endl;
      myOutputFileStream->flush();

  } // writeParticleVelocitiesToFile



  /// checks whether ratio of all particles have converged,
  /// i.e. move no longer
  bool IsConverged(float MinConvergenceRatio)
  {
		total_velVecLen = 0.0f;
    int NrParticlesConverged = 0;
    for (int i=0; i<Particles.size(); i++)
		{
			// get that particle
			Particle* p = Particles.at(i);

			// check for each dimension whether the particle has converged
			bool converged = true;
			for (int d=0; d<NrDims; d++)
			{
				if (p->vel[d] > ConvergenceThreshold[d])
				{				
					converged = false;
					break;
				}
      }
			
			// particle converged?
			if (converged)
					NrParticlesConverged++;
		}
		ratio_of_particles_already_converged = (float)NrParticlesConverged / (float)Particles.size();

		if (ratio_of_particles_already_converged >= MinConvergenceRatio)
				return true;
		else
				return false;

  } // IsConverged



	struct Particle
	{
		float* Loc;                         // particle location

		float* PersonalBestLoc;             // particle personal best location
		float  PersonalBestLocVal;          // value at particle personal best location

		float* vel;                         // velocities (per parameter search space dimension)
    float  velVecLen;                   // velocity vector length
	};




	int                         NrDims;

	float*                      MinParamVal;

	float*                      MaxParamVal;

	float*											ParamRange;

	float*											ConvergenceThreshold;

	vector<Particle*>           Particles;


	int                         OptStep;

	float*                      GlobalBestLoc;

	float                       GlobalBestLocVal;

	int													currParticleNr;


	float												total_velVecLen;

	float												ratio_of_particles_already_converged;


	static const int            MAX_NR_STEPS_TO_RECORD_OPT_PROGRESS = 5000;

	float                       OptimizationProgress[MAX_NR_STEPS_TO_RECORD_OPT_PROGRESS];



	float                       Gamma;                      // inertia weight

	float                       GammaStart;                 // inertia weight start value

	float                       GammaStop;                  // inertia weight stop value

	float                       GammaReductionRate;         // how much to reduce inertia weight after each optimization step?

	float                       Phi1;                       // acceleration personal best direction

	float                       Phi2;                       // acceleration global best direction


	float                       Xi;                         // constriction factor


	
	PSO_ObjectiveFunction*			yourObjectiveFunction;

	float												PARTICLE_HAS_CONVERGED_THRESHOLD;



	enum ExplorationVsExploitationControlType {InertiaWeight, ConstrictionCoefficient, RandomWalker};


	ExplorationVsExploitationControlType    PSOType;


  ofstream*                               myOutputFileStream;

}; // class PSO

PSO_ObjectiveFunction.h:

#pragma once


///
/// basis class for a Particle Swarm Optimizer objective function
///

class PSO_ObjectiveFunction
{
	public:

		PSO_ObjectiveFunction()
		{
			nrEvaluations = 0;
		}


		virtual float				evaluate(float* particleLocation) = 0; // pure virtual


		int									nrEvaluations;											   // we count how often the function has been evaluated


}; // PSO_ObjectiveFunction
 
public/particle_swarm_optimizer.txt · Last modified: 2016/03/24 13:56 by root · []
Recent changes RSS feed Powered by PHP Valid XHTML 1.0 Valid CSS Driven by DokuWiki