Particle swarm optimization

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.

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 acceptable resource usage. To provide a common basis for comparing implementations, the following test cases are recommended:

McCormick function - bowl-shaped, with a single minimum
- -1.5 < x1 < 4
- -3 < x2 < 4
- 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
- 0 < x1 < pi
- 0 < x2 < pi
- 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]]

C#

Translation of: java

<lang csharp>using System;

namespace ParticleSwarmOptimization {

   public struct Parameters {
       public double omega, phip, phig;

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

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

       public 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() {
           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;
       }

       public void Report(string testfunc) {
           Console.WriteLine("Test Function        : {0}", testfunc);
           Console.WriteLine("Iterations           : {0}", iter);
           Console.WriteLine("Global Best Position : {0}", string.Join(", ", gbpos));
           Console.WriteLine("Global Best Value    : {0}", gbval);
       }
   }

   class Program {
       public static State PsoInit(double[] min, double[] max, Parameters parameters, int nParticles) {
           var nDims = min.Length;
           double[][] pos = new double[nParticles][];
           for (int i = 0; i < nParticles; i++) {
               pos[i] = new double[min.Length];
               min.CopyTo(pos[i], 0);
           }
           double[][] vel = new double[nParticles][];
           for (int i = 0; i < nParticles; i++) {
               vel[i] = new double[nDims];
           }
           double[][] bpos = new double[nParticles][];
           for (int i = 0; i < nParticles; i++) {
               bpos[i] = new double[min.Length];
               min.CopyTo(bpos[i], 0);
           }
           double[] bval = new double[nParticles];
           for (int i = 0; i < nParticles; i++) {
               bval[i] = double.PositiveInfinity;
           }
           int iter = 0;
           double[] gbpos = new double[nDims];
           for (int i = 0; i < nDims; i++) {
               gbpos[i] = double.PositiveInfinity;
           }
           double gbval = double.PositiveInfinity;

           return new State(iter, gbpos, gbval, min, max, parameters, pos, vel, bpos, bval, nParticles, nDims);
       }

       static Random r = new Random();

       public static State Pso(Func<double[], double> fn, State y) {
           var 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] = new double[y.min.Length];
               y.min.CopyTo(bpos[i], 0);
           }
           double[] bval = new double[y.nParticles];
           double[] gbpos = new double[y.nDims];
           double gbval = double.PositiveInfinity;
           for (int j = 0; j < y.nParticles; j++) {
               // evaluate
               v[j] = fn.Invoke(y.pos[j]);
               // update
               if (v[j] < y.bval[j]) {
                   y.pos[j].CopyTo(bpos[j], 0);
                   bval[j] = v[j];
               }
               else {
                   y.bpos[j].CopyTo(bpos[j], 0);
                   bval[j] = y.bval[j];
               }
               if (bval[j] < gbval) {
                   gbval = bval[j];
                   bpos[j].CopyTo(gbpos, 0);
               }
           }
           double rg = r.NextDouble();
           double[][] pos = new double[y.nParticles][];
           double[][] vel = new double[y.nParticles][];
           for (int i = 0; i < y.nParticles; i++) {
               pos[i] = new double[y.nDims];
               vel[i] = new double[y.nDims];
           }
           for (int j = 0; j < y.nParticles; j++) {
               // migrate
               double rp = r.NextDouble();
               bool ok = true;
               for (int k = 0; k < y.nDims; k++) {
                   vel[j][k] = 0.0;
                   pos[j][k] = 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();
                   }
               }
           }
           var iter = 1 + y.iter;
           return new State(iter, gbpos, gbval, y.min, y.max, y.parameters, pos, vel, bpos, bval, y.nParticles, y.nDims);
       }

       public static State Iterate(Func<double[], double> fn, int n, State y) {
           State r = y;
           if (n == int.MaxValue) {
               State old = y;
               while (true) {
                   r = Pso(fn, r);
                   if (r.Equals(old)) break;
                   old = r;
               }
           }
           else {
               for (int i = 0; i < n; i++) {
                   r = Pso(fn, r);
               }
           }
           return r;
       }

       public static double Mccormick(double[] x) {
           var a = x[0];
           var b = x[1];
           return Math.Sin(a + b) + (a - b) * (a - b) + 1.0 + 2.5 * b - 1.5 * a;
       }

       public static double Michalewicz(double[] x) {
           int m = 10;
           int d = x.Length;
           double sum = 0.0;
           for (int i = 1; i < d; i++) {
               var j = x[i - 1];
               var k = Math.Sin(i * j * j / Math.PI);
               sum += Math.Sin(j) * Math.Pow(k, 2.0 * m);
           }
           return -sum;
       }

       static void Main(string[] args) {
           var 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(Mccormick, 40, state);
           state.Report("McCormick");
           Console.WriteLine("f(-.54719, -1.54719) : {0}", Mccormick(new double[] { -.54719, -1.54719 }));
           Console.WriteLine();

           state = PsoInit(
               new double[] { -0.0, -0.0 },
               new double[] { Math.PI, Math.PI },
               new Parameters(0.3, 0.3, 0.3),
               1000
               );
           state = Iterate(Michalewicz, 30, state);
           state.Report("Michalewicz (2D)");
           Console.WriteLine("f(2.20, 1.57)        : {0}", Michalewicz(new double[] { 2.20, 1.57 }));
       }
   }

}</lang>

Output:

Test Function        : McCormick
Iterations           : 40
Global Best Position : -0.546850526417689, -1.54649614884518
Global Best Value    : -1.91322235333426
f(-.54719, -1.54719) : -1.91322295488227

Test Function        : Michalewicz (2D)
Iterations           : 30
Global Best Position : 2.20290514143486, 2.20798457238775
Global Best Value    : -0.801303410096221
f(2.20, 1.57)        : -0.801166387820286

C++

Translation of: D

<lang cpp>#include <algorithm>

include <functional>
include <iostream>
include <random>
include <vector>

const auto PI = std::atan2(0, -1);

bool double_equals(double a, double b, double epsilon = 0.001) {

   return std::abs(a - b) < epsilon;

}

template <typename T> bool vector_equals(const std::vector<T> & lhs, const std::vector<T> & rhs) {

   if (lhs.size() != rhs.size()) {
       return false;
   }

   for (size_t i = 0; i < lhs.size(); i++) {
       if (!vector_equals(lhs[i], rhs[i])) {
           return false;
       }
   }

   return true;

}

template <typename T> bool vector_equals(const T & lhs, const T & rhs) {

   return lhs == rhs;

}

