can you change the zdt1 and zdt to the Multi-Objective Cost Function that i provided:

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