enda_main.m

由我收集或写出的GA源码

% Main program to run the ENDA_MOEA.

function ENDA_Main()
% Initialize the variables
sizeA = 100;   % The number of active populations
sizeP = 100;   % The number of passive populations
gen = 1000;       % The number of generations
M = 2;         % The number of objectives.       
N = 6;         % The number of decision variables.
V = 8;         % V is the sum of M and N.
xmin = 0;                       % Maximum of all the decision variables    
xmax = 1;                       % Minimum of all the decision variables

% Initialize the population
for i = 1 : sizeA
    % Initialize the decision variables
    for j = 1 : N
        hormonesA(i,j) = rand(1);                          % i.e f(i,j) = min + (max - min)*rand(1);
    end
end
for i = 1 : sizeP
    % Initialize the decision variables
    for j = 1 : N
        hormonesP(i,j) = rand(1);                         
    end
end

% Start the evolution process.The following are performed in each generation.Select the population using non-domination.Classy
% the hormones from A.Evaluate the crowding distance and preformance.Select parents.Geneate descendant.

for g= 1 : gen

% Function to evaluate the objective functions for the given input vector hormones.
    for j = 1 : sizeA
        hormonesA(j,N+1) = 1 - exp(-4*hormonesA(j,1))*(sin(6*pi*hormonesA(j,1)))^6;    % Objective function one
        sum = 0;
        for k = 2 : N
            sum = sum + hormonesA(j,k)/4;
        end
        gx = 1 + 9*(sum)^(0.25);                                                       % Intermediate function
        hormonesA(j,N+2) = gx*(1 - ((hormonesA(j,N+1))/(gx))^2);                       % Objective function anoter one
    end
    for j = 1 : sizeP
        hormonesP(j,N+1) = 1 - exp(-4*hormonesP(j,1))*(sin(6*pi*hormonesP(j,1)))^6; 
        sum = 0;
        for k = 2 : N
            sum = sum + hormonesP(j,k)/4;
        end
        gx = 1 + 9*(sum)^(0.25);                                      
        hormonesP(j,N+2) = gx*(1 - ((hormonesP(j,N+1))/(gx))^2); 
    end
    
% This function sort the current popultion based on non-domination. 
    hormonesU = [hormonesA;hormonesP];
    num = [];                                         % Number of individuals that dominate this individual
    p = [];                                           % hormonesP - non-dominated solution from U
    [sizeU,K] = size(hormonesU);
    for i = 1 : sizeU
        num(i) = 0;
        for j = 1 : sizeU
            dom_less = 0;dom_equal = 0;dom_more = 0;
            for k = 1 : M
                if (hormonesU(i,N + k) < hormonesU(j,N + k))          % p=u(i,N + k),q=u(j,N + k))
                    dom_less = dom_less + 1;             
                elseif (hormonesU(i,N + k) == hormonesU(j,N + k))
                    dom_equal = dom_equal + 1;
                else
                    dom_more = dom_more + 1;
                end
            end
           if dom_less == 0 & dom_equal ~= M
               num(i) = num(i) + 1;                                   % the number of q dominate p
           end
        end 
        if num(i) == 0                                                % the number of that dominate p
            p = [p  i];                                               % add u to the set P. i.e. P = P∪{u}
        end
    end
    hormonesP = hormonesU(p,:);
    [sizeP,K] = size(hormonesP);
        
% Classify the hormones from A according to P .Then calculate the value of p's class and performance value of a hormone a.
    [sizeP,K] = size(hormonesP);
    for  j = 1 : M
        minP = min(hormonesP(:,N + j));
        maxP = max(hormonesP(:,N + j));
        minA = min(hormonesA(:,N + j));
        maxA = max(hormonesA(:,N + j));
        fmin = min(abs(minP),abs(minA));
        fmax = max(abs(maxP),abs(maxA));
        if (fmax == fmin)
            if (fmax == 0)
                fun(j) = 1;
            else
                fun(j) = abs(fmax);
            end
        else
            fun(j) = max(abs(fmin),abs(fmax));
        end
        hormonesA(:,N+j) = hormonesA(:,N+j)/fun(j);
    end 
    for j = 1 : sizeA
        t = 1; 
        a = abs(hormonesA(j,N+1)-hormonesP(1,N+1)/fun(1))+abs(hormonesA(j,N+1)-hormonesP(1,N+2)/fun(1));
        for k = 2 : sizeP  
            b = abs(hormonesA(j,N+1)-hormonesP(k,N+1)/fun(2))+abs(hormonesA(j,N+1)-hormonesP(k,N+2)/fun(2));
            if b<=a
                a = b;t = k;
            end
        end
        hormonesA(j,V+1) = t;
    end
