MOEAD算法程序源代码

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(idealpoint, 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)

ncell = num2cell(subobjs(C));

[subproblems(C).optimal] = ncell{:}; [subproblems(C).optpoint] = deal(inds(i)); end end end

evaluate

function [v, x] = 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 structure. % 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 points. if isstruct(x)

v = prob.func(x.parameter); x.objective=v; else

v = prob.func(x);

end

gaussian_mutate

function 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 = length(x); lowend = domain(:,1); highend =domain(:,2); sigma = (highend-lowend)./20;

newparam = min(max(normrnd(x, sigma), lowend), highend); C = rand(parDim, 1)

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 subproblems 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.

neighbourindex = subproblems(index).neighbour;

%The random draw from the neighbours. nsize = length(neighbourindex); si = ones(1,3)*index;

si(1)=neighbourindex(ceil(rand*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)

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) si(3)=neighbourindex(ceil(rand*nsize)); end

%retrieve the individuals.

points = [subproblems(si).curpoint]; selectpoints = [points.parameter];

oldpoint = subproblems(index).curpoint.parameter; parDim = size(domain, 1);

jrandom = ceil(rand*parDim);

randomarray = rand(parDim, 1); deselect = randomarray

deselect(jrandom)=true;

newpoint = selectpoints(:,1)+params.F*(selectpoints(:,2)-selectpoints(:,3)); newpoint(~deselect)=oldpoint(~deselect);

%repair the new value.

newpoint=max(newpoint, domain(:,1)); newpoint=min(newpoint, domain(:,2));

ind = struct('parameter',newpoint,'objective',[], 'estimation',[]); %ind.parameter = newpoint;

%ind = realmutate(ind, domain, 1/parDim); ind = gaussian_mutate(ind, 1/parDim, domain);

%clear points selectpoints oldpoint randomarray deselect newpoint neighbourindex si;

end

get_structure

function str = get_structure( name )

%STRUCTURE 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 evaluate 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 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 name case 'individual'

str = struct('parameter',[],'objective'[],'estimation'[]); case 'subproblem'

str = struct('weight',[],'optimal',[],'curpoint',[],'optpoint',[]); case 'estimation'

str = struct(); otherwise end

init_weights

function 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=[]; for i=0:popsize if objDim==2

p=struct('weight',[],'neighbour',[],'optimal', Inf, 'optpoint',[], 'curpoint', []); weight=zeros(2,1); weight(1)=i/popsize;

weight(2)=(popsize-i)/popsize; p.weight=weight; subp=[subp p];

elseif objDim==3 %TODO end end

% weight = lhsdesign(popsize, objDim, 'criterion','maximin', 'iterations', 1000)'; % p=struct('weight',[],'neighbour',[],'optimal', Inf, 'optpoint',[], 'curpoint', []); % subp = repmat(p, popsize, 1); % cells = 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

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

[s,sindex]=sort(distanceMatrix(i,:)); subp(i).neighbour=sindex(1:niche)'; end end

moead

function 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 varargin. 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 'ws' or 'ts'.

starttime = clock;

%global variable definition.

global params idealpoint objDim parDim itrCounter; %set the random generator. rand('state',10);

%Set the algorithms parameters. paramIn = varargin;

[objDim, parDim, idealpoint, params, subproblems]=init(mop, paramIn);

itrCounter=1;

while ~terminate(itrCounter) tic;

subproblems = evolve(subproblems, mop, params);

disp(sprintf('iteration %u finished, time used: %u', itrCounter, toc)); itrCounter=itrCounter+1; end

%display the result.

pareto=[subproblems.curpoint]; pp=[pareto.objective]; scatter(pp(1,:), pp(2,:));

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;

idealp=ones(objDim,1)*inf;

%the default values for the parameters.

params.popsize=100;params.niche=30;params.iteration=100; params.dmethod='ts'; params.F = 0.5; params.CR = 0.5;

%handle the parameters, mainly about the popsize while length(propertyArgIn)>=2 prop = propertyArgIn{1}; val=propertyArgIn{2};

propertyArgIn=propertyArgIn(3:end);

switch prop case 'popsize'

params.popsize=val; case 'niche'

params.niche=val; case 'iteration'

params.iteration=val; case 'method'

params.dmethod=val; otherwise

warning('moea doesnot support the given parameters name'); end end

subproblems = init_weights(params.popsize, params.niche, objDim); params.popsize = length(subproblems);

%initial the subproblem's initital state. inds = randompoint(mop, params.popsize);

[V, INDS] = arrayfun(@evaluate, repmat(mop, size(inds)), inds, 'UniformOutput', 0);

