1 ## Copyright (C) 1996, 1997 R. Storn
2 ## Copyright (C) 2009-2010 Christian Fischer <cfischer@itm.uni-stuttgart.de>
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
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
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/>.
17 ## de_min: global optimisation using differential evolution
19 ## Usage: [x, obj_value, nfeval, convergence] = de_min(fcn, control)
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
29 ## fcn string : Name of function. Must return a real value
30 ## control vector : (Optional) Control variables, described below
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)
44 ## Control variable: (optional) may be named arguments (i.e. "name",value
45 ## ---------------- pairs), a struct, or a vector, where
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
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)
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.
108 ## Example to find the minimum of the Rosenbrock saddle:
109 ## ----------------------------------------------------
111 ## function result = f(x);
112 ## result = 100 * (x(2) - x(1)^2)^2 + (1 - x(1))^2;
116 ## ctl.XVmin = [-2 -2];
117 ## ctl.XVmax = [ 2 2];
118 ## [x, obj_value, nfeval, convergence] = de_min (@f, ctl);
120 ## Keywords: global-optimisation optimisation minimisation
122 function [bestmem, bestval, nfeval, convergence] = de_min(fcn, varargin)
129 NP = 0; # NP will be set later
139 ## ------------ Check input variables (ctl) --------------------------------
140 if nargin >= 2, # Read control arguments
143 ctl = struct (varargin{:});
145 ctl = varargin{va_arg_cnt++};
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
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)
178 if isfield (ctl, "maxiter") && ! isnan (ctl.maxiter)
179 maxiter = ctl.maxiter;
184 ## set dimension D and population size NP
190 ## -------- do a few sanity checks --------------------------------
191 if (length (XVmin) != length (XVmax))
192 error("Length of upper and lower bounds does not match.")
195 error("Population size NP must be bigger than 5.")
197 if ((F <= 0) || (F > 2))
198 error("Difference Factor F out of range (0,2].")
200 if ((CR < 0) || (CR > 1))
201 error("CR value out of range [0,1].")
204 error("maxiter must be positive.")
207 error("maxnfe must be positive.")
209 refresh = floor(abs(refresh));
211 ## ----- Initialize population and some arrays --------------------------
213 pop = zeros(NP,D); # initialize pop
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
218 pop(i,:) = XVmin + rand (1,D) .* (XVmax - XVmin);
221 ## initialise the weighting factors between 0.0 and 1.0
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
231 ## ------ Evaluate the best member after initialization ------------------
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
245 nfeval = nfeval + NP;
246 bestmemit = pop(ibest,:); # best member of current iteration
247 bestvalit = bestval; # best value of current iteration
249 bestmem = bestmemit; # best member ever
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.
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)
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
267 ind = randperm (4); # index pointer array
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);
274 rt = rem (rot +ind(3), NP);
276 rt = rem (rot + ind(4), NP);
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
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
289 mui = rand (NP, D) < CR; # mask for intermediate population
290 # all random numbers < CR are 1, 0 otherwise
293 st = strategy - 6; # binomial crossover
295 st = strategy; # exponential crossover
296 mui = sort (mui'); # transpose, collect 1's in each column
298 n = floor (rand * D);
300 rtd = rem (rotd + n, D);
301 mui(:,i) = mui(rtd+1,i); #rotate column i by n
304 mui = mui'; # transpose back
306 mpo = mui < 0.5; # inverse mask to mui
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
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
329 neigh_ind = i-nr:i+nr; # index range of neighbourhood
330 neigh_ind = neigh_ind + ((neigh_ind <= 0)-(neigh_ind > NP))*NP;
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
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
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
353 ui = popold.*mpo + ui.*mui;
355 ## enforce initial bounds/constraints if specified
358 ui(i,:) = max (ui(i,:), XVmin);
359 ui(i,:) = min (ui(i,:), XVmax);
363 ## ----- Select which vectors are allowed to enter the new population ------
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
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
378 endfor #---end for i = 1:NP
380 nfeval = nfeval + NP; # increase number of function evaluations
382 bestmemit = bestmem; # freeze the best member of this iteration for the
383 # coming iteration. This is needed for some of the
386 ## ---- Output section ----------------------------------------------------
389 if (rem (iter, refresh) == 0)
390 printf ('Iteration: %d, Best: %8.4e, Worst: %8.4e\n', ...
391 iter, bestval, max(val));
393 printf ('x(%d) = %e\n', n, bestmem(n));
399 endwhile #---end while ((iter < maxiter) ...
401 ## check that all variables are well within bounds/constraints
404 range = XVmax - XVmin;
405 if (ui(i,:) < XVmin + 0.01*range)
408 if (ui(i,:) > XVmax - 0.01*range)
413 ## create the convergence result
416 elseif (abs (max (val) - bestval) / max (1, abs (max (val))) <= tol)
418 elseif (boundsOK == 0)
420 elseif (iter >= maxiter)
422 elseif (nfeval >= maxnfe)
428 %!function result = f(x);
429 %! result = 100 * (x(2) - x(1)^2)^2 + (1 - x(1))^2;
434 %! ctl.XVmin = [-2 -2];
435 %! ctl.XVmax = [ 2 2];
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);
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
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)