From b79334b61f2fad32941edb6eea4eb28b4ab4bb21 Mon Sep 17 00:00:00 2001 From: "Federico E. Benelli" Date: Tue, 5 Mar 2024 13:41:15 -0300 Subject: [PATCH 1/4] minimal fpm structure, and a more versatile interface --- Fortran codes/PSOclassicG.f95 | 364 ------------------------------- Fortran codes/PSOclassicG_output | Bin 27016 -> 0 bytes src/PSOclassicG.f90 | 324 +++++++++++++++++++++++++++ 3 files changed, 324 insertions(+), 364 deletions(-) delete mode 100644 Fortran codes/PSOclassicG.f95 delete mode 100644 Fortran codes/PSOclassicG_output create mode 100644 src/PSOclassicG.f90 diff --git a/Fortran codes/PSOclassicG.f95 b/Fortran codes/PSOclassicG.f95 deleted file mode 100644 index a72b07e..0000000 --- a/Fortran codes/PSOclassicG.f95 +++ /dev/null @@ -1,364 +0,0 @@ -program PSOclassicG - -! Author: Mohammad Asif Zaman -! Fortran code (gfortran, G95) -! The fitness.f95 and timestamp.f95 files are expected to be in the same directory. - -! gfortran -c PSOclassicG.f95 -! gfortran PSOclassicG.f95 -o PSOclassicG.exe -! ./PSOclassicG.exe to run (in Cygwin or powershell terminal) - - -! Date: 13th Aug., 2012 - -! VTR = Value to reach. Required fitness function value for termination. -! NIR = Number of iterations required to obtain VTR. -! pbest = population best solution at each generation/iteration -! pbest_val = fitness value at pbset -! gbest = global best solution -! gbest_val = fitness value at gbest - - - -! -> Mathematical parameters -!real(kind(1.0d0)) , parameter :: pi = 3.141592653589793 -!complex, parameter :: j = (0,1) -!real(kind(1.0d0)), parameter :: eps = epsilon(eps) -! <- Mathematical parameters - - - -! May 2, 2020 -! version 1.2 -! - timestamp subroutine: date bug fixed -! - added Makefile -! - The make file works in both linux terminal and windows powershell as well -! - The code works fine in both linux and windows - - - -! -> File operation parameters -CHARACTER(len=18) :: str -character(len=29) :: command_Make -character(len=35) :: command_File -! <- File operation parameters - - -! -> PSO and fitness parameters -integer, parameter :: psize = 30, dimen = 30, max_iter = 5000, c1 = 2, c2 = 2 -real(kind(1.0d0)), parameter :: VTR = 1e-6 -integer :: iter = 1, overfly_counter = 0, NIR = 0, NFE = 0 -logical :: overfly_flag = .false. -real(kind(1.0d0)) :: xmin(dimen), xmax(dimen), vmin(dimen), vmax(dimen) -real(kind(1.0d0)) :: x(psize, dimen) = 0, v(psize, dimen) = 0 -real(kind(1.0d0)) :: pbest(psize, dimen), gbest(dimen) -real(kind(1.0d0)) :: pbest_val(psize) = 0 -real(kind(1.0d0)) :: gbest_val -real(kind(1.0d0)) :: w -real(kind(1.0d0)) :: gbest_val_store(max_iter) = 0, x11_store(max_iter) -! <- PSO and fitness parameters - - -integer :: loop1, loop2, loopG - - -real(kind(1.0d0)) :: time_start, time_end -real(kind(1.0d0)) :: temp1, temp2 - - - -! The following loop runs the PSO. Every time this loop starts, a new PSO run is initiated with new seeds. -loopG = 0 -! Do while (NIR == 0) ! Turn this loop on if you want to run PSO multiple times till it satisfies VTR - - ! until NIR reaches a finite value. If VTR is not achieved, this loop will continue. - - - loopG = loopG + 1 - print *, loopG - print *, NIR - - - NIR = 0 - NFE = 0 - overfly_counter = 0 - - call cpu_time(time_start) - - - ! -> Defining solution space - ! All dimensions are assigned same limit. But this can easily be changed by employing a loop - xmin = -30 - xmax = 30 - - vmin = 0.5 * xmin ! lower limit for velocity - vmax = 0.5 * xmax ! higher limit for velocity - ! <- Defining solution space - - - - call random_seed - call random_number(x) - call random_number(v) - - - ! -> Intializaing particle position and velocities - - ! The cost of the first particle is evaluated because it is needed as a stepping stool for evaluating gbest. - x(1,:) = x(1,:) * (xmax - xmin) + xmin ! fixing the location of the first particle - v(1,:) = v(1,:) * (vmax - vmin) + vmin - pbest(1,:) = x(1,:) - gbest = pbest(1,:) - call fitness(dimen, pbest(1,:), pbest_val(1)) ! evaluating the fitenss of the first particle - NFE = NFE + 1 - gbest_val = pbest_val(1) - - - ! Composite threaded loop (15th Aug., 2012) - ! As the population is adjusted, the corresponding fitness parameters of the population is assigned simultaneously. - - do loop1 = 2, psize - x(loop1,:) = x(loop1,:) * (xmax - xmin) + xmin - v(loop1,:) = v(loop1,:) * (vmax - vmin) + vmin - - pbest(loop1,:) = x(loop1,:) - - call fitness(dimen, x(loop1,:), pbest_val(loop1)) - NFE = NFE + 1 - if ( pbest_val(loop1) < gbest_val ) then ! Geared for minimization problem - gbest_val = pbest_val(loop1) - gbest = x(loop1,:) - end if - - end do - - - ! <- Intializaing particle position and velocities - - - - - - ! Main iteration loop. Creating the generations. - do iter = 1, max_iter - - w = 0.9 - (0.9 - 0.2) * iter /max_iter - - gbest_val_store(iter) = gbest_val - x11_store(iter) = x(1,1) - - - do loop1 = 1, psize - overfly_flag = .false. - do loop2 = 1, dimen - call random_number(temp1) - call random_number(temp2) - v(loop1, loop2) = w*v(loop1, loop2) + c1 * temp1 * (pbest(loop1,loop2) - x(loop1,loop2)) - v(loop1, loop2) = v(loop1, loop2) + c2* temp2 *(gbest(loop2) - x(loop1,loop2)) - - if ( v(loop1, loop2) > vmax(loop2) ) then - v(loop1, loop2) = vmax(loop2) - end if - - if ( v(loop1, loop2) < vmin(loop2) ) then - v(loop1, loop2) = vmin(loop2) - end if - - - x(loop1, loop2) = x(loop1, loop2) + v(loop1, loop2) - - if ( x(loop1, loop2) < xmin(loop2) .or. x(loop1, loop2) > xmax(loop2) ) then - overfly_flag = .true. - !print *, "fly" - overfly_counter = overfly_counter + 1 - end if - - end do - - ! Invisible boundary condition. Can be changed to reflecting or absorbing boundary conditions. - - if ( .not. overfly_flag ) then - - !print *, "No overfly" - call fitness(dimen, x(loop1,:), temp1) - if (NIR == 0) then - NFE = NFE + 1 - end if - - - ! Composite IF statement that adjusts the pbset and gbest parameters (15th Aug., 2012) - if ( temp1 < pbest_val(loop1) ) then - pbest_val(loop1) = temp1 - pbest(loop1,:) = x(loop1,:) - - if ( pbest_val(loop1) < gbest_val ) then - gbest_val = pbest_val(loop1) - gbest = pbest(loop1,:) - if ( NIR == 0 .and. gbest_val < VTR ) then - NIR = iter - end if - end if - - end if - - - end if - - end do ! Ending loop over population - - end do ! Ending main iteration loop - - - call cpu_time(time_end) - - - - - - call timestamp(str) - - - !do loop1 = 1, psize - !print "(4f10.3)", x(loop1,:), pbest_val(loop1) - !print *, x(loop1,:), pbest_val(loop1) - - !print "(1f10.3)", pbest_val(loop1) - !end do - - !print "(2f10.5)", x - !print * , "Best" - !print *, x(1,:) - !print "(4f10.6)", gbest, gbest_val - - print *, "======================================================" - print *, "Time stamp = ", str - print *, "======================================================" - - print *, "======================================================" - print *, "Algorithm parameters" - print *, "--------------------------------------------------------------------------------------------" - print *, "Problem dimension = ", dimen - print *, "Dynamic range = ", xmin(1) , "to", xmax(1) - print *, "Population size = ", psize - print *, "Maximum allowed iterations = ", max_iter - print *, "Defined value to reach (VTR) = ", VTR - print *, "======================================================" - - print *, - - - print *, "======================================================" - print *, "Performance parameters" - print *, "--------------------------------------------------------------------------------------------" - print *, "Total CPU time = ", time_end - time_start, "seconds" - print *, "Overfly counter = ", overfly_counter - print *, "No. of iterations required to reach VTR (NIR) = ", NIR - print *, "No. of function evaluations (NFE) = ", NFE - print *, "Global best value achieved after NIR = ", gbest_val_store(NIR) - print *, "Global best value achieved after max_iter = ", gbest_val - print *, - print *, "Global best solution =" - print *, - do loop1 = 1, dimen - print *, gbest(loop1) - end do - print *, - !call fitness(dimen, gbest, temp1) - !print *, temp1 - print *, "======================================================" - -!end do ! Ending NIR loop if active. - - - - -! -> CSV File operation 17th Aug., 2012. (M -> 18th Aug., 2012) - - -command_Make = ('mkdir ' // 'Data_') // str - -print *, "commandMake = ", command_Make - -call system(command_Make) -command_File = 'Data_' // str // '/' // 'converg.csv' ! the file name must be 11 characters, including . and extension -open (20, file = command_File) - do loop1 = 1, max_iter - write (20,*) gbest_val_store(loop1), ",", x11_store(loop1) - end do -close (20) - -command_File = 'Data_' // str // '/' // 'final_x.txt' ! the file name must be 11 characters, including . and extension -open (30, file = command_File) - do loop1 = 1, psize - - write (30,*) x(loop1,:), "," - - end do -close (30) - -! <- File operation 17th Aug., 2012. - - - -! -> Log file operation -command_File = 'Data_' // str // '/' // 'logfile.txt' -open (11, file = command_File) - -write (11,*) "======================================================" -write (11,*) "Time stamp = ", str -write (11,*) "======================================================" - -write (11,*) "======================================================" -write (11,*) "Algorithm parameters" -write (11,*) "--------------------------------------------------------------------------------------------" -write (11,*) "Problem dimension = ", dimen -write (11,*) "Dynamic range = ", xmin(1) , "to", xmax(1) -write (11,*) "Population size = ", psize -write (11,*) "Maximum allowed iterations = ", max_iter -write (11,*) "Defined value to reach (VTR) = ", VTR -write (11,*) "======================================================" - -write (11,*) - - -write (11,*) "======================================================" -write (11,*) "Performance parameters" -write (11,*) "---------------------------------------------------------------------------------------------" -write (11,*) "Total CPU time = ", time_end - time_start, "seconds" -write (11,*) "Overfly counter = ", overfly_counter -write (11,*) "No. of iterations required to reach VTR (NIR) = ", NIR -write (11,*) "No. of function evaluations (NFE) = ", NFE -write (11,*) "Global best value achieved after NIR = ", gbest_val_store(NIR) -write (11,*) "Global best value achieved after max_iter = ", gbest_val -write (11,*) -write (11,*) "Global best solution =" -write (11,*) -do loop1 = 1, dimen - write (11,*) gbest(loop1) -end do -write (11,*) -write (11,*) "======================================================" - -close(11) - -! <- Log file operation - - - - -contains - - include 'fitness.f95' - include 'timestamp.f95' - - - - - - - - - - - -end program PSOclassicG diff --git a/Fortran codes/PSOclassicG_output b/Fortran codes/PSOclassicG_output deleted file mode 100644 index 58f0de84da7e2bfe306330265d01ae7275ea56c1..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 27016 zcmeHweSB2ang7kpNW_rD*9EaUDrmGaCcH$zRc3&J8=63bK-m_Y49O%!lT11@iD0Wn zldz6)R(@S*q~>JI8ns-Qbkc_zt4HOZ*zxB zyY%zfKYo`_=A84K=lMR*dCqgsIrrR|!(IN0CD~b7ioCLw8x<<;$oDc5|5TKHK~@yB zxyqS}SGh#FNI3%(7yi7ghEkrB=`jH<(>%eKgESX^bF(>NZnl?YWO@-NLj%90g_M#$ z{ql+|Pb(^TQKpLGW(1NWTgNL|5ory}$nH>0VsrESxUO>@g;l+E@;!=9O}-GooPOP8-^nxr4y zD3Q(Q;g4c~>W{AUT%t_zt^KI=i4U*boBP20x?64`S^0YlS;NOse-jo5ym%+};J#_74@Y8=U<=(Xvy?!2 zb!DI~6bUsnMq{B!b>*U_a7(B%z@}$Xo+0Xh*CKx6}q>uqM*j(qLeY#scBi zP>X@0DQs+sg&LsPs#9YKMk2w@Y1Kev5`j9n!hjpy3>P&AqHP-tL=CHUBh0Ra57Y4K zg0WBl%?!jEn@xDJKy6bvY9xtTG8!?=XQ8$;3z%&lKQ0DU~(fB{& z$UzO~lhiDX#@B+G#&tLT+AeH7lO)#?0In<`A>>(o3DSMzQ>V!-eE9WTR(D<{I zDWZ<#IbODgcWgxsk+t;&E8t^tRSnBdz$fXGm)k%WDp*yqfll*{yh;poIfqi+TmzkS z$jfV>%Q>FvR0ACj(XJ{3{Va{DC~FP$Nd|h2fnH#sZ#2*g4fIw6UC!lHy~#j7+rZyp zpc}{Sdkyq+4g6aTbegN>)oY-iFNsL|4D<^O^oI=ezcSEw80eD?^qmI!6a#(GK)=vH z-)EptHPA;4bh(D4YPIu7p4#L3VzFCMy9Z;g`&IhPUP5%g=g_n4d zqw+?o92-sHn*KJGDfGrhSosYqQ^<`CvhsgXnL=%B2P?loWeTyeK30B~$}TEzW#wH| zrls3h2P^kenL=-@m6ac%GKJh&4J-dYDpRP9Rk8B_rSe27ds+D(s7xU@R>I2PrgA=& z-K=~cl_~Vbidgv`DpSafDXe@ul_}K5#!rASJwjy)v9VEBzJWpzDYUi;6 zHU4J%)SU!H-GBD_r>G(|3g#|`zkWCMdc8-M2dI9vdiPNVFmxgZ0BU?(9ZD`xv-Ya{ zPsGmG(9V)*g_kJncXRz6H#`XYmA122LonF_hpyl2+5%G6hfLE`3LYrb&IL-F>z%)c zoCfwj38Y5WkKGO#%CA!IJ~9_^rN*B;_<1T7fAl4MR>ZgOL-I@3pw8oG-Q`hx@?YyK z_bO^)^;R{0P|X@tv+z06Pa&zspI76<$t&R>wR3mYJHJwgyT>7aCfx+X-O3lK6h$_& zdB{@GA$trAjGE}C2vp%>0C!Cu;hMU$|r%=U^>^z#OuWw@|l^#v_Q*ZD%r&!JNdvh_iu@lDV)R!N6@$Z1kUI z;H`QE8ThSZ%!7j&4ZMOIXpY3wVBm*PY-`}NdIcHy$|uagj*JFg#SJtE<7qH(rEZ70 zi?7lv$iVq9u&29@2Dx4}QH4fTVj>w~6A3-tM#7F$J_(}f$V`AD)OyE_}Jm2GYo8XC!*lBH){5fN< z9E1q*LWEXJ%z6n=e8f*{rDPX}CnoT^ZQzNG*l8)1T+QJ-08hLSDY3#6AMw-DDmk6; zCx}Vy>8AO=pXg%p`YLL-sy*wEJ+oV`6ndx8Ot}rM1)I{ta-D-&ahAcb()&s(bq|_9 z+}(%VSHLsey%o6*mfMCLElFV)5h*Or9Z8&gTk}B4z#`wW;p1b-`Q|YkX$vtOXs?V-y|yG-Txfn5b&!r3HWb@7DdvTVNlOfZ+pq- zyz^jy>}|N6dvg9?NNr8?N!9ND(Aa^k%;{M#sPX;k(4k!Qs{I<_Wp(H<2rsMg58%=K z1FSi`ySoms4&;Sxtj|9vT(^t*Nd(0Gg1DACF84VmOHpy3-f@3|YStvFn)RksO%ai% zniGX;CTmfBgVsogR_C%-Q|vw0N_PR_JgGOKL`|sZSB&lCL9@efD?oyY_R+mgqR%26 zD9mWX=A>c9A}-O*bjJD6VHZfr!|ww*z8g0w-7n?wUHu(#^nU#Hr?z2CN>0HXh+&Qm zJ@XF6CaK-8#PZbm?~<>NV{URamAMs|6?fGlLwC(DS|~lcyB7uz`;X87-Fe56iwhrk z9oZiL5w&)NJ~5p|4pk()hjdSDQscXr6Eug$pIw&7&z9uhs6%fdn0}+iKU1?#V?5n+ zx~)1bp>}es{L(3>>MOe_e3W;@pAt`<$t-nfG)Hx@yYXaU^}u^9thfqIjD=QWJC8KP z`z|fc{vC@!+c;}KN^zDuTg2J%S#B0ay(&OyP}xJ6$gqg3FZ2wx*Jw&X&K31MB6&@s z)ev?w*$pT{FD=Fle$#HjeJr524HiB&*y%r#)jine8oLb7ti!I0VH^$mJ+8^f-Ei-@ zicM7P&S`&ZC<{7IOCy$twl=ndSFySP`gu>Pw$q}LAd*|7`FjNaoaq*9vASm4jsRj0KM z;ytUAYV8YZ=i7AObR8fu7AOvW;fnKF%Z-Y=_Ru8xXDH@#E5g@7@u@^pKI>||-5FfY zb!m}6;i1de`3fE~LUUJl?V*8-T726ANAX;PWwNJe$-^8rW!N?GhnI@6;+kv4*mgN( z?!f&h#+D=>TWax$(c^l(?^3069F)SY31TYk{gc)LFya))?ia`8l7I~ew&QgcMiJZb zdT)}ft@5}P42*?6u4Txv6@m{LyOFslL0ZTygp@D`vjPg{IartmMr$pXg z=~GMLmisYs8V&p(PA_&VynQ{cb7nB_U4|%3&_W8S%YzA$#ZG>YtN$9e5_cWJr_EF} z1=~bDuJ59fPe3e~dd}ZM+`F>k=ih-0o3k4DcFkPl0H256!=PV9G_BtU&Oxi`Sx=$B zus%y3rpe^w53C#BVyLu;>@NfRKrfb-1UT=oRzbyORJ7v1<&2+gz@J3;&kOwDzDf8A z-W+D5@mc+std_5Bq5W4>v}#W}YwywH|p_}829E9~&!b;i#%;4}Lxmh$L&<@cI@ z_!b5|(NLW2GAn+O_^b^(>a4g&hv}y+qJbKrcmq?MpxFX96AOFx0*AlynPw29=c9pE zA@*f1#vhE+#NO3gLv6U&&gl1?jec2&VfVs7L}>phyXU#uTrx{~Ynv`XihoReR&R|v zD_*a|^nZ!)yM*G|R>eOO6;pKyQhbe_;`f~u|G|KN2`PR^;QtaUAY=aw>=PCD=n|y( zKZwuT;t!k^->Soy2Zn{ZV!~f#hyRf?{#iOc!Q|V)eBuoAv~G>r%=K5$sr z|4jg5YVM{Pl&0zumae1ys_u>K-cRr(?Owjr07>6382Cg$7Mo3;XJ_(pXOrL8Vfamc zLg4Q+3`y)JFUC@sg6C(t1QiE~&+4nsofU7=Vfss`VMQt*a#vUtyM^L1U4j%ZwNw0s zv*MHie;O(F3dJwKrUe0?1f~haPwEn+_<7>9w)ljz;w?H%|COY8qfopafIY6Im%GK2 zJ6n7w?xkmTJ~wBEQJS=@*&JsDv3d%lT-^rZ3KLHs&bcgi`X= zd;etvL*F%V=!<{RVa$`modl>4t`D%{cCpg-BAW5#Am7e&;L68eK=<|2*P!Sf8cqGl zL&V+FO>Ao1kN0KWD$n*JTg9_|J>6@Oi~Aoce5`vNa`*ZlPV1HPom=pxZfJi2ln=Zn zqUnj>QVZ}#RV>lSGI2NlkX31c_6RHbxAJuh*l5V)uV5ya+O5D$z1{CnyIqWX?^EfH>PF&wd>mZ*i&bJHH%%t?e$GyO!bsoa?=yDp&HeqvM(mdpdk?XTz52L+ z2eR~P5f0voEWI>9b`V*53x({y!e6l0RP}NH2weOKsQ3C`NOSQNxL8d9F5XRVr-r+q zz$fYBk;8`tdP_-bZ$m^hAJlHEyYD!{_d-XA#X7J|wdq3()}?5f_Sk*$fpjm_${q!< zGVrMvZt+#a=MM*VnzQMGfDz zrU&q}khxrQafjw&3(!#j@sJtrZbfv87iJyEQdj`;USyfjaQ7yewQiEtpaQHLkzkbx ze?H16e8ga1uj}V=Gn~8u3S}&j1}*Ta4`>w>O(7I*#`(#1UE2r_wb^w}4?591HJaED(4M5B$eM&m=BCN6AoF$q^ zJ2a14`wF*F%6h{`c+hOP%qn4fv+rsxLGAOuJQ+%-S!u(>~K#`^!3vxphB% z!PL3|Dq6iqKS8oC7%bjv&xBuT#=qGPe~vT$1qOV!+aq_Y-raAq{{$7S_Ln*1Kc<`6 zZ+e>epO28^GVNBi z8^lbzRd;Z)8vXfA42A4_%?orfGP;$Rt-UeN+2~>&#=KjV`nA^EW^Wy4W@$^ePl(T| zc)qjZpXo5>-Ks5C#d5c*S(l&|-)X1#dS}IDI!wQ5w`!6>G0i#JZq-YA1u1sfDgJ9` z#ZT%m=7C|`^O~y(|Ho9cdh|wT{Pj9M!Q5;I<8y|&K!-7#`R;QjGfz;_YUV;`{Kxdx z@`q2kTa_^5|DzrLB4_+H2K@EZjg!m0^f-~6VaC7F4u7#T{>M6udAI7hUuE8{+QcUK z1W(d-tD*)-({9!6fHX{qJyf*1-0y7iA_G1>Fb>Gwsta{|ev`}Hsw`cCTz8J0;w8?C zpEcl{cB_8)pIXDXV!2y&uP#A~A0$3&iQ%uzEDW=NV7_S-K2lTXny2%;lQQlD(M5QNr!<}4ny@QFyo6{L17^XlR^ymt3%2B zy_&VI8Uw>(!9WlXLB&BR?IWdq$=~i_bKH8UUw06e?Iyr!OqM)UG|4UU5ZeiPWh;LB zwpR2ceveudzvFZCC7SJpOTc{6NK3Z?fOK4te_4LAH*S2eV` znL76kZeLSF7)Kp!Y<9QeT!7{fer+FBX8Z-Osz?}zHZ;5IU`0z52L-^0#hY7#&5gBg z94pWeLOB+$3b(d31?e0CceL?#qE!YrH8!_3yMs+l;oCxWZXC}LVGL1J!WXoSxj0nc z*n+zDU{hNNC~h3(P`lAxyrz2PRh&!d`(*Eq8>>PQoG#HEY^e2Uf(zR$e8r|TjjkaF{m*V(^O`h1MSW~#6zOgAp zSr(Xo!7!OIV_C3eNoYf5Fyd>CR0cOMYin87*5qqzSQTnrQ5#zxZeJX#Wme#a!&s=( zwMnGk6zQEJy<4PRBK=#D-Y3$8NVkgg??w6@kuC;yDu=0HADE74b6vaWX;8mz9OU^M zy*wRVFl*MV1&a6k{`&g*1)Rp?5c!(pmrfzRz6T$CM%D=yr9)nl*8hs}<`RH_%Ldsiy85$OfHQ>j6uUaa7fNH6;|l`4XluBJq4C1oqR}{>>jDPegOvZ#=P_(3A z^0LCSZp-UX7F>A!)w8BuMi9hDK3xu*uck)$3W~O5FPe0Q+L@E{r95Z{4Hxk>;jan{ ztv6;MW8^z9!^F1@)f11UQqPiJ$>#^&V3x*5ZFvHJ55p(dBh~oYISbkJD*i@LAENq& z1w{{JFD#h+cRBt7cOrLTLGk@AwV-6n8EV1Y&Iy$T-ble*UqOklpm|qX_-iuh;)re>qXiw(k&wW zj!3QTq4-_8XwiIk@#+m&>$JJ&cxHP_W|XzDQt6$gWuB7Rp3<3D@e23MlG3u0vXXgL z9J_K3?vAEWcey%ny8HI8%yUFp$OvmW4yv9gk`5ff}?L0_hWe!FQd483dj?qAjH)i@2 z<)AGcqa?mKX1#J0IhyJ<6yLeZ`zf8mb2Lt8De3mp_?!*@k6QYLkcsZL%nu$xx6T{1 zuE?U%T-Yh+iB*iA?$6r=j+}>3WY-paW-00Y@?+5L+Vv>tWam$WoouZI+S3fDNTDNe z=&E6#it-#%;vYHR%P4Hkh1^T{Bl~Hum9AX&sVJ`@B|2?|(k1Ei!vmtzrZ8QSejItC zdj;JsCN%nH7}4oxM|4R#eG7`{^b-@h#@Q$43((2VM;3V*g{}FJo6h*BC`l0~Qt!1A z=e`=o@BKbY6xy^WeEPBpU6M}wO+;6AbGpO{F*;t+iTh<5hXJp>|Jxb=If`206`VZj zV)XQQqguOuWNT?9{r>%DhLe83JWd87POSGS_fvZOyvA_S@6#$wyly4^zB>&5b4Z2d z`XmFr80&>|&|h>467cK)ot!LbV$CxBC(0`Z#XC|Yk zkE_EO_z!2Ge*(IlJ?ElbXXErgIX=s{x{S5gXgSZ%X8up7F9*HIreD6E0q3p^^u7%A z-5Ka(8R!$y5W9Yy3i{dEQxrND%y?ayfxlYtk6QSHjK4sk zpdZRW&&Gnf2;;yJi#-LP+qJhO1AR>f`X1abgrB!N$l4`8HZ#3b()4y@!1)(OFHq=Q zH@dE8pLovAz<(=iug&_THWG`*+Un~)wKzg^*`k$!it<&}I8#+Wc{LDgrqe)MLO7%o z2ag6C*p^kGE*6eN1HraUIGD4!wF%pZb)M_y%)`;3c1!{6DK+AxPde2!)Dnwq#_^ok zdkoaIH8*dDh@lh!Cmjl^vs$-8_07{+v*wm6K2M2f*6fm*N~N!SIZSV9jG$aLPid(S zQD&}!143y(4;x=>KM%y~N|a4V=FU^vDVeFX(iU6*yL41JN2za&p}|oP6qG5nuLg42 zY^6bR;h@oQgrUz;HkFoY#d)QgvihcA1A&3Uk4N>Cmd^9cx(q*10uJ9V2a%^qv;T&J> zNLkZiwP`?>(_c|aI7vU2B!rAICEKssAi`YFNPC}O-O72 zh=v0jvBS)cYSsdW0@mOa4M{(J*#aca_O+OmZi?wNWebKL>1a0#L=yudHLxKXP3to4 zoM$VyrrmT*aYiovz-TK=h(g8S4)X!a7H)Q+vECH@C}j(u;)!ll^D750y3sqZ<`Z-R2hH@iZ3UH5f|i9ge9G zKwTu%6eIzWYi)`t9@gU?WIYXG6k{R0Mexu}0H!c&zbCX&jQtzy03u6VlurU&LgsI& z1sl}R+auTzjuAX4>s?;1+b5V45~vr6?KY<^D)a08TS;>(5=1ycijrP!^W@d$YZ)s#g%C z{a7}zOYYZ+qOu$r+B1~$a=&G#Ajp*V9jWb79#1IjT7yq4AvJlqFViX%%6&=Bo$kMz zQ9*lxQor24IZ95!HJ>NiKB(lA>8+@zJw{oU`#ifw9k=O|4&zZUeo^;tI(`_gj<@@rztiDfdbW3{~5g*p_4?r`@%jdftLjE4hPNLK+ z@&4H&FZa=BFX0?#N=A`d+yAdZUiP1S4)Y56DvN;Ig7@zhdAXmqbs1-@u<%>p*y+C1 zeyLyXyNwEY`hPHtm(_pIpt1me6q|I(_XhI+C9(2b%LgD+h(9SW_Xmzu5?XqHS>+Fb zhfI|6+I}Am47jZP*79NSQ0z*1xvwt&-mGPC%?w!e=P@2US#@2mL>lUR@__- diff --git a/src/PSOclassicG.f90 b/src/PSOclassicG.f90 new file mode 100644 index 0000000..f5637ab --- /dev/null +++ b/src/PSOclassicG.f90 @@ -0,0 +1,324 @@ +module PSOclassicG + !! Author: Mohammad Asif Zaman + !! Fortran code (gfortran, G95) + !! The fitness.f95 and timestamp.f95 files are expected to be in the same directory. + !! Date: 13th Aug., 2012 + + !! VTR = Value to reach. Required fitness function value for termination. + !! NIR = Number of iterations required to obtain VTR. + !! pbest = population best solution at each generation/iteration + !! pbest_val = fitness value at pbset + !! gbest = global best solution + !! gbest_val = fitness value at gbest + ! + !! -> Mathematical parameters + !!real(pr) , parameter :: pi = 3.141592653589793 + !!complex, parameter :: j = (0,1) + !!real(pr), parameter :: eps = epsilon(eps) + !! <- Mathematical parameters + + !! May 2, 2020 + !! version 1.2 + !! - timestamp subroutine: date bug fixed + !! - added Makefile + !! - The make file works in both linux terminal and windows powershell as well + !! - The code works fine in both linux and windows + use iso_fortran_env, only: pr => real64 + implicit none + + ! -> File operation parameters + CHARACTER(len=18) :: str + character(len=:), allocatable :: command_Make + character(len=:), allocatable :: command_File + ! <- File operation parameters + + + integer :: loop1, loop2, loopG + + real(pr) :: time_start, time_end + real(pr) :: temp1, temp2 + +contains + + subroutine pso(fitness, xopt, verbose, psize) + interface + subroutine fitness(x, cost) + import pr + real(pr), intent(in) :: x(:) + real(pr), intent(out) :: cost + end subroutine + end interface + + real(pr), intent(in out) :: xopt(:) + integer, intent(in) :: psize + logical, intent(in) :: verbose + + integer :: dimen + real(pr) :: xmin(size(xopt)), xmax(size(xopt)) + real(pr) :: vmin(size(xopt)), vmax(size(xopt)) + real(pr) :: x(psize, size(xopt)) + real(pr) :: v(psize, size(xopt)) + real(pr) :: pbest(psize, size(xopt)), gbest(size(xopt)) + + ! -> PSO and fitness parameters + integer, parameter :: max_iter = 5000, c1 = 2, c2 = 2 + real(pr), parameter :: VTR = 1e-6 + integer :: iter = 1, overfly_counter = 0, NIR = 0, NFE = 0 + logical :: overfly_flag = .false. + real(pr) :: pbest_val(psize) + real(pr) :: gbest_val + real(pr) :: w + real(pr) :: gbest_val_store(max_iter), x11_store(max_iter) + ! <- PSO and fitness parameters + + pbest_val = 0 + gbest_val = 0 + gbest_val_store = 0 + + dimen = size(xopt) + x = 0 + v = 0 + + ! The following loop runs the PSO. Every time this loop starts, a new PSO run is initiated with new seeds. + loopG = 0 + ! Do while (NIR == 0) ! Turn this loop on if you want to run PSO multiple times till it satisfies VTR + + ! until NIR reaches a finite value. If VTR is not achieved, this loop will continue. + + loopG = loopG + 1 + print *, loopG + print *, NIR + + NIR = 0 + NFE = 0 + overfly_counter = 0 + + call cpu_time(time_start) + + ! -> Defining solution space + ! All dimensions are assigned same limit. + ! But this can easily be changed by employing a loop + xmin = -30 + xmax = 30 + + vmin = 0.5*xmin ! lower limit for velocity + vmax = 0.5*xmax ! higher limit for velocity + ! <- Defining solution space + + call random_seed + call random_number(x) + call random_number(v) + + ! -> Intializaing particle position and velocities + + ! The cost of the first particle is evaluated because it is needed + ! as a stepping stool for evaluating gbest. + ! fixing the location of the first particle + x(1, :) = x(1, :)*(xmax - xmin) + xmin + v(1, :) = v(1, :)*(vmax - vmin) + vmin + pbest(1, :) = x(1, :) + gbest = pbest(1, :) + + ! evaluating the fitenss of the first particle + call fitness(pbest(1, :), pbest_val(1)) + NFE = NFE + 1 + gbest_val = pbest_val(1) + + ! Composite threaded loop (15th Aug., 2012) + ! As the population is adjusted, the corresponding fitness parameters of the population is assigned simultaneously. + + do loop1 = 2, psize + x(loop1, :) = x(loop1, :)*(xmax - xmin) + xmin + v(loop1, :) = v(loop1, :)*(vmax - vmin) + vmin + + pbest(loop1, :) = x(loop1, :) + + call fitness(x(loop1, :), pbest_val(loop1)) + NFE = NFE + 1 + if (pbest_val(loop1) < gbest_val) then ! Geared for minimization problem + gbest_val = pbest_val(loop1) + gbest = x(loop1, :) + end if + + end do + + ! <- Intializaing particle position and velocities + + ! Main iteration loop. Creating the generations. + do iter = 1, max_iter + + w = 0.9 - (0.9 - 0.2)*iter/max_iter + + gbest_val_store(iter) = gbest_val + x11_store(iter) = x(1, 1) + + do loop1 = 1, psize + overfly_flag = .false. + do loop2 = 1, dimen + call random_number(temp1) + call random_number(temp2) + v(loop1, loop2) = w*v(loop1, loop2) + c1*temp1*(pbest(loop1, loop2) - x(loop1, loop2)) + v(loop1, loop2) = v(loop1, loop2) + c2*temp2*(gbest(loop2) - x(loop1, loop2)) + + if (v(loop1, loop2) > vmax(loop2)) then + v(loop1, loop2) = vmax(loop2) + end if + + if (v(loop1, loop2) < vmin(loop2)) then + v(loop1, loop2) = vmin(loop2) + end if + + x(loop1, loop2) = x(loop1, loop2) + v(loop1, loop2) + + if (x(loop1, loop2) < xmin(loop2) .or. x(loop1, loop2) > xmax(loop2)) then + overfly_flag = .true. + overfly_counter = overfly_counter + 1 + end if + + end do + + ! Invisible boundary condition. Can be changed to reflecting or absorbing boundary conditions. + if (.not. overfly_flag) then + + !print *, "No overfly" + call fitness(x(loop1, :), temp1) + if (NIR == 0) then + NFE = NFE + 1 + end if + + ! Composite IF statement that adjusts the pbset and gbest parameters (15th Aug., 2012) + if (temp1 < pbest_val(loop1)) then + pbest_val(loop1) = temp1 + pbest(loop1, :) = x(loop1, :) + + if (pbest_val(loop1) < gbest_val) then + gbest_val = pbest_val(loop1) + gbest = pbest(loop1, :) + if (NIR == 0 .and. gbest_val < VTR) then + NIR = iter + end if + end if + + end if + + end if + + end do ! Ending loop over population + + end do ! Ending main iteration loop + + call cpu_time(time_end) + + ! call timestamp(str) + + if (verbose) then + print *, "======================================================" + print *, "Time stamp = ", str + print *, "======================================================" + + print *, "======================================================" + print *, "Algorithm parameters" + print *, "--------------------------------------------------------------------------------------------" + print *, "Problem dimension = ", dimen + print *, "Dynamic range = ", xmin(1), "to", xmax(1) + print *, "Population size = ", psize + print *, "Maximum allowed iterations = ", max_iter + print *, "Defined value to reach (VTR) = ", VTR + print *, "======================================================" + + print *, "" + + print *, "======================================================" + print *, "Performance parameters" + print *, "--------------------------------------------------------------------------------------------" + print *, "Total CPU time = ", time_end - time_start, "seconds" + print *, "Overfly counter = ", overfly_counter + print *, "No. of iterations required to reach VTR (NIR) = ", NIR + print *, "No. of function evaluations (NFE) = ", NFE + print *, "Global best value achieved after NIR = ", gbest_val_store(NIR) + print *, "Global best value achieved after max_iter = ", gbest_val + print *, "" + print *, "Global best solution =" + print *, "" + do loop1 = 1, dimen + print *, gbest(loop1) + end do + print *, "" + !call fitness(dimen, gbest, temp1) + !print *, temp1 + print *, "======================================================" + end if + + !end do ! Ending NIR loop if active. + + ! -> CSV File operation 17th Aug., 2012. (M -> 18th Aug., 2012) + ! call make_csv + ! <- File operation 17th Aug., 2012. + + ! -> Log file operation + ! call make_log + ! <- Log file operation + contains + subroutine make_csv + integer :: file_unit + command_Make = ('mkdir '//'Data_')//str + call system(command_Make) + + command_File = 'Data_'//str//'/'//'converg.csv' + open (newunit=file_unit, file=command_File) + do loop1 = 1, max_iter + write (file_unit, *) gbest_val_store(loop1), ",", x11_store(loop1) + end do + close (file_unit) + + command_File = 'Data_'//str//'/'//'final_x.txt' + open (newunit=file_unit, file=command_File) + do loop1 = 1, psize + write (file_unit, *) x(loop1, :), "," + end do + close (file_unit) + end subroutine + + subroutine make_log + integer :: file_unit + command_File = 'Data_'//str//'/'//'logfile.txt' + open (newunit=file_unit, file=command_File) + + write (file_unit, *) "======================================================" + write (file_unit, *) "Time stamp = ", str + write (file_unit, *) "======================================================" + + write (file_unit, *) "======================================================" + write (file_unit, *) "Algorithm parameters" + write (file_unit, *) "--------------------------------------------------------------------------------------------" + write (file_unit, *) "Problem dimension = ", dimen + write (file_unit, *) "Dynamic range = ", xmin(1), "to", xmax(1) + write (file_unit, *) "Population size = ", psize + write (file_unit, *) "Maximum allowed iterations = ", max_iter + write (file_unit, *) "Defined value to reach (VTR) = ", VTR + write (file_unit, *) "======================================================" + + write (file_unit, *) + + write (file_unit, *) "======================================================" + write (file_unit, *) "Performance parameters" + write (file_unit, *) "---------------------------------------------------------------------------------------------" + write (file_unit, *) "Total CPU time = ", time_end - time_start, "seconds" + write (file_unit, *) "Overfly counter = ", overfly_counter + write (file_unit, *) "No. of iterations required to reach VTR (NIR) = ", NIR + write (file_unit, *) "No. of function evaluations (NFE) = ", NFE + write (file_unit, *) "Global best value achieved after NIR = ", gbest_val_store(NIR) + write (file_unit, *) "Global best value achieved after max_iter = ", gbest_val + write (file_unit, *) + write (file_unit, *) "Global best solution =" + write (file_unit, *) + do loop1 = 1, dimen + write (file_unit, *) gbest(loop1) + end do + write (file_unit, *) + write (file_unit, *) "======================================================" + + close (file_unit) + end subroutine make_log + end subroutine +end module PSOclassicG From 0b81a3656f32d9bc7d75f599c4aa9cd33093f955 Mon Sep 17 00:00:00 2001 From: "Federico E. Benelli" Date: Tue, 5 Mar 2024 13:46:29 -0300 Subject: [PATCH 2/4] example --- example/fitness.f90 | 31 +++++++++++++++++++++++++++++++ example/main.f90 | 11 +++++++++++ 2 files changed, 42 insertions(+) create mode 100644 example/fitness.f90 create mode 100644 example/main.f90 diff --git a/example/fitness.f90 b/example/fitness.f90 new file mode 100644 index 0000000..ba6ae2e --- /dev/null +++ b/example/fitness.f90 @@ -0,0 +1,31 @@ +module fitness_example + use iso_fortran_env, only: pr => real64 + implicit none + +contains + subroutine fitness(x, cost) + real(pr), intent(in) :: x(:) + real(pr), intent(out) :: cost + real(pr) :: temp1 + real(pr) :: temp2 + + integer :: loop1, D + + D = size(x) + + cost = 0 + temp1 = 0 + temp2 = 0 + + ! -> Ackley + do loop1 = 1, D + temp1 = temp1 + x(loop1)**2 + temp2 = temp2 + cos(2*3.141592653589793*x(loop1)) + end do + + temp1 = (temp1/D)**0.5 + temp2 = temp2/D + cost = -20*exp(-0.2*temp1) - exp(temp2) + 20.0 + exp(1.0) + ! <- Ackley + end subroutine fitness +end module diff --git a/example/main.f90 b/example/main.f90 new file mode 100644 index 0000000..c222973 --- /dev/null +++ b/example/main.f90 @@ -0,0 +1,11 @@ +program main + use fitness_example, only: fitness + use PSOclassicG, only: pso + use iso_fortran_env, only: pr => real64 + + real(pr) :: x(30), xmin(30), xmax(30) + + xmax = 30 + xmin = -30 + call pso(fitness, x, psize=30, verbose=.true., xmin=xmin, xmax=xmax) +end program main From be54fcdf5d8eb5a05c328b7073fdc407aaccb66c Mon Sep 17 00:00:00 2001 From: "Federico E. Benelli" Date: Tue, 5 Mar 2024 13:47:41 -0300 Subject: [PATCH 3/4] rework --- src/PSOclassicG.f90 | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/PSOclassicG.f90 b/src/PSOclassicG.f90 index f5637ab..d7fc0bf 100644 --- a/src/PSOclassicG.f90 +++ b/src/PSOclassicG.f90 @@ -40,7 +40,7 @@ module PSOclassicG contains - subroutine pso(fitness, xopt, verbose, psize) + subroutine pso(fitness, xopt, verbose, psize, xmax, xmin) interface subroutine fitness(x, cost) import pr @@ -49,12 +49,13 @@ subroutine fitness(x, cost) end subroutine end interface - real(pr), intent(in out) :: xopt(:) - integer, intent(in) :: psize - logical, intent(in) :: verbose + real(pr), intent(in out) :: xopt(:) !! Optimal values obtained + real(pr), intent(in) :: xmax(size(xopt)) !! Maxumim dimensions + real(pr), intent(in) :: xmin(size(xopt)) !! Minimum dimensions + integer, intent(in) :: psize !! Number of particles + logical, intent(in) :: verbose !! Print output integer :: dimen - real(pr) :: xmin(size(xopt)), xmax(size(xopt)) real(pr) :: vmin(size(xopt)), vmax(size(xopt)) real(pr) :: x(psize, size(xopt)) real(pr) :: v(psize, size(xopt)) @@ -98,9 +99,6 @@ subroutine fitness(x, cost) ! -> Defining solution space ! All dimensions are assigned same limit. ! But this can easily be changed by employing a loop - xmin = -30 - xmax = 30 - vmin = 0.5*xmin ! lower limit for velocity vmax = 0.5*xmax ! higher limit for velocity ! <- Defining solution space @@ -209,6 +207,8 @@ subroutine fitness(x, cost) call cpu_time(time_end) + xopt = gbest + ! call timestamp(str) if (verbose) then From 4f85d393c3126fa1fd1c6adb0cb0908533302e1c Mon Sep 17 00:00:00 2001 From: "Federico E. Benelli" Date: Tue, 5 Mar 2024 13:47:53 -0300 Subject: [PATCH 4/4] toml manifest --- fpm.toml | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) create mode 100644 fpm.toml diff --git a/fpm.toml b/fpm.toml new file mode 100644 index 0000000..43f2987 --- /dev/null +++ b/fpm.toml @@ -0,0 +1,18 @@ +name = "Particle-Swarm-Optimization-Fortran-95" +version = "0.1.0" +license = "license" +author = "Mohammad Asif Zaman" +maintainer = "http://web.stanford.edu/~zaman/" +copyright = "Copyright 2024, Mohammad Asif Zaman " +[build] +auto-executables = true +auto-tests = true +auto-examples = true +module-naming = false +[install] +library = false +[fortran] +implicit-typing = false +implicit-external = false +source-form = "free" +