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.
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 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
- 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:
Note this works with ALGOL 68 Genie versions 1, 2 and 3 on my Windows 11 laptop, however it runs out of memory on TIO.RUN, which is running (I believe) ALGOL 68 Genie version 2 under Linux. Recommend using version 3...
This problem appears to be data related as it occurs when the second optimisation is attempted. A possible work-around is indicated in the code.
BEGIN # particle swarm optimisation - translation of Go #
MODE FF = PROC([]REAL)REAL;
MODE PARAMETERS = STRUCT( REAL omega, phip, phig );
MODE STATE = STRUCT( INT iter
, REF[]REAL gbpos
, REAL gbval
, REF[]REAL min
, REF[]REAL max
, PARAMETERS params
, REF[,]REAL pos
, REF[,]REAL vel
, REF[,]REAL bpos
, REF[]REAL bval
, INT nparticles
, INT ndims
);
OP FMT = ( INT v )STRING: whole( v, 0 );
OP FMT = ( REAL v )STRING: fixed( v, -18, 15 );
OP FMT = ( []REAL v )STRING:
BEGIN
STRING result := "", prefix := "";
FOR i FROM LWB v TO UPB v DO
result +:= prefix + FMT v[ i ];
prefix := " "
OD;
"[" + result + "]"
END # FMT # ;
PROC println = ( STRING line )VOID: print( ( line, newline ) );
PRIO REPORT = 9;
OP REPORT = ( STATE s, STRING testfunc )VOID:
BEGIN
println( "Test Function : " + testfunc );
println( "Iterations : " + FMT iter OF s );
println( "Global Best Position : " + FMT gbpos OF s );
println( "Global Best Value : " + FMT gbval OF s )
END # REPORT # ;
PROC psoinit = ( []REAL min, max, PARAMETERS params, INT nparticles )REF STATE:
BEGIN
INT ndims = ( UPB min - LWB min ) + 1;
HEAP STATE result
:= STATE( # iter # 0
, # gbpos # HEAP[ 1 : ndims ]REAL
, # gbval # max real
, # min # HEAP[ 1 : ndims ]REAL := min[ AT 1 ]
, # max # HEAP[ 1 : ndims ]REAL := max[ AT 1 ]
, # params # params
, # pos # HEAP[ 1 : nparticles, 1 : ndims ]REAL
, # val # HEAP[ 1 : nparticles, 1 : ndims ]REAL
, # bpos # HEAP[ 1 : nparticles, 1 : ndims ]REAL
, # bval # HEAP[ 1 : nparticles ]REAL
, # nparticles # nparticles
, # ndims # ndims
);
FOR i TO nparticles DO
( pos OF result )[ i, : ] := min[ AT 1 ];
( bpos OF result )[ i, : ] := min[ AT 1 ];
( bval OF result )[ i ] := max real;
FOR j TO ndims DO ( vel OF result )[ i, j ] := 0 OD
OD;
FOR i TO ndims DO
( gbpos OF result )[ i ] := max real
OD;
result
END # psoinit # ;
PRIO PSO = 9;
OP PSO = ( REF STATE y, FF fn )REF STATE: # updates in-place the STATE y #
BEGIN
INT nparticles = nparticles OF y, ndims = ndims OF y;
[ 1 : nparticles ]REAL v;
PARAMETERS p = params OF y;
REAL omega = omega OF p, phip = phip OF p, phig = phig OF p;
REF[]REAL gbpos = gbpos OF y;
REF[]REAL min = min OF y;
REF[]REAL max = max OF y;
REF[,]REAL pos = pos OF y;
REF[,]REAL vel = vel OF y;
REF[,]REAL bpos = bpos OF y;
REF[]REAL bval = bval OF y;
FOR j TO nparticles DO
# evaluate #
v[ j ] := fn( pos[ j, : ] );
# update #
IF v[ j ] < bval[ j ] THEN
bpos[ j, : ] := pos[ j, : ]; # FOR TIO.RUN, replace this line with: #
# FOR k TO ndims DO bpos[ j, k ] := pos[ j, k ] OD; #
bval[ j ] := v[ j ]
FI;
IF bval[ j ] < gbval OF y THEN
gbval OF y := bval[ j ];
gbpos := bpos[ j, : ]
FI
OD;
REAL rg = random;
FOR j TO nparticles DO
# migrate #
REAL rp = random;
BOOL ok := TRUE;
FOR k TO ndims DO
vel[ j, k ] := omega * vel[ j, k ]
+ phip * rp * ( bpos[ j, k ] - pos[ j, k ] )
+ phig * rg * ( gbpos[ k ] - pos[ j, k ] );
pos[ j, k ] +:= vel[ j, k ];
IF ok THEN
ok := min[ k ] < pos[ j, k ] AND max[ k ] > pos[ j, k ]
FI
OD;
IF NOT ok THEN
FOR k TO ndims DO
pos[ j, k ] := min[ k ] + ( max[ k ] - min[ k ] ) * random
OD
FI
OD;
iter OF y +:= 1;
y
END # pso # ;
# performs n iterations of the STATE y, updating y in-place #
PROC iterate = ( FF fn, INT n, REF STATE y )REF STATE:
BEGIN
TO n DO y PSO fn OD;
y
END # iterate # ;
PROC mccormick = ( []REAL x )REAL:
BEGIN
REAL a = x[ LWB x ], b = x[ LWB x + 1 ];
sin( a + b ) + ( a - b ) * ( a - b ) + 1.0 + 2.5 * b - 1.5 * a
END # mccormick # ;
PROC michalewicz = ( []REAL x )REAL:
BEGIN
REAL sum := 0.0;
FOR i FROM LWB x TO UPB x DO
REAL j = x[ i ];
REAL k = sin( i * j * j / pi );
sum +:= sin( j ) * ( k ^ 20 )
OD;
- sum
END # michalewicz # ;
BEGIN
STATE st := psoinit( # min # ( -1.5, -3.0 )
, # max # ( 4.0, 4.0 )
, # params # ( 0.0, 0.6, 0.3 )
, # nparticles # 100
);
st := iterate( mccormick, 40, st );
st REPORT "McCormick";
println( "f(-.54719, -1.54719) : " + FMT mccormick( ( -0.54719, -1.54719 ) ) );
println( "" );
st := psoinit( # min # ( 0.0, 0.0 )
, # max # ( pi, pi )
, # params # ( 0.3, 0.3, 0.3 )
, # nparticles # 1000
);
st := iterate( michalewicz, 30, st );
st REPORT "Michalewicz (2D)";
println( "f(2.20, 1.57) : " + FMT michalewicz( ( 2.2, 1.57 ) ) )
END
END- Output:
Test Function : McCormick Iterations : 40 Global Best Position : [-0.547174623193133 -1.547119957495940] Global Best Value : -1.913222947617180 f(-.54719, -1.54719) : -1.913222954882270 Test Function : Michalewicz (2D) Iterations : 30 Global Best Position : [ 2.202913601885530 1.570798313997210] Global Best Value : -1.801303408877070 f(2.20, 1.57) : -1.801140718473820
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 }));
}
}
}
- 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
#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';
}
- 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
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;
ulong nParticles;
ulong 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, ulong 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]));
}
- 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
global target_fn omega phip phig iter xmin[] xmax[] pos[][] vel[][] bpos[][] bval[] gbpos[] gbval nparticles ndims .
#
func sinr x .
return sin (x * 180 / pi)
.
func evaluate x[] .
if target_fn = 1
# --- McCormick ---
a = x[1]
b = x[2]
s = sinr (a + b)
return s + (a - b) * (a - b) + 1 + 2.5 * b - 1.5 * a
elif target_fn = 2
# --- Michalewicz ---
m = 10
for i = 1 to len x[]
j = x[i]
k = sinr (i * j * j / pi)
sum += sinr j * pow k (2 * m)
.
return -sum
.
.
proc psoInit xmi[] xma[] o p g n .
xmin[] = xmi[]
xmax[] = xma[]
omega = o
phip = p
phig = g
nparticles = n
iter = 0
ndims = len xmin[]
gbval = 1e100
len gbpos[] ndims
for k = 1 to ndims : gbpos[k] = 1e100
len pos[][] nparticles
len vel[][] nparticles
len bpos[][] nparticles
len bval[] nparticles
for i = 1 to nparticles
len pos[i][] ndims
len vel[i][] ndims
len bpos[i][] ndims
for k = 1 to ndims
pos[i][k] = xmin[k] + (xmax[k] - xmin[k]) * randomf
vel[i][k] = 0
bpos[i][k] = pos[i][k]
.
bval[i] = 1e100
.
.
proc pso .
rg = randomf
for i = 1 to nparticles
v = evaluate pos[i][]
if v < bval[i]
bval[i] = v
bpos[i][] = pos[i][]
.
if bval[i] < gbval
gbval = bval[i]
gbpos[] = bpos[i][]
.
.
for i = 1 to nparticles
rp = randomf
for k = 1 to ndims
vel[i][k] = omega * vel[i][k] + phip * rp * (bpos[i][k] - pos[i][k]) + phig * rg * (gbpos[k] - pos[i][k])
pos[i][k] += vel[i][k]
if pos[i][k] <= xmin[k] or pos[i][k] >= xmax[k]
for k = 1 to ndims
pos[i][k] = xmin[k] + (xmax[k] - xmin[k]) * randomf
.
break 1
.
.
.
iter += 1
.
proc iterate n .
for i = 1 to n : pso
.
proc report name$ .
print "Test Function : " & name$
print "Iterations : " & iter
print "Global Best Position : " & gbpos[]
print "Global Best Value : " & gbval
.
numfmt 0 5
target_fn = 1
psoInit [ -1.5 -3.0 ] [ 4.0 4.0 ] 0.0 0.6 0.3 100
iterate 40
report "McCormick"
print evaluate [ -0.54719 -1.54719 ]
print ""
#
target_fn = 2
psoInit [ 0 0 ] [ pi pi ] 0.3 0.3 0.3 1000
iterate 30
report "Michalewicz (2D)"
print evaluate [ 2.2 1.57 ]- Output:
Test Function : McCormick Iterations : 40 Global Best Position : [ -0.54693 -1.54728 ] Global Best Value : -1.91322 -1.91322 Test Function : Michalewicz (2D) Iterations : 30 Global Best Position : [ 2.20291 1.57080 ] Global Best Value : -1.80130 -1.80114
! ==============================================================================
! Particle Swarm Optimization (PSO)
!
! PSO is a population-based metaheuristic that simulates the social behaviour
! of a flock of birds searching for food. Each candidate solution is called a
! "particle". Particles move through the search space, attracted both toward
! the best position they have personally visited (cognitive component) and
! toward the best position any particle in the swarm has ever visited (social
! component). No gradient information is required, so PSO works on
! discontinuous, non-differentiable, or black-box objective functions.
!
! Velocity update equation (applied per dimension j, per particle i):
!
! v_ij <- omega * v_ij
! + phi_p * r_p * (p_best_ij - x_ij) ! pull toward personal best
! + phi_g * r_g * (g_best_j - x_ij) ! pull toward global best
!
! x_ij <- x_ij + v_ij ! move particle
!
! where:
! omega - inertia weight: fraction of old velocity carried forward.
! High omega encourages global exploration; low omega (or zero)
! makes particles brake quickly and exploit known good regions.
! phi_p - cognitive coefficient: weight of the personal-best attraction.
! phi_g - social coefficient: weight of the global-best attraction.
! r_p, r_g - independent uniform random scalars in [0,1], redrawn each step.
!
! Compile: gfortran -O2 -o pso pso.f90
! ==============================================================================
module pso_module
use iso_fortran_env, only: real64
implicit none
! --------------------------------------------------------------------------
! Abstract interface that every objective function must satisfy.
! x - position vector of length n (the point being evaluated)
! n - number of dimensions
! f - scalar objective value; PSO minimises this quantity
!
! Declaring the interface here lets the compiler type-check any function
! passed to the pso() subroutine at compile time rather than silently
! accepting a mismatched argument.
! --------------------------------------------------------------------------
abstract interface
function obj_func_t(x, n) result(f)
import :: real64
integer, intent(in) :: n
real(kind=real64), intent(in) :: x(n)
real(kind=real64) :: f
end function obj_func_t
end interface
contains
! ==========================================================================
! subroutine pso
!
! Runs the Particle Swarm Optimization algorithm and returns the best
! position and value found.
!
! Arguments:
! n_dim - number of search-space dimensions
! n_particles - swarm size (more particles = broader initial coverage)
! n_iter - number of full swarm update cycles (iterations)
! omega - inertia weight (see header notes)
! phi_p - cognitive coefficient (personal-best attraction strength)
! phi_g - social coefficient (global-best attraction strength)
! lb(n_dim) - lower bounds for each dimension
! ub(n_dim) - upper bounds for each dimension
! obj_func - the objective function to minimise (procedure argument)
! best_pos - (out) position of the global best found
! best_val - (out) objective value at best_pos
! ==========================================================================
subroutine pso(n_dim, n_particles, n_iter, omega, phi_p, phi_g, &
lb, ub, obj_func, best_pos, best_val)
integer, intent(in) :: n_dim, n_particles, n_iter
real(kind=real64), intent(in) :: omega, phi_p, phi_g
real(kind=real64), intent(in) :: lb(n_dim), ub(n_dim)
real(kind=real64), intent(out) :: best_pos(n_dim), best_val
procedure(obj_func_t) :: obj_func
! pos(i,j) - current position of particle i in dimension j
! vel(i,j) - current velocity of particle i in dimension j
! p_best_pos(i,j)- best position particle i has ever personally visited
! p_best_val(i) - objective value at p_best_pos(i,:)
! g_best_pos(j) - best position any particle has ever visited (global)
! g_best_val - objective value at g_best_pos
real(kind=real64), allocatable :: pos(:,:), vel(:,:)
real(kind=real64), allocatable :: p_best_pos(:,:), p_best_val(:)
real(kind=real64) :: g_best_pos(n_dim), g_best_val
real(kind=real64) :: r_p ! random scalar for personal-best term
real(kind=real64) :: r_g ! random scalar for global-best term
real(kind=real64) :: fval ! objective value of current position
real(kind=real64) :: v_max ! per-dimension velocity clamp magnitude
integer :: i, j, k ! particle index, dimension index, iteration
! Heap-allocate swarm arrays so large swarms do not overflow the stack.
allocate(pos(n_particles, n_dim))
allocate(vel(n_particles, n_dim))
allocate(p_best_pos(n_particles, n_dim))
allocate(p_best_val(n_particles))
! Seed the RNG. Calling random_seed() with no argument lets the runtime
! choose a seed (often from the system clock), giving a different sequence
! each run. Replace with random_seed(put=...) for reproducibility.
call random_seed()
! ----------------------------------------------------------------------
! Phase 1: Initialisation
! Scatter particles uniformly at random within [lb, ub].
! Assign each particle a small random initial velocity (+/-10% of the
! dimension range) so they are not immediately stationary.
! Evaluate the objective at each starting position and record it as the
! particle's personal best (since it is the only position seen so far).
! ----------------------------------------------------------------------
do i = 1, n_particles
do j = 1, n_dim
! Uniform random position: lb + U[0,1] * (ub - lb)
call random_number(r_p)
pos(i, j) = lb(j) + r_p * (ub(j) - lb(j))
! Small random initial velocity centred on zero: U[-0.1, 0.1]*range
! r_p is in [0,1], so (r_p - 0.5) is in [-0.5, 0.5];
! multiplying by 0.2 gives [-0.1, 0.1] of the dimension range.
call random_number(r_p)
vel(i, j) = (ub(j) - lb(j)) * (r_p - 0.5_real64) * 0.2_real64
end do
! The initial position is trivially the personal best.
p_best_pos(i, :) = pos(i, :)
p_best_val(i) = obj_func(pos(i, :), n_dim)
end do
! ----------------------------------------------------------------------
! Phase 2: Identify the initial global best across the whole swarm.
! Start from particle 1, then scan the rest for anything better.
! ----------------------------------------------------------------------
g_best_val = p_best_val(1)
g_best_pos = p_best_pos(1, :)
do i = 2, n_particles
if (p_best_val(i) < g_best_val) then
g_best_val = p_best_val(i)
g_best_pos = p_best_pos(i, :)
end if
end do
! ----------------------------------------------------------------------
! Phase 3: Main PSO iteration loop
!
! Each iteration: for every particle in every dimension,
! (a) draw fresh random scalars r_p and r_g,
! (b) update velocity using the standard PSO formula,
! (c) clamp velocity so particles cannot travel more than v_max per
! step -- without this, particles with zero inertia (omega=0) can
! receive arbitrarily large kicks from phi_p and phi_g and overshoot
! the entire search space,
! (d) advance position and clip it to [lb, ub] if it strays out,
! (e) evaluate the objective at the new position,
! (f) update personal best and, if better still, global best.
! ----------------------------------------------------------------------
do k = 1, n_iter
do i = 1, n_particles
do j = 1, n_dim
! Fresh independent random scalars for this particle-dimension-
! iteration triplet. Using separate r_p and r_g breaks any
! correlation between the cognitive and social pulls.
call random_number(r_p)
call random_number(r_g)
! Velocity update: inertia + cognitive pull + social pull
vel(i, j) = omega * vel(i, j) &
+ phi_p * r_p * (p_best_pos(i, j) - pos(i, j)) &
+ phi_g * r_g * (g_best_pos(j) - pos(i, j))
! Velocity clamping: limit to +/- 20% of the dimension range.
! This prevents "velocity explosion" -- a failure mode where
! particles rocket out of the search space and never return,
! especially common when omega=0 removes damping.
v_max = 0.2_real64 * (ub(j) - lb(j))
vel(i, j) = max(-v_max, min(v_max, vel(i, j)))
! Advance position by the new velocity.
pos(i, j) = pos(i, j) + vel(i, j)
! Hard boundary enforcement: if the particle crosses a bound,
! clamp it to the boundary edge. The velocity is not reversed
! here; it will be redirected by the attraction terms on the
! next iteration once the particle is stuck on the wall.
pos(i, j) = max(lb(j), min(ub(j), pos(i, j)))
end do ! end dimension loop
! Evaluate the objective at the updated position.
fval = obj_func(pos(i, :), n_dim)
! Update personal best if this position is an improvement.
if (fval < p_best_val(i)) then
p_best_val(i) = fval
p_best_pos(i, :) = pos(i, :)
! A new personal best can only improve the global best,
! so check here rather than in a separate sweep.
if (fval < g_best_val) then
g_best_val = fval
g_best_pos = pos(i, :)
end if
end if
end do ! end particle loop
end do ! end iteration loop
! Return the global best found to the caller.
best_pos = g_best_pos
best_val = g_best_val
deallocate(pos, vel, p_best_pos, p_best_val)
end subroutine pso
end module pso_module
! ==============================================================================
! Standard benchmark objective functions used to exercise the PSO solver.
! Both conform to obj_func_t so they can be passed directly to pso().
! ==============================================================================
module test_functions
use iso_fortran_env, only: real64
implicit none
! pi is needed by Michalewicz; compute at compile time via acos identity.
real(kind=real64), parameter :: pi = acos(-1.0_real64)
contains
! --------------------------------------------------------------------------
! McCormick function (2-dimensional)
!
! f(x1, x2) = sin(x1 + x2) + (x1 - x2)^2 - 1.5*x1 + 2.5*x2 + 1
!
! The surface is bowl-shaped with a single global minimum and no
! significant local minima, making it a useful sanity-check: a correctly
! implemented PSO should always find it.
!
! Recommended search bounds: -1.5 <= x1 <= 4, -3 <= x2 <= 4
! Global minimum: f(-0.54719, -1.54719) = -1.9133285
! --------------------------------------------------------------------------
function mccormick(x, n) result(f)
integer, intent(in) :: n
real(kind=real64), intent(in) :: x(n)
real(kind=real64) :: f
! sin(x1+x2) gives the undulating ridge character of the surface.
! (x1-x2)^2 forms a parabolic valley running diagonally.
! The linear terms -1.5*x1 + 2.5*x2 tilt the surface and shift
! the minimum away from the origin.
f = sin(x(1) + x(2)) + (x(1) - x(2))**2 &
- 1.5_real64*x(1) + 2.5_real64*x(2) + 1.0_real64
end function mccormick
! --------------------------------------------------------------------------
! Michalewicz function (n-dimensional, parameter m controls steepness)
!
! f(x) = -sum_{i=1}^{n} sin(xi) * [ sin(i * xi^2 / pi) ]^(2*m)
!
! m=10 (used here) produces extremely steep ridges and deep, narrow
! valleys. The function has n! local minima, of which n are deep enough
! to trap a naive optimizer, making it a demanding multimodal test case.
!
! Recommended search bounds: 0 <= xi <= pi for each dimension.
!
! 2-D global minimum: f(2.2029, 1.5708) approx -1.8013
! x2 = pi/2 exactly (the peak of sin) -- the ridge is infinitely
! thin in the x2 direction when m is large, which is why a large
! swarm (1000 particles) is needed to land close enough to it.
! --------------------------------------------------------------------------
function michalewicz(x, n) result(f)
integer, intent(in) :: n
real(kind=real64), intent(in) :: x(n)
real(kind=real64) :: f
integer :: i
integer, parameter :: m = 10 ! steepness parameter; higher = sharper ridges
f = 0.0_real64
do i = 1, n
! Each term is negative (subtracted from zero) so the function is
! minimised where the absolute value of the sum is maximised.
! sin(xi) peaks at xi = pi/2; the power term (2*m=20) sharpens the
! valley around each i-dependent xi* location so that almost all the
! function value is concentrated in a tiny neighbourhood of xi*.
f = f - sin(x(i)) * (sin(real(i, real64) * x(i)**2 / pi))**(2*m)
end do
end function michalewicz
end module test_functions
! ==============================================================================
! Main program: run the two benchmark tests and report results.
! ==============================================================================
program particle_swarm
use iso_fortran_env, only: real64
use pso_module
use test_functions
implicit none
real(kind=real64) :: best_pos(2) ! 2-D position of the optimum found
real(kind=real64) :: best_val ! objective value at best_pos
real(kind=real64) :: lb(2), ub(2) ! search-space bounds for each test
print *, '================================================'
print *, ' Particle Swarm Optimization'
print *, '================================================'
print *
! ==========================================================================
! Test 1: McCormick function
!
! omega = 0 means particles carry NO inertia from the previous step: each
! velocity update is purely attraction-based. This is aggressive exploitation
! and works well for a single-basin function like McCormick where there is
! no danger of prematurely converging to a local minimum.
!
! phi_p = 0.6 (cognitive) > phi_g = 0.3 (social): particles trust their own
! history more than the swarm consensus, keeping them spread across the bowl
! until they individually close in on the minimum.
!
! 100 particles over 40 iterations is sufficient for a smooth 2-D landscape.
! ==========================================================================
print *, 'Test 1: McCormick function'
print *, ' f(x1,x2) = sin(x1+x2) + (x1-x2)^2 - 1.5*x1 + 2.5*x2 + 1'
print *, ' Bounds: -1.5 < x1 < 4, -3 < x2 < 4'
print *, ' omega=0.0 phi_p=0.6 phi_g=0.3 N=100 iter=40'
print *
lb = [-1.5_real64, -3.0_real64]
ub = [ 4.0_real64, 4.0_real64]
!subroutine pso(n_dim, n_particles, n_iter, omega, phi_p, phi_g, &
! lb, ub, obj_func, best_pos, best_val)
call pso(2, 500, 100, 0.0_real64, 0.6_real64, 0.3_real64, &
lb, ub, mccormick, best_pos, best_val)
print '(2x,a,f12.7)', 'Best x1 = ', best_pos(1)
print '(2x,a,f12.7)', 'Best x2 = ', best_pos(2)
print '(2x,a,f12.7)', 'Best f(x) = ', best_val
print *, ' Known minimum: f(-0.54719, -1.54719) = -1.9133285'
print *
! ==========================================================================
! Test 2: Michalewicz function
!
! omega = 0.3 retains some inertia, helping particles glide across the
! narrow inter-valley ridges rather than stalling on a slope.
!
! phi_p = phi_g = 0.3: balanced cognitive and social terms keep the swarm
! from collapsing onto a single discovered valley too quickly, which would
! miss the true global minimum hidden in a different (possibly narrower) one.
!
! 1000 particles are needed because the deep valleys at m=10 are so narrow
! that with fewer particles most runs fail to place any particle close enough
! to the correct valley for it to be discovered. 30 iterations suffice once
! the right valley is seeded, since the social pull then rapidly draws the
! whole swarm in.
! ==========================================================================
print *, 'Test 2: Michalewicz function (m=10)'
print *, ' f(x) = -sum_i sin(xi) * [sin(i*xi^2/pi)]^20'
print *, ' Bounds: 0 < x1 < pi, 0 < x2 < pi'
print *, ' omega=0.3 phi_p=0.3 phi_g=0.3 N=1000 iter=30'
print *
lb = [0.0_real64, 0.0_real64]
ub = [pi, pi]
call pso(2, 1000, 30, 0.3_real64, 0.3_real64, 0.3_real64, &
lb, ub, michalewicz, best_pos, best_val)
print '(2x,a,f12.7)', 'Best x1 = ', best_pos(1)
print '(2x,a,f12.7)', 'Best x2 = ', best_pos(2)
print '(2x,a,f12.7)', 'Best f(x) = ', best_val
print *, ' Known minimum: f(2.2029, 1.5708) approx -1.8013'
print *
print *, '================================================'
end program particle_swarm
================================================ Particle Swarm Optimization ================================================ Test 1: McCormick function f(x1,x2) = sin(x1+x2) + (x1-x2)^2 - 1.5*x1 + 2.5*x2 + 1 Bounds: -1.5 < x1 < 4, -3 < x2 < 4 omega=0.0 phi_p=0.6 phi_g=0.3 N=100 iter=40 Best x1 = -0.5471975 Best x2 = -1.5471976 Best f(x) = -1.9132230 Known minimum: f(-0.54719, -1.54719) = -1.9133285 Test 2: Michalewicz function (m=10) f(x) = -sum_i sin(xi) * [sin(i*xi^2/pi)]^20 Bounds: 0 < x1 < pi, 0 < x2 < pi omega=0.3 phi_p=0.3 phi_g=0.3 N=1000 iter=30 Best x1 = 2.2029054 Best x2 = 1.5707966 Best f(x) = -1.8013034 Known minimum: f(2.2029, 1.5708) approx -1.8013 ================================================
#define INF 1E+308
Type Parameters
omega As Double
phip As Double
phig As Double
End Type
Type State
iter As Integer
gbpos(Any) As Double
gbval As Double
minBound(Any) As Double
maxBound(Any) As Double
params As Parameters
posic(Any, Any) As Double
vel(Any, Any) As Double
bpos(Any, Any) As Double
bval(Any) As Double
nParticles As Integer
nDims As Integer
End Type
Function mccormick(x() As Double) As Double
Dim As Double a = x(0)
Dim As Double b = x(1)
Return Sin(a + b) + (a - b) * (a - b) + 1 + 2.5 * b - 1.5 * a
End Function
Function michalewicz(x() As Double) As Double
Const m As Integer = 10
Dim As Integer d = Ubound(x) + 1
Dim As Double sum = 0
For i As Integer = 1 To d
Dim As Double j = x(i - 1)
Dim As Double k = Sin(i * j * j / (4 * Atn(1)))
sum += Sin(j) * (k ^ (2 * m))
Next
Return -sum
End Function
Sub initState(Byref state As State, minBound() As Double, maxBound() As Double, params As Parameters, nParticles As Integer)
state.nDims = Ubound(minBound) + 1
state.nParticles = nParticles
Redim state.gbpos(state.nDims - 1)
Redim state.minBound(state.nDims - 1)
Redim state.maxBound(state.nDims - 1)
Redim state.posic(state.nParticles - 1, state.nDims - 1)
Redim state.vel(state.nParticles - 1, state.nDims - 1)
Redim state.bpos(state.nParticles - 1, state.nDims - 1)
Redim state.bval(state.nParticles - 1)
state.iter = 0
state.gbval = INF
state.params = params
For i As Integer = 0 To state.nDims - 1
state.minBound(i) = minBound(i)
state.maxBound(i) = maxBound(i)
state.gbpos(i) = INF
Next
Randomize Timer
For i As Integer = 0 To state.nParticles - 1
For j As Integer = 0 To state.nDims - 1
state.posic(i, j) = minBound(j)
state.vel(i, j) = 0
state.bpos(i, j) = minBound(j)
Next
state.bval(i) = INF
Next
End Sub
Sub psoIteration(Byref state As State, fn As Function(x() As Double) As Double)
Dim As Double v, rg = Rnd
For j As Integer = 0 To state.nParticles - 1
Dim As Double tempPos(state.nDims - 1)
For k As Integer = 0 To state.nDims - 1
tempPos(k) = state.posic(j, k)
Next
v = fn(tempPos())
If v < state.bval(j) Then
For k As Integer = 0 To state.nDims - 1
state.bpos(j, k) = state.posic(j, k)
Next
state.bval(j) = v
If v < state.gbval Then
state.gbval = v
For k As Integer = 0 To state.nDims - 1
state.gbpos(k) = state.bpos(j, k)
Next
End If
End If
Next
For j As Integer = 0 To state.nParticles - 1
Dim As Double rp = Rnd
Dim As Boolean ok = True
For k As Integer = 0 To state.nDims - 1
state.vel(j, k) = state.params.omega * state.vel(j, k) + _
state.params.phip * rp * (state.bpos(j, k) - state.posic(j, k)) + _
state.params.phig * rg * (state.gbpos(k) - state.posic(j, k))
state.posic(j, k) += state.vel(j, k)
If state.posic(j, k) <= state.minBound(k) Or state.posic(j, k) >= state.maxBound(k) Then
ok = False
End If
Next
If Not ok Then
For k As Integer = 0 To state.nDims - 1
state.posic(j, k) = state.minBound(k) + (state.maxBound(k) - state.minBound(k)) * Rnd
Next
End If
Next
state.iter += 1
End Sub
Sub report(state As State, testfunc As String)
Print Using !"Test Function : &\nIterations : &\nGlobal Best Position : [&, &]\nGlobal Best Value : &"; testfunc; state.iter; state.gbpos(0); state.gbpos(1); state.gbval
End Sub
' Main program
Dim As State state
Dim As Parameters params = Type(0.0, 0.6, 0.3)
Dim As Double minBound(1) = {-1.5, -3.0}
Dim As Double maxBound(1) = {4.0, 4.0}
initState(state, minBound(), maxBound(), params, 100)
For i As Integer = 1 To 40
psoIteration(state, @mccormick)
Next
report(state, "McCormick")
Dim As Double testArray(1) = {-0.54719, -1.54719}
Print "f(-0.54719, -1.54719) : " & mccormick(testArray())
Dim As Parameters params2 = Type(0.3, 0.3, 0.3)
Dim As Double minBound2(1) = {0.0, 0.0}
Dim As Double maxBound2(1) = {3.14159265359, 3.14159265359}
initState(state, minBound2(), maxBound2(), params2, 1000)
For i As Integer = 1 To 30
psoIteration(state, @michalewicz)
Next
Print
report(state, "Michalewicz (2D)")
Dim As Double testArray2(1) = {2.2, 1.57}
Print "f(2.20, 1.57) : " & michalewicz(testArray2())
Sleep
- Output:
Test Function : McCormick Iterations : 40 Global Best Position : [-0.547543880423963, -1.546859553291497] Global Best Value : -1.913222486647357 f(-0.54719, -1.54719) : -1.913222954882274 Test Function : Michalewicz (2D) Iterations : 30 Global Best Position : [2.202852310793197, 1.570815305240421] Global Best Value : -1.801303349496491 f(2.20, 1.57) : -1.801140718473825
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}))
}
- 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
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&{.
Apply to McCormick Function:
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
Apply to Michalewicz Function:
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
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}));
}
}
- 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
Translation of J.
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);
Example displayed result (random numbers are involved so there will be a bit of variance between repeated runs:
Iteration: 0
GlobalBestPosition: Infinity
GlobalBestValue: Infinity
Iteration: 40
GlobalBestPosition: -0.5134004259016365,-1.5512442672625184
GlobalBestValue: -1.9114053788600853
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]))
- 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)
// 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))}")
}
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
(* Particle Swarm Optimization in Mathematica *)
(* Constants *)
INFINITY = 10^100;
MAXINT = 2^31;
(* Parameters structure *)
Parameters[omega_, phip_, phig_] := <|
"omega" -> omega,
"phip" -> phip,
"phig" -> phig
|>;
(* State structure *)
State[iter_, gbpos_, gbval_, min_, max_, parameters_, pos_, vel_, bpos_, bval_, nParticles_, nDims_] := <|
"iter" -> iter,
"gbpos" -> gbpos,
"gbval" -> gbval,
"min" -> min,
"max" -> max,
"parameters" -> parameters,
"pos" -> pos,
"vel" -> vel,
"bpos" -> bpos,
"bval" -> bval,
"nParticles" -> nParticles,
"nDims" -> nDims
|>;
(* Report function *)
report[state_, testfunc_] := (
Print["Test Function : ", testfunc];
Print["Iterations : ", state["iter"]];
Print["Global Best Position : ", state["gbpos"]];
Print["Global Best Value : ", NumberForm[state["gbval"], 16]];
);
(* Random number generator (0 to 1) *)
uniform01[] := RandomReal[];
(* PSO initialization *)
psoInit[min_, max_, parameters_, nParticles_] := Module[
{nDims, pos, vel, bpos, bval, iter, gbpos, gbval},
nDims = Length[min];
pos = Table[min, {nParticles}];
vel = Table[ConstantArray[0.0, nDims], {nParticles}];
bpos = Table[min, {nParticles}];
bval = ConstantArray[INFINITY, nParticles];
iter = 0;
gbpos = ConstantArray[INFINITY, nDims];
gbval = INFINITY;
State[iter, gbpos, gbval, min, max, parameters, pos, vel, bpos, bval, nParticles, nDims]
];
(* Main PSO algorithm *)
pso[fn_, y_] := Module[
{p, v, bpos, bval, gbpos, gbval, rg, pos, vel, j, k, rp, ok},
p = y["parameters"];
v = ConstantArray[0.0, y["nParticles"]];
bpos = Table[y["min"], {y["nParticles"]}];
bval = ConstantArray[0.0, y["nParticles"]];
gbpos = ConstantArray[0.0, y["nDims"]];
gbval = INFINITY;
(* Evaluate and update best positions *)
For[j = 1, j <= y["nParticles"], j++,
(* Evaluate *)
v[[j]] = fn[y["pos"][[j]]];
(* Update personal best *)
If[v[[j]] < y["bval"][[j]],
bpos[[j]] = y["pos"][[j]];
bval[[j]] = v[[j]],
bpos[[j]] = y["bpos"][[j]];
bval[[j]] = y["bval"][[j]]
];
(* Update global best *)
If[bval[[j]] < gbval,
gbval = bval[[j]];
gbpos = bpos[[j]]
];
];
rg = uniform01[];
pos = Table[ConstantArray[0.0, y["nDims"]], {y["nParticles"]}];
vel = Table[ConstantArray[0.0, y["nDims"]], {y["nParticles"]}];
(* Update velocities and positions *)
For[j = 1, j <= y["nParticles"], j++,
rp = uniform01[];
ok = True;
For[k = 1, 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]];
];
(* Reset position if out of bounds *)
If[!ok,
For[k = 1, k <= y["nDims"], k++,
pos[[j, k]] = y["min"][[k]] + (y["max"][[k]] - y["min"][[k]]) * uniform01[];
];
];
];
iter = 1 + y["iter"];
State[iter, gbpos, gbval, y["min"], y["max"], y["parameters"], pos, vel, bpos, bval, y["nParticles"], y["nDims"]]
];
(* Iteration function *)
iterate[fn_, n_, y_] := Module[
{r, old},
r = y;
old = y;
If[n == MAXINT,
(* Iterate until convergence *)
While[True,
r = pso[fn, r];
If[r === old, Break[]];
old = r;
],
(* Fixed number of iterations *)
For[i = 1, i <= n, i++,
r = pso[fn, r];
];
];
r
];
(* Test functions *)
mccormick[x_] := Module[{a, b},
{a, b} = x;
Sin[a + b] + (a - b)^2 + 1.0 + 2.5 b - 1.5 a
];
michalewicz[x_] := Module[{m, d, sum, i, j, k},
m = 10;
d = Length[x];
sum = 0.0;
For[i = 1, i <= d, i++,
j = x[[i]];
k = Sin[i * j^2 / Pi];
sum += Sin[j] * k^(2.0 * m);
];
-sum
];
(* Main execution *)
main[] := (
(* Test McCormick function *)
state = psoInit[{-1.5, -3.0}, {4.0, 4.0}, Parameters[0.0, 0.6, 0.3], 100];
state = iterate[mccormick, 40, state];
report[state, "McCormick"];
Print["f(-.54719, -1.54719) : ", NumberForm[mccormick[{-0.54719, -1.54719}], 16]];
Print[];
(* Test Michalewicz function *)
state = psoInit[{0.0, 0.0}, {Pi, Pi}, Parameters[0.3, 0.3, 0.3], 1000];
state = iterate[michalewicz, 30, state];
report[state, "Michalewicz (2D)"];
Print["f(2.20, 1.57) : ", NumberForm[michalewicz[{2.2, 1.57}], 16]];
);
(* Run the main function *)
main[]
- Output:
Test Function : McCormick
Iterations : 40
Global Best Position : {-0.5471097794914574, -1.54684236508515}
Global Best Value : NumberForm[-1.9132227985012031, 16]
f(-.54719, -1.54719) : NumberForm[-1.913222954882274, 16]
Test Function : Michalewicz (2D)
Iterations : 30
Global Best Position : {2.202900184402327, 1.570795445580627}
Global Best Value : NumberForm[-1.8013034096043792, 16]
f(2.20, 1.57) : NumberForm[-1.801140718473825, 16]
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])
- 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
/* 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
- 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
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);
- 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
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
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 range(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 range(y.nParticles):
# migrate
rp = uniform01()
ok = True
vel[j] = [0.0] * (len(vel[j]))
pos[j] = [0.0] * (len(pos[j]))
for k in range(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 range(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 range(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 range(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()
- 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
#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 #:> <)
- 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
(formerly Perl 6)
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;
- Output:
McCormick 🌐 best position: -0.54714 -1.54710 🌐 best value: -1.91322 Michalewicz 🌐 best position: 2.20291 1.57080 🌐 best value: -1.80130
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.
/*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
- 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:
use std::f64::consts::PI;
use rand::prelude::*;
const EPSILON: f64 = 0.001;
fn double_equals(a: f64, b: f64) -> bool {
(a - b).abs() < EPSILON
}
fn vector_equals_f64(lhs: &[f64], rhs: &[f64]) -> bool {
lhs.len() == rhs.len() && lhs.iter().zip(rhs).all(|(a, b)| double_equals(*a, *b))
}
fn vector_equals_2d(lhs: &[Vec<f64>], rhs: &[Vec<f64>]) -> bool {
lhs.len() == rhs.len() && lhs.iter().zip(rhs).all(|(a, b)| vector_equals_f64(a, b))
}
#[derive(Debug, Clone)]
struct Parameters {
omega: f64,
phip: f64,
phig: f64,
}
impl PartialEq for Parameters {
fn eq(&self, other: &Self) -> bool {
double_equals(self.omega, other.omega)
&& double_equals(self.phip, other.phip)
&& double_equals(self.phig, other.phig)
}
}
#[derive(Debug, Clone)]
struct State {
iter: i32,
gbpos: Vec<f64>,
gbval: f64,
min: Vec<f64>,
max: Vec<f64>,
parameters: Parameters,
pos: Vec<Vec<f64>>,
vel: Vec<Vec<f64>>,
bpos: Vec<Vec<f64>>,
bval: Vec<f64>,
n_particles: usize,
n_dims: usize,
}
impl PartialEq for State {
fn eq(&self, other: &Self) -> bool {
self.iter == other.iter
&& vector_equals_f64(&self.gbpos, &other.gbpos)
&& double_equals(self.gbval, other.gbval)
&& vector_equals_f64(&self.min, &other.min)
&& vector_equals_f64(&self.max, &other.max)
&& self.parameters == other.parameters
&& vector_equals_2d(&self.pos, &other.pos)
&& vector_equals_2d(&self.vel, &other.vel)
&& vector_equals_2d(&self.bpos, &other.bpos)
&& vector_equals_f64(&self.bval, &other.bval)
&& self.n_particles == other.n_particles
&& self.n_dims == other.n_dims
}
}
impl State {
fn report(&self, test_func: &str) {
println!("Test Function : {}", test_func);
println!("Iterations : {}", self.iter);
println!("Global Best Position : {:?}", self.gbpos);
println!("Global Best Value : {}", self.gbval);
}
}
fn pso_init(min: Vec<f64>, max: Vec<f64>, parameters: Parameters, n_particles: usize) -> State {
let n_dims = min.len();
let pos = vec![min.clone(); n_particles];
let vel = vec![vec![0.0; n_dims]; n_particles];
let bpos = vec![min.clone(); n_particles];
let bval = vec![f64::INFINITY; n_particles];
let iter = 0;
let gbpos = vec![f64::INFINITY; n_dims];
let gbval = f64::INFINITY;
State {
iter,
gbpos,
gbval,
min,
max,
parameters,
pos,
vel,
bpos,
bval,
n_particles,
n_dims,
}
}
fn pso<F>(func: F, state: &State) -> State
where
F: Fn(&[f64]) -> f64,
{
let mut rng = rand::rng();
let p = &state.parameters;
let mut v = vec![0.0; state.n_particles];
let mut bpos = vec![state.min.clone(); state.n_particles];
let mut bval = vec![0.0; state.n_particles];
let mut gbpos = vec![0.0; state.n_dims];
let mut gbval = f64::INFINITY;
// Evaluate and update best positions
for j in 0..state.n_particles {
v[j] = func(&state.pos[j]);
if v[j] < state.bval[j] {
bpos[j] = state.pos[j].clone();
bval[j] = v[j];
} else {
bpos[j] = state.bpos[j].clone();
bval[j] = state.bval[j];
}
if bval[j] < gbval {
gbval = bval[j];
gbpos = bpos[j].clone();
}
}
let rg = rng.random::<f64>();
let mut pos = vec![vec![0.0; state.n_dims]; state.n_particles];
let mut vel = vec![vec![0.0; state.n_dims]; state.n_particles];
// Update particle positions and velocities
for j in 0..state.n_particles {
let rp = rng.random::<f64>();
let mut ok = true;
for k in 0..state.n_dims {
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 && state.min[k] < pos[j][k] && state.max[k] > pos[j][k];
}
if !ok {
for k in 0..state.n_dims {
pos[j][k] = state.min[k] + (state.max[k] - state.min[k]) * rng.random::<f64>();
}
}
}
let iter = state.iter + 1;
State {
iter,
gbpos,
gbval,
min: state.min.clone(),
max: state.max.clone(),
parameters: state.parameters.clone(),
pos,
vel,
bpos,
bval,
n_particles: state.n_particles,
n_dims: state.n_dims,
}
}
fn iterate<F>(func: F, n: Option<i32>, state: &State) -> State
where
F: Fn(&[f64]) -> f64,
{
let mut result = state.clone();
match n {
None => {
// Iterate until convergence
loop {
let old = result.clone();
result = pso(&func, &result);
if result == old {
break;
}
}
}
Some(iterations) => {
for _ in 0..iterations {
result = pso(&func, &result);
}
}
}
result
}
fn mccormick(x: &[f64]) -> f64 {
let a = x[0];
let b = x[1];
(a + b).sin() + (a - b).powi(2) + 1.0 + 2.5 * b - 1.5 * a
}
fn michalewicz(x: &[f64]) -> f64 {
let m = 10;
let d = x.len();
let mut sum = 0.0;
for i in 1..d {
let j = x[i - 1];
let k = ((i as f64) * j * j / PI).sin();
sum += j.sin() * k.powf(2.0 * m as f64);
}
-sum
}
fn main() {
// Test McCormick function
let mut state = pso_init(
vec![-1.5, -3.0],
vec![4.0, 4.0],
Parameters {
omega: 0.0,
phip: 0.6,
phig: 0.3,
},
100,
);
state = iterate(mccormick, Some(40), &state);
state.report("McCormick");
println!("f(-0.54719, -1.54719) : {}", mccormick(&[-0.54719, -1.54719]));
println!();
// Test Michalewicz function
state = pso_init(
vec![0.0, 0.0],
vec![PI, PI],
Parameters {
omega: 0.3,
phip: 0.3,
phig: 0.3,
},
1000,
);
state = iterate(michalewicz, Some(30), &state);
state.report("Michalewicz (2D)");
println!("f(2.20, 1.57) : {}", michalewicz(&[2.2, 1.57]));
}
- Output:
Test Function : McCormick Iterations : 40 Global Best Position : [-0.5476448555592675, -1.5475984929762063] Global Best Value : -1.9132226413204705 f(-0.54719, -1.54719) : -1.913222954882274 Test Function : Michalewicz (2D) Iterations : 30 Global Best Position : [2.2029078402505253, 2.53387403303021] Global Best Value : -0.801303410011067 f(2.20, 1.57) : -0.8011663878202854
import scala.util.Random
import scala.math._
object PSO {
case class Parameters(omega: Double, phip: Double, phig: Double)
case class State(
iter: Int,
gbpos: Array[Double],
gbval: Double,
min: Array[Double],
max: Array[Double],
parameters: Parameters,
pos: Array[Array[Double]],
vel: Array[Array[Double]],
bpos: Array[Array[Double]],
bval: Array[Double],
nParticles: Int,
nDims: Int
) {
def report(testFunc: String): Unit = {
println(s"Test Function : $testFunc")
println(s"Iterations : $iter")
println(s"Global Best Position : ${gbpos.mkString("[", ", ", "]")}")
println(f"Global Best value : $gbval%.15f")
}
}
private val random = new Random()
def psoInit(min: Array[Double], max: Array[Double], parameters: Parameters, nParticles: Int): State = {
val nDims = min.length
val pos = Array.fill(nParticles)(min.clone())
val vel = Array.fill(nParticles, nDims)(0.0)
val bpos = Array.fill(nParticles)(min.clone())
val bval = Array.fill(nParticles)(Double.PositiveInfinity)
val iter = 0
val gbpos = Array.fill(nDims)(Double.PositiveInfinity)
val gbval = Double.PositiveInfinity
State(iter, gbpos, gbval, min, max, parameters, pos, vel, bpos, bval, nParticles, nDims)
}
def pso(fn: Array[Double] => Double, y: State): State = {
val p = y.parameters
val v = Array.ofDim[Double](y.nParticles)
val bpos = Array.fill(y.nParticles)(y.min.clone())
val bval = Array.ofDim[Double](y.nParticles)
var gbpos = Array.ofDim[Double](y.nDims)
var gbval = Double.PositiveInfinity
// Evaluate and update best positions
for (j <- 0 until y.nParticles) {
v(j) = fn(y.pos(j))
if (v(j) < y.bval(j)) {
bpos(j) = y.pos(j).clone()
bval(j) = v(j)
} else {
bpos(j) = y.bpos(j).clone()
bval(j) = y.bval(j)
}
if (bval(j) < gbval) {
gbval = bval(j)
gbpos = bpos(j).clone()
}
}
val rg = random.nextDouble()
val pos = Array.ofDim[Double](y.nParticles, y.nDims)
val vel = Array.ofDim[Double](y.nParticles, y.nDims)
// Update particle positions and velocities
for (j <- 0 until y.nParticles) {
val rp = random.nextDouble()
var ok = true
for (k <- 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 <- 0 until y.nDims) {
pos(j)(k) = y.min(k) + (y.max(k) - y.min(k)) * random.nextDouble()
}
}
}
val iter = y.iter + 1
State(iter, gbpos, gbval, y.min, y.max, y.parameters, pos, vel, bpos, bval, y.nParticles, y.nDims)
}
def iterate(fn: Array[Double] => Double, n: Int, initialState: State): State = {
if (n == Int.MaxValue) {
// Infinite iteration until convergence
var current = initialState
var previous = initialState
do {
previous = current
current = pso(fn, current)
} while (current != previous)
current
} else {
// Fixed number of iterations
(0 until n).foldLeft(initialState)((state, _) => pso(fn, state))
}
}
def mccormick(x: Array[Double]): Double = {
val a = x(0)
val b = x(1)
sin(a + b) + (a - b) * (a - b) + 1.0 + 2.5 * b - 1.5 * a
}
def michalewicz(x: Array[Double]): Double = {
val m = 10
val d = x.length
val sum = (1 until d).map { i =>
val j = x(i - 1)
val k = sin(i * j * j / Pi)
sin(j) * pow(k, 2.0 * m)
}.sum
-sum
}
def main(args: Array[String]): Unit = {
// Test McCormick function
var state = psoInit(
Array(-1.5, -3.0),
Array(4.0, 4.0),
Parameters(0.0, 0.6, 0.3),
100
)
state = iterate(mccormick, 40, state)
state.report("McCormick")
println(f"f(-.54719, -1.54719) : ${mccormick(Array(-0.54719, -1.54719))}%.15f")
println()
// Test Michalewicz function
state = psoInit(
Array(0.0, 0.0),
Array(Pi, Pi),
Parameters(0.3, 3.0, 0.3),
1000
)
state = iterate(michalewicz, 30, state)
state.report("Michalewicz (2D)")
println(f"f(2.20, 1.57) : ${michalewicz(Array(2.20, 1.57))}%.15f")
}
}
- Output:
Test Function : McCormick Iterations : 40 Global Best Position : [-0.5475499889655079, -1.547524418209413] Global Best value : -1.913222754537433 f(-.54719, -1.54719) : -1.913222954882274 Test Function : Michalewicz (2D) Iterations : 30 Global Best Position : [2.2029055203309, 0.7707578361091618] Global Best value : -0.801303410098552 f(2.20, 1.57) : -0.801166387820286
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]))")
- 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
const std = @import("std");
const math = std.math;
const print = std.debug.print;
const ArrayList = std.ArrayList;
const Allocator = std.mem.Allocator;
const Random = std.Random;
const EPSILON: f64 = 0.001;
fn doubleEquals(a: f64, b: f64) bool {
return @abs(a - b) < EPSILON;
}
fn vectorEqualsF64(lhs: []const f64, rhs: []const f64) bool {
if (lhs.len != rhs.len) return false;
for (lhs, rhs) |a, b| {
if (!doubleEquals(a, b)) return false;
}
return true;
}
fn vectorEquals2D(lhs: []const []const f64, rhs: []const []const f64) bool {
if (lhs.len != rhs.len) return false;
for (lhs, rhs) |a, b| {
if (!vectorEqualsF64(a, b)) return false;
}
return true;
}
const Parameters = struct {
omega: f64,
phip: f64,
phig: f64,
fn equals(self: Parameters, other: Parameters) bool {
return doubleEquals(self.omega, other.omega) and
doubleEquals(self.phip, other.phip) and
doubleEquals(self.phig, other.phig);
}
};
const State = struct {
iter: i32,
gbpos: []f64,
gbval: f64,
min: []f64,
max: []f64,
parameters: Parameters,
pos: [][]f64,
vel: [][]f64,
bpos: [][]f64,
bval: []f64,
n_particles: usize,
n_dims: usize,
allocator: Allocator,
fn init(allocator: Allocator, min: []const f64, max: []const f64, parameters: Parameters, n_particles: usize) !State {
const n_dims = min.len;
// Allocate memory for all arrays
const gbpos = try allocator.alloc(f64, n_dims);
const min_copy = try allocator.alloc(f64, n_dims);
const max_copy = try allocator.alloc(f64, n_dims);
const pos = try allocator.alloc([]f64, n_particles);
const vel = try allocator.alloc([]f64, n_particles);
const bpos = try allocator.alloc([]f64, n_particles);
const bval = try allocator.alloc(f64, n_particles);
// Initialize arrays
for (gbpos) |*val| val.* = math.inf(f64);
@memcpy(min_copy, min);
@memcpy(max_copy, max);
for (0..n_particles) |i| {
pos[i] = try allocator.alloc(f64, n_dims);
vel[i] = try allocator.alloc(f64, n_dims);
bpos[i] = try allocator.alloc(f64, n_dims);
@memcpy(pos[i], min);
@memset(vel[i], 0.0);
@memcpy(bpos[i], min);
bval[i] = math.inf(f64);
}
return State{
.iter = 0,
.gbpos = gbpos,
.gbval = math.inf(f64),
.min = min_copy,
.max = max_copy,
.parameters = parameters,
.pos = pos,
.vel = vel,
.bpos = bpos,
.bval = bval,
.n_particles = n_particles,
.n_dims = n_dims,
.allocator = allocator,
};
}
fn deinit(self: *State) void {
self.allocator.free(self.gbpos);
self.allocator.free(self.min);
self.allocator.free(self.max);
self.allocator.free(self.bval);
for (self.pos) |p| self.allocator.free(p);
for (self.vel) |v| self.allocator.free(v);
for (self.bpos) |bp| self.allocator.free(bp);
self.allocator.free(self.pos);
self.allocator.free(self.vel);
self.allocator.free(self.bpos);
}
fn clone(self: *const State) !State {
var new_state = try State.init(self.allocator, self.min, self.max, self.parameters, self.n_particles);
new_state.iter = self.iter;
@memcpy(new_state.gbpos, self.gbpos);
new_state.gbval = self.gbval;
@memcpy(new_state.bval, self.bval);
for (0..self.n_particles) |i| {
@memcpy(new_state.pos[i], self.pos[i]);
@memcpy(new_state.vel[i], self.vel[i]);
@memcpy(new_state.bpos[i], self.bpos[i]);
}
return new_state;
}
fn equals(self: *const State, other: *const State) bool {
return self.iter == other.iter and
vectorEqualsF64(self.gbpos, other.gbpos) and
doubleEquals(self.gbval, other.gbval) and
vectorEqualsF64(self.min, other.min) and
vectorEqualsF64(self.max, other.max) and
self.parameters.equals(other.parameters) and
vectorEquals2D(self.pos, other.pos) and
vectorEquals2D(self.vel, other.vel) and
vectorEquals2D(self.bpos, other.bpos) and
vectorEqualsF64(self.bval, other.bval) and
self.n_particles == other.n_particles and
self.n_dims == other.n_dims;
}
fn report(self: *const State, test_func: []const u8) void {
print("Test Function : {s}\n", .{test_func});
print("Iterations : {d}\n", .{self.iter});
print("Global Best Position : [" , .{});
for (self.gbpos, 0..) |val, i| {
if (i > 0) print(", ", .{});
print("{d:.5}", .{val});
}
print("]\n", .{});
print("Global Best Value : {d}\n", .{self.gbval});
}
};
fn pso(func: fn([]const f64) f64, state: *const State, rng: *Random) !State {
const p = &state.parameters;
var v = try state.allocator.alloc(f64, state.n_particles);
defer state.allocator.free(v);
var bpos = try state.allocator.alloc([]f64, state.n_particles);
defer {
for (bpos) |bp| state.allocator.free(bp);
state.allocator.free(bpos);
}
var bval = try state.allocator.alloc(f64, state.n_particles);
defer state.allocator.free(bval);
const gbpos = try state.allocator.alloc(f64, state.n_dims);
defer state.allocator.free(gbpos);
var gbval: f64 = math.inf(f64);
// Initialize arrays
for (0..state.n_particles) |i| {
bpos[i] = try state.allocator.alloc(f64, state.n_dims);
@memcpy(bpos[i], state.min);
}
@memset(bval, 0.0);
@memset(gbpos, 0.0);
// Evaluate and update best positions
for (0..state.n_particles) |j| {
v[j] = func(state.pos[j]);
if (v[j] < state.bval[j]) {
@memcpy(bpos[j], state.pos[j]);
bval[j] = v[j];
} else {
@memcpy(bpos[j], state.bpos[j]);
bval[j] = state.bval[j];
}
if (bval[j] < gbval) {
gbval = bval[j];
@memcpy(gbpos, bpos[j]);
}
}
const rg = rng.float(f64);
var pos = try state.allocator.alloc([]f64, state.n_particles);
var vel = try state.allocator.alloc([]f64, state.n_particles);
for (0..state.n_particles) |i| {
pos[i] = try state.allocator.alloc(f64, state.n_dims);
vel[i] = try state.allocator.alloc(f64, state.n_dims);
@memset(pos[i], 0.0);
@memset(vel[i], 0.0);
}
// Update particle positions and velocities
for (0..state.n_particles) |j| {
const rp = rng.float(f64);
var ok = true;
for (0..state.n_dims) |k| {
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];
}
if (!ok) {
for (0..state.n_dims) |k| {
pos[j][k] = state.min[k] + (state.max[k] - state.min[k]) * rng.float(f64);
}
}
}
// Create new state
var new_state = State{
.iter = state.iter + 1,
.gbpos = try state.allocator.alloc(f64, state.n_dims),
.gbval = gbval,
.min = try state.allocator.alloc(f64, state.n_dims),
.max = try state.allocator.alloc(f64, state.n_dims),
.parameters = state.parameters,
.pos = pos,
.vel = vel,
.bpos = try state.allocator.alloc([]f64, state.n_particles),
.bval = try state.allocator.alloc(f64, state.n_particles),
.n_particles = state.n_particles,
.n_dims = state.n_dims,
.allocator = state.allocator,
};
@memcpy(new_state.gbpos, gbpos);
@memcpy(new_state.min, state.min);
@memcpy(new_state.max, state.max);
@memcpy(new_state.bval, bval);
for (0..state.n_particles) |i| {
new_state.bpos[i] = try state.allocator.alloc(f64, state.n_dims);
@memcpy(new_state.bpos[i], bpos[i]);
}
return new_state;
}
fn iterate(func: fn([]const f64) f64, n: ?i32, state: *State, rng: *Random) !State {
var result = try state.clone();
if (n) |iterations| {
for (0..@intCast(iterations)) |_| {
var old_result = result;
result = try pso(func, &result, rng);
old_result.deinit();
}
} else {
// Iterate until convergence
while (true) {
var old = try result.clone();
const new_result = try pso(func, &result, rng);
result.deinit();
result = new_result;
if (result.equals(&old)) {
old.deinit();
break;
}
old.deinit();
}
}
return result;
}
fn mccormick(x: []const f64) f64 {
const a = x[0];
const b = x[1];
return math.sin(a + b) + math.pow(f64, a - b, 2.0) + 1.0 + 2.5 * b - 1.5 * a;
}
fn michalewicz(x: []const f64) f64 {
const m = 10;
const d = x.len;
var sum: f64 = 0.0;
for (1..d) |i| {
const j = x[i - 1];
const k = math.sin(@as(f64, @floatFromInt(i)) * j * j / math.pi);
sum += math.sin(j) * math.pow(f64, k, 2.0 * @as(f64, @floatFromInt(m)));
}
return -sum;
}
pub fn main() !void {
var gpa = std.heap.GeneralPurposeAllocator(.{}){};
defer _ = gpa.deinit();
const allocator = gpa.allocator();
var prng = std.Random.DefaultPrng.init(@as(u64, @intCast(std.time.milliTimestamp())));
var rng = prng.random();
// Test McCormick function
var state = try State.init(
allocator,
&[_]f64{ -1.5, -3.0 },
&[_]f64{ 4.0, 4.0 },
Parameters{
.omega = 0.0,
.phip = 0.6,
.phig = 0.3,
},
100,
);
var result1 = try iterate(mccormick, 40, &state, &rng);
state.deinit();
result1.report("McCormick");
print("f(-0.54719, -1.54719) : {d}\n", .{mccormick(&[_]f64{ -0.54719, -1.54719 })});
print("\n", .{});
result1.deinit();
// Test Michalewicz function
var state2 = try State.init(
allocator,
&[_]f64{ 0.0, 0.0 },
&[_]f64{ math.pi, math.pi },
Parameters{
.omega = 0.3,
.phip = 0.3,
.phig = 0.3,
},
1000,
);
var result2 = try iterate(michalewicz, 30, &state2, &rng);
state2.deinit();
result2.report("Michalewicz (2D)");
print("f(2.20, 1.57) : {d}\n", .{michalewicz(&[_]f64{ 2.2, 1.57 })});
result2.deinit();
}
- Output:
Test Function : McCormick Iterations : 40 Global Best Position : [-0.54770, -1.54759] Global Best Value : -1.9132225947683208 f(-0.54719, -1.54719) : -1.913222954882274 Test Function : Michalewicz (2D) Iterations : 30 Global Best Position : [2.20291, 2.48195] Global Best Value : -0.8013034100853031 f(2.20, 1.57) : -0.8011663878202854