template <> bool vector_equals(const std::vector<double> & lhs, const std::vector<double> & rhs) {

   if (lhs.size() != rhs.size()) {
       return false;
   }

   for (size_t i = 0; i < lhs.size(); i++) {
       if (!double_equals(lhs[i], rhs[i])) {
           return false;
       }
   }

   return true;

}

template <typename T> std::ostream& operator<<(std::ostream & os, const std::vector<T> & v) {

   auto it = v.cbegin();
   auto end = v.cend();

   os << '[';
   if (it != end) {
       os << *it;
       it = std::next(it);
   }
   while (it != end) {
       os << ", " << *it;
       it = std::next(it);
   }
   return os << ']';

}

double uniform01() {

   static std::default_random_engine generator;
   static std::uniform_real_distribution<double> distribution(0.0, 1.0);
   return distribution(generator);

}

struct Parameters {

   double omega, phip, phig;

   bool operator==(const Parameters& rhs) {
       return double_equals(omega, rhs.omega)
           && double_equals(phip, rhs.phip)
           && double_equals(phig, rhs.phig);
   }

};

struct State {

   int iter;
   std::vector<double> gbpos;
   double gbval;
   std::vector<double> min;
   std::vector<double> max;
   Parameters parameters;
   std::vector<std::vector<double>> pos;
   std::vector<std::vector<double>> vel;
   std::vector<std::vector<double>> bpos;
   std::vector<double> bval;
   int nParticles;
   int nDims;

   bool operator==(const State& rhs) {
       return iter == rhs.iter
           && vector_equals(gbpos, rhs.gbpos)
           && double_equals(gbval, rhs.gbval)
           && vector_equals(min, rhs.min)
           && vector_equals(max, rhs.max)
           && parameters == rhs.parameters
           && vector_equals(pos, rhs.pos)
           && vector_equals(vel, rhs.vel)
           && vector_equals(bpos, rhs.bpos)
           && vector_equals(bval, rhs.bval)
           && nParticles == rhs.nParticles
           && nDims == rhs.nDims;
   }

   void report(const std::string& testFunc) {
       std::cout << "Test Function        : " << testFunc << '\n';
       std::cout << "Iterations           : " << iter << '\n';
       std::cout << "Global Best Position : " << gbpos << '\n';
       std::cout << "Global Best Value    : " << gbval << '\n';
   }

};

State psoInit(const std::vector<double> & min, const std::vector<double> & max, const Parameters & parameters, int nParticles) {

   int nDims = min.size();

   std::vector<std::vector<double>> pos(nParticles);
   for (int i = 0; i < nParticles; i++) {
       std::copy(min.cbegin(), min.cend(), std::back_inserter(pos[i]));
   }

   std::vector<std::vector<double>> vel(nParticles);
   for (int i = 0; i < nParticles; i++) {
       vel[i].resize(nDims);
   }

   std::vector<std::vector<double>> bpos(nParticles);
   for (int i = 0; i < nParticles; i++) {
       std::copy(min.cbegin(), min.cend(), std::back_inserter(bpos[i]));
   }

   std::vector<double> bval(nParticles, HUGE_VAL);

   auto iter = 0;

   std::vector<double> gbpos(nDims, HUGE_VAL);

   auto gbval = HUGE_VAL;

   return{ iter, gbpos, gbval, min, max, parameters, pos, vel, bpos, bval, nParticles, nDims };

}

State pso(const std::function<double(const std::vector<double>&)> & fn, const State & y) {

   auto p = y.parameters;

   std::vector<double> v(y.nParticles);

   std::vector<std::vector<double>> bpos(y.nParticles);
   for (int i = 0; i < y.nParticles; i++) {
       std::copy(y.min.cbegin(), y.min.cend(), std::back_inserter(bpos[i]));
   }

   std::vector<double> bval(y.nParticles);

   std::vector<double> gbpos(y.nDims);

   auto gbval = HUGE_VAL;

   for (int j = 0; j < y.nParticles; j++) {
       // 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();

   std::vector<std::vector<double>> pos(y.nParticles);
   for (size_t i = 0; i < pos.size(); i++) {
       pos[i].resize(y.nDims);
   }

   std::vector<std::vector<double>> vel(y.nParticles);
   for (size_t i = 0; i < vel.size(); i++) {
       vel[i].resize(y.nDims);
   }

   for (size_t j = 0; j < y.nParticles; j++) {
       // migrate
       auto rp = uniform01();
       bool ok = true;
       std::fill(vel[j].begin(), vel[j].end(), 0);
       std::fill(pos[j].begin(), pos[j].end(), 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]) * uniform01();
           }
       }
   }

   auto iter = 1 + y.iter;

   return { iter, gbpos, gbval, y.min, y.max, y.parameters, pos, vel, bpos, bval, y.nParticles, y.nDims };

}

State iterate(const std::function<double(const std::vector<double>&)> & fn, int n, const State & y) {

   State r(y);
   if (n == INT32_MAX) {
       State old(y);
       while (true) {
           r = pso(fn, r);
           if (r == old) {
               break;
           }
           old = r;
       }
   } else {
       for (int i = 0; i < n; i++) {
           r = pso(fn, r);
       }
   }
   return r;

}

