can you add one more multi optimization method to be compared to:
% compare_MOEMA_NSGAII_SmartBuilding.m
% Comprehensive script to compare MOEMA and NSGA-II optimizations for
% Smart Building Comfort (Temperature, Lux, IAQ) and Energy Consumption

%% ------------------ Initialization ------------------

clc;
clear;
close all;

%% ------------------ Data Initialization ------------------

% Data from UTHM 
% Columns: Time_instance, Temp_(F), Temp_(C), lux, IAQ
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
];

%% ------------------ Define Fitness Function ------------------

% Define the fitness function handle with dataSH
fitness_fcn = @(X) multi_objective_cost_function(X, dataSH);

%% ------------------ Optimization Parameters Setup ------------------

% Number of decision variables: 3 per period for 48 periods
dim = 144;

% Define lower and upper bounds for each variable
% Each period has T, L, IAQ variables with specific bounds
lb = zeros(1, dim);
ub = zeros(1, dim);

% Define bounds per variable type
% T: Temperature [20, 25.56]
% L: Light [720, 880]
% IAQ: IAQ [700, 880]

% Assign bounds for each variable in the decision vector
for k = 1:48
    idx = 3*(k-1) + 1;
    lb(idx) = 20;       % Temperature lower bound
    ub(idx) = 25.56;    % Temperature upper bound
    
    lb(idx+1) = 720;    % Light lower bound
    ub(idx+1) = 880;    % Light upper bound
    
    lb(idx+2) = 700;    % IAQ lower bound
    ub(idx+2) = 880;    % IAQ upper bound
end

%% ------------------ MOEMA Optimization ------------------

% Parameters for MOEMA
NPop_MOEMA = 500;  % Population size
T_MOEMA = 500;     % Number of iterations

% Run MOEMA
disp('Starting MOEMA Optimization...');
tic; % Start timer
[paretoFront_MOEMA, paretoSet_MOEMA, ~, ~, ~] = MOEMA(NPop_MOEMA, T_MOEMA, lb, ub, dim, fitness_fcn);
MOEMA_Time = toc; % Record execution time
disp('MOEMA Optimization Completed.');

%% ------------------ NSGA-II Optimization Using gamultiobj ------------------

% NSGA-II Parameters
population_size = 500;      % Number of individuals in each generation
max_generations = 500;      % Number of generations
crossover_fraction = 0.8;   % Fraction of offspring from crossover

% Define options for gamultiobj
options = optimoptions('gamultiobj', ...
    'PopulationSize', population_size, ...
    'MaxGenerations', max_generations, ...
    'CrossoverFraction', crossover_fraction, ...
    'ParetoFraction', 0.35, ...
    'Display', 'iter', ...
    'PlotFcn', {@gaplotpareto});

% Run NSGA-II
disp('Starting NSGA-II Optimization...');
tic; % Start timer
[paretoSet_NSGAII, paretoFront_NSGAII] = gamultiobj(fitness_fcn, dim, [], [], [], [], lb, ub, options);
NSGAII_Time = toc; % Record execution time
disp('NSGA-II Optimization Completed.');

%% ------------------ Compare Pareto Fronts ------------------

% Plotting both Pareto fronts together for comparison
figure;
hold on;

% Plot MOEMA Pareto Front
plot(paretoFront_MOEMA(:,1), paretoFront_MOEMA(:,2), 'ro', 'MarkerSize',6, 'DisplayName', 'MOEMA');

% Plot NSGA-II Pareto Front
plot(paretoFront_NSGAII(:,1), paretoFront_NSGAII(:,2), 'b+', 'MarkerSize',6, 'DisplayName', 'NSGA-II');

xlabel('-Average Comfort Index');
ylabel('Average Energy Consumption');
title('Pareto Front Comparison: MOEMA vs. NSGA-II');
legend('Location', 'best');
grid on;
set(gca, 'FontSize', 12);
hold off;

%% ------------------ Compare Temp, Lux, and IAQ ------------------

