أرغب في اعادة كتابة الكود التالي ليستخدم خوارزمية pso بدل الخوارزمية الجينية

clc; 
clear;
N = 300;
greedy_percent = 0.10;
random_percent = 0.90; 
alpha = 100; % 0.001
beta = 0.001;
generation = 1;
maximum_generation = 40;
num_elite_chromosmes = 1;
K_val = 4;
r = 0.8; 
coordinates = readcell('cordinates.csv');
[latitudes, longitudes] = convert_coordinates(coordinates(:,1), coordinates(:,2));
tbl_customers = readcell('c101customers35.csv');
tbl_customers(:, 2) = num2cell(latitudes);
tbl_customers(:, 3) = num2cell(longitudes);
writecell(tbl_customers, 'c101customers35.csv');
tbl_customers = readmatrix('c101customers35.csv');
total_time = tbl_customers(1, 6);
latitudes = tbl_customers(:, 2);
longitudes = tbl_customers(:, 3);
vertexs = [latitudes, longitudes];
n = size(vertexs, 1);
arr_distance_matrix = zeros(n);
for i = 1:n
    for j = 1:n
        arr_distance_matrix(i,j) = haversine_distance(vertexs(i,1), vertexs(i,2), vertexs(j,1), vertexs(j,2));
    end
end
columns_to_remove = [2, 3];
if width(tbl_customers) >= 8
    columns_to_remove = [columns_to_remove, 8];
end
if width(tbl_customers) >= 9
    columns_to_remove = [columns_to_remove, 9];
end
tbl_customers(:, columns_to_remove) = [];
arr_customers = tbl_customers;
num_customers = size(arr_customers, 1); % 26
number_random_chromosomes = floor(N*random_percent) - 1; 
number_greedy_chromosomes = N - number_random_chromosomes - 1; % 30 greedy chromosomes, N=300
number_nearest_nn_chromosomes = 1;
dataset_radius = (nanmax(arr_distance_matrix(:)) - nanmin(arr_distance_matrix(:)))/2;
total_capacity = 200;
all_customers_list = 1:num_customers;
radius = dataset_radius * (rand(1,1) * 0.5 + 0.5); 
[pop_random_routeList, pop_random_distanceList, pop_random_num_routes, pop_random_total_distance] = generate_random_strategy_routes(number_random_chromosomes, arr_customers, arr_distance_matrix, total_capacity, total_time, all_customers_list);
[pop_greedy_routeList, pop_greedy_distanceList,  pop_greedy_num_routes, pop_greedy_total_distance] = generate_greedy_strategy_routes(number_greedy_chromosomes, arr_customers, arr_distance_matrix, total_capacity, total_time, all_customers_list, radius);
[pop_nearest_routeList, pop_nearest_distanceList, pop_nearest_num_routes, pop_nearest_total_distance] = generate_nearest_neighbor_strategy_routes(number_nearest_nn_chromosomes, arr_customers, arr_distance_matrix, total_capacity, total_time, all_customers_list);
pop_chromosome_routeList_array = [pop_random_routeList; pop_greedy_routeList; pop_nearest_routeList];
pop_chromosome_distanceList_array = [pop_random_distanceList; pop_greedy_distanceList; pop_nearest_distanceList];
pop_num_routes = [pop_random_num_routes; pop_greedy_num_routes; pop_nearest_num_routes];
pop_total_distance = [pop_random_total_distance; pop_greedy_total_distance; pop_nearest_total_distance];
[pop_chromosome_routeList_array, pop_chromosome_distanceList_array, arr_pop_results] = initial_routing_phase2(N, alpha, beta, pop_chromosome_routeList_array, pop_chromosome_distanceList_array, arr_distance_matrix, arr_customers, total_time, total_capacity);
iter_time = [];
last_dist = 0;
stop_count = 0;
while generation <= maximum_generation
    tic;
    new_pop_chromosome_routeList_array = cell(N,1);
    new_pop_chromosome_distanceList_array = cell(N,1);
    arr_new_pop_results = zeros(N,3); % ['num_vehicles', 'distance', 'fitness']
    [arr_new_pop_results, new_pop_chromosome_routeList_array, new_pop_chromosome_distanceList_array] = full_recombination(alpha, beta, K_val,r, N,arr_pop_results, arr_new_pop_results, arr_customers, arr_distance_matrix, total_time, total_capacity, pop_chromosome_routeList_array, pop_chromosome_distanceList_array, new_pop_chromosome_routeList_array, new_pop_chromosome_distanceList_array);
    [new_pop_chromosome_routeList_array, new_pop_chromosome_distanceList_array, arr_new_pop_results] = mutation(alpha, beta, N, arr_new_pop_results, new_pop_chromosome_routeList_array, new_pop_chromosome_distanceList_array,arr_distance_matrix, arr_customers, total_time, total_capacity);
    [elite_chromosome_routeList, elite_chromosome_distanceList, elite_result] = elitism(pop_chromosome_routeList_array, pop_chromosome_distanceList_array, arr_pop_results);
    [~, worst_chromosome_number] = max(arr_new_pop_results(:,3));
    new_pop_chromosome_routeList_array{worst_chromosome_number} = elite_chromosome_routeList;
    new_pop_chromosome_distanceList_array{worst_chromosome_number} = elite_chromosome_distanceList;
    arr_new_pop_results(worst_chromosome_number,:) = elite_result(:);
    pop_chromosome_routeList_array = new_pop_chromosome_routeList_array;
    pop_chromosome_distanceList_array = new_pop_chromosome_distanceList_array;
    arr_pop_results = arr_new_pop_results;
    print_inf(pop_chromosome_routeList_array, pop_chromosome_distanceList_array, arr_pop_results, generation, vertexs,arr_distance_matrix);
    generation = generation + 1;
    [bestTotalDistance, ~] = min(arr_pop_results(:, 2));
    iter_time(generation) = toc;  
    if last_dist == bestTotalDistance
        stop_count = stop_count + 1;
        if stop_count > 5
            break;
        end
    else
        last_dist = bestTotalDistance;
        stop_count = 0;
    end 
