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. is the value is off, sa adjusts vm to the correct value).
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.
  • 99 - should not be seen; only used internally.

Calls

proc~~sa~~CallsGraph proc~sa simulated_annealing_module::simulated_annealing_type%sa proc~exprep simulated_annealing_module::exprep proc~sa->proc~exprep proc~func simulated_annealing_module::simulated_annealing_type%func proc~sa->proc~func proc~perturb_and_evaluate simulated_annealing_module::simulated_annealing_type%perturb_and_evaluate proc~sa->proc~perturb_and_evaluate proc~print_vector simulated_annealing_module::print_vector proc~sa->proc~print_vector proc~rand_init simulated_annealing_module::rand_init proc~sa->proc~rand_init proc~uniform_random_number simulated_annealing_module::uniform_random_number proc~sa->proc~uniform_random_number proc~uniform simulated_annealing_module::uniform proc~perturb_and_evaluate->proc~uniform proc~uniform->proc~uniform_random_number

Source Code

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

    implicit none

    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. is the value is off,
                                                         !! sa adjusts vm to the correct value).
    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.
                                                         !!  * 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)                    :: 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`
    real(wp),dimension(:),allocatable    :: vm_original   !! original input value of `vm` (size `n`)
    logical                     :: first !! indicates first function eval
    logical                     :: abort !! indicates the known solution has been found

    ! 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)
    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
        nacp  = 0

        ! the first function eval for a new main cycle
        first = k > 1

        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
                        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
                fopt = me%func(fopt)
                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
                end if
            end if

            !  if termination criteria is not met, prepare for another loop.
            t = abs(rt)*t
            do i = me%neps, 2, -1
                fstar(i) = fstar(i-1)
            end do
            f = fopt
            x = xopt

        end do main

    end do reset

    end subroutine sa