遗传算法结解决31个客户点的CVRP问题（存在车辆容量和行驶距离限制-MATLAB，引用他人的思想）

引用链接

针对31个客户的配送问题
引用链接中的代码采用了基于距离加入客户节点的编码方式，之后的迭代采用的是正常的遗传算法的锦标赛选择、顺序交叉、单点变异操作。我又严格按照引用链接中的引用链接（俄罗斯套娃，哈哈哈）中想法，选择自身逆序、当前最优和全局最优作为可能的一个父代进行交叉，且没有变异操作。得到的结果如下：①得到最优解的概率较低②程序运行时间较长，大于1分钟③针对7个和12个客户点的配送问题，用下面的代码就可以实现极大概率找到最优解
下面是我按照上面引用链接中的引用链接（俄罗斯套娃，哈哈哈）中想法进行的matlab编程，即将他的程序从python翻译到matlab。而上面引用链接中的代码是只采用了距离的编码。

代码

①主函数main

% 问题：找到精确解的概率较低
%编码：自然数编码，先随机选择一个客户点，再根据距离依次选择下一个客户点，直到所有客户都遍历到
%解码：根据客户需求量和车辆行驶距离进行分配车辆
%选择：单纯从自身的逆序、当前最优、全局最优选择另一个父代
%交叉：顺序交叉
%变异：无变异操作
tic;
clear;
clc;
close all;
%% 用importdata读取文件并计算各个点之间的距离
CVRP=importdata('c101.txt');
[row,~] = size(CVRP);
dist = zeros(row,row);
for i = 1:row-1
    for j = i+1:row
        dist(i,j) = sqrt(sum((CVRP(i,2:3)-CVRP(j,2:3)).^2));
    end
end
for i = 1:row
    dist(i,i) = inf;
end
dist = dist + dist';
row = row - 1;%客户点的数量

%% 参数设置
[N,D,volume,Dist,C0,C1,Iter,w,c1,c2] = constant(row);
GA = zeros(N,D+1);%GA数组位置申请
Pbest = zeros(N+1,D+1);%当前最优和全局最优
Path = inf(N+1,2*D+2);%定义路径的数组，前2*D存放具体路径，2*D+1、2*D+2分别为车辆书及路径长度

%% 初始化
[GA,Pbest,Path] = initialization(CVRP,dist,N,D,volume,Dist,C0,C1,GA,Pbest,Path);

%% 迭代
[GA,Pbest,Path,Gbest] = iteration(N,D,volume,Dist,C0,C1,Iter,w,c1,c2,CVRP,dist,GA,Pbest,Path);

%% 命令行显示
[value,row] = min(Path(:,2*D+2));
[~,col] = find(Path(row,1:2*D)==1);
[~,k] = size(col);
for i = 1:k-1
    disp(['车辆',num2str(i),':',num2str(Path(row,col(i):col(i+1)))]);
end
disp(['最短总路程：',num2str(Path(row,2*D+2))]);        %输入多个参数进行显示，除了用num2str显示数值，还需要将所有输入加上方括号！！

plot(Gbest);
draw(CVRP);
hold on;
grid on ;
for i = 1:k-1
    c = Path(row,col(i):col(i+1));
    plot(CVRP(c,2),CVRP(c,3));
end
toc;

②参数函数constant

function [N,D,volume,Dist,C0,C1,Iter,w,c1,c2] = constant(row)

N = 50;
D = row;
if D ==7                %车辆的容量
    volume = 1;
elseif D ==12
    volume = 6;
elseif D == 31
    volume = 120;
end
Dist = 250;             %车辆最远行驶距离
C0 = 30;                %车辆启动成本
C1 = 1;                 %车辆的行驶成本
Iter = 1000;             %迭代数

w = 0.2;
c1 = 0.4;
c2 = 0.4;%采用所谓的粒子群思想选择交叉的父代，以w/(w+c1+c2)、w+c1/(w+c1+c2)、w+c1+c2/(w+c1+c2)分别选择粒子自身逆序、自身最佳、全局最佳作为父代之一

③初始化函数initialization

function [GA,Pbest,Path] = initialization(CVRP,dist,N,D,volume,Dist,C0,C1,GA,Pbest,Path)

for i = 1:N
    dist_matrix = dist;
    dist_matrix(:,1) = inf;
    k = 1;%路径中节点的个数
    GA(i,k) = 1;
    while GA(i,k)==1
        GA(i,k) = randperm(D+1,1);%随机选择一个起始客户点
    end
    dist_matrix(:,GA(i,k)) = inf;
    while k<=30
        k = k + 1;
        [~,col] = min(dist_matrix(GA(i,k-1),:));
        GA(i,k) = col;
        dist_matrix(:,GA(i,k)) = inf;
    end
    Path(i,:) = pathdecode(CVRP,dist,D,volume,Dist,GA(i,1:D));
    GA(i,D+1) = fitness(Path(i,:),D,dist,C0,C1);%计算距离
    Path(i,2*D+2) = GA(i,D+1);