double mccormick(const std::vector<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(const std::vector<double> & x) {

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

}

int main() {

   auto state = psoInit(
       { -1.5, -3.0 },
       { 4.0, 4.0 },
       { 0.0, 0.6, 0.3 },
       100
   );
   state = iterate(mccormick, 40, state);
   state.report("McCormick");
   std::cout << "f(-0.54719, -1.54719) : " << mccormick({ -0.54719, -1.54719 }) << '\n';
   std::cout << '\n';

   state = psoInit(
       { 0.0, 0.0 },
       {PI, PI},
       { 0.3, 0.3, 0.3 },
       1000
   );
   state = iterate(michalewicz, 30, state);
   state.report("Michalewicz (2D)");
   std::cout << "f(2.20, 1.57)        : " << michalewicz({ 2.2, 1.57 }) << '\n';

}</lang>

Output:

Test Function        : McCormick
Iterations           : 40
Global Best Position : [-0.547284, -1.54737]
Global Best Value    : -1.91322
f(-0.54719, -1.54719) : -1.91322

Test Function        : Michalewicz (2D)
Iterations           : 30
Global Best Position : [2.20291, 1.2939]
Global Best Value    : -0.801303
f(2.20, 1.57)        : -0.801166

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

Go

Translation of: Kotlin

<lang go>package main

import (

   "fmt"
   "math"
   "math/rand"
   "time"

)

type ff = func([]float64) float64

type parameters struct{ omega, phip, phig float64 }

type state struct {

   iter       int
   gbpos      []float64
   gbval      float64
   min        []float64
   max        []float64
   params     parameters
   pos        [][]float64
   vel        [][]float64
   bpos       [][]float64
   bval       []float64
   nParticles int
   nDims      int

}

func (s state) report(testfunc string) {

   fmt.Println("Test Function        :", testfunc)
   fmt.Println("Iterations           :", s.iter)
   fmt.Println("Global Best Position :", s.gbpos)
   fmt.Println("Global Best Value    :", s.gbval)

}

func psoInit(min, max []float64, params parameters, nParticles int) *state {

   nDims := len(min)
   pos := make([][]float64, nParticles)
   vel := make([][]float64, nParticles)
   bpos := make([][]float64, nParticles)
   bval := make([]float64, nParticles)
   for i := 0; i < nParticles; i++ {
       pos[i] = min
       vel[i] = make([]float64, nDims)
       bpos[i] = min
       bval[i] = math.Inf(1)
   }
   iter := 0
   gbpos := make([]float64, nDims)
   for i := 0; i < nDims; i++ {
       gbpos[i] = math.Inf(1)
   }
   gbval := math.Inf(1)
   return &state{iter, gbpos, gbval, min, max, params,
       pos, vel, bpos, bval, nParticles, nDims}

}

func pso(fn ff, y *state) *state {

   p := y.params
   v := make([]float64, y.nParticles)
   bpos := make([][]float64, y.nParticles)
   bval := make([]float64, y.nParticles)
   gbpos := make([]float64, y.nDims)
   gbval := math.Inf(1)
   for j := 0; j < y.nParticles; j++ {
       // 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]
       }
   }
   rg := rand.Float64()
   pos := make([][]float64, y.nParticles)
   vel := make([][]float64, y.nParticles)
   for j := 0; j < y.nParticles; j++ {
       pos[j] = make([]float64, y.nDims)
       vel[j] = make([]float64, y.nDims)
       // migrate
       rp := rand.Float64()
       ok := true
       for z := 0; z < y.nDims; z++ {
           pos[j][z] = 0
           vel[j][z] = 0
       }
       for 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 k := 0; k < y.nDims; k++ {
               pos[j][k] = y.min[k] + (y.max[k]-y.min[k])*rand.Float64()
           }
       }
   }
   iter := 1 + y.iter
   return &state{iter, gbpos, gbval, y.min, y.max, y.params,
       pos, vel, bpos, bval, y.nParticles, y.nDims}

}

func iterate(fn ff, n int, y *state) *state {

   r := y
   for i := 0; i < n; i++ {
       r = pso(fn, r)
   }
   return r

}

func mccormick(x []float64) float64 {

   a, b := x[0], x[1]
   return math.Sin(a+b) + (a-b)*(a-b) + 1.0 + 2.5*b - 1.5*a

}

func michalewicz(x []float64) float64 {

   m := 10.0
   sum := 0.0
   for i := 1; i <= len(x); i++ {
       j := x[i-1]
       k := math.Sin(float64(i) * j * j / math.Pi)
       sum += math.Sin(j) * math.Pow(k, 2*m)
   }
   return -sum

}

