from this data, can you create a function to optimize the user comfort and energy consumption in smart building. I am trying to use this code to do multi objective optimization. can you try use my moema code. give me full code. and can you give me a plot for Temperature,lux,and IAQ before the optimization and after the optimization. 

here is the data:
% Data from UTHM
% 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 code:
% Main script to run MOEMA optimization
clc;
clear;
close all;
% Set problem parameters
NPop = 100; % Population size
T = 200; % Number of iterations
dim = 30; % Number of decision variables
lb = zeros(1, dim); % Lower bound
ub = ones(1, dim); % Upper bound

% Run MOEMA
[paretoFront, paretoSet, Trajectories, fitness_history, position_history] = MOEMA(NPop, T, lb, ub, dim, @zdt1);

% Plot the Pareto front
figure;
plot(paretoFront(:,1), paretoFront(:,2), 'ro');
xlabel('f1');
ylabel('f2');
title('Pareto Front for ZDT1 using MOEMA');

% Calculate and display metrics
disp(['Number of solutions in Pareto Front: ', num2str(size(paretoFront, 1))]);
disp(['Minimum f1: ', num2str(min(paretoFront(:,1)))]);
disp(['Minimum f2: ', num2str(min(paretoFront(:,2)))]);

function [paretoFront, paretoSet, Trajectories, fitness_history, position_history] = MOEMA(NPop, T, lb, ub, dim, FITNESSFCN)
% Multi-Objective Evolutionary Mating Algorithm (MOEMA)
% MOEMA attempts to solve multi-objective problems of the following form:
% min F(X) = [f1(X), f2(X), ..., fm(X)] subject to: lb <= X <= ub
%
% FITNESSFCN - function handle of the multi-objective fitness function
% lb - lower bounds on X
% ub - upper bounds on X
% T - number of iterations
% NPop - size of the population (class size)
% paretoFront - set of non-dominated solutions in objective space
% paretoSet - set of non-dominated solutions in decision space
% Trajectories - best solution trajectory over iterations
% fitness_history - history of fitness values
% position_history - history of population positions

% Preallocation
fitness_history = zeros(NPop, T, 2); % Assuming 2 objectives
position_history = zeros(NPop, T, dim);
Trajectories = zeros(dim, T);

% Generation of initial population
pop = initialization(NPop, dim, ub, lb);
pop_dad = pop(1:NPop/2, :); % Eqn (1) in the paper
pop_mum = pop(NPop/2+1:end, :); % Eqn (2) in the paper

% Evaluation of objective function for initialization
obj_dad = evaluate_population(pop_dad, FITNESSFCN);
obj_mum = evaluate_population(pop_mum, FITNESSFCN);
obj = [obj_dad; obj_mum];

% Initialize archive with non-dominated solutions
archive = update_archive([], [], pop, obj);

for gen = 1:T % Start iteration
k1 = randperm(NPop/2);
k2 = randperm(NPop/2);
k1x = pop_dad(k1, :);
k2x = pop_mum(k2, :);


for k = 1:NPop/2
    % Mating process
    Imates = 1 + (var(k1x(k, :)) - var(k2x(k, :))); % Eqn (4) in the paper
    p = randn(1, dim);  % Eqn (6) in the paper
    r = rand();
    
    if Imates >= 0  % Upper part of Eqn (5) in the paper
        temp1 = p .* k1x(k, :);
        temp2 = (1-p) .* k2x(k, :);
        pop_child = temp1 + temp2;
    else  % Lower part of Eqn (5) in the paper
        temp1 = p .* k2x(k, :);
        temp2 = (1-p) .* k1x(k, :);
        pop_child = temp1 + temp2;
    end
    
    % Impact of best environment (using random non-dominated solution)
    if ~isempty(archive.pop)
        best_pop = archive.pop(randi(size(archive.pop, 1)), :);
        K = rand(1, dim) < 0.8;  % Eqn (8) in the paper
        pop_child = best_pop .* (1-K) + pop_child .* K;  % Eqn (7) in the paper
    end
    
    % Bounding of the solution
    pop_child = max(min(ub, pop_child), lb);
    
    % Evaluate child
    obj_child = FITNESSFCN(pop_child);
    
    % Update population using non-dominated sorting
    [pop_dad, obj_dad] = update_population(pop_dad, obj_dad, pop_child, obj_child, k1(k));
    [pop_mum, obj_mum] = update_population(pop_mum, obj_mum, pop_child, obj_child, k2(k));
    
    % Exploration process
    if r < 0.2
        if ~isempty(archive.pop)
            temp1 = archive.pop(randi(size(archive.pop, 1)), :);
            pop_child = rand(1, dim) .* temp1;  % Eqn (9) in the paper
            pop_child = max(min(ub, pop_child), lb);
            obj_child = FITNESSFCN(pop_child);
            
            % Update population using non-dominated sorting
            [pop_dad, obj_dad] = update_population(pop_dad, obj_dad, pop_child, obj_child, randi(NPop/2));
            [pop_mum, obj_mum] = update_population(pop_mum, obj_mum, pop_child, obj_child, randi(NPop/2));
        end
    end
end

% Combine updated pool dad and mum
pop = [pop_dad; pop_mum];
obj = [obj_dad; obj_mum];

% Update archive
archive = update_archive(archive.pop, archive.obj, pop, obj);

% Record the history at each iteration
fitness_history(:, gen, :) = obj;
position_history(:, gen, :) = pop;

% Update Trajectories (using a random non-dominated solution)
if ~isempty(archive.pop)
    Trajectories(:, gen) = archive.pop(randi(size(archive.pop, 1)), :)';
end

disp(['Iteration #', num2str(gen), ' completed. Archive size: ', num2str(size(archive.pop, 1))]);
end

% Return Pareto front and Pareto set
paretoFront = archive.obj;
paretoSet = archive.pop;

end

% Helper Functions

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

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

function [pop, obj] = update_population(pop, obj, child, child_obj, index)
if dominates(child_obj, obj(index, :))
pop(index, :) = child;
obj(index, :) = child_obj;
elseif ~dominates(obj(index, :), child_obj)
if rand() < 0.5
pop(index, :) = child;
obj(index, :) = child_obj;
end
end
end

function archive = update_archive(archive_pop, archive_obj, pop, obj)
combined_pop = [archive_pop; pop];
combined_obj = [archive_obj; obj];


[is_dominated, ~] = pareto_dominance(combined_obj);
non_dominated = ~is_dominated;

archive.pop = combined_pop(non_dominated, :);
archive.obj = combined_obj(non_dominated, :);
end

function [is_dominated, domination_count] = pareto_dominance(objectives)
[n_solutions, n_objectives] = size(objectives);
is_dominated = false(n_solutions, 1);
domination_count = zeros(n_solutions, 1);


for i = 1:n_solutions
    for j = i+1:n_solutions
        if dominates(objectives(i, :), objectives(j, :))
            is_dominated(j) = true;
            domination_count(i) = domination_count(i) + 1;
        elseif dominates(objectives(j, :), objectives(i, :))
            is_dominated(i) = true;
            domination_count(j) = domination_count(j) + 1;
        end
    end
end
end

function result = dominates(obj1, obj2)
result = all(obj1 <= obj2) && any(obj1 < obj2);
end

function f = zdt1(x)
f1 = x(1);
g = 1 + 9 * sum(x(2:end)) / (length(x) - 1);
h = 1 - sqrt(f1 / g);
f2 = g * h;
f = [f1, f2];
end