end
Pbest = GA;
[~,row] = min(GA(:,D+1));
Pbest(N+1,:) = GA(row,:);
Path(N+1,:) = Path(row,:);

④解码函数pathdecode

function path = pathdecode(CVRP,dist,D,volume,Dist,x)
%解码
path = inf(1,2*D+2);%存放解码后的路径
k = 1;%路径中已有的节点个数
path(k) = 1;%第一个节点是配送中心
heavy = 0;%初始的车辆载重为0
numcar = 0;%车辆数初始为0
d = 0;%初始车辆的行驶距离
for j = 1:D
    k = k + 1;%已有的节点个数加1
    if (heavy+CVRP(x(j),4)<volume)&&(d + dist(path(k-1),x(j)) + dist(x(j),1)<Dist)
        path(k) = x(j);
        heavy = heavy + CVRP(x(j),4);%累计车辆的载重
        d = d + dist(path(k-1),x(j));%累计车辆的行驶距离
    else
        heavy = CVRP(x(j),4);%重置车辆的载货量
        path(k) = 1;%添加一个配送中心节点序号
        k = k + 1;
        path(k) = x(j);%下一个加入路径的节点
        d = dist(path(k-1),path(k));%重置车辆的行驶距离
        numcar = numcar + 1;%车辆数加1
    end
end
k = k + 1;
numcar = numcar + 1;
path(k) = 1;%车辆回到配送中心
path(2*D+1) = numcar;%记录车辆数

⑤适应度函数fitness

function y = fitness(path,D,dist,C0,C1)
%适应度值为车辆启动成本和车辆行驶成本的线性和
y = 0;%计算路径长度
pathlong = path(2*D+1)+1+D;
for i = 1:pathlong-1
    y = y + dist(path(i),path(i+1));
    %位置 1 处的索引无效。数组索引必须为正整数或逻辑值。原因是函数调用的参数输入问题或者确实是本行有问题!!!!
end
y = C1*y + C0*path(2*D+1);

⑥迭代函数iteration

function [GA,Pbest,Path,Gbest] = iteration(N,D,volume,Dist,C0,C1,Iter,w,c1,c2,CVRP,dist,GA,Pbest,Path)
Gbest = zeros(1,Iter);%存放每次迭代的最佳距离值
for iter = 1:Iter
    for i = 1:N
        GA(i,1:D) = crossover(N,D,w,c1,c2,GA(i,1:D),Pbest(i,1:D),Pbest(N+1,1:D));%交叉操作
        
        Path(i,:) = pathdecode(CVRP,dist,D,volume,Dist,GA(i,1:D));
        GA(i,D+1) = fitness(Path(i,:),D,dist,C0,C1);
        Path(i,2*D+2) = GA(i,D+1);
        if GA(i,D+1)<Pbest(i,D+1)
            Pbest(i,:) = GA(i,:);
        end
    end
    [value,row] = min(Pbest(:,D+1));
    if value<Pbest(N+1,D+1)
        Pbest(N+1,:) = Pbest(row,:);
    end
    Path(N+1,:) = pathdecode(CVRP,dist,D,volume,Dist,Pbest(N+1,1:D));
    Path(N+1,2*D+2) = fitness(Path(N+1,:),D,dist,C0,C1);

    Gbest(1,iter) = min(GA(:,D+1));
end

⑦交叉函数crossover

function x = crossover(N,D,w,c1,c2,x,pbest,gbest)
%顺序交叉
parent1 = x;
r = rand;
if r<w
    parent2 = fliplr(x);
elseif r<w+c1
    parent2 = pbest;
elseif r<w+c1+c2
    parent2 = gbest;
end
for i = 1:N
    point = randperm(D,2);
    point1 = min(point);
    point2 = max(point);%指定大小交叉点
    while max(point)-min(point)>29||max(point)==min(point)
        %交叉点之间的距离≥29则交叉不会产生变化
        point = randperm(D,2);%随机选择1-D范围内的两个不同的整数
        point1 = min(point);
        point2 = max(point);%指定大小交叉点
    end
    x(1,point1:point2) = parent1(point1:point2);%子代交叉点之间是父代1交叉点之间的序号
    parent2_1 = [parent2(point2+1:D) parent2(1:point2)];%父代2从交叉点2后一位开始循环一轮到该交叉点
    for cross = point1:point2
        parent2_1(parent2_1==x(1,cross)) = [];%将parent2_1中与子代中重复的序号去除
    end
    if point1 == 1
        x(1,point2+1:D) = parent2_1;%如果交叉点1是第一位，那直接赋值
    elseif point2 ==D
        x(1,1:point1-1) = parent2_1;%如果交叉点2是最后一位，也可以直接赋值
    else
        x(1,1:point1-1) = parent2_1(1:point1-1);%交叉点前的一段序号赋值
        x(1,point2+1:D) = parent2_1(point1:end);%交叉点后的一段序号赋值
    end
end

⑧绘图函数draw

function draw(CVRP)

figure(5);
% axis([0 110 0 110]);
title('配送路线');
hold on ;
grid on ;
plot(CVRP(:,2),CVRP(:,3),'r.','linewidth',10);%不用循环，
end

结果