can you try on zdt4? if yes give me full code
clc;
clear;
close all;

% Problem parameters
dim = 30;
lb = zeros(1, dim);
ub = ones(1, dim);

% Run MOEMA
NPop_MOEMA = 50;
T_MOEMA = 200;
[paretoFront_MOEMA, ~] = MOEMA(NPop_MOEMA, T_MOEMA, lb, ub, dim, @zdt2);

% Run NSGA-II
NPop_NSGA2 = 50;
T_NSGA2 = 200;
[paretoFront_NSGA2, ~] = NSGA2(NPop_NSGA2, T_NSGA2, lb, ub, dim);

% True Pareto front
x1 = linspace(0, 1, 100);
true_pareto_front = [x1', 1 - (x1.^2)'];

% Reference point for hypervolume calculation
reference_point = [1.1, 1.1];

% Calculate metrics for MOEMA
gd_MOEMA = calculate_gd(paretoFront_MOEMA, true_pareto_front);
spread_MOEMA = calculate_spread(paretoFront_MOEMA);
hv_MOEMA = calculate_hypervolume(paretoFront_MOEMA, reference_point);
convergence_MOEMA = calculate_convergence(paretoFront_MOEMA, true_pareto_front);

% Calculate metrics for NSGA-II
gd_NSGA2 = calculate_gd(paretoFront_NSGA2, true_pareto_front);
spread_NSGA2 = calculate_spread(paretoFront_NSGA2);
hv_NSGA2 = calculate_hypervolume(paretoFront_NSGA2, reference_point);
convergence_NSGA2 = calculate_convergence(paretoFront_NSGA2, true_pareto_front);

% Display Comparison
disp('Comparison Table:');
disp('-----------------------------------------------------------');
disp('| Metric | MOEMA | NSGA-II |');
disp('-----------------------------------------------------------');
fprintf('| GD | %-14.4f | %-14.4f |\n', gd_MOEMA, gd_NSGA2);
fprintf('| Spread | %-14.4f | %-14.4f |\n', spread_MOEMA, spread_NSGA2);
fprintf('| Hypervolume | %-14.4f | %-14.4f |\n', hv_MOEMA, hv_NSGA2);
fprintf('| Convergence | %-14.4f | %-14.4f |\n', convergence_MOEMA, convergence_NSGA2);
disp('-----------------------------------------------------------');

% Plotting Results
figure;
plot(paretoFront_MOEMA(:,1), paretoFront_MOEMA(:,2), 'ro', 'DisplayName', 'MOEMA');
hold on;
plot(paretoFront_NSGA2(:,1), paretoFront_NSGA2(:,2), 'bo', 'DisplayName', 'NSGA-II');
plot(true_pareto_front(:,1), true_pareto_front(:,2), 'k-', 'DisplayName', 'True Pareto Front');
xlabel('f1');
ylabel('f2');
title('Comparison of MOEMA and NSGA-II on zdt2');
legend('show');
grid on;

% Plotting Results
figure;

% MOEMA vs True Pareto Front
subplot(1, 2, 1); % 1 row, 2 columns, first plot
plot(paretoFront_MOEMA(:,1), paretoFront_MOEMA(:,2), 'ro', 'DisplayName', 'MOEMA');
hold on;
plot(true_pareto_front(:,1), true_pareto_front(:,2), 'k-', 'DisplayName', 'True Pareto Front');
xlabel('f1');
ylabel('f2');
title('MOEMA vs True Pareto Front for zdt2');
legend('show');
grid on;

% NSGA-II vs True Pareto Front
subplot(1, 2, 2); % 1 row, 2 columns, second plot
plot(paretoFront_NSGA2(:,1), paretoFront_NSGA2(:,2), 'bo', 'DisplayName', 'NSGA-II');
hold on;
plot(true_pareto_front(:,1), true_pareto_front(:,2), 'k-', 'DisplayName', 'True Pareto Front');
xlabel('f1');
ylabel('f2');
title('NSGA-II vs True Pareto Front for zdt2');
legend('show');
grid on;

% Function Definitions

