% 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 = 100;       % Population size
maxGen = 2000;       % Maximum number of generations
mutationRate = 0.05; % 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

% 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
    
    % 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

% Final Solution
bestSolution = pop(1, :);
bestFitness = fitnessFunction(bestSolution, t, a, c, V, travelTime);

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

% 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