【VRP问题】基于遗传算法求解带时间窗多机场多飞机车辆路径问题Matlab代码

1 简介

某地区拥有M个机场，机场分布在不同的地理位置，每个机场均拥有K_m(m=1,2,…,M)架无人机。机场m(m=1,2,…,M)的第k(k=1,2,…,K_m)架飞机的载弹量为Q_mk，平均速度V_mk。在T₀时刻，在机场周边分布N个任务目标，且某些（个）任务对执行时间有较严格的要求，任务i的执行时间窗为[ET_i ,LT_i]，执行任务i需要的弹药数量为W_i(i=0,1,2,…,N)。其中ET_i表示任务目标i的最早执行时间，LT_i表示任务目标i最晚执行时间。p_mk1表示提前完成损失系数，即机场m的无人机k每提前单位时间到达的时间成本；m_k2p表示延期完成惩罚系数，即机场m的无人机k每推迟单位时间到达的时间成本。当该地区收到对多个目标进行攻击的任务时，可调用该地区所有的无人机，且目标的数量、位置及时间窗是确定的；军用机场的数量、位置，飞机的类型是已知的。

该问题不仅需要确定每个目标应该由哪个机场以及无人机去执行，同时还要确定执行任务的先后顺序，且在不违反各个任务时间窗限制，飞机航程限制、火力威胁、探测威胁等约束条件下满足代价最小化。具体符号如表所示：

根据上述描述，问题有如下的假设：

a) 机场的位置及无人机种类、无人机数量、无人机时间窗已知。任务目标位置、完成任务所需弹药量、任务执行时间窗已知；

b) 每架无人机从机场起飞出发执行任务最终又回到原机场；

c) 每个机场调度的无人机数量不能超过其所拥有的最大无人机数；

d) 每次任务中每架无人机能调度一次，且不得超过其最大载弹量；

e) 每个任务目标必须且只能由一架无人机完成；

f) 假设无人机飞行的高度有航管中心统一分配，不考虑飞机碰撞等特殊状况。

构建数学模型时，同时考虑车辆时间窗和客户软时间窗约束，根据违反客户要求时间的长短施以相应的惩罚，从而在目标函数增加相应的时间成本。因此，带时间窗和任务软时间窗的多车场车辆路径问题的数学模型为：

目标函数：

2 部分代码

clc;

clear all;

close all;

global popsize

global D

global dis_airport

global point_num

global plane_number

global all_point

global belong_to_airport

global demand

global time_window

global speed

global service_time

global capacity_plane

global danger_area_axis

global danger_radius

global threat_matrix

global penatly_coefficient_dis

global penatly_coefficient_load

global penatly_coefficient_early_time

global penatly_coefficient_late_time

global airport_scope

global plotting_sacle    %比例尺

global danger_num 

maxgan=400; %迭代次数

popsize=100;  %种群大小

%各种惩罚系数

penatly_coefficient_dis = 100;%距离惩罚系数

penatly_coefficient_load = 1000;%负载惩罚系数

penatly_coefficient_early_time =100;%最早时间惩罚系数

penatly_coefficient_late_time = 100;%最晚时间惩罚系数

airport_scope = 300; %机场辐射范围

plotting_sacle=2;    %比例尺等于1

%%读取数据

file='data.xlsx';

data_1=xlsread(file,'机场信息');

data_2=xlsread(file,'飞机信息');

data_3=xlsread(file,'目标点信息');

data_4=xlsread(file,'威胁区域');

airport_num = size(data_1,1);   %飞机场数量

coor_airport = data_1(1:airport_num,1:2);%飞机场经纬度

scope_airport = data_1(1:2,3)/1000;%飞机场辐射范围

plane_number = size(data_2,1);%飞机数量

capacity_plane = data_2(1:plane_number,2);%各个飞机容量

speed = data_2(1:plane_number,3)*3.6;%飞机速度 km/h

belong_to_airport = data_2(1:plane_number,4);%飞机所属机场

point_num = length(data_3);%点的数量

coor_point = data_3(1:point_num,1:2);%点的经纬度

demand = data_3(1:point_num,3);%点的需求

time_window = data_3(1:point_num,4:5);%时间窗

service_time = data_3(1:point_num,6);

danger_num = size(data_4,1);%危险区域数量

coor_danger_point=data_4(1:danger_num,1:2);%危险区域经纬度

danger_radius = data_4(1:danger_num,3:4)/1000;%危险区域半径 km

[n,m]=size(coor_point);

X=zeros(n+airport_num,m);   %横纵坐标矩阵

min_long=floor(min(coor_point(:,1)));

min_lat=floor(min(coor_point(:,2)));

for i =1:n

   X(i,1)=computeDistanceLatLon(coor_point(i,1),min_lat,min_long,min_lat)/1000;%单位是km

   X(i,2)=computeDistanceLatLon(min_long,coor_point(i,2),min_long,min_lat)/1000;

end

for i=1:airport_num

   X(n+i,1)=computeDistanceLatLon(coor_airport(i,1),min_lat,min_long,min_lat)/1000;%单位是米

   X(n+i,2)=computeDistanceLatLon(min_long,coor_airport(i,2),min_long,min_lat)/1000;

end

danger_area_axis=[danger_num,2];

for i = 1:danger_num

   danger_area_axis(i,1)=computeDistanceLatLon(coor_danger_point(i,1),min_lat,min_long,min_lat)/1000;%单位是km

   danger_area_axis(i,2)=computeDistanceLatLon(min_long,coor_danger_point(i,2),min_long,min_lat)/1000;

end

D=calculate_distance(X);  %距离矩阵

