I. Introduction to VRP

1 Basic principles of VRPVehicle Routing Problem (VRP) is one of the most important research problems in operations research. VRP is concerned with the path planning of a supplier and K points of sale, which can be summarized as: for a series of delivery points and receiving points, the organization calls certain vehicles, arranges appropriate driving routes, and makes the vehicles pass through them in an orderly manner, under specified constraints (e.g. The demand and quantity of goods, delivery and delivery time, vehicle capacity limit, mileage limit, driving time limit, etc.), and strive to achieve certain goals (such as the shortest total mileage of empty vehicles, the lowest total transportation cost, vehicles arrive at a certain time, the use of the smallest number of vehicles, etc.). The legend of VRP is as follows: 2 Problem Attributes and FaQsThe characteristics of vehicle routing problem are relatively complex and generally include four attributes: (1) Address characteristics include the number of parking lots, demand types and operation requirements. (2) Vehicle characteristics include: vehicle number, carrying weight constraints, carrying variety constraints, running route constraints, and working time constraints. (3) Other characteristics of the problem. (4) The objective function may be to minimize the total cost, or minimize the maximum operating cost, or maximize just-in-time operations.

3 Common problems are as follows:(1) Traveling salesman Problem (2) Vehicle Routing problem with Capacity Constraint (CVRP) The model was difficult to extend to other scenarios of VRP and the execution paths of specific vehicles were not known, so the model continued to be improved. (3) Vehicle routing problem with Time Windows (VRP with Time Windows) is a vehicle routing problem with Time Windows (VRP with Time Windows, VRPTW), considering that demand points have requirements on vehicle arrival Time. The Vehicle Routing problem with Time Windows (VRPTW) is a time window constraint for the customer to be visited on the VRP. In the VRPTW problem, in addition to the driving cost, the cost function also includes the waiting time caused by arriving at a customer early and the service time required by the customer. In VRPTW, vehicles must not only meet the constraints of VRP problems, but also meet the Time Window constraints of demand points. The Time Window constraints of demand points can be divided into two types, one is Hard Time Window, which requires that vehicles must arrive within the Time Window. Early arrival must wait, while late arrival is rejected. The other is the Soft Time Window. It is not necessary to arrive in the Time Window, but to arrive outside the Time Window must be punished. The biggest difference between the Soft Time Window and the hard Time Window is to replace waiting and rejection with punishment. Model 2(refer to 2017 A generalized Formulation for Vehicle Routing problems) : This model was formulated with two-dimensional decision variables (4) Collection and distribution problem (5) Reference for multi-yard vehicle routing problem (2005 Lim, Genetic Algorithm for multi-yard vehicle routing problem _ Zou Tong, 1996 Renaud)Since vehicles are homogeneous, the modeling here does not include vehicle dimensions in the variables. (6) priority constraint vehicle route problem (7) compatibility constraint vehicle route problem (8) random demand vehicle route problem

4 solutions (1) mathematical analysis method (2) human-computer interaction method (3) group and then arrange route method (4) group and then arrange route method (5) save or insert method (6) improve or exchange method (7) mathematical programming approximation method (8) heuristic algorithm

5 Comparison between VRP and VRPTW

2. Ant colony algorithm introduction

1 Origin and development history of Ant Colony Algorithm (ACA) Marco Dorigo et al. found that ant colonies can quickly find targets by secreting biohormones called pheromones to communicate foraging information when searching for food. Therefore, in his doctoral dissertation in 1991, he first systematically proposed a new intelligent optimization algorithm “Ant System” (AS for short) based on Ant population. Later, the proposer and many researchers made various improvements to the algorithm and applied it to a wider range of fields. Figure coloring problem, secondary assignment problem, workpiece sequencing problem, vehicle routing problem, job shop scheduling problem, network routing problem, large-scale integrated circuit design and so on. In recent years, M.Dorigo et al. further developed Ant algorithm into a general Optimization technology “Ant Colony Optimization (ACO)”, and called all the algorithms conforming to ACO framework “ACO algorithm”.

Specifically, individual ants start looking for food without first telling them where it is. When an ant finds food, it releases a volatile pheromone (called a pheromone, it evaporates over time and its concentration indicates how far the path is) into the environment that other ants sense as a guide. Generally, when there are pheromones on multiple paths, the ant will preferentially choose the path with high pheromone concentration, so that the pheromone concentration of the path with high pheromone concentration will be higher, forming a positive feedback. Some ants do not repeat the same route as others. They take a different route. If the alternative path is shorter than the original one, gradually more ants are attracted to the shorter path. Finally, after some time of running, there may be a shortest path repeated by most ants. In the end, the path with the highest pheromone concentration is the optimal one selected by the ant. Compared with other algorithms, ant colony algorithm is a relatively young algorithm, characteristics, such as distributed computing center control and asynchronous indirect communication between individuals, and is easy to be combined with other optimization algorithms, through many healthyenterprise continuously explore, to today has developed a variety of improved ant colony algorithm, but the principle of ant colony algorithm is still the main.

