can you combine and compare all these methods in one code. i want you to create a table to compare all the method. and please display the original electricity comsumption and the optimized electricity comsumption and original comfort and optimized comfort. here is all the method i want you to combine and compare. show me full code

Moema:
% run_MOEMA_MultiObjective.m
% Comprehensive script to run MOEMA optimization for Smart Building Comfort 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
];

% ------------------ MOEMA Parameters ------------------

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

% Population size and number of iterations
NPop = 300;  % Population size
T = 300;     % Number of iterations

% 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

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

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

% ------------------ Run MOEMA ------------------

[paretoFront, paretoSet, Trajectories, fitness_history, position_history] = MOEMA(NPop, T, lb, ub, dim, fitness_fcn);

% ------------------ Plot the Pareto Front ------------------

figure;
plot(paretoFront(:,1), paretoFront(:,2), 'ro', 'MarkerSize',8, 'LineWidth',1.5);
xlabel('-Average Comfort Index');
ylabel('Average Energy Consumption');
title('Pareto Front for Smart Building Optimization using MOEMA');
grid on;
set(gca, 'FontSize', 12);

% ------------------ Plot Settings Before and After Optimization ------------------

% Select a Pareto solution to compare (e.g., with minimum energy consumption)
[~, minEC_idx] = min(paretoFront(:,2));
optimized_solution = paretoSet(minEC_idx, :);  % 1x144 vector

