CAN YOU CHANGE MY NSGA2 CODE 1 STYLE OF WRITING TO BE MORE LIKE MY NSGA2 CODE 2? PLEASE DONT CHANGE THE OPERATOR. YOU JUST NEED TO CHANGE THE WRITING STYLE LIKE IN CODE 1 USE evaluate_objective CHANGE TO LIKE CODE 2 JUST USE FUNCTION ZDT1? DONT CHANGE THE OPERATOR IN THE FIRST CODE.

FIRST CODE:
% NSGA-II optimization
% Main script to run NSGA-II optimization
clc; 
clear; 
close all;

% Set problem parameters
NPop = 100;  % Population size
T = 500;     % Number of iterations
M = 2;
dim= 30;
lb = zeros(1, dim);  % Lower bound
ub = ones(1, dim);   % Upper bound

% Run NSGA-II
[paretoFront, paretoSet, Trajectories, fitness_history, position_history] = NSGA2(NPop, T, lb, ub, dim, M);

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

% 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] = NSGA2(NPop, T, lb, ub, dim, M)
% Non-dominated Sorting Genetic Algorithm II (NSGA-II)
% NSGA-II attempts to solve multi-objective problems of the following form:
%  min F(X) = [f1(X), f2(X), ..., fm(X)]  subject to: lb <= X <= ub
%                  
%  lb           - lower bounds on X
%  ub           - upper bounds on X
%  T            - number of iterations
%  NPop         - size of the population
%  M            - number of objectives
%  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, M);
position_history = zeros(NPop, T, dim);
Trajectories = zeros(dim, T);

% Initialize population
chromosome = initialize_variables(NPop, M, dim, lb, ub);

% Sort the initialized population
chromosome = non_domination_sort_mod(chromosome, M, dim);

% Start the evolution process
for i = 1:T
    % Selection
    pool = round(NPop/2);
    tour = 2;
    parent_chromosome = tournament_selection(chromosome, pool, tour);
    
    % Genetic operators (crossover and mutation)
    mu = 20;
    mum = 20;
    offspring_chromosome = genetic_operator(parent_chromosome, M, dim, mu, mum, lb, ub);
    
    % Merge populations
    [main_pop,~] = size(chromosome);
    [offspring_pop,~] = size(offspring_chromosome);
    intermediate_chromosome(1:main_pop,:) = chromosome;
    intermediate_chromosome(main_pop + 1:main_pop + offspring_pop,1:M+dim) = offspring_chromosome;
    
    % Non-dominated sorting
    intermediate_chromosome = non_domination_sort_mod(intermediate_chromosome, M, dim);
    
    % Selection for next generation
    chromosome = replace_chromosome(intermediate_chromosome, M, dim, NPop);
    
    % Record history
    fitness_history(:,i,:) = chromosome(:,dim+1:dim+M);
    position_history(:,i,:) = chromosome(:,1:dim);
    
    % Update Trajectories (using the best ranked solution)
    [~, best_idx] = min(chromosome(:,dim+M+1));
    Trajectories(:,i) = chromosome(best_idx,1:dim)';
    
    if ~mod(i,50)
        disp(['Iteration #', num2str(i), ' completed']);
    end
end

% Extract Pareto front and set
[~, unique_idx] = unique(chromosome(:,dim+1:dim+M), 'rows');
paretoFront = chromosome(unique_idx,dim+1:dim+M);
paretoSet = chromosome(unique_idx,1:dim);

end


function f = evaluate_objective(x, M, V)
    f = [];
    f(1) = x(1);
    sum = 0;
    for i = 2:V
        sum = sum + x(i);
    end
    g = 1 + 9 * sum / (V - 1);
    h = 1 - sqrt(f(1) / g);
    f(2) = g * h;


end

function f = initialize_variables(N, M, V, min_range, max_range)

min = min_range;
max = max_range;

K = M + V;

for i = 1 : N
    for j = 1 : V
        f(i,j) = min(j) + (max(j) - min(j))*rand(1);
    end
    f(i,V + 1: K) = evaluate_objective(f(i,:), M, V);
end
end



function f = non_domination_sort_mod(x, M, V)

%% function f = non_domination_sort_mod(x, M, V)
[N, m] = size(x);
clear m

% Initialize the front number to 1.
front = 1;

% There is nothing to this assignment, used only to manipulate easily in
% MATLAB.
F(front).f = [];
individual = [];