function [paretoFront, paretoSet, Trajectories, fitness_history, position_history] = MOEMA(NPop, T, lb, ub, dim, FITNESSFCN)
fitness_history = zeros(NPop, T, 2);
position_history = zeros(NPop, T, dim);
Trajectories = zeros(dim, T);
pop = initialization(NPop, dim, ub, lb);
pop_dad = pop(1:NPop/2, :);
pop_mum = pop(NPop/2+1:end, :);
obj_dad = evaluate_population(pop_dad, FITNESSFCN);
obj_mum = evaluate_population(pop_mum, FITNESSFCN);
obj = [obj_dad; obj_mum];
archive = update_archive([], [], pop, obj);


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
        Imates = 1 + (var(k1x(k,:)) - var(k2x(k,:)));
        p = randn(1, dim);
        r = rand();
        
        if Imates >= 0
            temp1 = p .* k1x(k,:);
            temp2 = (1-p) .* k2x(k,:);
            pop_child = temp1 + temp2;
        else
            temp1 = p .* k2x(k,:);
            temp2 = (1-p) .* k1x(k,:);
            pop_child = temp1 + temp2;
        end
        
        if ~isempty(archive.pop)
            best_pop = archive.pop(randi(size(archive.pop, 1)), :);
            K = rand(1, dim) < 0.8;
            pop_child = best_pop .* (1-K) + pop_child .* K;
        end
        
        pop_child = max(min(ub, pop_child), lb);
        obj_child = FITNESSFCN(pop_child);
        [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));
        
        if r < 0.2
            if ~isempty(archive.pop)
                temp1 = archive.pop(randi(size(archive.pop, 1)), :);
                pop_child = rand(1, dim) .* temp1;
                pop_child = max(min(ub, pop_child), lb);
                obj_child = FITNESSFCN(pop_child);
                [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
    
    pop = [pop_dad; pop_mum];
    obj = [obj_dad; obj_mum];
    archive = update_archive(archive.pop, archive.obj, pop, obj);
    fitness_history(:, gen, :) = obj;
    position_history(:, gen, :) = pop;
    
    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

paretoFront = archive.obj;
paretoSet = archive.pop;
end

function [paretoFront, paretoSet] = NSGA2(NPop, T, lb, ub, dim)
pop = initialization(NPop, dim, ub, lb);
fitness = evaluate_population(pop, @zdt2);


for gen = 1:T
    [fronts, ranks] = non_dominated_sort(fitness);
    mating_pool = tournament_selection(pop, ranks, NPop);
    offspring = crossover(mating_pool, 0.9, dim);
    offspring = mutation(offspring, 1/dim, dim);
    offspring_fitness = evaluate_population(offspring, @zdt2);
    comb_pop = [pop; offspring];
    comb_fitness = [fitness; offspring_fitness];
    [fronts, ~] = non_dominated_sort(comb_fitness);
    [pop, fitness] = select_new_population(comb_pop, comb_fitness, fronts, NPop);
    disp(['Generation #', num2str(gen), ' completed.']);
end

paretoFront = fitness;
paretoSet = pop;
end

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

function fitness = evaluate_population(pop, FITNESSFCN)
fitness = zeros(size(pop, 1), 2);
for i = 1:size(pop, 1)
fitness(i, :) = FITNESSFCN(pop(i, :));
end
end

function [fronts, rank] = non_dominated_sort(fitness)
num_individuals = size(fitness, 1);
domination_count = zeros(num_individuals, 1);
dominated_set = cell(num_individuals, 1);
rank = zeros(num_individuals, 1);


fronts = cell(1);
for i = 1:num_individuals
    for j = 1:num_individuals
        if i ~= j
            if dominates(fitness(i, :), fitness(j, :))
                dominated_set{i} = [dominated_set{i}, j];
            elseif dominates(fitness(j, :), fitness(i, :))
                domination_count(i) = domination_count(i) + 1;
            end
        end
    end
    
    if domination_count(i) == 0
        rank(i) = 1;
        fronts{1} = [fronts{1}, i];
    end
end

current_front = 1;
while ~isempty(fronts{current_front})
    next_front = [];
    for i = fronts{current_front}
        for j = dominated_set{i}
            domination_count(j) = domination_count(j) - 1;
            if domination_count(j) == 0
                rank(j) = current_front + 1;
                next_front = [next_front, j];
            end
        end
    end
    current_front = current_front + 1;
    fronts{current_front} = next_front;
end
fronts = fronts(1:current_front-1);
end

function is_dominated = dominates(x, y)
is_dominated = all(x <= y) && any(x < y);
end

function selected_pop = tournament_selection(pop, ranks, num_selected)
selected_pop = zeros(num_selected, size(pop, 2));
for i = 1:num_selected
candidates = randperm(size(pop, 1), 2);
if ranks(candidates(1)) < ranks(candidates(2))
selected_pop(i, :) = pop(candidates(1), :);
else
selected_pop(i, :) = pop(candidates(2), :);
end
end
end

function offspring = crossover(pop, crossover_prob, num_vars)
offspring = zeros(size(pop));
for i = 1:2:size(pop, 1)
if rand < crossover_prob
parent1 = pop(i, :);
parent2 = pop(i+1, :);
cross_point = randi([1, num_vars-1], 1, 1);
offspring(i, :) = [parent1(1:cross_point), parent2(cross_point+1:end)];
offspring(i+1, :) = [parent2(1:cross_point), parent1(cross_point+1:end)];
else
offspring(i, :) = pop(i, :);
offspring(i+1, :) = pop(i+1, :);
end
end
end

function mutated_pop = mutation(pop, mutation_prob, num_vars)
eta_m = 20;
mutated_pop = pop;
for i = 1:size(pop, 1)
for j = 1:num_vars
if rand < mutation_prob
r = rand;
if r < 0.5
delta = (2 * r)^(1/(eta_m + 1)) - 1;
else
delta = 1 - (2 * (1 - r))^(1/(eta_m + 1));
end
mutated_pop(i, j) = pop(i, j) + delta;
mutated_pop(i, j) = min(max(mutated_pop(i, j), 0), 1);
end
end
end
end

function [pop, fitness] = select_new_population(comb_pop, comb_fitness, fronts, NPop)
pop = [];
fitness = [];
for i = 1:length(fronts)
if size(pop, 1) + length(fronts{i}) > NPop
break;
end
pop = [pop; comb_pop(fronts{i}, :)];
fitness = [fitness; comb_fitness(fronts{i}, :)];
end


remaining = NPop - size(pop, 1);
if remaining > 0
    selected_front = fronts{i};
    distances = crowding_distance(comb_fitness(selected_front, :));
    [~, sorted_indices] = sort(distances, 'descend');
    selected = selected_front(sorted_indices(1:remaining));
    pop = [pop; comb_pop(selected, :)];
    fitness = [fitness; comb_fitness(selected, :)];
end
end

function distances = crowding_distance(fitness)
num_individuals = size(fitness, 1);
num_objectives = size(fitness, 2);
distances = zeros(num_individuals, 1);


for m = 1:num_objectives
    [sorted_values, sorted_idx] = sort(fitness(:, m));
    distances(sorted_idx(1)) = inf;
    distances(sorted_idx(end)) = inf;
    for i = 2:num_individuals-1
        if sorted_values(end) ~= sorted_values(1)
            distances(sorted_idx(i)) = distances(sorted_idx(i)) + ...
                (sorted_values(i+1) - sorted_values(i-1)) / (sorted_values(end) - sorted_values(1));
        end
    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 [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 f = zdt2(x)
f1 = x(1);
g = 1 + 9 * sum(x(2:end)) / (length(x) - 1);
h = 1 - (f1 / g)^2;
f2 = g * h;
f = [f1, f2];
end

function spread = calculate_spread(front)
sorted_front = sortrows(front);
distances = sqrt(sum(diff(sorted_front).^2, 2));
spread = std(distances) / mean(distances);
end

function gd = calculate_gd(approx_front, true_front)
sum_of_distances = 0;
for i = 1:size(approx_front, 1)
distances = sqrt(sum((true_front - approx_front(i, :)).^2, 2));
min_distance = min(distances);
sum_of_distances = sum_of_distances + min_distance^2;
end
gd = sqrt(sum_of_distances) / size(approx_front, 1);
end

function hv = calculate_hypervolume(front, reference_point)
sorted_front = sortrows(front);
hypervolume = 0;
for i = 1:size(sorted_front, 1)
if i == 1
width = reference_point(1) - sorted_front(i, 1);
else
width = sorted_front(i-1, 1) - sorted_front(i, 1);
end
height = reference_point(2) - sorted_front(i, 2);
hypervolume = hypervolume + width * height;
end
hv = hypervolume;
end

function convergence = calculate_convergence(approx_front, true_front)
distances = zeros(size(approx_front, 1), 1);
for i = 1:size(approx_front, 1)
dist = sqrt(sum((true_front - approx_front(i, :)).^2, 2));
distances(i) = min(dist);
end
convergence = mean(distances);
end