% Select Pareto solutions with minimum energy consumption
[~, minEC_idx_MOEMA] = min(paretoFront_MOEMA(:,2));
optimized_solution_MOEMA = paretoSet_MOEMA(minEC_idx_MOEMA, :);  % 1x144 vector

[~, minEC_idx_NSGAII] = min(paretoFront_NSGAII(:,2));
optimized_solution_NSGAII = paretoSet_NSGAII(minEC_idx_NSGAII, :);  % 1x144 vector

% Extract optimized settings for each period
optimized_settings_MOEMA = reshape(optimized_solution_MOEMA, 3, 48)';  % 48x3 matrix: [T, L, IAQ]
optimized_settings_NSGAII = reshape(optimized_solution_NSGAII, 3, 48)';  % 48x3 matrix: [T, L, IAQ]

% Original settings from dataSH
original_temp = dataSH(:,3);   % Temperature in C
original_lux = dataSH(:,4);    % lux
original_IAQ = dataSH(:,5);    % IAQ

% Optimized settings
optimized_temp_MOEMA = optimized_settings_MOEMA(:,1);
optimized_lux_MOEMA = optimized_settings_MOEMA(:,2);
optimized_IAQ_MOEMA = optimized_settings_MOEMA(:,3);

optimized_temp_NSGAII = optimized_settings_NSGAII(:,1);
optimized_lux_NSGAII = optimized_settings_NSGAII(:,2);
optimized_IAQ_NSGAII = optimized_settings_NSGAII(:,3);

% Desired values
desired_temp = 22.78;  % As per cost function
desired_lux = 800;     % As per cost function
desired_IAQ = 800;     % As per cost function

% Create a time vector for 48 periods
time = dataSH(:,1);  % Assuming time instance is from 1 to 48

% ------------- Temperature Comparison -------------

figure;
plot(time, original_temp, 'k-', 'LineWidth',1.5); hold on;
plot(time, optimized_temp_MOEMA, 'r--', 'LineWidth',1.5);
plot(time, optimized_temp_NSGAII, 'b-.', 'LineWidth',1.5);
yline(desired_temp, 'g:', 'LineWidth',1.5);
xlabel('Time Instance');
ylabel('Temperature (°C)');
title('Temperature Before and After Optimization');
legend('Original', 'MOEMA Optimized', 'NSGA-II Optimized', 'Desired', 'Location', 'best');
grid on;
set(gcf, 'Position', [100, 100, 800, 400]);  % Set window size
hold off;

% ------------- Light Level Comparison -------------

figure;
plot(time, original_lux, 'k-', 'LineWidth',1.5); hold on;
plot(time, optimized_lux_MOEMA, 'r--', 'LineWidth',1.5);
plot(time, optimized_lux_NSGAII, 'b-.', 'LineWidth',1.5);
yline(desired_lux, 'g:', 'LineWidth',1.5);
xlabel('Time Instance');
ylabel('Light Level (lux)');
title('Light Level Before and After Optimization');
legend('Original', 'MOEMA Optimized', 'NSGA-II Optimized', 'Desired', 'Location', 'best');
grid on;
set(gcf, 'Position', [200, 200, 800, 400]);  % Set window size
hold off;

% ------------- IAQ Comparison -------------

figure;
plot(time, original_IAQ, 'k-', 'LineWidth',1.5); hold on;
plot(time, optimized_IAQ_MOEMA, 'r--', 'LineWidth',1.5);
plot(time, optimized_IAQ_NSGAII, 'b-.', 'LineWidth',1.5);
yline(desired_IAQ, 'g:', 'LineWidth',1.5);
xlabel('Time Instance');
ylabel('IAQ');
title('Indoor Air Quality (IAQ) Before and After Optimization');
legend('Original', 'MOEMA Optimized', 'NSGA-II Optimized', 'Desired', 'Location', 'best');
grid on;
set(gcf, 'Position', [300, 300, 800, 400]);  % Set window size
hold off;

%% ------------------ Compare Metrics ------------------

% Metrics for MOEMA
num_pareto_MOEMA = size(paretoFront_MOEMA, 1);
min_CI_MOEMA = min(paretoFront_MOEMA(:,1));   % Since it's -CI_avg
min_EC_MOEMA = min(paretoFront_MOEMA(:,2));