end

function distance = distance_calculator(route, arr_distance_matrix)
    r = route(1:end-1);
    c = route(2:end);
    distance = sum(arr_distance_matrix(sub2ind(size(arr_distance_matrix), r, c)));
end

function arr_feasible_customers = feasible_neighbours(arr_customers, arr_distance_matrix, available_customers_list, total_time, total_capacity, curr_time, curr_capacity, curr_customer)
    num_available_customers = length(available_customers_list);
    arr_feasible_customers = zeros(num_available_customers, 11);
    arr_feasible_customers(:, 1:5) = arr_customers(available_customers_list, :);
    arr_feasible_customers(:, 2) = arr_feasible_customers(:, 2) + curr_capacity;
    arr_feasible_customers(:, 6) = arr_distance_matrix(curr_customer, available_customers_list);
    arr_feasible_customers(:, 7) = arr_feasible_customers(:, 6) + curr_time;
    arr_feasible_customers(:, 8) = max(0, arr_feasible_customers(:, 3) - arr_feasible_customers(:, 7));
    arr_feasible_customers(:, 9) = arr_feasible_customers(:, 7) + arr_feasible_customers(:, 8);
    arr_feasible_customers(:, 10) = arr_feasible_customers(:, 9) + arr_feasible_customers(:, 5);
    arr_feasible_customers(:, 11) = arr_distance_matrix(available_customers_list, 1) + arr_feasible_customers(:, 10);
    feasible_mask = (arr_feasible_customers(:, 2) <= total_capacity) ...
        & (arr_feasible_customers(:, 11) <= total_time) ...
        & (arr_feasible_customers(:, 9) >= arr_feasible_customers(:, 3)) ...
        & (arr_feasible_customers(:, 9) <= arr_feasible_customers(:, 4)) ;
    arr_feasible_customers = arr_feasible_customers(feasible_mask, :);
end


function [pop_chromosome_routeList_array, pop_chromosome_distanceList_array, pop_num_routes, pop_total_distance] = generate_greedy_strategy_routes(number_greedy_chromosomes, arr_customers, arr_distance_matrix, total_capacity, total_time, all_customers_list, radius)
    pop_chromosome_routeList_array = cell(number_greedy_chromosomes, 1);
    pop_chromosome_distanceList_array = cell(number_greedy_chromosomes, 1);
    pop_num_routes = zeros(number_greedy_chromosomes, 1);
    pop_total_distance = zeros(number_greedy_chromosomes, 1);
    for i = 1:number_greedy_chromosomes
        available_customers_list = setdiff(all_customers_list, 1);
        chromosome_route_list = {};
        chromosome_distance_list = [];
        route_list = [1];
        curr_distance = 0;
        curr_capacity = 0;
        curr_time = 0;
        curr_customer = 1;
        while ~isempty(available_customers_list)
            arr_feasible_neighbours = feasible_neighbours(arr_customers, arr_distance_matrix, available_customers_list, total_time, total_capacity, curr_time, curr_capacity, curr_customer);   
            if ~isempty(arr_feasible_neighbours)                
                if sum(arr_feasible_neighbours(:, 6) <= radius) ~= 0
                   arr_feasible_neighbours = arr_feasible_neighbours(arr_feasible_neighbours(:, 6) <= radius, :);
                end
                feasible_neighbours_list = arr_feasible_neighbours(:, 1);
                curr_customer = feasible_neighbours_list(randi(length(feasible_neighbours_list)));
                curr_customer_data = arr_feasible_neighbours(arr_feasible_neighbours(:,1) == curr_customer, :);
                route_list = [route_list, curr_customer];
                curr_distance = curr_distance + curr_customer_data(1, 6);
                curr_capacity = curr_customer_data(1, 2);
                curr_time = curr_customer_data(1, 10);
                available_customers_list = setdiff(available_customers_list, curr_customer);
            else
                route_list = [route_list, 1];
                curr_distance = curr_distance + arr_distance_matrix(curr_customer, 1);
                chromosome_route_list{end+1} = route_list;
                chromosome_distance_list = [chromosome_distance_list, curr_distance];
                route_list = [1];
                curr_distance = 0;
                curr_capacity = 0;
                curr_time = 0;
                curr_customer = 1;
            end
        end
        route_list = [route_list, 1];
        curr_distance = curr_distance + arr_distance_matrix(curr_customer, 1);
        chromosome_route_list{end+1} = route_list;
        chromosome_distance_list = [chromosome_distance_list, curr_distance];
        pop_chromosome_routeList_array{i} = chromosome_route_list;
        pop_chromosome_distanceList_array{i} = chromosome_distance_list;
        pop_num_routes(i) = length(chromosome_route_list);
        pop_total_distance(i) = sum(chromosome_distance_list);
    end
