Particle swarm optimization: Difference between revisions
(→{{header|J}}: tacit version of McCormick Function) |
(→{{header|J}}: Use local assignment within verb) |
||
Line 15: | Line 15: | ||
</p> |
</p> |
||
=={{header|J}}== |
=={{header|J}}== |
||
<lang J> |
<lang J>pso_init =: 3 : 0 |
||
pso_init =: 3 : 0 |
|||
'Min Max parameters nParticles' =. y |
'Min Max parameters nParticles' =. y |
||
smoutput 4 2 $ 'Min';Min;'Max';Max;'omega, phip, phig';parameters;'nParticles';nParticles |
smoutput 4 2 $ 'Min';Min;'Max';Max;'omega, phip, phig';parameters;'nParticles';nParticles |
||
Line 36: | Line 35: | ||
NB. evaluate |
NB. evaluate |
||
val = |
val =. pso_function"1 pos |
||
NB. update |
NB. update |
||
better = |
better =. val < bval |
||
bpos = |
bpos =. (better * pos) + ((1-better) * bpos0) |
||
bval = |
bval =. pso_function"1 bpos |
||
gbval = |
gbval =. <./,bval |
||
gmask = |
gmask =. gbval = bval |
||
gindex = |
gindex =. +/(gmask*(i.#gmask)) |
||
gbpos = |
gbpos =. gindex { bpos |
||
NB. migrate |
NB. migrate |
||
omega = |
omega =. 0{parameters |
||
phip = |
phip =. 1{parameters |
||
phig = |
phig =. 2{parameters |
||
rp = |
rp =. 1e6%~?(#pos)$1e6 |
||
rg = |
rg =. 1e6%~?1e6 |
||
vel = |
vel =. (omega*vel) + (phip * rp * (bpos-pos)) + (phig * rg * (gbpos -"1 1 pos)) |
||
pos = |
pos =. pos + vel |
||
NB. reset out-of-bounds particles |
NB. reset out-of-bounds particles |
||
pmask = |
pmask =. Min <"1 pos +. pos <"1 Max |
||
rnd = |
rnd =. (1e6%~ ? ($pos) $ 1e6) |
||
Range = |
Range =. Max - Min |
||
newpos = |
newpos =. |: Min + Range * |: rnd |
||
pos = |
pos =. (pmask * pos) + ((1-pmask) * newpos) |
||
iter = |
iter =. >: iter |
||
NB. new state |
NB. new state |
||
iter;gbpos;gbval;Min;Max;parameters;pos;vel;bpos;bval |
iter;gbpos;gbval;Min;Max;parameters;pos;vel;bpos;bval |
||
⚫ | |||
) |
|||
⚫ | |||
Apply to McCormick Function: |
Apply to McCormick Function: |
||
<lang J> |
<lang J> load'trig' |
||
load'trig' |
|||
pso_function =: 3 : 0 |
pso_function =: 3 : 0 |
||
(sin (0{y)+(1{y)) + (((0{y) - (1{y))^2) + (_1.5 * (0{y)) + (2.5 * (1{y)) + 1 |
(sin (0{y)+(1{y)) + (((0{y) - (1{y))^2) + (_1.5 * (0{y)) + (2.5 * (1{y)) + 1 |
||
Line 99: | Line 95: | ||
├──────────────────┼──────────────────┤ |
├──────────────────┼──────────────────┤ |
||
│GlobalBestValue │_1.91322 │ |
│GlobalBestValue │_1.91322 │ |
||
└──────────────────┴──────────────────┘ |
└──────────────────┴──────────────────┘</lang> |
||
</lang> |
|||
=={{header|ooRexx}}== |
=={{header|ooRexx}}== |
Revision as of 05:24, 3 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 McCormick function, and possibly other standard or well-known optimization test cases.
References:
J
<lang J>pso_init =: 3 : 0
'Min Max parameters nParticles' =. y smoutput 4 2 $ 'Min';Min;'Max';Max;'omega, phip, phig';parameters;'nParticles';nParticles Range =. Max - Min nDims =. #Min searchSpaceBounds =. |: (2,nDims) $ Min, Max rnd =. (1e6%~ ? (nParticles,nDims) $ 1e6) pos =. |: Min + Range * |: rnd bpos =. pos bval =. (#pos) $ _ vel =. ($pos) $ 0 0;_;_;Min;Max;parameters;pos;vel;bpos;bval NB. initial state
)
pso =: 3 : 0
NB. previous state 'iter gbpos gbval Min Max parameters pos vel bpos0 bval' =: y
NB. evaluate val =. pso_function"1 pos
NB. update better =. val < bval bpos =. (better * pos) + ((1-better) * bpos0) bval =. pso_function"1 bpos gbval =. <./,bval gmask =. gbval = bval gindex =. +/(gmask*(i.#gmask)) gbpos =. gindex { bpos
NB. migrate omega =. 0{parameters phip =. 1{parameters phig =. 2{parameters rp =. 1e6%~?(#pos)$1e6 rg =. 1e6%~?1e6 vel =. (omega*vel) + (phip * rp * (bpos-pos)) + (phig * rg * (gbpos -"1 1 pos)) pos =. pos + vel
NB. reset out-of-bounds particles pmask =. Min <"1 pos +. pos <"1 Max rnd =. (1e6%~ ? ($pos) $ 1e6) Range =. Max - Min newpos =. |: Min + Range * |: rnd pos =. (pmask * pos) + ((1-pmask) * newpos) iter =. >: iter
NB. new state iter;gbpos;gbval;Min;Max;parameters;pos;vel;bpos;bval
)</lang> Apply to McCormick Function: <lang J> load'trig'
pso_function =: 3 : 0 (sin (0{y)+(1{y)) + (((0{y) - (1{y))^2) + (_1.5 * (0{y)) + (2.5 * (1{y)) + 1 )
pso_function =: sin@(+/) + *:@(-/) + 1 + _1.5 2.5 +/@:* ] NB. tacit version
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 =: pso^:40 state smoutput |: 2 3 $ 'Iteration';'GlobalBestPosition';'GlobalBestValue';iter;gbpos;gbval
┌──────────────────┬──────────────────┐ │Iteration │40 │ ├──────────────────┼──────────────────┤ │GlobalBestPosition│_0.547599 _1.54788│ ├──────────────────┼──────────────────┤ │GlobalBestValue │_1.91322 │ └──────────────────┴──────────────────┘</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 16 digits).
The numeric precision is only limited to the number of decimal digits in the pi variable (in this case, 77). <lang rexx>/*REXX pgm calc. Particle Swarm Optimization as it migrates through a solution*/ numeric digits length(pi()); sDig=16 /*SDIG: the number of displayed digits.*/ parse arg x y d p . /*obtain optional arguments from the CL*/ if x== | x==',' then x= -0.5 /*X not defined? Then use the default.*/ if y== | y==',' then y= -1.5 /*Y " " " " " " */ if d== | d==',' then d= 1 /*D " " " " " " */ if p== | p==',' then p= 1e12 /*P " " " " " " */ minF=p /*P the number of particles: 1 billion*/ say center('X', sDig+3, '═') center('Y', sDig+3, '═') center('D', sDig+3, '═') call refine x,y
do r=1 for 10; d=d*.5 call refine minX, minY end /*r*/
say say 'Which is better (less) than the global minimum at:' say ' f(-.54719, -1.54719) ───► ' fmt(f(-.54719, -1.54719)) say 'The published global minimum is: -1.9133' exit /*────────────────────────────────────────────────────────────────────────────*/ 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 /*────────────────────────────────────────────────────────────────────────────*/ fmt: parse arg ?; ?=format(?,,sDig) /*format number with Sdig decimal digs.*/ L=length(?); if pos(.,?)\==0 then ?=strip(strip(?,'T',0),'T',.);return left(?,L) /*────────────────────────────────────────────────────────────────────────────*/ 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
/*──────────────────────────────────one─liner subroutines──────────────────────────────────────*/ f: procedure: parse arg a,b; return sin(a+b) + (a-b)**2 - 1.5*a + 2.5*b + 1 pi: pi=3.1415926535897932384626433832795028841971693993751058209749445923078164062862; return pi r2r: return arg(1) // (pi()*2) /*normalize radians ───► a unit circle.*/</lang> output when using the default inputs:
═════════X═════════ ═════════Y═════════ ═════════D═════════ -1.5 -2.5 -1.2431975046920717 -1 -2 -1.6411200080598672 -0.5 -1.5 -1.9092974268256817 -0.5625 -1.5625 -1.912819789818452 -0.5625 -1.546875 -1.9128819293954732 -0.546875 -1.546875 -1.9132227747573614 -0.54736328125 -1.54736328125 -1.9132229074107836 Which is better (less) than the global minimum at: f(-.54719, -1.54719) ───► -1.9132229548822736 The published global minimum is: -1.9133