sa Subroutine

private subroutine sa(me, x, rt, t, vm, xopt, fopt, nacc, nfcnev, ier)

Continuous simulated annealing global optimization algorithm

Brief Overview

sa tries to find the global optimum of an n dimensional function. it moves both up and downhill and as the optimization process proceeds, it focuses on the most promising area.

to start, it randomly chooses a trial point within the step length vm (a vector of length n) of the user selected starting point. the function is evaluated at this trial point and its value is compared to its value at the initial point.

in a maximization problem, all uphill moves are accepted and the algorithm continues from that trial point. downhill moves may be accepted; the decision is made by the metropolis criteria. it uses t (temperature) and the size of the downhill move in a probabilistic manner. the smaller t and the size of the downhill move are, the more likely that move will be accepted. if the trial is accepted, the algorithm moves on from that point. if it is rejected, another point is chosen instead for a trial evaluation.

each element of vm periodically adjusted so that half of all function evaluations in that direction are accepted.

a fall in t is imposed upon the system with the rt variable by t(i+1) = rt*t(i) where i is the ith iteration. thus, as t declines, downhill moves are less likely to be accepted and the percentage of rejections rise. given the scheme for the selection for vm, vm falls. thus, as t declines, vm falls and sa focuses upon the most promising area for optimization.

the parameter t is crucial in using sa successfully. it influences vm, the step length over which the algorithm searches for optima. for a small initial t, the step length may be too small; thus not enough of the function might be evaluated to find the global optima. the user should carefully examine vm in the intermediate output (set iprint = 1) to make sure that vm is appropriate. the relationship between the initial temperature and the resulting step length is function dependent.

to determine the starting temperature that is consistent with optimizing a function, it is worthwhile to run a trial run first. set rt = 1.5 and t = 1.0. with rt > 1.0, the temperature increases and vm rises as well. then select the t that produces a large enough vm.

Notes

  • as far as possible, the parameters here have the same name as in the description of the algorithm on pp. 266-8 of corana et al.
  • the suggested input values generally come from corana et al. to drastically reduce runtime, see goffe et al., pp. 90-1 for suggestions on choosing the appropriate rt and nt.

History

  • differences compared to version 2.0:
    1. if a trial is out of bounds, a point is randomly selected from lb(i) to ub(i). unlike in version 2.0, this trial is evaluated and is counted in acceptances and rejections. all corresponding documentation was changed as well.
  • differences compared to version 3.0:
    1. if vm(i) > (ub(i) - lb(i)), vm is set to ub(i) - lb(i). the idea is that if t is high relative to lb & ub, most points will be accepted, causing vm to rise. but, in this situation, vm has little meaning; particularly if vm is larger than the acceptable region. setting vm to this size still allows all parts of the allowable region to be selected.
  • differences compared to version 3.1:
    1. test made to see if the initial temperature is positive.
    2. write statements prettied up.
    3. references to paper updated.

Type Bound

simulated_annealing_type

Arguments

Type IntentOptional Attributes Name
class(simulated_annealing_type), intent(inout) :: me
real(kind=wp), intent(inout), dimension(me%n) :: x

on input: the starting values for the variables of the function to be optimized. [Will be replaced by final point]
real(kind=wp), intent(in) :: rt

the temperature reduction factor. the value suggested by corana et al. is .85. see goffe et al. for more advice.
real(kind=wp), intent(inout) :: t

on input, the initial temperature. see goffe et al. for advice. on output, the final temperature. Note that if t=0, then all downhill steps will be rejected
real(kind=wp), intent(inout), dimension(me%n) :: vm