dis_airport = D(1:point_num,point_num+1:point_num+airport_num);

threat_matrix=cell(airport_num+point_num,airport_num+point_num);

for i = 1:airport_num+point_num

   for j = 1:airport_num+point_num

       if i == j 

           continue;

       else

           threat_matrix{i,j}=threat_value(i,j,X);

       end

   end

end

%%获得初始解

all_point = point_num+plane_number; %飞机数量和需求点数量的总和

pop=zeros(popsize,all_point);   

first_route=[point_num+1,1:point_num];

for k=2:plane_number

   r = randperm(length(first_route),1);

   first_route=[first_route(1:r),point_num+k,first_route(r+1:length(first_route))];

end

pop(1,:)=first_route;    

for k=2:popsize

   start_point = randperm(plane_number,1)+point_num;

   p = [1:start_point-1,start_point+1:all_point];

   pop(k,1)=start_point;

   ord = randperm(all_point-1);

   for j = 1:all_point-1

       pop(k,j+1)=p(ord(j));

   end

end

%%开始迭代

bestfitvalue=inf; %最优函数值

bestroute=[];%最优路线

distHistory = zeros(maxgan,1);%历史最优

for iter=1:maxgan  %迭代次数

   fprintf('第%d轮迭代\n',iter); 

  [fitvalue,efficiency]=func(pop);

  [mindist,ind]=min(fitvalue);

   if mindist<bestfitvalue && efficiency(ind)==1   %更新最优解

       bestfitvalue=mindist;

       bestroute=pop(ind,:);

   end

   distHistory(iter)=bestfitvalue;   

   %%交叉变异

   randomOrder = randperm(popsize);

   newPop=zeros(popsize,all_point);

   for p = 4:4:popsize

       tmpPop=zeros(4,all_point);

       rtes = pop(randomOrder(p-3:p),:);

       dists = fitvalue(randomOrder(p-3:p));

      [ignore,idx] = min(dists); 

       bestOf4Route = rtes(idx,:);

       routeInsertionPoints = sort(ceil(rand(1,2)*(all_point-1))+1);

       I = routeInsertionPoints(1);

       J = routeInsertionPoints(2);

       for k = 1:4 % Mutate the Best to get Three New Routes

           tmpPop(k,:) = bestOf4Route;

           switch k

               case 1 

                   tmpPop(k,[I J]) = tmpPop(k,[J I]);

               case 2 % Flip

                   tmpPop(k,I:J) = tmpPop(k,J:-1:I);

               case 3 % Swap

                   tmpPop(k,[I J]) = tmpPop(k,[J I]);

               case 4 % Slide

                   tmpPop(k,I:J) = tmpPop(k,[I+1:J I]);

               otherwise % Do Nothing

           end

       end

       newPop(p-3:p,:) = tmpPop;

   end

   parent_pop = [pop;newPop];

   parent_pop=unique(parent_pop,'rows');

  [par_val,whe]=func(parent_pop);

  [fit,index_new]=sort(par_val);

   for u=1:popsize

       pop(u,:)=parent_pop(index_new(u),:);

   end

end

bestfitvalue

bestroute

figure(1);

plot(1:maxgan,distHistory);

hold on 

aa=zeros(plane_number+1,1);

aa(plane_number+1)=all_point+1;

aa(1:plane_number) = find(bestroute>point_num);

route = cell(1,plane_number);

for i = 1:plane_number

   route{i}=bestroute(aa(i):(aa(i+1)-1));

end

color=['k','r','b','g','m','c','k','r','b','g','m','c','k','r','b','g','m','c'];

figure(2)

%各个点的散点图

hold off

unit_long = 1000/computeDistanceLatLon(min_long,min_lat,min_long+1,min_lat);

unit_lat = 1000/computeDistanceLatLon(min_long,min_lat,min_long,min_lat+1);

kk=[coor_point;coor_airport];

xlswrite('result.xlsx',bestfitvalue,'sheet1','A1:A1');

xlswrite('result.xlsx',bestroute,'sheet1','A2:BB2');

all=zeros(100,300);

index_re=1;

for i = 1:plane_number

   ll = route{i};

   for j = 1:length(ll)-1

       rou = adjustment(transform(ll(j)),ll(j+1),X);

       if size(rou,1)==2

           change = [kk(transform(ll(j)),:);kk(ll(j+1),:)];

       else

           %change=rou.*[unit_long,unit_lat]+[min_long,min_lat];

           change=[rou(:,1)*unit_long,rou(:,2)*unit_lat];

           change=[change(:,1)+min_long,change(:,2)+min_lat];

       end        

       all(index_re,1:2)=[ll(j),ll(j+1)];        

       all(index_re+1:index_re+2,1:size(change,1))=change';   

       index_re=index_re+3;

   end

   if j

       rou = adjustment(ll(j+1),transform(ll(1)),X);

       if size(rou,1)==2

           change=[kk(ll(j+1),:);kk(transform(ll(1)),:)];

       else

           %change=rou.*[unit_long,unit_lat]+[min_long,min_lat];

           change=[rou(:,1)*unit_long,rou(:,2)*unit_lat];

           change=[change(:,1)+min_long,change(:,2)+min_lat];

       end

       all(index_re,1:2)=[ll(j+1),ll(1)];        

       all(index_re+1:index_re+2,1:size(change,1))=change';   

       index_re=index_re+3;

   end    

end

xlswrite('result.xlsx',all,'sheet2')

3 仿真结果

4 参考文献

[1]张露. "基于改进遗传算法求解带时间窗车辆路径规划问题." 中国物流与采购 14(2020):4.

部分理论引用网络文献，若有侵权联系博主删除。