% Extract optimized settings for each period
optimized_settings = reshape(optimized_solution, 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 = optimized_settings(:,1);
optimized_lux = optimized_settings(:,2);
optimized_IAQ = optimized_settings(:,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 Plot
figure;
plot(time, original_temp, 'b-', 'LineWidth',1.5); hold on;
plot(time, optimized_temp, 'r--', 'LineWidth',1.5);
yline(desired_temp, 'g:', 'LineWidth',1.5);
xlabel('Time Instance');
ylabel('Temperature (°C)');
title('Temperature Before and After Optimization');
legend('Original', 'Optimized', 'Desired', 'Location', 'best');
grid on;
set(gcf, 'Position', [100, 100, 800, 400]);  % Set window size

% Light Level Plot
figure;
plot(time, original_lux, 'b-', 'LineWidth',1.5); hold on;
plot(time, optimized_lux, 'r--', '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', 'Optimized', 'Desired', 'Location', 'best');
grid on;
set(gcf, 'Position', [200, 200, 800, 400]);  % Set window size

% IAQ Plot
figure;
plot(time, original_IAQ, 'b-', 'LineWidth',1.5); hold on;
plot(time, optimized_IAQ, 'r--', 'LineWidth',1.5);
yline(desired_IAQ, 'g:', 'LineWidth',1.5);
xlabel('Time Instance');
ylabel('IAQ');
title('IAQ Before and After Optimization');
legend('Original', 'Optimized', 'Desired', 'Location', 'best');
grid on;
set(gcf, 'Position', [300, 300, 800, 400]);  % Set window size

% ------------------ Display Metrics ------------------

disp(['Number of solutions in Pareto Front: ', num2str(size(paretoFront, 1))]);
disp(['Minimum -Average Comfort Index: ', num2str(min(paretoFront(:,1)),4)]);
disp(['Minimum Average Energy Consumption: ', num2str(min(paretoFront(:,2)),4)]);

%% ------------------ 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);
            
            % Exploration process
            if r < 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 decide where to place the exploration child
                    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
        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
        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

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

Mopso:
% run_MOPSO_MultiObjective.m
% Comprehensive script to run MOPSO optimization for Smart Building Comfort and Energy Consumption

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

clc;
clear;
close all;

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

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

% ------------------ MOPSO Parameters ------------------

params.Np = 300;        % Population size
params.Nr = 150;        % Repository size
params.maxgen = 300;    % Maximum number of generations
params.W = 0.5;         % Inertia weight
params.C1 = 1.5;        % Personal acceleration coefficient
params.C2 = 1.5;        % Social acceleration coefficient
params.ngrid = 20;      % Number of grids in each dimension
params.maxvel = 5;      % Maximum velocity
params.u_mut = 0.1;     % Mutation rate

% Dimensionality and bounds
dim = 144;
lb = zeros(1, dim);
ub = zeros(1, dim);

% Define variable-specific bounds
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

% Custom fitness function including shared data
MultiObj.fun = @(pos) cell2mat(arrayfun(@(i) multi_objective_cost_function(pos(i,:), dataSH), 1:size(pos,1), 'UniformOutput', false)');
MultiObj.nVar = dim;
MultiObj.var_min = lb(:)';
MultiObj.var_max = ub(:)';

% ------------------ Run MOPSO ------------------

REP = MOPSO(params, MultiObj);

% ------------------ Plot the Pareto Front ------------------

figure;
plot(REP.pos_fit(:,1), REP.pos_fit(:,2), 'ro', 'MarkerSize',8, 'LineWidth',1.5);
xlabel('-Average Comfort Index');
ylabel('Average Energy Consumption');
title('Pareto Front for Smart Building Optimization using MOPSO');
grid on;
set(gca, 'FontSize', 12);
% Add to the script after the MOPSO function has finished running:

% ------------------ Plot Settings Before and After Optimization ------------------

% Select a Pareto solution to compare (e.g., with minimum energy consumption)
[~, minEC_idx] = min(REP.pos_fit(:,2));
optimized_solution = REP.pos(minEC_idx, :);  % 1x144 vector

% Extract optimized settings for each period
optimized_settings = reshape(optimized_solution, 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 = optimized_settings(:,1);
optimized_lux = optimized_settings(:,2);
optimized_IAQ = optimized_settings(:,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 Plot
figure;
plot(time, original_temp, 'b-', 'LineWidth', 1.5); hold on;
plot(time, optimized_temp, 'r--', 'LineWidth', 1.5);
yline(desired_temp, 'g:', 'LineWidth', 1.5);
xlabel('Time Instance');
ylabel('Temperature (°C)');
title('Temperature Before and After Optimization');
legend('Original', 'Optimized', 'Desired', 'Location', 'best');
grid on;
set(gcf, 'Position', [100, 100, 800, 400]);  % Set window size

% Light Level Plot
figure;
plot(time, original_lux, 'b-', 'LineWidth', 1.5); hold on;
plot(time, optimized_lux, 'r--', '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', 'Optimized', 'Desired', 'Location', 'best');
grid on;
set(gcf, 'Position', [200, 200, 800, 400]);  % Set window size

% IAQ Plot
figure;
plot(time, original_IAQ, 'b-', 'LineWidth', 1.5); hold on;
plot(time, optimized_IAQ, 'r--', 'LineWidth', 1.5);
yline(desired_IAQ, 'g:', 'LineWidth', 1.5);
xlabel('Time Instance');
ylabel('IAQ');
title('IAQ Before and After Optimization');
legend('Original', 'Optimized', 'Desired', 'Location', 'best');
grid on;
set(gcf, 'Position', [300, 300, 800, 400]);  % Set window size
%% ------------------ Functions ------------------

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

function REP = MOPSO(params, MultiObj)
    % Extract parameters
    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;

    % Initialize repository
    DOMINATED = checkDomination(POS_fit);
    REP.pos = POS(~DOMINATED, :);
    REP.pos_fit = POS_fit(~DOMINATED, :);
    REP = updateGrid(REP, ngrid);

    % Set velocity bounds
    maxvel = (var_max - var_min) .* maxvel ./ 100;

    % Main optimization loop
    for gen = 1:maxgen
        % Select leader
        h = selectLeader(REP);
        
        % Update velocities 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 positions
        POS_fit = fun(POS);
        
        % Update 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 personal bests
        is_better = dominates(POS_fit, PBEST_fit);
        PBEST_fit(is_better, :) = POS_fit(is_better, :);
        PBEST(is_better, :) = POS(is_better, :);
    end
end

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 = dominates(fitness(all_perm(:, 1), :), fitness(all_perm(:, 2), :));
    dom_vector(unique(all_perm(d == 1, 2))) = 1;
end

function d = dominates(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)
        [~, idx] = sort(REP.pos_fit(:, m), 'ascend');
        m_up = [REP.pos_fit(idx(2:end), m); Inf];
        m_down = [Inf; REP.pos_fit(idx(1:end-1), m)];
        distance = m_up - m_down;
        distance = distance / (max(REP.pos_fit(:, m)) - min(REP.pos_fit(:, m)));
        [~, orig_idx] = sort(idx, 'ascend');
        crowding = crowding + distance(orig_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), floor(Np/3)];
    else
        sub_sizes = [floor(Np/3), round(Np/3), round(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

MOEAD:
% run_MOEAD_MultiObjective.m
% Comprehensive script to run MOEA/D optimization for Smart Building Comfort 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
];

%% ------------------ MOEA/D Parameters ------------------

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

% Number of objectives
nObj = 2;

% Population size and number of iterations
nPop = 300;          % Population size (number of weight vectors / subproblems)
MaxIt = 300;         % Maximum number of iterations
T = 20;              % Neighborhood size

% Crossover and Mutation Parameters
pCrossover = 1;      % Probability of crossover
pMutation = 1/dim;   % Probability of mutation
etaC = 20;           % Crossover distribution index
etaM = 20;           % Mutation distribution index

% Define lower and upper bounds for each variable
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]
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

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

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

%% ------------------ MOEA/D Execution ------------------

% Execute MOEA/D
[paretoFront, paretoSet, fitness_history, position_history] = MOEAD(dim, nObj, nPop, MaxIt, T, lb, ub, fitness_fcn, ...
    pCrossover, pMutation, etaC, etaM);

%% ------------------ Plot the Pareto Front ------------------

figure;
plot(paretoFront(:,1), paretoFront(:,2), 'ro', 'MarkerSize',8, 'LineWidth',1.5);
xlabel('-Average Comfort Index');
ylabel('Average Energy Consumption');
title('Pareto Front for Smart Building Optimization using MOEA/D');
grid on;
set(gca, 'FontSize', 12);

%% ------------------ Plot Settings Before and After Optimization ------------------

% Select a Pareto solution to compare (e.g., with minimum energy consumption)
[~, minEC_idx] = min(paretoFront(:,2));
optimized_solution = paretoSet(minEC_idx, :); % 1x144 vector

% Extract optimized settings for each period
optimized_settings = reshape(optimized_solution, 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 = optimized_settings(:,1);
optimized_lux = optimized_settings(:,2);
optimized_IAQ = optimized_settings(:,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 Plot
figure;
plot(time, original_temp, 'b-', 'LineWidth',1.5); hold on;
plot(time, optimized_temp, 'r--', 'LineWidth',1.5);
yline(desired_temp, 'g:', 'LineWidth',1.5);
xlabel('Time Instance');
ylabel('Temperature (°C)');
title('Temperature Before and After Optimization');
legend('Original', 'Optimized', 'Desired', 'Location', 'best');
grid on;
set(gcf, 'Position', [100, 100, 800, 400]); % Set window size

% Light Level Plot
figure;
plot(time, original_lux, 'b-', 'LineWidth',1.5); hold on;
plot(time, optimized_lux, 'r--', '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', 'Optimized', 'Desired', 'Location', 'best');
grid on;
set(gcf, 'Position', [200, 200, 800, 400]); % Set window size

% IAQ Plot
figure;
plot(time, original_IAQ, 'b-', 'LineWidth',1.5); hold on;
plot(time, optimized_IAQ, 'r--', 'LineWidth',1.5);
yline(desired_IAQ, 'g:', 'LineWidth',1.5);
xlabel('Time Instance');
ylabel('IAQ');
title('IAQ Before and After Optimization');
legend('Original', 'Optimized', 'Desired', 'Location', 'best');
grid on;
set(gcf, 'Position', [300, 300, 800, 400]); % Set window size

%% ------------------ Display Metrics ------------------

disp(['Number of solutions in Pareto Front: ', num2str(size(paretoFront, 1))]);
disp(['Minimum -Average Comfort Index: ', num2str(min(paretoFront(:,1)),4)]);
disp(['Minimum Average Energy Consumption: ', num2str(min(paretoFront(:,2)),4)]);

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

    % Reshape X into 48x3 matrix
    settings = reshape(X, 3, 48)';
    Tk = settings(:,1); % Temperature settings (48x1)
    Lk = settings(:,2); % Light settings (48x1)
    Ik = settings(:,3); % IAQ settings (48x1)

    % --- Calculate Comfort Index Components ---
    Ti = 1 - ((Tk - Td) ./ Tc').^2;
    Li = 1 - ((Lk - Ld) ./ Lc').^2;
    Ii = 1 - ((Ik - Id) ./ Ic').^2;
    CI = p1 * Ti + p2 * Li + p3 * Ii;
    CI_avg = mean(CI);

    % --- Calculate Energy Consumption ---
    ET = pT * abs(Tk - Tc');
    EL = pL * abs(Lk - Lc');
    EI = pI * abs(Ik - Ic');
    EC = ET + EL + EI;
    EC_avg = mean(EC);

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

% ---------- MOEAD Optimization Function ----------
function [paretoFront, paretoSet, fitness_history, position_history] = MOEAD(dim, nObj, nPop, MaxIt, T, lb, ub, FITNESSFCN, ...
    pCrossover, pMutation, etaC, etaM)
% MOEAD: Multi-Objective Evolutionary Algorithm based on Decomposition
%
% Inputs:
%   dim         - Number of decision variables
%   nObj        - Number of objectives
%   nPop        - Population size (number of weight vectors / subproblems)
%   MaxIt       - Maximum number of iterations
%   T           - Neighborhood size
%   lb          - Lower bounds (1 x dim)
%   ub          - Upper bounds (1 x dim)
%   FITNESSFCN  - Fitness function handle
%   pCrossover  - Crossover probability
%   pMutation   - Mutation probability
%   etaC        - Crossover distribution index
%   etaM        - Mutation distribution index
%
% Outputs:
%   paretoFront     - Non-dominated solutions (objectives)
%   paretoSet       - Non-dominated solutions (decision variables)
%   fitness_history - History of fitness values (ideal point z over iterations)
%   position_history- History of population positions

    %% Step 1: Generate Weight Vectors
    [W, nPop] = GenerateUniformWeightVectors(nPop, nObj);

    %% Step 2: Initialize Neighborhood
    Dist = squareform(pdist(W));
    [~, sortedIndices] = sort(Dist, 2, 'ascend');
    Neighborhood = sortedIndices(:, 2:T+1); % Exclude self (first)

    %% Step 3: Initialize Population
    Pop = struct('Position', [], 'Cost', [], 'G', []);
    for i = 1:nPop
        Pop(i).Position = unifrnd(lb, ub);
        Pop(i).Cost = FITNESSFCN(Pop(i).Position);
    end

    % Initialize Ideal Point as the min over all solutions
    allCosts = arrayfun(@(p) p.Cost, Pop, 'UniformOutput', false);
    allCosts = cell2mat(allCosts');
    z = min(allCosts, [], 1); % 1 x nObj

    % Calculate initial decomposed cost
    for i = 1:nPop
        Pop(i).G = DecomposedCost(Pop(i).Cost, W(i,:), z);
    end

    %% Step 4: Initialize Archive
    Archive = UpdateParetoArchive(Pop, dim, nObj, lb, ub, FITNESSFCN, nPop);

    %% Initialize History Variables
    fitness_history = zeros(nObj, MaxIt);
    position_history = cell(nPop, MaxIt);

    %% Step 5: Main Loop
    for it = 1:MaxIt
        for i = 1:nPop
            % Mating selection based on neighborhood
            neighbors = Neighborhood(i, :);
            selector = randperm(length(neighbors),2);
            P1 = Pop(neighbors(selector(1))).Position;
            P2 = Pop(neighbors(selector(2))).Position;

            % Crossover
            C = SimulatedBinaryCrossover(P1, P2, pCrossover, etaC, lb, ub);

            % Mutation
            C = PolynomialMutation(C, pMutation, etaM, lb, ub);

            % Evaluation
            ChildCost = FITNESSFCN(C);

            % Update Ideal Point
            z = min(z, ChildCost);

            % Decomposition
            g_child = DecomposedCost(ChildCost, W(i,:), z);

            % Update the population
            for j = 1:length(Neighborhood(i,:))
                nj = Neighborhood(i,j);
                if g_child < Pop(nj).G
                    Pop(nj).Position = C;
                    Pop(nj).Cost = ChildCost;
                    Pop(nj).G = DecomposedCost(ChildCost, W(nj,:), z);
                end
            end
        end

        % Update Archive with current population
        Archive = UpdateParetoArchive(Pop, dim, nObj, lb, ub, FITNESSFCN, nPop);

        % Record Fitness and Positions
        for p = 1:nPop
            position_history{p, it} = Pop(p).Position;
        end
        % Record ideal point
        fitness_history(:, it) = z;

        % Display iteration information
        disp(['Iteration ' num2str(it) ': Archive size = ' num2str(size(Archive,1))]);
    end

    %% Step 6: Extract Pareto Front and Set
    ParetoCost = cell2mat({Archive.Cost}');
    ParetoSet = cell2mat({Archive.Position}');

    paretoFront = ParetoCost;
    paretoSet = ParetoSet;
end

% ---------- Generate Uniform Weight Vectors ----------
function [W, nPop_out] = GenerateUniformWeightVectors(nPop, nObj)
% GenerateUniformWeightVectors Generates uniformly distributed weight vectors
% for bi-objective optimization using a triangular lattice.

    if nObj ~= 2
        error('Only 2-objective problems are supported in this implementation.');
    end

    % Calculate the number of divisions (H)
    H = 1;
    while (H*(H+1)/2 < nPop)
        H = H + 1;
    end

    % Generate weight vectors
    W = [];
    for i = 0:H
        for j = 0:H-i
            w1 = i / H;
            w2 = j / H;
            if w1 + w2 <=1
                W = [W; w1 w2];
            end
        end
    end

    % If more weight vectors than needed, randomly select
    if size(W,1) > nPop
        idx = randperm(size(W,1), nPop);
        W = W(idx, :);
    end

    nPop_out = size(W,1);
end

% ---------- Decomposed Cost Calculation using Tchebycheff Approach ----------
function g = DecomposedCost(fx, lambda, z)
% DecomposedCost Calculates the Tchebycheff decomposition for scalarizing multiple objectives.
%
% Inputs:
%   fx     - Objective vector of the solution (1 x nObj)
%   lambda - Weight vector (1 x nObj)
%   z      - Ideal point (1 x nObj)
%
% Output:
%   g      - Scalarized fitness value

    % Ensure lambda and fx are row vectors
    lambda = lambda(:)';
    fx = fx(:)';
    z = z(:)';

    g = max(lambda .* abs((fx - z)));
end

% ---------- Simulated Binary Crossover (SBX) ----------
function child = SimulatedBinaryCrossover(parent1, parent2, pCrossover, etaC, lb, ub)
% SimulatedBinaryCrossover Performs SBX crossover on two parents.
%
% Inputs:
%   parent1    - First parent solution (1 x dim)
%   parent2    - Second parent solution (1 x dim)
%   pCrossover - Crossover probability
%   etaC       - Crossover distribution index
%   lb         - Lower bounds (1 x dim)
%   ub         - Upper bounds (1 x dim)
%
% Output:
%   child      - Offspring solution (1 x dim)

    child = parent1;
    if rand <= pCrossover
        for i = 1:length(parent1)
            if rand <= 0.5
                if abs(parent1(i) - parent2(i)) > 1e-14
                    x1 = min(parent1(i), parent2(i));
                    x2 = max(parent1(i), parent2(i));
                    rand_val = rand;
                    beta_q = 1.0 + (2.0*(x1 - lb(i))/(x2 - x1));
                    alpha = 2.0 - beta_q^-(etaC+1);
                    if rand_val <= 1.0/alpha
                        betaq = (rand_val * alpha)^(1.0/(etaC+1));
                    else
                        betaq = (1.0/(2.0 - rand_val*alpha))^(1.0/(etaC+1));
                    end
                    c1 = 0.5*((x1 + x2) - betaq*(x2 - x1));

                    beta_q = 1.0 + (2.0*(ub(i) - x2)/(x2 - x1));
                    alpha = 2.0 - beta_q^-(etaC+1);
                    if rand_val <= 1.0/alpha
                        betaq = (rand_val * alpha)^(1.0/(etaC+1));
                    else
                        betaq = (1.0/(2.0 - rand_val*alpha))^(1.0/(etaC+1));
                    end
                    c2 = 0.5*((x1 + x2) + betaq*(x2 - x1));

                    if rand <= 0.5
                        child(i) = c2;
                    else
                        child(i) = c1;
                    end
                else
                    child(i) = parent1(i);
                end
            else
                child(i) = parent1(i);
            end
        end
    end
end

% ---------- Polynomial Mutation ----------
function mutant = PolynomialMutation(child, pMutation, etaM, lb, ub)
% PolynomialMutation Applies polynomial mutation to a solution.
%
% Inputs:
%   child      - Solution to mutate (1 x dim)
%   pMutation  - Mutation probability per gene
%   etaM       - Mutation distribution index
%   lb         - Lower bounds (1 x dim)
%   ub         - Upper bounds (1 x dim)
%
% Output:
%   mutant     - Mutated solution (1 x dim)

    mutant = child;
    for i = 1:length(child)
        if rand <= pMutation
            delta1 = (child(i) - lb(i)) / (ub(i) - lb(i));
            delta2 = (ub(i) - child(i)) / (ub(i) - lb(i));
            rand_val = rand;
            mut_pow = 1.0/(etaM+1.0);
            if rand_val < 0.5
                xy = 1.0 - delta1;
                val = 2.0*rand_val + (1.0 - 2.0*rand_val)*(xy)^(etaM+1);
                delta_q = val^(mut_pow) - 1.0;
            else
                xy = 1.0 - delta2;
                val = 2.0*(1.0 - rand_val) + 2.0*(rand_val - 0.5)*(xy)^(etaM+1);
                delta_q = 1.0 - val^(mut_pow);
            end
            mutant(i) = child(i) + delta_q*(ub(i) - lb(i));
            mutant(i) = min(max(mutant(i), lb(i)), ub(i));
        end
    end
end

% ---------- Update Pareto Archive ----------
function Archive = UpdateParetoArchive(Pop, dim, nObj, lb, ub, FITNESSFCN, nPop)
% UpdateParetoArchive Maintains an archive of non-dominated solutions.
%
% Inputs:
%   Pop        - Current population structure
%   dim        - Number of decision variables
%   nObj       - Number of objectives
%   lb         - Lower bounds
%   ub         - Upper bounds
%   FITNESSFCN - Fitness function handle
%   nPop       - Population size
%
% Output:
%   Archive    - Updated archive of non-dominated solutions (struct array)

    % Extract costs and positions
    allCosts = arrayfun(@(p) p.Cost, Pop, 'UniformOutput', false);
    allCosts = cell2mat(allCosts');
    allPositions = arrayfun(@(p) p.Position, Pop, 'UniformOutput', false);
    allPositions = cell2mat(allPositions')';

    % Combine costs and positions
    combined = [allCosts, allPositions'];

    % Perform non-dominated sorting
    [isDominated, ~] = pareto_dominance(combined(:,1:nObj));

    % Keep non-dominated solutions
    nonDominated = combined(~isDominated, :);

    % If archive size exceeds, randomly select
    if size(nonDominated,1) > nPop
        indices = randperm(size(nonDominated,1), nPop);
        nonDominated = nonDominated(indices, :);
    end

    % Create Archive as struct array
    Archive = struct('Cost', num2cell(nonDominated(:,1:nObj),2), ...
                    'Position', num2cell(nonDominated(:,nObj+1:end),2));
end

% ---------- Pareto Dominance Function ----------
function [is_dominated, dominated_count] = pareto_dominance(objectives)
% pareto_dominance Determines domination relationships among solutions.
%
% Inputs:
%   objectives - Matrix of objective values (nSolutions x nObjectives)
%
% Outputs:
%   is_dominated     - Logical array indicating if a solution is dominated
%   dominated_count  - Number of times a solution is dominated

    nSolutions = size(objectives,1);
    is_dominated = false(nSolutions,1);
    dominated_count = zeros(nSolutions,1);

    for i = 1:nSolutions
        for j = 1:nSolutions
            if i ~= j
                if dominates(objectives(j,:), objectives(i,:))
                    is_dominated(i) = true;
                    dominated_count(j) = dominated_count(j) +1;
                end
            end
        end
    end
end

% ---------- Dominates Function ----------
function flag = dominates(obj1, obj2)
% dominates Returns true if obj1 dominates obj2.

    flag = all(obj1 <= obj2) && any(obj1 < obj2);
end