%% Non-Dominated sort. 
for i = 1 : N
    % Number of individuals that dominate this individual
    individual(i).n = 0;
    % Individuals which this individual dominate
    individual(i).p = [];
    for j = 1 : N
        dom_less = 0;
        dom_equal = 0;
        dom_more = 0;
        for k = 1 : M
            if (x(i,V + k) < x(j,V + k))
                dom_less = dom_less + 1;
            elseif (x(i,V + k) == x(j,V + k))
                dom_equal = dom_equal + 1;
            else
                dom_more = dom_more + 1;
            end
        end
        if dom_less == 0 && dom_equal ~= M
            individual(i).n = individual(i).n + 1;
        elseif dom_more == 0 && dom_equal ~= M
            individual(i).p = [individual(i).p j];
        end
    end   
    if individual(i).n == 0
        x(i,M + V + 1) = 1;
        F(front).f = [F(front).f i];
    end
end
% Find the subsequent fronts
while ~isempty(F(front).f)
   Q = [];
   for i = 1 : length(F(front).f)
       if ~isempty(individual(F(front).f(i)).p)
        	for j = 1 : length(individual(F(front).f(i)).p)
            	individual(individual(F(front).f(i)).p(j)).n = ...
                	individual(individual(F(front).f(i)).p(j)).n - 1;
        	   	if individual(individual(F(front).f(i)).p(j)).n == 0
               		x(individual(F(front).f(i)).p(j),M + V + 1) = ...
                        front + 1;
                    Q = [Q individual(F(front).f(i)).p(j)];
                end
            end
       end
   end
   front =  front + 1;
   F(front).f = Q;
end

[temp,index_of_fronts] = sort(x(:,M + V + 1));
for i = 1 : length(index_of_fronts)
    sorted_based_on_front(i,:) = x(index_of_fronts(i),:);
end
current_index = 0;

%% Crowding distance
% Find the crowding distance for each individual in each front
for front = 1 : (length(F) - 1)
%    objective = [];
    distance = 0;
    y = [];
    previous_index = current_index + 1;
    for i = 1 : length(F(front).f)
        y(i,:) = sorted_based_on_front(current_index + i,:);
    end
    current_index = current_index + i;
    % Sort each individual based on the objective
    sorted_based_on_objective = [];
    for i = 1 : M
        [sorted_based_on_objective, index_of_objectives] = ...
            sort(y(:,V + i));
        sorted_based_on_objective = [];
        for j = 1 : length(index_of_objectives)
            sorted_based_on_objective(j,:) = y(index_of_objectives(j),:);
        end
        f_max = ...
            sorted_based_on_objective(length(index_of_objectives), V + i);
        f_min = sorted_based_on_objective(1, V + i);
        y(index_of_objectives(length(index_of_objectives)),M + V + 1 + i)...
            = Inf;
        y(index_of_objectives(1),M + V + 1 + i) = Inf;
         for j = 2 : length(index_of_objectives) - 1
            next_obj  = sorted_based_on_objective(j + 1,V + i);
            previous_obj  = sorted_based_on_objective(j - 1,V + i);
            if (f_max - f_min == 0)
                y(index_of_objectives(j),M + V + 1 + i) = Inf;
            else
                y(index_of_objectives(j),M + V + 1 + i) = ...
                     (next_obj - previous_obj)/(f_max - f_min);
            end
         end
    end
    distance = [];
    distance(:,1) = zeros(length(F(front).f),1);
    for i = 1 : M
        distance(:,1) = distance(:,1) + y(:,M + V + 1 + i);
    end
    y(:,M + V + 2) = distance;
    y = y(:,1 : M + V + 2);
    z(previous_index:current_index,:) = y;
end
f = z();
end

function f  = replace_chromosome(intermediate_chromosome, M, V,pop)

%% function f  = replace_chromosome(intermediate_chromosome,pro,pop)

[N, m] = size(intermediate_chromosome);

% Get the index for the population sort based on the rank
[temp,index] = sort(intermediate_chromosome(:,M + V + 1));

clear temp m

% Now sort the individuals based on the index
for i = 1 : N
    sorted_chromosome(i,:) = intermediate_chromosome(index(i),:);
end

% Find the maximum rank in the current population
max_rank = max(intermediate_chromosome(:,M + V + 1));

