求解带时间窗的车辆路径问题（matlab实现）

直接上代码

在一个目录下把下面的文件都建好直接就能跑起来了
后面我会把用到的数学模型和算法详细介绍、根据不同实际问题给出变种问题的求解思路，先去给老板写java了555555555555

main.m

%%
Iteration = 1000;   %迭代轮数
number = 51;      %订单数
route_number = 20;  %车辆数
capacity = 30;  %车辆容量
tensure = 50;  %禁忌时长
Tabu = zeros(number+1,route_number);  %禁忌表
Tabu_new = zeros(number+1,1);  %禁忌表，使用新路径
distance = zeros(number+1,number+1);  %距离矩阵
num = zeros(route_number,1);  %路径中的节点数
bag = zeros(route_number,1);  %路径中的需求量
val_route = zeros(route_number,1); %单条路径的路程成本
val_time = zeros(route_number,1);%单条路径的超时成本（依赖时间窗）
%%
filename = 'R101.txt';  
[s.id,s.x,s.y,s.demand,s.early,s.last,s.ser] = textread(filename,'%n%n%n%n%n%n%n');  %读取订单数据
for i = 1:number
    for j = 1:number
        distance(i,j) = ((s.x(i)-s.x(j))^2+(s.y(i)-s.y(j))^2)^0.5;  %初始化距离矩阵
    end
end
%%
route=cell(route_number,1);%初始化路径，c1数量
r = zeros(number);%节点归属数组
for i=1:route_number
    route{i,1}=[route{i,1},1,1];%插入起始点与终点
    num(i,1)=2;%更新路径中节点个数为2
end
temp=randperm(number);%将节点打乱，存放在temp
for k = 1:number
    c = temp(k);
    for i = 1:route_number %寻找可行的路径插入节点c
        if(c~=1 && bag(i)+s.demand(c) < capacity)
            for j= 2:num(i,1)
                if s.last(c) < s.early(route{i,1}(j))
                    route{i,1}=[route{i,1}(1:j-1),c,route{i,1}(j:num(i,1))];%节点插入路径
                    bag(i)=bag(i)+s.demand(c);%i路径增加对应的需求
                    r(c)=i;%确定节点归属哪条路径
                    num(i,1)=num(i,1)+1;
                end
                if r(c)~=0 %如果节点已有归属，退出for循环
                    break;
                end
            end
        end
        if r(c)~=0 %如果节点已有归属，退出for循环
            break;
        end
    end
end
for i=1:route_number
    if num(i,1)>2 %路径中存在订单
        size=num(i,1);
        t=0;
        for j=1:size-1
            val_route(i,1)=val_route(i,1)+distance(route{i,1}(j),route{i,1}(j+1));%路径i的成本：第j与第j+1个节点的距离之和，j<size
            t=t+distance(route{i,1}(j),route{i,1}(j+1))+s.ser(route{i,1}(j+1));
            if t>s.last(route{i,1}(j+1))
                val_time(i,1)=val_time(i,1)+t-s.last(route{i,1}(j+1));
            else
                t=s.last(route{i,1}(j+1));
            end
        end
    end
end
%%
[val_log,val_test,val_gol,route,num,val_time,val_route]=shift(r,s,route,num,val_route,val_time,Iteration,number,route_number,distance,Tabu,tensure,bag,capacity,Tabu_new);

shift.m

