% Genetic Algorithm for Minimizing Network Traversal Time

clc;
clear;
close all;

% Problem Parameters
N = 17;                                                % Number of roads
t = 1.5 * ones(1, N);                                  % Fixed time for each road
a = [1.25 * ones(1, 5), 1.5 * ones(1, 5), ones(1, 7)]; % Weighting factor
c = [
    54.13, 21.56, 34.08, 49.19, 33.03, 21.84, 29.96, 24.87, 47.24, 33.97, ...
    26.89, 32.76, 39.98, 37.12, 53.83, 61.65, 59.73];  % Road capacities
V = 100;                                               % Incoming vehicle rate

% Travel Time Function
travelTime = @(xi, ti, ai, ci) ti + ai * xi / (1 - xi / ci);
% Normalization Function (Infinite norm normalized to S)
normalizeSum  = @(x, S) (x ./ sum(x)) * S;     % Ensure sum of single row x equals S
normalizeSum2 = @(x, S) (x ./ sum(x, 2)) * S;  % Ensure sum of each row of 2D matrix x equals S

% Genetic Algorithm Parameters
popSize = 36;       % Population size
maxGen = 2000;       % Maximum number of generations
mutationRate = 0.1; % Mutation probability

% Initialize Population
pop = rand(popSize, N) .* c;    % Random initial solutions (0 <= x <= c)
pop = normalizeSum2(pop, V);    % Ensure sum of each solution equals V

newPop = zeros(popSize, N);     % Pre-allocate new population buffer
bestFitness = zeros(maxGen, 1); % Result array

% Genetic Algorithm Execution
for gen = 1:maxGen
    % Fitness Calculation
    fitness = arrayfun(@(i) fitnessFunction(pop(i, :), t, a, c, V, travelTime), 1:popSize);

    % Selection
    [~, idx] = sort(fitness); % Sort based on fitness (ascending order)
    pop = pop(idx, :);        % Retain the best solutions
    
    % Keep the best chromosome
    bestFitness(gen) = fitnessFunction(pop(1, :), t, a, c, V, travelTime);

    % Crossover
    newPop(1:popSize/2, :) = pop(1:popSize/2, :); % Retain top half
    for i = 1:popSize/2
        parent1 = newPop(randi(popSize/2), :);
        parent2 = newPop(randi(popSize/2), :);
        crossPoint = randi(N);
        child = [parent1(1:crossPoint), parent2(crossPoint+1:end)];
        child = normalizeSum(child, V);
        newPop(popSize/2 + i, :) = child;
    end
    
    % Mutation
    for i = 1:popSize
        if rand < mutationRate
            mutationIdx = randi(N);
            newPop(i, mutationIdx) = rand * c(mutationIdx);
            newPop(i, :) = normalizeSum(newPop(i, :), V);
        end
    end
    
    % Replacement
    pop = newPop;
end

% Results
bestSolution = pop(1, :);
disp('Best Solution [veh/min]:');
disp(bestSolution);
disp(['Best Objective Value: ', num2str(bestFitness(end)), ' [min]']);

figure('Name', 'Time over generations', 'NumberTitle', 'off');
set(gcf, 'Position', [100, 100, 960, 640]); % Set the figure size
plot(1:maxGen, bestFitness, '-b', 'LineWidth', 1);

% Customize the plot
title(['Population = ', num2str(popSize), '   -   Mutation = ', num2str(mutationRate)], 'Interpreter', 'latex', 'FontSize', 16); % Title of the plot
xlabel('Generations') ;
ylabel('T_{total}');

% save the figure
print(gcf, ['figures/constV_pop_', num2str(popSize), 'mut_', num2str(mutationRate), '.png'], '-dpng', '-r300');


% Fitness Function
function T_total = fitnessFunction(x, t, a, c, V, travelTime)
    if abs(sum(x) - V) > 1e-6 || any(x < 0) || any(x > c)
        T_total = inf; % Infeasible solutions
        return;
    end
    
    T = arrayfun(@(xi, ti, ai, ci) travelTime(xi, ti, ai, ci), x, t, a, c); % Apply function to all elements
    T_total = sum(T .* x); % Total traversal time
end