can you add  the comparison Plot  Before and After optimization for temp, iaq and lux. please display the maximum, minimum and averange. please show full 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
];

% 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);
for k = 1:48
    idx = 3*(k-1) + 1;
    lb(idx) = 20;       % Temp: lower bound
    ub(idx) = 25.56;    % Temp: 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 MOEMA
[paretoFrontMOEMA, ~, ~, ~, ~] = MOEMA(NPop, T, lb, ub, dim, @(X) multi_objective_cost_function(X, dataSH));

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

% Plot Comparisons
figure;
hold on;
plot(paretoFrontMOEMA(:,1), paretoFrontMOEMA(:,2), 'bs', 'DisplayName', 'MOEMA');
plot(paretoFrontNSGA2(:,1), paretoFrontNSGA2(:,2), 'ro', 'DisplayName', 'NSGA-II');
xlabel('-Average Comfort Index');
ylabel('Average Energy Consumption');
title('Comparison of Pareto Fronts: MOEMA vs NSGA-II');
legend('show');
grid on;


% Calculate performance metrics
referenceSet = [paretoFrontMOEMA; paretoFrontNSGA2];
[~, uniqueIdx] = unique(referenceSet, 'rows');
referenceSet = referenceSet(uniqueIdx, :);

% Hypervolume
hvMOEMA = hypervolume(paretoFrontMOEMA, max(referenceSet, [], 1) * 1.1);
hvNSGA2 = hypervolume(paretoFrontNSGA2, max(referenceSet, [], 1) * 1.1);

% Generational Distance
gdMOEMA = generational_distance(paretoFrontMOEMA, referenceSet);
gdNSGA2 = generational_distance(paretoFrontNSGA2, referenceSet);

% Inverted Generational Distance
igdMOEMA = inverted_generational_distance(paretoFrontMOEMA, referenceSet);
igdNSGA2 = inverted_generational_distance(paretoFrontNSGA2, referenceSet);

% Spread (Diversity)
spreadMOEMA = spread_metric(paretoFrontMOEMA);
spreadNSGA2 = spread_metric(paretoFrontNSGA2);

% Display results
disp('Performance Metrics:');
disp(['MOEMA Hypervolume: ', num2str(hvMOEMA)]);
disp(['NSGA-II Hypervolume: ', num2str(hvNSGA2)]);
disp(['MOEMA Generational Distance: ', num2str(gdMOEMA)]);
disp(['NSGA-II Generational Distance: ', num2str(gdNSGA2)]);
disp(['MOEMA Inverted Generational Distance: ', num2str(igdMOEMA)]);
disp(['NSGA-II Inverted Generational Distance: ', num2str(igdNSGA2)]);
disp(['MOEMA Spread: ', num2str(spreadMOEMA)]);
disp(['NSGA-II Spread: ', num2str(spreadNSGA2)]);

% Plot Comparisons
figure;
hold on;
scatter(paretoFrontMOEMA(:,1), paretoFrontMOEMA(:,2), 50, 'bs', 'filled', 'DisplayName', 'MOEMA');
scatter(paretoFrontNSGA2(:,1), paretoFrontNSGA2(:,2), 50, 'ro', 'filled', 'DisplayName', 'NSGA-II');
xlabel('-Average Comfort Index');
ylabel('Average Energy Consumption');
title('Comparison of Pareto Fronts: MOEMA vs NSGA-II');
legend('show');
grid on;



% Functions

