遗传算法求解带时间窗的车辆路径问题 (VRPTW)

%clear, clc, 
%close all;

%% input data
% test_date = importdata('test.xlsx');
test_date = importdata('c101.txt');
% tw1 means the earliest start time, tw2 means the latest start time
depot_time_window1 = test_date(1,5);
depot_time_window2 = test_date(1,6);
time_window1 = test_date(2:end, 5); % tw, start time
time_window2 = test_date(2:end, 6); % tw, end time
service_time = test_date(2:end, 7); % service time 
vertexs = test_date(:, 2:3);
customer_location = vertexs(2:end, :);
customer_number = size(customer_location, 1);

robot_number = 5; % number of vehicles

% dist matrix 距离矩阵
distance_pair = pdist(vertexs); % pairwise distance
distance_matrix = squareform(distance_pair); 

%% GA parameters 
population = 100;            %种群内个体数（一般20-100个）
generation = 1;              %开始种群代数
maximum_generation = 50;     %终止进化代数（一般100-500代）
probability_crossover = 0.9; %交叉概率（一般0.4-0.99）
probability_mutation = 0.05; %变异概率（一般0.0001-0.1）
rate_gap = 0.9; 
length_chrom = customer_number + robot_number - 1;  %染色体长度：染色体长度为顾客数量加上机器人数量减一，表示所有顾客和机器人的访问顺序。

%% initialize population
init_chrom = createInitChrom(customer_number, time_window1);
%调用 createInitChrom 函数生成一个初始染色体 init_chrom。
chroms = createInitPopulation(population, length_chrom, customer_number, init_chrom);
%调用 createInitPopulation 函数生成一个大小为 population 的种群 chroms，其中每个染色体的长度为 length_chrom，初始染色体为 init_chrom。

% display some values of the initial random population 这段代码用于显示初始随机种群中某个染色体的解码结果
% 调用 decode 函数将染色体解码成多条路径，并计算路径的总行驶距离和违规情况。然后，它会将解码结果输出到控制台，包括车辆数、总行驶距离、违规路径数和违规顾客数等信息。
[VC,NV,TD,violate_num,violate_cus] = decode(chroms(1,:),customer_number,time_window1,time_window2,depot_time_window2,service_time,distance_matrix);
disp('The initial population:');
disp(['Number of Robots: ', num2str(NV), ', Total Distance: ', num2str(TD), ', Number of violated Path: ', num2str(violate_num), ', Number of Violated Customer: ', num2str(violate_cus)]);
fprintf('\n');

% calsulate the value of objective function 
%调用 calObj 函数来计算每个染色体的目标函数值，并将计算结果存储在 ObjV 变量中
ObjV = calObj(chroms,customer_number,time_window1,time_window2,depot_time_window2,service_time,distance_matrix);
% the minimal value
%使用 min 函数找到 ObjV 中的最小值，存储在 pre_objective_value 变量中。
pre_objective_value=min(ObjV);
% output the objective function value of current generation
disp(['Generation ', num2str(generation), ': Objective Function Value = ', num2str(min(ObjV))]);


 % Display the change of the objective function value     画图一：最优值迭代图
figure;                         % figure 函数用于创建一个新的图形窗口
hold on;   %函数用于保持当前图形窗口,并允许在其上添加新的图形元素
box on;    %用于在图形窗口中显示坐标轴的边框。
xlim([0, maximum_generation]); %用于设置横轴的范围
title('Optimization Process');  %优化过程
xlabel('Generation');
ylabel('The Optimal Value'); %最优值

%定义了一些变量，用于记录迭代次数、上一次的目标函数值和停止计数器
iter_time = [];
last_dist = 0;
stop_count = 0;


%% the loop
while generation <= maximum_generation
    
    tic;
    
    % calculate the fitness value
    ObjV = calObj(chroms,customer_number,time_window1,time_window2,depot_time_window2,service_time,distance_matrix);
    % draw the line
%     line([generation - 1, generation], [pre_objective_value, min(ObjV)]); pause(0.0001);
    pre_objective_value = min(ObjV);
    FitnV = Fitness(ObjV);
    % select
    SelCh = Select(chroms,FitnV,rate_gap);
    % crossover
    SelCh = Recombin(SelCh,probability_crossover);
    % mutation
    SelCh = Mutate(SelCh,probability_mutation);
    
    
    % local search
    SelCh = neighborhoodSearch(SelCh, customer_number, time_window1, time_window2, depot_time_window2, service_time, distance_matrix);
    % 
    chroms = Reins(chroms,SelCh,ObjV);
    % Delete duplicate chromsomes 
    chroms = DealRepeat(chroms);
    ObjV = calObj(chroms,customer_number,time_window1,time_window2,depot_time_window2,service_time,distance_matrix);
    [minObjV,minInd]=min(ObjV);
    disp(['Generation:' num2str(generation)]);
    [bestVC,bestNV,bestTD,best_vionum,best_viocus]=decode(chroms(minInd(1),:),customer_number,time_window1,time_window2,depot_time_window2,service_time,distance_matrix);
    disp(['Number of Robots: ', num2str(bestNV), ',Objective Function Value : ', num2str(min(ObjV)),', Total Distance: ', num2str(bestTD), ', Number of violated Path: ', num2str(best_vionum), ', Number of Violated Customer: ', num2str(best_viocus)]);
    fprintf('\n');
    generation = generation + 1;
    
    iter_time(generation) = toc;
    
    if last_dist == bestTD
        stop_count = stop_count + 1;
        if stop_count > 30
            break;
        end
    else
        last_dist = bestTD;
        stop_count = 0;
    end
    
    
end
ObjV = calObj(chroms,customer_number,time_window1,time_window2,depot_time_window2,service_time,distance_matrix);
[minObjV,minInd]=min(ObjV);

disp('Optimal Result:')
bestChrom=chroms(minInd(1),:);
[bestVC,bestNV,bestTD,best_vionum,best_viocus]=decode(bestChrom, customer_number, time_window1, time_window2, depot_time_window2, service_time, distance_matrix);
disp(['Number of Robots: ', num2str(bestNV), ', Total Distance: ', num2str(bestTD), ', Number of violated Path: ', num2str(best_vionum), ', Number of Violated Customer: ', num2str(best_viocus)]);

flag = Judge(bestVC, time_window1, time_window2, depot_time_window2, service_time, distance_matrix);

% check the missing element in the final result
missing_element = judgeFullElement(bestVC, customer_number);

drawMap(bestVC, vertexs);

编码函数:

很抱歉，代码中并没有提供编码函数的具体实现，需要根据上下文推测编码方式。一般来说，遗传算法中的编码方式有二进制编码、实数编码、排列编码等多种方式，具体应该根据问题的特点选择合适的编码方式。

推测：

基于代码中 length_chrom = customer_number + robot_number - 1 的定义，以及 chroms(1,:) 传递给 decode 函数，可以推测代码使用的编码方式可能是排列编码，即用一个数组表示所有顾客和机器人的访问顺序。

需要进一步了解具体实现才能确定编码方式。