end


function [pop_chromosome_routeList_array, pop_chromosome_distanceList_array, pop_num_routes, pop_total_distance] = generate_nearest_neighbor_strategy_routes(number_nearest_nn_chromosomes, arr_customers, arr_distance_matrix, total_capacity, total_time, all_customers_list)
    pop_chromosome_routeList_array = cell(number_nearest_nn_chromosomes, 1);
    pop_chromosome_distanceList_array = cell(number_nearest_nn_chromosomes, 1);
    pop_num_routes = zeros(number_nearest_nn_chromosomes, 1);
    pop_total_distance = zeros(number_nearest_nn_chromosomes, 1);
    for i = 1:number_nearest_nn_chromosomes
        available_customers_list = setdiff(all_customers_list, 1);
        chromosome_route_list = {};
        chromosome_distance_list = [];
        route_list = [1];
        curr_distance = 0;
        curr_capacity = 0;
        curr_time = 0;
        curr_customer = 1;
        while ~isempty(available_customers_list)
            arr_feasible_neighbours = feasible_neighbours(arr_customers, arr_distance_matrix, available_customers_list, total_time, total_capacity, curr_time, curr_capacity, curr_customer);   
            if ~isempty(arr_feasible_neighbours) 
                [~, closest_idx] = min(arr_feasible_neighbours(:,6));
                curr_customer = arr_feasible_neighbours(closest_idx, 1); 
                curr_customer_data = arr_feasible_neighbours(arr_feasible_neighbours(:,1) == curr_customer, :);
                route_list = [route_list, curr_customer];
                curr_distance = curr_distance + curr_customer_data(1, 6);
                curr_capacity = curr_customer_data(1, 2);
                curr_time = curr_customer_data(1, 10);
                available_customers_list = setdiff(available_customers_list, curr_customer);
            else
                route_list = [route_list, 1];
                curr_distance = curr_distance + arr_distance_matrix(curr_customer, 1);
                chromosome_route_list{end+1} = route_list;
                chromosome_distance_list = [chromosome_distance_list, curr_distance];
                route_list = [1];
                curr_distance = 0;
                curr_capacity = 0;
                curr_time = 0;
                curr_customer = 1;
            end
        end
        route_list = [route_list, 1];
        curr_distance = curr_distance + arr_distance_matrix(curr_customer, 1);
        chromosome_route_list{end+1} = route_list;
        chromosome_distance_list = [chromosome_distance_list, curr_distance];
        pop_chromosome_routeList_array{i} = chromosome_route_list;
        pop_chromosome_distanceList_array{i} = chromosome_distance_list;
        pop_num_routes(i) = length(chromosome_route_list);
        pop_total_distance(i) = sum(chromosome_distance_list);
    end
end


