how to improve this code:
% MOEMA_PerformanceMetrics_Full.m
% Enhanced implementation with NSGA-inspired crowding distance and Levy flights
% Using zdt2 benchmark function

clc;
clear;
close all;

%% Parameters
NPop = 500;              % Population size
T = 1000;               % Number of iterations
dim = 30;              % Number of decision variables
lb = zeros(1, dim);    % Lower bound
ub = ones(1, dim);     % Upper bound
num_runs = 5;          % Number of runs for statistical evaluation

% Levy flight parameters
beta = 1.5;            % Levy flight parameter (1 < beta < 2)
sigma_u = (gamma(1+beta)*sin(pi*beta/2)/(gamma((1+beta)/2)*beta*2^((beta-1)/2)))^(1/beta);
sigma_v = 1;

% 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];

%% 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, @zdt2, beta, sigma_u, sigma_v);
    
    % 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
%% Plot Results per Run
for run = 1:num_runs
    % Plot Pareto Front for each run
    figure;
    plot(squeeze(fitness_history(:, run, 1)), squeeze(fitness_history(:, run, 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 - Run #' num2str(run)]);
    legend('Obtained Pareto Front', 'True Pareto Front');
    grid on;
    hold off;
    
    % Display metrics for this run
    fprintf('\nRun #%d Metrics:\n', run);
    fprintf('  GD          = %.4f\n', GD_values(run));
    fprintf('  IGD         = %.4f\n', IGD_values(run));
    fprintf('  Spacing     = %.4f\n', Spacing_values(run));
    fprintf('  Spread      = %.4f\n', Spread_values(run));
    fprintf('  Hypervolume = %.4f\n', Hypervolume_values(run));
end

%% Calculate Statistics and Plot Results
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
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 zdt2 using Enhanced MOEMA');
legend('Obtained Pareto Front', 'True Pareto Front');
grid on;

%% Main MOEMA Function
function [paretoFront, paretoSet, Trajectories, fitness_history, position_history] = ...
    MOEMA(NPop, T, lb, ub, dim, FITNESSFCN, beta, sigma_u, sigma_v)
    
    % Initialize storage
    fitness_history = zeros(NPop, T, 2);
    position_history = zeros(NPop, T, dim);
    Trajectories = zeros(dim, T);
    
    % Initialize population
    pop = initialization(NPop, dim, ub, lb);
    pop_dad = pop(1:NPop/2, :);
    pop_mum = pop(NPop/2+1:end, :);
    
    % Evaluate initial population
    obj_dad = evaluate_population(pop_dad, FITNESSFCN);
    obj_mum = evaluate_population(pop_mum, FITNESSFCN);
    obj = [obj_dad; obj_mum];
    
    % Initialize archive with crowding distance
    archive = update_archive_with_crowding([], [], pop, obj);
    
    % Main loop
    for gen = 1:T
        k1 = randperm(NPop/2);
        k2 = randperm(NPop/2);
        k1x = pop_dad(k1, :);
        k2x = pop_mum(k2, :);
        
        for k = 1:NPop/2
            % Enhanced mating process with crowding distance consideration
            [I_mates, selected_mate] = select_mate_with_crowding(k1x(k,:), k2x(k,:), archive);
            
            % Generate offspring using selected mate
            p = randn(1, dim);
            if I_mates >= 0
                pop_child = p .* k1x(k,:) + (1-p) .* selected_mate;
            else
                pop_child = p .* selected_mate + (1-p) .* k1x(k,:);
            end
            
            % Levy flight exploration
            if rand() < 0.25
                pop_child = levy_flight_exploration(pop_child, archive, beta, sigma_u, sigma_v, lb, ub);
            end
            
            % Bound solutions
            pop_child = max(min(ub, pop_child), lb);
            
            % Evaluate child
            obj_child = FITNESSFCN(pop_child);
            
            % Update populations with crowding distance consideration
            [pop_dad, obj_dad] = update_population_with_crowding(pop_dad, obj_dad, pop_child, obj_child, k1(k));
            [pop_mum, obj_mum] = update_population_with_crowding(pop_mum, obj_mum, pop_child, obj_child, k2(k));
        end
        
        % Combine populations
        pop = [pop_dad; pop_mum];
        obj = [obj_dad; obj_mum];
        
        % Update archive with crowding distance
        archive = update_archive_with_crowding(archive.pop, archive.obj, pop, obj);
        
        % Record history
        fitness_history(:, gen, :) = obj;
        position_history(:, gen, :) = pop;
        
        if ~isempty(archive.pop)
            Trajectories(:, gen) = archive.pop(randi(size(archive.pop, 1)), :)';
        end
    end
    
    paretoFront = archive.obj;
    paretoSet = archive.pop;
end

%% Enhanced Crowding Distance Functions
function crowding_distances = calculate_crowding_distance(objectives)
    [n_solutions, n_objectives] = size(objectives);
    
    if n_solutions <= 2
        crowding_distances = inf(n_solutions, 1);
        return;
    end
    
    crowding_distances = zeros(n_solutions, 1);
    
    for m = 1:n_objectives
        [sorted_obj, indices] = sort(objectives(:,m));
        
        % Set boundary points to infinity
        crowding_distances(indices(1)) = inf;
        crowding_distances(indices(end)) = inf;
        
        % Calculate crowding distance for intermediate points
        obj_range = sorted_obj(end) - sorted_obj(1);
        if obj_range == 0
            obj_range = 1;
        end
        
        for i = 2:(n_solutions-1)
            crowding_distances(indices(i)) = crowding_distances(indices(i)) + ...
                (sorted_obj(i+1) - sorted_obj(i-1)) / obj_range;
        end
    end
end

function archive = update_archive_with_crowding(archive_pop, archive_obj, pop, obj)
    combined_pop = [archive_pop; pop];
    combined_obj = [archive_obj; obj];
    
    % Get non-dominated solutions
    [is_dominated, ~] = pareto_dominance(combined_obj);
    non_dominated = ~is_dominated;
    
    % Calculate crowding distances for non-dominated solutions
    selected_obj = combined_obj(non_dominated, :);
    crowding_distances = calculate_crowding_distance(selected_obj);
    
    % Sort by crowding distance if too many solutions
    if sum(non_dominated) > 100  % Maximum archive size
        [~, sorted_indices] = sort(crowding_distances, 'descend');
        selected_indices = sorted_indices(1:100);
        
        archive.pop = combined_pop(non_dominated, :);
        archive.obj = selected_obj;
        archive.crowding = crowding_distances;
        
        archive.pop = archive.pop(selected_indices, :);
        archive.obj = archive.obj(selected_indices, :);
        archive.crowding = archive.crowding(selected_indices);
    else
        archive.pop = combined_pop(non_dominated, :);
        archive.obj = selected_obj;
        archive.crowding = crowding_distances;
    end
end

%% Levy Flight Functions
function child = levy_flight_exploration(child, archive, beta, sigma_u, sigma_v, lb, ub)
    % Generate Levy flight step
    u = normrnd(0, sigma_u, [1, length(child)]);
    v = normrnd(0, sigma_v, [1, length(child)]);
    step = u./abs(v).^(1/beta);
    
    % Select random archive solution as reference
    if ~isempty(archive.pop)
        ref_solution = archive.pop(randi(size(archive.pop, 1)), :);
        
        % Apply Levy flight step
        child = child + step .* (ref_solution - child);
        
        % Bound the solution
        child = max(min(ub, child), lb);
    end
end

%% Enhanced Selection and Update Functions
function [I_mates, selected_mate] = select_mate_with_crowding(solution1, solution2, archive)
    if isempty(archive.pop)
        I_mates = 1 + (var(solution1) - var(solution2));
        selected_mate = solution2;
        return;
    end
    
    % Calculate diversity contribution
    distances1 = pdist2(solution1, archive.pop);
    distances2 = pdist2(solution2, archive.pop);
    
    diversity1 = min(distances1);
    diversity2 = min(distances2);
    
    % Consider both variance and diversity
    I_mates = 1 + (var(solution1) - var(solution2)) + 0.5*(diversity1 - diversity2);
    
    % Select mate based on crowding distance
    if I_mates >= 0
        selected_mate = solution2;
    else
        selected_mate = solution1;
    end
end

function [pop, obj] = update_population_with_crowding(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)
        % Calculate crowding distance for both solutions
        combined_obj = [obj; child_obj];
        crowding_distances = calculate_crowding_distance(combined_obj);
        
        % Replace if child has better crowding distance
        if crowding_distances(end) > crowding_distances(index)
            pop(index, :) = child;
            obj(index, :) = child_obj;
        end
    end
end

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

function obj = evaluate_population(pop, FITNESSFCN)
    num_sol = size(pop, 1);
    obj = zeros(num_sol, 2);
    parfor i = 1:num_sol
        obj(i, :) = FITNESSFCN(pop(i, :));
    end
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);
    
    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 true_pareto = get_true_pareto()
    % Generate the true Pareto front for ZDT2
    num_points = 100;
    f1 = linspace(0, 1, num_points)';
    f2 = 1 - (f1).^2;  % ZDT2 has a concave Pareto front
    true_pareto = [f1, f2];
end

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

function f = zdt2(x)
    % ZDT2 Test Problem
    f1 = x(1);
    g = 1 + 9 * sum(x(2:end)) / (length(x) - 1);
    h = 1 - (f1/g)^2;  % This is the key difference from zdt2
    f2 = g * h;
    f = [f1, f2];
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 = 2;
    igd = (sum(distances.^p)/size(true_pf, 1))^(1/p);
end

function sp = spacing_metric(obtained_pf)
    if size(obtained_pf, 1) < 2
        sp = Inf;
        return;
    end
    
    % Calculate distances between consecutive solutions
    distances = zeros(size(obtained_pf, 1), 1);
    for i = 1:size(obtained_pf, 1)
        min_dist = Inf;
        for j = 1:size(obtained_pf, 1)
            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
    
    % Calculate spacing metric
    mean_dist = mean(distances);
    sp = sqrt(sum((distances - mean_dist).^2) / (size(obtained_pf, 1) - 1));
end

function delta = spread_metric(obtained_pf, true_pf)
    if size(obtained_pf, 1) < 2
        delta = Inf;
        return;
    end
    
    % Find extreme solutions
    f1_min_idx = find(obtained_pf(:,1) == min(obtained_pf(:,1)), 1);
    f1_max_idx = find(obtained_pf(:,1) == max(obtained_pf(:,1)), 1);
    
    % Calculate distances between consecutive solutions
    distances = zeros(size(obtained_pf, 1)-1, 1);
    for i = 1:(size(obtained_pf, 1)-1)
        distances(i) = sqrt(sum((obtained_pf(i+1,:) - obtained_pf(i,:)).^2));
    end
    
    % Calculate mean distance
    mean_dist = mean(distances);
    
    % Calculate distances to extreme solutions
    df = sqrt(sum((obtained_pf(f1_min_idx,:) - obtained_pf(f1_max_idx,:)).^2));
    
    % Calculate spread metric
    delta = (df + sum(abs(distances - mean_dist))) / (df + size(obtained_pf, 1) * mean_dist);
end

function hv = calculate_hypervolume(obtained_pf, reference_point)
    if isempty(obtained_pf)
        hv = 0;
        return;
    end
    
    % Remove dominated points
    [is_dominated, ~] = pareto_dominance(obtained_pf);
    non_dominated = obtained_pf(~is_dominated, :);
    
    % Sort points by first objective
    [~, sort_idx] = sort(non_dominated(:,1));
    sorted_pf = non_dominated(sort_idx, :);
    
    % Calculate hypervolume
    hv = 0;
    for i = 1:size(sorted_pf, 1)
        if i == 1
            width = reference_point(1) - sorted_pf(i,1);
        else
            width = sorted_pf(i-1,1) - sorted_pf(i,1);
        end
        
        if i == size(sorted_pf, 1)
            height = reference_point(2) - sorted_pf(i,2);
        else
            height = reference_point(2) - sorted_pf(i,2);
        end
        
        hv = hv + width * height;
    end
end