% ---------- MOEMA Optimization Function ----------
function [paretoFront, paretoSet, Trajectories, fitness_history, position_history] = MOEMA(NPop, T, lb, ub, dim, FITNESSFCN)
    fitness_history = zeros(NPop, T, 2);
    position_history = zeros(NPop, T, dim);
    Trajectories = zeros(dim, T);
    
    pop = initialization(NPop, dim, ub, lb);
    pop_dad = pop(1:floor(NPop/2), :);
    pop_mum = pop(floor(NPop/2)+1:end, :);
    
    obj_dad = evaluate_population(pop_dad, FITNESSFCN);
    obj_mum = evaluate_population(pop_mum, FITNESSFCN);
    obj = [obj_dad; obj_mum];
    
    archive = update_archive([], [], pop, obj);
    
    for gen = 1:T
        k1 = randperm(size(pop_dad,1));
        k2 = randperm(size(pop_mum,1));
        k1x = pop_dad(k1, :);
        k2x = pop_mum(k2, :);
        
        for k = 1:size(k1x,1)
            imates = randn(1, dim);
            p = rand(1, dim);
            r = rand();
            
            if mean(imates) >= 0
                child = p .* k1x(k, :) + (1 - p) .* k2x(k, :);
            else
                child = p .* k2x(k, :) + (1 - p) .* k1x(k, :);
            end
            
            if ~isempty(archive.pop)
                best_pop = archive.pop(randi(size(archive.pop,1)), :);
                K = rand(1, dim) < 0.8;
                child = best_pop .* (1 - K) + child .* K;
            end
            
            child = max(min(child, ub), lb);
            obj_child = FITNESSFCN(child);
            
            [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);
            
            if r < 0.2 && ~isempty(archive.pop)
                temp1 = archive.pop(randi(size(archive.pop,1)), :);
                temp_child = rand(1, dim) .* temp1;
                temp_child = max(min(temp_child, ub), lb);
                obj_temp_child = FITNESSFCN(temp_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
        
        pop = [pop_dad; pop_mum];
        obj = [obj_dad; obj_mum];
        archive = update_archive(archive.pop, archive.obj, pop, obj);
        
        fitness_history(:, gen, :) = obj;
        position_history(:, gen, :) = pop;
        
        if ~isempty(archive.pop)
            Trajectories(:, gen) = archive.pop(randi(size(archive.pop,1)), :)';
        end
        
        disp(['MOEMA Iteration #', num2str(gen), ' completed. Archive size: ', num2str(size(archive.pop,1))]);
    end
    
    paretoFront = archive.obj;
    paretoSet = archive.pop;
end

% ---------- NSGA-II Optimization Function ----------
function [paretoFront, paretoSet] = NSGA2(NPop, T, lb, ub, dim, FITNESSFCN)
    pop = initialization(NPop, dim, ub, lb);
    fitness = evaluate_population(pop, FITNESSFCN);

    for gen = 1:T
        [fronts, ranks] = non_dominated_sort(fitness);
        mating_pool = tournament_selection(pop, ranks, NPop);
        offspring = crossover(mating_pool, 0.9, dim);
        offspring = mutation(offspring, 1/dim, dim);
        offspring_fitness = evaluate_population(offspring, FITNESSFCN);
        
        comb_pop = [pop; offspring];
        comb_fitness = [fitness; offspring_fitness];
        
        [fronts, ~] = non_dominated_sort(comb_fitness);
        [pop, fitness] = select_new_population(comb_pop, comb_fitness, fronts, NPop);
        
        disp(['NSGA-II Generation #', num2str(gen), ' completed.']);
    end

    paretoFront = fitness;
    paretoSet = pop;
end

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

% ---------- Evaluate Population Function ----------
function obj = evaluate_population(pop, FITNESSFCN)
    num_individuals = size(pop,1);
    obj = zeros(num_individuals, 2); 
    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)
    if dominates(child_obj, obj(index, :))
        pop(index, :) = child;
        obj(index, :) = child_obj;
    elseif ~dominates(obj(index, :), child_obj)
        if rand() < 0.5
            pop(index, :) = child;
            obj(index, :) = child_obj;
        end
    end
end

% ---------- Update Archive Function ----------
function archive = update_archive(archive_pop, archive_obj, pop, obj)
    combined_pop = [archive_pop; pop];
    combined_obj = [archive_obj; obj];
    
    [is_dominated, ~] = pareto_dominance(combined_obj);
    non_dominated = ~is_dominated;
    
    archive.pop = combined_pop(non_dominated, :);
    archive.obj = combined_obj(non_dominated, :);
end

% ---------- Pareto Dominance Function ----------
function [is_dominated, domination_count] = pareto_dominance(objectives)
    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)
    result = all(obj1 <= obj2) && any(obj1 < obj2);
end

% ---------- Non-Dominated Sort Function ----------
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

% ---------- Tournament Selection Function ----------
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

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

% ---------- Mutation Function ----------
function mutated_pop = mutation(pop, mutation_prob, num_vars)
    eta_m = 20;
    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

% ---------- Select New Population Function ----------
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

% ---------- Crowding Distance Function ----------
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

% ---------- Multi-Objective Cost Function ----------
function objectives = multi_objective_cost_function(X, dataSH)
    if length(X) ~= 144
        error('Decision vector X must have 144 elements (3 variables per 48 periods).');
    end
    
    p1 = 1/3; 
    p2 = 1/3; 
    p3 = 1/3;
    
    Td = 22.78;
    Ld = 800;
    Id = 800;
    
    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);
        
        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;
        
        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

% Add performance metric functions
function hv = hypervolume(paretoFront, referencePoint)
    hv = 0;
    sortedFront = sortrows(paretoFront);
    for i = 1:size(sortedFront, 1)
        if i == size(sortedFront, 1)
            hv = hv + prod(referencePoint - sortedFront(i, :));
        else
            hv = hv + prod([sortedFront(i+1, 1), referencePoint(2)] - sortedFront(i, :));
        end
    end
end

function gd = generational_distance(approximationSet, referenceSet)
    distances = pdist2(approximationSet, referenceSet);
    gd = sqrt(sum(min(distances, [], 2).^2)) / size(approximationSet, 1);
end

function igd = inverted_generational_distance(approximationSet, referenceSet)
    distances = pdist2(referenceSet, approximationSet);
    igd = sum(min(distances, [], 2)) / size(referenceSet, 1);
end

function s = spread_metric(paretoFront)
    sortedFront = sortrows(paretoFront);
    distances = sqrt(sum(diff(sortedFront).^2, 2));
    d = mean(distances);
    df = sqrt(sum((sortedFront(1,:) - sortedFront(end,:)).^2));
    s = (df + sum(abs(distances - d))) / (df + size(paretoFront, 1) * d);
end