Rosetta Code

Particle swarm optimization

Revision as of 17:44, 15 July 2018 by Robbie (talk | contribs) (→‎{{header|Java}})

Particle Swarm Optimization (PSO) is an optimization method in which multiple candidate solutions ('particles') migrate through the solution space under the influence of local and global best known positions. PSO does not require that the objective function be differentiable and can optimize over very large problem spaces, but is not guaranteed to converge. The method should be demonstrated by application of the functions recommended below, and possibly other standard or well-known optimization test cases.

Particle swarm optimization is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

The goal of parameter selection is to ensure that the global minimum is discriminated from any local minima, and that the minimum is accurately determined, and that convergence is achieved with acceptible resource usage. To provide a common basis for comparing implementations, the following test cases are recommended:

  • McCormick function - bowl-shaped, with a single minimum
      function parameters and bounds (recommended):
    • -1.5 < x1 < 4
    • -3 < x2 < 4
      search parameters (suggested):
    • omega = 0
    • phi p = 0.6
    • phi g = 0.3
    • number of particles = 100
    • number of iterations = 40
  • Michalewicz function - steep ridges and valleys, with multiple minima
      function parameters and bounds (recommended):
    • 0 < x1 < pi
    • 0 < x2 < pi
      search parameters (suggested):
    • omega = 0.3
    • phi p = 0.3
    • phi g = 0.3
    • number of particles = 1000
    • number of iterations = 30

References:

  • [Particle Swarm Optimization[1]]
  • [Virtual Library of Optimization Test Functions[2]]

D

Translation of: Kotlin

<lang D>import std.math; import std.random; import std.stdio;

alias Func = double function(double[]);

struct Parameters { 

   double omega, phip, phig;

}

struct State { 

   int iter;
   double[] gbpos;
   double gbval;
   double[] min;
   double[] max;
   Parameters parameters;
   double[][] pos;
   double[][] vel;
   double[][] bpos;
   double[] bval;
   int nParticles;
   int nDims;

   void report(string testfunc) {
       writeln("Test Function        : ", testfunc);
       writeln("Iterations           : ", iter);
       writefln("Global Best Position : [%(%.16f, %)]", gbpos);
       writefln("Global Best Value    : %.16f", gbval);
   }

}

State psoInit(double[] min, double[] max, Parameters parameters, int nParticles) { 

   auto nDims = min.length;
   double[][] pos;
   pos.length = nParticles;
   pos[] = min;
   double[][] vel;
   vel.length = nParticles;
   for (int i; i<nParticles; i++) vel[i].length = nDims;
   double[][] bpos;
   bpos.length = nParticles;
   bpos[] = min;
   double[] bval;
   bval.length = nParticles;
   bval[] = double.infinity;
   auto iter = 0;
   double[] gbpos;
   gbpos.length = nDims;
   gbpos[] = double.infinity;
   auto gbval = double.infinity;
   return State(iter, gbpos, gbval, min, max, parameters, pos, vel, bpos, bval, nParticles, nDims);

}

State pso(Func fn, State y) { 

   auto p = y.parameters;
   double[] v;
   v.length = y.nParticles;
   double[][] bpos;
   bpos.length = y.nParticles;
   bpos[] = y.min;
   double[] bval;
   bval.length = y.nParticles;
   double[] gbpos;
   gbpos.length = y.nDims;
   auto gbval = double.infinity;
   foreach (j; 0..y.nParticles) {
       // evaluate
       v[j] = fn(y.pos[j]);
       // update
       if (v[j] < y.bval[j]) {
           bpos[j] = y.pos[j];
           bval[j] = v[j];
       } else {
           bpos[j] = y.bpos[j];
           bval[j] = y.bval[j];
       }
       if (bval[j] < gbval) {
           gbval = bval[j];
           gbpos = bpos[j];
       }
   }
   auto rg = uniform01();
   double[][] pos;
   pos.length = y.nParticles;
   for (int i; i<pos.length; i++) pos[i].length = y.nDims;
   double[][] vel;
   vel.length = y.nParticles;
   for (int i; i<vel.length; i++) vel[i].length = y.nDims;
   foreach (j; 0..y.nParticles) {
       // migrate
       auto rp = uniform01();
       bool ok = true;
       vel[j][] = 0;
       pos[j][] = 0;
       foreach (k; 0..y.nDims) {
           vel[j][k] = p.omega * y.vel[j][k] +
                       p.phip * rp * (bpos[j][k] - y.pos[j][k]) +
                       p.phig * rg * (gbpos[k] - y.pos[j][k]);
           pos[j][k] = y.pos[j][k] + vel[j][k];
           ok = ok && y.min[k] < pos[j][k] && y.max[k] > pos[j][k];
       }
       if (!ok) {
           foreach (k; 0..y.nDims) {
               pos[j][k] = y.min[k] + (y.max[k] - y.min[k]) * uniform01();
           }
       }
   }
   auto iter = 1 + y.iter;
   return State(iter, gbpos, gbval, y.min, y.max, y.parameters, pos, vel, bpos, bval, y.nParticles, y.nDims);

}

State iterate(Func fn, int n, State y) { 

   auto r = y;
   auto old = y;
   if (n == int.max) {
       while (true) {
           r = pso(fn, r);
           if (r == old) break;
           old = r;
       }
   } else {
       foreach (_; 0..n) r = pso(fn, r);
   }
   return r;

}

