I want to convert my NSGA 2 algorithm into matlab code.

To convert and implement the NSGA-II (Non-dominated Sorting Genetic Algorithm II) in MATLAB, you have two primary methods: use MATLAB's native multi-objective solver gamultiobj() (which is built on the NSGA-II algorithm with multi-threading and constraint handling), or construct a custom standalone NSGA-II script featuring fast non-dominated sorting, crowding distance assignment, Simulated Binary Crossover (SBX), and elitist selection.

Method 1: Native MATLAB NSGA-II Implementation via gamultiobj (Recommended)

MATLAB's Global Optimization Toolbox implements the NSGA-II framework directly inside gamultiobj(). Below is a complete script solving the standard bi-objective benchmark problem (Kursawe function):

% =========================================================================
% Bi-Objective Optimization using MATLAB's Native NSGA-II Solver
% =========================================================================
clc;
clear;
close all;

% Step 1: Define Multi-Objective Fitness Function Handle (2 Objectives)
% Objective 1: f1(x) = sum(-10 * exp(-0.2 * sqrt(x_i^2 + x_{i+1}^2)))
% Objective 2: f2(x) = sum(abs(x_i)^0.8 + 5 * sin(x_i^3))
fitnessFunction = @(x) [ ...
    sum(-10 * exp(-0.2 * sqrt(x(1:end-1).^2 + x(2:end).^2))), ...
    sum(abs(x).^0.8 + 5 * sin(x.^3)) ...
];

% Step 2: Problem Dimensions and Decision Variable Bounds
nvars = 3;                  % Number of design parameters
lb = -5 * ones(1, nvars);   % Lower bounds
ub =  5 * ones(1, nvars);   % Upper bounds

% Step 3: Configure NSGA-II Genetic Algorithm Options
options = optimoptions('gamultiobj', ...
    'PopulationSize', 120, ...              % Population size (N)
    'MaxGenerations', 150, ...              % Maximum generations
    'ParetoFraction', 0.4, ...              % Fraction of population on Pareto front
    'CrossoverFraction', 0.85, ...          % Crossover probability
    'PlotFcn', @gaplotpareto, ...           % Real-time Pareto front plotting
    'Display', 'iter', ...                  % Command window logging
    'UseParallel', false);                  % Set true for multi-core acceleration

% Step 4: Execute NSGA-II Optimization
fprintf('Starting NSGA-II Evolutionary Optimization...\n');
[x_opt, fval_opt, exitflag, output] = gamultiobj(fitnessFunction, nvars, [], [], [], [], lb, ub, [], options);

% Step 5: Plot Final Non-Dominated Pareto Optimal Front
figure('Name', 'NSGA-II Pareto Optimal Front', 'Color', 'w');
scatter(fval_opt(:,1), fval_opt(:,2), 45, 'filled', 'MarkerFaceColor', [0.85, 0.32, 0.09]);
grid on;
xlabel('Objective 1: f_1(x)', 'FontSize', 12, 'FontWeight', 'bold');
ylabel('Objective 2: f_2(x)', 'FontSize', 12, 'FontWeight', 'bold');
title(['NSGA-II Pareto Front (Total Solutions: ' num2str(size(fval_opt,1)) ')'], 'FontSize', 13, 'FontWeight', 'bold');

Method 2: Standalone Custom NSGA-II Algorithm from Scratch

If you want full algorithmic control without requiring the Global Optimization Toolbox, the following script executes the core NSGA-II loop:

% =========================================================================
% Custom NSGA-II Algorithm (Fast Non-Dominated Sorting + Crowding Distance)
% =========================================================================
clc;
clear;

% Algorithm Parameters
pop_size = 60;      % Population size
max_gen = 50;       % Generations
nvars = 2;          % Decision variables
lb = [0, 0];        % Lower bounds
ub = [1, 1];        % Upper bounds

% Step 1: Initialize Random Population
pop = lb + (ub - lb) .* rand(pop_size, nvars);

% Bi-Objective Function (ZDT1 Benchmark)
eval_obj = @(x) [x(1), (1 + 9*x(2)) * (1 - sqrt(x(1) / (1 + 9*x(2))))];

