! ! ---------------------------------- ! Sample program of PSO optimization ! ---------------------------------- ! ! Advanced example ! ! 1. Separate library independent, algorithm independent routines ! library independent: sharable with other optimization libraries ! algorithm independent: sharable with SCE-UA module ! 2. Use low-level primitive routine ! 3. Use user-defined callback function for ! 1) history saving ! 2) convergence checking as exit criterion ! !---------------------------------------------------------------------- ! Module for Common Data (library dependent, algorithm independent) !---------------------------------------------------------------------- module common_data integer, parameter, public :: itmax = 100000 integer, public :: it0(itmax) ! history of iteration count integer, public :: iev0(itmax) ! history of evaluation count real, public :: ev0(itmax) ! history of evaluation value integer, public :: nt = 0 ! number of history data end module common_data !---------------------------------------------------------------------- ! Main Program (library independent, algorithm independent) !---------------------------------------------------------------------- program test3p use common_data integer, parameter :: nmax = 100 integer :: n, i real :: xinit(1:nmax), xmin(1:nmax), xmax(1:nmax) real :: x(1:nmax), r !====== Initialize parameters ====== n = 2 ! number of parameter xmin(1) = 0.0 ! lower bound of x(1) xmax(1) = 12.0 ! upper bound of x(1) xmin(2) = 0.0 ! lower bound of x(2) xmax(2) = 12.0 ! upper bound of x(2) xinit(1) = 2.43 ! initial value of x(1) xinit(2) = 7.57 ! initial value of x(2) !====== Optimize ====== call optimize(n, xmin, xmax, xinit, x, r) !====== Print result ====== print * do i=1, n print *, 'P', i, x(i) enddo print *, 'E', r print * !====== Print history ====== do i=1, nt print *, '---', it0(i), ev0(i), iev0(i) end do print * end program test3p !---------------------------------------------------------------------- ! Objective Function (library independent, algorithm independent) ! ! If you want to share some data with other routines, ! make a module which contains public variables like "common_data", ! referable from other routines. ! ! However, when you enable OpenMP, this function is parallelized. ! Never update shared data without locking from this function. ! ! When OpenMP is enabled, you can write to only local (automatic) ! variables allocated in this function like "t1" and "t2". ! !---------------------------------------------------------------------- function f(x, y) result(z) real, intent(in) :: x, y real :: z real :: t1, t2 t1 = x**4 - 24. * x**3 + 193. * x**2 - 570. * x + 400. t2 = y**4 - 21. * y**3 + 151. * y**2 - 411. * y + 280. z = -3. / 2. * (t1 + t2) end function f !---------------------------------------------------------------------- ! Evaluation Function (library dependent, algorithm independent) !---------------------------------------------------------------------- function evtest(x, ev) result(r) real, intent(in) :: x(:) ! parameter list real, intent(out) :: ev ! evaluation value integer :: r ! number of constrained condition real :: f ev = f(x(1), x(2)) r = 0 end function evtest !---------------------------------------------------------------------- ! Callback Function (library dependent, algorithm independent) ! ! This function is not parallelized during OpenMP enabled, ! you can update shared data safely. ! ! If you don't want to use user-defined callback function, ! please use pre-defined function "pso_cbdummy" in module "opti_pso". ! !---------------------------------------------------------------------- function mycallback(it, iev, ev, x) result(r) use common_data integer, intent(in) :: it ! current iteration counter integer, intent(in) :: iev ! current evaluation counter real, intent(in) :: ev ! current evaluation value real, intent(in) :: x(:) ! current parameter value integer :: r ! continuation flag (>=0:cont., <0:exit) real, parameter :: eps = 1.e-5 !====== Update history ====== nt = nt + 1 it0(nt) = it iev0(nt) = iev ev0(nt) = ev !====== Check Convergence ====== if(nt > 2 .and. abs(ev - ev0(max(nt - 2, 1))) < eps) then r = -1 else r = 1 end if !====== Print intermediate status ====== print *, it, ev, '(', x(1), ',', x(2), ')' end function mycallback !---------------------------------------------------------------------- ! Optimization Routine (library dependent, algorithm dependent) !---------------------------------------------------------------------- subroutine optimize(n, xmin, xmax, xinit, x, r) use opti_pso integer, intent(in) :: n real, intent(in) :: xmin(1:n), xmax(1:n), xinit(1:n) real, intent(out) :: x(1:n), r real :: w, c1, c2, c3, vmax, vrand integer :: np, it, ie, ir, ic, flag external :: evtest, mycallback ! interface declaration is better !====== Setup PSO parameters ====== np = 7 ! population size (number of particle) w = 0.9 ! inertia parameter c1 = 2.0 ! self intention parameter c2 = 2.0 ! swarm intention parameter c3 = 1.e-7 ! random search parameter (f_lrand) vmax = 0.5 ! limit of vector (f_vlimit) vrand = 0.15 ! velocity parametebtion parameter (f_vrand) it = 1000 ! maximum iteration count (<=0 means unlimitted) ie = 0 ! maximum evaluation coun (<=0 means unlimitted) ir = 0 ! interval of message printing (<=0 means no message) ic = 50 ! interval of callback (<=0 means no callback) flag = 0 ! control flag flag = flag + f_quiet ! restrain messages flag = flag + f_initp ! use initial parameters flag = flag + f_lrand ! local random search algorithm flag = flag + f_vlimit ! velocity limit algorithm flag = flag + f_vrand ! velocity perturbation algorithm flag = flag + f_gcenter ! evaluate center of gravity ! flag = flag + f_gjump ! move to center of gravity (unrecommended) ! flag = flag + f_qmcinit ! initialize by Quasi-Monte Carlo call pso(n, xmin, xmax, xinit, np, w, c1, c2, c3, vmax, vrand, & & it, ie, ir, ic, flag, pso_maximize, evtest, mycallback, x, r) end subroutine optimize