can you try on zdt2:
clc; 
clear; 
close all;

%% ======================= Parameter Initialization =======================

% Problem Parameters for ZDT3
NPop = 50;          % Population size
T = 200;             % Number of iterations
dim = 30;            % Number of decision variables
lb = zeros(1, dim);  % Lower bound
ub = ones(1, dim);   % Upper bound

% Multiple Runs Parameters
num_runs = 10;       % Define the number of runs

%% ======================= Metrics Storage Initialization =======================

% Initialize metrics storage
GD_values = zeros(num_runs,1);
IGD_values = zeros(num_runs,1);
Convergence_values = zeros(num_runs,1);
Diversity_values = zeros(num_runs,1);

% Arrays to store GD statistics
best_GD_values = zeros(num_runs,1);
average_GD_values = zeros(num_runs,1);
median_GD_values = zeros(num_runs,1);
worst_GD_values = zeros(num_runs,1);

% Arrays to store IGD statistics
best_IGD_values = zeros(num_runs,1);
average_IGD_values = zeros(num_runs,1);
median_IGD_values = zeros(num_runs,1);
worst_IGD_values = zeros(num_runs,1);

% Initialize cell array to store Pareto fronts from all runs
paretoFront_allRuns = cell(num_runs,1);

%% ======================= Define True Pareto Front for ZDT3 =======================

% User-Provided True Pareto Front
x1 = linspace(0, 1, 500);
x2 = 1 - sqrt(x1) - x1.*sin(10*pi*x1);
true_pareto_front = [x1', x2'];

% Optional: Remove any points with x2 < 0 if desired
% Uncomment the following lines to exclude invalid points
% valid_indices = x2 >= 0;
% true_pareto_front = [x1(valid_indices)', x2(valid_indices)'];

%% ======================= Multiple Runs of MOEMA =======================

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, @zdt3);
    
    % Store Pareto front of the current run
    paretoFront_allRuns{run} = paretoFront;
    
    % Calculate Generational Distance (GD)
    GD = generational_distance(paretoFront, true_pareto_front);
    
    % Calculate Inverted Generational Distance (IGD) using user-provided function
    IGD = inverted_generational_distance(paretoFront, true_pareto_front);
    
    % Calculate Convergence Metric (using GD as an example)
    Convergence = GD; % You can define a different convergence metric if needed
    
    % Calculate Diversity Metric
    Diversity = diversity_metric(paretoFront, true_pareto_front);
    
    % Store metrics
    GD_values(run) = GD;
    IGD_values(run) = IGD;
    Convergence_values(run) = Convergence;
    Diversity_values(run) = Diversity;
    
    % Store GD statistics
    best_GD_values(run) = min(GD_values(1:run));
    average_GD_values(run) = mean(GD_values(1:run));
    median_GD_values(run) = median(GD_values(1:run));
    worst_GD_values(run) = max(GD_values(1:run));
    
    % Store IGD statistics
    best_IGD_values(run) = min(IGD_values(1:run));
    average_IGD_values(run) = mean(IGD_values(1:run));
    median_IGD_values(run) = median(IGD_values(1:run));
    worst_IGD_values(run) = max(IGD_values(1:run));
    
    % Optionally, you can plot each run's Pareto front
    % Uncomment the following lines if you wish to visualize each run
    % figure;
    % plot(paretoFront(:,1), paretoFront(:,2), 'ro', 'DisplayName', 'MOEMA');
    % hold on;
    % plot(true_pareto_front(:,1), true_pareto_front(:,2), 'b-', 'DisplayName', 'True Pareto Front');
    % xlabel('f1');
    % ylabel('f2');
    % title(sprintf('MOEMA vs True Pareto Front for ZDT3 (Run #%d)', run));
    % legend('show');
    % grid on;
    
end

%% ======================= Identify and Display the Best Pareto Front =======================

% Define which metric to use for selecting the best run
% Options: 'GD', 'IGD', 'Convergence', 'Diversity'
best_metric = 'GD';  % Change this to your preferred metric ('GD', 'IGD', etc.)

