Matlab使用遗传算法（GA）求解应急物资配送路径（VRP）问题

1 应急物资配送路径问题描述与建模

本文研究要解决的应急物资配送路径优化问题有：单个调配中心，多个受灾点，已知调配中心和各受灾点的位置，以及各受灾点的需求量、卸货时间等。要求每个受灾点只能由一辆车进行配送，且完成配送任务后需重返配送中心。其目标是优化配送车辆的路径，使得在时间窗约束下总成本最低，应急物资配送路线示意图如图所示。

模型假设
(1)假设每个受灾点的应急物资需求量不超过运输车辆的最大载重量，且每个受灾点仅能被一辆车服务。
(2)假设运输车辆的初始位置为各调配中心，且完成配送任务后需重返调配中心。
(3)假设调配中心的应急生活物资数量足以满足各受灾点的需求。
(4)假设运输车辆为同一类车型，具有相同的载重量。

决策变量如下：

目标函数应急物资配送车辆的总运输费用如下：

其中c为应急物资配送车辆的运输耗油总费用，c2为每辆车的固定成本，ct1和ct2分别是装卸货物的时间。

GA求解

求解流程图如下：

其中主要确定的就是编码方式，初始种群的产生方式，适应值计算函数，选择策略，交叉策略和变异策略，接下来仔细说明。

整数编码

对于VRP问题，传统的编码方式主要是利用整数编码，加入受灾点有N个，运输车辆有K辆，我们可以使用1-N个整数来代表受灾点，K-1（首尾不加0）或者K+1（首尾加0）个0来表示运输车辆的路径。本文选用N+K-1的方式进行编码。染色体示例如下：

其中有6辆车和6个受灾点。路径分别为：1，0，3-4，2，6，5。两个0相邻表示该辆车没有被使用。

种群初始化

种群的初始化采用随机初始化的方式。代码如下：其中

function [outputs] = create_generation(size, len, city_num, demand)
    % size 为种群的数量   len为染色体的长度  city_num为受灾城市的数量  demand为受灾点的需求
    for i=1:size
        while 1
            zero_pop = zeros(1,len);
            rand_city = randperm(city_num);
            city_index = randperm(len, city_num);
            zero_pop(city_index) = rand_city;
            if check_constrain(zero_pop, demand)  % 检查是否满足对应的约束条件，这里可以改为自己的约束条件
                break;
            end
        end
        outputs(i,:) = zero_pop;
    end
end

主要是通过随机生成城市的插入位置，将对应的城市插入到全零的数组中即可，注意的是初始化的时候加入了一个约束检查，确保产生的解都是满足约束条件的，这里可以根据问题的不同换成自己对应的约束条件。

适应值计算

本次的需求是最小化成本也就是常规的运输费用等等。主要涉及的是一个染色体如何进行解码的问题。同时考虑到后面的交叉和变异会产生出不满足约束条件的个体，这里我们选用惩罚值进行限制，也就是将不满足约束条件的适应值一个惩罚，如将其适应值变为最低的哪个，下次选择就不会进行选择到了。

本次计算的流程：先将染色体两端的0去除，因为这个就相当于空车。接下来进行遍历。如果遇到两个连续的0也进行跳过，此外每次经过城市对应的计算相应的参数如运输距离、油费等，加入适应值最好即可完成计算。代码如下：