for gen = 1:max_gen
    % Evaluate Fitness for Current Population
    objs = zeros(pop_size, 2);
    for i = 1:pop_size
        objs(i, :) = eval_obj(pop(i, :));
    end
    
    % Step 2: Fast Non-Dominated Sorting
    ranks = non_dominated_sort(objs);
    
    % Step 3: Crowding Distance Calculation
    crowd_dist = calc_crowding_distance(objs, ranks);
    
    % Step 4: Binary Tournament Selection, Crossover (SBX), and Mutation
    offspring = zeros(pop_size, nvars);
    for i = 1:2:pop_size
        % Tournament selection
        p1 = tournament_select(ranks, crowd_dist);
        p2 = tournament_select(ranks, crowd_dist);
        
        % Simulated Binary Crossover (SBX)
        [c1, c2] = sbx_crossover(pop(p1,:), pop(p2,:), lb, ub);
        
        % Polynomial Mutation
        offspring(i, :)   = poly_mutation(c1, lb, ub);
        if i+1 <= pop_size
            offspring(i+1, :) = poly_mutation(c2, lb, ub);
        end
    end
    
    % Step 5: Elitist (2N -> N) Combination
    combined_pop = [pop; offspring];
    combined_objs = zeros(2 * pop_size, 2);
    for i = 1:(2 * pop_size)
        combined_objs(i, :) = eval_obj(combined_pop(i, :));
    end
    
    c_ranks = non_dominated_sort(combined_objs);
    c_dist = calc_crowding_distance(combined_objs, c_ranks);
    
    % Select Top N based on Rank and Crowding Distance
    [~, sort_idx] = sortrows([c_ranks, -c_dist], [1, 2]);
    pop = combined_pop(sort_idx(1:pop_size), :);
end

disp('NSGA-II Execution Completed Successfully.');

% =========================================================================
% Local Auxiliary Functions for NSGA-II
% =========================================================================
function ranks = non_dominated_sort(objs)
    N = size(objs, 1);
    ranks = zeros(N, 1);
    dominated_count = zeros(N, 1);
    dominates_list = cell(N, 1);
    
    for p = 1:N
        for q = 1:N
            if all(objs(p,:) <= objs(q,:)) && any(objs(p,:) < objs(q,:))
                dominates_list{p} = [dominates_list{p}, q];
            elseif all(objs(q,:) <= objs(p,:)) && any(objs(q,:) < objs(p,:))
                dominated_count(p) = dominated_count(p) + 1;
            end
        end
    end
    
    front = find(dominated_count == 0);
    current_rank = 1;
    while ~isempty(front)
        ranks(front) = current_rank;
        next_front = [];
        for i = 1:length(front)
            p = front(i);
            for j = 1:length(dominates_list{p})
                q = dominates_list{p}(j);
                dominated_count(q) = dominated_count(q) - 1;
                if dominated_count(q) == 0
                    next_front = [next_front, q];
                end
            end
        end
        current_rank = current_rank + 1;
        front = next_front;
    end
end

function dist = calc_crowding_distance(objs, ranks)
    N = size(objs, 1);
    dist = zeros(N, 1);
    num_obj = size(objs, 2);
    max_rank = max(ranks);
    
    for r = 1:max_rank
        idx = find(ranks == r);
        if length(idx) <= 2
            dist(idx) = inf;
            continue;
        end
        
        for m = 1:num_obj
            [sorted_vals, sort_order] = sort(objs(idx, m));
            dist(idx(sort_order(1))) = inf;
            dist(idx(sort_order(end))) = inf;
            
            val_range = sorted_vals(end) - sorted_vals(1);
            if val_range == 0, val_range = eps; end
            
            for k = 2:(length(idx) - 1)
                dist(idx(sort_order(k))) = dist(idx(sort_order(k))) + ...
                    (sorted_vals(k+1) - sorted_vals(k-1)) / val_range;
            end
        end
    end
end

function idx = tournament_select(ranks, dist)
    cand = randperm(length(ranks), 2);
    if ranks(cand(1)) < ranks(cand(2))
        idx = cand(1);
    elseif ranks(cand(2)) < ranks(cand(1))
        idx = cand(2);
    else
        idx = cand(1);
        if dist(cand(2)) > dist(cand(1)), idx = cand(2); end
    end
end

function [c1, c2] = sbx_crossover(p1, p2, lb, ub)
    eta_c = 20;
    u = rand(size(p1));
    beta = zeros(size(p1));
    beta(u <= 0.5) = (2 * u(u <= 0.5)).^(1 / (eta_c + 1));
    beta(u > 0.5)  = (1 ./ (2 * (1 - u(u > 0.5)))).^(1 / (eta_c + 1));
    c1 = min(max(0.5 * ((1 + beta).*p1 + (1 - beta).*p2), lb), ub);
    c2 = min(max(0.5 * ((1 - beta).*p1 + (1 + beta).*p2), lb), ub);
end

function mutant = poly_mutation(p, lb, ub)
    eta_m = 20;
    r = rand(size(p));
    delta = zeros(size(p));
    delta(r < 0.5) = (2 * r(r < 0.5)).^(1 / (eta_m + 1)) - 1;
    delta(r >= 0.5) = 1 - (2 * (1 - r(r >= 0.5))).^(1 / (eta_m + 1));
    mutant = min(max(p + delta .* (ub - lb), lb), ub);
end

Key NSGA-II Algorithmic Mechanisms

Practical Optimization Tips

• For Speed & Parallelization: Use gamultiobj() with 'UseParallel', true when your objective function evaluates complex Simulink models or FEA simulations.
• Vectorized Evaluation: Enable 'UseVectorized', true to evaluate entire populations in a single matrix operation rather than iterating one by one.
• Pareto Front Coverage: Set 'ParetoFraction' between 0.35 and 0.5 to ensure an even distribution of compromise solutions along conflicting objectives.