double mccormick(double[] x) { 

   auto a = x[0];
   auto b = x[1];
   return sin(a + b) + (a - b) * (a - b) + 1.0 + 2.5 * b - 1.5 * a;

}

double michalewicz(double[] x) { 

   auto m = 10;
   auto d = x.length;
   auto sum = 0.0;
   foreach (i; 1..d) {
       auto j = x[i - 1];
       auto k = sin(i * j * j / PI);
       sum += sin(j) * k^^(2.0*m);
   }
   return -sum;

}

void main() { 

   auto state = psoInit(
       [-1.5, -3.0],
       [4.0, 4.0],
       Parameters(0.0, 0.6, 0.3),
       100
   );
   state = iterate(&mccormick, 40, state);
   state.report("McCormick");
   writefln("f(-.54719, -1.54719) : %.16f", mccormick([-.54719, -1.54719]));
   writeln;
   state = psoInit(
       [0.0, 0.0],
       [PI, PI],
       Parameters(0.3, 0.3, 0.3),
       1000
   );
   state = iterate(&michalewicz, 30, state);
   state.report("Michalewicz (2D)");
   writefln("f(2.20, 1.57)        : %.16f", michalewicz([2.2, 1.57]));

}</lang>

Output:
Test Function        : McCormick
Iterations           : 40
Global Best Position : [-0.5673174452967942, -1.5373177402652800]
Global Best Value    : -1.9122776571457756
f(-.54719, -1.54719) : -1.9132229548822735

Test Function        : Michalewicz (2D)
Iterations           : 30
Global Best Position : [2.1907380816516597, 1.5608474620076016]
Global Best Value    : -1.7949374368688056
f(2.20, 1.57)        : -1.8011407184738251

J

<lang J>load 'format/printf'

pso_init =: verb define 

  'Min Max parameters nParticles' =. y
  'Min: %j\nMax: %j\nomega, phip, phig: %j\nnParticles: %j\n' printf Min;Max;parameters;nParticles
  nDims =. #Min
  pos =. Min +"1 (Max - Min) *"1 (nParticles,nDims) ?@$ 0
  bpos =. pos
  bval =. (#pos) $ _
  vel  =. ($pos) $ 0
  0;_;_;Min;Max;parameters;pos;vel;bpos;bval      NB. initial state

)