2 solving principle of ant colony algorithm Based on the above description of ant colony foraging behavior, the algorithm is mainly to simulate the foraging behavior of the following aspects: (1) the simulation diagram contains two kinds of pheromones in the scene, said a home, a place, said food and both pheromone volatilization in at a certain rate. (2) Each ant can perceive information only in a small part of the area around it. Ants searching for food, if within the scope of the perception, can directly in the past, if is beyond the scope of awareness, will be toward more than pheromones, ants can have a small probability pheromone many places don’t go, and instead, it is very important to the small probability event, represents a way of innovation, is very important to find a better solution. (3) Ants return to the nest using the same rules as when they find food. (4) When ants move, they will first follow the guidance of pheromone. If there is no guidance of pheromone, they will follow the direction of their moving inertia, but there is a certain probability of changing direction. Ants can also remember the path they have already walked, so as to avoid repeating the same place. (5) The ants leave the most pheromones when they find food, and then the farther away they are from the food, the less pheromone they leave. The rules for finding nest pheromones are the same as for food. Ant colony algorithm has the following characteristics: positive feedback algorithm, concurrency algorithm, strong robustness, probabilistic global search, does not rely on strict mathematical properties, long search time, easy to stop phenomenon. Ant transfer probability formula:In the formula, is the probability of ant K transferring from city I to city J; α and β were the relative importance of pheromones and heuristic factors, respectively. Is the pheromone quantity on edge (I, j); Is the heuristic factor; The next step for Ant K allows the selection of cities. The above formula is the pheromone update formula in the ant system, and is the pheromone quantity on the edge (I,j). ρ is the evaporation coefficient of pheromone, 0<ρ<1; Is the pheromone quantity left by the KTH ant on the edge (I,j) in this iteration; Q is a normal coefficient; Is the path length of the k ant during this tour. In ant system, pheromone update formula is:3. Solving steps of ant colony algorithm: (1) Initialization parameters At the beginning of calculation, it is necessary to initialize related parameters, such as ant colony size (ant number) m, pheromone importance factor α, heuristic function importance factor β, pheromone will emit money ρ, total pheromone release Q, maximum iteration times iter_max, initial value of iteration times iter=1. (2) construct solution space and randomly place each ant at different starting points. For each ant k (k=1,2,3… M), calculate the next city to be visited according to (2-1) until all ants have visited all cities. (3) update information su to calculate the path length of each ant Lk(k=1,2… , m), record the optimal solution (shortest path) in the current iteration number. At the same time, pheromone concentration on the connection path of each city was updated according to Equations (2-2) and (2-3). (4) Determine whether to terminate if iter<iter_max, set iter=iter+1 to clear the record table of ant paths and return to Step 2; Otherwise, the calculation is terminated and the optimal solution is output. (5) Determine whether to terminate if iter<iter_max, set iter=iter+1 to clear the record table of ant paths and return to Step 2; Otherwise, the calculation is terminated and the optimal solution is output. 3. Determine whether to terminate. If iter<iter_max, set iter=iter+1 to clear the record table of ant paths and return to Step 2. Otherwise, the calculation is terminated and the optimal solution is output.

Three, some source code

clear; clc; close all;
tic
%% input 
c101 = importdata('c101.txt');
% c101 = importdata('my_test_data.xlsx');
% depot_time_window1 = c101(1.5); % time window of depot
% depot_time_window2 = c101(1.6);

depot_time_window1 = TimeTrans(c101(1.5)); % time window of depot
depot_time_window2 = TimeTrans(c101(1.6));
vertexs = c101(:,2:3); 
customer = vertexs(2:end,:); % customer locations
customer_number = size(customer,1);
% vehicle_number = 25;
% time_window1 = c101(2:end,5);
% time_window2 = c101(2:end,6);

time_window1 = TimeTrans(c101(2:end,5));
time_window2 = TimeTrans(c101(2:end,6));

width = time_window2-time_window1; % width of time window
service_time = c101(2:end,7); 
h = pdist(vertexs);
dist = squareform(h); % distance matrix
%% initialize the parameters
ant_number = floor(customer_number * 1.5);                                                  % number of ants
alpha = 4;                                                        % parameter for pheromone
beta = 5;                                                         % paremeter for heuristic information
gamma = 2;                                                        % parameter for waiting time
delta = 3;                                                        % parameter for width of time window
r0 = 0.5;                                                         % a constant to control the movement of ants
rho = 0.85;                                                       % pheromone evaporation rate
Q = 5;                                                            % a constant to influence the update of pheromene
Eta = 1./dist;                                                    % heuristic function
iter = 1;                                                         % initial iteration number
iter_max = 200;                                                    % maximum iteration number

