here is a code for comparing moema and moea-d. can you add mopso to the code and compare it with others.
moema and moea-d code:
% Main script to run MOEMA and MOEA/D 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_MOEMA, ~, ~, ~, ~] = MOEMA(NPop, T, lb, ub, dim, @zdt1);

% Run MOEA/D
nVar = dim;
VarMin = 0;
VarMax = 1;
nPop_MOEAD = 200;
nArchive = 50;
MaxIt = T;
T_MOEAD = 15;
CostFunction = @(x) ZDT(x);
crossover_params.gamma = 0.5;
crossover_params.VarMin = VarMin;
crossover_params.VarMax = VarMax;
[paretoFront_MOEAD, ~, ~, ~] = MOEAD(nVar, VarMin, VarMax, nPop_MOEAD, nArchive, MaxIt, T_MOEAD, CostFunction, crossover_params);

% Plot the Pareto fronts
figure;
plot(paretoFront_MOEMA(:,1), paretoFront_MOEMA(:,2), 'ro', 'DisplayName', 'MOEMA');
hold on;
plot(paretoFront_MOEAD(1,:), paretoFront_MOEAD(2,:), 'bx', 'DisplayName', 'MOEA/D');
xlabel('f1');
ylabel('f2');
title('Pareto Fronts for ZDT1 using MOEMA and MOEA/D');
legend;

% Generate true ZDT1 Pareto front
f1 = 0:0.01:1;
g = 1;
h = 1 - sqrt(f1 ./ g);
f2 = g .* h;
plot(f1, f2, 'k-', 'DisplayName', 'True Pareto Front');
legend;

% Create a table for comparison
MOEMA_solutions = size(paretoFront_MOEMA, 1);
MOEAD_solutions = size(paretoFront_MOEAD, 2);
MOEMA_min_f1 = min(paretoFront_MOEMA(:,1));
MOEAD_min_f1 = min(paretoFront_MOEAD(1,:));
MOEMA_min_f2 = min(paretoFront_MOEMA(:,2));
MOEAD_min_f2 = min(paretoFront_MOEAD(2,:));

table = array2table([MOEMA_solutions; MOEAD_solutions; MOEMA_min_f1; MOEAD_min_f1; MOEMA_min_f2; MOEAD_min_f2]);
table.Properties.VariableNames = {'Value'};
table.Properties.RowNames = {'MOEMA Solutions', 'MOEA/D Solutions', 'MOEMA Min f1', 'MOEA/D Min f1', 'MOEMA Min f2', 'MOEA/D Min f2'};
disp(table);

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