v = cell2mat(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 subproblems = evolve(subproblems, mop, params) global idealpoint;

for i=1:length(subproblems)

%new point generation using genetic operations, and evaluate it. ind = genetic_op(subproblems, i, mop.domain, params); [obj,ind] = evaluate(mop, ind); %update the idealpoint.

idealpoint = min(idealpoint, obj);

%update the neighbours.

neighbourindex = subproblems(i).neighbour;

subproblems(neighbourindex)=update(subproblems(neighbourindex),ind, idealpoint);

%clear ind obj neighbourindex neighbours;

clear ind obj neighbourindex; end end

function subp =update(subp, ind, idealpoint)

newobj=subobjective([subp.weight], ind.objective, idealpoint, 'te'); oops = [subp.curpoint];

oldobj=subobjective([subp.weight], [oops.objective], idealpoint, 'te' );

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

randompoint

function 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(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(repmat(estimation, prob.od, 1)); end

realmutate

function 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, val, xy; % double eta_m = id_mu;

eta_m=20;

numVariables = size(domains,1); if (isstruct(ind)) a = ind.parameter; else a = ind; end

for j = 1:numVariables if (rand() <= rate) y = a(j);

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

rnd = rand();

mut_pow = 1.0 / (eta_m + 1.0); if (rnd <= 0.5) xy = 1.0 - delta1;

val = 2.0 * rnd + (1.0 - 2.0 * rnd) * (xy^(eta_m + 1.0)); deltaq = (val^mut_pow) - 1.0; else

xy = 1.0 - delta2;

val = 2.0 * (1.0 - rnd) + 2.0 * (rnd - 0.5) * (xy^ (eta_m + 1.0)); deltaq = 1.0 - (val^mut_pow); end

y = y + deltaq * (yu - yl);

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

a(j) = y; end end

if isstruct(ind) ind.parameter = a; else ind = a; end end

subobjective

function obj = subobjective(weight, ind, idealpoint, method)

%SUBOBJECTIVE function evaluate a point's objective with a given method of Tcomposition.

% Two method are implemented by far is Weighted-Sum and Tchebesheff. % weight: is the decomposition weight.(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. %

% 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 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 vector of subobjective is return in both case.

if (nargin==2)

obj = ws(weight, ind); elseif (nargin==3)

obj = te(weight, ind, idealpoint); else

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

obj= te(weight, ind, idealpoint); end end end

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;

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 weight size, or equals 1'); end end

testmop

function mop = testmop( testname, dimension ) %Get test multi-objective problems from a given name.

% The method get testing or benchmark problems for Multi-Objective % Optimization. The implemented problems included ZDT, OKA, KNO. % User get the corresponding test problem, which is an instance of class % mop, by passing the problem name and optional dimension parameters.

mop=struct('name',[],'od',[],'pd',[],'domain',[],'func',[]); switch lower(testname) case 'kno1'

mop=kno1(mop); case 'zdt1'

mop=zdt1(mop, dimension); otherwise

error('Undefined test problem name'); end end

%KNO1 function generator function p=kno1(p) p.name='KNO1'; p.od = 2; p.pd = 2;

p.domain= [0 3;0 3]; p.func = @evaluate;

%KNO1 evaluation function. function y = evaluate(x) y=zeros(2,1);

c = x(1)+x(2);

f = 9-(3*sin(2.5*c^0.5) + 3*sin(4*c) + 5 *sin(2*c+2)); g = (pi/2.0)*(x(1)-x(2)+3.0)/6.0; y(1)= 20-(f*cos(g)); y(2)= 20-(f*sin(g)); end end

%ZDT1 function generator function p=zdt1(p,dim) p.name='ZDT1'; p.pd=dim; p.od=2;

p.domain=[zeros(dim,1) ones(dim,1)]; p.func=@evaluate;

%KNO1 evaluation function. function y=evaluate(x) y=zeros(2,1); y(1) = x(1);

su = sum(x)-x(1); g = 1 + 9 * su / (dim - 1); y(2) =g*(1 - sqrt(x(1) / g)); end end

c = x(1)+x(2);

f = 9-(3*sin(2.5*c^0.5) + 3*sin(4*c) + 5 *sin(2*c+2)); g = (pi/2.0)*(x(1)-x(2)+3.0)/6.0; y(1)= 20-(f*cos(g)); y(2)= 20-(f*sin(g)); end end

%ZDT1 function generator function p=zdt1(p,dim) p.name='ZDT1'; p.pd=dim; p.od=2;

p.domain=[zeros(dim,1) ones(dim,1)]; p.func=@evaluate;

%KNO1 evaluation function. function y=evaluate(x) y=zeros(2,1); y(1) = x(1);

su = sum(x)-x(1); g = 1 + 9 * su / (dim - 1); y(2) =g*(1 - sqrt(x(1) / g)); end end

本文来源：https://www.bwwdw.com/article/npir.html