% Metrics for NSGA-II
num_pareto_NSGAII = size(paretoFront_NSGAII, 1);
min_CI_NSGAII = min(paretoFront_NSGAII(:,1)); % Since it's -CI_avg
min_EC_NSGAII = min(paretoFront_NSGAII(:,2));

% Display Metrics
fprintf('\n--- Optimization Metrics Comparison ---\n');
fprintf('MOEMA:\n');
fprintf('Number of Pareto Front Solutions: %d\n', num_pareto_MOEMA);
fprintf('Maximum Average Comfort Index: %.4f\n', -min_CI_MOEMA);
fprintf('Minimum Average Energy Consumption: %.4f\n', min_EC_MOEMA);
fprintf('Execution Time: %.2f seconds\n', MOEMA_Time);

fprintf('\nNSGA-II:\n');
fprintf('Number of Pareto Front Solutions: %d\n', num_pareto_NSGAII);
fprintf('Maximum Average Comfort Index: %.4f\n', -min_CI_NSGAII);
fprintf('Minimum Average Energy Consumption: %.4f\n', min_EC_NSGAII);
fprintf('Execution Time: %.2f seconds\n', NSGAII_Time);

%% ------------------ Calculate Original Metrics ------------------

% Construct the original decision vector from dataSH
original_decision = zeros(1, dim);
for k = 1:48
    original_decision(3*(k-1)+1) = dataSH(k, 3); % Temperature (°C)
    original_decision(3*(k-1)+2) = dataSH(k, 4); % Light (lux)
    original_decision(3*(k-1)+3) = dataSH(k, 5); % IAQ
end

% Evaluate the original metrics using the multi-objective cost function
original_objectives = multi_objective_cost_function(original_decision, dataSH);
original_CI_avg = -original_objectives(1); % Convert back to positive Comfort Index
original_EC_avg = original_objectives(2);  % Electricity Consumption

%% ------------------ Extract Optimized Metrics ------------------

% Metrics for MOEMA
max_CI_MOEMA = -min(paretoFront_MOEMA(:,1));    % Maximum Average Comfort Index
min_EC_MOEMA = min(paretoFront_MOEMA(:,2));     % Minimum Average Energy Consumption

% Metrics for NSGA-II
max_CI_NSGAII = -min(paretoFront_NSGAII(:,1));  % Maximum Average Comfort Index
min_EC_NSGAII = min(paretoFront_NSGAII(:,2));   % Minimum Average Energy Consumption

%% ------------------ Create Comparison Table ------------------

% Define the row names as a column vector
Rows = {'Original'; 'MOEMA Optimized'; 'NSGA-II Optimized'};

% Define the metrics
Avg_Comfort_Index = [original_CI_avg; max_CI_MOEMA; max_CI_NSGAII];
Avg_Electricity_Consumption = [original_EC_avg; min_EC_MOEMA; min_EC_NSGAII];
Optimization_Time = [NaN; MOEMA_Time; NSGAII_Time]; % No transpose

% Calculate Improvements
Improvement_CI_MOEMA = ((max_CI_MOEMA - original_CI_avg) / original_CI_avg) * 100;
Improvement_CI_NSGAII = ((max_CI_NSGAII - original_CI_avg) / original_CI_avg) * 100;

Improvement_EC_MOEMA = ((original_EC_avg - min_EC_MOEMA) / original_EC_avg) * 100;
Improvement_EC_NSGAII = ((original_EC_avg - min_EC_NSGAII) / original_EC_avg) * 100;

% Handle potential infinite or NaN improvements
Improvement_CI_MOEMA = finiteness_handling(Improvement_CI_MOEMA);
Improvement_CI_NSGAII = finiteness_handling(Improvement_CI_NSGAII);
Improvement_EC_MOEMA = finiteness_handling(Improvement_EC_MOEMA);
Improvement_EC_NSGAII = finiteness_handling(Improvement_EC_NSGAII);