function [pop_chromosome_routeList_array, pop_chromosome_distanceList_array, pop_num_routes, pop_total_distance] = generate_random_strategy_routes(number_random_chromosomes, arr_customers, arr_distance_matrix, total_capacity, total_time, all_customers_list)
    pop_chromosome_routeList_array = cell(number_random_chromosomes, 1);
    pop_chromosome_distanceList_array = cell(number_random_chromosomes, 1);
    pop_num_routes = zeros(number_random_chromosomes, 1);
    pop_total_distance = zeros(number_random_chromosomes, 1);
    for i = 1:number_random_chromosomes
        available_customers_list = setdiff(all_customers_list, 1);
        chromosome_route_list = {};
        chromosome_distance_list = [];
        route_list = [1];
        curr_distance = 0;
        curr_capacity = 0;
        curr_time = 0;
        curr_customer = 1;
        while ~isempty(available_customers_list)         
            arr_feasible_neighbours = feasible_neighbours(arr_customers, arr_distance_matrix, available_customers_list, total_time, total_capacity, curr_time, curr_capacity, curr_customer);   
            if ~isempty(arr_feasible_neighbours)              
                feasible_neighbours_list = arr_feasible_neighbours(:,1);
                curr_customer = datasample(feasible_neighbours_list, 1);                
                curr_customer_data = arr_feasible_neighbours(arr_feasible_neighbours(:,1) == curr_customer, :);
                route_list = [route_list, curr_customer];
                curr_distance = curr_distance + curr_customer_data(1, 6);
                curr_capacity = curr_customer_data(1, 2);
                curr_time = curr_customer_data(1, 10);
                available_customers_list = setdiff(available_customers_list, curr_customer);
            else
                route_list = [route_list, 1];
                curr_distance = curr_distance + arr_distance_matrix(curr_customer, 1);
                chromosome_route_list{end+1} = route_list;
                chromosome_distance_list = [chromosome_distance_list, curr_distance];   
                route_list = [1];
                curr_distance = 0;
                curr_capacity = 0;
                curr_time = 0;
                curr_customer = 1;
            end
        end        
        route_list = [route_list, 1];
        curr_distance = curr_distance + arr_distance_matrix(curr_customer, 1);
        chromosome_route_list{end+1} = route_list;
        chromosome_distance_list = [chromosome_distance_list, curr_distance];
        pop_chromosome_routeList_array{i} = chromosome_route_list;
        pop_chromosome_distanceList_array{i} = chromosome_distance_list;
        pop_num_routes(i) = length(chromosome_route_list);
        pop_total_distance(i) = sum(chromosome_distance_list);
    end
end

function [pop_chromosome_routeList_array, pop_chromosome_distanceList_array, arr_pop_results] = initial_routing_phase2(N, alpha, beta, pop_chromosome_routeList_array, pop_chromosome_distanceList_array, arr_distance_matrix, arr_customers, total_time, total_capacity)
    results_list_routeList = cell(N, 1);
    results_list_distanceList = cell(N, 1);
    results_num_vehicles = zeros(N,1);
    results_total_distance = zeros(N,1);
    results_fitness = zeros(N,1);
    arr_pop_results = zeros(N, 3);
    for i = 1:N
        [results_list_routeList{i}, results_list_distanceList{i}, results_num_vehicles(i), results_total_distance(i), results_fitness(i)] = routing_phase2(alpha, beta, pop_chromosome_routeList_array{i}, pop_chromosome_distanceList_array{i}, total_capacity, total_time, arr_distance_matrix, arr_customers);
    end
    for i = 1:N
        pop_chromosome_routeList_array{i} = results_list_routeList{i};
        pop_chromosome_distanceList_array{i} = results_list_distanceList{i};
        arr_pop_results(i, :) = [results_num_vehicles(i), results_total_distance(i), results_fitness(i)];
    end

end


function [pop_chromosome_routeList_array, pop_chromosome_distanceList_array, num_vehicles, total_distance, fitness] = routing_phase2(alpha, beta, pop_chromosome_routeList_array, pop_chromosome_distanceList_array, total_capacity, total_time, arr_distance_matrix, arr_customers)
    chr_index_list = 1:length(pop_chromosome_routeList_array)+1;
    chr_index_list(end) = 1;    
    for i = 1:length(chr_index_list)-1
        r1_index = chr_index_list(i);
        r2_index = chr_index_list(i+1);
        route_r1 = pop_chromosome_routeList_array{r1_index};
        route_r2 = pop_chromosome_routeList_array{r2_index};
        dist_r1 = pop_chromosome_distanceList_array(r1_index);
        dist_r2 = pop_chromosome_distanceList_array(r2_index);
        customer_r1 = route_r1(end-1);
        route_r2 = [route_r2(1), customer_r1, route_r2(2:end)];
        [route_validity_status, new_r2_dist] = route_validity(route_r2, arr_distance_matrix, arr_customers, total_time, total_capacity);  
        if route_validity_status == true
            if length(route_r1) > 3
                new_r1_dist = dist_r1 - arr_distance_matrix(route_r1(end-1),1) + arr_distance_matrix(route_r1(end-2),1);          
                if (dist_r1 + dist_r2) > (new_r1_dist + new_r2_dist)
                    route_r1(end-1) = [];    
                    pop_chromosome_routeList_array{r1_index} = route_r1;
                    pop_chromosome_routeList_array{r2_index} = route_r2;
                    pop_chromosome_distanceList_array(r1_index) = new_r1_dist;
                    pop_chromosome_distanceList_array(r2_index) = new_r2_dist;
                end
            else
                pop_chromosome_routeList_array(r1_index) = [];
                pop_chromosome_routeList_array{r2_index} = route_r2;
                pop_chromosome_distanceList_array(r1_index) = [];
                pop_chromosome_distanceList_array(r2_index) = new_r2_dist;
            end
        end
        total_distance = sum(pop_chromosome_distanceList_array);  
        num_vehicles = length(pop_chromosome_routeList_array);
        fitness = (alpha * num_vehicles) + (beta * total_distance);
    end
end