func main() {

   rand.Seed(time.Now().UnixNano())
   st := psoInit(
       []float64{-1.5, -3.0},
       []float64{4.0, 4.0},
       parameters{0.0, 0.6, 0.3},
       100,
   )
   st = iterate(mccormick, 40, st)
   st.report("McCormick")
   fmt.Println("f(-.54719, -1.54719) :", mccormick([]float64{-.54719, -1.54719}))
   fmt.Println()
   st = psoInit(
       []float64{0.0, 0.0},
       []float64{math.Pi, math.Pi},
       parameters{0.3, 0.3, 0.3},
       1000,
   )
   st = iterate(michalewicz, 30, st)
   st.report("Michalewicz (2D)")
   fmt.Println("f(2.20, 1.57)        :", michalewicz([]float64{2.2, 1.57}))

</lang>

Output:

Sample output:

Test Function        : McCormick
Iterations           : 40
Global Best Position : [-0.5473437041724806 -1.5464923165739348]
Global Best Value    : -1.9132220947578635
f(-.54719, -1.54719) : -1.913222954882274

Test Function        : Michalewicz (2D)
Iterations           : 30
Global Best Position : [2.2029051150895165 1.570796212894911]
Global Best Value    : -1.8013034100953598
f(2.20, 1.57)        : -1.8011407184738253

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>

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

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>

Julia

<lang julia>using Optim

const mcclow = [-1.5, -3.0] const mccupp = [4.0, 4.0] const miclow = [0.0, 0.0] const micupp = Float64.([pi, pi]) const npar = [100, 1000] const x0 = [0.0, 0.0]

michalewicz(x, m=10) = -sum(i -> sin(x[i]) * (i * sin( x[i]^2/pi))^(2*m), 1:length(x))

mccormick(x) = sin(x[1] + x[2]) + (x[1] - x[2])^2 - 1.5 * x[1] + 2.5 * x[2] + 1

println(optimize(mccormick, x0, ParticleSwarm(;lower=mcclow, upper=mccupp, n_particles=npar[1]))) @time optimize(mccormick, x0, ParticleSwarm(;lower=mcclow, upper=mccupp, n_particles=npar[1]))

println(optimize(michalewicz, x0, ParticleSwarm(;lower=miclow, upper=micupp, n_particles=npar[2]))) @time optimize(michalewicz, x0, ParticleSwarm(;lower=miclow, upper=micupp, n_particles=npar[2]))

</lang>

Output:

Results of Optimization Algorithm
 * Algorithm: Particle Swarm
 * Starting Point: [0.0,0.0]
 * Minimizer: [-0.5471975503990738,-1.5471975447742121]
 * Minimum: -1.913223e+00
 * Iterations: 1000
 * Convergence: false
   * |x - x'| ≤ 0.0e+00: false
     |x - x'| = NaN
   * |f(x) - f(x')| ≤ 0.0e+00 |f(x)|: false
     |f(x) - f(x')| = NaN |f(x)|
   * |g(x)| ≤ 1.0e-08: false
     |g(x)| = NaN
   * Stopped by an increasing objective: false
   * Reached Maximum Number of Iterations: true
 * Objective Calls: 101001
 * Gradient Calls: 0
0.087319 seconds (228.91 k allocations: 12.098 MiB, 59.41% gc time)

Results of Optimization Algorithm
 * Algorithm: Particle Swarm
 * Starting Point: [0.0,0.0]
 * Minimizer: [2.202905520771759,1.5707963264041795]
 * Minimum: -1.801303e+00
 * Iterations: 1000
 * Convergence: false
   * |x - x'| ≤ 0.0e+00: false
     |x - x'| = NaN
   * |f(x) - f(x')| ≤ 0.0e+00 |f(x)|: false
     |f(x) - f(x')| = NaN |f(x)|
   * |g(x)| ≤ 1.0e-08: false
     |g(x)| = NaN
   * Stopped by an increasing objective: false
   * Reached Maximum Number of Iterations: true
 * Objective Calls: 1001001
 * Gradient Calls: 0
2.312291 seconds (3.52 M allocations: 153.253 MiB, 0.49% gc time)

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

Nim

Translation of: Kotlin

<lang Nim>import math, random, sequtils, sugar

type

 Func = seq[float] -> float

 Parameters = tuple[omega, phip, phig: float]

 State = object
   iter: int
   gbpos: seq[float]
   gbval: float
   min, max: seq[float]
   parameters: Parameters
   pos, vel, bpos: seq[seq[float]]
   bval: seq[float]
   nParticles, nDims: int

func initState(min, max: seq[float]; parameters: Parameters; nParticles: int): State =

 let nDims = min.len
 State(iter: 0,
       gbpos: repeat(Inf, nDims),
       gbval: Inf,
       min: min,
       max: max,
       parameters: parameters,
       pos: repeat(min, nParticles),
       vel: newSeqWith(nParticles,
       newSeq[float](nDims)),
       bpos: repeat(min, nParticles),
       bval: repeat(Inf, nParticles),
       nParticles: nParticles,
       nDims: nDims)

proc report(state: State; testFunc: string) =

 echo "Test Function:        ", testfunc
 echo "Iterations:           ", state.iter
 echo "Global Best Position: ", state.gbpos
 echo "Global Best Value:    ", state.gbval

proc pso(fn: Func; y: State): State =

 let p = y.parameters
 var v = newSeq[float](y.nParticles)
 var bpos  = repeat(y.min, y.nParticles)
 var bval  = newSeq[float](y.nParticles)
 var gbpos = newSeq[float](y.nDims)
 var gbval = Inf

 for j in 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]

 let rg = rand(1.0)
 var pos = newSeqWith(y.nParticles, newSeq[float](y.nDims))
 var vel = newSeqWith(y.nParticles, newSeq[float](y.nDims))
 for j in 0..<y.nParticles:
   # migrate.
   let rp = rand(1.0)
   var ok = true
   for k in 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 and y.min[k] < pos[j][k] and y.max[k] > pos[j][k]
   if not ok:
     for k in 0..<y.nDims:
       pos[j][k]= y.min[k] + (y.max[k] - y.min[k]) * rand(1.0)

 result = State(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)

proc iterate(fn: Func; n: int; y: State): State =

 result = y
 if n == int.high:
   while true:
     let old = result
     result = pso(fn, result)
     if result == old: break
 else:
   for _ in 1..n:
     result = pso(fn, result)

func mccormick(x: seq[float]): float =

 let a = x[0]
 let b = x[1]
 result = sin(a + b) + (a - b) * (a - b) + 1.0 + 2.5 * b - 1.5 * a

func michalewicz(x: seq[float]): float =

 const M = 10
 for i in 1..x.len:
   let j = x[i - 1]
   let k = sin(i.toFloat * j * j / PI)
   result -= sin(j) * k^(2 * M)

randomize()

var state = initState(min = @[-1.5, -3],

                     max = @[4.0, 4.0],
                     parameters = (0.0, 0.6, 0.3),
                     nParticles = 100)

state = iterate(mccormick, 40, state) state.report("McCormick") echo "f(-.54719, -1.54719): ", mccormick(@[-0.54719, -1.54719]) echo() state = initState(min = @[0.0, 0.0],

                 max = @[PI, PI],
                 parameters = (0.3, 0.3, 0.3),
                 nParticles = 1000)

state = iterate(michalewicz, 30, state) state.report("Michalewicz (2D)") echo "f(2.20, 1.57): ", michalewicz(@[2.2, 1.57])</lang>

Output:

Test Function:        McCormick
Iterations:           40
Global Best Position: @[-0.5470347980396687, -1.547176688676891]
Global Best Value:    -1.913222920248667
f(-.54719, -1.54719): -1.913222954882274

Test Function:        Michalewicz (2D)
Iterations:           30
Global Best Position: @[2.202898715299719, 1.570804023976923]
Global Best Value:    -1.801303406946448
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

Perl

Translation of: Raku

<lang perl>use strict; use warnings; use feature 'say';

use constant PI => 2 * atan2(1, 0); use constant Inf => 1e10;

sub pso_init {

   my(%y) = @_;
   my $d = @{$y{'min'}};
   my $n = $y{'n'};

   $y{'gbval'} = Inf;
   $y{'gbpos'} = [(Inf) x $d];
   $y{'bval'}  = [(Inf) x $n];
   $y{'bpos'}  = [($y{'min'}) x $n];
   $y{'pos'}   = [($y{'min'}) x $n];
   $y{'vel'}   = [([(0) x $d]) x $n];

%y

}

sub pso {

   my($fn, %y) = @_;
   my $p      = $y{'p'};
   my $n      = $y{'n'};
   my $d      = @{$y{'min'}};
   my @bpos   = ($y{'min'}) x $n;
   my $gbval  = Inf;
   my $rand_g = rand;
   my (@pos, @vel, @v, @gbpos, @bval);

   for my $j (0 .. $n-1) {
       $v[$j] = &$fn(@{$y{'pos'}[$j]}); # evaluate

       # 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) {
           @gbpos = @{$bpos[$j]};
           $gbval = $bval[$j];
       }
   }

   # migrate
   for my $j (0 .. $n-1) {
       my $rand_p = rand;
       my $ok = 1;
       for my $k (0 .. $d-1) {
           $vel[$j][$k] = $$p{'omega'} * $y{'vel'}[$j][$k]
                      + $$p{'phi_p'} * $rand_p * ($bpos[$j][$k] - $y{'pos'}[$j][$k])
                      + $$p{'phi_g'} * $rand_g * ($gbpos[$k]    - $y{'pos'}[$j][$k]);
           $pos[$j][$k] = $y{'pos'}[$j][$k] + $vel[$j][$k];
           $ok = ($y{'min'}[$k] < $pos[$j][$k]) && ($pos[$j][$k] < $y{'max'}[$k]) && $ok;
       }
       next if $ok;
       $pos[$j][$_] = $y{'min'}[$_] + ($y{'max'}[$_] - $y{'min'}[$_]) * rand for 0 .. $d-1;
   }
   return {gbpos => \@gbpos, gbval => $gbval, bpos => \@bpos, bval => \@bval, pos => \@pos, vel => \@vel,
              min => $y{'min'}, max => $y{'max'}, p=> $y{'p'}, n => $n};

}

sub report {

   my($function_name, %state) = @_;
   say $function_name;
   say 'Global best position: ' . sprintf "%.5f %.5f", @{$state{'gbpos'}};
   say 'Global best value:    ' . sprintf "%.5f",      $state{'gbval'};
   say ;

}

McCormick

sub mccormick {

   my($a,$b) = @_;
   sin($a+$b) + ($a-$b)**2 + (1 + 2.5*$b - 1.5*$a)

}

my %state = pso_init( (

   min => [-1.5, -3],
   max => [4, 4],
   n   => 100,
   p   => {omega => 0, phi_p => 0.6, phi_g => 0.3},

) ); %state = %{pso(\&mccormick, %state)} for 1 .. 40; report('McCormick', %state);

Michalewicz

sub michalewicz {

   my(@x) = @_;
   my $sum;
   my $m = 10;
   for my $i (1..@x) {
       my $j = $x[$i - 1];
       my $k = sin($i * $j**2/PI);
       $sum += sin($j) * $k**(2*$m)
   }
   -$sum

}

%state = pso_init( (

   min => [0, 0],
   max => [PI, PI],
   n   => 1000,
   p   => {omega => 0.3, phi_p => 0.3, phi_g => 0.3},

) ); %state = %{pso(\&michalewicz, %state)} for 1 .. 30; report('Michalewicz', %state);</lang>

Output:

McCormick
Global best position: -0.54675 -1.54665
Global best value:    -1.91322

Michalewicz
Global best position: 2.20293 1.57080
Global best value:    -1.80130

Phix

Translation of: Kotlin

with javascript_semantics
enum OMEGA, PHIP, PHIG
enum ITER, GBPOS, GBVAL, MIN, MAX, PARAMS,
     POS, VEL, BPOS, BVAL, NPARTICLES, NDIMS 
 
constant inf = 1e308*1e308
 
constant fmt = """
Test Function        : %s
Iterations           : %d
Global Best Position : %s
Global Best Value    : %f
"""
 
procedure report(sequence state, string testfunc)
    printf(1,fmt,{testfunc,state[ITER],sprint(state[GBPOS]),state[GBVAL]})
end procedure
 
function psoInit(sequence mins, maxs, params, integer nParticles)
    integer nDims = length(mins), iter=0
    atom gbval = inf
    sequence gbpos = repeat(inf,nDims),
             pos = repeat(mins,nParticles),
             vel = repeat(repeat(0,nDims),nParticles),
             bpos = repeat(mins,nParticles),
             bval = repeat(inf,nParticles)
    return {iter,gbpos,gbval,mins,maxs,params,pos,vel,bpos,bval,nParticles,nDims}
end function
 
function pso(integer fn, sequence state)
    integer particles = state[NPARTICLES],
            dims = state[NDIMS]
    sequence p = state[PARAMS],
             v = repeat(0,particles),
             bpos = repeat(state[MIN],particles),
             bval = repeat(0,particles),
             gbpos = repeat(0,dims)
    atom gbval = inf
    for j=1 to particles do
        -- evaluate
        v[j] = fn(state[POS][j])
        -- update
        if v[j] < state[BVAL][j] then
            bpos[j] = state[POS][j]
            bval[j] = v[j]
        else
            bpos[j] = state[BPOS][j]
            bval[j] = state[BVAL][j]
        end if
        if bval[j] < gbval then
            gbval = bval[j]
            gbpos = bpos[j]
        end if
    end for
    atom rg = rnd()
    sequence pos = repeat(repeat(0,dims),particles),
             vel = repeat(repeat(0,dims),particles)
    for j=1 to particles do
        -- migrate
        atom rp = rnd()
        bool ok = true
        vel[j] = repeat(0,dims)
        pos[j] = repeat(0,dims)
        for k=1 to dims do
            vel[j][k] = p[OMEGA] * state[VEL][j][k] +
                        p[PHIP] * rp * (bpos[j][k] - state[POS][j][k]) +
                        p[PHIG] * rg * (gbpos[k] - state[POS][j][k])
            pos[j][k] = state[POS][j][k] + vel[j][k]
            ok = ok and state[MIN][k] < pos[j][k] and state[MAX][k] > pos[j][k]
        end for
        if not ok then
            for k=1 to dims do
                pos[j][k]= state[MIN][k] + (state[MAX][k] - state[MIN][k]) * rnd()
            end for
        end if
    end for
    integer iter = 1 + state[ITER]
    return {iter, gbpos, gbval, state[MIN], state[MAX], state[PARAMS],
            pos, vel, bpos, bval, particles, dims}
end function
 
function iterate(integer fn, n, sequence state)
    sequence r = state,
             old = state
    if n=-1 then
        while true do
            r = pso(fn, r)
            if (r == old) then exit end if
            old = r
        end while
    else
        for i=1 to n do
            r = pso(fn, r)
        end for
    end if
    return r
end function
 
function mccormick(sequence x)
    atom {a, b} = x
    return sin(a + b) + (a - b) * (a - b) + 1.0 + 2.5 * b - 1.5 * a
end function
 
function michalewicz(sequence x)
    integer m = 10,
            d = length(x)
    atom total = 0.0
    for i=1 to d do
        atom j = x[i],
             k = sin(i * j * j / PI)
        total += sin(j) * power(k, 2.0 * m)
    end for
    return -total
end function
 
procedure main()
    sequence mins = {-1.5, -3.0},
             maxs = {4.0, 4.0},
             params = {0.0, 0.6, 0.3}
    integer nParticles = 100
    sequence state = psoInit(mins,maxs,params,nParticles)
    state = iterate(mccormick, 40, state)
    report(state,"McCormick")
    atom {x,y} = state[GBPOS]
    printf(1,"f(%.4f, %.4f)  : %f\n\n",{x,y,mccormick({x,y})})
 
    mins = {0.0, 0.0}
    maxs = {PI, PI}
    params = {0.3, 0.3, 0.3}
    nParticles = 1000
    state = psoInit(mins,maxs,params,nParticles)
    state = iterate(michalewicz, 30, state)
    report(state,"Michalewicz (2D)")
    {x,y} = state[GBPOS]
    printf(1,"f(%.5f, %.5f)  : %f\n\n",{x,y,michalewicz({x,y})})
end procedure
main()

Output:

Test Function        : McCormick
Iterations           : 40
Global Best Position : {-0.5471808566,-1.547021879}
Global Best Value    : -1.913223
f(-0.5472, -1.5470)  : -1.913223

Test Function        : Michalewicz (2D)
Iterations           : 30
Global Best Position : {2.202905614,1.570796293}
Global Best Value    : -1.801303
f(2.20291, 1.57080)  : -1.801303

Python

Translation of: D

<lang python>import math import random

INFINITY = 1 << 127 MAX_INT = 1 << 31

class Parameters:

   def __init__(self, omega, phip, phig):
       self.omega = omega
       self.phip = phip
       self.phig = phig

class State:

   def __init__(self, iter, gbpos, gbval, min, max, parameters, pos, vel, bpos, bval, nParticles, nDims):
       self.iter = iter
       self.gbpos = gbpos
       self.gbval = gbval
       self.min = min
       self.max = max
       self.parameters = parameters
       self.pos = pos
       self.vel = vel
       self.bpos = bpos
       self.bval = bval
       self.nParticles = nParticles
       self.nDims = nDims

   def report(self, testfunc):
       print "Test Function :", testfunc
       print "Iterations    :", self.iter
       print "Global Best Position :", self.gbpos
       print "Global Best Value    : %.16f" % self.gbval

def uniform01():

   v = random.random()
   assert 0.0 <= v and v < 1.0
   return v

def psoInit(min, max, parameters, nParticles):

   nDims = len(min)
   pos = [min[:]] * nParticles
   vel = [[0.0] * nDims] * nParticles
   bpos = [min[:]] * nParticles
   bval = [INFINITY] * nParticles
   iter = 0
   gbpos = [INFINITY] * nDims
   gbval = INFINITY
   return State(iter, gbpos, gbval, min, max, parameters, pos, vel, bpos, bval, nParticles, nDims);

def pso(fn, y):

   p = y.parameters
   v = [0.0] * (y.nParticles)
   bpos = [y.min[:]] * (y.nParticles)
   bval = [0.0] * (y.nParticles)
   gbpos = [0.0] * (y.nDims)
   gbval = INFINITY
   for j in xrange(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][:]
   rg = uniform01()
   pos = [[None] * (y.nDims)] * (y.nParticles)
   vel = [[None] * (y.nDims)] * (y.nParticles)
   for j in xrange(0, y.nParticles):
       # migrate
       rp = uniform01()
       ok = True
       vel[j] = [0.0] * (len(vel[j]))
       pos[j] = [0.0] * (len(pos[j]))
       for k in xrange(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 and y.min[k] < pos[j][k] and y.max[k] > pos[j][k]
       if not ok:
           for k in xrange(0, y.nDims):
               pos[j][k] = y.min[k] + (y.max[k] - y.min[k]) * uniform01()
   iter = 1 + y.iter
   return State(iter, gbpos, gbval, y.min, y.max, y.parameters, pos, vel, bpos, bval, y.nParticles, y.nDims);

def iterate(fn, n, y):

   r = y
   old = y
   if n == MAX_INT:
       while True:
           r = pso(fn, r)
           if r == old:
               break
           old = r
   else:
       for _ in xrange(0, n):
           r = pso(fn, r)
   return r

def mccormick(x):

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

def michalewicz(x):

   m = 10
   d = len(x)
   sum = 0.0
   for i in xrange(1, d):
       j = x[i - 1]
       k = math.sin(i * j * j / math.pi)
       sum += math.sin(j) * k ** (2.0 * m)
   return -sum

def main():

   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")
   print "f(-.54719, -1.54719) : %.16f" % (mccormick([-.54719, -1.54719]))

   print

   state = psoInit([0.0, 0.0], [math.pi, math.pi], Parameters(0.3, 0.3, 0.3), 1000)
   state = iterate(michalewicz, 30, state)
   state.report("Michalewicz (2D)")
   print "f(2.20, 1.57)        : %.16f" % (michalewicz([2.2, 1.57]))

main()</lang>

Output:

Test Function : McCormick
Iterations    : 40
Global Best Position : [-0.5471069930124911, -1.5471582891466962]
Global Best Value    : -1.9132229450518705
f(-.54719, -1.54719) : -1.9132229548822739

Test Function : Michalewicz (2D)
Iterations    : 30
Global Best Position : [2.2029052187108036, 0.9404640520657541]
Global Best Value    : -0.8013034100970750
f(2.20, 1.57)        : -0.8011663878202856

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

Raku

(formerly Perl 6)

Translation of: J (via Javascript, Kotlin, D, Python)

<lang perl6>sub pso-init (%y) {

   my $d = @(%y{'min'});
   my $n = %y{'n'};

   %y{'gbval'} = Inf;
   %y{'gbpos'} = [Inf xx $d];
   %y{'bval'}  = [Inf xx $n];
   %y{'bpos'}  = [%y{'min'} xx $n];
   %y{'pos'}   = [%y{'min'} xx $n];
   %y{'vel'}   = [[0 xx $d] xx $n];

%y;

}

sub pso (&fn, %y) {

   my %p      = %y{'p'};
   my $n      = %y{'n'};
   my $d      = @(%y{'min'});
   my @bpos   = %y{'min'} xx $n;
   my $gbval  = Inf;
   my $rand-g = rand;
   my (@pos, @vel, @v, @gbpos, @bval);

   for 0 ..^ $n -> \j {
       @v[j] = &fn(%y{'pos'}[j]); # evaluate

       # 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];
       }
   }

   # migrate
   for 0 ..^ $n -> \j {
       my $rand-p = rand;
       my $ok = True;
       for 0 ..^ $d -> \k {
           @vel[j;k] = %p{'ω'} × %y{'vel'}[j;k]
                      + %p{'φ-p'} × $rand-p × (@bpos[j;k] - %y{'pos'}[j;k])
                      + %p{'φ-g'} × $rand-g × (@gbpos[k]  - %y{'pos'}[j;k]);
           @pos[j;k] = %y{'pos'}[j;k] + @vel[j;k];
           $ok = %y{'min'}[k] < @pos[j;k] < %y{'max'}[k] if $ok;
       }
       next if $ok;
       @pos[j;$_] = %y{'min'}[$_] + (%y{'max'}[$_] - %y{'min'}[$_]) × rand for 0 ..^ $d;
   }

   return {gbpos => @gbpos, gbval => $gbval, bpos => @bpos, bval => @bval, pos => @pos, vel => @vel,
           min => %y{'min'}, max => %y{'max'}, p=> %y{'p'}, n => $n}

}

sub report ($function-name, %state) {

   say $function-name;
   say '🌐 best position: ' ~ %state{'gbpos'}.fmt('%.5f');
   say '🌐 best value:    ' ~ %state{'gbval'}.fmt('%.5f');
   say ;

}

sub mccormick (@x) {

   my ($a,$b) = @x;
   sin($a+$b) + ($a-$b)**2 + (1 + 2.5×$b - 1.5×$a)

}

my %state = pso-init( {

   min => [-1.5, -3],
   max => [4, 4],
   n   => 100,
   p   => {ω=> 0, φ-p=> 0.6, φ-g=> 0.3},

} ); %state = pso(&mccormick, %state) for 1 .. 40; report 'McCormick', %state;

sub michalewicz (@x) {

   my $sum;
   my $m = 10;
   for 1..@x -> $i {
       my $j = @x[$i-1];
       my $k = sin($i × $j**2/π);
       $sum += sin($j) × $k**(2×$m)
   }
   -$sum

}

%state = pso-init( {

   min => [0, 0],
   max => [π, π],
   n   => 1000,
   p   => {ω=> 0.3, φ-p=> 0.3, φ-g=> 0.3},

} ); %state = pso(&michalewicz, %state) for 1 .. 30; report 'Michalewicz', %state;</lang>

Output:

McCormick
🌐 best position: -0.54714 -1.54710
🌐 best value:    -1.91322

Michalewicz
🌐 best position: 2.20291 1.57080
🌐 best value:    -1.80130

REXX

Translation of: ooRexx

This REXX version uses a large numeric digits (the number of decimal digits in pi), 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, 110).

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. A little extra code was added to show a title and align the output columns.

The refinement loop is stopped when the calculation of the function value stabilizes.

Note that REXX uses decimal floating point, not binary. <lang rexx>/*REXX program calculates Particle Swarm Optimization as it migrates through a solution.*/ numeric digits length( pi() ) - length(.) /*use the number of decimal digs in pi.*/ 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 /* " " " " " " */ old= /*#part: 1e12 ≡ is one trillion. */ minF= #part /*the minimum for the function (#part).*/ show= sDigs + 3 /*adjust number decimal digits for show*/ say "══iteration══" center('X',show,"═") center('Y',show,"═") center('D',show,"═")

= 0; call refine x, y /*#: REFINE iterations; invoke REFINE. */

               do  until refine(minX, minY)     /*perform until the mix is  "refined". */
               d=d  * .2                        /*decrease the difference in the mix. .*/
               end   /*until*/                  /* [↑]  stop refining if no difference.*/

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

                 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
                   #= # + 1;     say center(#,13) new            /*bump iter.; show new*/
                   minF= f;    minX= x;     minY= y;    old= new /*assign new values.  */
                   end   /*y*/
                 end     /*x*/
       return 0

/*──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────*/ f: procedure: parse arg a,b; return sin(a+b) + (a-b)**2 - 1.5*a + 2.5*b + 1 fmt: ?= format(arg(1), , sDigs); L= length(?); if pos(., ?)\==0 then ?= strip( strip(?, 'T', 0), "T", .); return left(?,L) pi: pi=3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865; return pi r2r: return arg(1) // ( pi() * 2) /*normalize radians ───► a unit circle.*/ sin: procedure; arg x; x= r2r(x); z=x; xx= x*x; do k=2 by 2 until p=z; p=z; x= -x* xx/ (k*(k+1)); z= z+x; end; return z</lang>

output when using the default input:

══iteration══ ═════════════X══════════════ ═════════════Y══════════════ ═════════════D══════════════
      1       -1.5                         -2.5                         -1.2431975046920717486273609
      2       -1                           -2                           -1.6411200080598672221007448
      3       -0.5                         -1.5                         -1.9092974268256816953960199
      4       -0.54                        -1.54                        -1.9131329795075164948766768
      5       -0.548                       -1.548                       -1.9132218400165267634506035
      6       -0.548                       -1.5472                      -1.9132220344928294065568196
      7       -0.5472                      -1.5472                      -1.9132229549706499208388746
      8       -0.5472                      -1.54719872                  -1.9132229549737311254290577
      9       -0.54719872                  -1.54719872                  -1.9132229549786702369612333
     10       -0.54719872                  -1.54719744                  -1.91322295497891365438682
     11       -0.54719744                  -1.54719744                  -1.9132229549810149766572388
     12       -0.5471975424                -1.5471975424                -1.9132229549810362588916172
     13       -0.54719755264               -1.54719755264               -1.9132229549810363893093655
     14       -0.547197550592              -1.547197550592              -1.9132229549810363922848065
     15       -0.5471975514112             -1.5471975514112             -1.9132229549810363928381695
     16       -0.5471975510016             -1.5471975510016             -1.9132229549810363928520779
     17       -0.54719755116544            -1.54719755116544            -1.9132229549810363929162561
     18       -0.547197551198208           -1.547197551198208           -1.9132229549810363929179331
     19       -0.547197551198208           -1.54719755119755264         -1.9132229549810363929179344
     20       -0.54719755119755264         -1.54719755119755264         -1.9132229549810363929179361
     21       -0.54719755119755264         -1.54719755119689728         -1.9132229549810363929179365
     22       -0.54719755119689728         -1.54719755119689728         -1.9132229549810363929179375
     23       -0.54719755119689728         -1.547197551196766208        -1.9132229549810363929179375
     24       -0.547197551196766208        -1.547197551196766208        -1.9132229549810363929179376
     25       -0.547197551196766208        -1.547197551196635136        -1.9132229549810363929179376
     26       -0.547197551196635136        -1.547197551196635136        -1.9132229549810363929179376
     27       -0.547197551196635136        -1.5471975511966089216       -1.9132229549810363929179376
     28       -0.5471975511966089216       -1.5471975511966089216       -1.9132229549810363929179376
     29       -0.5471975511966089216       -1.54719755119660367872      -1.9132229549810363929179376
     30       -0.54719755119660367872      -1.54719755119660367872      -1.9132229549810363929179376
     31       -0.54719755119660367872      -1.54719755119659843584      -1.9132229549810363929179376
     32       -0.54719755119659843584      -1.54719755119659843584      -1.9132229549810363929179376
     33       -0.547197551196597387264     -1.547197551196597387264     -1.9132229549810363929179376
     34       -0.5471975511965978066944    -1.5471975511965978066944    -1.9132229549810363929179376
     35       -0.5471975511965978066944    -1.54719755119659776475136   -1.9132229549810363929179376
     36       -0.54719755119659776475136   -1.54719755119659776475136   -1.9132229549810363929179376
     37       -0.54719755119659776475136   -1.547197551196597756362752  -1.9132229549810363929179376
     38       -0.547197551196597756362752  -1.547197551196597756362752  -1.9132229549810363929179376
     39       -0.547197551196597756362752  -1.547197551196597747974144  -1.9132229549810363929179376
     40       -0.547197551196597747974144  -1.547197551196597747974144  -1.9132229549810363929179376
     41       -0.547197551196597747974144  -1.5471975511965977462964224 -1.9132229549810363929179376
     42       -0.5471975511965977462964224 -1.5471975511965977462964224 -1.9132229549810363929179376
     43       -0.5471975511965977462964224 -1.5471975511965977462293135 -1.9132229549810363929179376
     44       -0.5471975511965977462293135 -1.5471975511965977462293135 -1.9132229549810363929179376
     45       -0.5471975511965977462293135 -1.5471975511965977461622047 -1.9132229549810363929179376
     46       -0.5471975511965977461622047 -1.5471975511965977461622047 -1.9132229549810363929179376
     47       -0.5471975511965977461487829 -1.5471975511965977461487829 -1.9132229549810363929179376
     48       -0.5471975511965977461541516 -1.5471975511965977461541516 -1.9132229549810363929179376
     49       -0.547197551196597746154259  -1.547197551196597746154259  -1.9132229549810363929179376
     50       -0.547197551196597746154259  -1.5471975511965977461542375 -1.9132229549810363929179376
     51       -0.5471975511965977461542375 -1.5471975511965977461542375 -1.9132229549810363929179376
     52       -0.5471975511965977461542375 -1.547197551196597746154216  -1.9132229549810363929179376
     53       -0.547197551196597746154216  -1.547197551196597746154216  -1.9132229549810363929179376
     54       -0.547197551196597746154216  -1.5471975511965977461542152 -1.9132229549810363929179376
     55       -0.5471975511965977461542152 -1.5471975511965977461542152 -1.9132229549810363929179376
     56       -0.5471975511965977461542152 -1.5471975511965977461542143 -1.9132229549810363929179376
     57       -0.5471975511965977461542143 -1.5471975511965977461542143 -1.9132229549810363929179376
     58       -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

Wren

Translation of: Kotlin

Library: Wren-dynamic

<lang ecmascript>import "random" for Random import "/dynamic" for Tuple

var Parameters = Tuple.create("Parameters", ["omega", "phip", "phig"])

var fields = [

   "iter", "gbpos", "gbval", "min", "max", "parameters",
   "pos", "vel", "bpos", "bval", "nParticles", "nDims"

]

var State = Tuple.create("State", fields)

var report = Fn.new { |state, testfunc|

   System.print("Test Function         : %(testfunc)")
   System.print("Iterations            : %(state.iter)")
   System.print("Global Best Position  : %(state.gbpos)")
   System.print("Global Best Value     : %(state.gbval)")

}

var psoInit = Fn.new { |min, max, parameters, nParticles|

   var nDims = min.count
   var pos   = List.filled(nParticles, null)
   var vel   = List.filled(nParticles, null)
   var bpos  = List.filled(nParticles, null)
   var bval  = List.filled(nParticles, 1/0)
   for (i in 0...nParticles) {
       pos[i]  = min.toList
       vel[i]  = List.filled(nDims, 0)
       bpos[i] = min.toList
   }
   var iter  = 0
   var gbpos = List.filled(nDims, 1/0 )
   var gbval = 1/0
   return State.new(iter, gbpos, gbval, min, max, parameters,
                    pos, vel, bpos, bval, nParticles, nDims)

}

var r = Random.new()

var pso = Fn.new { |fn, y|

   var p = y.parameters
   var v = List.filled(y.nParticles, 0)
   var bpos  = List.filled(y.nParticles, null)
   for (i in 0...y.nParticles) bpos[i] = y.min.toList
   var bval  = List.filled(y.nParticles, 0)
   var gbpos = List.filled(y.nDims, 0)
   var gbval = 1/0
   for (j in 0...y.nParticles) {
       // evaluate
       v[j] = fn.call(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]
       }
   }
   var rg = r.float()
   var pos = List.filled(y.nParticles, null)
   var vel = List.filled(y.nParticles, null)
   for (i in 0...y.nParticles) {
       pos[i] = List.filled(y.nDims, 0)
       vel[i] = List.filled(y.nDims, 0)
   }
   for (j in 0...y.nParticles) {
       // migrate
       var rp = r.float()
       var ok = true
       for (k in 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) {
           for (k in 0...y.nDims) {
               pos[j][k]= y.min[k] + (y.max[k] - y.min[k]) * r.float()
           }
       }
   }
   var iter = 1 + y.iter
   return State.new(
       iter, gbpos, gbval, y.min, y.max, y.parameters,
       pos, vel, bpos, bval, y.nParticles, y.nDims
   )

}

var iterate = Fn.new { |fn, n, y|

   var r = y
   var old = y
   if (n == 2147483647) {
       while (true) {
           r = pso.call(fn, r)
           if (r == old) break
           old = r
       }
   } else {
       for (i in 1..n) r = pso.call(fn, r)
   }
   return r

}

var mccormick = Fn.new { |x|

   var a = x[0]
   var b = x[1]
   return (a + b).sin + (a - b) * (a - b) + 1 + 2.5 * b - 1.5 * a

}

var michalewicz = Fn.new { |x|

   var m = 10
   var d = x.count
   var sum = 0
   for (i in 1..d) {
       var j = x[i - 1]
       var k = (i * j * j / Num.pi).sin
       sum = sum + j.sin * k.pow(2 * m)
   }
   return -sum

}

var state = psoInit.call([-1.5, -3], [4, 4], Parameters.new(0, 0.6, 0.3), 100) state = iterate.call(mccormick, 40, state) report.call(state, "McCormick") System.print("f(-0.54719, -1.54719) : %(mccormick.call([-0.54719, -1.54719]))") System.print() state = psoInit.call([0, 0], [Num.pi, Num.pi], Parameters.new(0.3, 0.3, 0.3), 1000) state = iterate.call(michalewicz, 30, state) report.call(state, "Michalewicz (2D)") System.print("f(2.20, 1.57) : %(michalewicz.call([2.2, 1.57]))")</lang>

Output:

Sample run:

Test Function         : McCormick
Iterations            : 40
Global Best Position  : [-0.54763537556709, -1.5469760587453]
Global Best Value     : -1.9132225000184
f(-0.54719, -1.54719) : -1.9132229548823

Test Function         : Michalewicz (2D)
Iterations            : 30
Global Best Position  : [2.2029075565418, 1.570796180786]
Global Best Value     : -1.8013034100303
f(2.20, 1.57)         : -1.8011407184738