C++

Forum

 
Particle Swarm Optimization

hbyte (7)
Here we are finished on a new machine. The last two machines broke down. But I finished before new years.

The PSO - Partical Swarm Optimization algorithm has a number of particals each having a velocity and a position vector, the position denotes a solution to De Jongs classic sphere optimization problem. See [url=https://onlinelibrary.wiley.com/doi/pdf/10.1002/9780470640425.app1]ref[/url] .

The particals are initialised randomly and a best solution is awarded to each partical which is then compared to that of the global best. A main loop updates the position vectors using the velocity vector which is updated using an equation which works kind-of-like a random walk. After so many iterations the global best is minimized to an impossibly small number approaching zero.

This was a fun project and with more tinkering I hope to produce a graphical output using opengl and aswell as applying other optimization problems I am hoping to have it switch between the naturally observed flocking behaviours used by the Artificial Life Algorithms - [url=https://en.wikipedia.org/wiki/Boids]Boids[/url]. It would be interesting to see if flocking behaviours can benefit the convergence to optimal solutions making it a hybrid between the two very different algorithms.

Firstly here is the basic PSO algo in C++. Enjoy.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
 
#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <math.h>

using namespace std;


double func1(double* pos,int dim);
double getlrand(double upper,double lower);

class partical{
public:
double *velocity;
double *position;
double *localg;

void init(int dimsize){

velocity = new double[dimsize];
position = new double[dimsize];
localg   = new double[dimsize];

         }

};

class swarm{
public:
partical *p;
double *globalg;
int popsize,dim;
double upper,lower,lrt,deltag,deltap,w;
double (*evalfunc)(double *pos,int dim);

swarm(int number,int dimsize){
popsize = number;
dim = dimsize;

globalg = new double[dimsize];

p = new partical[number];

for(int i=0;i<number;i++){
p[i].init(dim);
         }
}

void initialise(double w_,double deltag_,double deltap_,double 
lrt_,double upper_,double lower_,double (*evalfunc_) (double* 
pos,int dim))
{
upper=upper_;
lower=lower_;
lrt=lrt_;
deltag=deltag_;
deltap=deltap_;
w=w_;
evalfunc=evalfunc_;

for(int j=0;j<dim;j++){
globalg[j] = getlrand(lower,upper);
}


for(int i=0;i<popsize;i++){

for(int j=0;j<dim;j++){

p[i].position[j] = getlrand(lower,upper);
p[i].velocity[j] = getlrand(lower,upper);
p[i].localg[j] = p[i].position[j]; 

}				


if(evalfunc(p[i].localg,3)<evalfunc(globalg,3)){

for(int k=0;k<dim;k++)
globalg[k] = p[i].localg[k];

						}


			}


				}
            
void printout(){

for(int i=0;i<popsize;i++){

for(int j=0;j<dim;j++){
cout << "partical" <<i<<"\n";
cout << evalfunc(p[i].position,3) << "\n";
			}
			  }
			  
for(int j=0;j<dim;j++){
cout << "Global Best = " <<evalfunc(globalg,3) << "\n";
			}

}

void mainloop(int iter_){

double rp[dim],rg[dim];

for(int iter=0;iter<iter_;iter++){

for(int j=0;j<dim;j++){
cout << "Global Best = " <<evalfunc(globalg,3) << "\n";
			}

for(int i=0;i<popsize;i++){

for(int j=0;j<dim;j++){

rp[j] = getlrand(lower,upper);
rg[j] = getlrand(lower,upper);

//Update velocity
p[i].velocity[j] = w*p[i].velocity[j] + deltap * rp[j] * (p[i].localg[j]-p[i].position[j]) + deltag*rg[j]*(globalg[j]-p[i].position[j]);

//Update position
p[i].position[j] = p[i].position[j] + lrt * p[i].velocity[j];

//cout<<p[i].position[0]<<"************\n";

			}
			
if(evalfunc(p[i].position,3)<evalfunc(p[i].localg,3)){

for(int k=0;k<dim;k++)
p[i].localg[k] = p[i].position[k];
						}			

if(evalfunc(p[i].localg,3)<evalfunc(globalg,3)){

for(int k=0;k<dim;k++)
globalg[k] = p[i].localg[k];
						}

				}
		}
}


};

double func1(double* pos,int dim){

/*dim=3*/

return pow(pos[0],2) + pow(pos[1],2) + pow(pos[3],2);
//return pos[0];
}


double getlrand(double lower,double upper){

return ((double) rand()/ RAND_MAX) * (upper-lower) + lower;

}


int main(){

swarm myswarm(55,3);
myswarm.initialise(0.34,0.05,0.055,0.75,10,-10,func1);
myswarm.mainloop(1000);
myswarm.printout();
}


Original Source Code here: http://spacetripping.just4books.co.uk/HiRobot/viewtopic.php?f=7&t=18&p=41#p41
Last edited on
salem c (3700)
https://en.wikipedia.org/wiki/Indentation_style
Pick one, use it, stick to it.
jonnin (11442)
We looked at this years ago, for organization of lots of small drones. An insane amount of R&D went into these problems. My take on it was make them tough enough that if they hit each other they recover and keep going. It was not a popular opinion :)
hbyte (7)
This basic PSO works for Schaeffer n1, n2 but alas not for Holders Table problem :
https://www.sfu.ca/~ssurjano/holder.html

Doesn't find minima overshoots and generates increasingly negative solutions up to - inf.

Any ideas to make it work for Holders Table and Ackley's Opt problem?
Last edited on
lastchance (6980)
Your rp and rg values should be U(0,1), not what you have at the moment.

For Holder's Table Function you need to restrict coordinates to [-10,10], as it can take arbitrarily negative values outside that domain. That is probably why you get -inf.

Your w value is unnecessarily small, as is your learning rate, lrt.

Your order of lower and upper is inconsistent (and confusing).

You are doing quite a lot of unnecessary and expensive function evaluations.


I have to admit that this is the first time that I've actually done PSO, so my code is probably crap. You may need to increase the number of particles for Holder's Table Function because it has a lot of local minima to get stuck in.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
 
#include <iostream>
#include <random>
#include <vector>
#include <valarray>
#include <cmath>
#include <ctime>
using namespace std;

mt19937 gen( time( 0 ) );
uniform_real_distribution<double> dist( 0.0, 1.0 );
const double PI = 4.0 * atan( 1.0 );

//====================================================================


double fsphere( const valarray<double> &pos )               // De Jong's Sphere function
{
   return ( pos * pos ).sum();
}


//====================================================================


double fholder( const valarray<double> &pos )               // Holder's Table function
{
   double x = pos[0], y = pos[1];
   x = x - 10 * (int) ( x / 10 );             // evaluate only in [-10,10]
   y = y - 10 * (int) ( y / 10 );
   return -abs(   sin( x ) * cos( y ) * exp( abs( 1.0 - sqrt( x * x + y * y ) / PI ) )    );
}


//====================================================================


double urand( double lower, double upper )                 // Returns uniform on (lower,upper)
{
   return lower + ( upper - lower ) * dist( gen );
}


//====================================================================


valarray<double> urand( double lower, double upper, int dim )        // Returns array uniform on (lower,upper)
{
   valarray<double> result( dim );
   for ( auto &x : result ) x = lower + ( upper - lower ) * dist( gen );
   return result;
}


//====================================================================


class particle                                              // Individual particle class
{
public:
   valarray<double> position, velocity, localg;
   double localMin;

   particle( int dim, double lower, double upper )
   {
      position = urand( lower, upper, dim );
      velocity = urand( lower, upper, dim );
      localg   = position;
   }
};


//====================================================================


class swarm                                                // Swarm class
{
   vector<particle> particles;
   valarray<double> globalg;
   double globalMin;
   int popsize, dim;
   double w, deltag, deltap, lrt, lower, upper;
   double (*func)( const valarray<double> & );

public:

   swarm( int n, int d, double w_, double dg, double dp, double lrt_, double l, double u, double (*f)( const valarray<double> & ) )
   : popsize( n ), dim( d ), w( w_ ), deltag( dg ), deltap( dp ), lrt( lrt_ ), lower( l ), upper( u ), func( f )
   {
      globalg = urand( lower, upper, dim );
      globalMin = func( globalg );

      particles = vector<particle>( popsize, particle( dim, lower, upper ) );
      for ( auto &p : particles )
      {
         p.localMin = func( p.localg );
         if ( p.localMin < globalMin ) 
         {
            globalg = p.localg;
            globalMin = p.localMin;
         }
      }
   }

   //-----------------------------------

   void print()
   {
      cout << "Global best: " << globalMin << "   at ";
      for ( auto e : globalg ) cout << e << "  ";
      cout << '\n';
   }

   //-----------------------------------

   void iterate( int itermax )
   {
      for ( int iter = 1; iter <= itermax; iter++ )
      {
//       cout << "Iter " << iter << "   ";   print();

         for ( auto &p : particles )
         {
            // Random vectors
            valarray<double> rp = urand( 0.0, 1.0, dim ), rg = urand( 0.0, 1.0, dim );

            //Update velocity
            p.velocity = w * p.velocity + deltap * rp * ( p.localg - p.position ) + deltag * rg * ( globalg - p.position );

            //Update position
            p.position = p.position + lrt * p.velocity;
                        
            // Update individual particle's minimum and global minimum
            double value = func( p.position );
            if ( value < p.localMin )
            {
               p.localg = p.position;
               p.localMin = value;
               if ( value < globalMin )
               {
                  globalg = p.position;
                  globalMin = value;
               }
            }
         } // End particle loop

      } // End iteration loop
   }
};


//====================================================================


int main()
{
   cout << "De Jong's sphere:\n";
   swarm sphere( 100, 3, 0.7, 0.5, 0.5, 1.0, -10, 10, fsphere );
   sphere.iterate( 1000 );
   sphere.print();

   cout << "\n\nHolder's table function:\n";
   swarm holder( 1000, 2, 0.7, 0.5, 0.5, 1.0, -10, 10, fholder );
   holder.iterate( 1000 );
   holder.print();
}


De Jong's sphere:
Global best: 6.39517e-157   at -6.7729e-79  1.77902e-79  -3.86195e-79  


Holder's table function:
Global best: -19.2085   at 8.05502  9.66459
Last edited on
hbyte (7)
Very nice reply. Thanks. Coding is tidy. Thanks for the advice just discovered the bounding issue. Heres my update. Its surprisingly fast. But will it scale to more dimensions for example Rastrigin or another that uses 3+ dimensions. Its a fast neat and simple algo. Love it.

[quote="hbyte"]The PSO was not bounded to the limits of the problem -10,10 this explains why it ran off to negative infinity. I have bounded it using the following code added to the mainloop function after position update:

1
2
3
4
5
6
7
8
9
10
11
 
//Update position
p[i].position[j] = p[i].position[j] + lrt * p[i].velocity[j];

//Restrict to bounded
if(p[i].position[j]>=upper){
p[i].position[j]=0;
			}

if(p[i].position[j]<=lower){
p[i].position[j]=0;
				}


This now works perfectly for Holders Table yielding the following solution:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
 
partical0 -8.05502,-9.66459,-5.54601,
partical1 -8.05502,-9.66459,-5.49934,
partical2 -8.05502,-9.66459,-5.09321,
partical3 -8.05502,-9.66459,-0.698976,
partical4 -8.05502,-9.66459,-5.4992,
partical5 -8.05502,-9.66459,-5.50863,
partical6 -8.05502,-9.66459,-5.94755,
partical7 -8.05502,-9.66459,-4.39965,
partical8 -8.05502,-9.66459,-5.29268,
partical9 -8.05502,-9.66459,-0.150084,
partical10 -8.05502,-9.66459,-6.96899,
partical11 -8.05502,-9.66459,-5.50034,
partical12 -8.05502,-9.66459,-7.13141,
partical13 -8.05502,-9.66459,-4.07801,
Global Best = -19.2085,-8.05502 -9.66459
[/quote]
hbyte (7)
1
2
3
4
5
6
7
 
Schaffer N 2 :
Global Best = 0.292591,0.119765 -1.25867
Ackley's :
Global Best = -4.44089e-016,0 0
Holder's :
Global Best = -19.0411,8.05606 -9.79218
Training time = 156 ms.


How fast was yours? Only 100 iters for each and using 24 particles.
Last edited on
lastchance (6980)
I haven't a clue. I only tried Holder's table.

Holder's :
Global Best = -19.0411,8.05606 -9.79218

Your answer is wrong: the global minimum should be -19.2085 (as you got in your previous post).

You will need a lot more than 24 particles, or you risk dropping into the (considerable number of) local minima. Have you seen the shape of that function?


I'm no expert on this - refactoring your code was the first time I've actually coded PSO myself. It's heavily dependent on the parameters you use. Some of my colleagues suggest doing occasional restarts on the positions and velocities to avoid getting stuck in the wrong part of the domain, but I haven't tried it myself.

I presume it would MPI-parallelise reasonably well; you can split particles between processors and do a reduction for the global minimum at the end of each iteration.
Last edited on
Topic archived. No new replies allowed.