% Start adding each front based on rank and crowing distance until the
% whole population is filled.
previous_index = 0;
for i = 1 : max_rank
    % Get the index for current rank i.e the last the last element in the
    % sorted_chromosome with rank i. 
    current_index = max(find(sorted_chromosome(:,M + V + 1) == i));
    % Check to see if the population is filled if all the individuals with
    % rank i is added to the population. 
    if current_index > pop
        % If so then find the number of individuals with in with current
        % rank i.
        remaining = pop - previous_index;
        % Get information about the individuals in the current rank i.
        temp_pop = ...
            sorted_chromosome(previous_index + 1 : current_index, :);
        % Sort the individuals with rank i in the descending order based on
        % the crowding distance.
        [temp_sort,temp_sort_index] = ...
            sort(temp_pop(:, M + V + 2),'descend');
        % Start filling individuals into the population in descending order
        % until the population is filled.
        for j = 1 : remaining
            f(previous_index + j,:) = temp_pop(temp_sort_index(j),:);
        end
        return;
    elseif current_index < pop
        % Add all the individuals with rank i into the population.
        f(previous_index + 1 : current_index, :) = ...
            sorted_chromosome(previous_index + 1 : current_index, :);
    else
        % Add all the individuals with rank i into the population.
        f(previous_index + 1 : current_index, :) = ...
            sorted_chromosome(previous_index + 1 : current_index, :);
        return;
    end
    % Get the index for the last added individual.
    previous_index = current_index;
end
end

function f = tournament_selection(chromosome, pool_size, tour_size)

%% function tournament_selection(chromosome, pool_size, tour_size) 

%% Tournament selection process

% Get the size of chromosome. The number of chromosome is not important
% while the number of elements in chromosome are important.
[pop, variables] = size(chromosome);
% The peunltimate element contains the information about rank.
rank = variables - 1;
% The last element contains information about crowding distance.
distance = variables;

% Until the mating pool is filled, perform tournament selection
for i = 1 : pool_size
    % Select n individuals at random, where n = tour_size
    for j = 1 : tour_size
        % Select an individual at random
        candidate(j) = round(pop*rand(1));
        % Make sure that the array starts from one. 
        if candidate(j) == 0
            candidate(j) = 1;
        end
        if j > 1
            % Make sure that same candidate is not choosen.
            while ~isempty(find(candidate(1 : j - 1) == candidate(j)))
                candidate(j) = round(pop*rand(1));
                if candidate(j) == 0
                    candidate(j) = 1;
                end
            end
        end
    end
    % Collect information about the selected candidates.
    for j = 1 : tour_size
        c_obj_rank(j) = chromosome(candidate(j),rank);
        c_obj_distance(j) = chromosome(candidate(j),distance);
    end
    % Find the candidate with the least rank
    min_candidate = ...
        find(c_obj_rank == min(c_obj_rank));
    % If more than one candiate have the least rank then find the candidate
    % within that group having the maximum crowding distance.
    if length(min_candidate) ~= 1
        max_candidate = ...
        find(c_obj_distance(min_candidate) == max(c_obj_distance(min_candidate)));
        % If a few individuals have the least rank and have maximum crowding
        % distance, select only one individual (not at random). 
        if length(max_candidate) ~= 1
            max_candidate = max_candidate(1);
        end
        % Add the selected individual to the mating pool
        f(i,:) = chromosome(candidate(min_candidate(max_candidate)),:);
    else
        % Add the selected individual to the mating pool
        f(i,:) = chromosome(candidate(min_candidate(1)),:);
    end
end
end

function f  = genetic_operator(parent_chromosome, M, V, mu, mum, l_limit, u_limit)

[N,m] = size(parent_chromosome);

clear m
p = 1;
% Flags used to set if crossover and mutation were actually performed. 
was_crossover = 0;
was_mutation = 0;