%% Main MOEA/D Function
function [paretoFront, paretoSet, fitness_history, position_history] = MOEAD(nVar, VarMin, VarMax, nPop, nArchive, MaxIt, T, CostFunction, crossover_params)

    % Problem Definition
    VarSize = [nVar, 1];   % Decision Variables Matrix Size
    nObj = numel(CostFunction(unifrnd(VarMin, VarMax, VarSize)));  % Number of Objectives

    % Create Sub-problems
    sp = CreateSubProblems(nObj, nPop, T);

    % Empty Individual Template
    empty_individual.Position = [];
    empty_individual.Cost = [];
    empty_individual.g = [];
    empty_individual.IsDominated = [];

    % Initialize Goal Point
    z = zeros(nObj, 1);

    % Create Initial Population
    pop = repmat(empty_individual, nPop, 1);
    for i = 1:nPop
        pop(i).Position = unifrnd(VarMin, VarMax, VarSize);
        pop(i).Cost = CostFunction(pop(i).Position);
        z = min(z, pop(i).Cost);  % Update goal point
    end

    for i = 1:nPop
        pop(i).g = DecomposedCost(pop(i), z, sp(i).lambda);
    end

    % Determine Population Domination Status
    pop = DetermineDomination(pop);

    % Initialize Estimated Pareto Front (EP)
    EP = pop(~[pop.IsDominated]);

    % History variables for fitness and positions
    fitness_history = zeros(nObj, MaxIt);
    position_history = cell(nPop, MaxIt);

    % Main Loop
    for it = 1:MaxIt
        for i = 1:nPop
            % Reproduction (Crossover)
            K = randsample(T, 2);
            
            j1 = sp(i).Neighbors(K(1));
            p1 = pop(j1);
            
            j2 = sp(i).Neighbors(K(2));
            p2 = pop(j2);
            
            y = empty_individual;
            y.Position = Crossover(p1.Position, p2.Position, crossover_params);
            
            y.Cost = CostFunction(y.Position);
            
            z = min(z, y.Cost);  % Update goal point
            
            % Update neighboring solutions
            for j = sp(i).Neighbors
                y.g = DecomposedCost(y, z, sp(j).lambda);
                if y.g <= pop(j).g
                    pop(j) = y;
                end
            end
        end
        
        % Determine Population Domination Status
        pop = DetermineDomination(pop);
        
        ndpop = pop(~[pop.IsDominated]);
        
        EP = [EP; ndpop]; %#ok<AGROW>
        EP = DetermineDomination(EP);
        EP = EP(~[EP.IsDominated]);
        
        if numel(EP) > nArchive
            Extra = numel(EP) - nArchive;
            ToBeDeleted = randsample(numel(EP), Extra);
            EP(ToBeDeleted) = [];
        end

        % Store fitness data and positions for this iteration
        fitness_history(:, it) = mean([pop.Cost], 2);  % Store average fitness
        for i = 1:nPop
            position_history{i, it} = pop(i).Position;
        end
        
        % Display Iteration Information
        disp(['Iteration ' num2str(it) ': Number of Pareto Solutions = ' num2str(numel(EP))]);
    end

    % Prepare Results
    paretoFront = [EP.Cost];
    paretoSet = [EP.Position];

end

%% Helper Functions

% ZDT1 Cost Function
function z = ZDT(x)
    n = numel(x);
    f1 = x(1);
    g = 1 + 9 / (n-1) * sum(x(2:end));
    h = 1 - sqrt(f1 / g);
    f2 = g * h;
    z = [f1; f2];
end

% Plot Pareto Front
function PlotCosts(paretoFront)
    plot(paretoFront(1,:), paretoFront(2,:), 'x');
    xlabel('1^{st} Objective');
    ylabel('2^{nd} Objective');
    grid on;
end

% Check if one solution dominates another
function b = Dominates(x, y)
    if isfield(x, 'Cost')
        x = x.Cost;
    end
    if isfield(y, 'Cost')
        y = y.Cost;
    end
    b = all(x <= y) && any(x < y);
end

% Determine Domination in the Population
function pop = DetermineDomination(pop)
    nPop = numel(pop);
    for i = 1:nPop
        pop(i).IsDominated = false;
    end
    for i = 1:nPop
        for j = i+1:nPop
            if Dominates(pop(i), pop(j))
                pop(j).IsDominated = true;
            elseif Dominates(pop(j), pop(i))
                pop(i).IsDominated = true;
            end
        end
    end
end

% Decomposed Cost Calculation using Tchebycheff Approach
function g = DecomposedCost(individual, z, lambda)
    if isfield(individual, 'Cost')
        fx = individual.Cost;
    else
        fx = individual;
    end
    g = max(lambda .* abs(fx - z));
end

% Crossover Operation
function y = Crossover(x1, x2, params)
    gamma = params.gamma;
    VarMin = params.VarMin;
    VarMax = params.VarMax;
    alpha = unifrnd(-gamma, 1+gamma, size(x1));
    y = alpha .* x1 + (1 - alpha) .* x2;
    y = min(max(y, VarMin), VarMax);
end

% Create Sub-problems
function sp = CreateSubProblems(nObj, nPop, T)
    empty_sp.lambda = [];
    empty_sp.Neighbors = [];
    sp = repmat(empty_sp, nPop, 1);

    for i = 1:nPop
        lambda = rand(nObj, 1);
        lambda = lambda / norm(lambda);
        sp(i).lambda = lambda;
    end

    LAMBDA = [sp.lambda]';
    D = pdist2(LAMBDA, LAMBDA);  % Distance between lambda vectors

    for i = 1:nPop
        [~, SO] = sort(D(i, :));
        sp(i).Neighbors = SO(1:T);
    end