Tau = ones(customer_number+1,customer_number+1);                  % a matrix to store pheromone
Table = zeros(ant_number,customer_number);                        % a matrix to save the route
Route_best = zeros(iter_max,customer_number);                     % the best route
Cost_best = zeros(iter_max,1);                                    % the cost of best route



iter_time = [];
last_dist = 0;
stop_count = 0;



%% find the best route 
while iter <= iter_max
    %tic;
    % ConstructAntSolutions
    for i = 1:ant_number
        for j = 1:customer_number r = rand; np = NextPoint(i,Table,Tau,Eta,alpha,beta,gamma,delta,r,r0,time_window1,time_window2,width,service_time,depot_time_window2,di st); Table(i,j) = np; end end %% calculate the costfor each ant
    cost = zeros(ant_number,1);
    NV = zeros(ant_number,1);
    TD = zeros(ant_number,1);
    for i=1:ant_number
        VC = decode(Table(i,:),time_window1,time_window2,depot_time_window2,service_time,dist);
        [cost(i,1),NV(i,1),TD(i,1)] = CostFun(VC,dist);
    end
    %% find the minimal cost and the best route
    if iter == 1
        [min_Cost,min_index] = min(cost);
        Cost_best(iter) = min_Cost;
        Route_best(iter,:) = Table(min_index,:);
    else
        % compare the min_cost in this iteration with the last iter
        [min_Cost,min_index] = min(cost);
        Cost_best(iter) = min(Cost_best(iter - 1),min_Cost); 
        if Cost_best(iter) == min_Cost
            Route_best(iter,:) = Table(min_index,:);
        else
            Route_best(iter,:) = Route_best((iter- 1), :); end end %% update the pheromene bestR = Route_best(iter,:); % find out the best route [bestVC,bestNV,bestTD] = decode(bestR,time_window1,time_window2,depot_time_window2,service_time,dist); Tau = updateTau(Tau,bestR,rho,Q,time_window1,time_window2,depot_time_window2,service_time,dist); % %print 
    disp(['Iterration: ',num2str(iter)])
    disp(['Number of Robots: ',num2str(bestNV),', Total Distance: ',num2str(bestTD)]);
    fprintf('\n')
    %
    iter = iter+1;
    Table = zeros(ant_number,customer_number);
    
    %iter_time(iter) = 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
%% draw
bestRoute=Route_best(iter- 1, :); [bestVC,NV,TD]=decode(bestRoute,time_window1,time_window2,depot_time_window2,service_time,dist); draw_Best(bestVC,vertexs); figure(2)
plot(1:iter_max,Cost_best,'b')
xlabel('Iteration')
ylabel('Cost')
title('Change of Cost')
%% check the constraints, 1 == no violation
flag = Check(bestVC,time_window1,time_window2,depot_time_window2,service_time,dist)

toc

function draw_Best(VC,vertexs)
customer=vertexs(2:end,:);  
NV=size(VC,1); 
figure
hold on;box on
title('Map of Best Solution')
hold on;
C=linspecer(NV);
for i=1:NV
    part_seq=VC{i};
    len=length(part_seq);
    for j=0:len
        if j==0
            fprintf('%s'.'Route',num2str(i),'：');
            fprintf('%d->'.0);
            c1=customer(part_seq(1), :); plot([vertexs(1.1),c1(1)],[vertexs(1.2),c1(2)].The '-'.'color',C(i,:),'linewidth'.1);
        elseif j==len
            fprintf('%d->',part_seq(j));
            fprintf('%d'.0);
            fprintf('\n');
            c_len=customer(part_seq(len),:);
            plot([c_len(1),vertexs(1.1)],[c_len(2),vertexs(1.2)].The '-'.'color',C(i,:),'linewidth'.1);
        else
            fprintf('%d->',part_seq(j));
            c_pre=customer(part_seq(j),:);
            c_lastone=customer(part_seq(j+1), :); plot([c_pre(1),c_lastone(1)],[c_pre(2),c_lastone(2)].The '-'.'color',C(i,:),'linewidth'.1);
        end
    end
end
plot(customer(:,1),customer(:,2),'ro'.'linewidth'.1); hold on; plot(vertexs(1.1),vertexs(1.2),'s'.'linewidth'.2.'MarkerEdgeColor'.'b'.'MarkerFaceColor'.'b'.'MarkerSize'.10);
end

Copy the code

4. Operation results

Matlab version and references

1 matlab version 2014A

[1] Yang Baoyang, YU Jizhou, Yang Shan. Intelligent Optimization Algorithm and Its MATLAB Example (2nd Edition) [M]. Publishing House of Electronics Industry, 2016. [2] ZHANG Yan, WU Shuigen. MATLAB Optimization Algorithm source code [M]. Tsinghua University Press, 2017.