the step length vector. on input it should encompass the region of interest given the starting value x. for point x(i), the next trial point is selected is from x(i) - vm(i) to x(i) + vm(i). since vm is adjusted so that about half of all points are accepted, the input value is not very important (i.e. if the value is off, sa adjusts vm to the correct value). note: if vm=0.0, then it is set to abs(xu-xl)`.
real(kind=wp), intent(out), dimension(me%n) :: xopt

the variables that optimize the function.
real(kind=wp), intent(out) :: fopt

the optimal value of the function.
integer, intent(out) :: nacc

the number of accepted function evaluations.
integer, intent(out) :: nfcnev

the total number of function evaluations. in a minor point, note that the first evaluation is not used in the core of the algorithm; it simply initializes the algorithm.
integer, intent(out) :: ier

the error return number:

  • 0 - normal return; termination criteria achieved.
  • 1 - number of function evaluations (nfcnev) is greater than the maximum number (maxevl).
  • 3 - the initial temperature is not positive.
  • 4 - user stop in problem function.
  • 5 - step sizes collapsed below vm_min threshold.
  • 99 - should not be seen; only used internally.

Calls

proc~~sa~~CallsGraph proc~sa simulated_annealing_type%sa proc~exprep exprep proc~sa->proc~exprep proc~func simulated_annealing_type%func proc~sa->proc~func proc~perturb_and_evaluate simulated_annealing_type%perturb_and_evaluate proc~sa->proc~perturb_and_evaluate proc~print_vector print_vector proc~sa->proc~print_vector proc~rand_init rand_init proc~sa->proc~rand_init proc~uniform_random_number uniform_random_number proc~sa->proc~uniform_random_number proc~perturb_variable simulated_annealing_type%perturb_variable proc~perturb_and_evaluate->proc~perturb_variable proc~bipareto bipareto proc~perturb_variable->proc~bipareto proc~cauchy cauchy proc~perturb_variable->proc~cauchy proc~triangular_dist triangular_dist proc~perturb_variable->proc~triangular_dist proc~truncated_normal truncated_normal proc~perturb_variable->proc~truncated_normal proc~uniform uniform proc~perturb_variable->proc~uniform proc~bipareto->proc~uniform_random_number proc~bipareto->proc~uniform proc~cauchy->proc~uniform_random_number proc~triangular_dist->proc~uniform_random_number proc~truncated_normal->proc~uniform proc~normal normal proc~truncated_normal->proc~normal proc~uniform->proc~uniform_random_number proc~normal->proc~uniform_random_number

Called by

proc~~sa~~CalledByGraph proc~sa simulated_annealing_type%sa proc~solve_simulated_annealing solve_simulated_annealing proc~solve_simulated_annealing->proc~sa

Source Code

   subroutine sa(me, x, rt, t, vm, xopt, fopt, nacc, nfcnev, ier)

      class(simulated_annealing_type),intent(inout) :: me
      real(wp), dimension(me%n), intent(inout)  :: x       !! on input: the starting values for the variables of the
                                                           !! function to be optimized. [Will be replaced by final point]
      real(wp), intent(in)                      :: rt      !! the temperature reduction factor.  the value suggested by
                                                           !! corana et al. is .85. see goffe et al. for more advice.
      real(wp), intent(inout)                   :: t       !! on input, the initial temperature. see goffe et al. for advice.
                                                           !! on output, the final temperature. Note that if `t=0`, then
                                                           !! all downhill steps will be rejected
      real(wp), dimension(me%n), intent(inout)  :: vm      !! the step length vector. on input it should encompass the region of
                                                           !! interest given the starting value x.  for point x(i), the next
                                                           !! trial point is selected is from x(i) - vm(i)  to  x(i) + vm(i).
                                                           !! since vm is adjusted so that about half of all points are accepted,
                                                           !! the input value is not very important (i.e. if the value is off,
                                                           !! sa adjusts vm to the correct value).
                                                           !! note: if `vm=0.0`, then it is set to abs(xu-xl)`.
      real(wp), dimension(me%n), intent(out)  :: xopt      !! the variables that optimize the function.
      real(wp), intent(out)                   :: fopt      !! the optimal value of the function.
      integer, intent(out)                    :: nacc      !! the number of accepted function evaluations.
      integer, intent(out)                    :: nfcnev    !! the total number of function evaluations. in a minor
                                                           !! point, note that the first evaluation is not used in the
                                                           !! core of the algorithm; it simply initializes the
                                                           !! algorithm.
      integer, intent(out)                    :: ier       !! the error return number:
                                                           !!
                                                           !!  * 0 - normal return; termination criteria achieved.
                                                           !!  * 1 - number of function evaluations (nfcnev) is
                                                           !!    greater than the maximum number (maxevl).
                                                           !!  * 3 - the initial temperature is not positive.
                                                           !!  * 4 - user stop in problem function.
                                                           !!  * 5 - step sizes collapsed below `vm_min` threshold.
                                                           !!  * 99 - should not be seen; only used internally.

      real(wp),dimension(:),allocatable :: xp           !! perturbed `x` vector (size `n`)
      real(wp),dimension(:),allocatable :: fstar        !! array of optimals found so far (size `neps`)
      real(wp),dimension(:),allocatable :: vm_original  !! original input value of `vm` (size `n`)
      real(wp)  :: f             !! function value
      real(wp)  :: fp            !! function value
      real(wp)  :: p             !! for metropolis criteria
      real(wp)  :: pp            !! for metropolis criteria
      real(wp)  :: ratio         !! ratio of `nacp` to `ns`
      integer   :: nup           !! number of uphill steps
      integer   :: ndown         !! number of accepted downhill steps
      integer   :: nrej          !! number of rejected downhill steps
      integer   :: nnew          !! new optimal for this temperature
      integer   :: i, j, m, k    !! counter
      integer   :: totmov        !! total moves
      integer   :: nacp          !! number of accepts steps in a complete `ns` cycle
      logical   :: quit          !! if the termination criteria was achieved
      integer   :: unit          !! unit for printing
      real(wp)  :: t_original    !! original input value of `t`
      logical   :: first         !! indicates first function eval
      logical   :: abort         !! indicates the known solution has been found
      real(wp)  :: t0            !! initial temperature (for non-geometric cooling schedules)
      integer   :: temp_iter     !! temperature iteration counter (for non-geometric cooling schedules)

      ! initialize:
      allocate(xp(me%n))
      allocate(fstar(me%neps))
      allocate(vm_original(me%n))
      unit        = me%iunit
      nfcnev      = 0
      ier         = 99
      xopt        = x
      fopt        = huge(1.0_wp)
      vm          = abs(vm)
      where (vm == 0.0_wp)
         vm = abs(me%ub - me%lb)
      end where
      t_original  = t
      vm_original = vm
      abort       = .false.

      ! if the initial temperature is not positive, notify the user and
      ! return to the calling routine.
      if (t < 0.0_wp) then
         if (me%iprint>0) write(unit,'(A)') 'Error: The initial temperature must be non-negative.'
         ier = 3
         return
      end if

      ! if the initial value is out of bounds, then just put
      ! the violated ones on the nearest bound.
      where (x > me%ub)
         x = me%ub
      else where (x < me%lb)
         x = me%lb
      end where

      !  initialize the random number generator
      call rand_init(me%iseed1,me%iseed2)

      ! evaluate the function with input x and return value as f.
      ! keep trying if necessary until there is a valid function
      ! evaluation.
      call me%perturb_and_evaluate(x,vm,xp,f,nfcnev,ier,first=.true.)
      select case (ier)
       case(1,4)
         return
      end select

      f    = me%func(f)
      xopt = xp ! it may have been perturbed above
      fopt = f
      if (me%iprint >= 1) then
         write(unit, '(A)') '  '
         call print_vector(unit,xp,me%n,'initial x')
         write(unit, '(A,1X,G25.18)')  '  initial f: ', me%func(f)
      end if

      reset: do k = 1, me%n_resets

         if (me%iprint >= 1) then
            write(unit,'(A)')       ''
            write(unit,'(A)')       '=========================================='
            write(unit,'(A,1X,I5)') ' main loop ', k
            write(unit,'(A)')       '=========================================='
            write(unit,'(A)')       ''
         end if

         ! restart the main loop
         fstar     = huge(1.0_wp)
         t         = t_original
         vm        = vm_original
         nacc      = 0
         temp_iter = 0      ! reset temperature iteration counter
         nacp      = 0
         first     = k > 1  ! the first function eval for a new main cycle

         main : do
            ! start the main loop. note that it terminates if (i) the algorithm
            ! succesfully optimizes the function or (ii) there are too many
            ! function evaluations (more than maxevl).
            nup = 0
            nrej = 0
            nnew = 0
            ndown = 0

            nt_loop : do m = 1, me%nt
               ns_loop : do j = 1, me%ns

                  ! evaluate the function with the trial point xp and return as fp.
                  call me%perturb_and_evaluate(x,vm,xp,fp,nfcnev,ier,first=first)
                  first = .false.
                  select case (ier)
                   case(1,4)
                     x    = xopt
                     fopt = me%func(fopt)
                     return
                  end select

                  fp = me%func(fp)
                  if (me%iprint >= 3) then
                     write(unit,'(A)') ' '
                     call print_vector(unit,x,me%n,'current x')
                     write(unit,'(A,G25.18)') ' current f: ', me%func(f)
                     call print_vector(unit,xp,me%n,'trial x')
                     write(unit,'(A,G25.18)') ' resulting f: ', me%func(fp)
                  end if

                  if (fp >= f) then
                     ! function value increased: accept the new point

                     if (me%iprint >= 3) write(unit,'(A)') '  point accepted'
                     x = xp
                     f = fp
                     nacc = nacc + 1
                     nacp = nacp + 1
                     nup = nup + 1

                     ! if greater than any other point, record as new optimum.
                     if (fp > fopt) then
                        if (me%iprint >= 3) write(unit,'(A)') '  new optimum'
                        xopt = xp
                        fopt = fp
                        nnew = nnew + 1
                        if (me%ireport == 2 .or. me%ireport == 3) then ! report this value to the user
                           call me%report(xopt, me%func(fopt), istat=2) ! convert f back to user's sign if necessary
                        end if
                     end if

                  else
                     ! if the point is lower, use the metropolis criteria to decide on
                     ! acceptance or rejection.
                     ! NOTE: if t=0, all downhill steps are rejected (monotonic)

                     if (t/=0.0_wp) then
                        p = exprep((fp - f)/t)
                        pp = uniform_random_number()
                        if (pp < p) then
                           if (me%iprint >= 3) then
                              if (me%maximize) then
                                 write(unit,'(A)')  '  though lower, point accepted'
                              else
                                 write(unit,'(A)')  '  though higher, point accepted'
                              end if
                           end if
                           x = xp
                           f = fp
                           nacc = nacc + 1
                           nacp = nacp + 1
                           ndown = ndown + 1
                           cycle  ! accepted
                        end if
                     end if

                     ! rejected:
                     nrej = nrej + 1
                     if (me%iprint >= 3) then
                        if (me%maximize) then
                           write(unit,'(A)') '  lower point rejected'
                        else
                           write(unit,'(A)') '  higher point rejected'
                        end if
                     end if

                  end if

               end do ns_loop ! j - ns loop

               select case (me%step_mode)
                case(1)
                  ! adjust vm so that approximately half of all evaluations are accepted.
                  ratio = real(nacp,wp) /real(me%ns,wp)
                  do i = 1, me%n
                     if (ratio > 0.6_wp) then
                        vm(i) = vm(i)*(1.0_wp + me%c(i)*(ratio - 0.6_wp)/0.4_wp)
                     else if (ratio < 0.4_wp) then
                        vm(i) = vm(i)/(1.0_wp + me%c(i)*((0.4_wp - ratio)/0.4_wp))
                     end if
                     if (vm(i) > (me%ub(i)-me%lb(i))) then
                        vm(i) = me%ub(i) - me%lb(i)
                     end if
                  end do
                  if (me%iprint >= 2) then
                     write(unit,'(/A)') '---------------------------------------------------'
                     write(unit,'(A)')  ' intermediate results after step length adjustment '
                     write(unit,'(A/)') '---------------------------------------------------'
                     call print_vector(unit, vm,   me%n, 'new step length (vm)')
                     call print_vector(unit, xopt, me%n, 'current optimal x')
                     call print_vector(unit, x,    me%n, 'current x')
                     write(unit,'(A)') ' '
                  end if
                case (2)
                  ! keep vm as is
                case (3)
                  ! use the constant factor:
                  vm = abs(me%vms) * vm
                case default
                  error stop ' invalid value of step_mode'
               end select

               nacp = 0

            end do nt_loop ! m - nt loop

            if (me%iprint >= 1) then
               totmov = nup + ndown + nrej
               write(unit,'(/A)')      '--------------------------------------------------------'
               write(unit,'(A)')       ' intermediate results before next temperature reduction '
               write(unit,'(A/)')      '--------------------------------------------------------'
               write(unit,'(A,G12.5)') '  current temperature:            ', t
               if (me%maximize) then
                  write(unit,'(A,G25.18)') '  max function value so far:  ', fopt
                  write(unit,'(A,I8)')     '  total moves:                ', totmov
                  write(unit,'(A,I8)')     '     uphill:                  ', nup
                  write(unit,'(A,I8)')     '     accepted downhill:       ', ndown
                  write(unit,'(A,I8)')     '     rejected downhill:       ', nrej
                  write(unit,'(A,I8)')     '  new maxima this temperature:', nnew
               else
                  write(unit,'(A,G25.18)') '  min function value so far:  ', -fopt
                  write(unit,'(A,I8)')     '  total moves:                ', totmov
                  write(unit,'(A,I8)')     '     downhill:                ', nup
                  write(unit,'(A,I8)')     '     accepted uphill:         ', ndown
                  write(unit,'(A,I8)')     '     rejected uphill:         ', nrej
                  write(unit,'(A,I8)')     '  new minima this temperature:', nnew
               end if
               call print_vector(unit,xopt, me%n, 'current optimal x')
               call print_vector(unit,vm,   me%n, 'step length (vm)')
               write(unit, '(A)') ' '
            end if

            ! check termination criteria.
            quit = .false.
            ! first, check `f` if the known value is available
            if (me%optimal_f_specified) then
               quit = abs(me%func(fopt) - me%optimal_f) <= abs(me%optimal_f_tol)
               if (quit) abort = .true.
            end if
            ! check the convergence tolerances:
            if (.not. quit) then
               fstar(1) = f
               quit = ((fopt-f)<=me%eps) .and. (all(abs(f-fstar)<=me%eps))
            end if

            ! terminate sa if appropriate.
            if (quit) then
               x = xopt
               ier = 0
               if (me%iprint >= 1) then
                  write(unit,'(/A)') '----------------------------------------------'
                  write(unit,'(A)')  '  sa achieved termination criteria. ier = 0.  '
                  write(unit,'(A/)') '----------------------------------------------'
               end if
               if (abort) then
                  exit reset ! the solution has been found, so can ignore
                  ! the reset loop if that is being used
               else
                  exit main ! exit the main loop but continue resets
               end if
            end if

            ! check if step sizes have collapsed
            if (maxval(vm) < me%vm_min) then
               if (me%iprint >= 1) then
                  write(unit,'(/A)') '-----------------------------------------------------'
                  write(unit,'(A)')  '  Step sizes collapsed below vm_min threshold'
                  write(unit,'(A/)') '-----------------------------------------------------'
               end if
               x = xopt
               ier = 5  ! error code for "step collapse"
               exit main  ! exit main loop to try next reset if available
            end if

            !  if termination criteria is not met, prepare for another loop.
            temp_iter = temp_iter + 1

            if (t_original /= 0.0_wp) then
               ! apply the selected cooling schedule
               select case (me%cooling_schedule)
               case (1)
                  ! geometric (classical): T(k+1) = rt * T(k)
                  t = abs(rt)*t
               case (2)
                  ! fast annealing (Cauchy): T(k) = T0 / (1 + k)
                  t = t_original / real(1 + temp_iter, wp)
               case (3)
                  ! huang: T(k) = T0 / (1 + c*k)^n
                  t = t_original / (1.0_wp + me%cooling_param * real(temp_iter, wp))**me%n
               case (4)
                  ! boltzmann: T(k) = T0 / log(1 + k + e)
                  t = t_original / log(1.0_wp + real(temp_iter, wp) + exp(1.0_wp))
               case (5)
                  ! logarithmic: T(k) = T0 / (1 + c*log(1+k))
                  t = t_original / (1.0_wp + me%cooling_param * log(1.0_wp + real(temp_iter, wp)))
               case default
                  ! fallback to geometric
                  t = abs(rt)*t
               end select
            end if

            do i = me%neps, 2, -1
               fstar(i) = fstar(i-1)
            end do
            f = fopt
            x = xopt

         end do main

      end do reset

      ! convert fopt back to user space before returning
      fopt = me%func(fopt)

   end subroutine sa