-
Notifications
You must be signed in to change notification settings - Fork 4
/
PSO.m
121 lines (108 loc) · 6.06 KB
/
PSO.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
(* ::Package:: *)
(************************************************************************)
(* This file was generated automatically by the Mathematica front end. *)
(* It contains Initialization cells from a Notebook file, which *)
(* typically will have the same name as this file except ending in *)
(* ".nb" instead of ".m". *)
(* *)
(* This file is intended to be loaded into the Mathematica kernel using *)
(* the package loading commands Get or Needs. Doing so is equivalent *)
(* to using the Evaluate Initialization Cells menu command in the front *)
(* end. *)
(* *)
(* DO NOT EDIT THIS FILE. This entire file is regenerated *)
(* automatically each time the parent Notebook file is saved in the *)
(* Mathematica front end. Any changes you make to this file will be *)
(* overwritten. *)
(************************************************************************)
(* ::Input:: *)
ClearAll[initialiseParticles];
initialiseParticles[objFun_,nparts_Integer,bndLo_List,bndUp_List,r_Integer:1]:=
Module[{ndims,positions,velocities,neighbours,idx,bestPositions,bestNeighbourPositions,vals},
ndims=Length[bndLo];
(*Print[ndims];*)
positions=Transpose@(RandomReal[#,nparts]&/@Transpose[{bndLo,bndUp}]);
(*Print[positions];*)
velocities=Transpose@(RandomReal[#,nparts]&/@Transpose[{bndLo,bndUp}]);
(*Print[velocities];*)
neighbours[poss_]:=Transpose@(RotateLeft[poss,#]&/@Range[-r,r]);
(*Print[neighbours[positions]];*)
bestPositions=positions;
(*for maximisation, replace Ordering[#,1] with Ordering[#,1,Greater]*)
vals=objFun/@positions;
(*Print[vals];*)
(*idx=Flatten@(Ordering[#,1]&/@Map[objFun,neighbours[positions],{2}]);*)
idx=Flatten@Ordering[#,1]&/@neighbours[vals];
(*Print[idx];*)
(*also:
bestNeighbourPositions=neighbours[positions]\[LeftDoubleBracket]#\[RightDoubleBracket]&@@Transpose[{Range[nparts],idx}]
*)
bestNeighbourPositions=Table[neighbours[positions][[i,idx[[i]]]],{i,1,nparts}];
(*Print[bestNeighbourPositions];*)
Return[{positions,velocities,neighbours,bestPositions,bestNeighbourPositions,vals}];
]
(* ::Input:: *)
psoStep[objFun_,bndLo_List,bndUp_List,positions_List,velocities_List,neighbours_,prevBestPositions_List,prevBestNeighbourPositions_List,vals_List,\[Omega]_Real: 0.5/Log[2.],c_Real:(0.5+Log[2.])]:=
Module[{nparts,ndims,bestPositions,bestNeighbourPositions,idx,\[Phi]1,\[Phi]2,newPositions,newVelocities,newVals},
{nparts,ndims}=Dimensions[positions];
(*Print[{nparts,ndims}];*)
(*replace < with > for maximisation*)
bestPositions=MapThread[If[objFun[#1]<objFun[#2],#1,#2]&,{positions,prevBestPositions}];
(*Print[bestPositions];*)
(*idx=Flatten@(Ordering[#,1]&/@Map[objFun,neighbours[positions],{2}]);*)
idx=Flatten@(Ordering[#,1]&/@(neighbours[vals]));
(*Print[idx];*)
bestNeighbourPositions=Table[neighbours[positions][[i,idx[[i]]]],{i,nparts}];
(*Print[bestNeighbourPositions];*)
\[Phi]1=RandomReal[{0,c},nparts];
(*Print[\[Phi]1];*)
\[Phi]2=RandomReal[{0,c},nparts];
(*Print[\[Phi]2];*)
(*velocities=Table[\[Omega] velocities\[LeftDoubleBracket]i\[RightDoubleBracket]+\[CurlyPhi]1\[LeftDoubleBracket]i\[RightDoubleBracket](bestPositions\[LeftDoubleBracket]i\[RightDoubleBracket]-positions\[LeftDoubleBracket]i\[RightDoubleBracket])+\[CurlyPhi]2\[LeftDoubleBracket]i\[RightDoubleBracket](bestNeighbourPositions\[LeftDoubleBracket]i\[RightDoubleBracket]-positions\[LeftDoubleBracket]i\[RightDoubleBracket]),{i,1,nparts}];*)
newVelocities=MapThread[\[Omega] #2+#5(#3-#1)+#6(#4-#1)&,
{positions,velocities,bestPositions,bestNeighbourPositions,\[Phi]1,\[Phi]2}];
(*Print[newVelocities];*)
newPositions=MapThread[#1+#2&,{positions,newVelocities}];
(*Print[newPositions];*)
(*lower confinement*)
{newPositions,newVelocities}=Transpose@Table[Transpose@MapThread[If[#2<#1,{#1,-0.5*#3},{#2,#3}]&,{bndLo,newPositions[[i]],newVelocities[[i]]}],{i,nparts}];
(*upper confinement*)
{newPositions,newVelocities}=Transpose@Table[Transpose@MapThread[If[#2>#1,{#1,-0.5*#3},{#2,#3}]&,{bndUp,newPositions[[i]],newVelocities[[i]]}],{i,nparts}];
(*Print[newPositions];*)
(*Print[newVelocities];*)
newVals=objFun/@newPositions;
Return[{newPositions,newVelocities,bestPositions,bestNeighbourPositions,newVals}];
]
(* ::Input:: *)
ClearAll[pso];
pso::usage="pso[objFun_,nparts_Integer,bndLo_List,bndUp_List,niter_Integer:100,r_Integer:1]
nparts: the number of particles
bndLo: a list of the lower bounds, one for each dimension
bndUp: a list of the upper bounds, one for each dimension
niter: Number of iterations
r: length of the toroidal neighbour scheme e.g. r=1 gives (i-1 mod nparts, i, i+1 mod nparts) as the neighbours of particle i";
Options[pso]=
pso[objFun_,nparts_Integer,bndLo_List,bndUp_List,niter_Integer:100,r_Integer:1]:=
Module[{positions,velocities,neighbours,bestPositions,bestNeighbourPositions,vals,iterCnt,globalBestIdx,globalBestPosition,globalBestValue},
{positions,velocities,neighbours,bestPositions,bestNeighbourPositions,vals}=initialiseParticles[objFun,nparts,bndLo,bndUp,r];
For[iterCnt=1,iterCnt<niter,iterCnt++,
{positions,velocities,bestPositions,bestNeighbourPositions,vals}=psoStep[objFun,bndLo,bndUp,positions,velocities,neighbours,bestPositions,bestNeighbourPositions,vals]];
globalBestIdx=First@Ordering[bestPositions,1];
globalBestPosition=bestPositions[[globalBestIdx]];
globalBestValue=objFun[globalBestPosition];
Return[{globalBestValue,globalBestPosition}];
]
(* ::Input:: *)
ClearAll[PSO];
Options[PSO]={"nparts"->32,"bndLo"->-10,"bndUp"->10,"niter"->100,"r"->1,"RandomSeed"->0};
PSO[objFun_,x_,opts:OptionsPattern[]]:=Block[{solution},
SeedRandom[OptionValue["RandomSeed"]];
solution=pso[objFun,
OptionValue["nparts"],
Table[OptionValue["bndLo"],{Length[x]}],
Table[OptionValue["bndUp"],{Length[x]}],
OptionValue["niter"],
OptionValue["r"]];
SeedRandom[];
solution
]