function [fitn, valide] = fitness(population, d_t, city, need)
% population就是种群  d_t是每一个城市的卸货时间  city是每一个城市的坐标  need是需求
    v = 60;  %运输车辆在未发生灾难前平均速度
    c1 = 7.64;  % 运输车辆的油价
    c2 = 300; % 运输车辆正常磨损维护费用
    r = 0.18;  % 运输车辆耗油量
    sc = 100;  % 时间成本
    %% 计算成本
    fitn = [];
    valide = [];
    
    for pop=population'
        if ~check_constrain(pop', need)  % 如果不满足约束 直接进行跳过
            valide(end+1) = false;
            fitn(end+1) = 0;  % 适应值为0
            continue;
        end
        valide(end+1) = true;
        C = 0;  % 运输耗油成本
        Ct1 = 2*sum(d_t)*sc;  % 装卸时间成本 固定
        Ct2 = 0;  % 运输时间成本
        C2 = 0;  % 车辆维修成本
        % 事先去除首尾的多余的0只保留一个0，因为没有意义
        first_nonzero = find(pop, 1, 'first');
        last_nonzero = find(pop, 1, 'last');
        pop = [0;pop(first_nonzero:last_nonzero);0];
        car_nums = 0;
        for i=1:length(pop)
            if pop(i)==0 
                if ((i == 1) || (pop(i-1)==0)) % 两个零连一起没有意义
                    continue;
                end
                car_nums = car_nums + 1;               
            end
            %  在这里计算相应的成本
            distance = pdist([city(pop(i-1)+1,:); city(pop(i)+1,:)], 'euclidean');
            C = C + distance*r*c1;
            Ct2 = Ct2 + distance*sc/v;
        end
        fitn(end+1) = 1/(C + Ct1 +Ct2 + car_nums*c2);
    end

选择操作

选择操作主要用于从父代中选择出用于交叉和变异的子代个体，这里我们采用的是轮盘赌注选择的方法，很经典具体不介绍了，具体看代码；

function SelCh = Select(chroms,FitnV,GGAP)
% chroms种群   FitnV适应值列表  GGAP选择父代的比例
    N = size(chroms,1); 
    NSel = max(floor(N * GGAP + .5), 2);   % 计算选择多少父代
    % Calculate the cumulative fitness values
    CumFitnV = cumsum(FitnV) / sum(FitnV);
   % Select chromosomes using roulette wheel selection
    SelCh = zeros(NSel, size(chroms,2));
    for i = 1:NSel
        r = rand;
        ix = find(CumFitnV >= r, 1, 'first');
        SelCh(i,:) = chroms(ix,:);
    end
end

交叉操作

由于问题的复杂性，在进行交叉的时候很容易产生不符合约束的子代，此时可以通过修复算子进行修复（将交叉后的子代按照一定的规则进行修改使其符合约束），本次为了简单没有单独编写复杂的修复算子，而是将修复算子直接定义为重新初始化一个新的个体。

此外交叉的具体方式采用顺序交叉的方式，先随机选择两个交叉点，交换父代的中间的这段，然后将剩余的按照另一个父代的出现顺序依次插入。

function [a_org,b_org]=Crossover(a_org,b_org)  % 只管交叉不用管是否符合约束。
    a_index = find(a_org~=0);
    b_index = find(b_org~=0);
    a = a_org(a_index);
    b = b_org(b_index);
    L=length(a);
    while 1
        % generate a random number in the range [1:L]
        % 1. choose 2 points
        r1 = randsrc(1,1,[1:L]);
        r2 = randsrc(1,1,[1:L]);
        if r1 ~= r2
            % start and end
            s = min([r1,r2]);
            e = max([r1,r2]);
            % choose 2 parts of chromosome
            a0 = [b(s:e),a];
            b0 = [a(s:e),b];
            for i = 1:length(a0)
                aindex = find(a0==a0(i));
                bindex = find(b0==b0(i));
                % delete the same index
                if length(aindex) > 1
                    a0(aindex(2))=[];
                end
                
                if length(bindex) > 1
                    b0(bindex(2)) = [];
                end
                if i == length(a)
                    break
                end
            end
            a = [a0(e-s+1:e-s+s), a0(1:e-s), a0(e-s+s+1:end)];
            b = [b0(e-s+1:e-s+s), b0(1:e-s), b0(e-s+s+1:end)];
            break
        end
    end
    a_org(a_index) = a;
    b_org(b_index) = b;
end

变异操作

变异采用简单的两点交换，也就是随机选择两个点位，交换俩个基因的位置即可。

function SelCh = Mutate(SelCh,Pm)  % 两点交叉
    [NSel,L] = size(SelCh);
    for i = 1:NSel
        if Pm >= rand
            % Random permutation of integers
            R = randperm(L);
            % choose 2 random number in ith chromosome
            SelCh(i,R(1:2)) = SelCh(i,R(2:-1:1));
        end
    end
end

至此大概的流程基本都已经介绍完毕，如果想要做对应的更改只需要简单的修改对应的函数步骤约束检查和适应度计算即可。

实验

采用20个受灾点和20辆车进行实验。参数设定如下

结果如下：

适应度不断下降。最终得到的路线图如下：

需要代码的同学可以咸鱼搜breeze002联系，需要帮助也可以找我。