
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]

讀者回應 ( 0 意見 )
訂閱發佈留言 (Atom)
發佈留言
Please leave your name and tell me what you thought about this site. Comments, suggestions and views are welcomed.
如果這篇文章對你有幫助,那請留個訊息給我~