USER
can you plot all of pareto plot and display the metrics for every run:
% MOEMA_PerformanceMetrics_Full.m
% Main script to run MOEMA optimization and evaluate performance metrics
clc;
clear;
close all;
%% Parameters
NPop = 50; % Population size
T = 50; % Number of iterations
dim = 10; % Number of decision variables
lb = [-5*ones(1,dim-1), 0]; % Lower bounds: [-5,...,-5, 0]
ub = [5*ones(1,dim-1), 1]; % Upper bounds: [5,...,5, 1]
num_runs = 10; % Number of runs for statistical evaluation
% Initialize metrics storage
GD_values = zeros(num_runs, 1);
IGD_values = zeros(num_runs, 1);
Spacing_values = zeros(num_runs, 1);
Spread_values = zeros(num_runs, 1);
Hypervolume_values = zeros(num_runs, 1);
% Reference point for hypervolume calculation
reference_point = [1.1, 1.1]; % Adjust based on your problem
%% Run MOEMA Multiple Times
for run = 1:num_runs
fprintf('Run #%d/%d\n', run, num_runs);
% Run MOEMA
[paretoFront, paretoSet, Trajectories, fitness_history, position_history] = MOEMA(NPop, T, lb, ub, dim, @zdt4);
% Get True Pareto Front
true_pareto = get_true_pareto();
% Calculate all metrics
[GD_values(run), IGD_values(run), Spacing_values(run), Spread_values(run)] = ...
calculate_metrics(paretoFront, true_pareto);
% Calculate hypervolume
Hypervolume_values(run) = calculate_hypervolume(paretoFront, reference_point);
end
%% Calculate Statistics
metrics = {GD_values, IGD_values, Spacing_values, Spread_values, Hypervolume_values};
metric_names = {'GD', 'IGD', 'Spacing', 'Spread', 'Hypervolume'};
% Print results
fprintf('\nPerformance Metrics over %d runs:\n', num_runs);
for i = 1:length(metrics)
fprintf('\n%s Statistics:\n', metric_names{i});
fprintf(' Best = %.4f\n', min(metrics{i}));
fprintf(' Average = %.4f\n', mean(metrics{i}));
fprintf(' Median = %.4f\n', median(metrics{i}));
fprintf(' Worst = %.4f\n', max(metrics{i}));
fprintf(' Std Dev = %.4f\n', std(metrics{i}));
end
%% Plot Results
% Plot Pareto Front
figure;
plot(paretoFront(:,1), paretoFront(:,2), 'ro', 'MarkerSize', 6, 'MarkerFaceColor', 'r');
hold on;
plot(true_pareto(:,1), true_pareto(:,2), 'b--', 'LineWidth', 1.5);
xlabel('f1');
ylabel('f2');
title('Pareto Front for zdt4 using MOEMA');
legend('Obtained Pareto Front', 'True Pareto Front');
grid on;
hold off;
% Plot the Pareto front
figure;
plot(paretoFront(:,1), paretoFront(:,2), 'ro');
xlabel('f1');
ylabel('f2');
title('Pareto Front for ZDT4 using MOEMA');
grid on;
hold off;
% Calculate and display metrics
disp(['Number of solutions in Pareto Front: ', num2str(size(paretoFront, 1))]);
disp(['Minimum f1: ', num2str(min(paretoFront(:,1)))]);
disp(['Minimum f2: ', num2str(min(paretoFront(:,2)))]);
function [paretoFront, paretoSet, Trajectories, fitness_history, position_history] = MOEMA(NPop, T, lb, ub, dim, FITNESSFCN)
% Multi-Objective Evolutionary Mating Algorithm (MOEMA)
% MOEMA attempts to solve multi-objective problems of the following form:
% min F(X) = [f1(X), f2(X), ..., fm(X)] subject to: lb <= X <= ub
%
% FITNESSFCN - function handle of the multi-objective fitness function
% lb - lower bounds on X
% ub - upper bounds on X
% T - number of iterations
% NPop - size of the population (class size)
% paretoFront - set of non-dominated solutions in objective space
% paretoSet - set of non-dominated solutions in decision space
% Trajectories - best solution trajectory over iterations
% fitness_history - history of fitness values
% position_history - history of population positions
% Preallocation
fitness_history = zeros(NPop, T, 2); % Assuming 2 objectives
position_history = zeros(NPop, T, dim);
Trajectories = zeros(dim, T);
% Generation of initial population
pop = initialization(NPop, dim, ub, lb);
pop_dad = pop(1:NPop/2, :); % Eqn (1) in the paper
pop_mum = pop(NPop/2+1:end, :); % Eqn (2) in the paper
% Evaluation of objective function for initialization
obj_dad = evaluate_population(pop_dad, FITNESSFCN);
obj_mum = evaluate_population(pop_mum, FITNESSFCN);
obj = [obj_dad; obj_mum];
% Initialize archive with non-dominated solutions
archive = update_archive([], [], pop, obj);
for gen = 1:T % Start iteration
k1 = randperm(NPop/2);
k2 = randperm(NPop/2);
k1x = pop_dad(k1, :);
k2x = pop_mum(k2, :);
for k = 1:NPop/2
% Mating process
Imates = 1 + (var(k1x(k, :)) - var(k2x(k, :))); % Eqn (4) in the paper
p = randn(1, dim); % Eqn (6) in the paper
r = rand();
if Imates >= 0 % Upper part of Eqn (5) in the paper
temp1 = p .* k1x(k, :);
temp2 = (1-p) .* k2x(k, :);
pop_child = temp1 + temp2;
else % Lower part of Eqn (5) in the paper
temp1 = p .* k2x(k, :);
temp2 = (1-p) .* k1x(k, :);
pop_child = temp1 + temp2;
end
% Impact of best environment (using random non-dominated solution)
if ~isempty(archive.pop)
best_pop = archive.pop(randi(size(archive.pop, 1)), :);
K = rand(1, dim) < 0.8; % Eqn (8) in the paper
pop_child = best_pop .* (1-K) + pop_child .* K; % Eqn (7) in the paper
end
% Bounding of the solution
pop_child = max(min(ub, pop_child), lb);
% Evaluate child
obj_child = FITNESSFCN(pop_child);
% Update population using non-dominated sorting
[pop_dad, obj_dad] = update_population(pop_dad, obj_dad, pop_child, obj_child, k1(k));
[pop_mum, obj_mum] = update_population(pop_mum, obj_mum, pop_child, obj_child, k2(k));
% Exploration process
if r < 0.5
if ~isempty(archive.pop)
temp1 = archive.pop(randi(size(archive.pop, 1)), :);
pop_child = rand(1, dim) .* temp1; % Eqn (9) in the paper
pop_child = max(min(ub, pop_child), lb);
obj_child = FITNESSFCN(pop_child);
% Update population using non-dominated sorting
[pop_dad, obj_dad] = update_population(pop_dad, obj_dad, pop_child, obj_child, randi(NPop/2));
[pop_mum, obj_mum] = update_population(pop_mum, obj_mum, pop_child, obj_child, randi(NPop/2));
end
end
end
% Combine updated pool dad and mum
pop = [pop_dad; pop_mum];
obj = [obj_dad; obj_mum];
% Update archive
archive = update_archive(archive.pop, archive.obj, pop, obj);
% Record the history at each iteration
fitness_history(:, gen, :) = obj;
position_history(:, gen, :) = pop;
% Update Trajectories (using a random non-dominated solution)
if ~isempty(archive.pop)
Trajectories(:, gen) = archive.pop(randi(size(archive.pop, 1)), :)';
end
disp(['Iteration #', num2str(gen), ' completed. Archive size: ', num2str(size(archive.pop, 1))]);
end
% Return Pareto front and Pareto set
paretoFront = archive.obj;
paretoSet = archive.pop;
end
% Helper Functions
function pop = initialization(NPop, dim, ub, lb)
pop = rand(NPop, dim) .* (ub - lb) + lb;
end
function obj = evaluate_population(pop, FITNESSFCN)
obj = zeros(size(pop, 1), 2); % Assuming 2 objectives
for i = 1:size(pop, 1)
obj(i, :) = FITNESSFCN(pop(i, :));
end
end
function [pop, obj] = update_population(pop, obj, child, child_obj, index)
if dominates(child_obj, obj(index, :))
pop(index, :) = child;
obj(index, :) = child_obj;
elseif ~dominates(obj(index, :), child_obj)
if rand() < 0.5
pop(index, :) = child;
obj(index, :) = child_obj;
end
end
end
function archive = update_archive(archive_pop, archive_obj, pop, obj)
combined_pop = [archive_pop; pop];
combined_obj = [archive_obj; obj];
[is_dominated, ~] = pareto_dominance(combined_obj);
non_dominated = ~is_dominated;
archive.pop = combined_pop(non_dominated, :);
archive.obj = combined_obj(non_dominated, :);
end
function [is_dominated, domination_count] = pareto_dominance(objectives)
[n_solutions, n_objectives] = size(objectives);
is_dominated = false(n_solutions, 1);
domination_count = zeros(n_solutions, 1);
for i = 1:n_solutions
for j = i+1:n_solutions
if dominates(objectives(i, :), objectives(j, :))
is_dominated(j) = true;
domination_count(i) = domination_count(i) + 1;
elseif dominates(objectives(j, :), objectives(i, :))
is_dominated(i) = true;
domination_count(j) = domination_count(j) + 1;
end
end
end
end
function result = dominates(obj1, obj2)
result = all(obj1 <= obj2) && any(obj1 < obj2);
end
%% Performance Metrics Functions
function [GD, IGD, Spacing, Spread] = calculate_metrics(obtained_pf, true_pf)
GD = generational_distance(obtained_pf, true_pf);
IGD = inverted_generational_distance(obtained_pf, true_pf);
Spacing = spacing_metric(obtained_pf);
Spread = spread_metric(obtained_pf, true_pf);
end
function gd = generational_distance(obtained_pf, true_pf)
if isempty(obtained_pf)
gd = Inf;
return;
end
distances = zeros(size(obtained_pf, 1), 1);
for i = 1:size(obtained_pf, 1)
dist_to_true = sqrt(sum((true_pf - obtained_pf(i,:)).^2, 2));
distances(i) = min(dist_to_true);
end
p = 2;
gd = (sum(distances.^p)/size(obtained_pf, 1))^(1/p);
end
function igd = inverted_generational_distance(obtained_pf, true_pf)
if isempty(obtained_pf)
igd = Inf;
return;
end
distances = zeros(size(true_pf, 1), 1);
for i = 1:size(true_pf, 1)
dist_to_obtained = sqrt(sum((obtained_pf - true_pf(i,:)).^2, 2));
distances(i) = min(dist_to_obtained);
end
p = 1;
igd = (sum(distances.^p)/size(true_pf, 1))^(1/p);
end
function spacing = spacing_metric(obtained_pf)
if size(obtained_pf, 1) < 2
spacing = Inf;
return;
end
n = size(obtained_pf, 1);
distances = zeros(n, 1);
for i = 1:n
min_dist = Inf;
for j = 1:n
if i ~= j
dist = sqrt(sum((obtained_pf(i,:) - obtained_pf(j,:)).^2));
min_dist = min(min_dist, dist);
end
end
distances(i) = min_dist;
end
d_bar = mean(distances);
spacing = sqrt(sum((distances - d_bar).^2)/(n-1));
end
function spread = spread_metric(obtained_pf, true_pf)
if size(obtained_pf, 1) < 2
spread = Inf;
return;
end
[~, idx] = sort(obtained_pf(:,1));
obtained_pf = obtained_pf(idx,:);
df = zeros(1, 2);
for i = 1:2
obtained_extreme = obtained_pf(i,:);
true_extreme = true_pf(i,:);
df(i) = sqrt(sum((obtained_extreme - true_extreme).^2));
end
distances = sqrt(sum(diff(obtained_pf).^2, 2));
d_bar = mean(distances);
spread = (sum(df) + sum(abs(distances - d_bar))) / (sum(df) + size(obtained_pf,1) * d_bar);
end
function hypervolume = calculate_hypervolume(obtained_pf, reference_point)
if isempty(obtained_pf)
hypervolume = 0;
return;
end
[obtained_pf, ~] = sortrows(obtained_pf);
hypervolume = 0;
for i = 1:size(obtained_pf, 1)
if i == 1
width = reference_point(1) - obtained_pf(i,1);
else
width = obtained_pf(i,1) - obtained_pf(i-1,1);
end
height = reference_point(2) - obtained_pf(i,2);
hypervolume = hypervolume + width * height;
end
end
% ------------------- True Pareto Front Function -------------------
function true_pareto = get_true_pareto()
% Generate the true Pareto front for ZDT4
num_points = 100; % Number of points in the true Pareto front
f1 = linspace(0, 1, num_points)';
f2 = 1 - sqrt(f1);
true_pareto = [f1, f2];
end
function f = zdt4(x)
f1 = x(1);
g = 1 + 10*(length(x)-1);
for i = 2:length(x)
g = g + (x(i)^2 - 10*cos(4*pi*x(i)));
end
h = 1 - sqrt(f1/g);
f2 = g*h;
f = [f1, f2];
end