(list of modules)(lisf of routines)
    1   program optimisation
    2   
    3     ! This program is meant to solve an optimization problem: it finds
    4     ! a local/global minimum of a function E (objective function).
    5     ! The optimization is always done under a basic constraint
    6     ! on the range of each variable.
    7     ! Additional constraints can be added, of the form H_i(x) = 0
    8     ! Several methods are available, with tunable parameters.
    9     ! Local optimization methds include Nelder-Mead ('simplex'),
   10     ! steepest descent ('steepes') and BFGS ('bfgs'), while
   11     ! global optimization methods include simulated annealing
   12     ! ('annealing') and a genetic algorithm ('genetic').
   13     ! Note that steepest descent and BFGS require that partial
   14     ! derivates of E are also given.
   15     ! Constrained optimization uses either a penalty method or
   16     ! Lagrange multipliers (and an unconstrained optimization is
   17     ! solved iteratively).
   18   
   19     ! An input file ('optim.in') governs the computation.
   20     ! The result is written to screen
   21     ! Convergence history is written in file 'optim.out'
   22   
   23     ! References:
   24     !    - Numerical Recipes (Nelder-Mead)
   25     !    - Numerical Mathematics, A. Quarteroni, Springer (constrained)
   26     !    - Wikipedia (BFGS, line search)
   27   
   28     ! Author: L. Klinger (2006)
   29   
   30     ! History:
   31   
   32   
   33     use knuth
   34   
   35     implicit none
   36   
   37     integer, parameter :: KD = selected_real_kind(14,200)
   38   
   39     ! this hold the configuration for the optimization
   40     ! task (problem size, algorithm to use, numerical parameters,...)
   41     type oparam
   42        integer           :: N, M, n_cons
   43        character(len=16) :: method ! simplex,genetic,steepest,annealing,BFG
   44        character(len=16) :: cons_method ! penalty, multiplier
   45        integer           :: itmax, GA_popsize
   46        integer           :: itmaxcons, seed
   47        real(KD)          :: tol, ctol, ls_sigma, ls_rho, simplex_lambda
   48        real(KD)          :: prob_mut, prob_cross, SA_beta
   49        real(KD)          :: calpha, cbeta, clambda
   50        logical           :: use_cons
   51        real(KD), dimension(:), pointer :: pmin, pmax
   52     end type oparam
   53   
   54   
   55     ! global variable to hold constraint
   56     logical  :: constrained, use_lagrange_multipliers
   57     real(KD) :: cons_alpha
   58     integer  :: n_cons
   59     real(KD), dimension(:), pointer :: lagrange_multipliers
   60     real(KD), dimension(:), pointer :: p
   61   
   62     real         :: t1, t2
   63     integer      :: iter, fec, out
   64     type(oparam) :: op
   65   
   66   
   67   
   68     ! initialize global variables
   69     constrained = .false.
   70     n_cons      = 1
   71     cons_alpha  = 0.0d0
   72     fec         = 0      ! function evaluation counte
   73     out         = 11     ! unit for convergence histor
   74   
   75     ! initialize computation by (mainly) reading the config file
   76     ! and filling OP
   77     call init( op, p )
   78     call randomize( op%seed )
   79     open( unit = out, file = 'optim.out' )
   80   
   81     ! computational part
   82     call cpu_time( t1 )
   83     if ( op%use_cons ) then
   84        call COptim( op, p, iter )
   85     else
   86        call UOptim( op, p, iter )
   87     end if
   88     call cpu_time( t2 )
   89   
   90     ! display results
   91     print *, 'Best solution found: at'
   92     do n_cons = 1, size(p)
   93        write( *, '(es20.10)' ) p(n_cons)
   94     end do
   95     print *, 'objective function is ', E(p)
   96   
   97     print *
   98     print *, 'Number of function evaluations: ', fec
   99     print *, 'Number of iterations:           ', iter
  100     print *, 'Elapsed CPU time:               ', t2 - t1
  101   
  102     close( out )
  103   
  104     ! clean-up
  105     if ( associated( op%pmin ) ) deallocate( op%pmin )
  106     if ( associated( op%pmax ) ) deallocate( op%pmax )
  107     if ( associated( p ) )       deallocate( p )
  108   
  109   
  110   contains
  111   
  112   !==========================================================================
  113   function E( x )
  114   !==========================================================================
  115   ! The objective function.
  116   
  117     real(KD), dimension(:), intent(in) :: x
  118     real(KD) :: E
  119   
  120     E = 1.0d2 * (x(2)-x(1)**2)**2 + (1.0d0 - x(1))**2
  121   !  E = 0.5d0 * x(1) + 0.25d0 * x(2) + 1.0d0
  122   
  123     if ( constrained ) then
  124        E = E + 0.5d0 * cons_alpha * NormL2( H( x, n_cons ) )
  125        if ( use_lagrange_multipliers ) then
  126           E = E + dot_product( lagrange_multipliers, H(x,n_cons) )
  127        end if
  128     end if
  129   
  130     fec = fec + 1
  131   
  132   end function E
  133   
  134   
  135   !==========================================================================
  136   function GradE( x )
  137   !==========================================================================
  138   ! The gradient of the objective function (used for steepest descent and
  139   ! BFSG algorithms)
  140   
  141     real(KD), dimension(:), intent(in) :: x
  142     real(KD), dimension(size(x))       :: GradE
  143   
  144     integer :: i
  145   
  146     do i = 1, size(x)
  147        GradE(1) = -( 4.0d2*x(1)*(x(2)-x(1)**2) + 2.0d0 *(1.0d0-x(1)))
  148        GradE(2) = 2.0d2 * (x(2) - x(1)**2)
  149   !     GradE(1) = 0.5d0
  150   !     GradE(2) = 0.25d0
  151     end do
  152   
  153   end function GradE
  154   
  155   
  156   !==========================================================================
  157   function H( x, m )
  158   !==========================================================================
  159   ! constaints in the form H(x) == 0
  160   
  161     real(KD), dimension(:), intent(in) :: x
  162     integer,                intent(in) :: m ! the number of constaint
  163     real(KD), dimension(m) :: H
  164   
  165     H = 0.0d0
  166   
  167     H(1) = (x(1)+2.0d0)**2 + (x(2)+2.0d0)**2 - 0.25d0
  168     !H(2) = ...
  169     !H(3) = ...
  170   
  171   end function H
  172   
  173   
  174   
  175   !**************************************************************************
  176   !**************************************************************************
  177   !**
  178   !**  Constrained optim (penalty method or Lagrange multipliers one)
  179   !**
  180   !**************************************************************************
  181   !**************************************************************************
  182   
  183   
  184   !==========================================================================
  185   subroutine COptim( params, sol, iter )
  186   !==========================================================================
  187   ! constrained optimization driver. Uses penalty method or Lagrange
  188   ! multipliers. SOL, on input is an initial guess (unconstrained), and
  189   ! on output, is the results. ITER returns then number of outer iterations 
  190   ! used.
  191   
  192     type(oparam),           intent(in)    :: params
  193     real(KD), dimension(:), intent(inout) :: sol
  194     integer,                intent(out)   :: iter
  195   
  196     real(KD) :: alpha, beta, worst_cons, tol, lambda0
  197     integer  :: uiter, meth
  198   
  199     real(KD), dimension(size(sol)) :: x
  200   
  201     meth = 1
  202     if ( trim(params%cons_method) == 'multipliers' ) meth = 2
  203   
  204     alpha   = op%calpha
  205     beta    = op%cbeta
  206     n_cons  = op%n_cons
  207     iter    = 0
  208     tol     = op%ctol
  209     lambda0 = op%clambda
  210   
  211     constrained = .true.
  212     cons_alpha  = alpha
  213     use_lagrange_multipliers = .false.
  214   
  215     if ( meth == 2 ) then
  216        use_lagrange_multipliers = .true.
  217        if ( associated( lagrange_multipliers ) ) deallocate( lagrange_multipliers )
  218        allocate( lagrange_multipliers(n_cons) )
  219        lagrange_multipliers = lambda0
  220     end if     
  221   
  222     x = sol
  223   
  224     ! worst_cons tells how well the current X satisfies the constraint
  225     worst_cons = maxval( H( x, n_cons ) )
  226   
  227     do while ( worst_cons > tol )
  228   
  229        call UOptim( params, x, uiter )
  230        worst_cons = maxval( H( x, n_cons ) )
  231   
  232        write( *, '(i3,6es18.10)' ) iter, worst_cons, x, cons_alpha
  233        write( out, '(6es20.10)' ) x(1:min(op%n,6))
  234   
  235        iter = iter + 1
  236        if ( iter > params%itmaxcons ) exit
  237   
  238        cons_alpha = cons_alpha * beta
  239   
  240        if ( meth == 2 ) then
  241           lagrange_multipliers = lagrange_multipliers + cons_alpha * H( x, n_cons )
  242        end if
  243   
  244     end do
  245   
  246     sol = x
  247     constrained = .false.
  248   
  249     if ( meth == 2 ) deallocate( lagrange_multipliers )
  250   
  251   end subroutine COptim
  252   
  253   
  254   
  255   
  256   !**************************************************************************
  257   !**************************************************************************
  258   !**
  259   !**  Unconstrained optimization (wrapper)
  260   !**
  261   !**************************************************************************
  262   !**************************************************************************
  263   
  264   
  265   !==========================================================================
  266   subroutine UOptim( params, sol, iter )
  267   !==========================================================================
  268   ! solves an Unconstrained optimization problem (minimizing E), according
  269   ! to PARAMS, to get SOL and ITER (number of iterations used). On input,
  270   ! SOL is an initial guess (not relevant for genetic and annealing)
  271   
  272     type(oparam),           intent(in)    :: params
  273     real(KD), dimension(:), intent(inout) :: sol
  274     integer,                intent(out)   :: iter
  275   
  276     select case ( trim(params%method) )
  277     case( 'simplex' )
  278        call Simplex_optim( params, sol, iter )
  279     case( 'genetic' )
  280        call GA_optim( params, sol, iter )
  281     case( 'steepest' )
  282        call SD_optim( params, sol, iter )
  283     case( 'annealing' )
  284        call SA_optim( params, sol, iter )
  285     case( 'BFGS' )
  286        call BFGS_optim( params, sol, iter )
  287     case default
  288        print *, 'optimization method "', trim(params%method), '" unknown.'
  289        sol  = 1.0D+33
  290        iter = 0
  291     end select
  292   
  293   end subroutine UOptim
  294   
  295   
  296   
  297   !**************************************************************************
  298   !**************************************************************************
  299   !**
  300   !**  BFGS
  301   !**
  302   !**************************************************************************
  303   !**************************************************************************
  304   
  305   
  306   !==========================================================================
  307   subroutine BFGS_optim( params, optim, iter )
  308   !==========================================================================
  309   ! finds a local minimum with BFGS method
  310   
  311     type(oparam),           intent(in)    :: params
  312     real(KD), dimension(:), intent(inout) :: optim
  313     integer,                intent(out)   :: iter
  314   
  315     real(KD), dimension(size(optim)) :: p0, y, p, grad, s
  316     real(KD), dimension(size(optim),size(optim)) :: B, A
  317   
  318     real(KD) :: alpha, tra1, tra2, f, fold
  319   
  320     integer :: i, j
  321   
  322     ! initalize B
  323     B = 0.0d0
  324     do i = 1, size(optim)
  325        B(i,i) = 1.0d0
  326     end do
  327   
  328     p    = optim
  329     fold = E(p)
  330     iter = 0
  331   
  332     do 
  333        grad = GradE( p )
  334        A = B ! make a copy of B since it will be overwritten by solv
  335        call solve( A, s, -grad )
  336        call line_search( p, s, alpha, params )
  337        p = p + alpha * s
  338   
  339        if ( params%n_cons == 0 ) write( out, '(6es20.10)' ) p(1:min(params%n,6))
  340   
  341        y = GradE( p ) - grad
  342   
  343        if ( modulus(y) > 1.0D-9 ) then
  344           tra1 = 1.0d0 / dot_product( y, s )
  345           tra2 = 1.0d0 / dot_product( s, matmul( B, s ) )
  346        
  347           do j = 1, size(optim)
  348              A(:,j) = s * s(j)
  349           end do
  350           A = matmul( A, B )
  351           A = matmul( B, A )
  352           do i = 1, size(optim)
  353              do j = 1, size(optim)
  354                 B(i,j) = B(i,j) + y(i)*y(j) * tra1 - &
  355                      A(i,j) * tra2
  356              end do
  357           end do
  358        end if
  359        f = E(p)
  360        iter = iter + 1
  361        if ( iter > params%itmax .or. abs(f-fold) < params%tol ) exit
  362        fold = f
  363     end do
  364     
  365     optim = p
  366   
  367   end subroutine BFGS_optim
  368   
  369   
  370   !==========================================================================
  371   subroutine solve( A, x, b )
  372   !==========================================================================
  373   ! wrapper for lapack solve
  374   
  375     real(KD), dimension(:,:), intent(in)    :: A
  376     real(KD), dimension(:),   intent(in)    :: b
  377     real(KD), dimension(:),   intent(inout) :: x
  378   
  379     integer :: info
  380     integer, dimension(size(b)) :: ipiv
  381   
  382     x = b
  383     call DGESV( size(A,1), 1, A, size(A,1), ipiv, x, size(b), info ) 
  384     if ( info /= 0 ) then
  385        x = 1.0D+33
  386        print *, 'info = ', info
  387     end if
  388   
  389   end subroutine solve
  390   
  391   
  392   !==========================================================================
  393   subroutine SD_optim( params, optim, iter )
  394   !==========================================================================
  395   ! uses steepest descent to perform optimizsation of function E
  396   
  397     type(oparam),           intent(in)    :: params
  398     real(KD), dimension(:), intent(inout) :: optim
  399     integer,                intent(out)   :: iter
  400   
  401     real(KD), dimension(size(optim)) :: p0, y, p, grad, s
  402   
  403     real(KD) :: alpha, tra1, tra2
  404   
  405     integer :: i, j
  406   
  407     p = optim
  408     tra1 = E(p)
  409     iter = 0
  410   
  411     do 
  412        grad = GradE( p )
  413        s = -grad;
  414        call line_search( p, s, alpha, params )
  415        p = p + alpha * s
  416   
  417        if ( params%n_cons == 0 ) write( out, '(6es20.10)' ) p(1:min(params%n,6))
  418   
  419        tra2 = E(p)
  420        iter = iter + 1
  421        if ( abs(tra2-tra1) < params%tol .or. iter > params%itmax ) exit
  422        tra1 = tra2
  423     end do
  424     
  425     optim = p
  426   
  427   end subroutine SD_optim
  428   
  429   
  430   !==========================================================================
  431   subroutine line_search( x, s, alpha, params )
  432   !==========================================================================
  433   
  434     real(KD), dimension(:), intent(in)    :: x
  435     real(KD), dimension(:), intent(inout) :: s
  436     real(KD),               intent(out)   :: alpha
  437     type(oparam),           intent(in)    :: params
  438   
  439     integer  :: j
  440     real(KD) :: sigma, rho
  441   
  442     ! make sure the search won't end outside the allowed domain
  443     call largest_alpha( x, s, alpha, params%pmin, params%pmax )
  444   
  445     if ( alpha == 0.0d0 ) then
  446        ! project s onto the constraint
  447        do j = 1, size(x)
  448           if ( x(j) == params%pmin(j) .or. x(j) == params%pmax(j) ) then
  449              s(j) = 0.0d0
  450              exit
  451           end if
  452        end do
  453        call largest_alpha( x, s, alpha, params%pmin, params%pmax )
  454     end if
  455   
  456     if ( alpha < 0.0d0 ) alpha = -alpha
  457   
  458     sigma = params%ls_sigma !1.0D-
  459     rho   = params%ls_rho   !0.5d
  460     alpha = min( 1.0, alpha )
  461   
  462     do j = 1, 20
  463        if ( E(x) - E(x + alpha * s) >= -sigma * alpha * &
  464             dot_product( GradE( x ), s ) ) return
  465        alpha = alpha * rho
  466     end do
  467   
  468   end subroutine line_search
  469   
  470   
  471   !==========================================================================
  472   subroutine largest_alpha( x, s, alpha, pmin, pmax )
  473   !==========================================================================
  474   ! finds the largest value of ALPHA, so that X + ALPHA*S remains in the
  475   ! limits of the range of all variables
  476   
  477     real(KD), dimension(:), intent(in)  :: x, s
  478     real(KD),               intent(out) :: alpha
  479     real(KD), dimension(:), intent(in)  :: pmin, pmax
  480   
  481     integer  :: i
  482     real(KD) :: a
  483   
  484     alpha = 1.0D+33
  485     do i = 1, size(x)
  486        ! find alpha for this variable
  487        if ( abs(s(i)) < 1.0D-9 ) cycle
  488        a = (pmax(i) - x(i)) / s(i)
  489        if ( a < 0.0d0 ) a = (pmin(i) - x(i)) / s(i)
  490        if ( a < alpha ) alpha = a
  491     end do
  492   
  493   end subroutine largest_alpha
  494   
  495   
  496   
  497   
  498   
  499   
  500   
  501   
  502   
  503   !**************************************************************************
  504   !**************************************************************************
  505   !**
  506   !**  Nelder-Mead
  507   !**
  508   !**************************************************************************
  509   !**************************************************************************
  510   
  511   
  512   
  513   
  514   !==========================================================================
  515   function simplex_try( p, y, psum, ndim, ihi, fac, pmin, pmax )
  516   !==========================================================================
  517   
  518     integer,  intent(in)                    :: ihi, ndim
  519     real(KD), intent(in)                    :: fac
  520     real(KD), dimension(:,:), intent(inout) :: p
  521     real(KD), dimension(:),   intent(inout) :: psum, y
  522     real(KD), dimension(:),   intent(in)    :: pmin, pmax
  523     real(KD) :: simplex_try
  524   
  525     integer :: j
  526     real(KD) :: fac1, fac2, ytry
  527     real(KD), dimension(ndim) :: ptry
  528   
  529     fac1 = (1.0d0 - fac)/real(ndim,KD)
  530     fac2 = fac1 - fac
  531     do j = 1, ndim
  532        ptry(j) = psum(j) * fac1 - p(ihi,j)*fac2
  533        ! if ptry is outside of the domain, put it back it
  534        if ( ptry(j) < pmin(j) ) ptry(j) = pmin(j)
  535        if ( ptry(j) > pmax(j) ) ptry(j) = pmax(j)
  536     end do
  537   
  538     ytry = E(ptry)
  539     if ( ytry < y(ihi) ) then
  540        y(ihi) = ytry
  541        do j = 1, ndim
  542           psum(j) = psum(j) - p(ihi,j) + ptry(j)
  543           p(ihi,j) = ptry(j)
  544        end do
  545     end if
  546     simplex_try = ytry
  547   
  548   end function simplex_try
  549   
  550   
  551   !==========================================================================
  552   subroutine Simplex_optim( params, optim, iter )
  553   !==========================================================================
  554   ! Performs Nelder-Mead optimization (simplex method)
  555   
  556     type(oparam),           intent(in)    :: params
  557     real(KD), dimension(:), intent(inout) :: optim
  558     integer,                intent(out)   :: iter
  559     
  560     real(KD), dimension(size(optim)+1,size(optim)) :: p
  561     real(KD), dimension(size(optim))               :: x
  562     real(KD), dimension(size(optim)+1)             :: y
  563   
  564     real(KD), parameter :: TINY = 1.0d-10
  565   
  566     real(KD) :: ftol
  567     integer  :: ndim, mp, np
  568   
  569     real(KD) :: rtol, lsum, swap, ysave, ytry, amotry
  570     real(KD), dimension(size(optim)) :: psum, cg
  571     integer :: i, ihi, ilo, inhi, j, m, n, kk
  572   
  573     ftol = params%tol
  574     ndim = size(optim)
  575   
  576     call init_simplex( params%pmin, params%pmax, p, y, optim, params )
  577   
  578     iter = 0
  579   1 continue
  580     do n = 1, ndim
  581        lsum = 0.0d0
  582        do m = 1, ndim + 1
  583           lsum = lsum + p(m,n)
  584        end do
  585        psum(n) = lsum
  586     end do
  587   
  588   2 continue
  589     ilo = 1
  590     if ( y(1) > y(2) ) then
  591        ihi = 1 ; inhi = 2
  592     else
  593        ihi = 2 ; inhi = 1
  594     end if
  595     
  596     do i = 1, ndim + 1
  597        if ( y(i) <= y(ilo) ) ilo = i
  598        if ( y(i) >  y(ihi) ) then
  599           inhi = ihi
  600           ihi = i
  601        else if ( y(i) > y(inhi) ) then
  602           if ( i /= ihi ) inhi = i
  603        end if
  604     end do
  605   
  606     rtol = 2.0d0 * abs( y(ihi)-y(ilo) ) / ( abs(y(ihi)) + abs( y(ilo) ) + TINY )
  607     if ( rtol < ftol ) then
  608        swap = y(1) ; y(1) = y(ilo) ; y(ilo) = swap
  609        do n = 1, ndim
  610           swap = p(1,n) ; p(1,n) = p(ilo,n) ; p(ilo,n) = swap
  611        end do
  612        optim = 0.0d0
  613        do n = 1, size(p,2)
  614           optim(n) = sum( p(:,n) ) / real(size(p,1),KD)
  615        end do
  616        return
  617     end if
  618   
  619     if ( iter >= params%itmax ) then
  620        return
  621     end if
  622   
  623     iter = iter + 2
  624     
  625     ytry = simplex_try( p, y, psum, ndim, ihi, -1.0d0, params%pmin, params%pmax )
  626     if ( ytry <= y(ilo) ) then
  627        ytry = simplex_try( p, y, psum, ndim, ihi, 2.0d0, params%pmin, params%pmax )
  628     else if ( ytry >= y(inhi) ) then
  629        ysave = y(ihi)
  630        ytry = simplex_try( p, y, psum, ndim, ihi, 0.5d0, params%pmin, params%pmax )
  631        if ( ytry >= ysave ) then
  632           do i = 1, ndim + 1
  633              if ( i /= ilo ) then
  634                 do j = 1, ndim
  635                    psum(j) = 0.5d0 * ( p(i,j) + p(ilo,j) )
  636                    p(i,j) = psum(j)
  637                 end do
  638                 y(i) = E(psum)
  639              end if
  640           end do
  641           iter = iter + ndim
  642           goto 1
  643        end if
  644     else
  645        iter = iter - 1
  646     end if
  647   
  648     if ( params%n_cons == 0 ) then
  649        call simplex_cg( p, cg )
  650        write( out, '(6es20.10)' ) cg(1:min(params%n,6))
  651     end if
  652   
  653     goto 2
  654   
  655   end subroutine Simplex_optim
  656   
  657   
  658   
  659   !==========================================================================
  660   subroutine simplex_cg( p, c )
  661   !==========================================================================
  662   ! returns C the centre of the simplex P
  663   
  664     real(KD), dimension(:,:), intent(in)  :: p
  665     real(KD), dimension(:),   intent(out) :: c
  666     
  667     integer :: i
  668   
  669     c = 0.0d0
  670     do i = 1, size(p,2)
  671        c(i) = sum( p(:,i) ) / real(size(p,1),KD)
  672     end do
  673   
  674   end subroutine simplex_cg
  675   
  676   
  677   
  678   !==========================================================================
  679   subroutine init_simplex( pmin, pmax, p, y, p0, op ) 
  680   !==========================================================================
  681   ! initializes a simplex for the Nelder-Mead optimization algorithm
  682   ! PMIN and PMAX are used to devise a characteristic length
  683   
  684     real(KD), dimension(:),   intent(in)  :: pmin, pmax, p0
  685     real(KD), dimension(:,:), intent(out) :: p
  686     real(KD), dimension(:),   intent(out) :: y
  687     type(oparam),             intent(in)  :: op
  688   
  689     integer  :: i
  690     real(KD) :: lambda
  691   
  692     p(1,:) = p0
  693     y(1)   = E( p0 )
  694     do i = 1, size(p0)
  695        lambda = (pmax(i)-pmin(i)) * op%simplex_lambda
  696        p(1+i,:) = p0
  697        p(1+i,i) = p(1+i,i) + lambda
  698        y(1+i)   = E( p(1+i,:) )
  699     end do
  700   
  701   end subroutine init_simplex
  702   
  703   
  704   
  705   
  706   !**************************************************************************
  707   !**************************************************************************
  708   !**
  709   !**  Genetic algorithm
  710   !**
  711   !**************************************************************************
  712   !**************************************************************************
  713   
  714   !==========================================================================
  715   subroutine GA_optim( params, optim, iter )
  716   !==========================================================================
  717   ! uses a genetic algorithm to perform optimizsation of function E
  718   
  719     type(oparam),           intent(in)    :: params
  720     real(KD), dimension(:), intent(inout) :: optim
  721     integer,                intent(out)   :: iter
  722     
  723     real(KD), dimension(size(optim),params%GA_popsize) :: pop
  724     real(KD), dimension(params%GA_popsize) :: score
  725   
  726     real(KD), dimension(size(optim)) :: p0
  727     
  728   
  729     integer :: i, j, ibest, k, l, ip
  730     integer, dimension(1) :: iloc
  731     real(KD) :: best_score
  732   
  733     ! initialize population
  734     best_score = 1.0d20
  735     do i = 1, params%GA_popsize
  736        call initvar( pop(:,i), params%pmin, params%pmax )
  737        score(i) = E( pop(:,i) )
  738        if ( score(i) < best_score ) then
  739           best_score = score(i)
  740           ibest = i
  741        end if
  742     end do
  743   
  744     ! perform optimization
  745     ip = 1
  746     do i = 1, params%itmax
  747        ! mutate
  748        do j = 1, params%GA_popsize
  749           if ( rRan() <= params%prob_mut ) then
  750              call mutate( pop(:,j), params%pmin, params%pmax )
  751              score(j) = E( pop(:,j) )
  752              if ( score(j) < best_score ) then
  753                 best_score = score(j)
  754                 ibest = j
  755              end if
  756           end if
  757        end do
  758        ! cross-over
  759        do j = 1, params%GA_popsize
  760           if ( rRan() <= params%prob_cross ) then
  761              k = iranrange( 1, params%GA_popsize )
  762              if ( k /= j ) then
  763                 call crossover( pop(:,j), pop(:,k), params%pmin, params%pmax )
  764                 score(j) = E( pop(:,j) )
  765                 if ( score(j) < best_score ) then
  766                    best_score = score(j)
  767                    ibest = j
  768                 end if
  769                 score(k) = E( pop(:,k) )
  770                 if ( score(k) < best_score ) then
  771                    best_score = score(k)
  772                    ibest = k
  773                 end if
  774              end if
  775           end if
  776        end do
  777        ! tournament
  778        do j = 1, params%GA_popsize
  779           k = iranrange( 1, params%GA_popsize )
  780           l = iranrange( 1, params%GA_popsize )
  781           do while ( k == l )
  782              l = iranrange( 1, params%GA_popsize)
  783           end do
  784           if ( score(k) < score(l) ) then
  785              score(l) = score(k)
  786              pop(:,l) = pop(:,k)
  787           else
  788              score(k) = score(l)
  789              pop(:,k) = pop(:,l)
  790           end if
  791        end do
  792        iloc = minloc( score )
  793        ibest = iloc(1)
  794        if ( params%n_cons == 0 ) then
  795           if ( ip == i ) then
  796              write( out, '(6es20.10)' ) pop(1:min(6,params%n),ibest)
  797              ip = 2 * ip
  798           end if
  799        end if
  800     end do
  801   
  802     iter = params%itmax
  803     optim = pop(:,ibest)
  804   
  805   end subroutine GA_optim
  806   
  807   
  808   
  809   
  810   !==========================================================================
  811   subroutine mutate( var, pmin, pmax )
  812   !==========================================================================
  813   ! updates the variables VAR by changing one, so as to satisfy constraints
  814   ! on their range
  815   
  816     real(KD), dimension(:), intent(inout) :: var
  817     real(KD), dimension(:), intent(in)    :: pmin, pmax
  818   
  819     integer :: i
  820   
  821     i = iranrange( 1, size(var) )
  822   
  823     var(i) = pmin(i) + rRan() * (pmax(i) - pmin(i))
  824   
  825   end subroutine mutate
  826   
  827   
  828   
  829   !==========================================================================
  830   subroutine crossover( var1, var2, pmin, pmax )
  831   !==========================================================================
  832   ! updates the variables VAR1 and VAR2 by exchanging some of them
  833   
  834   
  835     real(KD), dimension(:), intent(inout) :: var1, var2
  836     real(KD), dimension(:), intent(in)    :: pmin, pmax
  837   
  838     integer  :: i1, i2, j
  839     real(KD) :: tra
  840   
  841     i1 = iranrange( 1, size(var1) )
  842     i2 = iranrange( 1, size(var1) )
  843     
  844     if ( i1 > i2 ) then
  845        j = i1 ; i1 = i2 ; i2 = j
  846     end if
  847   
  848     do j = i1, i2
  849        tra = var1(j) ; var1(j) = var2(j) ; var2(j) = tra
  850     end do
  851   
  852   end subroutine crossover
  853   
  854   
  855   !==========================================================================
  856   subroutine initvar( var, pmin, pmax )
  857   !==========================================================================
  858   ! updates the variables VAR by changing one, so as to satisfy constraints
  859   ! on their range
  860   
  861     real(KD), dimension(:), intent(inout) :: var
  862     real(KD), dimension(:), intent(in)    :: pmin, pmax
  863   
  864     integer :: i
  865   
  866     do i = 1, size(var)
  867        var(i) = pmin(i) + rRan() * (pmax(i) - pmin(i))
  868     end do
  869   
  870   end subroutine initvar
  871   
  872   
  873   
  874   
  875   
  876   !**************************************************************************
  877   !**************************************************************************
  878   !**
  879   !**  Simulated annealing
  880   !**
  881   !**************************************************************************
  882   !**************************************************************************
  883   
  884   
  885   
  886   !==========================================================================
  887   subroutine SA_optim( params, optim, iter )
  888   !==========================================================================
  889   ! uses simulated annealing algorithm to perform optimizsation of function E
  890   
  891     type(oparam),           intent(in)    :: params
  892     real(KD), dimension(:), intent(inout) :: optim
  893     integer,                intent(out)   :: iter
  894   
  895     real(KD), dimension(size(optim)) :: path, bestp, new_path
  896     real(KD) :: cost, beta, delta, delta_cost, best
  897   
  898     integer :: i, j, c
  899   
  900     iter = params%itmax
  901     cost = E(optim)
  902   
  903     path = optim
  904   
  905     beta = 1.0D2 ; j = 1 ; c = 0 ; best = 1.0D10
  906     delta = (1.0D4/1.0D2)**(1.0D0/real(iter,KD))
  907   
  908     do i = 1, iter
  909        beta = beta * delta
  910        call SA_mutate( cost, delta_cost, path, new_path, params%pmin, params%pmax )
  911        if ( delta_cost < 0.0d0 .or. rRan() < exp( -beta*delta_cost ) ) then
  912           cost = cost + delta_cost
  913           path = new_path
  914           c = c + 1
  915           if ( cost < best ) then
  916              best  = cost
  917              bestp = path
  918           end if
  919        end if
  920        if ( params%n_cons == 0 ) then
  921           if ( mod(i,iter/100) == 0 ) &
  922                write( out, '(6es20.10)' ) bestp(1:min(6,params%n)) 
  923        end if
  924     end do
  925     
  926     optim = bestp
  927   
  928   end subroutine SA_optim
  929   
  930   
  931   !==========================================================================
  932   subroutine SA_mutate( cost, delta_cost, path, new_path, pmin, pmax )
  933   !==========================================================================
  934   
  935     real(KD), intent(in)  :: cost
  936     real(KD), intent(out) :: delta_cost
  937     real(KD), dimension(:), intent(in)  :: path
  938     real(KD), dimension(:), intent(out) :: new_path
  939     real(KD), dimension(:), intent(in)  :: pmin, pmax
  940   
  941     integer :: i
  942   
  943     new_path = path
  944     do i = 1, size(path)
  945        if ( rRan() > 0.5d0 ) then
  946           new_path(i) = pmin(i) + rRan() * (pmax(i) - pmin(i))
  947        end if
  948     end do
  949     delta_cost = E(new_path) - cost
  950   
  951   end subroutine SA_mutate
  952   
  953   
  954   
  955   !**************************************************************************
  956   !**************************************************************************
  957   !**
  958   !**  Utility routines
  959   !**
  960   !**************************************************************************
  961   !**************************************************************************
  962   
  963   
  964   !==========================================================================
  965   function NormL2( H )
  966   !==========================================================================
  967   
  968     real(KD), dimension(:), intent(in) :: H
  969     real(KD) :: NormL2
  970   
  971     integer :: i
  972   
  973     NormL2 = 0.0d0
  974     do i = 1, size(H)
  975        NormL2 = NormL2 + H(i)*H(i)
  976     end do
  977   
  978   end function NormL2
  979   
  980   
  981   !==========================================================================
  982   function modulus( x )
  983   !==========================================================================
  984   
  985     real(KD), dimension(:), intent(in) :: x
  986     real(KD) :: modulus
  987   
  988     modulus = sqrt(NormL2(x))
  989   
  990   end function modulus
  991   
  992   
  993   
  994   !==========================================================================
  995   subroutine init( op, p )
  996   !==========================================================================
  997   ! initialize an optimization by reading an input file
  998   
  999     type(oparam),           intent(out) :: op
 1000     real(KD), dimension(:), pointer     :: p
 1001   
 1002     logical            :: fexist
 1003     character(len=128) :: line
 1004     integer            :: i
 1005   
 1006     real(KD), dimension(:), pointer :: pmin, pmax
 1007   
 1008     inquire( file = 'optim.in', exist = fexist )
 1009     if ( .not. fexist ) stop 'ERR: cannot read optim.in.'
 1010   
 1011     open( unit = 9, file = 'optim.in' )
 1012     read( 9, '(a)' ) line
 1013     read( 9, * ) op%n ! number of variable
 1014     if ( op%n <= 0 ) stop 'ERR: invalid number of variables.'
 1015   
 1016     read( 9, '(a)' ) line
 1017     read( 9, * ) op%m ! number of parameter
 1018     if ( op%m < 0 ) stop 'ERR: invalid number of parameters.'
 1019   
 1020   
 1021     read( 9, '(a)' ) line
 1022     read( 9, * ) op%n_cons ! number of constraint
 1023     op%use_cons = .false.
 1024     if ( op%n_cons > 0 ) op%use_cons = .true.
 1025   
 1026     read( 9, '(a)' ) line
 1027     allocate( pmin(op%n), pmax(op%n) )
 1028     do i = 1, op%n
 1029        read( 9, * ) pmin(i), pmax(i)
 1030     end do
 1031     op%pmin => pmin
 1032     op%pmax => pmax
 1033   
 1034     ! initial solution
 1035     read( 9, '(a)' ) line
 1036     allocate( p(op%N) )
 1037     do i = 1, op%N
 1038        read( 9, * ) p(i)
 1039     end do
 1040   
 1041     ! for the unconstrained problem: method, itmax, tol
 1042     read( 9, '(a)' ) line
 1043     read( 9, '(a)' ) op%method
 1044     read( 9, * ) op%itmax
 1045     read( 9, * ) op%tol
 1046   
 1047     ! for the constrained problem: method, citmax, ctol
 1048     read( 9, '(a)' ) line
 1049     read( 9, '(a)' ) op%cons_method
 1050     read( 9, * ) op%itmaxcons
 1051     read( 9, * ) op%ctol
 1052   
 1053     ! Now come various parameters
 1054   
 1055     ! line search: sigma and rho
 1056     read( 9, '(a)' ) line
 1057     read( 9, * ) op%ls_sigma
 1058     read( 9, * ) op%ls_rho
 1059   
 1060     ! simplex: scaling (lambda) for initial simplex
 1061     read( 9, '(a)' ) line
 1062     read( 9, * ) op%simplex_lambda
 1063   
 1064     ! genetic algorithm: population size, probabilities of mutation and crossover
 1065     read( 9, '(a)' ) line
 1066     read( 9, * ) op%GA_popsize
 1067     read( 9, * ) op%prob_mut
 1068     read( 9, * ) op%prob_cross
 1069   
 1070     ! simulated annealing: initial beta
 1071     read( 9, '(a)' ) line
 1072     read( 9, * ) op%SA_beta
 1073   
 1074     ! constrained optimization: alpha, beta, lambda_0
 1075     read( 9, '(a)' ) line
 1076     read( 9, * ) op%calpha
 1077     read( 9, * ) op%cbeta
 1078     read( 9, * ) op%clambda
 1079   
 1080     ! seed for the PRNG
 1081     read( 9, '(a)' ) line
 1082     read( 9, * ) op%seed
 1083     close( 9 )
 1084   
 1085     ! display these values to screen
 1086     print *
 1087     print *, 'O P T I M I Z A T I O N  summary'
 1088     print *
 1089     write( *, '(a32,i5)' ) 'problem size', op%N
 1090     write( *, '(a32,i5)' ) 'number of parameters', op%M
 1091     write( *, '(a32,i5)' ) 'number of constraints', op%n_cons
 1092   
 1093     print *
 1094     print *, 'Initial solution:'
 1095     write( *, '(10es15.4)' ) p(1:min(op%n,10))
 1096   
 1097     print *
 1098     print *, 'Solution range:'
 1099     write( *, '(10es15.4)' ) pmin(1:min(op%n,10))
 1100     write( *, '(10es15.4)' ) pmax(1:min(op%n,10))
 1101   
 1102     print *
 1103     print *, 'Problem/method: '
 1104     if ( op%n_cons > 0 ) then
 1105        print *, '   constrained'
 1106        print *, '  ', op%cons_method
 1107        print *, '   unconstrained solver: ', op%method
 1108     else
 1109        print *, '   unconstrained'
 1110        print *, '   method: ', op%method
 1111     end if
 1112   
 1113     print *
 1114     print *, 'Numerical parameters'
 1115     write( *, '(a32,i8)' ) 'number of iterations', op%itmax
 1116     write( *, '(a32,es10.2)' ) 'tolerance (unconstrained)', op%tol
 1117   
 1118     if ( op%n_cons > 0 ) then
 1119          write( *, '(a32,es10.2)' ) 'tolerance (constrained)', op%ctol
 1120          write( *, '(a32,es10.2)' ) 'alpha (constr.)', op%calpha
 1121          write( *, '(a32,es10.2)' ) 'beta (constr.)', op%cbeta
 1122          if ( index( op%cons_method, 'multipl' ) > 0 ) then
 1123             write( *, '(a32,es10.2)' ) 'initial multiplier', op%clambda
 1124          end if
 1125          print *
 1126     end if
 1127   
 1128     if ( index( op%method, 'BFGS' ) > 0 .or. &
 1129          index( op%method, 'steepest' ) > 0 ) then
 1130        write( *, '(a32,es10.2)' ) 'line_search, sigma', op%ls_sigma
 1131        write( *, '(a32,es10.2)' ) 'line_search, rho', op%ls_rho
 1132     end if
 1133     if ( index( op%method, 'simplex' ) > 0 ) then
 1134        write( *, '(a32,es10.2)' ) 'simplex init, lambda', op%simplex_lambda
 1135     end if
 1136     if ( index( op%method, 'genetic' ) > 0 ) then
 1137        write( *, '(a32,i8)' ) 'population size', op%GA_popsize
 1138        write( *, '(a32,es10.2)' ) 'mutation probability', op%prob_mut
 1139        write( *, '(a32,es10.2)' ) 'cross-over probability', op%prob_cross
 1140     end if
 1141     if ( index( op%method, 'annealing' ) > 0 ) then
 1142        write( *, '(a32,es10.2)' ) 'initial beta', op%SA_beta
 1143     end if
 1144   
 1145     print *
 1146   
 1147   end subroutine init
 1148   
 1149   end program optimisation
List of routines:
BFGS_optim
COptim
GA_optim
SA_mutate
SA_optim
SD_optim
Simplex_optim
UOptim
crossover
init
init_simplex
initvar
largest_alpha
line_search
mutate
simplex_cg
solve
E
GradE
H
NormL2
modulus
simplex_try
(list of modules)