Particle swarm optimization

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.

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 2D Rosenbrock function, and possibly other standard or well-known optimization test cases.

References:

• [Particle Swarm Optimization[1]]
• [Rosenbrock function[2]]
• [Test functions for optimization[3]]

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
  )

  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│GlobalBestPosition│GlobalBestValue│ ├─────────┼──────────────────┼───────────────┤ │40 │_0.547804 _1.54794│_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