function [route_validity_status, curr_dist] = route_validity(route, arr_distance_matrix, arr_customers, total_time, total_capacity)
    route_validity_status = true;
    service_time = arr_customers(route(2), 5);
    curr_dist = arr_distance_matrix(route(1), route(2));
    curr_capacity = arr_customers(route(2), 2); 
    ready_time = arr_customers(route(2), 3);
    curr_time = arr_distance_matrix(route(1), route(2)); 
    waiting_time = max(0, ready_time - curr_time);
    curr_time = curr_time + waiting_time + service_time; 
    route_customers = route(2:end-1);
    for i = 1:length(route_customers)-1
        curr_dist = curr_dist + arr_distance_matrix(route_customers(i), route_customers(i+1));
        curr_time = curr_time + arr_distance_matrix(route_customers(i), route_customers(i+1));
        ready_time = arr_customers(route_customers(i+1), 3);
        due_time = arr_customers(route_customers(i+1), 4);
        curr_capacity = curr_capacity + arr_customers(route_customers(i+1), 2);
        if (curr_time<=due_time && curr_capacity<=total_capacity)
            waiting_time = max(0, ready_time - curr_time);
            if ((curr_time + waiting_time>=ready_time) && (curr_time + waiting_time + service_time <= total_time))
                curr_time = curr_time + waiting_time + service_time;
            else
                route_validity_status = false;
                break;
            end
        else
            route_validity_status = false;
            break;
        end
    end
    if route_validity_status == true
        curr_dist = curr_dist + arr_distance_matrix(route(end-1), route(end));
        curr_time = curr_time + arr_distance_matrix(route(end-1), route(end));
        if curr_time>total_time
            route_validity_status = false;
        end
    end
end

function [curr_chromosome, curr_chromosome_distances, curr_chromosome_total_distance] = insertion(RP, chromosome, chromosome_distances, arr_customers, arr_distance_matrix, total_capacity, total_time)
    curr_chromosome = chromosome;
    curr_chromosome_distances = chromosome_distances;
    for customer_to_insert = RP
        arr_tracker = cell(length(curr_chromosome), 6);
        arr_tracker(:, 1) = chromosome;
        arr_tracker(:, 2) = num2cell(chromosome_distances);
        arr_tracker(:, 3) = num2cell(cellfun(@(x) false, chromosome)); % false indicating customer not inserted yet
        arr_tracker(:, 4) = num2cell(cellfun(@(x) nan, chromosome));  % NaN indicating no route update yet
        for i = 1:length(chromosome)
            route = arr_tracker{i, 1};
            route_distance = arr_tracker{i, 2};
            [inserted_status, updated_route, updated_distance, update_route_distance_eval] = insertion_per_customer(customer_to_insert, route, route_distance, arr_customers, arr_distance_matrix, total_capacity, total_time);
            if inserted_status
                arr_tracker{i, 3} = true;
                arr_tracker{i, 4} = updated_route;
                arr_tracker{i, 5} = updated_distance;
                arr_tracker{i, 6} = update_route_distance_eval;
            end
        end
        updated_routes_indices = find(cell2mat(arr_tracker(:,3)));
        if ~isempty(updated_routes_indices)
            [~, min_route_idx] = min(cell2mat(arr_tracker(updated_routes_indices,6)));
            best_route_idx = updated_routes_indices(min_route_idx);
            curr_chromosome{best_route_idx} = arr_tracker{best_route_idx,4};
            curr_chromosome_distances(best_route_idx) = arr_tracker{best_route_idx,5};
        else
            additional_route = [1, customer_to_insert, 1];
            curr_chromosome{end+1} = additional_route;
            curr_chromosome_distances(end+1) = distance_calculator(additional_route, arr_distance_matrix);
        end

        chromosome = curr_chromosome;
        chromosome_distances = curr_chromosome_distances;
        curr_chromosome_total_distance = sum(curr_chromosome_distances);
    end
end

function [inserted_status, updated_route, updated_distance, update_route_distance_eval] = insertion_per_customer(customer_to_insert, route, route_distance, arr_customers, arr_distance_matrix, total_capacity, total_time)
    inserted_status = false;
    arr_positions_to_insert = 2:length(route);
    update_route_distance_eval = Inf;
    updated_route = route;
    updated_distance = route_distance;
    for positions_to_insert = arr_positions_to_insert
        new_route = [route(1:positions_to_insert-1) customer_to_insert route(positions_to_insert:end)];
        [route_validity_status, curr_dist] = route_validity(new_route, arr_distance_matrix, arr_customers, total_time, total_capacity);
        if route_validity_status
            inserted_status = true;
            distance_eval = curr_dist - route_distance;
            if distance_eval < update_route_distance_eval
                updated_route = new_route;
                update_route_distance_eval = distance_eval;
                updated_distance = curr_dist;
            end
        end
    end
end