switch best_metric
    case 'GD'
        [~, best_run] = min(GD_values);
    case 'IGD'
        [~, best_run] = min(IGD_values);
    case 'Convergence'
        [~, best_run] = min(Convergence_values);
    case 'Diversity'
        [~, best_run] = max(Diversity_values); % Assuming higher diversity is better
    otherwise
        error('Invalid metric selected. Choose from ''GD'', ''IGD'', ''Convergence'', or ''Diversity''.');
end

% Retrieve the best Pareto front
best_paretoFront = paretoFront_allRuns{best_run};

% Plot the best Pareto front alongside the true Pareto front
figure;
plot(best_paretoFront(:,1), best_paretoFront(:,2), 'ro', 'DisplayName', sprintf('Best MOEMA Run #%d (%s)', best_run, best_metric));
hold on;
plot(true_pareto_front(:,1), true_pareto_front(:,2), 'b-', 'DisplayName', 'True Pareto Front');
xlabel('f1');
ylabel('f2');
title(sprintf('Best MOEMA Pareto Front vs True Pareto Front based on %s', best_metric));
legend('show');
grid on;

fprintf('\nBest Pareto Front is from Run #%d based on %s.\n', best_run, best_metric);

%% ======================= Compute Final Statistics =======================

% Compute Final Statistics for GD
final_best_GD = min(GD_values);
final_average_GD = mean(GD_values);
final_median_GD = median(GD_values);
final_worst_GD = max(GD_values);

% Compute Final Statistics for IGD
final_best_IGD = min(IGD_values);
final_average_IGD = mean(IGD_values);
final_median_IGD = median(IGD_values);
final_worst_IGD = max(IGD_values);

%% ======================= Display Metrics =======================

fprintf('\nPerformance Metrics over %d runs:\n', num_runs);

% Generational Distance (GD) Statistics
fprintf('Generational Distance (GD):\n');
fprintf('  Best    = %.4f\n', final_best_GD);
fprintf('  Average = %.4f\n', final_average_GD);
fprintf('  Median  = %.4f\n', final_median_GD);
fprintf('  Worst   = %.4f\n', final_worst_GD);

% Inverted Generational Distance (IGD) Statistics
fprintf('Inverted Generational Distance (IGD):\n');
fprintf('  Best    = %.4f\n', final_best_IGD);
fprintf('  Average = %.4f\n', final_average_IGD);
fprintf('  Median  = %.4f\n', final_median_IGD);
fprintf('  Worst   = %.4f\n', final_worst_IGD);

% Convergence Metric Statistics
fprintf('Convergence Metric (GD as example):\n');
fprintf('  Mean = %.4f, Std = %.4f\n', mean(Convergence_values), std(Convergence_values));

% Diversity Metric Statistics
fprintf('Diversity Metric:\n');
fprintf('  Mean = %.4f, Std = %.4f\n', mean(Diversity_values), std(Diversity_values));

%% ======================= Optional: Plot Metrics =======================

% Plot GD over Runs
figure;
plot(1:num_runs, GD_values, 'b-o', 'LineWidth', 1.5);
xlabel('Run');
ylabel('Generational Distance (GD)');
title('GD over Multiple Runs');
grid on;

% Plot IGD over Runs
figure;
plot(1:num_runs, IGD_values, 'r-o', 'LineWidth', 1.5);
xlabel('Run');
ylabel('Inverted Generational Distance (IGD)');
title('IGD over Multiple Runs');
grid on;

% Plot Convergence over Runs
figure;
plot(1:num_runs, Convergence_values, 'g-o', 'LineWidth', 1.5);
xlabel('Run');
ylabel('Convergence Metric (GD)');
title('Convergence over Multiple Runs');
grid on;

% Plot Diversity over Runs
figure;
plot(1:num_runs, Diversity_values, 'm-o', 'LineWidth', 1.5);
xlabel('Run');
ylabel('Diversity Metric');
title('Diversity over Multiple Runs');
grid on;

%% ======================= Display Aggregate Comparison Table =======================

