octave_packages/optim-1.2.0/de_min.m

   1 ## Copyright (C) 1996, 1997 R. Storn
   2 ## Copyright (C) 2009-2010 Christian Fischer <cfischer@itm.uni-stuttgart.de>
   3 ##
   4 ## This program is free software; you can redistribute it and/or modify it under
   5 ## the terms of the GNU General Public License as published by the Free Software
   6 ## Foundation; either version 3 of the License, or (at your option) any later
   7 ## version.
   8 ##
   9 ## This program is distributed in the hope that it will be useful, but WITHOUT
  10 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  11 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
  12 ## details.
  13 ##
  14 ## You should have received a copy of the GNU General Public License along with
  15 ## this program; if not, see <http://www.gnu.org/licenses/>.
  16
  17 ## de_min: global optimisation using differential evolution
  18 ##
  19 ## Usage: [x, obj_value, nfeval, convergence] = de_min(fcn, control)
  20 ##
  21 ## minimization of a user-supplied function with respect to x(1:D),
  22 ## using the differential evolution (DE) method based on an algorithm
  23 ## by  Rainer Storn (http://www.icsi.berkeley.edu/~storn/code.html)
  24 ## See: http://www.softcomputing.net/tevc2009_1.pdf
  25 ##
  26 ##
  27 ## Arguments:
  28 ## ---------------
  29 ## fcn        string : Name of function. Must return a real value
  30 ## control    vector : (Optional) Control variables, described below
  31 ##         or struct
  32 ##
  33 ## Returned values:
  34 ## ----------------
  35 ## x          vector : parameter vector of best solution
  36 ## obj_value  scalar : objective function value of best solution
  37 ## nfeval     scalar : number of function evaluations
  38 ## convergence       : 1 = best below value to reach (VTR)
  39 ##                     0 = population has reached defined quality (tol)
  40 ##                    -1 = some values are close to constraints/boundaries
  41 ##                    -2 = max number of iterations reached (maxiter)
  42 ##                    -3 = max number of functions evaluations reached (maxnfe)
  43 ##
  44 ## Control variable:   (optional) may be named arguments (i.e. "name",value
  45 ## ----------------    pairs), a struct, or a vector, where
  46 ##                     NaN's are ignored.
  47 ##
  48 ## XVmin        : vector of lower bounds of initial population
  49 ##                *** note: by default these are no constraints ***
  50 ## XVmax        : vector of upper bounds of initial population
  51 ## constr       : 1 -> enforce the bounds not just for the initial population
  52 ## const        : data vector (remains fixed during the minimization)
  53 ## NP           : number of population members
  54 ## F            : difference factor from interval [0, 2]
  55 ## CR           : crossover probability constant from interval [0, 1]
  56 ## strategy     : 1 --> DE/best/1/exp           7 --> DE/best/1/bin
  57 ##                2 --> DE/rand/1/exp           8 --> DE/rand/1/bin
  58 ##                3 --> DE/target-to-best/1/exp 9 --> DE/target-to-best/1/bin
  59 ##                4 --> DE/best/2/exp           10--> DE/best/2/bin
  60 ##                5 --> DE/rand/2/exp           11--> DE/rand/2/bin
  61 ##                6 --> DEGL/SAW/exp            else  DEGL/SAW/bin
  62 ## refresh      : intermediate output will be produced after "refresh"
  63 ##                iterations. No intermediate output will be produced
  64 ##                if refresh is < 1
  65 ## VTR          : Stopping criterion: "Value To Reach"
  66 ##                de_min will stop when obj_value <= VTR.
  67 ##                Use this if you know which value you expect.
  68 ## tol          : Stopping criterion: "tolerance"
  69 ##                stops if (best-worst)/max(1,worst) < tol
  70 ##                This stops basically if the whole population is "good".
  71 ## maxnfe       : maximum number of function evaluations
  72 ## maxiter      : maximum number of iterations (generations)
  73 ##
  74 ##       The algorithm seems to work well only if [XVmin,XVmax] covers the
  75 ##       region where the global minimum is expected.
  76 ##       DE is also somewhat sensitive to the choice of the
  77 ##       difference factor F. A good initial guess is to choose F from
  78 ##       interval [0.5, 1], e.g. 0.8.
  79 ##       CR, the crossover probability constant from interval [0, 1]
  80 ##       helps to maintain the diversity of the population and is
  81 ##       rather uncritical but affects strongly the convergence speed.
  82 ##       If the parameters are correlated, high values of CR work better.
  83 ##       The reverse is true for no correlation.
  84 ##       Experiments suggest that /bin likes to have a slightly
  85 ##       larger CR than /exp.
  86 ##       The number of population members NP is also not very critical. A
  87 ##       good initial guess is 10*D. Depending on the difficulty of the
  88 ##       problem NP can be lower than 10*D or must be higher than 10*D
  89 ##       to achieve convergence.
  90 ##
  91 ## Default Values:
  92 ## ---------------
  93 ## XVmin = [-2];
  94 ## XVmax = [ 2];
  95 ## constr= 0;
  96 ## const = [];
  97 ## NP    = 10 *D
  98 ## F     = 0.8;
  99 ## CR    = 0.9;
 100 ## strategy = 12;
 101 ## refresh  = 0;
 102 ## VTR   = -Inf;
 103 ## tol   = 1.e-3;
 104 ## maxnfe  = 1e6;
 105 ## maxiter = 1000;
 106 ##
 107 ##
 108 ## Example to find the minimum of the Rosenbrock saddle:
 109 ## ----------------------------------------------------
 110 ## Define f as:
 111 ##                    function result = f(x);
 112 ##                      result = 100 * (x(2) - x(1)^2)^2 + (1 - x(1))^2;
 113 ##                    end
 114 ## Then type:
 115 ##
 116 ##      ctl.XVmin = [-2 -2];
 117 ##      ctl.XVmax = [ 2  2];
 118 ##      [x, obj_value, nfeval, convergence] = de_min (@f, ctl);
 119 ##
 120 ## Keywords: global-optimisation optimisation minimisation
 121
 122 function [bestmem, bestval, nfeval, convergence] = de_min(fcn, varargin)
 123
 124 ## Default options
 125 XVmin = [-2 ];
 126 XVmax = [ 2 ];
 127 constr= 0;
 128 const = [];
 129 NP    = 0;         # NP will be set later
 130 F     = 0.8;
 131 CR    = 0.9;
 132 strategy = 12;
 133 refresh  = 0;
 134 VTR   = -Inf;
 135 tol   = 1.e-3;
 136 maxnfe  = 1e6;
 137 maxiter = 1000;
 138
 139 ## ------------ Check input variables (ctl) --------------------------------
 140 if nargin >= 2,                 # Read control arguments
 141   va_arg_cnt = 1;
 142   if nargin > 2,
 143     ctl = struct (varargin{:});
 144   else
 145     ctl = varargin{va_arg_cnt++};
 146   end
 147   if isnumeric (ctl)
 148     if length (ctl)>=1 && !isnan (ctl(1)), XVmin  = ctl(1); end
 149     if length (ctl)>=2 && !isnan (ctl(2)), XVmax  = ctl(2); end
 150     if length (ctl)>=3 && !isnan (ctl(3)), constr = ctl(3); end
 151     if length (ctl)>=4 && !isnan (ctl(4)), const  = ctl(4); end
 152     if length (ctl)>=5 && !isnan (ctl(5)), NP  = ctl(5); end
 153     if length (ctl)>=6 && !isnan (ctl(6)), F   = ctl(6); end
 154     if length (ctl)>=7 && !isnan (ctl(7)), CR  = ctl(7); end
 155     if length (ctl)>=8 && !isnan (ctl(8)), strategy = ctl(8); end
 156     if length (ctl)>=9 && !isnan (ctl(9)), refresh  = ctl(9); end
 157     if length (ctl)>=10&& !isnan (ctl(10)), VTR   = ctl(10); end
 158     if length (ctl)>=11&& !isnan (ctl(11)), tol   = ctl(11); end
 159     if length (ctl)>=12&& !isnan (ctl(12)), maxnfe  = ctl(12); end
 160     if length (ctl)>=13&& !isnan (ctl(13)), maxiter = ctl(13); end
 161   else
 162     if isfield (ctl,"XVmin") && !isnan (ctl.XVmin), XVmin=ctl.XVmin; end
 163     if isfield (ctl,"XVmax") && !isnan (ctl.XVmax), XVmax=ctl.XVmax; end
 164     if isfield (ctl,"constr")&& !isnan (ctl.constr), constr=ctl.constr; end
 165     if isfield (ctl,"const") && !isnan (ctl.const), const=ctl.const; end
 166     if isfield (ctl, "NP" ) && ! isnan (ctl.NP  ), NP  = ctl.NP  ; end
 167     if isfield (ctl, "F"  ) && ! isnan (ctl.F   ), F   = ctl.F   ; end
 168     if isfield (ctl, "CR" ) && ! isnan (ctl.CR  ), CR  = ctl.CR  ; end
 169     if isfield (ctl, "strategy") && ! isnan (ctl.strategy),
 170       strategy = ctl.strategy  ; end
 171     if isfield (ctl, "refresh") && ! isnan (ctl.refresh),
 172       refresh  = ctl.refresh   ; end
 173     if isfield (ctl, "VTR") && ! isnan (ctl.VTR ), VTR = ctl.VTR ; end
 174     if isfield (ctl, "tol") && ! isnan (ctl.tol ), tol = ctl.tol ; end
 175     if isfield (ctl, "maxnfe")  && ! isnan (ctl.maxnfe)
 176       maxnfe = ctl.maxnfe;
 177     end
 178     if isfield (ctl, "maxiter") && ! isnan (ctl.maxiter)
 179       maxiter = ctl.maxiter;
 180     end
 181   end
 182 end
 183
 184 ## set dimension D and population size NP
 185 D  = length (XVmin);
 186 if (NP == 0);
 187   NP = 10 * D;
 188 end
 189
 190 ## -------- do a few sanity checks --------------------------------
 191 if (length (XVmin) != length (XVmax))
 192   error("Length of upper and lower bounds does not match.")
 193 end
 194 if (NP < 5)
 195   error("Population size NP must be bigger than 5.")
 196 end
 197 if ((F <= 0) || (F > 2))
 198   error("Difference Factor F out of range (0,2].")
 199 end
 200 if ((CR < 0) || (CR > 1))
 201   error("CR value out of range [0,1].")
 202 end
 203 if (maxiter <= 0)
 204   error("maxiter must be positive.")
 205 end
 206 if (maxnfe <= 0)
 207   error("maxnfe must be positive.")
 208 end
 209 refresh = floor(abs(refresh));
 210
 211 ## ----- Initialize population and some arrays --------------------------
 212
 213 pop = zeros(NP,D);  # initialize pop
 214
 215 ## pop is a matrix of size NPxD. It will be initialized with
 216 ## random values between the min and max values of the parameters
 217 for i = 1:NP
 218    pop(i,:) = XVmin + rand (1,D) .* (XVmax - XVmin);
 219 end
 220
 221 ## initialise the weighting factors between 0.0 and 1.0
 222 w  = rand (NP,1);
 223 wi = w;
 224
 225 popold    = zeros (size (pop));   # toggle population
 226 val       = zeros (1, NP);        # create and reset the "cost array"
 227 bestmem   = zeros (1, D);         # best population member ever
 228 bestmemit = zeros (1 ,D);         # best population member in iteration
 229 nfeval    = 0;                    # number of function evaluations
 230
 231 ## ------ Evaluate the best member after initialization ------------------
 232
 233 ibest   = 1;                      # start with first population member
 234 val(1)  = feval (fcn, [pop(ibest,:) const]);
 235 bestval = val(1);                 # best objective function value so far
 236 bestw   = w(1);                   # weighting of best design so far
 237 for i = 2:NP                      # check the remaining members
 238   val(i) = feval (fcn, [pop(i,:) const]);
 239   if (val(i) < bestval)           # if member is better
 240      ibest   = i;                 # save its location
 241      bestval = val(i);
 242      bestw   = w(i);
 243   end
 244 end
 245 nfeval  = nfeval + NP;
 246 bestmemit = pop(ibest,:);         # best member of current iteration
 247 bestvalit = bestval;              # best value of current iteration
 248
 249 bestmem = bestmemit;              # best member ever
 250
 251 ## ------ DE - Minimization ---------------------------------------
 252 ## popold is the population which has to compete. It is static
 253 ## through one iteration. pop is the newly emerging population.
 254
 255 bm_n= zeros (NP, D);            # initialize bestmember matrix in neighbourh.
 256 lpm1= zeros (NP, D);            # initialize local population matrix 1
 257 lpm1= zeros (NP, D);            # initialize local population matrix 2
 258 rot = 0:1:NP-1;                 # rotating index array (size NP)
 259 rotd= 0:1:D-1;                  # rotating index array (size D)
 260
 261 iter = 1;
 262 while ((iter < maxiter) && (nfeval < maxnfe) &&  (bestval > VTR)  && ...
 263        ((abs (max (val) - bestval) / max (1, abs (max (val))) > tol)))
 264   popold = pop;                   # save the old population
 265   wold   = w;                     # save the old weighting factors
 266
 267   ind = randperm (4);             # index pointer array
 268
 269   a1  = randperm (NP);            # shuffle locations of vectors
 270   rt = rem (rot + ind(1), NP);    # rotate indices by ind(1) positions
 271   a2  = a1(rt+1);                 # rotate vector locations
 272   rt = rem (rot + ind(2), NP);
 273   a3  = a2(rt+1);
 274   rt = rem (rot  +ind(3), NP);
 275   a4  = a3(rt+1);
 276   rt = rem (rot + ind(4), NP);
 277   a5  = a4(rt+1);
 278
 279   pm1 = popold(a1,:);             # shuffled population 1
 280   pm2 = popold(a2,:);             # shuffled population 2
 281   pm3 = popold(a3,:);             # shuffled population 3
 282   w1  = wold(a4);                 # shuffled weightings 1
 283   w2  = wold(a5);                 # shuffled weightings 2
 284
 285   bm = repmat (bestmemit, NP, 1); # population filled with the best member
 286                                   # of the last iteration
 287   bw = repmat (bestw, NP, 1);     # the same for the weighting of the best
 288
 289   mui = rand (NP, D) < CR;        # mask for intermediate population
 290                                   # all random numbers < CR are 1, 0 otherwise
 291
 292   if (strategy > 6)
 293     st = strategy - 6;                # binomial crossover
 294   else
 295     st = strategy;                        # exponential crossover
 296     mui = sort (mui');            # transpose, collect 1's in each column
 297     for i = 1:NP
 298       n = floor (rand * D);
 299       if (n > 0)
 300          rtd = rem (rotd + n, D);
 301          mui(:,i) = mui(rtd+1,i); #rotate column i by n
 302       endif
 303     endfor
 304     mui = mui';                                   # transpose back
 305   endif
 306   mpo = mui < 0.5;                # inverse mask to mui
 307
 308   if (st == 1)                            # DE/best/1
 309     ui = bm + F*(pm1 - pm2);              # differential variation
 310   elseif (st == 2)                        # DE/rand/1
 311     ui = pm3 + F*(pm1 - pm2);             # differential variation
 312   elseif (st == 3)                        # DE/target-to-best/1
 313     ui = popold + F*(bm-popold) + F*(pm1 - pm2);
 314   elseif (st == 4)                        # DE/best/2
 315     pm4 = popold(a4,:);                 # shuffled population 4
 316     pm5 = popold(a5,:);                 # shuffled population 5
 317     ui = bm + F*(pm1 - pm2 + pm3 - pm4);  # differential variation
 318   elseif (st == 5)                        # DE/rand/2
 319     pm4 = popold(a4,:);                 # shuffled population 4
 320     pm5 = popold(a5,:);                 # shuffled population 5
 321     ui = pm5 + F*(pm1 - pm2 + pm3 - pm4); # differential variation
 322   else                                    # DEGL/SAW
 323     ## The DEGL/SAW method is more complicated.
 324         ## We have to generate a neighbourhood first.
 325         ## The neighbourhood size is 10% of NP and at least 1.
 326         nr = max (1, ceil ((0.1*NP -1)/2));  # neighbourhood radius
 327         ## FIXME: I don't know how to vectorise this. - if possible
 328         for i = 1:NP
 329           neigh_ind = i-nr:i+nr;              # index range of neighbourhood
 330           neigh_ind = neigh_ind + ((neigh_ind <= 0)-(neigh_ind > NP))*NP;
 331                                           # do wrap around
 332           [x, ix] = min (val(neigh_ind));     # find the local best and its index
 333           bm_n(i,:) = popold(neigh_ind(ix),:); # copy the data from the local best
 334           neigh_ind(nr+1) = [];                 # remove "i"
 335           pq = neigh_ind(randperm (length (neigh_ind)));
 336                                         # permutation of the remaining ind.
 337           lpm1(i,:) = popold(pq(1),:);        # create the local pop member matrix
 338           lpm2(i,:) = popold(pq(2),:);        # for the random point p,q
 339         endfor
 340     ## calculate the new weihting factors
 341     wi = wold + F*(bw - wold) + F*(w1 - w2); # use DE/target-to-best/1/nocross
 342                                           # for optimisation of weightings
 343     ## fix bounds for weightings
 344     o = ones (NP, 1);
 345     wi = sort ([0.05*o, wi, 0.95*o],2)(:,2); # sort and take the second column
 346     ## fill weighting matrix
 347     wm = repmat (wi, 1, D);
 348     li = popold + F*(bm_n- popold) + F*(lpm1 - lpm2);
 349     gi = popold + F*(bm  - popold) + F*(pm1 - pm2);
 350     ui = wm.*gi + (1-wm).*li;             # combine global and local part
 351   endif
 352   ## crossover
 353   ui = popold.*mpo + ui.*mui;
 354
 355   ## enforce initial bounds/constraints if specified
 356   if (constr == 1)
 357     for i = 1:NP
 358       ui(i,:) = max (ui(i,:), XVmin);
 359       ui(i,:) = min (ui(i,:), XVmax);
 360     end
 361   end
 362
 363   ## ----- Select which vectors are allowed to enter the new population ------
 364   for i = 1:NP
 365     tempval = feval (fcn, [ui(i,:) const]);   # check cost of competitor
 366     if (tempval <= val(i))  # if competitor is better
 367        pop(i,:) = ui(i,:);  # replace old vector with new one
 368        val(i)   = tempval;  # save value in "cost array"
 369        w(i)     = wi(i);    # save the weighting factor
 370
 371        ## we update bestval only in case of success to save time
 372        if (tempval <= bestval)     # if competitor better than the best one ever
 373           bestval = tempval;      # new best value
 374           bestmem = ui(i,:);      # new best parameter vector ever
 375           bestw   = wi(i);        # save best weighting
 376        end
 377     end
 378   endfor #---end for i = 1:NP
 379
 380   nfeval  = nfeval + NP;     # increase number of function evaluations
 381
 382   bestmemit = bestmem;       # freeze the best member of this iteration for the
 383                              # coming iteration. This is needed for some of the
 384                              # strategies.
 385
 386   ## ---- Output section ----------------------------------------------------
 387
 388   if (refresh > 0)
 389     if (rem (iter, refresh) == 0)
 390        printf ('Iteration: %d,  Best: %8.4e,  Worst: %8.4e\n', ...
 391                 iter, bestval, max(val));
 392        for n = 1:D
 393          printf ('x(%d) = %e\n', n, bestmem(n));
 394        end
 395     end
 396   end
 397
 398   iter = iter + 1;
 399 endwhile #---end while ((iter < maxiter) ...
 400
 401 ## check that all variables are well within bounds/constraints
 402 boundsOK = 1;
 403 for i = 1:NP
 404   range = XVmax - XVmin;
 405   if (ui(i,:) < XVmin + 0.01*range)
 406     boundsOK = 0;
 407   end
 408   if (ui(i,:) > XVmax - 0.01*range)
 409     boundsOK = 0;
 410   end
 411 end
 412
 413 ## create the convergence result
 414 if (bestval <= VTR)
 415   convergence = 1;
 416 elseif (abs (max (val) - bestval) / max (1, abs (max (val))) <= tol)
 417   convergence = 0;
 418 elseif (boundsOK == 0)
 419   convergence = -1;
 420 elseif (iter >= maxiter)
 421   convergence = -2;
 422 elseif (nfeval >= maxnfe)
 423   convergence = -3;
 424 end
 425
 426 endfunction
 427
 428 %!function result = f(x);
 429 %!  result = 100 * (x(2) - x(1)^2)^2 + (1 - x(1))^2;
 430 %!test
 431 %! tol = 1.0e-4;
 432 %! ctl.tol   = 0.0;
 433 %! ctl.VTR   = 1.0e-6;
 434 %! ctl.XVmin = [-2 -2];
 435 %! ctl.XVmax = [ 2  2];
 436 %! rand("state", 11)
 437 %! [x, obj_value, nfeval, convergence] = de_min (@f, ctl);
 438 %! assert (convergence == 1);
 439 %! assert (f(x) == obj_value);
 440 %! assert (obj_value < ctl.VTR);
 441
 442 %!demo
 443 %! ## define a simple example function
 444 %! f = @(x) peaks(x(1), x(2));
 445 %! ## plot the function to see where the minimum might be
 446 %! peaks()
 447 %! ## first we set the region where we expect the minimum
 448 %! ctl.XVmin = [-3 -3];
 449 %! ctl.XVmax = [ 3  3];
 450 %! ## and solve it with de_min
 451 %! [x, obj_value, nfeval, convergence] = de_min (f, ctl)