matlab的差分算法实现.m

一个计算Difference Algorithm的程序！

%   CR = 0.5; 
%   strategy = 7;
%   refresh = 10; 
%
% Cost function:   function result = f(x,y);
%   has to be defined by the user and is minimized
%      w.r. to x(1:D).
%
% Example to find the minimum of the Rosenbrock saddle:
% ----------------------------------------------------
% Define f.m as:
% function result = f(x,y);
% result = 100*(x(2)-x(1)^2)^2+(1-x(1))^2;
% end
% Then type:
%
%   VTR = 1.e-6;
%   D = 2; 
%   XVmin = [-2 -2]; 
%   XVmax = [2 2]; 
%   [bestmem,bestval,nfeval] = devec3("f",VTR,D,XVmin,XVmax);
%
% The same example with a more complete argument list is handled in 
% run1.m
%
% About devec3.m
% --------------
% Differential Evolution for MATLAB
% Copyright (C) 1996, 1997 R. Storn
% International Computer Science Institute (ICSI)
% 1947 Center Street, Suite 600
% Berkeley, CA 94704
% E-mail: [e-mail]
% WWW: http://http.icsi.berkeley.edu/~storn
%
% devec is a vectorized variant of DE which, however, has a
% propertiy which differs from the original version of DE:
% 1) The random selection of vectors is performed by shuffling the
% population array. Hence a certain vector can't be chosen twice
% in the same term of the perturbation expression.
%
% Due to the vectorized expressions devec3 executes fairly fast
% in MATLAB's interpreter environment.
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 1, or (at your option)
% any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details. A copy of the GNU 
% General Public License can be obtained from the 
% Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.

%-----Check input variables---------------------------------------------
err=[];
if nargin<1, error('devec3 1st argument must be function name'); else 
if exist(fname)<1; err(1,length(err)+1)=1; end; end;
if nargin<2, VTR = 1.e-6; else 
if length(VTR)~=1; err(1,length(err)+1)=2; end; end;
if nargin<3, D = 2; else
if length(D)~=1; err(1,length(err)+1)=3; end; end; 
if nargin<4, XVmin = [-2 -2];else
if length(XVmin)~=D; err(1,length(err)+1)=4; end; end; 
if nargin<5, XVmax = [2 2]; else
if length(XVmax)~=D; err(1,length(err)+1)=5; end; end; 
if nargin<6, y=[]; end; 
if nargin<7, NP = 10*D; else
if length(NP)~=1; err(1,length(err)+1)=7; end; end; 
if nargin<8, itermax = 200; else
if length(itermax)~=1; err(1,length(err)+1)=8; end; end; 
if nargin<9, F = 0.8; else
if length(F)~=1; err(1,length(err)+1)=9; end; end;
if nargin<10, CR = 0.5; else
if length(CR)~=1; err(1,length(err)+1)=10; end; end; 
if nargin<11, strategy = 7; else
if length(strategy)~=1; err(1,length(err)+1)=11; end; end;
if nargin<12, refresh = 10; else
if length(refresh)~=1; err(1,length(err)+1)=12; end; end; 
if length(err)>0
fprintf(stdout,'error in parameter %d\n', err);
usage('devec3 (string,scalar,scalar,vector,vector,any,integer,integer,scalar,scalar,integer,integer)');   
end

if (NP < 5)
NP=5;
fprintf(1,' NP increased to minimal value 5\n');
end
if ((CR < 0) | (CR > 1))
CR=0.5;
fprintf(1,'CR should be from interval [0,1]; set to default value 0.5\n');
end
if (itermax <= 0)
itermax = 200;
fprintf(1,'itermax should be > 0; set to default value 200\n');
end
refresh = floor(refresh);

%-----Initialize population and some arrays-------------------------------

pop = zeros(NP,D); %initialize pop to gain speed

%----pop is a matrix of size NPxD. It will be initialized-------------
%----with random values between the min and max values of the---------
%----parameters-------------------------------------------------------

for i=1:NP
pop(i,:) = XVmin + rand(1,D).*(XVmax - XVmin);
end