end


mopso code:
% Main script to run MOPSO on ZDT1

% Parameters
params.Np = 200;        % Population size
params.Nr = 100;        % Repository size
params.maxgen = 100;    % Maximum number of generations
params.W = 0.4;         % Inertia weight
params.C1 = 2;          % Individual confidence factor
params.C2 = 2;          % Swarm confidence factor
params.ngrid = 20;      % Number of grids in each dimension
params.maxvel = 5;      % Maximum velocity in percentage
params.u_mut = 0.5;     % Uniform mutation percentage

% Define ZDT1 problem
MultiObj.nVar = 30;
MultiObj.var_min = zeros(MultiObj.nVar, 1);
MultiObj.var_max = ones(MultiObj.nVar, 1);
MultiObj.fun = @(POS) cell2mat(arrayfun(@(i) zdt1(POS(i,:)), 1:size(POS, 1), 'UniformOutput', false)');

% Generate true Pareto front
MultiObj.truePF = generateZDT1TruePF(1000);

% Run MOPSO
REP = MOPSO(params, MultiObj);

% Display results
disp('Repository fitness values (POS_fit) and positions (POS):');
display(REP.pos_fit);
display(REP.pos);

% Plot final results with the true Pareto front
figure;
plot(REP.pos_fit(:,1), REP.pos_fit(:,2), 'ro', 'DisplayName', 'MOPSO Solutions');
hold on;
plot(MultiObj.truePF(:,1), MultiObj.truePF(:,2), 'b-', 'DisplayName', 'True Pareto Front');
xlabel('f1');
ylabel('f2');
title('MOPSO Solutions vs True Pareto Front for ZDT1');
legend;
grid on;

% ZDT1 Objective Function
function y = zdt1(x)
    n = length(x);
    f1 = x(1);
    g = 1 + (9/(n - 1)) * sum(x(2:end));
    f2 = g * (1 - sqrt(f1 / g));
    y = [f1, f2];
end

% Generate True Pareto Front for ZDT1
function truePF = generateZDT1TruePF(numPoints)
    f1 = linspace(0, 1, numPoints)';
    f2 = 1 - sqrt(f1);
    truePF = [f1, f2];
end

% MOPSO Function
function REP = MOPSO(params, MultiObj)
    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;
    DOMINATED = checkDomination(POS_fit);
    REP.pos = POS(~DOMINATED, :);
    REP.pos_fit = POS_fit(~DOMINATED, :);
    REP = updateGrid(REP, ngrid);
    maxvel = (var_max - var_min) .* maxvel ./ 100;
    gen = 1;

    % Main MOPSO loop
    while gen <= maxgen
        % Select leader
        h = selectLeader(REP);
        
        % Update speeds 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 the population
        POS_fit = fun(POS);
        
        % Update the 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 the best positions found so far for each particle
        pos_best = dominates(POS_fit, PBEST_fit);
        best_pos = ~dominates(PBEST_fit, POS_fit);
        best_pos(rand(Np, 1) >= 0.5) = 0;
        PBEST_fit(pos_best, :) = POS_fit(pos_best, :);
        PBEST(pos_best, :) = POS(pos_best, :);
        
        % Update generation
        gen = gen + 1;
    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), :));
    dominated_particles = unique(all_perm(d == 1, 2));
    dom_vector(dominated_particles) = 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)
        [m_fit, idx] = sort(REP.pos_fit(:, m), 'ascend');
        m_up = [m_fit(2:end); Inf];
        m_down = [Inf; m_fit(1:end-1)];
        distance = (m_up - m_down) ./ (max(m_fit) - min(m_fit));
        [~, idx] = sort(idx, 'ascend');
        crowding = crowding + distance(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), round(Np/3)];
    else
        sub_sizes = [round(Np/3), round(Np/3), floor(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