% Define Improvement vectors as column vectors
Improvement_Comfort_Index = [NaN; Improvement_CI_MOEMA; Improvement_CI_NSGAII];
Improvement_Electricity_Consumption = [NaN; Improvement_EC_MOEMA; Improvement_EC_NSGAII];

% Create the table ensuring all variables have the same number of rows
Comparison_Table = table(Rows, Avg_Comfort_Index, Improvement_Comfort_Index, ...
    Avg_Electricity_Consumption, Improvement_Electricity_Consumption, Optimization_Time, ...
    'VariableNames', {'Scenario', 'Average_Comfort_Index', 'CI_Improvement_Percentage', ...
    'Average_Electricity_Consumption', 'EC_Improvement_Percentage', 'Optimization_Time_sec'});

% Display the table
disp('--- Comfort Index and Electricity Consumption Comparison ---');
disp(Comparison_Table);

%% ------------------ Additional Comparisons ------------------

% ------------- Bar Chart Comparison -------------

figure('Position', [400, 400, 1000, 500]);

% Subplot for Average Comfort Index
subplot(1,2,1);
bar(Avg_Comfort_Index, 'FaceColor', [0.2 0.6 0.5]);
set(gca, 'XTickLabel', Rows, 'FontSize', 12);
ylabel('Average Comfort Index');
title('Comparison of Average Comfort Index');
grid on;

% Subplot for Average Electricity Consumption
subplot(1,2,2);
bar(Avg_Electricity_Consumption, 'FaceColor', [0.8 0.3 0.3]);
set(gca, 'XTickLabel', Rows, 'FontSize', 12);
ylabel('Average Electricity Consumption');
title('Comparison of Average Electricity Consumption');
grid on;

sgtitle('Optimization Results Comparison');

% ------------- Pie Chart for Improvement Percentage -------------

figure('Position', [600, 100, 800, 600]);

% Comfort Index Improvement
subplot(2,2,1);
if all(isfinite([Improvement_CI_MOEMA, Improvement_CI_NSGAII]))
    pie([Improvement_CI_MOEMA, Improvement_CI_NSGAII], {'MOEMA', 'NSGA-II'});
    title('Average Comfort Index Improvement (%)');
else
    disp('Non-finite values detected in Comfort Index improvements. Pie chart not generated for this metric.');
end

% Electricity Consumption Improvement
subplot(2,2,2);
if all(isfinite([Improvement_EC_MOEMA, Improvement_EC_NSGAII]))
    pie([Improvement_EC_MOEMA, Improvement_EC_NSGAII], {'MOEMA', 'NSGA-II'});
    title('Average Electricity Consumption Improvement (%)');
else
    disp('Non-finite values detected in Electricity Consumption improvements. Pie chart not generated for this metric.');
end

% Optimization Time Comparison
subplot(2,2,3);
bar([MOEMA_Time, NSGAII_Time], 'FaceColor', [0.4 0.4 0.7]);
set(gca, 'XTickLabel', {'MOEMA', 'NSGA-II'}, 'FontSize', 12);
ylabel('Execution Time (sec)');
title('Optimization Time Comparison');
grid on;

% Handle potential overlapped titles
sgtitle('Additional Optimization Comparisons');

% ------------- Pareto Fronts Scatter Plot -------------

figure;
hold on;

% Scatter plot with EC vs CI for MOEMA
scatter(paretoFront_MOEMA(:,2), -paretoFront_MOEMA(:,1), 50, 'ro', 'filled', 'DisplayName', 'MOEMA');

% Scatter plot with EC vs CI for NSGA-II
scatter(paretoFront_NSGAII(:,2), -paretoFront_NSGAII(:,1), 50, 'b+', 'DisplayName', 'NSGA-II');

xlabel('Average Energy Consumption');
ylabel('Average Comfort Index');
title('Pareto Fronts: MOEMA vs. NSGA-II');
legend('Location', 'best');
grid on;
set(gca, 'FontSize', 12);
hold off;

%% ------------------ Determine Best Method ------------------