% Calculate the value of a's class.
    for j = 1 : sizeP
        num = 0;
        for k = 1 : sizeA
            if hormonesA(k,V+1) == j
                num = num + 1;
            end
        end
        if num == 0
            val = 0;
        else
            val =1 / num;                            % Calculate the value of a's class.
        end
        for k = 1 : sizeA
            if hormonesA(k,V+1) == j
                hormonesA(k,V+2) = val;
            end
        end
    end    
    % Calculate the performance value of a hormone.
    ndnum = [];                                      % Number of solutions that dominated by the hormone a.  
    for a = 1 : sizeA
        ndnum = 0;
        for j = 1 : sizeA
            dom_less = 0;
            dom_equal = 0;
            dom_more = 0;
            for k = 1 : M
                if (hormonesA(a,N + k) < hormonesA(j,N + k))
                    dom_less = dom_less + 1;
                elseif (hormonesA(a,N + k) == hormonesA(j,N + k))
                    dom_equal = dom_equal + 1;
                else
                    dom_more = dom_more + 1;
                end
            end
            if dom_more == 0 & dom_equal ~= M
                ndnum = ndnum + 1;              % the number of p dominate q                
            end
            hormonesA(:,V+3) = ndnum/sizeA;       % the performance value of a hormone a.
        end 
    end
    
    
 % This function select the first parent of the descendant. First parent is selected by using the tournament selection
 % operator and the second is selected among the neighborhood of the first one, by using the proportional selection.
    firstsize = sizeA;
    for j = 1 : firstsize
        for k = 1 : 2
            a(k) = round(sizeA*rand(1));
            if a(k) == 0
                a(k) = 1;
            end
            if k > 1
                while ~isempty(find(a(1 : k - 1) == a(k)))
                    a(k) = round(sizeA*rand(1));
                    if a(k) == 0
                        a(k) = 1;
                    end
                end
            end
        end
        if hormonesA(a(1),V+2) > hormonesA(a(2),V+2)
            y = hormonesA(a(1),:);
        elseif hormonesA(a(1),V+2) == hormonesA(a(2),V+2)
            if hormonesA(a(1),N+1) + hormonesA(a(1),N+2) <= hormonesA(a(2),N+1) + hormonesA(a(2),N+2)
                y = hormonesA(a(1),:);
            else
                y = hormonesA(a(2),:);
            end
        else
            y = hormonesA(a(2),:);
        end
       first(j,:) =y;                                              % This function generate the first parent.
    end
    for j = 1 : firstsize
        second = [];y = [];
        x = first(j,:);
        for k = 1 : sizeA
            if hormonesA(k,V + 1) == x(V + 1)
                second = [second;hormonesA(k,:)];
            end
        end
        [sizesecond,K] = size(second);
        if sizesecond == 1
            y = rand(1,V+3); 
        else
            a = find(second(:,V+3) == max(second(:,V+3)));
            [sizea,K] = size(a);
            l = round(sizea*rand(1));
            if l == 0
                l = 1;
            end
            y = second(a(l),:);                         % This function generate the second parent.
        end
        p = rand(1);q = rand(1);
        if p > 0.5
            z(1) = (1 - p)*x(1)+p*y(1);
        else
            z(1) = rand(1);
        end    
        for k = 2 : N
            z(k) = (1 - q)*x(k)+q*y(k);
        end
        for k = 1 : N
            if z(k) > xmax
                z(k) = xmax;
            elseif z(k) < xmin
                z(k) = xmin;
            end
        end
        descendant(j,:) = z;
    end
    hormonesA = descendant;    
    gen = gen + 1;
end

% Visualize.The following is used to visualize the result for the given problem.
plot(hormonesP(:,N + 1),hormonesP(:,N + 2),'*');
title('MOP using ENDA');
xlabel('f(x_1)');
ylabel('f(x_2)');