pso =: adverb define 

  NB. previous state
  'iter gbpos gbval Min Max parameters pos vel bpos0 bval' =. y 

  NB. evaluate
  val    =. u"1 pos

  NB. update
  better =. val < bval
  bpos   =. (better # pos) (I. better)} bpos0
  bval   =. u"1 bpos
  gbval  =. <./ bval
  gbpos  =. bpos {~ (i. <./) bval

  NB. migrate
  'omega phip phig' =. parameters
  rp  =. (#pos) ?@$ 0
  rg  =. ? 0
  vel =. (omega*vel) + (phip * rp * bpos - pos) + (phig * rg * gbpos -"1 pos)
  pos =. pos + vel

  NB. reset out-of-bounds particles
  bad    =. +./"1 (Min >"1 pos) ,. (pos >"1 Max)
  newpos =. Min +"1 (Max-Min) *"1 ((+/bad),#Min) ?@$ 0
  pos    =. newpos (I. bad)} pos
  iter   =. >: iter

  NB. new state
  iter;gbpos;gbval;Min;Max;parameters;pos;vel;bpos;bval

)

reportState=: 'Iteration: %j\nGlobalBestPosition: %j\nGlobalBestValue: %j\n' printf 3&{.</lang> Apply to McCormick Function:<lang J> require 'trig' 

  mccormick =: sin@(+/) + *:@(-/) + 1 + _1.5 2.5 +/@:* ]

  state =: pso_init _1.5 _3 ; 4 4 ; 0 0.6 0.3; 100

Min: _1.5 _3 Max: 4 4 omega, phip, phig: 0 0.6 0.3 nParticles: 100 

  state =: (mccormick pso)^:40 state
  reportState state

Iteration: 40 GlobalBestPosition: _0.547399 _1.54698 GlobalBestValue: _1.91322</lang> Apply to Michalewicz Function: <lang J> michalewicz =: 3 : '- +/ (sin y) * 20 ^~ sin (>: i. #y) * (*:y) % pi' 

  michalewicz =: [: -@(+/) sin * 20 ^~ sin@(pi %~ >:@i.@# * *:)  NB. tacit equivalent
   
  state =: pso_init 0 0 ; (pi,pi) ; 0.3 0.3 0.3; 1000

Min: 0 0 Max: 3.14159 3.14159 omega, phip, phig: 0.3 0.3 0.3 nParticles: 1000 

  state =: (michalewicz pso)^:30 state
  reportState state

Iteration: 30 GlobalBestPosition: 2.20296 1.57083 GlobalBestValue: _1.8013</lang>

JavaScript

Translation of J.

<lang JavaScript>function pso_init(y) { 

 var nDims= y.min.length;
 var pos=[], vel=[], bpos=[], bval=[];
 for (var j= 0; j<y.nParticles; j++) {
   pos[j]= bpos[j]= y.min;
   var v= []; for (var k= 0; k<nDims; k++) v[k]= 0;
   vel[j]= v;
   bval[j]= Infinity}
 return {

iter: 0, gbpos: Infinity, gbval: Infinity, min: y.min, max: y.max, parameters: y.parameters, pos: pos, vel: vel, bpos: bpos, bval: bval, 

       nParticles: y.nParticles,
       nDims: nDims}

}

function pso(fn, state) { 

 var y= state;
 var p= y.parameters;
 var val=[], bpos=[], bval=[], gbval= Infinity, gbpos=[];
 for (var j= 0; j<y.nParticles; j++) {
   // evaluate
   val[j]= fn.apply(null, y.pos[j]);
   // update
   if (val[j] < y.bval[j]) {
     bpos[j]= y.pos[j];
     bval[j]= val[j];
   } else {
     bpos[j]= y.bpos[j];
     bval[j]= y.bval[j]}
   if (bval[j] < gbval) {
     gbval= bval[j];
     gbpos= bpos[j]}}
 var rg= Math.random(), vel=[], pos=[];
 for (var j= 0; j<y.nParticles; j++) {
   // migrate
   var rp= Math.random(), ok= true;
   vel[j]= [];
   pos[j]= [];
   for (var k= 0; k < y.nDims; k++) {
     vel[j][k]= p.omega*y.vel[j][k] + p.phip*rp*(bpos[j]-y.pos[j]) + p.phig*rg*(gbpos-y.pos[j]);
     pos[j][k]= y.pos[j]+vel[j][k];
     ok= ok && y.min[k]<pos[j][k] && y.max>pos[j][k];}
   if (!ok)
     for (var k= 0; k < y.nDims; k++)
       pos[j][k]= y.min[k] + (y.max[k]-y.min[k])*Math.random()}
 return {

iter: 1+y.iter, gbpos: gbpos, gbval: gbval, min: y.min, max: y.max, parameters: y.parameters, pos: pos, vel: vel, bpos: bpos, bval: bval, 

       nParticles: y.nParticles,
       nDims: y.nDims}

}

function display(text) { 

 if (document) {
   var o= document.getElementById('o');
   if (!o) {
     o= document.createElement('pre');
     o.id= 'o';
     document.body.appendChild(o)}
   o.innerHTML+= text+'\n';
   window.scrollTo(0,document.body.scrollHeight);
 }
 if (console.log) console.log(text)

}

function reportState(state) { 

 var y= state;
 display();
 display('Iteration: '+y.iter);
 display('GlobalBestPosition: '+y.gbpos);
 display('GlobalBestValue: '+y.gbval);

}

function repeat(fn, n, y) { 

 var r=y, old= y;
 if (Infinity == n)
   while ((r= fn(r)) != old) old= r;
 else
   for (var j= 0; j<n; j++) r= fn(r);
 return r

}

function mccormick(a,b) { 

 return Math.sin(a+b) + Math.pow(a-b,2) + (1 + 2.5*b - 1.5*a)

}

state= pso_init({ 

 min: [-1.5,-3], max:[4,4],
 parameters: {omega: 0, phip: 0.6, phig: 0.3},
 nParticles: 100});

reportState(state);

state= repeat(function(y){return pso(mccormick,y)}, 40, state);

reportState(state);</lang>

Example displayed result (random numbers are involved so there will be a bit of variance between repeated runs:

<lang Javascript> Iteration: 0 GlobalBestPosition: Infinity GlobalBestValue: Infinity

Iteration: 40 GlobalBestPosition: -0.5134004259016365,-1.5512442672625184 GlobalBestValue: -1.9114053788600853</lang>

Java

Translation of: Kotlin

<lang java>import java.util.Arrays; import java.util.Objects; import java.util.Random; import java.util.function.Function;

public class App { 

   static class Parameters {
       double omega;
       double phip;
       double phig;

       Parameters(double omega, double phip, double phig) {
           this.omega = omega;
           this.phip = phip;
           this.phig = phig;
       }
   }

   static class State {
       int iter;
       double[] gbpos;
       double gbval;
       double[] min;
       double[] max;
       Parameters parameters;
       double[][] pos;
       double[][] vel;
       double[][] bpos;
       double[] bval;
       int nParticles;
       int nDims;

       State(int iter, double[] gbpos, double gbval, double[] min, double[] max, Parameters parameters, double[][] pos, double[][] vel, double[][] bpos, double[] bval, int nParticles, int nDims) {
           this.iter = iter;
           this.gbpos = gbpos;
           this.gbval = gbval;
           this.min = min;
           this.max = max;
           this.parameters = parameters;
           this.pos = pos;
           this.vel = vel;
           this.bpos = bpos;
           this.bval = bval;
           this.nParticles = nParticles;
           this.nDims = nDims;
       }

       void report(String testfunc) {
           System.out.printf("Test Function        : %s\n", testfunc);
           System.out.printf("Iterations           : %d\n", iter);
           System.out.printf("Global Best Position : %s\n", Arrays.toString(gbpos));
           System.out.printf("Global Best value    : %.15f\n", gbval);
       }
   }

   private static State psoInit(double[] min, double[] max, Parameters parameters, int nParticles) {
       int nDims = min.length;
       double[][] pos = new double[nParticles][];
       for (int i = 0; i < nParticles; ++i) {
           pos[i] = min.clone();
       }
       double[][] vel = new double[nParticles][nDims];
       double[][] bpos = new double[nParticles][];
       for (int i = 0; i < nParticles; ++i) {
           bpos[i] = min.clone();
       }
       double[] bval = new double[nParticles];
       for (int i = 0; i < bval.length; ++i) {
           bval[i] = Double.POSITIVE_INFINITY;
       }
       int iter = 0;
       double[] gbpos = new double[nDims];
       for (int i = 0; i < gbpos.length; ++i) {
           gbpos[i] = Double.POSITIVE_INFINITY;
       }
       double gbval = Double.POSITIVE_INFINITY;
       return new State(iter, gbpos, gbval, min, max, parameters, pos, vel, bpos, bval, nParticles, nDims);
   }

   private static Random r = new Random();

   private static State pso(Function<double[], Double> fn, State y) {
       Parameters p = y.parameters;
       double[] v = new double[y.nParticles];
       double[][] bpos = new double[y.nParticles][];
       for (int i = 0; i < y.nParticles; ++i) {
           bpos[i] = y.min.clone();
       }
       double[] bval = new double[y.nParticles];
       double[] gbpos = new double[y.nDims];
       double gbval = Double.POSITIVE_INFINITY;
       for (int j = 0; j < y.nParticles; ++j) {
           // evaluate
           v[j] = fn.apply(y.pos[j]);
           // update
           if (v[j] < y.bval[j]) {
               bpos[j] = y.pos[j];
               bval[j] = v[j];
           } else {
               bpos[j] = y.bpos[j];
               bval[j] = y.bval[j];
           }
           if (bval[j] < gbval) {
               gbval = bval[j];
               gbpos = bpos[j];
           }
       }
       double rg = r.nextDouble();
       double[][] pos = new double[y.nParticles][y.nDims];
       double[][] vel = new double[y.nParticles][y.nDims];
       for (int j = 0; j < y.nParticles; ++j) {
           // migrate
           double rp = r.nextDouble();
           boolean ok = true;
           Arrays.fill(vel[j], 0.0);
           Arrays.fill(pos[j], 0.0);
           for (int k = 0; k < y.nDims; ++k) {
               vel[j][k] = p.omega * y.vel[j][k] +
                   p.phip * rp * (bpos[j][k] - y.pos[j][k]) +
                   p.phig * rg * (gbpos[k] - y.pos[j][k]);
               pos[j][k] = y.pos[j][k] + vel[j][k];
               ok = ok && y.min[k] < pos[j][k] && y.max[k] > pos[j][k];
           }
           if (!ok) {
               for (int k = 0; k < y.nDims; ++k) {
                   pos[j][k] = y.min[k] + (y.max[k] - y.min[k]) * r.nextDouble();
               }
           }
       }
       int iter = 1 + y.iter;
       return new State(
           iter, gbpos, gbval, y.min, y.max, y.parameters,
           pos, vel, bpos, bval, y.nParticles, y.nDims
       );
   }

   private static State iterate(Function<double[], Double> fn, int n, State y) {
       State r = y;
       if (n == Integer.MAX_VALUE) {
           State old = y;
           while (true) {
               r = pso(fn, r);
               if (Objects.equals(r, old)) break;
               old = r;
           }
       } else {
           for (int i = 0; i < n; ++i) {
               r = pso(fn, r);
           }
       }
       return r;
   }

   private static double mccormick(double[] x) {
       double a = x[0];
       double b = x[1];
       return Math.sin(a + b) + (a - b) * (a - b) + 1.0 + 2.5 * b - 1.5 * a;
   }

   private static double michalewicz(double[] x) {
       int m = 10;
       int d = x.length;
       double sum = 0.0;
       for (int i = 1; i < d; ++i) {
           double j = x[i - 1];
           double k = Math.sin(i * j * j / Math.PI);
           sum += Math.sin(j) * Math.pow(k, 2.0 * m);
       }
       return -sum;
   }

   public static void main(String[] args) {
       State state = psoInit(
           new double[]{-1.5, -3.0},
           new double[]{4.0, 4.0},
           new Parameters(0.0, 0.6, 0.3),
           100
       );
       state = iterate(App::mccormick, 40, state);
       state.report("McCormick");
       System.out.printf("f(-.54719, -1.54719) : %.15f\n", mccormick(new double[]{-.54719, -1.54719}));
       System.out.println();

       state = psoInit(
           new double[]{0.0, 0.0},
           new double[]{Math.PI, Math.PI},
           new Parameters(0.3, 3.0, 0.3),
           1000
       );
       state = iterate(App::michalewicz, 30, state);
       state.report("Michalewicz (2D)");
       System.out.printf("f(2.20, 1.57)        : %.15f\n", michalewicz(new double[]{2.20, 1.57}));
   }

}</lang>

Output:
Test Function        : McCormick
Iterations           : 40
Global Best Position : [-0.5468738679864172, -1.547048532862534]
Global Best value    : -1.913222827709136
f(-.54719, -1.54719) : -1.913222954882274

Test Function        : Michalewicz (2D)
Iterations           : 30
Global Best Position : [2.2029055320517994, 1.832848319327826]
Global Best value    : -0.801303410098550
f(2.20, 1.57)        : -0.801166387820286

Kotlin

Translation of: JavaScript

<lang scala>// version 1.1.51

import java.util.Random

typealias Func = (DoubleArray) -> Double

class Parameters(val omega: Double, val phip: Double, val phig: Double)

class State( 

   val iter: Int,
   val gbpos: DoubleArray,
   val gbval: Double,
   val min: DoubleArray,
   val max: DoubleArray,
   val parameters: Parameters,
   val pos: Array<DoubleArray>,
   val vel: Array<DoubleArray>,
   val bpos: Array<DoubleArray>,
   val bval: DoubleArray,
   val nParticles: Int,
   val nDims: Int

) { 

   fun report(testfunc: String) {
       println("Test Function        : $testfunc")
       println("Iterations           : $iter")
       println("Global Best Position : ${gbpos.asList()}")
       println("Global Best Value    : $gbval")
   }

}

fun psoInit( 

   min: DoubleArray,
   max: DoubleArray,
   parameters: Parameters,
   nParticles: Int

): State { 

   val nDims = min.size
   val pos   = Array(nParticles) { min }
   val vel   = Array(nParticles) { DoubleArray(nDims) }
   val bpos  = Array(nParticles) { min }
   val bval  = DoubleArray(nParticles) { Double.POSITIVE_INFINITY}
   val iter  = 0
   val gbpos = DoubleArray(nDims) { Double.POSITIVE_INFINITY }
   val gbval = Double.POSITIVE_INFINITY
   return State(iter, gbpos, gbval, min, max, parameters,
                pos, vel, bpos, bval, nParticles, nDims)

}

val r = Random()

fun pso(fn: Func, y: State): State { 

   val p = y.parameters
   val v = DoubleArray(y.nParticles)
   val bpos  = Array(y.nParticles) { y.min }
   val bval  = DoubleArray(y.nParticles)
   var gbpos = DoubleArray(y.nDims)
   var gbval = Double.POSITIVE_INFINITY
   for (j in 0 until y.nParticles) {
       // evaluate
       v[j] = fn(y.pos[j])
       // update
       if (v[j] < y.bval[j]) {
           bpos[j] = y.pos[j]
           bval[j] = v[j]
       }
       else {
           bpos[j] = y.bpos[j]
           bval[j] = y.bval[j]
       }
       if (bval[j] < gbval) {
           gbval = bval[j]
           gbpos = bpos[j]
       }
   }
   val rg = r.nextDouble()
   val pos = Array(y.nParticles) { DoubleArray(y.nDims) }
   val vel = Array(y.nParticles) { DoubleArray(y.nDims) }
   for (j in 0 until y.nParticles) {
       // migrate
       val rp = r.nextDouble()
       var ok = true
       vel[j].fill(0.0)
       pos[j].fill(0.0)
       for (k in 0 until y.nDims) {
           vel[j][k] = p.omega * y.vel[j][k] +
                       p.phip * rp * (bpos[j][k] - y.pos[j][k]) +
                       p.phig * rg * (gbpos[k] - y.pos[j][k])
           pos[j][k] = y.pos[j][k] + vel[j][k]
           ok = ok && y.min[k] < pos[j][k] && y.max[k] > pos[j][k]
       }
       if (!ok) {
           for (k in 0 until y.nDims) {
               pos[j][k]= y.min[k] + (y.max[k] - y.min[k]) * r.nextDouble()
           }
       }
   }
   val iter = 1 + y.iter
   return State(
       iter, gbpos, gbval, y.min, y.max, y.parameters,
       pos, vel, bpos, bval, y.nParticles, y.nDims
   )

}

fun iterate(fn: Func, n: Int, y: State): State { 

   var r = y
   var old = y
   if (n == Int.MAX_VALUE) {
       while (true) {
           r = pso(fn, r)
           if (r == old) break
           old = r
       }
   }
   else {
       repeat(n) { r = pso(fn, r) }
   }
   return r

}

fun mccormick(x: DoubleArray): Double { 

   val (a, b) = x
   return Math.sin(a + b) + (a - b) * (a - b) + 1.0 + 2.5 * b - 1.5 * a

}

fun michalewicz(x: DoubleArray): Double { 

   val m = 10
   val d = x.size
   var sum = 0.0
   for (i in 1..d) {
       val j = x[i - 1]
       val k = Math.sin(i * j * j / Math.PI)
       sum += Math.sin(j) * Math.pow(k, 2.0 * m)
   }
   return -sum

}

fun main(args: Array<String>) { 

   var state = psoInit(
       min = doubleArrayOf(-1.5, -3.0),
       max = doubleArrayOf(4.0, 4.0),
       parameters = Parameters(0.0, 0.6, 0.3),
       nParticles = 100
   )
   state = iterate(::mccormick, 40, state)
   state.report("McCormick")
   println("f(-.54719, -1.54719) : ${mccormick(doubleArrayOf(-.54719, -1.54719))}")
   println()
   state = psoInit(
       min = doubleArrayOf(0.0, 0.0),
       max = doubleArrayOf(Math.PI, Math.PI),
       parameters = Parameters(0.3, 0.3, 0.3),
       nParticles = 1000
   )
   state = iterate(::michalewicz, 30, state)
   state.report("Michalewicz (2D)")
   println("f(2.20, 1.57)        : ${michalewicz(doubleArrayOf(2.2, 1.57))}")

}</lang>

Sample output: 

Test Function        : McCormick
Iterations           : 40
Global Best Position : [-0.5471015946082899, -1.5471991634200966]
Global Best Value    : -1.913222941607108
f(-.54719, -1.54719) : -1.913222954882274

Test Function        : Michalewicz (2D)
Iterations           : 30
Global Best Position : [2.202908690715102, 1.5707970450218895]
Global Best Value    : -1.8013034099142804
f(2.20, 1.57)        : -1.801140718473825

ooRexx

<lang oorexx>/* REXX ---------------------------------------------------------------

  • Test for McCormick function
  • --------------------------------------------------------------------*/

Numeric Digits 16 Parse Value '-.5 -1.5 1' With x y d fmin=1e9 Call refine x,y Do r=1 To 10 

 d=d/5
 Call refine xmin,ymin
 End

Say 'which is better (less) than' Say ' f(-.54719,-1.54719)='f(-.54719,-1.54719) Say 'and differs from published -1.9133' Exit

refine: Parse Arg xx,yy Do x=xx-d To xx+d By d/2 

 Do y=yy-d To yy+d By d/2
   f=f(x,y)
   If f<fmin Then Do
     Say x y f
     fmin=f
     xmin=x
     ymin=y
     End
   End
 End

Return

f: Parse Arg x,y res=rxcalcsin(x+y,16,'R')+(x-y)**2-1.5*x+2.5*y+1 Return res

requires rxmath library</lang>
Output:
-1.5 -2.5 -1.243197504692072
-1.0 -2.0 -1.641120008059867
-0.5 -1.5 -1.909297426825682
-0.54 -1.54 -1.913132979507516
-0.548 -1.548 -1.913221840016527
-0.5480 -1.5472 -1.913222034492829
-0.5472 -1.5472 -1.913222954970650
-0.54720000 -1.54719872 -1.913222954973731
-0.54719872 -1.54719872 -1.913222954978670
-0.54719872 -1.54719744 -1.913222954978914
-0.54719744 -1.54719744 -1.913222954981015
-0.5471975424 -1.5471975424 -1.913222954981036
which is better (less) than
        f(-.54719,-1.54719)=-1.913222954882273
and differs from published  -1.9133

Racket

<lang racket>#lang racket/base (require racket/list racket/math)

(define (unbox-into-cycle s) 

 (if (box? s) (in-cycle (in-value (unbox s))) s))
Tries to "maximise" function > (so if you want a minimum, set #
> to <, IYSWIM)

(define (PSO f particles iterations hi lo #:ω ω #:φ_p φ_p #:φ_g φ_g #:> (>? >)) 

 (define dimensions (procedure-arity f))
 (unless (exact-nonnegative-integer? dimensions)
   (raise-argument-error 'PSO "function of fixed arity" 1 f))
 
 (define-values (x v)
   (for/lists (x v)
     ((_ particles))
     (for/lists (xi vi)
       ((d (in-range dimensions))
        (h (unbox-into-cycle hi))
        (l (unbox-into-cycle lo)))
       (define h-l (- h l))
       (values (+ l (* (random) h-l)) (+ (- h-l) (* 2 (random) h-l))))))
 
 (define (particle-step x_i v_i p_i g)
   (for/lists (x_i+ v_i+)
     ((x_id (in-list x_i))
      (v_id (in-list v_i))
      (p_id (in-list p_i))
      (g_d (in-list g)))
     (define v_id+ (+ (* ω v_id)
                      (* φ_p (random) (- p_id x_id))
                      (* φ_g (random) (- g_d x_id))))
     (values (+ x_id v_id+) v_id+)))
 
 (define (call-f args) (apply f args))
 (define g0 (argmax call-f x))
 (define-values (_X _V _P _P. G G.)
   (for/fold ; because of g and g., we can't use for/lists
    ((X x) (V v) (P x) (P. (map call-f x)) (g g0) (g. (apply f g0)))
    ((_ iterations))
     (for/fold
      ((x+ null) (v+ null) (p+ null) (p.+ null) (g+ g) (g.+ g.))                
      ((x_i (in-list X))
       (v_i (in-list V))
       (p_i (in-list P))
       (p._i (in-list P.)))        
       (define-values (x_i+ v_i+) (particle-step x_i v_i p_i g+))        
       (let* ((x._i+ (apply f x_i+))
              (new-p_i? (>? x._i+ p._i))
              (new-g? (>? x._i+ g.+)))
         (values (cons x_i+ x+)
                 (cons v_i+ v+)
                 (cons (if new-p_i? x_i+ p_i) p+)
                 (cons (if new-p_i? x._i+ p._i) p.+)
                 (if new-g? x_i+ g+)
                 (if new-g? x._i+ g.+))))))
 (values G G.))

(define (McCormick x1 x2) 

 (+ (sin (+ x1 x2)) (sqr (- x1 x2)) (* -1.5 x1) (* 2.5 x2) 1))

(define (Michalewitz d #:m (m 10)) 

 (define 2m (* 2 m))
 (define /pi (/ pi))
 (define (f . xx)
   (let Σ ((s 0) (i 1) (xx xx))
     (if (null? xx)
         (- s)
         (let ((x (car xx)))
           (Σ (+ s (* (sin x) (expt (sin (* i (sqr x) /pi)) 2m))) (+ i 1) (cdr xx))))))
 (procedure-reduce-arity f d))

(displayln "McCormick [-1.993] @ (-0.54719, -1.54719)") (PSO McCormick 1000 100 #(-1.5 -3) #(4 4) #:ω 0 #:φ_p 0.6 #:φ_g 0.3 #:> <) (displayln "Michalewitz 2d [-1.8013] @ (2.20, 1.57)") (PSO (Michalewitz 2) 1000 30 (box 0) (box pi) #:ω 0.3 #:φ_p 0.3 #:φ_g 0.3 #:> <) (displayln "Michalewitz 5d [-4.687658]") (PSO (Michalewitz 5) 1000 30 (box 0) (box pi) #:ω 0.3 #:φ_p 0.3 #:φ_g 0.3 #:> <) (displayln "Michalewitz 10d [-9.66015]") (PSO (Michalewitz 10) 1000 30 (box 0) (box pi) #:ω 0.3 #:φ_p 0.3 #:φ_g 0.3 #:> <)</lang>

Output:

Here is a sample run, the particles roll downhill quite nicely for McCormick, but there's a lot of space to search with the 10-dimensional Michalewitz; so YMMV with that one! 

McCormick [-1.993] @ (-0.54719, -1.54719)
'(-0.5471975539492846 -1.547197548223612)
-1.9132229549810367
Michalewitz 2d [-1.8013] @ (2.20, 1.57)
'(2.20290527060906 1.5707963523178217)
-1.8013034100975123
Michalewitz 5d [-4.687658]
'(2.188617053067511
  1.571283730996248
  1.2884975345181757
  1.9194689579781514
  1.7202092563763838)
-4.680722049442259
Michalewitz 10d [-9.66015]
'(1.359756739301337
  2.7216986742916007
  1.2823734619604734
  1.097509491839529
  2.2225042675789752
  0.9162856379217913
  1.8753760783453128
  0.7909979596555162
  0.46574677476493
  1.8558804696523914)
-6.432092623300999

REXX

Translation of: ooRexx

This REXX version uses a large   numeric digits   (but only displays 25 digits).

Classic REXX doesn't have a   sine   function, so a RYO version is included here.

The numeric precision is only limited to the number of decimal digits defined in the   pi   variable   (in this case,   118).

This REXX version supports the specifying of   X,   Y,   and   D,   as well as the number of particles,   and the number of decimal digits to be displayed.

The refinement loop is stopped when the function value stabilizes.

Note that REXX used decimal floating point, not binary. <lang rexx>/*REXX program calculates Particle Swarm Optimization as it migrates through a solution.*/ numeric digits length(pi()) - 1 /*sDigs: is the # of displayed digits.*/ parse arg x y d #part sDigs . /*obtain optional arguments from the CL*/ if x== | x=="," then x= -0.5 /*Not specified? Then use the default.*/ if y== | y=="," then y= -1.5 /* " " " " " " */ if d== | d=="," then d= 1 /* " " " " " " */ if #part== | #part=="," then #part=1e12 /* " " " " " " */ if sDigs== | sDigs=="," then sDigs= 25 /* " " " " " " */ minF=#part; old= /*number of particles is one billion. */ say center('X', sDigs+3, "═") center('Y', sDigs+3, "═") center('D', sDigs+3, "═") call refine x,y 

               do  until refine(minX,minY);   d=d*.2         /*increase the difference.*/
               end   /*until ···*/              /* [↑]  stop refining if no difference.*/

say indent=1 + (sDigs+3) * 2 /*compute the indentation for alignment*/ say right('The global minimum for f(-.54719, -1.54719) ───► ', indent) fmt(f(-.54719, -1.54719)) say right('The published global minimum is:' , indent) fmt( -1.9133 ) exit /*stick a fork in it, we're all done. */ /*──────────────────────────────────────────────────────────────────────────────────────*/ refine: parse arg xx,yy; h=d * .5 

                    do   x=xx-d  to xx+d  by h
                      do y=yy-d  to yy+d  by h;    f=f(x,y);   if f>=minF  then iterate
                      new=fmt(x) fmt(y) fmt(f);                if new=old  then return 1
                      say new;   minF=f;  minX=x;  minY=y;     old=new
                      end   /*y*/
                    end     /*x*/
       return 0

/*──────────────────────────────────────────────────────────────────────────────────one─liner subroutines───────────────────────────────*/ f: procedure: parse arg a,b; return sin(a+b) + (a-b)**2 - 1.5*a + 2.5*b + 1 fmt: parse arg ?;  ?=format(?,,sDigs); L=length(?); if pos(.,?)\==0 then ?=strip(strip(?,'T',0),"T",.); return left(?,L) pi: pi=3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282307; return pi r2r: return arg(1) // (pi()*2) /*normalize radians ───► a unit circle.*/ sin: procedure; parse arg x; x=r2r(x); numeric fuzz 5; z=x; _=x; q=x*x; do k=2 by 2 until p=z; p=z; _=-_*q/(k*(k+1)); z=z+_; end; return z</lang> output   when using the default inputs: 

═════════════X══════════════ ═════════════Y══════════════ ═════════════D══════════════
-1.5                         -2.5                         -1.2431975046920717486273609
-1                           -2                           -1.6411200080598672221007448
-0.5                         -1.5                         -1.9092974268256816953960199
-0.54                        -1.54                        -1.9131329795075164948766768
-0.548                       -1.548                       -1.9132218400165267634506035
-0.548                       -1.5472                      -1.9132220344928294065568196
-0.5472                      -1.5472                      -1.9132229549706499208388746
-0.5472                      -1.54719872                  -1.9132229549737311254290577
-0.54719872                  -1.54719872                  -1.9132229549786702369612333
-0.54719872                  -1.54719744                  -1.91322295497891365438682
-0.54719744                  -1.54719744                  -1.9132229549810149766572388
-0.5471975424                -1.5471975424                -1.9132229549810362588916172
-0.54719755264               -1.54719755264               -1.9132229549810363893093655
-0.547197550592              -1.547197550592              -1.9132229549810363922848065
-0.5471975514112             -1.5471975514112             -1.9132229549810363928381695
-0.5471975510016             -1.5471975510016             -1.9132229549810363928520779
-0.54719755116544            -1.54719755116544            -1.9132229549810363929162561
-0.547197551198208           -1.547197551198208           -1.9132229549810363929179331
-0.547197551198208           -1.54719755119755264         -1.9132229549810363929179344
-0.54719755119755264         -1.54719755119755264         -1.9132229549810363929179361
-0.54719755119755264         -1.54719755119689728         -1.9132229549810363929179365
-0.54719755119689728         -1.54719755119689728         -1.9132229549810363929179375
-0.54719755119689728         -1.547197551196766208        -1.9132229549810363929179375
-0.547197551196766208        -1.547197551196766208        -1.9132229549810363929179376
-0.547197551196766208        -1.547197551196635136        -1.9132229549810363929179376
-0.547197551196635136        -1.547197551196635136        -1.9132229549810363929179376
-0.547197551196635136        -1.5471975511966089216       -1.9132229549810363929179376
-0.5471975511966089216       -1.5471975511966089216       -1.9132229549810363929179376
-0.5471975511966089216       -1.54719755119660367872      -1.9132229549810363929179376
-0.54719755119660367872      -1.54719755119660367872      -1.9132229549810363929179376
-0.54719755119660367872      -1.54719755119659843584      -1.9132229549810363929179376
-0.54719755119659843584      -1.54719755119659843584      -1.9132229549810363929179376
-0.547197551196597387264     -1.547197551196597387264     -1.9132229549810363929179376
-0.5471975511965978066944    -1.5471975511965978066944    -1.9132229549810363929179376
-0.5471975511965978066944    -1.54719755119659776475136   -1.9132229549810363929179376
-0.54719755119659776475136   -1.54719755119659776475136   -1.9132229549810363929179376
-0.54719755119659776475136   -1.547197551196597756362752  -1.9132229549810363929179376
-0.547197551196597756362752  -1.547197551196597756362752  -1.9132229549810363929179376
-0.547197551196597756362752  -1.547197551196597747974144  -1.9132229549810363929179376
-0.547197551196597747974144  -1.547197551196597747974144  -1.9132229549810363929179376
-0.547197551196597747974144  -1.5471975511965977462964224 -1.9132229549810363929179376
-0.5471975511965977462964224 -1.5471975511965977462964224 -1.9132229549810363929179376
-0.5471975511965977462964224 -1.5471975511965977462293135 -1.9132229549810363929179376
-0.5471975511965977462293135 -1.5471975511965977462293135 -1.9132229549810363929179376
-0.5471975511965977462293135 -1.5471975511965977461622047 -1.9132229549810363929179376
-0.5471975511965977461622047 -1.5471975511965977461622047 -1.9132229549810363929179376
-0.5471975511965977461487829 -1.5471975511965977461487829 -1.9132229549810363929179376
-0.5471975511965977461541516 -1.5471975511965977461541516 -1.9132229549810363929179376
-0.547197551196597746154259  -1.547197551196597746154259  -1.9132229549810363929179376
-0.547197551196597746154259  -1.5471975511965977461542375 -1.9132229549810363929179376
-0.5471975511965977461542375 -1.5471975511965977461542375 -1.9132229549810363929179376
-0.5471975511965977461542375 -1.547197551196597746154216  -1.9132229549810363929179376
-0.547197551196597746154216  -1.547197551196597746154216  -1.9132229549810363929179376
-0.547197551196597746154216  -1.5471975511965977461542152 -1.9132229549810363929179376
-0.5471975511965977461542152 -1.5471975511965977461542152 -1.9132229549810363929179376
-0.5471975511965977461542152 -1.5471975511965977461542143 -1.9132229549810363929179376
-0.5471975511965977461542143 -1.5471975511965977461542143 -1.9132229549810363929179376
-0.5471975511965977461542145 -1.5471975511965977461542145 -1.9132229549810363929179376

      The global minimum for  f(-.54719, -1.54719)  ───►  -1.9132229548822735814541188
                         The published global minimum is: -1.9133

Output note:   the published global minimum (referenced above, as well as the function's arguments) can be found at:

  http://www.sfu.ca/~ssurjano/mccorm.html



Retrieved from "https://rosettacode.org/wiki/Particle_swarm_optimization?oldid=259775"