MOEAD算法程序源代碼

1、eval_update     function idealpoint, subproblems= eval_update(idealpoint, subproblems, inds)%EvaluationUpdate Summary of this function goes here%   Detailed explanation goes here objs = inds.objective;weights = subproblems.weight; idealpoint = min(id

2、ealpoint, min(objs,2);for i=1:length(inds)    subobjs = subobjective(weights, objs(:,i), idealpoint, 'ws');     %update the values.    C = subobjs<subproblems.optimal;    if any(C)    

3、    ncell = num2cell(subobjs(C);        subproblems(C).optimal = ncell:;           subproblems(C).optpoint = deal(inds(i);    endendend evaluatefunction v, x =

4、evaluate( prob, x )%EVALUATE function evaluate an individual structure of a vector point with%the given multiobjective problem. %   Detailed explanation goes here%   prob: is the multiobjective problem.%   x: is a vector point, or a individual struct

5、ure.%   v: is the result objectives evaluated by the mop.%   x: if x is a individual structure, then x's objective field is modified%   with the evaluated value and pass back. %   TODO, need to refine it to operate on a vector of p

6、oints.    if isstruct(x)        v = prob.func(x.parameter);        x.objective=v;    else        v = prob.func(x);end gaussian_mutatefu

7、nction ind = gaussian_mutate( ind, prob, domain)%GAUSSIAN_MUTATE Summary of this function goes here%   Detailed explanation goes here if isstruct(ind)    x = ind.parameter;else    x  = ind;end    parDim = lengt

8、h(x);   lowend  = domain(:,1);   highend =domain(:,2);   sigma = (highend-lowend)./20;     newparam = min(max(normrnd(x, sigma), lowend), highend);   C = rand(parDim, 1)<prob;   x(C) = ne

9、wparam(C);  if isstruct(ind)    ind.parameter = x;else    ind = x;end genetic_op  function ind=genetic_op(subproblems, index, domain, params)%GENETICOP function implemented the DE operation to generate a new%individual from a subproble

10、ms and its neighbours. %   subproblems: is all the subproblems.%   index: the index of the subproblem need to handle.%   domain: the domain of the origional multiobjective problem.%   ind: is an individual structure.   &

11、#160;neighbourindex = subproblems(index).neighbour;       %The random draw from the neighbours.    nsize = length(neighbourindex);    si = ones(1,3)*index;       si(1)=neighbourindex(ceil(ra

12、nd*nsize);    while si(1)=index        si(1)=neighbourindex(ceil(rand*nsize);    end       si(2)=neighbourindex(ceil(rand*nsize);    while si(2)=index | si(2)=si(1)&

13、#160;       si(2)=neighbourindex(ceil(rand*nsize);    end       si(3)=neighbourindex(ceil(rand*nsize);    while si(3)=index | si(3)=si(2) | si(3)=si(1)      

14、60; si(3)=neighbourindex(ceil(rand*nsize);    end        %retrieve the individuals.    points = subproblems(si).curpoint;    selectpoints = points.parameter;      &#

15、160;oldpoint = subproblems(index).curpoint.parameter;    parDim = size(domain, 1);       jrandom = ceil(rand*parDim);       randomarray = rand(parDim, 1);    deselect = randomarray<params

16、.CR;    deselect(jrandom)=true;    newpoint = selectpoints(:,1)+params.F*(selectpoints(:,2)-selectpoints(:,3);    newpoint(deselect)=oldpoint(deselect);       %repair the new value.    

17、newpoint=max(newpoint, domain(:,1);    newpoint=min(newpoint, domain(:,2);       ind = struct('parameter',newpoint,'objective', 'estimation',);    %ind.parameter = newpoint;    %ind

18、 = realmutate(ind, domain, 1/parDim);    ind = gaussian_mutate(ind, 1/parDim, domain);       %clear points selectpoints oldpoint randomarray deselect newpoint neighbourindex si;end  get_structurefunction str = get_structure( name )%STR

19、UCTURE Summary of this function goes here% Structure used in this toolbox.% individual structure:% parameter: the parameter space point of the individual. it's a column-wise% vector.% objective: the objective space point of the individual. it's column-wise% vector. It only have value after e

20、valuate function is called upon the% individual.% estimation: Also a structure array of the individual. It's not used in% MOEA/D but used in MOEA/D/GP. For every objective, the field contains the% estimation from the GP model.% estimation structure:% obj: the estimated mean.% std: the estimated

21、standard deviation for the mean.% subproblem structure:% weight: the decomposition weight for the subproblem.% optimal: the current optimal value of the current structure.% curpoiont: the current individual of the subproblem.% optpoint: the point that gain the optimal on the subproblem.% switch

22、 name    case 'individual'        str = struct('parameter','objective','estimation');    case 'subproblem'        str = struct('weight&#

23、39;,'optimal','curpoint','optpoint',);    case 'estimation'        str = struct();                    otherw

24、ise       end  init_weightsfunction subp=init_weights(popsize, niche, objDim)% init_weights function initialize a pupulation of subproblems structure% with the generated decomposition weight and the neighbourhood% relationship.    subp

25、=;    for i=0:popsize        if objDim=2            p=struct('weight','neighbour','optimal', Inf, 'optpoint', 'curpoint', );  

26、0;         weight=zeros(2,1);            weight(1)=i/popsize;            weight(2)=(popsize-i)/popsize;     &#

27、160;      p.weight=weight;            subp=subp p;        elseif objDim=3        %TODO       &#

28、160;end    end % weight = lhsdesign(popsize, objDim, 'criterion','maximin', 'iterations', 1000)'% p=struct('weight','neighbour','optimal', Inf, 'optpoint', 'curpoint', );% subp = repmat(p, popsize, 1);% cell

29、s = num2cell(weight);% subp.weight=cells:;     %Set up the neighbourhood.    leng=length(subp);    distanceMatrix=zeros(leng, leng);    for i=1:leng        for j=i+1:leng 

30、;           A=subp(i).weight;B=subp(j).weight;            distanceMatrix(i,j)=(A-B)'*(A-B);            distanceMatrix(j,i

31、)=distanceMatrix(i,j);        end        s,sindex=sort(distanceMatrix(i,:);        subp(i).neighbour=sindex(1:niche)'    end  end   

32、;      moeadfunction pareto = moead( mop, varargin)%MOEAD runs moea/d algorithms for the given mop.%   Detailed explanation goes here%   The mop must to be minimizing.%   The parameters of the algorithms can be set through va

33、rargin. including%   popsize: The subproblem's size.%   niche: the neighboursize, must less then the popsize.%   iteration: the total iteration of the moead algorithms before finish.%   method: the decomposition method, the value can be

34、 'ws' or 'ts'.     starttime = clock;    %global variable definition.    global params idealpoint objDim parDim itrCounter;    %set the random generator.    rand('state',1

35、0);       %Set the algorithms parameters.    paramIn = varargin;    objDim, parDim, idealpoint, params, subproblems=init(mop, paramIn);       itrCounter=1;    while termi

36、nate(itrCounter)        tic;        subproblems = evolve(subproblems, mop, params);        disp(sprintf('iteration %u finished, time used: %u', itrCounter, toc);  

37、;      itrCounter=itrCounter+1;    end       %display the result.    pareto=subproblems.curpoint;    pp=pareto.objective;    scatter(pp(1,:), pp(2,:); 

38、;   disp(sprintf('total time used %u', etime(clock, starttime);end function objDim, parDim, idealp, params, subproblems=init(mop, propertyArgIn)%Set up the initial setting for the MOEA/D.    objDim=mop.od;    parDim=mop.pd; 

39、60;     idealp=ones(objDim,1)*inf;       %the default values for the parameters.    params.popsize=100;params.niche=30;params.iteration=100;    params.dmethod='ts'    param

40、s.F = 0.5;    params.CR = 0.5;       %handle the parameters, mainly about the popsize    while length(propertyArgIn)>=2        prop = propertyArgIn1;     

41、60;  val=propertyArgIn2;        propertyArgIn=propertyArgIn(3:end);         switch prop            case 'popsize'    

42、            params.popsize=val;            case 'niche'                params.niche=val;

43、0;           case 'iteration'                params.iteration=val;            case 'method

44、9;                params.dmethod=val;            otherwise                warnin

45、g('moea doesnot support the given parameters name');        end    end       subproblems = init_weights(params.popsize, params.niche, objDim);    params.popsize = length(subprob

46、lems);       %initial the subproblem's initital state.    inds = randompoint(mop, params.popsize);    V, INDS = arrayfun(evaluate, repmat(mop, size(inds), inds, 'UniformOutput', 0);    v = cell

47、2mat(V);    idealp = min(idealp, min(v,2);        %indcells = mat2cell(INDS, 1, ones(1,params.popsize);    subproblems.curpoint = INDS:;    clear inds INDS V indcells;end   function sub

48、problems = evolve(subproblems, mop, params)    global idealpoint;      for i=1:length(subproblems)        %new point generation using genetic operations, and evaluate it.      

49、0; ind = genetic_op(subproblems, i, mop.domain, params);        obj,ind = evaluate(mop, ind);        %update the idealpoint.        idealpoint = min(idealpoint, obj); 

50、60;              %update the neighbours.        neighbourindex = subproblems(i).neighbour;        subproblems(neighbourindex)=update(subproblems(neigh

51、bourindex),ind, idealpoint);        %clear ind obj neighbourindex neighbours;                clear ind obj neighbourindex;    endend function subp =update(

52、subp, ind, idealpoint)    newobj=subobjective(subp.weight, ind.objective,  idealpoint, 'te');    oops = subp.curpoint;    oldobj=subobjective(subp.weight, oops.objective, idealpoint, 'te' );    

53、    C = newobj < oldobj;    subp(C).curpoint= deal(ind);    clear C newobj oops oldobj;end function y =terminate(itrcounter)    global params;    y = itrcounter>params.iteration;end &#

54、160; randompointfunction ind = randompoint(prob, n)%RANDOMNEW to generate n new point randomly from the mop problem given. if (nargin=1)    n=1;end randarray = rand(prob.pd, n);lowend = prob.domain(:,1);span = prob.domain(:,2)-lowend;point = randarray.*(span(:,ones

55、(1, n)+ lowend(:,ones(1,n);cellpoints = num2cell(point, 1); indiv = struct('parameter','objective', 'estimation', );ind = repmat(indiv, 1, n);ind.parameter = cellpoints:; % estimation = struct('obj', NaN ,'std', NaN);% ind.estimation = deal(repm

56、at(estimation, prob.od, 1);end  realmutatefunction ind = realmutate(ind, domains, rate)%REALMUTATE Summary of this function goes here%   Detailed explanation goes here    % double rnd, delta1, delta2, mut_pow, deltaq;   % double y, yl, yu,

57、val, xy;   % double eta_m = id_mu;    eta_m=20;   numVariables = size(domains,1);   if (isstruct(ind)       a = ind.parameter;   else       a = ind;&#

58、160;  end   for j = 1:numVariables        if (rand() <= rate)            y = a(j);               

59、        yl = domains(j,1);            yu = domains(j,2);            delta1 = (y - yl) / (yu - yl);       

60、     delta2 = (yu - y) / (yu - yl);             rnd = rand();            mut_pow = 1.0 / (eta_m + 1.0);        

61、;    if (rnd <= 0.5)                xy = 1.0 - delta1;                val = 2.0 * rnd + (1.0 - 2.0 * rnd) * (xy(eta_m

62、 + 1.0);                deltaq = (valmut_pow) - 1.0;            else               &#

63、160;xy = 1.0 - delta2;                val = 2.0 * (1.0 - rnd) + 2.0 * (rnd - 0.5) * (xy (eta_m + 1.0);                deltaq = 1.0 - (valmut

64、_pow);            end                       y = y + deltaq * (yu - yl);         

65、   if (y < yl)                y = yl;            end            if (y > yu)

66、0;               y = yu;            end            a(j) = y;       

67、60;       end   end   if isstruct(ind)       ind.parameter = a;   else       ind = a;   endend      subob

68、jectivefunction obj = subobjective(weight, ind, idealpoint, method)%SUBOBJECTIVE function evaluate a point's objective with a given method of%decomposition. %   Two method are implemented by far is Weighted-Sum and Tchebesheff.%   weight: is the decomposition w

69、eight.(column wise vector).%   ind: is the individual point(column wise vector).%   idealpoint: the idealpoint for Tchebesheff decomposition.%   method: is the decomposition method, the default is 'te' when is%   omitted.% 

70、0;%   weight and ind can also be matrix. in which have two scenairos:%   When weight is a matrix, then it's treated as a column wise set of%   weights. in that case, if ind is a size 1 column vector, then the%   subobjective is computed

71、 with every weight and the ind; if ind is also%   a matrix of the same size as weight, then the subobjective is computed%   in a column-to-column, with each column of weight computed against the%   corresponding column of ind.%   A row vect

72、or of subobjective is return in both case.     if (nargin=2)        obj = ws(weight, ind);    elseif (nargin=3)        obj = te(weight, ind, idealpoint);   

73、0;else        if strcmp(method, 'ws')            obj=ws(weight, ind);        elseif strcmp(method, 'te')      

74、0;     obj=te(weight, ind, idealpoint);        else            obj= te(weight, ind, idealpoint);        end    ende

75、nd function obj = ws(weight, ind)    obj = (weight'*ind)'end function obj = te(weight, ind, idealpoint)    s = size(weight, 2);    indsize = size(ind,2);       weight(weight = 0)=0.00001;

76、       if indsize=s        part2 = abs(ind-idealpoint(:,ones(1, indsize);        obj = max(weight.*part2);    elseif indsize =1        part2 = abs(ind-idealpoint);        obj = max(weight.*part2(:,ones(1, s);      else        error('individual size must be same as we