fprintf('\nAggregate Comparison Table:\n');
disp('------------------------------------------------------------');
disp('| Metric          | Best      | Average    | Median   | Worst   |');
disp('------------------------------------------------------------');
fprintf('| GD              | %-9.4f | %-10.4f | %-8.4f | %-7.4f |\n', final_best_GD, final_average_GD, final_median_GD, final_worst_GD);
fprintf('| IGD             | %-9.4f | %-10.4f | %-8.4f | %-7.4f |\n', final_best_IGD, final_average_IGD, final_median_IGD, final_worst_IGD);
fprintf('| Convergence     | %-9.4f | %-10.4f | %-8.4f | %-7.4f |\n', mean(Convergence_values), std(Convergence_values), median(Convergence_values), max(Convergence_values));
fprintf('| Diversity       | %-9.4f | %-10.4f | %-8.4f | %-7.4f |\n', mean(Diversity_values), std(Diversity_values), median(Diversity_values), max(Diversity_values));
disp('------------------------------------------------------------');

%% ======================= Function Definitions =======================

% ======================= MOEMA Algorithm =======================
function [paretoFront, paretoSet, Trajectories, fitness_history, position_history] = MOEMA(NPop, T, lb, ub, dim, FITNESSFCN)
    % Multi-Objective Evolutionary Mating Algorithm (MOEMA)
    % Attempts to solve multi-objective optimization problems
    
    % 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, :);                 % First half
    pop_mum = pop(NPop/2+1:end, :);             % Second half
    
    % 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, k);
            [pop_mum, obj_mum] = update_population(pop_mum, obj_mum, pop_child, obj_child, 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

% ======================= Initialization Function =======================
function pop = initialization(NPop, dim, ub, lb)
    pop = rand(NPop, dim) .* (ub - lb) + lb;
end

% ======================= Evaluate Population =======================
function obj = evaluate_population(pop, FITNESSFCN)
    obj = zeros(size(pop, 1), 2);  % Assuming 2 objectives
    parfor i = 1:size(pop, 1)      % Use parfor for faster computation if Parallel Toolbox is available
        obj(i, :) = FITNESSFCN(pop(i, :));
    end
end

% ======================= Update Population Function =======================
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

% ======================= Update Archive Function =======================
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

% ======================= Pareto Dominance Function =======================
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

% ======================= Dominance Check Function =======================
function result = dominates(obj1, obj2)
    result = all(obj1 <= obj2) && any(obj1 < obj2);
end

% ======================= ZDT3 Fitness Function =======================
function f = zdt3(x)
    f1 = x(1);
    g = 1 + 9 * sum(x(2:end)) / (length(x) - 1);
    h = 1 - sqrt(f1 / g) - (f1 / g) * sin(10 * pi * f1);
    f2 = g * h;
    f = [f1, f2];
end

% ======================= Generational Distance (GD) Function =======================
function gd = generational_distance(obtained_pf, true_pf)
    % Calculate Generational Distance (GD)
    distances = zeros(size(obtained_pf,1),1);
    for i = 1:size(obtained_pf,1)
        % Find the minimum distance from obtained_pf(i,:) to true_pf
        distances(i) = min( sqrt( sum( (true_pf - obtained_pf(i,:)).^2, 2 ) ) );
    end
    gd = mean(distances);
end

% ======================= Inverted Generational Distance (IGD) Function =======================
function igd = inverted_generational_distance(obtained_pf, true_pf)
    % Calculate Inverted Generational Distance (IGD)
    distances = sqrt(sum((true_pf - closest_point(true_pf, obtained_pf)).^2, 2));
    igd = mean(distances);
end

% ======================= Diversity Metric Function =======================
function spread = diversity_metric(obtained_pf, true_pf)
    % Calculate Diversity Metric (Spread)
    sorted_pf = sortrows(obtained_pf, 1); % Sort based on f1
    num_obtained = size(sorted_pf, 1);
    
    % Handle edge case with insufficient points
    if num_obtained < 2
        spread = Inf;
        return;
    end
    
    % Calculate gaps between consecutive points
    gaps = sqrt(sum(diff(sorted_pf).^2, 2));
    spread = sum(gaps) / (sqrt(2) * (num_obtained - 1));
end

% ======================= Closest Point Helper Function =======================
function closest = closest_point(src, target)
    % Find the closest point in target for each point in src
    num_src = size(src, 1);
    closest = zeros(size(src));
    
    for i = 1:num_src
        distances = sqrt(sum((target - src(i,:)).^2, 2));
        [~, idx] = min(distances);
        closest(i,:) = target(idx,:);
    end
end