can you change the use of zdt1 and zdt in my code to be change to Multi-Objective Cost Function. 
Data:
% Data from UTHM
dataSH = [
    1.0000   60.0000   15.5600  711.0000  610.0000
    2.0000   65.0000   18.3300  717.0000  670.0000
    3.0000   67.0000   19.4400  718.0000  650.0000
    4.0000   66.0000   18.8900  717.0000  640.0000
    5.0000   67.0000   19.4400  700.0000  600.0000
    6.0000   63.0000   17.2200  713.0000  620.0000
    7.0000   62.0000   16.6700  724.0000  670.0000
    8.0000   81.0000   27.2200  710.0000  920.0000
    9.0000   80.0000   26.6700  914.0000  920.0000
   10.0000   81.0000   27.2200  900.0000  960.0000
   11.0000   83.0000   28.3300  897.0000  900.0000
   12.0000   80.0000   26.6700  895.0000  930.0000
   13.0000   79.0000   26.1100  904.0000  950.0000
   14.0000   82.0000   27.7800  903.0000  960.0000
   15.0000   81.0000   27.2200  906.0000  970.0000
   16.0000   80.0000   26.6700  898.0000  925.0000
   17.0000   66.0000   18.8900  728.0000  610.0000
   18.0000   63.0000   17.2200  703.0000  670.0000
   19.0000   66.0000   18.8900  714.0000  630.0000
   20.0000   66.0000   18.8900  721.0000  645.0000
   21.0000   66.0000   18.8900  722.0000  650.0000
   22.0000   67.0000   19.4400  703.0000  600.0000
   23.0000   66.0000   18.8900  701.0000  670.0000
   24.0000   81.0000   27.2200  728.0000  980.0000
   25.0000   81.0000   27.2200  891.0000  930.0000
   26.0000   82.0000   27.7800  906.0000  948.0000
   27.0000   80.0000   26.6700  901.0000  965.0000
   28.0000   81.0000   27.2200  901.0000  916.0000
   29.0000   79.0000   26.1100  915.0000  900.0000
   30.0000   83.0000   28.3300  914.0000  960.0000
   31.0000   79.0000   26.1100  890.0000  970.0000
   32.0000   81.0000   27.2200  912.0000  930.0000
   33.0000   65.0000   18.3300  717.0000  610.0000
   34.0000   67.0000   19.4400  714.0000  620.0000
   35.0000   65.0000   18.3300  724.0000  670.0000
   36.0000   67.0000   19.4400  726.0000  680.0000
   37.0000   66.0000   18.8900  723.0000  650.0000
   38.0000   64.0000   17.7800  726.0000  660.0000
   39.0000   65.0000   18.3300  708.0000  640.0000
   40.0000   80.0000   26.6700  712.0000  900.0000
   41.0000   83.0000   28.3300  917.0000  910.0000
   42.0000   79.0000   26.1100  905.0000  920.0000
   43.0000   82.0000   27.7800  893.0000  980.0000
   44.0000   80.0000   26.6700  920.0000  940.0000
   45.0000   81.0000   27.2200  913.0000  950.0000
   46.0000   81.0000   27.2200  905.0000  940.0000
   47.0000   82.0000   27.7800  896.0000  970.0000
   48.0000   80.0000   26.6700  898.0000  980.0000
];
Multi-Objective Cost Function:
% ---------- Multi-Objective Cost Function ----------
function objectives = multi_objective_cost_function(X, dataSH)
    % MULTI_OBJECTIVE_COST_FUNCTION Evaluates user comfort and energy consumption.
    %
    % This function is designed for multi-objective optimization with two objectives:
    %   1. Maximize User Comfort (to be minimized as -CI_avg)
    %   2. Minimize Energy Consumption (EC_avg)
    %
    % Inputs:
    %   X      - A vector of decision variables for 48 periods, structured as:
    %            [T1, L1, IAQ1, T2, L2, IAQ2, ..., T48, L48, IAQ48]
    %   dataSH - A 48x5 matrix where:
    %            dataSH(:,3) = Temperature coefficients (Tc)
    %            dataSH(:,4) = Light coefficients (Lc)
    %            dataSH(:,5) = IAQ coefficients (IAQc)
    %
    % Output:
    %   objectives - A 2-element vector:
    %                [ -Average Comfort Index, Average Energy Consumption ]
    
    % Ensure that X has the correct length
    if length(X) ~= 144
        error('Decision vector X must have 144 elements (3 variables per 48 periods).');
    end
    
    % Weight coefficients for comfort components
    p1 = 1/3; % Temperature Weight
    p2 = 1/3; % Light Weight
    p3 = 1/3; % IAQ Weight
    
    % Desired setpoints for comfort
    Td = 22.78; % Desired Temperature (°C)
    Ld = 800;   % Desired Light Level (lux)
    Id = 800;   % Desired IAQ Level
    
    % Energy consumption weights
    pT = 9; % Temperature Energy Consumption Weight
    pL = 1; % Light Energy Consumption Weight
    pI = 1; % IAQ Energy Consumption Weight
    
    % Extract coefficients from dataSH for all 48 periods
    Tc = dataSH(:,3)'; % Temperature coefficients (1x48 vector)
    Lc = dataSH(:,4)'; % Light coefficients (1x48 vector)
    Ic = dataSH(:,5)'; % IAQ coefficients (1x48 vector)
    
    % Initialize vectors to store metrics for each period
    CI = zeros(48,1);    % Comfort Index for each period
    EC = zeros(48,1);    % Energy Consumption for each period
    
    % Loop through each of the 48 periods to compute metrics
    for k = 1:48
        % Extract decision variables for period k
        Tk = X(3*k-2); % Temperature setting
        Lk = X(3*k-1); % Light setting
        Ik = X(3*k);   % IAQ setting
        
        % --- Calculate Comfort Index Components ---
        Ti = 1 - ((Tk - Td) / Tc(k))^2;
        Li = 1 - ((Lk - Ld) / Lc(k))^2;
        Ii = 1 - ((Ik - Id) / Ic(k))^2;
        CI(k) = p1 * Ti + p2 * Li + p3 * Ii;
        
        % --- Calculate Energy Consumption ---
        ET = pT * abs(Tk - Tc(k));
        EL = pL * abs(Lk - Lc(k));
        EI = pI * abs(Ik - Ic(k));
        EC(k) = ET + EL + EI;
    end
    
    % --- Aggregate Metrics Over All Periods ---
    CI_avg = mean(CI);   % Average Comfort Index
    EC_avg = mean(EC);   % Average Energy Consumption
    
    % --- Prepare Objectives for Multi-Objective Optimization ---
    % Objective 1: Maximize User Comfort -> Minimize (-CI_avg)
    % Objective 2: Minimize Energy Consumption -> Minimize (EC_avg)
    objectives = [-CI_avg, EC_avg];