function [pop_chromosome_routeList_array, pop_chromosome_distanceList_array, arr_pop_results] = mutation(alpha, beta, N, arr_pop_results, pop_chromosome_routeList_array, pop_chromosome_distanceList_array, arr_distance_matrix, arr_customers, total_time, total_capacity)
    for i = 1:N
        [route_validity_status, updated_chromosome, updated_chromosome_distance] = mutation_chromosome(pop_chromosome_routeList_array{i}, pop_chromosome_distanceList_array{i}, arr_distance_matrix, arr_customers, total_time, total_capacity);
        if route_validity_status
            pop_chromosome_routeList_array{i} = updated_chromosome;
            pop_chromosome_distanceList_array{i} = updated_chromosome_distance;
            arr_pop_results(i, 1) = length(updated_chromosome_distance);
            arr_pop_results(i, 2) = sum(updated_chromosome_distance);
            arr_pop_results(i, 3) = alpha * length(updated_chromosome_distance) + beta * sum(updated_chromosome_distance);
        end
    end
end


function [route_validity_status, updated_chromosome, updated_chromosome_distance] = mutation_chromosome(chromosome, chromosome_distance, arr_distance_matrix, arr_customers, total_time, total_capacity, service_time)    
    route_validity_status = false;
    list_len_routes = cellfun(@length, chromosome);
    feasible_indices = find(list_len_routes >= 5);
    if ~isempty(feasible_indices) % Check if there are feasible routes
        if length(feasible_indices) == 1
           rand_route_no = feasible_indices;
        else
           rand_route_no = randsample(feasible_indices, 1);
        end
        route = chromosome{rand_route_no};
        len_rand_route = length(route);
        mutation_len_choices = [2, 3];
        mutation_len = randsample(mutation_len_choices, 1);
        start_index = 2;     
        if len_rand_route - 2 > mutation_len 
            start_index = randi([2, len_rand_route - mutation_len]);
        end
        temp_list = route(start_index:start_index + mutation_len - 1);
        reverse_temp = fliplr(temp_list);
        route(start_index:start_index + mutation_len - 1) = reverse_temp;
        [route_validity_status, dist] = route_validity(route, arr_distance_matrix, arr_customers, total_time, total_capacity);
        if route_validity_status && dist <= chromosome_distance(rand_route_no)
            chromosome{rand_route_no} = route;
            chromosome_distance(rand_route_no) = dist;
        end
    end    
    updated_chromosome = chromosome;
    updated_chromosome_distance = chromosome_distance;
end


function [chromosome1_routes, chromosome2_routes, chromosome1_distances, chromosome2_distances, chromosome1_total_distance, chromosome2_total_distance] = recombination(K_val, r, arr_pop_results, N, pop_chromosome_routeList_array,pop_chromosome_distanceList_array, arr_customers, arr_distance_matrix, total_capacity, total_time)
    parent_chromosome_no1 = tournament_selection(K_val, r, arr_pop_results, N);
    parent_chromosome_no2 = tournament_selection(K_val, r, arr_pop_results, N);
    while parent_chromosome_no1 == parent_chromosome_no2
        parent_chromosome_no2 = tournament_selection(K_val, r, arr_pop_results, N);
    end
    chromosome1_routes = pop_chromosome_routeList_array{parent_chromosome_no1};
    chromosome2_routes = pop_chromosome_routeList_array{parent_chromosome_no2};
    chromosome1_distances = pop_chromosome_distanceList_array{parent_chromosome_no1};
    chromosome2_distances = pop_chromosome_distanceList_array{parent_chromosome_no2};
    chromosome1_total_distance = arr_pop_results(parent_chromosome_no1, 2);
    chromosome2_total_distance = arr_pop_results(parent_chromosome_no2, 2);
    rand_route_no1 = randi([1, length(chromosome1_routes)]);
    rand_route_no2 = randi([1, length(chromosome2_routes)]);
    RP1 = chromosome1_routes{rand_route_no1}(2:end-1);
    RP2 = chromosome2_routes{rand_route_no2}(2:end-1);
    [chromosome1_routes, chromosome1_distances] = removal_crossover_customers(RP2, chromosome1_routes, chromosome1_distances, arr_distance_matrix);
    [chromosome2_routes, chromosome2_distances] = removal_crossover_customers(RP1, chromosome2_routes, chromosome2_distances, arr_distance_matrix);
    [chromosome1_routes, chromosome1_distances, chromosome1_total_distance] = insertion(RP2, chromosome1_routes, chromosome1_distances, arr_customers, arr_distance_matrix, total_capacity, total_time);
    [chromosome2_routes, chromosome2_distances, chromosome2_total_distance] = insertion(RP1, chromosome2_routes, chromosome2_distances, arr_customers, arr_distance_matrix, total_capacity, total_time);

end