function [val_log,val_test,val_gol,route,num,val_time,val_route]=shift(r,s,route,num,val_route,val_time,Iteration,number,route_number,distance,Tabu,tensure,bag,capacity,Tabu_new) %r节点所属路径、s为节点信息、route为路径、num为每条路径的节点个数、val_route为每条路径的成本
    val_log=[]; %用来记录每一轮找到的最优候选解
    val_test=[];%用来验证运算正确
    val_gol=99999; %用来记录找到的全局最优解
    for ite=1:Iteration
        val_cur=99999;%用来保存当前轮中的最有候选解
        for i=2:number
            cur_r=r(i,1);%得到i节点所在路径cur_r
            former=find(route{cur_r,1}==i);%记录取出的节点先前的位置
            route1=val_route(cur_r,1);%记录被取节点路径的原先成本
            time1=val_time(cur_r,1);%记录被取节点路径的原先时间成本
            val_route(cur_r,1)=val_route(cur_r,1)+distance(route{cur_r,1}(former-1),route{cur_r,1}(former+1))-distance(route{cur_r,1}(former-1),i)-distance(route{cur_r,1}(former+1),i);%计算路径被取出节点后的成本
            route{cur_r,1}(former)=[];%删除i节点
            val_time(cur_r,1)=getTime(route{cur_r,1},s,distance);
            num(cur_r,1)=num(cur_r,1)-1;%更新cur_r号路径所剩节点数
            for j=1:route_number
                if j~=cur_r && (num(j,1)>2 || (num(j,1)==2 && Tabu_new(i,1) < ite)) && Tabu(i,j) < ite%排除禁止新路径和节点原属路径
                    route2=val_route(j,1);%保存j路径的成本，后面会实施插入删除节点操作，方便复原
                    time2=val_time(j,1);%保存cur_r和j路径的时间成本
                    for k=2:num(j,1) %尝试将i节点插入j路径的k位置
                        val_route(j,1)=val_route(j,1)+distance(route{j,1}(k-1),i)+distance(route{j,1}(k),i)-distance(route{j,1}(k-1),route{j,1}(k));
                        route{j,1}=[route{j,1}(1:k-1),i,route{j,1}(k:num(j,1))];%插入i节点
                        num(j,1)=num(j,1)+1;%更新j路径节点数
                        val_time(j,1)=getTime(route{j,1},s,distance);
                        now=sum(val_route)+sum(val_time);%计算插入后的总成本
                        flag=overflow(bag,capacity,j,s.demand(i));%判断是否满足背包约束，满足返回0，否则返回1
                        if now<=val_cur && flag==0
                            val_cur=now;besti=i;bestr=j;bestk=k;%记录最优候选解信息，被取节点i，插入路径j，插入位置k
                        end
                        route{j,1}(k)=[];%删除i节点
                        num(j,1)=num(j,1)-1;%更新j号路径所剩节点数
                        val_route(j,1)=route2;%路径成本复原
                        val_time(j,1)=time2;%路径时间成本复原
                        if flag==1%如果超过背包约束，不再尝试后续位置插入
                            break;
                        end
                    end
                end
            end
            val_route(cur_r,1)=route1;%路径成本复原
            val_time(cur_r,1)=time1;%路径时间成本复原
            route{cur_r,1}=[route{cur_r,1}(1:former-1),i,route{cur_r,1}(former:num(cur_r,1))];%i返回原先位置
            num(cur_r,1)=num(cur_r,1)+1;%更新cur_r路径节点数
        end
        %下面开始用besti，bestr，bestk将路径变成最优候选解
        cur_r=r(besti,1);
        former=find(route{cur_r,1}==besti);
        i=besti;
        val_route(cur_r,1)=val_route(cur_r,1)+distance(route{cur_r,1}(former-1),route{cur_r,1}(former+1))-distance(route{cur_r,1}(former-1),i)-distance(route{cur_r,1}(former+1),i);
        route{cur_r,1}(former)=[];%删除i节点
        num(cur_r,1)=num(cur_r,1)-1;
        bag(cur_r,1)=bag(cur_r,1)-s.demand(i);
        j=bestr;k=bestk;
        val_route(j,1)=val_route(j,1)+distance(route{j,1}(k-1),i)+distance(route{j,1}(k),i)-distance(route{j,1}(k-1),route{j,1}(k));
        route{j,1}=[route{j,1}(1:k-1),i,route{j,1}(k:num(j,1))];%插入i节点
        num(j,1)=num(j,1)+1;
        bag(j,1)=bag(j,1)+s.demand(i);
        val_time(cur_r,1)=getTime(route{cur_r,1},s,distance);
        val_time(j,1)=getTime(route{j,1},s,distance);
        r(i,1)=j;%更新i节点的归属为j路径
        %更新禁忌表
        Tabu(i,j)=Tabu(i,j)+tensure;
        Tabu(j,i)=Tabu(i,j);
        if num(j,1)==3%使用新路径的禁忌表更新
            Tabu_new(i,1)=Tabu_new(i,1)+tensure;
        end
        %下面开始结算最优解
        val_log=[val_log,val_cur];  %记录每一轮的最优候选解，最后log与test应该相同，如果不同说明中间计算的代码有问题！
        val_test=[val_test,sum(val_route)+sum(val_time)];
        if val_cur < val_gol
            val_gol=val_cur%输出找到的最优解
        end     
    end
end

overflow.m

function f = overflow(bag,capacity,j,i)
    temp=bag(j,1)+i;
    if temp <= capacity
        f=0;
    else
        f=1;
    end
end

getTime.m

function f=getTime(route,s,distance)
    t=0;f=0;
    size=numel(route);
    for i=1:size-1
        t=t+distance(route(i),route(i+1))+s.ser(route(i+1));
        if t>s.last(route(i+1))
            f=f+t-s.last(route(i+1));
        else
            t=s.last(route(i+1));
        end
    end
end

R101.txt

序号、xy坐标、需求量、时间窗、服务时间

1  15  15  0  300  300  0
2  8  3  3  205  217  1
3  23  24  3  59  79  1
4  23  14  6  103  123  1
5  23  5  5  35  51  1
6  26  28  2  97  114  1
7  15  5  3  117  134  1
8  5  26  1  199  216  1
9  5  18  3  127  139  1
10  5  4  9  36  46  1
11  29  15  3  220  239  1
12  25  23  8  42  53  1
13  25  26  4  113  127  1
14  16  17  4  71  86  1
15  9  19  8  112  130  1
16  22  19  3  201  216  1
17  14  17  4  37  53  1
18  12  19  7  149  163  1
19  8  26  7  29  39  1
20  7  12  4  124  142  1
21  18  12  9  144  161  1
22  3  5  8  217  232  1
23  10  22  2  141  160  1
24  15  5  6  76  90  1
25  7  12  8  53  68  1
26  28  0  1  104  114  1
27  13  10  0  179  196  0
28  5  2  0  27  42  0
29  27  24  0  90  104  0
30  28  3  0  14  27  0
31  14  26  0  62  74  0
32  28  29  0  92  102  0
33  13  1  0  188  198  0
34  21  29  0  98  112  0
35  29  26  0  14  28  0
36  23  26  0  199  218  0
37  14  27  0  17  35  0
38  29  26  0  87  103  0
39  27  2  0  48  63  0
40  9  16  0  83  102  0
41  25  26  0  190  207  0
42  22  6  0  17  29  0
43  22  19  0  128  139  0
44  8  22  0  16  27  0
45  0  24  0  109  124  0
46  28  9  0  120  140  0
47  13  7  0  196  209  0
48  6  16  0  110  120  0
49  27  16  0  42  62  0
50  16  23  0  33  50  0
51  10  19  0  83  95  0