Stats

Popular Posts

Followers

Mathematica 教學:Particle swarm optimization

戴忠淵 於 2012年7月8日星期日 下午3:06 發表


Global Optimization by PSO: A Mathematica Program. Download  myPSO.m here.

PSOMax[fun_,cons_,varrange_List,OptionsPattern[]]:=
Module[
{k=Length[Variables[varrange]],
i=0,
dims=Length[Variables[varrange]],
\[DoubleStruckCapitalI],
iter,
lr=varrange[[All,{2,3}]],
vr,
initpso,
mypso,
\[CurlyPhi]1=2.05,
\[CurlyPhi]2=2.05,
TPPSO,
soltemp,
anstemp,
constrains,
method,
psotemp={},
tol,
std,
varobj},

Options[
PSOMax]={"Velocity"->
Table[{-(varrange[[z,3]]-varrange[[z,2]]),
varrange[[z,3]]-varrange[[z,2]]},{z,dims}],
"InitialPoints"->10dims,
"MaxInterations"->1000,
"Tollerance"->5*10^-3,
"Method"->"ConstrictionFactor",
"Parallelization"->False,
"PSOPlot"->False};

\[DoubleStruckCapitalI]=OptionValue["InitialPoints"];
iter=OptionValue["MaxInterations"];
method=OptionValue["Method"];
vr=OptionValue["Velocity"];
tol=OptionValue["Tollerance"];

constrains=(If[StringTake[ToString[Head[cons[[#1]]]],1]=="G",
Max[0,cons[[#1,2]]-cons[[#1,1]]]^2,
Max[0,cons[[#1,1]]-cons[[#1,2]]]^2]&)/@
Range[Length[cons]];


varobj=Transpose@{Evaluate@varrange[[All,1]],_Real&/@
Range[Length@Evaluate[varrange[[All,1]]]]};

TPPSO=Compile[Evaluate@varrange[[All,1]],
Evaluate[Simplify[fun-10^8Total[constrains]]]];

initpso[pnum_Integer,locationRanges_List,velocityRanges_List,fitFct_]:=
Module[{locations,velocities,fit,pb,gb,bb},
locations=Table[(Random[Real,lr[[#1]]]&)/@Range[dims],{\[DoubleStruckCapitalI]}];
velocities=Table[(Random[Real,vr[[#1]]]&)/@Range[dims],{\[DoubleStruckCapitalI]}];
fit=Apply[TPPSO,locations,{1}];
pb=Transpose[{fit,locations}];
gb=SortBy[Transpose[{fit,locations}],First][[-1,2]];
{locations,velocities,fit,pb,gb}];

mypso[temp_List]:=
Module[{bb,pb=temp[[4]][[All,2]],gb,pbtemp,gbtemp,vrtemp,
locations,fit,fittemp,particle=temp[[1]],
pbest=temp[[4]][[All,1]]},

fittemp=temp[[3]];gbtemp=TPPSO@@temp[[-1]];

gb=If[gbtemp>=Max[fittemp],temp[[-1]],
SortBy[Transpose[{fittemp,temp[[1]]}],First][[-1,2]]];

vrtemp=MapThread[If[method=="InertiaWeight",
((0.4+0.5*(iter-i)/iter)#2+\[CurlyPhi]1*Random[]*(#3-#1)+\[CurlyPhi]2*Random[]*(#4-#1)),
2/Abs[2-\[CurlyPhi]1-\[CurlyPhi]2-Sqrt[(-4+\[CurlyPhi]1+\[CurlyPhi]2)(\[CurlyPhi]1+\
\[CurlyPhi]2)]](#2+\[CurlyPhi]1*Random[]*(#3-#1)+\[CurlyPhi]2*Random[]*(#4-#1))]&,
{temp[[1]],temp[[2]],pb,Table[gb,{\[DoubleStruckCapitalI]}]}];

vrtemp=(MapThread[Min[Max[#1,First[#2]],Last[#2]]&,{#1,vr}]&)/@vrtemp;
locations=MapThread[#1+#2&,{temp[[1]],vrtemp}];
locations=(MapThread[Min[Max[#1,First[#2]],Last[#2]]&,{#1,lr}]&)/@locations;

If[OptionValue["Parallelization"]==True,DistributeDefinitions[TPPSO]];

fit=If[OptionValue["Parallelization"]==True,
Parallelize@Apply[TPPSO,locations,{1}],
Apply[TPPSO,locations,{1}]];

pb=(If[fit[[#1]]>pbest[[#1]],{fit[[#1]],locations[[#1]]},
temp[[4,#1]]]&)/@Range[Length[fittemp]];

N[{locations,vrtemp,fit,pb,gb}]];

soltemp=NestWhile[(i++;psotemp=Insert[psotemp,#[[1]],-1];
std=StandardDeviation[#1[[3]]];mypso[#1])&,
initpso[\[DoubleStruckCapitalI],lr,vr,TPPSO],std>tol&,2,iter][[-1]];

If[OptionValue["PSOPlot"]==True&&dims==2,{TPPSO@@soltemp,
Apply[Rule,Transpose[{varrange[[All,1]],soltemp}],{1}],
ToExpression["Iterations->"<>ToString[i]],"Standard Deviation->"<>ToString[std],
Manipulate[
ContourPlot[fun,Evaluate@varrange[[1]],Evaluate@varrange[[2]],
Epilog->{PointSize[0.02],Red,Point[psotemp[[z]]]}],{{z,
1,"Generation"},1,i,1}]},{TPPSO@@soltemp,
Apply[Rule,Transpose[{varrange[[All,1]],soltemp}],{1}],
ToExpression["Iterations->"<>ToString[i]],"Standard Deviation->"<>ToString[std]}]
];

Usage: PSOMax[function, constraints, variable range, options]

Show[
Plot3D[x^2-(y-1)^2,{x,0,5},{y,-2,2},PlotStyle->None],
Plot3D[x^2-(y-1)^2,{x,1,5},{y,-2,2},
RegionFunction->Function[{x,y},x>1&&x^2+y^2-4<=0], Filling->Bottom]]

AbsoluteTiming@PSOMax[x^2-(y-1)^2,
{x>1,x^2+y^2-4<=0},{{x,1,5},{y,-2,2}}, "Velocity"->{{-3,3},{-3,3}},
"Method"->"InertiaWeight","PSOPlot"->True]

AbsoluteTiming@
PSOMax[x^2-(y-1)^2,
{x>1,x^2+y^2-4<=0},{{x,1,5},{y,-2,2}}, "Velocity"->{{-3,3},{-3,3}},"InitialPoints"->50,
"MaxInterations"->1000,"Parallelization"->True,
"PSOPlot"->True]

Options[PSOMax]

Tags: ,

讀者回應 ( 0 意見 )

發佈留言

Please leave your name and tell me what you thought about this site. Comments, suggestions and views are welcomed.

如果這篇文章對你有幫助,那請留個訊息給我~