can you display the IGD, SP values, and MS values for this code? i want the Best, Avg , Worst, Stddev and Median for each of IGD, SP values, and MS values for this moema code.

% Main script to run MOEMA optimization
clc;
clear;
close all;
% Set problem parameters
NPop = 50; % Population size
T = 50; % Number of iterations
dim = 30; % Number of decision variables
lb = zeros(1, dim); % Lower bound
ub = ones(1, dim); % Upper bound

% Run MOEMA
[paretoFront, paretoSet, Trajectories, fitness_history, position_history] = MOEMA(NPop, T, lb, ub, dim, @zdt1);

% Plot the Pareto front
figure;
plot(paretoFront(:,1), paretoFront(:,2), 'ro');
xlabel('f1');
ylabel('f2');
title('Pareto Front for ZDT1 using MOEMA');

% 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.2
        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, ~] = size(objectives);
is_dominated = false(n_solutions, 1);
domination_count = zeros(n_solutions, 1); % Ensure this variable is initialized


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; % Track how many solutions 'i' dominates
        elseif dominates(objectives(j, :), objectives(i, :))
            is_dominated(i) = true;
            domination_count(j) = domination_count(j) + 1; % Track how many solutions 'j' dominates
        end
    end
end
end

function result = dominates(obj1, obj2)
result = all(obj1 <= obj2) && any(obj1 < obj2);
end

function f = zdt1(x)
f1 = x(1);
g = 1 + 9 * sum(x(2:end)) / (length(x) - 1);
h = 1 - sqrt(f1 / g);
f2 = g * h;
f = [f1, f2];
end