popold = zeros(size(pop)); % toggle population
val = zeros(1,NP); % create and reset the "cost array"
bestmem = zeros(1,D); % best population member ever
bestmemit = zeros(1,D); % best population member in iteration
nfeval = 0; % number of function evaluations

%------Evaluate the best member after initialization----------------------

ibest = 1; % start with first population member
val(1) = feval(fname,pop(ibest,:),y);
bestval = val(1); % best objective function value so far
nfeval = nfeval + 1;
for i=2:NP % check the remaining members
val(i) = feval(fname,pop(i,:),y);
nfeval = nfeval + 1;
if (val(i) < bestval) % if member is better
ibest = i; % save its location
bestval = val(i);
end 
end
bestmemit = pop(ibest,:); % best member of current iteration
bestvalit = bestval; % best value of current iteration

bestmem = bestmemit; % best member ever

%------DE-Minimization---------------------------------------------
%------popold is the population which has to compete. It is--------
%------static through one iteration. pop is the newly--------------
%------emerging population.----------------------------------------

pm1 = zeros(NP,D); % initialize population matrix 1
pm2 = zeros(NP,D); % initialize population matrix 2
pm3 = zeros(NP,D); % initialize population matrix 3
pm4 = zeros(NP,D); % initialize population matrix 4
pm5 = zeros(NP,D); % initialize population matrix 5
bm = zeros(NP,D); % initialize bestmember matrix
ui = zeros(NP,D); % intermediate population of perturbed vectors
mui = zeros(NP,D); % mask for intermediate population
mpo = zeros(NP,D); % mask for old population
rot = (0:1:NP-1); % rotating index array (size NP)
rotd= (0:1:D-1); % rotating index array (size D)
rt = zeros(NP); % another rotating index array
rtd = zeros(D); % rotating index array for exponential crossover
a1 = zeros(NP); % index array
a2 = zeros(NP); % index array
a3 = zeros(NP); % index array
a4 = zeros(NP); % index array
a5 = zeros(NP); % index array
ind = zeros(4);

iter = 1;
while ((iter < itermax) & (bestval > VTR))
popold = pop; % save the old population

ind = randperm(4); % index pointer array

a1 = randperm(NP); % shuffle locations of vectors
rt = rem(rot+ind(1),NP); % rotate indices by ind(1) positions
a2 = a1(rt+1); % rotate vector locations
rt = rem(rot+ind(2),NP);
a3 = a2(rt+1); 
rt = rem(rot+ind(3),NP);
a4 = a3(rt+1); 
rt = rem(rot+ind(4),NP);
a5 = a4(rt+1); 

pm1 = popold(a1,:); % shuffled population 1
pm2 = popold(a2,:); % shuffled population 2
pm3 = popold(a3,:); % shuffled population 3
pm4 = popold(a4,:); % shuffled population 4
pm5 = popold(a5,:); % shuffled population 5

for i=1:NP % population filled with the best member
bm(i,:) = bestmemit; % of the last iteration
end

mui = rand(NP,D) < CR; % all random numbers < CR are 1, 0 otherwise