% Criteria for Best Method:
% - Higher Comfort Index Improvement
% - Greater Electricity Consumption Reduction
% - Shorter Optimization Time

% Initialize scores
MOEMA_Score = 0;
NSGAII_Score = 0;

% Assign scores based on improvements

% For Comfort Index
if Improvement_CI_MOEMA > Improvement_CI_NSGAII
    MOEMA_Score = MOEMA_Score + 1;
elseif Improvement_CI_MOEMA < Improvement_CI_NSGAII
    NSGAII_Score = NSGAII_Score + 1;
end

% For Electricity Consumption
if Improvement_EC_MOEMA > Improvement_EC_NSGAII
    MOEMA_Score = MOEMA_Score + 1;
elseif Improvement_EC_MOEMA < Improvement_EC_NSGAII
    NSGAII_Score = NSGAII_Score + 1;
end

% For Optimization Time (prefer shorter time)
if MOEMA_Time < NSGAII_Time
    MOEMA_Score = MOEMA_Score + 1;
elseif MOEMA_Time > NSGAII_Time
    NSGAII_Score = NSGAII_Score + 1;
end

% Display the best method based on score
if MOEMA_Score > NSGAII_Score
    best_method = 'MOEMA';
elseif NSGAII_Score > MOEMA_Score
    best_method = 'NSGA-II';
else
    best_method = 'Both methods perform equally well';
end

fprintf('\n--- Best Optimization Method ---\n');
fprintf('MOEMA Score: %d\n', MOEMA_Score);
fprintf('NSGA-II Score: %d\n', NSGAII_Score);
fprintf('Best Method: %s\n', best_method);

%% ------------------ Helper Function ------------------

% ---------- Finite Handling Function ----------
function val = finiteness_handling(val)
    % FINITENESS_HANDLING Replaces non-finite values with 0
    if ~isfinite(val)
        val = 0;
    end
end

%% ------------------ Functions ------------------

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

% ---------- MOEMA Optimization Function ----------
function [paretoFront, paretoSet, Trajectories, fitness_history, position_history] = MOEMA(NPop, T, lb, ub, dim, FITNESSFCN)
    % MOEMA: Multi-Objective Evolutionary Mating Algorithm
    % 
    % Solves multi-objective optimization problems of the form:
    %   min F(X) = [f1(X), f2(X)]  subject to: lb <= X <= ub
    % 
    % Inputs:
    %   NPop      - Population size
    %   T         - Number of iterations
    %   lb        - Lower bounds (1 x dim)
    %   ub        - Upper bounds (1 x dim)
    %   dim       - Number of decision variables
    %   FITNESSFCN - Function handle for the multi-objective fitness function
    % 
    % Outputs:
    %   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:floor(NPop/2), :);      % Parent subset 1
    pop_mum = pop(floor(NPop/2)+1:end, :);  % Parent subset 2
    
    % Evaluation of objective function for initial population
    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);
    
    % Evolutionary Loop
    for gen = 1:T
        % Shuffle parents
        k1 = randperm(size(pop_dad,1));
        k2 = randperm(size(pop_mum,1));
        k1x = pop_dad(k1, :);
        k2x = pop_mum(k2, :);
        
        % Generate children
        for k = 1:size(k1x,1)
            % Mating process
            imates = randn(1, dim);  % Mating factor
            p = rand(1, dim);        % Crossover probability
            r = rand();
            
            if mean(imates) >= 0  % Decide mating based on imates
                child = p .* k1x(k, :) + (1 - p) .* k2x(k, :);
            else
                child = p .* k2x(k, :) + (1 - p) .* k1x(k, :);
            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;  % Modification probability
                child = best_pop .* (1 - K) + child .* K;
            end
            
            % Ensure child is within bounds
            child = max(min(child, ub), lb);
            
            % Evaluate child
            obj_child = FITNESSFCN(child);
            
            % Update population using non-dominated sorting
            [pop_dad, obj_dad] = update_population(pop_dad, obj_dad, child, obj_child, k);
            [pop_mum, obj_mum] = update_population(pop_mum, obj_mum, child, obj_child, k);
        end
        
        % Exploration process: Introduce diversity
        if rand() < 0.2
            if ~isempty(archive.pop)
                temp1 = archive.pop(randi(size(archive.pop,1)), :);
                temp_child = rand(1, dim) .* temp1;  % Exploration child
                temp_child = max(min(temp_child, ub), lb);
                obj_temp_child = FITNESSFCN(temp_child);
                
                % Randomly replace individuals in both parents
                replace_idx_dad = randi(size(pop_dad,1));
                replace_idx_mum = randi(size(pop_mum,1));
                
                pop_dad(replace_idx_dad, :) = temp_child;
                obj_dad(replace_idx_dad, :) = obj_temp_child;
                
                pop_mum(replace_idx_mum, :) = temp_child;
                obj_mum(replace_idx_mum, :) = obj_temp_child;
            end
        end
        
        % Combine updated populations
        pop = [pop_dad; pop_mum];
        obj = [obj_dad; obj_mum];
        
        % Update archive with non-dominated solutions
        archive = update_archive(archive.pop, archive.obj, pop, obj);
        
        % Record fitness and positions
        fitness_history(:, gen, :) = obj;
        position_history(:, gen, :) = pop;
        
        % Update Trajectories with a random non-dominated solution
        if ~isempty(archive.pop)
            Trajectories(:, gen) = archive.pop(randi(size(archive.pop,1)), :)';
        end
        
        % Display iteration status every 50 generations
        if mod(gen, 50) == 0 || gen == 1 || gen == T
            disp(['MOEMA Iteration #', num2str(gen), ' completed. Archive size: ', num2str(size(archive.pop,1))]);
        end
    end
    
    % Return Pareto front and Pareto set
    paretoFront = archive.obj;
    paretoSet = archive.pop;