end


My code:
% Main script to run MOEMA, MOEA/D, and MOPSO optimization
clc; 
clear; 
close all;

% Set problem parameters
NPop = 100;  % 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

% Run MOEMA
[paretoFront_MOEMA, ~, ~, ~, ~] = MOEMA(NPop, T, lb, ub, dim, @zdt1);

% Run MOEA/D
nVar = dim;
VarMin = 0;
VarMax = 1;
nPop_MOEAD = 200;
nArchive = 50;
MaxIt = T;
T_MOEAD = 15;
CostFunction = @(x) ZDT(x);
crossover_params.gamma = 0.5;
crossover_params.VarMin = VarMin;
crossover_params.VarMax = VarMax;
[paretoFront_MOEAD, ~, ~, ~] = MOEAD(nVar, VarMin, VarMax, nPop_MOEAD, nArchive, MaxIt, T_MOEAD, CostFunction, crossover_params);

% Run MOPSO
params.Np = 200;        % Population size
params.Nr = 100;        % Repository size
params.maxgen = T;      % Maximum number of generations
params.W = 0.4;         % Inertia weight
params.C1 = 2;          % Individual confidence factor
params.C2 = 2;          % Swarm confidence factor
params.ngrid = 20;      % Number of grids in each dimension
params.maxvel = 5;      % Maximum velocity in percentage
params.u_mut = 0.5;     % Uniform mutation percentage