if (strategy > 5)
st = strategy-5;     % binomial crossover
else
st = strategy;     % exponential crossover
mui=sort(mui');   % transpose, collect 1's in each column
for i=1:NP
n=floor(rand*D);
if n > 0
rtd = rem(rotd+n,D);
mui(:,i) = mui(rtd+1,i); %rotate column i by n
end
end
mui = mui';       % transpose back
end
mpo = mui < 0.5; % inverse mask to mui

if (st == 1) % DE/best/1
ui = bm + F*(pm1 - pm2); % differential variation
ui = popold.*mpo + ui.*mui; % crossover
elseif (st == 2) % DE/rand/1
ui = pm3 + F*(pm1 - pm2); % differential variation
ui = popold.*mpo + ui.*mui; % crossover
elseif (st == 3) % DE/rand-to-best/1
ui = popold + F*(bm-popold) + F*(pm1 - pm2); 
ui = popold.*mpo + ui.*mui; % crossover
elseif (st == 4) % DE/best/2
ui = bm + F*(pm1 - pm2 + pm3 - pm4); % differential variation
ui = popold.*mpo + ui.*mui; % crossover
elseif (st == 5) % DE/rand/2
ui = pm5 + F*(pm1 - pm2 + pm3 - pm4); % differential variation
ui = popold.*mpo + ui.*mui; % crossover
end

%-----Select which vectors are allowed to enter the new population------------
for i=1:NP
tempval = feval(fname,ui(i,:),y); % check cost of competitor
nfeval = nfeval + 1;
if (tempval <= val(i)) % if competitor is better than value in "cost array"
pop(i,:) = ui(i,:); % replace old vector with new one (for new iteration)
val(i) = tempval; % save value in "cost array"

%----we update bestval only in case of success to save time-----------
if (tempval < bestval) % if competitor better than the best one ever
bestval = tempval; % new best value
bestmem = ui(i,:); % new best parameter vector ever
end
end
end %---end for imember=1:NP

bestmemit = bestmem; % freeze the best member of this iteration for the coming 
% iteration. This is needed for some of the strategies.

%----Output section----------------------------------------------------------

if (refresh > 0)
if (rem(iter,refresh) == 0)
fprintf(1,'Iteration: %d, Best: %f, F: %f, CR: %f, NP: %d\n',iter,bestval,F,CR,NP);
for n=1:D
fprintf(1,'best(%d) = %f\n',n,bestmem(n));
end
end
end

iter = iter + 1;
end %---end while ((iter < itermax) ...


function result = rosen(x,y);
% Objective function for Differential Evolution
%
% Input Arguments: 
% ---------------
% x : parameter vector to be optimized
% y : data vector (remains fixed during optimization)
%
% Output Arguments:
% ----------------
% result : objective function value
%

%-----Rosenbrock's saddle------------------
result = 100*(x(2)-x(1)^2)^2+(1-x(1))^2;


% Initialization and run of differential evolution optimizer.
% A simpler version with fewer explicit parameters is in run0.m
%
% Here for Rosenbrock's function
% Change relevant entries to adapt to your personal applications
%
% The file ofunc.m must also be changed 
% to return the objective function
%

% VTR    "Value To Reach" (stop when ofunc < VTR)
    VTR = 1.e-6; 

% D    number of parameters of the objective function 
    D = 2; 

% XVmin,XVmax vector of lower and bounds of initial population
%     the algorithm seems to work well only if [XVmin,XVmax] 
%     covers the region where the global minimum is expected
% *** note: these are no bound constraints!! ***
    XVmin = [-2 -2]; 
    XVmax = [2 2];

% y    problem data vector (remains fixed during optimization)
    y=[]; 

% NP number of population members
    NP = 15; 

% itermax maximum number of iterations (generations)
    itermax = 200; 

% F DE-stepsize F ex [0, 2]
    F = 0.8; 

% CR crossover probabililty constant ex [0, 1]
    CR = 0.8; 

% strategy 1 --> DE/best/1/exp 6 --> DE/best/1/bin
% 2 --> DE/rand/1/exp 7 --> DE/rand/1/bin
% 3 --> DE/rand-to-best/1/exp 8 --> DE/rand-to-best/1/bin
% 4 --> DE/best/2/exp 9 --> DE/best/2/bin
% 5 --> DE/rand/2/exp else DE/rand/2/bin

    strategy = 7

% refresh intermediate output will be produced after "refresh"
% iterations. No intermediate output will be produced
% if refresh is < 1
    refresh = 10; 

[x,f,nf] = devec3('rosen',VTR,D,XVmin,XVmax,y,NP,itermax,F,CR,strategy,refresh)



% Initialization and run of differential evolution optimizer
% a more complex version with more explicit parameters is in run.m
%
% Here for Rosenbrock's function
% Change relevant entries to adapt to your personal applications
%
% The file ofunc.m must also be changed, 
% to return the objective function
%

% VTR    "Value To Reach" (stop when ofunc < VTR)
    VTR = 1.e-6; 

% D    number of parameters of the objective function 
    D = 2; 

% XVmin,XVmax vector of lower and upper bounds of initial population
%     the algorithm seems to work well only if [XVmin,XVmax] 
%     covers the region where the global minimum is expected
% *** note: these are no bound constraints!! ***
    XVmin = [-2 -2]; 
    XVmax = [2 2];