end

% ---------- Initialization Function ----------
function pop = initialization(NPop, dim, ub, lb)
    % INITIALIZATION Generates an initial population within bounds
    pop = rand(NPop, dim) .* (ub - lb) + lb;
end

% ---------- Evaluate Population Function ----------
function obj = evaluate_population(pop, FITNESSFCN)
    % EVALUATE_POPULATION Evaluates the objective function for each individual
    num_individuals = size(pop,1);
    obj = zeros(num_individuals, 2);  % Assuming 2 objectives
    parfor i = 1:num_individuals
        obj(i, :) = FITNESSFCN(pop(i, :));
    end
end

% ---------- Update Population Function ----------
function [pop, obj] = update_population(pop, obj, child, child_obj, index)
    % UPDATE_POPULATION Updates the population based on dominance
    if dominates(child_obj, obj(index, :))
        pop(index, :) = child;
        obj(index, :) = child_obj;
    elseif dominates(obj(index, :), child_obj)
        % Do nothing, parent dominates child
    else
        % Neither dominates; randomly decide to replace
        if rand() < 0.5
            pop(index, :) = child;
            obj(index, :) = child_obj;
        end
    end
end

% ---------- Update Archive Function ----------
function archive = update_archive(archive_pop, archive_obj, pop, obj)
    % UPDATE_ARCHIVE Updates the archive with non-dominated solutions
    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

% ---------- Pareto Dominance Function ----------
function [is_dominated, domination_count] = pareto_dominance(objectives)
    % PARETO_DOMINANCE Determines domination relationships among solutions
    % 
    % Outputs:
    %   is_dominated - Logical array indicating if a solution is dominated
    %   domination_count - Number of times a solution is dominated
    
    n_solutions = size(objectives,1);
    is_dominated = false(n_solutions,1);
    domination_count = zeros(n_solutions, 1);
    
    for i = 1:n_solutions
        for j = 1:n_solutions
            if i ~= j
                if dominates(objectives(j,:), objectives(i,:))
                    is_dominated(i) = true;
                    domination_count(j) = domination_count(j) + 1;
                end
            end
        end
    end
end

% ---------- Dominates Function ----------
function result = dominates(obj1, obj2)
    % DOMINATES Checks if obj1 dominates obj2
    result = all(obj1 <= obj2) && any(obj1 < obj2);
end