MultiObj.nVar = dim;
MultiObj.var_min = lb';
MultiObj.var_max = ub';
MultiObj.fun = @(POS) cell2mat(arrayfun(@(i) zdt1(POS(i,:)), 1:size(POS, 1), 'UniformOutput', false)');

REP_MOPSO = MOPSO(params, MultiObj);
paretoFront_MOPSO = REP_MOPSO.pos_fit;

% Plot the Pareto fronts
figure;
plot(paretoFront_MOEMA(:,1), paretoFront_MOEMA(:,2), 'ro', 'DisplayName', 'MOEMA');
hold on;
plot(paretoFront_MOEAD(1,:), paretoFront_MOEAD(2,:), 'bx', 'DisplayName', 'MOEA/D');
plot(paretoFront_MOPSO(:,1), paretoFront_MOPSO(:,2), 'g+', 'DisplayName', 'MOPSO');
xlabel('f1');
ylabel('f2');
title('Pareto Fronts for ZDT1 using MOEMA, MOEA/D, and MOPSO');
legend;

% Generate true ZDT1 Pareto front
f1 = 0:0.01:1;
g = 1;
h = 1 - sqrt(f1 ./ g);
f2 = g .* h;
plot(f1, f2, 'k-', 'DisplayName', 'True Pareto Front');
legend;

% Create a table for comparison
MOEMA_solutions = size(paretoFront_MOEMA, 1);
MOEAD_solutions = size(paretoFront_MOEAD, 2);
MOPSO_solutions = size(paretoFront_MOPSO, 1);
MOEMA_min_f1 = min(paretoFront_MOEMA(:,1));
MOEAD_min_f1 = min(paretoFront_MOEAD(1,:));
MOPSO_min_f1 = min(paretoFront_MOPSO(:,1));
MOEMA_min_f2 = min(paretoFront_MOEMA(:,2));
MOEAD_min_f2 = min(paretoFront_MOEAD(2,:));
MOPSO_min_f2 = min(paretoFront_MOPSO(:,2));

table = array2table([MOEMA_solutions; MOEAD_solutions; MOPSO_solutions; MOEMA_min_f1; MOEAD_min_f1; MOPSO_min_f1; MOEMA_min_f2; MOEAD_min_f2; MOPSO_min_f2]);
table.Properties.VariableNames = {'Value'};
table.Properties.RowNames = {'MOEMA Solutions', 'MOEA/D Solutions', 'MOPSO Solutions', 'MOEMA Min f1', 'MOEA/D Min f1', 'MOPSO Min f1', 'MOEMA Min f2', 'MOEA/D Min f2', 'MOPSO Min f2'};
disp(table);

% Main script continuation
% Calculate performance metrics
truePF = [f1', f2'];

[GD_MOEMA, IGD_MOEMA, Spread_MOEMA, HV_MOEMA] = calculateMetrics(paretoFront_MOEMA, truePF);
[GD_MOEAD, IGD_MOEAD, Spread_MOEAD, HV_MOEAD] = calculateMetrics(paretoFront_MOEAD', truePF);
[GD_MOPSO, IGD_MOPSO, Spread_MOPSO, HV_MOPSO] = calculateMetrics(paretoFront_MOPSO, truePF);

% Display performance metrics
disp('Performance Metrics:');
metrics_table = array2table([GD_MOEMA, IGD_MOEMA, Spread_MOEMA, HV_MOEMA;
                             GD_MOEAD, IGD_MOEAD, Spread_MOEAD, HV_MOEAD;
                             GD_MOPSO, IGD_MOPSO, Spread_MOPSO, HV_MOPSO]);
metrics_table.Properties.VariableNames = {'GD', 'IGD', 'Spread', 'HV'};
metrics_table.Properties.RowNames = {'MOEMA', 'MOEA/D', 'MOPSO'};
disp(metrics_table);

% Plotting
figure;
subplot(2,2,1);
plot(paretoFront_MOEMA(:,1), paretoFront_MOEMA(:,2), 'ro', 'DisplayName', 'MOEMA');
hold on;
plot(paretoFront_MOEAD(1,:), paretoFront_MOEAD(2,:), 'bx', 'DisplayName', 'MOEA/D');
plot(paretoFront_MOPSO(:,1), paretoFront_MOPSO(:,2), 'g+', 'DisplayName', 'MOPSO');
plot(f1, f2, 'k-', 'DisplayName', 'True Pareto Front');
xlabel('f1'); ylabel('f2');
title('Pareto Fronts Comparison');
legend('Location', 'southwest');

subplot(2,2,2);
bar([GD_MOEMA, GD_MOEAD, GD_MOPSO]);
title('Generation Distance (GD)');
set(gca, 'XTickLabel', {'MOEMA', 'MOEA/D', 'MOPSO'});

subplot(2,2,3);
bar([IGD_MOEMA, IGD_MOEAD, IGD_MOPSO]);
title('Inverted Generation Distance (IGD)');
set(gca, 'XTickLabel', {'MOEMA', 'MOEA/D', 'MOPSO'});

subplot(2,2,4);
bar([HV_MOEMA, HV_MOEAD, HV_MOPSO]);
title('Hypervolume (HV)');
set(gca, 'XTickLabel', {'MOEMA', 'MOEA/D', 'MOPSO'});

% Save results
save('optimization_results.mat', 'paretoFront_MOEMA', 'paretoFront_MOEAD', 'paretoFront_MOPSO', ...
     'GD_MOEMA', 'IGD_MOEMA', 'Spread_MOEMA', 'HV_MOEMA', ...
     'GD_MOEAD', 'IGD_MOEAD', 'Spread_MOEAD', 'HV_MOEAD', ...
     'GD_MOPSO', 'IGD_MOPSO', 'Spread_MOPSO', 'HV_MOPSO');

disp('Optimization completed. Results saved to optimization_results.mat');

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, 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 result = dominates(obj1, obj2)
    result = all(obj1 <= obj2) && any(obj1 < obj2);
end


%% Main MOEA/D Function
function [paretoFront, paretoSet, fitness_history, position_history] = MOEAD(nVar, VarMin, VarMax, nPop, nArchive, MaxIt, T, CostFunction, crossover_params)

    % Problem Definition
    VarSize = [nVar, 1];   % Decision Variables Matrix Size
    nObj = numel(CostFunction(unifrnd(VarMin, VarMax, VarSize)));  % Number of Objectives

    % Create Sub-problems
    sp = CreateSubProblems(nObj, nPop, T);

    % Empty Individual Template
    empty_individual.Position = [];
    empty_individual.Cost = [];
    empty_individual.g = [];
    empty_individual.IsDominated = [];

    % Initialize Goal Point
    z = zeros(nObj, 1);

    % Create Initial Population
    pop = repmat(empty_individual, nPop, 1);
    for i = 1:nPop
        pop(i).Position = unifrnd(VarMin, VarMax, VarSize);
        pop(i).Cost = CostFunction(pop(i).Position);
        z = min(z, pop(i).Cost);  % Update goal point
    end

    for i = 1:nPop
        pop(i).g = DecomposedCost(pop(i), z, sp(i).lambda);
    end

    % Determine Population Domination Status
    pop = DetermineDomination(pop);

    % Initialize Estimated Pareto Front (EP)
    EP = pop(~[pop.IsDominated]);

    % History variables for fitness and positions
    fitness_history = zeros(nObj, MaxIt);
    position_history = cell(nPop, MaxIt);

    % Main Loop
    for it = 1:MaxIt
        for i = 1:nPop
            % Reproduction (Crossover)
            K = randsample(T, 2);
            
            j1 = sp(i).Neighbors(K(1));
            p1 = pop(j1);
            
            j2 = sp(i).Neighbors(K(2));
            p2 = pop(j2);
            
            y = empty_individual;
            y.Position = Crossover(p1.Position, p2.Position, crossover_params);
            
            y.Cost = CostFunction(y.Position);
            
            z = min(z, y.Cost);  % Update goal point
            
            % Update neighboring solutions
            for j = sp(i).Neighbors
                y.g = DecomposedCost(y, z, sp(j).lambda);
                if y.g <= pop(j).g
                    pop(j) = y;
                end
            end
        end
        
        % Determine Population Domination Status
        pop = DetermineDomination(pop);
        
        ndpop = pop(~[pop.IsDominated]);
        
        EP = [EP; ndpop]; %#ok<AGROW>
        EP = DetermineDomination(EP);
        EP = EP(~[EP.IsDominated]);
        
        if numel(EP) > nArchive
            Extra = numel(EP) - nArchive;
            ToBeDeleted = randsample(numel(EP), Extra);
            EP(ToBeDeleted) = [];
        end

        % Store fitness data and positions for this iteration
        fitness_history(:, it) = mean([pop.Cost], 2);  % Store average fitness
        for i = 1:nPop
            position_history{i, it} = pop(i).Position;
        end
        
        % Display Iteration Information
        disp(['Iteration ' num2str(it) ': Number of Pareto Solutions = ' num2str(numel(EP))]);
    end

    % Prepare Results
    paretoFront = [EP.Cost];
    paretoSet = [EP.Position];

end

%% Helper Functions

% Plot Pareto Front
function PlotCosts(paretoFront)
    plot(paretoFront(1,:), paretoFront(2,:), 'x');
    xlabel('1^{st} Objective');
    ylabel('2^{nd} Objective');
    grid on;
end

% Check if one solution dominates another
function b = Dominates(x, y)
    if isfield(x, 'Cost')
        x = x.Cost;
    end
    if isfield(y, 'Cost')
        y = y.Cost;
    end
    b = all(x <= y) && any(x < y);
end

% Determine Domination in the Population
function pop = DetermineDomination(pop)
    nPop = numel(pop);
    for i = 1:nPop
        pop(i).IsDominated = false;
    end
    for i = 1:nPop
        for j = i+1:nPop
            if Dominates(pop(i), pop(j))
                pop(j).IsDominated = true;
            elseif Dominates(pop(j), pop(i))
                pop(i).IsDominated = true;
            end
        end
    end
end

% Decomposed Cost Calculation using Tchebycheff Approach
function g = DecomposedCost(individual, z, lambda)
    if isfield(individual, 'Cost')
        fx = individual.Cost;
    else
        fx = individual;
    end
    g = max(lambda .* abs(fx - z));
end

% Crossover Operation
function y = Crossover(x1, x2, params)
    gamma = params.gamma;
    VarMin = params.VarMin;
    VarMax = params.VarMax;
    alpha = unifrnd(-gamma, 1+gamma, size(x1));
    y = alpha .* x1 + (1 - alpha) .* x2;
    y = min(max(y, VarMin), VarMax);
end

% Create Sub-problems
function sp = CreateSubProblems(nObj, nPop, T)
    empty_sp.lambda = [];
    empty_sp.Neighbors = [];
    sp = repmat(empty_sp, nPop, 1);

    for i = 1:nPop
        lambda = rand(nObj, 1);
        lambda = lambda / norm(lambda);
        sp(i).lambda = lambda;
    end

    LAMBDA = [sp.lambda]';
    D = pdist2(LAMBDA, LAMBDA);  % Distance between lambda vectors

    for i = 1:nPop
        [~, SO] = sort(D(i, :));
        sp(i).Neighbors = SO(1:T);
    end
end

% MOPSO Function
function REP = MOPSO(params, MultiObj)
    Np = params.Np;
    Nr = params.Nr;
    maxgen = params.maxgen;
    W = params.W;
    C1 = params.C1;
    C2 = params.C2;
    ngrid = params.ngrid;
    maxvel = params.maxvel;
    u_mut = params.u_mut;
    fun = MultiObj.fun;
    nVar = MultiObj.nVar;
    var_min = MultiObj.var_min';
    var_max = MultiObj.var_max';
    
    % Initialize population
    POS = repmat((var_max - var_min), Np, 1) .* rand(Np, nVar) + repmat(var_min, Np, 1);
    VEL = zeros(Np, nVar);
    POS_fit = fun(POS);
    PBEST = POS;
    PBEST_fit = POS_fit;
    DOMINATED = checkDomination(POS_fit);
    REP.pos = POS(~DOMINATED, :);
    REP.pos_fit = POS_fit(~DOMINATED, :);
    REP = updateGrid(REP, ngrid);
    maxvel = (var_max - var_min) .* maxvel ./ 100;
    gen = 1;

    % Main MOPSO loop
    while gen <= maxgen
        % Select leader
        h = selectLeader(REP);
        
        % Update speeds and positions
        VEL = W .* VEL + C1 * rand(Np, nVar) .* (PBEST - POS) ...
              + C2 * rand(Np, nVar) .* (repmat(REP.pos(h,:), Np, 1) - POS);
        POS = POS + VEL;
        
        % Perform mutation
        POS = mutation(POS, gen, maxgen, Np, var_max, var_min, nVar, u_mut);
        
        % Check boundaries
        [POS, VEL] = checkBoundaries(POS, VEL, maxvel, var_max, var_min);       
        
        % Evaluate the population
        POS_fit = fun(POS);
        
        % Update the repository
        REP = updateRepository(REP, POS, POS_fit, ngrid);
        if size(REP.pos, 1) > Nr
            REP = deleteFromRepository(REP, size(REP.pos, 1) - Nr, ngrid);
        end
        
        % Update the best positions found so far for each particle
        pos_best = dominate(POS_fit, PBEST_fit);
        best_pos = ~dominate(PBEST_fit, POS_fit);
        best_pos(rand(Np, 1) >= 0.5) = 0;
        PBEST_fit(pos_best, :) = POS_fit(pos_best, :);
        PBEST(pos_best, :) = POS(pos_best, :);
        
        % Update generation
        gen = gen + 1;
    end
end

% Helper functions for MOPSO
function REP = updateRepository(REP, POS, POS_fit, ngrid)
    DOMINATED = checkDomination(POS_fit);
    REP.pos = [REP.pos; POS(~DOMINATED, :)];
    REP.pos_fit = [REP.pos_fit; POS_fit(~DOMINATED, :)];
    DOMINATED = checkDomination(REP.pos_fit);
    REP.pos_fit = REP.pos_fit(~DOMINATED, :);
    REP.pos = REP.pos(~DOMINATED, :);
    REP = updateGrid(REP, ngrid);
end

function [POS, VEL] = checkBoundaries(POS, VEL, maxvel, var_max, var_min)
    Np = size(POS, 1);
    MAXLIM = repmat(var_max, Np, 1);
    MINLIM = repmat(var_min, Np, 1);
    MAXVEL = repmat(maxvel, Np, 1);
    MINVEL = repmat(-maxvel, Np, 1);

    VEL(VEL > MAXVEL) = MAXVEL(VEL > MAXVEL);
    VEL(VEL < MINVEL) = MINVEL(VEL < MINVEL);
    VEL(POS > MAXLIM) = -VEL(POS > MAXLIM);
    POS(POS > MAXLIM) = MAXLIM(POS > MAXLIM);
    VEL(POS < MINLIM) = -VEL(POS < MINLIM);
    POS(POS < MINLIM) = MINLIM(POS < MINLIM);
end

function dom_vector = checkDomination(fitness)
    Np = size(fitness, 1);
    dom_vector = zeros(Np, 1);
    all_perm = nchoosek(1:Np, 2);
    all_perm = [all_perm; [all_perm(:, 2), all_perm(:, 1)]];

    d = dominate(fitness(all_perm(:, 1), :), fitness(all_perm(:, 2), :));
    dominated_particles = unique(all_perm(d == 1, 2));
    dom_vector(dominated_particles) = 1;
end

function d = dominate(x, y)
    d = all(x <= y, 2) & any(x < y, 2);
end

function REP = updateGrid(REP, ngrid)
    ndim = size(REP.pos_fit, 2);
    REP.hypercube_limits = zeros(ngrid + 1, ndim);
    for dim = 1:ndim
        REP.hypercube_limits(:, dim) = linspace(min(REP.pos_fit(:, dim)), max(REP.pos_fit(:, dim)), ngrid + 1)';
    end
    
    npar = size(REP.pos_fit, 1);
    REP.grid_idx = zeros(npar, 1);
    REP.grid_subidx = zeros(npar, ndim);
    for n = 1:npar
        idnames = [];
        for d = 1:ndim
            REP.grid_subidx(n, d) = find(REP.pos_fit(n, d) <= REP.hypercube_limits(:, d)', 1, 'first') - 1;
            if REP.grid_subidx(n, d) == 0, REP.grid_subidx(n, d) = 1; end
            idnames = [idnames ',' num2str(REP.grid_subidx(n, d))];
        end
        REP.grid_idx(n) = eval(['sub2ind(ngrid.*ones(1,ndim)' idnames ');']);
    end
    
    REP.quality = zeros(ngrid, 2);
    ids = unique(REP.grid_idx);
    for i = 1:length(ids)
        REP.quality(i, 1) = ids(i);
        REP.quality(i, 2) = 10 / sum(REP.grid_idx == ids(i));
    end
end

function selected = selectLeader(REP)
    prob = cumsum(REP.quality(:, 2));
    sel_hyp = REP.quality(find(rand(1, 1) * max(prob) <= prob, 1, 'first'), 1);
    
    idx = 1:length(REP.grid_idx);
    selected = idx(REP.grid_idx == sel_hyp);
    selected = selected(randi(length(selected)));
end

function REP = deleteFromRepository(REP, n_extra, ngrid)
    crowding = zeros(size(REP.pos, 1), 1);
    for m = 1:size(REP.pos_fit, 2)
        [m_fit, idx] = sort(REP.pos_fit(:, m), 'ascend');
        m_up = [m_fit(2:end); Inf];
        m_down = [Inf; m_fit(1:end-1)];
        distance = (m_up - m_down) ./ (max(m_fit) - min(m_fit));
        [~, idx] = sort(idx, 'ascend');
        crowding = crowding + distance(idx);
    end
    crowding(isnan(crowding)) = Inf;
    
    [~, del_idx] = sort(crowding, 'ascend');
    del_idx = del_idx(1:n_extra);
    REP.pos(del_idx, :) = [];
    REP.pos_fit(del_idx, :) = [];
    REP = updateGrid(REP, ngrid); 
end

function POS = mutation(POS, gen, maxgen, Np, var_max, var_min, nVar, u_mut)
    fract = Np / 3 - floor(Np / 3);
    if fract < 0.5
        sub_sizes = [ceil(Np/3), round(Np/3), round(Np/3)];
    else
        sub_sizes = [round(Np/3), round(Np/3), floor(Np/3)];
    end
    cum_sizes = cumsum(sub_sizes);
    
    nmut = round(u_mut * sub_sizes(2));
    if nmut > 0
        idx = cum_sizes(1) + randperm(sub_sizes(2), nmut);
        POS(idx, :) = repmat((var_max - var_min), nmut, 1) .* rand(nmut, nVar) + repmat(var_min, nmut, 1);
    end
    
    per_mut = (1 - gen / maxgen)^(5 * nVar);
    nmut = round(per_mut * sub_sizes(3));
    if nmut > 0
        idx = cum_sizes(2) + randperm(sub_sizes(3), nmut);
        POS(idx, :) = repmat((var_max - var_min), nmut, 1) .* rand(nmut, nVar) + repmat(var_min, nmut, 1);
    end
end

% ZDT1 Objective Function
function y = zdt1(x)
    n = length(x);
    f1 = x(1);
    g = 1 + (9/(n - 1)) * sum(x(2:end));
    f2 = g * (1 - sqrt(f1 / g));
    y = [f1, f2];
end

% ZDT Function for MOEA/D
function z = ZDT(x)
    n = numel(x);
    f1 = x(1);
    g = 1 + 9 / (n-1) * sum(x(2:end));
    h = 1 - sqrt(f1 / g);
    f2 = g * h;
    z = [f1; f2];
end

% Additional helper functions for MOEMA and MOEA/D (if needed)
% ... (Include any additional helper functions here)

% Performance Metrics
function [GD, IGD, Spread, HV] = calculateMetrics(paretoFront, truePF)
    % Generation Distance (GD)
    GD = mean(min(pdist2(paretoFront, truePF), [], 2));
    
    % Inverted Generation Distance (IGD)
    IGD = mean(min(pdist2(truePF, paretoFront), [], 2));
    
    % Spread
    d = min(pdist2(paretoFront, paretoFront), [], 2);
    d_mean = mean(d);
    d_extreme = max(min(pdist2(paretoFront, [0, 1; 1, 0]), [], 2));
    Spread = (d_extreme + sum(abs(d - d_mean))) / (d_extreme + length(d) * d_mean);
    
    % Hypervolume (HV)
    ref_point = max(truePF, [], 1) * 1.1;
    HV = hypervolume(paretoFront, ref_point);
end

% Hypervolume calculation
function hv = hypervolume(points, ref_point)
    sorted_points = sortrows(points);
    hv = 0;
    for i = 1:size(sorted_points, 1)
        if i == 1
            area = (ref_point(1) - sorted_points(i, 1)) * (ref_point(2) - sorted_points(i, 2));
        else
            area = (sorted_points(i, 1) - sorted_points(i-1, 1)) * (ref_point(2) - sorted_points(i, 2));
        end
        hv = hv + area;
    end
end