Particle Swarm Optimization

initialPSOPopulation[individuals_Integer,locationRanges_List,velocityRanges_List,fitFct_]:=
Module[{locations,velocities,fitnesses,locationsBestPrivate,locationsBestNeighbourhood, f, neighbour
},
(* calculate the fitness *)

f[x__]:=Max[Re[Map[fitFct,{x}]]];


(* get neigbours *)

neighbour[func__,lst_,r_Integer:1]:=MapThread[func,Map[RotateRight[lst,#]&,Range[-r,r]]];

(*initialize the locations *)
locations=Table[Table[Random[Real,locationRanges[[d]]],{d,1,Length[locationRanges]}],{individuals}];

velocities=Table[Table[Random[Real,velocityRanges[[d]]],{d,1,Length[velocityRanges]}],{individuals}];

fitnesses=Map[fitFct,locations];

locationsBestPrivate=locations;

locationsBestNeighbourhood=neighbour[f,locationsBestPrivate, 1];

{locations,
velocities,
fitnesses,
locationsBestPrivate,
locationsBestNeighbourhood}
]

psoStep[psoPopulation_List,locationRanges_List,velocityRanges_List,fitFct_,φ1_Real:1.0,φ2_Real:1.0]:=

Module[{locations,velocities,fitnesses,locationsBestPrivate,locationsBestNeighbourhood,f, neighbour,n},

f[x__]:=Max[Re[Map[fitFct,{x}]]];


neighbour[func__,lst_,r_Integer:1]:=MapThread[func,Map[RotateRight[lst,#]&,Range[-r,r]]];

locationsBestPrivate=MapThread[If[fitFct[#1]>fitFct[#2],#1,#2]&,{psoPopulation[[4]],psoPopulation[[1]]}];
(* Print["locationsBestPrivate: ",locationsBestPrivate]; *)

locationsBestNeighbourhood=bestNeighbourPerformer[locationsBestPrivate,1,fitFct];
(* Print["locationsBestNeighbourhood: ",locationsBestNeighbourhood]; *)

velocities=MapThread[
#2+φ1(#3-#1)+φ2(#4-#1)&,
{psoPopulation[[1]],psoPopulation[[2]],locationsBestPrivate,locationsBestNeighbourhood}];
(*Print["velocities: ",velocities]; *)

(* Enforce constraints *)
velocities=Map[MapThread[Min[Max[#1,First[#2]],Last[#2]]&,{#,velocityRanges}]&,velocities];

(*Print["Constrained velocities: ",velocities];*)

locations=MapThread[#1+#2&,{psoPopulation[[1]],velocities}];
(*Print["locations: ",locations];*)

(* Enforce constraints *)
locations=Map[MapThread[Min[Max[#1,First[#2]],Last[#2]]&,{#,locationRanges}]&,locations];

(*Print["Constrained locations: ",locations];*)

fitnesses=Map[fitFct,locations];

{locations,
velocities,
fitnesses,
locationsBestPrivate,
locationsBestNeighbourhood}
]

bestNeighbourPerformer[lst_,r_Integer:1,fitFct_]:=
Module[{neighbour,bestFits,bestFitsIndices,neighboursList,neighbourFits,bestNeighbourLocations},
neighbour[func_,list_,rr_Integer:1]:=MapThread[func,Map[RotateRight[list,#]&,Range[-rr,rr]]];

neighboursList=neighbour[List,lst,r];
neighbourFits=Map[Map[fitFct,#]&,neighboursList];

bestFits=Map[Max,neighbourFits];
bestFitsIndices=MapThread[First[Position[#1,#2,1]]&,{neighbourFits,bestFits}];

bestNeighbourLocations=MapThread[Part[#1,First[#2]]&,{neighboursList,bestFitsIndices}];

(* {neighboursList,neighbourFits,bestFits,bestFitsIndices,bestNeighbourLocations} *)

bestNeighbourLocations
]

plotPSO[psoPopulation_List,fitFct_,dimension_Integer]:=
Module[{secondPlot},
If[dimension ≥2 ,

If[dimension ==2,
secondPlot=fitPlot3D;
MapIndexed[Show[{secondPlot,Graphics3D[{RGBColor[0.8,0,0],PointSize[0.02],Map[Point[Flatten[{#1[[1]],#1[[2]],fitFct[#1]+0.5}]]&,#]}]},AspectRatio->Automatic,
PlotLabel->("PSO"),DisplayFunction->$DisplayFunction]&,Map[First[#]&,psoPopulation]],

Print ["Multi-dimensional is available for Contour and Density Plots"];


],
secondPlot=fitPlot2D;
MapIndexed[Show[{secondPlot,Graphics[Map[{RGBColor[0.8,0,0],Disk[{First[#1],fitFct[#1]},.15]}&,#1]]},AspectRatio->Automatic,
Frame->True,
FrameLabel->{None,None,("PSO"),None},DisplayFunction->$DisplayFunction]&,Map[First[#]&,psoPopulation]]

]
]

fitness2D[{x__}]:=-(((1.2+((Cos[ (((Apply[Plus,{x}]))^2)])/(1+0.01 (((Apply[Plus,{x}]))^2)))))^(-1))

g[v_List]:=(50 Sin[√Plus@@((((#1 ))^2)&)/@v])/(√Plus@@((((#1 ))^2)&)/@v)-√Plus@@((((#1 ))^2)&)/@v

fitness3D[v_List]:=g[v+{11,9}]+g[v-{11,3}]+g[v+{6,-9}]

Manipulate[
Switch[funSelection, 
func1,
  𝕀=20;
       λ={{-8,8}};
       μ={{-mu,mu}};
  ℱ = fitness2D;
      fitPlot2D=Plot[ℱ[{x}],{x,-10,10},
PlotRange->{Automatic,{-5,0}},DisplayFunction->Identity];
 pop=initialPSOPopulation[𝕀,λ,μ,ℱ],

func2,
  𝕀=20;
      λ={{-20,20},{-20,20}};
      μ={{-mu,mu},{-mu,mu}};
 ℱ = fitness3D;
      fitPlot3D=Plot3D[ℱ[{x,y}],{x,-20,20},{y,-20,20},PlotPoints->20,Boxed->False,Mesh->False,Axes->{True,True,False},AxesEdge->{{-1,-1},{1,-1},{-1,-1}},AxesLabel->{"x","y",None},ViewPoint->{1.3,-2.4,1.5},
DisplayFunction->Identity];
 pop=initialPSOPopulation[𝕀,λ,μ,ℱ]
 ];
Dynamic[
If[running,
pop=psoStep[pop,λ,μ,ℱ,φ1,φ2];
ListAnimate[plotPSO[{pop},ℱ,Length[λ]], AnimatorElements->None],
ListAnimate[plotPSO[{pop},ℱ,Length[λ]], AnimatorElements->None]
]
],

{{running,True,"Run"},{False,True}},
{{φ1,0.1,"Exploitation"},0.0,10.0},
{{φ2,1.0,"Exploration"},0.0,10.0},
{{mu,1.0,"Velocity"},0.0,10.0},
{{funSelection, func1, "Function"},{func1-> "2-D", func2->"3-D"}},
SaveDefinitions->True

]

Particle Swarm Optimization (PSO) is an useful technique for parameter optimization in continuous, multi-dimensional search spaces. PSO is mainly inspired by the social behaviour patterns of organisms that live and interact within large groups such as flocks, swarms, or herds. The groups of particles are attracted towards each other, with the particles with higher fitness excerting their attractive forces.

This demonstration provides an exploratory platform for PSO and demonstrates its performance on two-dimensional and three-dimensional continuous search problems. A thorough understanding of the algorithm can be experienced by dynamically changing the values of the inherent parameters by the sliders provided.

In a Particle Swarm Optimizer search is performed by a 'population' of individuals i with location vector (Subscript[(OverVector[p]), i]), which represent a potential solution to the optimization problem. The criterion for optimization is provided through a fitness function ℱ. By applying ℱ to the solution vectors their relative 'fitnesses' can be calculated. Hence, the fitness of individual i is defined by ℱ((Subscript[(OverVector[p]), i])). The fitness of an individual can also be interpreted as its attractive force for other individuals. That is, the fitter an individual is the more 'followers' it will have.

In addition to its location (Subscript[(OverVector[p]), i]) an individual also keeps track of its 'velocity', (Subscript[(OverVector[ν]), i]), that determines in which direction of the search space and how far the individual will move in the next time step. The evolution of the overall population of search points is carried out in discrete time steps, where the potential solutions 'fly through problem space', where the flight paths are determined by both the inertia of the particles themselves and their attractive forces towards each other.

The following scheme gives an outline of a basic Particle Swarm Optimizer.

1. Initialize the particle population by stochastically assigning locations P=((Subscript[(OverVector[p]), 0]),…,(Subscript[(OverVector[p]), i]),…,(Subscript[(OverVector[p]), N])) and velocities V=((Subscript[(OverVector[ν]), 0]),…,(Subscript[(OverVector[ν]), i]),…,(Subscript[(OverVector[ν]), N])).

2. Evaluate the fitness of all particles: ℱ(P)=(ℱ((Subscript[(OverVector[p]), 0])),…,ℱ((Subscript[(OverVector[p]), i])),…,ℱ((Subscript[(OverVector[p]), N]))).

3. Keep track of the locations where each individual had its highest fitness so far:

	P=((Subscript[(OverVector[p]), 0]^best),…,(Subscript[(OverVector[p]), i]^best),…,(Subscript[(OverVector[p]), N]^best)).

4. Keep track of the position with the global best fitness: 

	(Subscript[(OverVector[p]), g]^best)=(max_((OverVector[p])∈P))(ℱ((OverVector[p])) ).

5. Modify the particle velocities based on the previous best and global best positions:
	
	((Subscript[(OverVector[ν]), i])^new)=(Subscript[(OverVector[ν]), i])+(Subscript[φ, 1])((Subscript[(OverVector[p]), i]^best)-(Subscript[(OverVector[p]), i]))+(Subscript[φ, 2])((Subscript[(OverVector[p]), g]^best)-(Subscript[(OverVector[p]), i])) for 1≤i≤N.

6. Update the particle locations:

	((Subscript[(OverVector[p]), i])^new)=(Subscript[(OverVector[p]), i])+((Subscript[(OverVector[ν]), i])^new) for 1≤i≤N.

7. Terminate if some condition is met (e.g., global maximum found or sufficient number of iteration steps).

8. Go to step 2.

Contributed by: Namrata Khemka and Christian Jacob

After work by: Khemka N., Jacob C., Exploratory Toolkit for Evolutionary and Swarm-based Optimization Algorithms, Proceedings of the Eight International Mathematica Symposium (IMS'06).