can you compare these 2 method into 1 script code. please also compare the pareto code, compare the Plot Before and After Optimization for tempreture, lux and iaq. please give me full  code. 
moema code:
% 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

ngsa ii code:
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
dim = 144; % (3 variables per period for 48 periods)
NPop = 300; % Population size
T = 300; % Number of iterations

% Define bounds
lb = zeros(1, dim);
ub = zeros(1, dim);

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

% Run NSGA-II
[paretoFront, paretoSet] = NSGA2(NPop, T, lb, ub, dim, @(X) multi_objective_cost_function(X, dataSH));

% ------------------ Plot the Pareto Front ------------------
% Plot the results
figure;
plot(paretoFront(:,1), paretoFront(:,2), 'ro', 'DisplayName', 'NSGA-II');
xlabel('-Average Comfort Index');
ylabel('Average Energy Consumption');
title('NSGA-II Pareto Front for Smart Building Optimization');
legend('show');
grid on;

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

function [paretoFront, paretoSet] = NSGA2(NPop, T, lb, ub, dim, FITNESSFCN)
% Adapted NSGA-II implementation for specific function


% Initialize population
pop = initialization(NPop, dim, ub, lb);
fitness = evaluate_population(pop, FITNESSFCN);

for gen = 1:T
    % Non-dominated sorting
    [fronts, ranks] = non_dominated_sort(fitness);
    
    % Create mating pool using tournament selection
    mating_pool = tournament_selection(pop, ranks, NPop);
    
    % Generate offspring using crossover and mutation
    offspring = crossover(mating_pool, 0.9, dim);
    offspring = mutation(offspring, 1/dim, dim);
    
    % Evaluate offspring
    offspring_fitness = evaluate_population(offspring, FITNESSFCN);
    
    % Combine population and offspring
    comb_pop = [pop; offspring];
    comb_fitness = [fitness; offspring_fitness];
    
    % Non-dominated sorting on combined population
    [fronts, ~] = non_dominated_sort(comb_fitness);
    
    % Create new population
    [pop, fitness] = select_new_population(comb_pop, comb_fitness, fronts, NPop);
    
    % Display generation progress
    disp(['Generation #', num2str(gen), ' completed.']);
end

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

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

function fitness = evaluate_population(pop, FITNESSFCN)
fitness = zeros(size(pop, 1), 2); % Assuming 2 objectives
for i = 1:size(pop, 1)
fitness(i, :) = FITNESSFCN(pop(i, :));
end
end

function [fronts, rank] = non_dominated_sort(fitness)
num_individuals = size(fitness, 1);
domination_count = zeros(num_individuals, 1);
dominated_set = cell(num_individuals, 1);
rank = zeros(num_individuals, 1);


fronts = cell(1);
for i = 1:num_individuals
    for j = 1:num_individuals
        if i ~= j
            if dominates(fitness(i, :), fitness(j, :))
                dominated_set{i} = [dominated_set{i}, j];
            elseif dominates(fitness(j, :), fitness(i, :))
                domination_count(i) = domination_count(i) + 1;
            end
        end
    end
    
    if domination_count(i) == 0
        rank(i) = 1;
        fronts{1} = [fronts{1}, i];
    end
end

current_front = 1;
while ~isempty(fronts{current_front})
    next_front = [];
    for i = fronts{current_front}
        for j = dominated_set{i}
            domination_count(j) = domination_count(j) - 1;
            if domination_count(j) == 0
                rank(j) = current_front + 1;
                next_front = [next_front, j];
            end
        end
    end
    current_front = current_front + 1;
    fronts{current_front} = next_front;
end
fronts = fronts(1:current_front-1);
end

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

function selected_pop = tournament_selection(pop, ranks, num_selected)
selected_pop = zeros(num_selected, size(pop, 2));
for i = 1:num_selected
candidates = randperm(size(pop, 1), 2);
if ranks(candidates(1)) < ranks(candidates(2))
selected_pop(i, :) = pop(candidates(1), :);
else
selected_pop(i, :) = pop(candidates(2), :);
end
end
end

function offspring = crossover(pop, crossover_prob, num_vars)
offspring = zeros(size(pop));
for i = 1:2:size(pop, 1)
if rand < crossover_prob
parent1 = pop(i, :);
parent2 = pop(i+1, :);
cross_point = randi([1, num_vars-1], 1, 1);
offspring(i, :) = [parent1(1:cross_point), parent2(cross_point+1:end)];
offspring(i+1, :) = [parent2(1:cross_point), parent1(cross_point+1:end)];
else
offspring(i, :) = pop(i, :);
offspring(i+1, :) = pop(i+1, :);
end
end
end

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

function [pop, fitness] = select_new_population(comb_pop, comb_fitness, fronts, NPop)
pop = [];
fitness = [];
for i = 1:length(fronts)
if size(pop, 1) + length(fronts{i}) > NPop
break;
end
pop = [pop; comb_pop(fronts{i}, :)];
fitness = [fitness; comb_fitness(fronts{i}, :)];
end


remaining = NPop - size(pop, 1);
if remaining > 0
    selected_front = fronts{i};
    distances = crowding_distance(comb_fitness(selected_front, :));
    [~, sorted_indices] = sort(distances, 'descend');
    selected = selected_front(sorted_indices(1:remaining));
    pop = [pop; comb_pop(selected, :)];
    fitness = [fitness; comb_fitness(selected, :)];
end
end

function distances = crowding_distance(fitness)
num_individuals = size(fitness, 1);
num_objectives = size(fitness, 2);
distances = zeros(num_individuals, 1);


for m = 1:num_objectives
    [sorted_values, sorted_idx] = sort(fitness(:, m));
    distances(sorted_idx(1)) = inf;
    distances(sorted_idx(end)) = inf;
    for i = 2:num_individuals-1
        if sorted_values(end) ~= sorted_values(1)
            distances(sorted_idx(i)) = distances(sorted_idx(i)) + ...
                (sorted_values(i+1) - sorted_values(i-1)) / (sorted_values(end) - sorted_values(1));
        end
    end
end
end

% ---------- MOC Function ----------
function objectives = multi_objective_cost_function(X, dataSH)
% 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
p1 = 1/3; % Temperature Weight
p2 = 1/3; % Light Weight
p3 = 1/3; % IAQ Weight

% Desired setpoints
Td = 22.78;
Ld = 800;
Id = 800;

% Energy consumption weights
pT = 9;
pL = 1;
pI = 1;

Tc = dataSH(:,3)';
Lc = dataSH(:,4)';
Ic = dataSH(:,5)';

CI = zeros(48,1);
EC = zeros(48,1);

for k = 1:48
    Tk = X(3*k-2);
    Lk = X(3*k-1);
    Ik = X(3*k);
    
    % Comfort Index
    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;
    
    % 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

CI_avg = mean(CI);
EC_avg = mean(EC);

objectives = [-CI_avg, EC_avg];
end