function [chromosome, chromosome_distances] = removal_crossover_customers(RP, chromosome, chromosome_distances, arr_distance_matrix)
    for i = 1:length(RP)
        customer = RP(i);
        for j = 1:length(chromosome)
            if ismember(customer, chromosome{j})
                if length(chromosome{j}) > 3
                    chromosome{j}(chromosome{j} == customer) = [];
                    chromosome_distances(j) = distance_calculator(chromosome{j}, arr_distance_matrix);
                    break; 
                elseif length(chromosome{j}) == 3
                    chromosome(j) = [];
                    chromosome_distances(j) = [];
                    j = j - 1; 
                    break; 
                end
            end
        end
    end
end


function parent_chromosome_no = tournament_selection(K_val, r, arr_pop_results, N)
    rand_chromosome_list = randperm(N, K_val);
    arr_tournament_set = arr_pop_results(rand_chromosome_list, :);    
    ran = rand();
    if ran < r
        [~, idx] = min(arr_tournament_set(:, 2));
        parent_chromosome_no = rand_chromosome_list(idx);
    else
        parent_chromosome_no = rand_chromosome_list(randi(K_val));
    end
end


function [elite_chromosome_routeList, elite_chromosome_distanceList, elite_result] = elitism(pop_chromosome_routeList_array, pop_chromosome_distanceList_array, arr_pop_results)
    [~, min_fitness_chromosome_index] = min(arr_pop_results(:, 3));
    elite_chromosome_routeList = pop_chromosome_routeList_array{min_fitness_chromosome_index};
    elite_chromosome_distanceList = pop_chromosome_distanceList_array{min_fitness_chromosome_index};
    elite_result = arr_pop_results(min_fitness_chromosome_index, :);

end


function print_inf(pop_chromosome_routeList_array, pop_chromosome_distanceList_array, arr_pop_results, gen,vertexs,arr_distance_matrix)   
    [bestFitness, bestIndex] = min(arr_pop_results(:, 3));
    bestRoute = pop_chromosome_routeList_array{bestIndex};
    bestDistance = pop_chromosome_distanceList_array{bestIndex};
     
    num_routes = numel(bestRoute);
    initial_routes = cell(num_routes, 1);
    for i = 1:num_routes
        initial_routes{i, 1} = bestRoute{i};
    end   
    chromosome = [];

    for i = 1 : num_routes
            chromosome = [chromosome, initial_routes{i},0];
    end
    best_routes = decodeChromosome(chromosome, num_routes);
    disp(['Generations: ', num2str(gen)]);
    fprintf('Best Fitness: %f\n', bestFitness);
    fprintf('Best Total Distance: %f\n', min(arr_pop_results(:, 2)));
    fprintf('Best Number of Tracks: %d\n', min(arr_pop_results(:, 1)));
    for i = 1:length(bestDistance)
        fprintf('Corresponding Distance of Each Route: %f, [', bestDistance(i));
        fprintf('%d ', bestRoute{i});
        fprintf(']\n');
    end
    drawMap(best_routes,vertexs,arr_distance_matrix);
    pause(0.5);
    fprintf('\n===================================================================================\n\n');
end

function [arr_new_pop_results, new_pop_chromosome_routeList_array, new_pop_chromosome_distanceList_array] = full_recombination(alpha, beta, K_val, r, N, arr_pop_results, arr_new_pop_results, arr_customers, arr_distance_matrix, total_time, total_capacity, pop_chromosome_routeList_array, pop_chromosome_distanceList_array, new_pop_chromosome_routeList_array, new_pop_chromosome_distanceList_array) 
    for i = 1:floor(N/2)
        [chromosome1_routes, chromosome2_routes, chromosome1_distances, chromosome2_distances, chromosome1_total_distance, chromosome2_total_distance] = recombination(K_val, r, arr_pop_results, N, pop_chromosome_routeList_array, pop_chromosome_distanceList_array, arr_customers, arr_distance_matrix, total_capacity, total_time);
        new_pop_chromosome_routeList_array{i*2 - 1} = chromosome1_routes;
        new_pop_chromosome_routeList_array{i*2} = chromosome2_routes;  
        new_pop_chromosome_distanceList_array{i*2 - 1} = chromosome1_distances;
        new_pop_chromosome_distanceList_array{i*2} = chromosome2_distances;  
        arr_new_pop_results(i*2 - 1, :) = [length(chromosome1_routes), chromosome1_total_distance, alpha*length(chromosome1_routes) + beta*chromosome1_total_distance];
        arr_new_pop_results(i*2, :) = [length(chromosome2_routes), chromosome2_total_distance, alpha*length(chromosome2_routes) + beta*chromosome2_total_distance];
    end
    
end

function [lat_dec, lon_dec] = convert_coordinates(lat, lon)
    n = length(lat);
    lat_dec = zeros(n, 1);
    lon_dec = zeros(n, 1);
    for i = 1:n
        lat_parts = sscanf(lat{i}, '%d° %d''%d"N');
        lon_parts = sscanf(lon{i}, '%d° %d''%d"E');
        lat_dec(i) = lat_parts(1) + lat_parts(2)/60 + lat_parts(3)/3600;
        lon_dec(i) = lon_parts(1) + lon_parts(2)/60 + lon_parts(3)/3600;
    end