for i = 1 : N
    % With 90 % probability perform crossover
    if rand(1) < 0.9
        % Initialize the children to be null vector.
        child_1 = [];
        child_2 = [];
        % Select the first parent
        parent_1 = round(N*rand(1));
        if parent_1 < 1
            parent_1 = 1;
        end
        % Select the second parent
        parent_2 = round(N*rand(1));
        if parent_2 < 1
            parent_2 = 1;
        end
        % Make sure both the parents are not the same. 
        while isequal(parent_chromosome(parent_1,:),parent_chromosome(parent_2,:))
            parent_2 = round(N*rand(1));
            if parent_2 < 1
                parent_2 = 1;
            end
        end
        % Get the chromosome information for each randomnly selected
        % parents
        parent_1 = parent_chromosome(parent_1,:);
        parent_2 = parent_chromosome(parent_2,:);
        % Perform corssover for each decision variable in the chromosome.
        for j = 1 : V
            % SBX (Simulated Binary Crossover).
            % For more information about SBX refer the enclosed pdf file.
            % Generate a random number
            u(j) = rand(1);
            if u(j) <= 0.5
                bq(j) = (2*u(j))^(1/(mu+1));
            else
                bq(j) = (1/(2*(1 - u(j))))^(1/(mu+1));
            end
            % Generate the jth element of first child
            child_1(j) = ...
                0.5*(((1 + bq(j))*parent_1(j)) + (1 - bq(j))*parent_2(j));
            % Generate the jth element of second child
            child_2(j) = ...
                0.5*(((1 - bq(j))*parent_1(j)) + (1 + bq(j))*parent_2(j));
            % Make sure that the generated element is within the specified
            % decision space else set it to the appropriate extrema.
            if child_1(j) > u_limit(j)
                child_1(j) = u_limit(j);
            elseif child_1(j) < l_limit(j)
                child_1(j) = l_limit(j);
            end
            if child_2(j) > u_limit(j)
                child_2(j) = u_limit(j);
            elseif child_2(j) < l_limit(j)
                child_2(j) = l_limit(j);
            end
        end
        % Evaluate the objective function for the offsprings and as before
        % concatenate the offspring chromosome with objective value.
        child_1(:,V + 1: M + V) = evaluate_objective(child_1, M, V);
        child_2(:,V + 1: M + V) = evaluate_objective(child_2, M, V);
        % Set the crossover flag. When crossover is performed two children
        % are generate, while when mutation is performed only only child is
        % generated.
        was_crossover = 1;
        was_mutation = 0;
    % With 10 % probability perform mutation. Mutation is based on
    % polynomial mutation. 
    else
        % Select at random the parent.
        parent_3 = round(N*rand(1));
        if parent_3 < 1
            parent_3 = 1;
        end
        % Get the chromosome information for the randomnly selected parent.
        child_3 = parent_chromosome(parent_3,:);
        % Perform mutation on eact element of the selected parent.
        for j = 1 : V
           r(j) = rand(1);
           if r(j) < 0.5
               delta(j) = (2*r(j))^(1/(mum+1)) - 1;
           else
               delta(j) = 1 - (2*(1 - r(j)))^(1/(mum+1));
           end
           % Generate the corresponding child element.
           child_3(j) = child_3(j) + delta(j);
           % Make sure that the generated element is within the decision
           % space.
           if child_3(j) > u_limit(j)
               child_3(j) = u_limit(j);
           elseif child_3(j) < l_limit(j)
               child_3(j) = l_limit(j);
           end
        end
        % Evaluate the objective function for the offspring and as before
        % concatenate the offspring chromosome with objective value.    
        child_3(:,V + 1: M + V) = evaluate_objective(child_3, M, V);
        % Set the mutation flag
        was_mutation = 1;
        was_crossover = 0;
    end
    % Keep proper count and appropriately fill the child variable with all
    % the generated children for the particular generation.
    if was_crossover
        child(p,:) = child_1;
        child(p+1,:) = child_2;
        was_cossover = 0;
        p = p + 2;
    elseif was_mutation
        child(p,:) = child_3(1,1 : M + V);
        was_mutation = 0;
        p = p + 1;
    end
end
f = child;
end


SECOND CODE:
clc;
clear;
close all;

% Set problem parameters
NPop = 150;  % 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
num_runs = 10;  % Number of runs for statistical evaluation

% Storage for metrics
GD_values = zeros(num_runs, 1);
IGD_values = zeros(num_runs, 1);
Spacing_values = zeros(num_runs, 1);
Spread_values = zeros(num_runs, 1);
DU_values = zeros(num_runs, 1);
MS_values = zeros(num_runs, 1);  % Added MS storage

% store best pareto
paretoFronts = cell(num_runs, 1);
paretoSets = cell(num_runs, 1);

