hey how can i include this switch MultiObjFnc into my moema code:


switch MultiObjFnc code:

% Multi-objective function
%MultiObjFnc = 'Schaffer';
%MultiObjFnc = 'Kursawe';
%MultiObjFnc = 'Poloni';
%MultiObjFnc = 'Viennet2';
%MultiObjFnc = 'Viennet3';
MultiObjFnc = 'ZDT1';
%MultiObjFnc = 'ZDT2';
%MultiObjFnc = 'ZDT3';
%MultiObjFnc = 'ZDT6';

switch MultiObjFnc
    case 'Schaffer'         % Schaffer
        MultiObj.fun = @(x) [x(:).^2, (x(:)-2).^2];
        MultiObj.nVar = 1;
        MultiObj.var_min = -5;
        MultiObj.var_max = 5;
        load('Schaffer.mat');
        MultiObj.truePF = PF;
    case 'Kursawe'          % Kursawe 
        MultiObj.fun = @(x) [-10.*(exp(-0.2.*sqrt(x(:,1).^2+x(:,2).^2)) + exp(-0.2.*sqrt(x(:,2).^2+x(:,3).^2))), ...
                             sum(abs(x).^0.8 + 5.*sin(x.^3),2)];
        MultiObj.nVar = 3;
        MultiObj.var_min = -5.*ones(1,MultiObj.nVar);
        MultiObj.var_max = 5.*ones(1,MultiObj.nVar);
        load('Kursawe.mat');
        MultiObj.truePF = PF;
    case 'Poloni'           % Poloni's two-objective
        A1 = 0.5*sin(1)-2*cos(1)+sin(2)-1.5*cos(2);
        A2 = 1.5*sin(1)-cos(1)+2*sin(2)-0.5*cos(2);
        B1 = @(x,y) 0.5.*sin(x)-2.*cos(x)+sin(y)-1.5.*cos(y);
        B2 = @(x,y) 1.5.*sin(x)-cos(x)+2.*sin(y)-0.5.*cos(y);
        f1 = @(x,y) 1+(A1-B1(x,y)).^2+(A2-B2(x,y)).^2;
        f2 = @(x,y) (x+3).^2+(y+1).^2;
        MultiObj.fun = @(x) [f1(x(:,1),x(:,2)), f2(x(:,1),x(:,2))];
        MultiObj.nVar = 2;
        MultiObj.var_min = -pi.*ones(1,MultiObj.nVar);
        MultiObj.var_max = pi.*ones(1,MultiObj.nVar);
    case 'Viennet2'         % Viennet2
        f1 = @(x,y) 0.5.*(x-2).^2+(1/13).*(y+1).^2+3;
        f2 = @(x,y) (1/36).*(x+y-3).^2+(1/8).*(-x+y+2).^2-17;
        f3 = @(x,y) (1/175).*(x+2.*y-1).^2+(1/17).*(2.*y-x).^2-13;
        MultiObj.fun = @(x) [f1(x(:,1),x(:,2)), f2(x(:,1),x(:,2)), f3(x(:,1),x(:,2))];
        MultiObj.nVar = 2;
        MultiObj.var_min = [-4, -4];
        MultiObj.var_max = [4, 4];
        load('Viennet2.mat');
        MultiObj.truePF = PF;
    case 'Viennet3'         % Viennet3
        f1 = @(x,y) 0.5.*(x.^2+y.^2)+sin(x.^2+y.^2);
        f2 = @(x,y) (1/8).*(3.*x-2.*y+4).^2 + (1/27).*(x-y+1).^2 +15;
        f3 = @(x,y) (1./(x.^2+y.^2+1))-1.1.*exp(-(x.^2+y.^2));
        MultiObj.fun = @(x) [f1(x(:,1),x(:,2)), f2(x(:,1),x(:,2)), f3(x(:,1),x(:,2))];
        MultiObj.nVar = 2;
        MultiObj.var_min = [-3, -10];
        MultiObj.var_max = [10, 3];
        load('Viennet3.mat');
        MultiObj.truePF = PF;
    case 'ZDT1'             % ZDT1 (convex)
        g = @(x) 1+9.*sum(x(:,2:end),2)./(size(x,2)-1);
        MultiObj.fun = @(x) [x(:,1), g(x).*(1-sqrt(x(:,1)./g(x)))];
        MultiObj.nVar = 30; 
        MultiObj.var_min = zeros(1,MultiObj.nVar);
        MultiObj.var_max = ones(1,MultiObj.nVar);
        load('ZDT1.mat');
        MultiObj.truePF = PF;
    case 'ZDT2'             % ZDT2 (non-convex)
        f = @(x) x(:,1);
        g = @(x) 1+9.*sum(x(:,2:end),2)./(size(x,2)-1);
        h = @(x) 1-(f(x)./g(x)).^2;
        MultiObj.fun = @(x) [f(x), g(x).*h(x)];
        MultiObj.nVar = 30; 
        MultiObj.var_min = zeros(1,MultiObj.nVar);
        MultiObj.var_max = ones(1,MultiObj.nVar);
        load('ZDT2.mat');
        MultiObj.truePF = PF;
    case 'ZDT3'             % ZDT3 (discrete)
        f = @(x) x(:,1);
        g  = @(x) 1+(9/size(x,2)-1).*sum(x(:,2:end),2);
        h  = @(x) 1 - sqrt(f(x)./g(x)) - (f(x)./g(x)).*sin(10.*pi.*f(x));
        MultiObj.fun = @(x) [f(x), g(x).*h(x)];
        MultiObj.nVar = 30;
        MultiObj.var_min = 0.*ones(1,MultiObj.nVar);
        MultiObj.var_max = 1.*ones(1,MultiObj.nVar);
        load('ZDT3.mat');
        MultiObj.truePF = PF;
    case 'ZDT6'             % ZDT6 (non-uniform)
        f = @(x) 1 - exp(-4.*x(:,1)).*sin(6.*pi.*x(:,1));
        g = @(x) 1 + 9.*(sum(x(:,2:end),2)./(size(x,2)-1)).^0.25;
        h = @(x) 1 - (f(x)./g(x)).^2;
        MultiObj.fun = @(x) [f(x), g(x).*h(x)];
        MultiObj.nVar = 10;
        MultiObj.var_min = 0.*ones(1,MultiObj.nVar);
        MultiObj.var_max = 1.*ones(1,MultiObj.nVar);
        load('ZDT6.mat');
        MultiObj.truePF = PF;
end


moema code:
% Main script to run MOEMA optimization
clc; 
clear; 
close all;
% Set problem parameters
NPop = 50;  % Population size
T = 50;     % 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