end

function d = haversine_distance(lat1, lon1, lat2, lon2)
    R = 6371; 
    lat1 = deg2rad(lat1);
    lon1 = deg2rad(lon1);
    lat2 = deg2rad(lat2);
    lon2 = deg2rad(lon2);
    dlat = lat2 - lat1;
    dlon = lon2 - lon1;
    a = sin(dlat/2)^2 + cos(lat1) * cos(lat2) * sin(dlon/2)^2;
    c = 2 * atan2(sqrt(a), sqrt(1-a));
    d = R * c;
end

function [routes] = decodeChromosome(chromosome, num_routes)
    routes = cell(num_routes, 1);
    route_count = 1;
    split_locations = find(chromosome == 0); 
    for i = 1:length(split_locations)
        if i == 1
            current_route = chromosome(1:split_locations(i));
            current_route(current_route == chromosome(split_locations(i))) = [];
        else
            current_route = chromosome(split_locations(i - 1):split_locations(i));
            current_route(current_route == chromosome(split_locations(i - 1))) = [];
            current_route(current_route == chromosome(split_locations(i))) = [];
        end
        routes{route_count} = current_route; 
        route_count = route_count + 1;
    end
    if isempty(split_locations) % If there are no separation points, consider the entire chromosome as one route
        current_route = chromosome;
    end
    routes{route_count} = current_route;
end

function drawMap(best_routes, vertexs, arr_distance_matrix)
    clf;
    customer = vertexs(1:end,:);
    num_vehicules = size(best_routes, 1);
    hold on;
    box on;
    title('Optimal Roadmap')
    hold on;
    C = linspecer(num_vehicules); % give each route a unique color
    % Store colors for each route outside the loop
    route_colors = cell(num_vehicules, 1);
    for i = 1:num_vehicules
        route_colors{i} = C(i, :);
    end
    for i = 1:num_vehicules-1
        part_seq = best_routes{i};
        len = length(part_seq);
        color = route_colors{i}; % Get color for current route
        for j = 1:len
            if j == 1
                fprintf('Route %d:', i);
                fprintf('%d->', part_seq(1));
                c1 = customer(part_seq(2), :);
                plot([vertexs(1, 1), c1(1)], [vertexs(1, 2), c1(2)], '-', 'color', color, 'linewidth', 1);
                % Plot the distance
                mid_x = (vertexs(1, 1) + c1(1)) / 2;
                mid_y = (vertexs(1, 2) + c1(2)) / 2;
                text(mid_x, mid_y, num2str(arr_distance_matrix(1, part_seq(2))), 'Color', color, 'FontSize', 8, 'FontWeight', 'bold');
            elseif j == len
                fprintf('%d->', part_seq(j));
                fprintf('\n');
                c_len = customer(part_seq(len), :);
                plot([c_len(1), vertexs(1, 1)], [c_len(2), vertexs(1, 2)], '-', 'color', color, 'linewidth', 1);
                % Plot the distance
                mid_x = (c_len(1) + vertexs(1, 1)) / 2;
                mid_y = (c_len(2) + vertexs(1, 2)) / 2;
                text(mid_x, mid_y, num2str(arr_distance_matrix(part_seq(len), 1)), 'Color', color, 'FontSize', 8, 'FontWeight', 'bold');
            else
                fprintf('%d->', part_seq(j));
                c_pre = customer(part_seq(j), :);
                c_lastone = customer(part_seq(j + 1), :);
                plot([c_pre(1), c_lastone(1)], [c_pre(2), c_lastone(2)], '-', 'color', color, 'linewidth', 1);
                % Plot the distance
                mid_x = (c_pre(1) + c_lastone(1)) / 2;
                mid_y = (c_pre(2) + c_lastone(2)) / 2;
                text(mid_x, mid_y, num2str(arr_distance_matrix(part_seq(j), part_seq(j + 1))), 'Color', color, 'FontSize', 8, 'FontWeight', 'bold');
            end
        end
    end
    % Plot customer points
    plot(customer(:, 1), customer(:, 2), 'ro', 'linewidth', 1);
    % Add customer numbers
    for i = 2:size(customer, 1)  % Start from 2 to skip the depot
        text(customer(i, 1), customer(i, 2), num2str(i-1), 'HorizontalAlignment', 'center', 'VerticalAlignment', 'bottom', 'Color', 'k', 'FontWeight', 'bold');
    end
    hold on;
    % Plot depot
    plot(vertexs(1, 1), vertexs(1, 2), 's', 'linewidth', 2, 'MarkerEdgeColor', 'b', 'MarkerFaceColor', 'b', 'MarkerSize', 10);
    % Add depot label
    text(vertexs(1, 1), vertexs(1, 2), 'Depot', 'HorizontalAlignment', 'center', 'VerticalAlignment', 'top', 'Color', 'b', 'FontWeight', 'bold');
    
    % Enable zoom functionality
    zoom on;
end