%% Run NSGA-II Multiple Times
% Run NSGA-II Multiple Times
for run = 1:num_runs
    fprintf('Run #%d/%d\n', run, num_runs);
    
    % Run NSGA-II
    [paretoFront, paretoSet] = NSGA2(NPop, T, lb, ub, dim, @zdt1);
    
    % Store Pareto front and set
    paretoFronts{run} = paretoFront;
    paretoSets{run} = paretoSet;
    
    % Get True Pareto Front
    true_pareto = get_true_pareto();

    % Inside the main script
    [GD_values(run), IGD_values(run), Spacing_values(run), Spread_values(run), ...
        DU_values(run), MS_values(run)] = calculate_metrics(paretoFront, true_pareto);
end
metrics = {GD_values, IGD_values, Spacing_values, Spread_values, DU_values, MS_values};
metric_names = {'GD', 'IGD', 'Spacing', 'Spread', 'DU', 'MS'};

% Find the best run based on the chosen metric (e.g., IGD)
[~, best_run_idx] = min(IGD_values);
best_paretoFront = paretoFronts{best_run_idx};
best_paretoSet = paretoSets{best_run_idx};

% 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 the best Pareto Front
figure;
plot(best_paretoFront(:,1), best_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 ZDT1 (NGSAII)');
legend('Obtained Pareto Front', 'True Pareto Front');
grid on;
hold off;

% Display metrics for the best run
fprintf('\nBest Run #%d Metrics:\n', best_run_idx);
fprintf('GD    = %.4f\n', GD_values(best_run_idx));
fprintf('IGD   = %.4f\n', IGD_values(best_run_idx));
fprintf('MS    = %.4f\n', MS_values(best_run_idx));

% 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 ZDT1 (NGSAII)');
legend('Obtained Pareto Front', 'True Pareto Front');
grid on;
hold off;

% Plot the results
figure;
plot(paretoFront(:,1), paretoFront(:,2), 'ro', 'DisplayName', 'NSGA-II');
xlabel('f1');
ylabel('f2');
title('Pareto Front for ZDT1 (NGSAII)');
legend('show');
grid on;

function [paretoFront, paretoSet] = NSGA2(NPop, T, lb, ub, dim, fitnessFcn)
    % Initialize population
    pop = initialization(NPop, dim, ub, lb);
    fitness = evaluatePopulation(pop, fitnessFcn);

    for gen = 1:T
        % Non-dominated sorting
        [fronts, ranks] = nonDominatedSort(fitness);
        
        % Create mating pool using tournament selection
        matingPool = tournamentSelection(pop, ranks, NPop);
        
        % Generate offspring using SBX crossover and polynomial mutation
        offspring = crossover(matingPool, 0.9, dim);
        offspring = mutation(offspring, 1/dim, dim);
        
        % Evaluate offspring
        offspringFitness = evaluatePopulation(offspring, fitnessFcn);
        
        % Combine population and offspring
        combinedPop = [pop; offspring];
        combinedFitness = [fitness; offspringFitness];
        
        % Non-dominated sorting on combined population
        [fronts, ~] = nonDominatedSort(combinedFitness);
        
        % Create new population
        [pop, fitness] = selectNewPopulation(combinedPop, combinedFitness, fronts, NPop);
        
        % Display generation progress
        if mod(gen, 10) == 0
            fprintf('Generation #%d/%d completed.\n', gen, T);
        end
    end

    % Return last population as Pareto front
    paretoFront = fitness;
    paretoSet = pop;
end

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

function fitness = evaluatePopulation(pop, fitnessFcn)
    numIndividuals = size(pop, 1);
    numObjectives = 2; % This should be set based on the problem
    fitness = zeros(numIndividuals, numObjectives);
    for i = 1:numIndividuals
        fitness(i, :) = fitnessFcn(pop(i, :));
    end
end

function [fronts, rank] = nonDominatedSort(fitness)
    numIndividuals = size(fitness, 1);
    dominationCount = zeros(numIndividuals, 1);
    dominatedSet = cell(numIndividuals, 1);
    rank = zeros(numIndividuals, 1);
    fronts = cell(1);

    for i = 1:numIndividuals
        for j = 1:numIndividuals
            if i ~= j
                if dominates(fitness(i, :), fitness(j, :))
                    dominatedSet{i}(end+1) = j;
                elseif dominates(fitness(j, :), fitness(i, :))
                    dominationCount(i) = dominationCount(i) + 1;
                end
            end
        end
        
        if dominationCount(i) == 0
            rank(i) = 1;
            fronts{1}(end+1) = i;
        end
    end
    
    currentFront = 1;
    while ~isempty(fronts{currentFront})
        nextFront = [];
        for i = fronts{currentFront}
            for j = dominatedSet{i}
                dominationCount(j) = dominationCount(j) - 1;
                if dominationCount(j) == 0
                    rank(j) = currentFront + 1;
                    nextFront(end+1) = j;
                end
            end
        end
        currentFront = currentFront + 1;
        fronts{currentFront} = nextFront;
    end
    fronts = fronts(1:currentFront-1);
end

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

function selectedPop = tournamentSelection(pop, ranks, numSelected)
    selectedPop = zeros(numSelected, size(pop, 2));
    for i = 1:numSelected
        candidates = randperm(size(pop, 1), 2);
        if ranks(candidates(1)) < ranks(candidates(2))
            selectedPop(i, :) = pop(candidates(1), :);
        else
            selectedPop(i, :) = pop(candidates(2), :);
        end
    end
end

function offspring = crossover(pop, crossoverProb, numVars)
    etaC = 20; % Distribution index for crossover
    numIndividuals = size(pop, 1);
    offspring = zeros(numIndividuals, numVars);
    
    for i = 1:2:numIndividuals
        if rand < crossoverProb
            parent1 = pop(i, :);
            parent2 = pop(i+1, :);
            
            child1 = zeros(1, numVars);
            child2 = zeros(1, numVars);
            
            for j = 1:numVars
                u = rand;
                if u <= 0.5
                    beta = (2*u)^(1/(etaC+1));
                else
                    beta = (1/(2*(1-u)))^(1/(etaC+1));
                end
                
                child1(j) = 0.5 * ((1+beta) * parent1(j) + (1-beta) * parent2(j));
                child2(j) = 0.5 * ((1-beta) * parent1(j) + (1+beta) * parent2(j));
                
                % Handle bounds
                child1(j) = min(max(child1(j), 0), 1);
                child2(j) = min(max(child2(j), 0), 1);
            end
            
            offspring(i, :) = child1;
            offspring(i+1, :) = child2;
        else
            offspring(i, :) = pop(i, :);
            offspring(i+1, :) = pop(i+1, :);
        end
    end
end

function mutatedPop = mutation(pop, mutationProb, numVars)
    etaM = 20; % Mutation distribution index
    mutatedPop = pop;
    numIndividuals = size(pop, 1);
    
    for i = 1:numIndividuals
        for j = 1:numVars
            if rand < mutationProb
                r = rand;
                if r < 0.5
                    delta = (2*r)^(1/(etaM + 1)) - 1;
                else
                    delta = 1 - (2*(1-r))^(1/(etaM + 1));
                end
                mutatedPop(i, j) = pop(i, j) + delta;
                mutatedPop(i, j) = min(max(mutatedPop(i, j), 0), 1);
            end
        end
    end
end

function [pop, fitness] = selectNewPopulation(combinedPop, combinedFitness, fronts, NPop)
    pop = [];
    fitness = [];
    for i = 1:length(fronts)
        if size(pop, 1) + length(fronts{i}) > NPop
            break;
        end
        pop = [pop; combinedPop(fronts{i}, :)];
        fitness = [fitness; combinedFitness(fronts{i}, :)];
    end
    
    remaining = NPop - size(pop, 1);
    if remaining > 0
        selectedFront = fronts{i};
        distances = crowdingDistance(combinedFitness(selectedFront, :));
        [~, sortedIndices] = sort(distances, 'descend');
        selected = selectedFront(sortedIndices(1:remaining));
        pop = [pop; combinedPop(selected, :)];
        fitness = [fitness; combinedFitness(selected, :)];
    end
end

function distances = crowdingDistance(fitness)
    numIndividuals = size(fitness, 1);
    numObjectives = size(fitness, 2);
    distances = zeros(numIndividuals, 1);
    
    for m = 1:numObjectives
        [sortedValues, sortedIdx] = sort(fitness(:, m));
        distances(sortedIdx(1)) = inf;
        distances(sortedIdx(end)) = inf;
        for i = 2:numIndividuals-1
            range = sortedValues(end) - sortedValues(1);
            if range > 0
                distances(sortedIdx(i)) = distances(sortedIdx(i)) + ...
                    (sortedValues(i+1) - sortedValues(i-1)) / range;
            end
        end
    end
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

function true_pareto = get_true_pareto()
    num_points = 100;
    f1 = linspace(0, 1, num_points)';
    f2 = 1 - sqrt(f1);
    true_pareto = [f1, f2];
end