Particle swarm optimization: Difference between revisions
(→{{header|REXX}}: added/changed whitespace and comments, add more variables that can be specified, increased the number of iterations and shown digits, increased the precision, and other refinements.) |
(Javascript draft: TODO - implement and run the Michalewicz function) |
||
| Line 115: | Line 115: | ||
GlobalBestPosition: 2.20296 1.57083 |
GlobalBestPosition: 2.20296 1.57083 |
||
GlobalBestValue: _1.8013</lang> |
GlobalBestValue: _1.8013</lang> |
||
=={{header|Javascript}}== |
|||
Translation of [[Particle_Swarm_Optimization#J|J]]. |
|||
<lang Javascript>function pso_init(y) { |
|||
var nDims= y.min.length; |
|||
var pos=[], vel=[], bpos=[], bval=[]; |
|||
for (var j= 0; j<y.nParticles; j++) { |
|||
pos[j]= bpos[j]= y.min; |
|||
var v= []; for (var k= 0; k<nDims; k++) v[k]= 0; |
|||
vel[j]= v; |
|||
bval[j]= Infinity} |
|||
return { |
|||
iter: 0, |
|||
gbpos: Infinity, |
|||
gbval: Infinity, |
|||
min: y.min, |
|||
max: y.max, |
|||
parameters: y.parameters, |
|||
pos: pos, |
|||
vel: vel, |
|||
bpos: bpos, |
|||
bval: bval, |
|||
nParticles: y.nParticles, |
|||
nDims: nDims} |
|||
} |
|||
function pso(fn, state) { |
|||
var y= state; |
|||
var p= y.parameters; |
|||
var val=[], bpos=[], bval=[], gbval= Infinity, gbpos=[] |
|||
for (var j= 0; j<y.nParticles; j++) { |
|||
// evaluate |
|||
val[j]= fn.apply(null, y.pos[j]); |
|||
// update |
|||
if (val[j] < y.bval[j]) { |
|||
bpos[j]= y.pos[j]; |
|||
bval[j]= val[j]; |
|||
} else { |
|||
bpos[j]= y.bpos[j]; |
|||
bval[j]= y.bval[j]} |
|||
if (bval[j] < gbval) { |
|||
gbval= bval[j]; |
|||
gbpos= bpos[j]}} |
|||
var rg= Math.random(), vel=[], pos=[]; |
|||
for (var j= 0; j<y.nParticles; j++) { |
|||
// migrate |
|||
var rp= Math.random(), ok= true; |
|||
vel[j]= []; |
|||
pos[j]= []; |
|||
for (var k= 0; k < y.nDims; k++) { |
|||
vel[j][k]= p.omega*y.vel[j][k] + p.phip*rp*(bpos[j]-y.pos[j]) + p.phig*rg*(gbpos-y.pos[j]); |
|||
pos[j][k]= y.pos[j]+vel[j][k]; |
|||
ok= ok && y.min>pos[j][k] || y.max<pos[j][k];} |
|||
if (!ok) |
|||
for (var k= 0; k < y.nDims; k++) |
|||
pos[j][k]= y.min + (y.max-y.min)*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,2], max:[4,4], |
|||
parameters: {omega: 0, phip: 0.6, phig: 0.3}, |
|||
nParticles: 100}); |
|||
reportState(state); |
|||
state= repeat(function(y){return pso(mccormick,y)}, 40, state); |
|||
reportState(state);</lang> |
|||
Example displayed result (random numbers are involved so there will be a bit of variance between repeated runs: |
|||
<lang Javascript> |
|||
Iteration: 0 |
|||
GlobalBestPosition: Infinity |
|||
GlobalBestValue: Infinity |
|||
Iteration: 40 |
|||
GlobalBestPosition: -1.5,2 |
|||
GlobalBestValue: 20.979425538604204</lang> |
|||
=={{header|ooRexx}}== |
=={{header|ooRexx}}== |
||
Revision as of 19:53, 15 August 2015
Particle Swarm Optimization (PSO) is an optimization method in which multiple candidate solutions ('particles') migrate through the solution space under the influence of local and global best known positions. PSO does not require that the objective function be differentiable and can optimize over very large problem spaces, but is not guaranteed to converge. The method should be demonstrated by application of the functions recommended below, and possibly other standard or well-known optimization test cases.
The goal of parameter selection is to ensure that the global minimum is discriminated from any local minima, and that the minimum is accurately determined, and that convergence is achieved with acceptible resource usage. To provide a common basis for comparing implementations, the following test cases are recommended:
- McCormick function - bowl-shaped, with a single minimum
- function parameters and bounds (recommended):
- -1.5 < x1 < 4
- -3 < x2 < 4
- search parameters (suggested):
- omega = 0
- phi p = 0.6
- phi g = 0.3
- number of particles = 100
- number of iterations = 40
- Michalewicz function - steep ridges and valleys, with multiple minima
- function parameters and bounds (recommended):
- 0 < x1 < pi
- 0 < x2 < pi
- search parameters (suggested):
- omega = 0.3
- phi p = 0.3
- phi g = 0.3
- number of particles = 1000
- number of iterations = 30
References:
J
<lang J>load 'format/printf'
pso_init =: verb define
'Min Max parameters nParticles' =. y 'Min: %j\nMax: %j\nomega, phip, phig: %j\nnParticles: %j\n' printf Min;Max;parameters;nParticles nDims =. #Min pos =. Min +"1 (Max - Min) *"1 (nParticles,nDims) ?@$ 0 bpos =. pos bval =. (#pos) $ _ vel =. ($pos) $ 0 0;_;_;Min;Max;parameters;pos;vel;bpos;bval NB. initial state
)
pso =: adverb define
NB. previous state 'iter gbpos gbval Min Max parameters pos vel bpos0 bval' =. y
NB. evaluate val =. u"1 pos
NB. update
better =. val < bval
bpos =. (better # pos) (I. better)} bpos0
bval =. u"1 bpos
gbval =. <./ bval
gbpos =. bpos {~ (i. <./) bval
NB. migrate 'omega phip phig' =. parameters rp =. (#pos) ?@$ 0 rg =. ? 0 vel =. (omega*vel) + (phip * rp * bpos - pos) + (phig * rg * gbpos -"1 pos) pos =. pos + vel
NB. reset out-of-bounds particles bad =. +./"1 (Min >"1 pos) ,. (pos >"1 Max) newpos =. Min +"1 (Max-Min) *"1 ((+/bad),#Min) ?@$ 0 pos =. newpos (I. bad)} pos iter =. >: iter
NB. new state iter;gbpos;gbval;Min;Max;parameters;pos;vel;bpos;bval
)
reportState=: 'Iteration: %j\nGlobalBestPosition: %j\nGlobalBestValue: %j\n' printf 3&{.</lang> Apply to McCormick Function:<lang J> require 'trig'
mccormick =: sin@(+/) + *:@(-/) + 1 + _1.5 2.5 +/@:* ]
state =: pso_init _1.5 _3 ; 4 4 ; 0 0.6 0.3; 100
Min: _1.5 _3 Max: 4 4 omega, phip, phig: 0 0.6 0.3 nParticles: 100
state =: (mccormick pso)^:40 state reportState state
Iteration: 40 GlobalBestPosition: _0.547399 _1.54698 GlobalBestValue: _1.91322</lang> Apply to Michalewicz Function: <lang J> michalewicz =: 3 : '- +/ (sin y) * 20 ^~ sin (>: i. #y) * (*:y) % pi'
michalewicz =: [: -@(+/) sin * 20 ^~ sin@(pi %~ >:@i.@# * *:) NB. tacit equivalent state =: pso_init 0 0 ; (pi,pi) ; 0.3 0.3 0.3; 1000
Min: 0 0 Max: 3.14159 3.14159 omega, phip, phig: 0.3 0.3 0.3 nParticles: 1000
state =: (michalewicz pso)^:30 state reportState state
Iteration: 30 GlobalBestPosition: 2.20296 1.57083 GlobalBestValue: _1.8013</lang>
JavaScript
Translation of J.
<lang Javascript>function pso_init(y) {
var nDims= y.min.length;
var pos=[], vel=[], bpos=[], bval=[];
for (var j= 0; j<y.nParticles; j++) {
pos[j]= bpos[j]= y.min;
var v= []; for (var k= 0; k<nDims; k++) v[k]= 0;
vel[j]= v;
bval[j]= Infinity}
return {
iter: 0, gbpos: Infinity, gbval: Infinity, min: y.min, max: y.max, parameters: y.parameters, pos: pos, vel: vel, bpos: bpos, bval: bval,
nParticles: y.nParticles,
nDims: nDims}
}
function pso(fn, state) {
var y= state;
var p= y.parameters;
var val=[], bpos=[], bval=[], gbval= Infinity, gbpos=[]
for (var j= 0; j<y.nParticles; j++) {
// evaluate
val[j]= fn.apply(null, y.pos[j]);
// update
if (val[j] < y.bval[j]) {
bpos[j]= y.pos[j];
bval[j]= val[j];
} else {
bpos[j]= y.bpos[j];
bval[j]= y.bval[j]}
if (bval[j] < gbval) {
gbval= bval[j];
gbpos= bpos[j]}}
var rg= Math.random(), vel=[], pos=[];
for (var j= 0; j<y.nParticles; j++) {
// migrate
var rp= Math.random(), ok= true;
vel[j]= [];
pos[j]= [];
for (var k= 0; k < y.nDims; k++) {
vel[j][k]= p.omega*y.vel[j][k] + p.phip*rp*(bpos[j]-y.pos[j]) + p.phig*rg*(gbpos-y.pos[j]);
pos[j][k]= y.pos[j]+vel[j][k];
ok= ok && y.min>pos[j][k] || y.max<pos[j][k];}
if (!ok)
for (var k= 0; k < y.nDims; k++)
pos[j][k]= y.min + (y.max-y.min)*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,2], max:[4,4],
parameters: {omega: 0, phip: 0.6, phig: 0.3},
nParticles: 100});
reportState(state);
state= repeat(function(y){return pso(mccormick,y)}, 40, state);
reportState(state);</lang>
Example displayed result (random numbers are involved so there will be a bit of variance between repeated runs:
<lang Javascript> Iteration: 0 GlobalBestPosition: Infinity GlobalBestValue: Infinity
Iteration: 40 GlobalBestPosition: -1.5,2 GlobalBestValue: 20.979425538604204</lang>
ooRexx
<lang oorexx>/* REXX ---------------------------------------------------------------
- Test for McCormick function
- --------------------------------------------------------------------*/
Numeric Digits 16 Parse Value '-.5 -1.5 1' With x y d fmin=1e9 Call refine x,y Do r=1 To 10
d=d/5 Call refine xmin,ymin End
Say 'which is better (less) than' Say ' f(-.54719,-1.54719)='f(-.54719,-1.54719) Say 'and differs from published -1.9133' Exit
refine: Parse Arg xx,yy Do x=xx-d To xx+d By d/2
Do y=yy-d To yy+d By d/2
f=f(x,y)
If f<fmin Then Do
Say x y f
fmin=f
xmin=x
ymin=y
End
End
End
Return
f: Parse Arg x,y res=rxcalcsin(x+y,16,'R')+(x-y)**2-1.5*x+2.5*y+1 Return res
- requires rxmath library</lang>
- Output:
-1.5 -2.5 -1.243197504692072
-1.0 -2.0 -1.641120008059867
-0.5 -1.5 -1.909297426825682
-0.54 -1.54 -1.913132979507516
-0.548 -1.548 -1.913221840016527
-0.5480 -1.5472 -1.913222034492829
-0.5472 -1.5472 -1.913222954970650
-0.54720000 -1.54719872 -1.913222954973731
-0.54719872 -1.54719872 -1.913222954978670
-0.54719872 -1.54719744 -1.913222954978914
-0.54719744 -1.54719744 -1.913222954981015
-0.5471975424 -1.5471975424 -1.913222954981036
which is better (less) than
f(-.54719,-1.54719)=-1.913222954882273
and differs from published -1.9133
REXX
This REXX version uses a large numeric digits (but only displays 25 digits).
Classic REXX doesn't have a sine function, so a RYO version is included here.
The numeric precision is only limited to the number of decimal digits defined in the pi variable (in this case, 100).
This REXX version supports the specifying of X, Y, and D, as well as the number of particles, the number of times the
computation loop is performed, and the number of decimal digits to be displayed.
The refinement loop is stopped when the function value stabilizes, or the limit of iterations is reached. <lang rexx>/*REXX pgm calc. Particle Swarm Optimization as it migrates through a solution*/ numeric digits length(pi()) /*sDigs: is the # of displayed digits.*/ parse arg x y d #part times sDigs . /*obtain optional arguments from the CL*/ if x== | x==',' then x= -0.5 /*is X not defined?*/ if y== | y==',' then y= -1.5 /* " Y " " */ if d== | d==',' then d= 1 /* " D " " */ if #part== | #part==',' then #part=1e12 /* " the # particles " " */ if times== | times==',' then times= 40 /* " the # of times " " */ if sDigs== | sDigs==',' then sDigs= 25 /* " the # of digits " " */ minF=#part /*number of particles is one billion. */ say center('X',sDigs+3,'═') center('Y',sDigs+3,'═') center('D',sDigs+3,'═') call refine x,y
do stuff=1 for times until old=f; d=d*.2; old=f
call refine minX, minY
end /*stuff*/ /* [↑] stop refining after TIMES, or */
say /* when the value of F stabilizes.*/ indent=1 + 2*(sDigs+3) /*compute the indentation for alignment*/ say right('The global minimum for f(-.54719, -1.54719) ───► ', indent) fmt(f(-.54719, -1.54719)) say right('The published global minimum is:' , indent) fmt( -1.9133 ) exit /*stick a fork in it, we're all done. */ /*────────────────────────────────────────────────────────────────────────────*/ refine: parse arg xx,yy; dh=d * 0.5
do x=xx-d to xx+d by dh
do y=yy-d to yy+d by dh; f=f(x,y); if f>=minF then iterate
say fmt(x) fmt(y) fmt(f); minF=f; minX=x; minY=y
end /*y*/
end /*x*/
return /*──────────────────────────────────────────────────────────────────────────────────one─liner subroutines───────────────────────────────*/ f: procedure: parse arg a,b; return sin(a+b) + (a-b)**2 - 1.5*a + 2.5*b + 1 fmt: parse arg ?; ?=format(?,,sDigs); L=length(?); if pos(.,?)\==0 then ?=strip(strip(?,'T',0),'T',.); return left(?,L) pi: pi=3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068; return pi r2r: return arg(1) // (pi()*2) /*normalize radians ───► a unit circle.*/ sin: procedure; parse arg x; x=r2r(x); numeric fuzz 5; z=x; _=x; q=x*x; do k=2 by 2 until p=z; p=z; _=-_*q/(k*(k+1)); z=z+_; end; return z</lang> output when using the default inputs:
═════════════X══════════════ ═════════════Y══════════════ ═════════════D══════════════
-1.5 -2.5 -1.2431975046920717486273609
-1 -2 -1.6411200080598672221007448
-0.5 -1.5 -1.9092974268256816953960199
-0.54 -1.54 -1.9131329795075164948766768
-0.548 -1.548 -1.9132218400165267634506035
-0.548 -1.5472 -1.9132220344928294065568196
-0.5472 -1.5472 -1.9132229549706499208388746
-0.5472 -1.54719872 -1.9132229549737311254290577
-0.54719872 -1.54719872 -1.9132229549786702369612333
-0.54719872 -1.54719744 -1.91322295497891365438682
-0.54719744 -1.54719744 -1.9132229549810149766572388
-0.5471975424 -1.5471975424 -1.9132229549810362588916172
-0.54719755264 -1.54719755264 -1.9132229549810363893093655
-0.547197550592 -1.547197550592 -1.9132229549810363922848065
-0.5471975514112 -1.5471975514112 -1.9132229549810363928381695
-0.5471975510016 -1.5471975510016 -1.9132229549810363928520779
-0.54719755116544 -1.54719755116544 -1.9132229549810363929162561
-0.547197551198208 -1.547197551198208 -1.9132229549810363929179331
-0.547197551198208 -1.54719755119755264 -1.9132229549810363929179344
-0.54719755119755264 -1.54719755119755264 -1.9132229549810363929179361
-0.54719755119755264 -1.54719755119689728 -1.9132229549810363929179365
-0.54719755119689728 -1.54719755119689728 -1.9132229549810363929179375
-0.54719755119689728 -1.547197551196766208 -1.9132229549810363929179375
-0.547197551196766208 -1.547197551196766208 -1.9132229549810363929179376
-0.547197551196766208 -1.547197551196635136 -1.9132229549810363929179376
-0.547197551196635136 -1.547197551196635136 -1.9132229549810363929179376
-0.547197551196635136 -1.5471975511966089216 -1.9132229549810363929179376
-0.5471975511966089216 -1.5471975511966089216 -1.9132229549810363929179376
-0.5471975511966089216 -1.54719755119660367872 -1.9132229549810363929179376
-0.54719755119660367872 -1.54719755119660367872 -1.9132229549810363929179376
-0.54719755119660367872 -1.54719755119659843584 -1.9132229549810363929179376
-0.54719755119659843584 -1.54719755119659843584 -1.9132229549810363929179376
-0.547197551196597387264 -1.547197551196597387264 -1.9132229549810363929179376
-0.5471975511965978066944 -1.5471975511965978066944 -1.9132229549810363929179376
-0.5471975511965978066944 -1.54719755119659776475136 -1.9132229549810363929179376
-0.54719755119659776475136 -1.54719755119659776475136 -1.9132229549810363929179376
-0.54719755119659776475136 -1.547197551196597756362752 -1.9132229549810363929179376
-0.547197551196597756362752 -1.547197551196597756362752 -1.9132229549810363929179376
-0.547197551196597756362752 -1.547197551196597747974144 -1.9132229549810363929179376
-0.547197551196597747974144 -1.547197551196597747974144 -1.9132229549810363929179376
-0.547197551196597747974144 -1.5471975511965977462964224 -1.9132229549810363929179376
-0.5471975511965977462964224 -1.5471975511965977462964224 -1.9132229549810363929179376
-0.5471975511965977462964224 -1.5471975511965977462293135 -1.9132229549810363929179376
-0.5471975511965977462293135 -1.5471975511965977462293135 -1.9132229549810363929179376
-0.5471975511965977462293135 -1.5471975511965977461622047 -1.9132229549810363929179376
-0.5471975511965977461622047 -1.5471975511965977461622047 -1.9132229549810363929179376
-0.5471975511965977461487829 -1.5471975511965977461487829 -1.9132229549810363929179376
-0.5471975511965977461541516 -1.5471975511965977461541516 -1.9132229549810363929179376
-0.547197551196597746154259 -1.547197551196597746154259 -1.9132229549810363929179376
-0.547197551196597746154259 -1.5471975511965977461542375 -1.9132229549810363929179376
-0.5471975511965977461542375 -1.5471975511965977461542375 -1.9132229549810363929179376
-0.5471975511965977461542375 -1.547197551196597746154216 -1.9132229549810363929179376
-0.547197551196597746154216 -1.547197551196597746154216 -1.9132229549810363929179376
-0.547197551196597746154216 -1.5471975511965977461542152 -1.9132229549810363929179376
-0.5471975511965977461542152 -1.5471975511965977461542152 -1.9132229549810363929179376
-0.5471975511965977461542152 -1.5471975511965977461542143 -1.9132229549810363929179376
-0.5471975511965977461542143 -1.5471975511965977461542143 -1.9132229549810363929179376
-0.5471975511965977461542145 -1.5471975511965977461542145 -1.9132229549810363929179376
-0.5471975511965977461542145 -1.5471975511965977461542145 -1.9132229549810363929179376
-0.5471975511965977461542145 -1.5471975511965977461542145 -1.9132229549810363929179376
-0.5471975511965977461542145 -1.5471975511965977461542145 -1.9132229549810363929179376
-0.5471975511965977461542145 -1.5471975511965977461542145 -1.9132229549810363929179376
-0.5471975511965977461542145 -1.5471975511965977461542145 -1.9132229549810363929179376
-0.5471975511965977461542145 -1.5471975511965977461542145 -1.9132229549810363929179376
The global minimum for f(-.54719, -1.54719) ───► -1.9132229548822735814541188
The published global minimum is: -1.9133
Output note: the published global minimum (referenced above, as well as the function's arguments) can be found at: