cvrp.gms : Capacitated Vehicle Routing Problem

Description

This model tries to find the optimal route for a fleet of vehicles with limited capacities
by minimizing the total distance travelled while fulfilling the demands of a set of customers.
The model employs two methods to eliminate subtours in the route, viz, i) miller-tucker-zemlin,
ii) dantzig-fulkerson-johnson. The data is from VRPLIB and embedded code Python is used to read the data.
Augerat 1995 - Set A | a-n32-k05, http://www.vrp-rep.org/datasets/item/2014-0000.html


Large Model of Type : MIP


Category : GAMS Model library


Main file : cvrp.gms   includes :  a-n32-k05.xml

$title Capacitated Vehicle Routing Problem (CVRP,SEQ=435)

$onText
This model tries to find the optimal route for a fleet of vehicles with limited capacities
by minimizing the total distance travelled while fulfilling the demands of a set of customers.
The model employs two methods to eliminate subtours in the route, viz, i) miller-tucker-zemlin,
ii) dantzig-fulkerson-johnson. The data is from VRPLIB and embedded code Python is used to read the data.
Augerat 1995 - Set A | a-n32-k05, http://www.vrp-rep.org/datasets/item/2014-0000.html


This formulation is described in detail in:
G. B. Dantzig, J. H. Ramser, (1959) The Truck Dispatching Problem. Management Science 6(1):80-91.

Keywords: mixed integer linear programming, capacitated vehicle routing problem,
          miller-tucker-zemlin, dantzig-fulkerson-johnson
$offText


$eolCom !!

$if not set TOURELIMINATE   $set TOURELIMINATE  mtz !! subtour elimination methods [mtz, dfj]
$if not set TIMELIMIT       $set TIMELIMIT      10  !! total time limit
$if not set DEBUG           $set DEBUG          0   !! DEBUG=1 activates additional output

* Vehicle Flow Formulation
Set ii              superset of customer nodes
    i(ii)           subset of customer nodes
    kk              superset of vehicles     
    k(kk)           subset of vehicles 
    depot(ii)       depot nodes
    ijk(ii,ii,kk)   allowed arcs for vehicles;
     
Alias (ii,jj),(i,i2,j);

Parameter
    demand(ii)       demand of consumer at node
    distance(ii,jj)  distance between nodes
    capacity(kk<)    capacity of each vehicle;


* Reading the data from the xml file
$onEmbeddedCode Python:
import pandas as pd
from math import dist, ceil

nodes = pd.read_xml('a-n32-k05.xml', xpath='./network/nodes/*', parser='etree')
nodeId = nodes['id'].values
depotId = nodes[nodes['type']== 0]['id'].values
capacity = pd.read_xml('a-n32-k05.xml', xpath='./fleet/vehicle_profile', parser='etree')["capacity"].values[0]
demandFrame = pd.read_xml('a-n32-k05.xml', xpath='./requests/request', parser='etree')
demand = demandFrame[['node', 'quantity']].values
nvehicle = ceil(demandFrame.quantity.sum()/capacity)

distArray = nodes[['cx', 'cy']].to_numpy()
df = list()
for i in range(len(distArray)):
    for j in range(len(distArray)):
        df.append(('n'+str(i+1) ,'n'+str(j+1), dist(distArray[i], distArray[j])))

ii = [('n'+str(i)) for i in nodeId]
kk = [('k'+str(i)) for i in range(1, nvehicle+1)]
depot = [('n'+str(i)) for i in depotId]
demand = [('n'+str(int(i[0])), i[1]) for i in demand]
capacity = [('k'+str(i), float(capacity)) for i in range(1, nvehicle+1)]

gams.set('ii', ii)
gams.set('demand', demand)
gams.set('distance', df)
gams.set('depot', depot)
gams.set('capacity', capacity)
$offEmbeddedCode ii, demand, distance, depot, capacity

* use only 16 out of 32 nodes, full problem becomes computationally challenging

i(ii)$[ord(ii) <= 16]   = yes;     !! select first 16 nodes includind depot (n1)
k(kk)$[ord(kk) <= 3]    = yes;     !! select first 3 out of 5 vehicles


Binary Variable
    X(ii,jj,kk) 'arc from node i to node j is used by vehicle k';

Free variable
    TOTDIST     'total distance';

Equation
    eq_tot_dist             'Objective function'
    eq_node_balance(ii,kk)  'Every vehicle must leave the node once entered'
    eq_enter_once(ii)       'A node can only be visited once'
    eq_leave_depot(kk,ii)   'All the vehicles should leave the depot'
    eq_capacity(kk)         'Capacity limit of each vehicle must be met';
    

eq_tot_dist..                       TOTDIST =e= sum(ijk(i,j,k), distance(i,j)*X(i,j,k));

eq_node_balance(j,k)..              sum(i$ijk(i,j,k), X(i,j,k)) =E= sum(i$ijk(j,i,k), X(j,i,k));

eq_enter_once(j)$(not depot(j))..   sum((i,k)$ijk(i,j,k), X(i,j,k)) =E= 1;

eq_leave_depot(k,i)$depot(i)..      sum(j$(not depot(j) and ijk(i,j,k)), X(i,j,k)) =E= 1;

eq_capacity(k)..                    sum((i,j)$(not depot(j) and ijk(i,j,k)), demand(j)*X(i,j,k)) =L= capacity(k);


option limrow=0, limcol=0, solprint=off, solvelink=5;
$ifE %DEBUG%>1 option limrow=1e6, limcol=1e6, solprint=on;

* populate allowed arcs, no loops allowed
ijk(i,j,k)$[not sameas(i,j)] = yes;


$ifThenI.TOURELIMINATE %TOURELIMINATE%==mtz
  !! Miller - Tucker - Zemlin subtour elimination
  Positive Variable P(ii) 'MTZ variable that determines the order of a tour and puts a bound on the current load';
  
  !! bounds for MTZ Variable P
  P.fx(ii)$depot(ii)  = 0;
  P.up(ii)            = card(jj)-1;
  
  Equation eq_MTZ(ii,jj) 'Miller, Tucker and Zemlin subtour elimination';
  
  eq_MTZ(i,j)$[not sameas(i,j)].. P(i) - P(j) =l= card(j) - card(j)*sum(k$ijk(i,j,k), X(i,j,k)) - 1 + card(j)$(depot(j));
  
  Model vrp_mtz 'Miller - Tucker - Zemlin subtour elimination' /all/;

  vrp_mtz.reslim      = %TIMELIMIT%;

  Solve vrp_mtz min TOTDIST use MIP;
  display x.l;

$elseIfI.TOURELIMINATE %TOURELIMINATE%==dfj
  !! Dantzig - Fulkerson - Johnson subtour elimination
  Set stCut           possible subtour elimination cuts   /c1*c1000/ 
      a(stCut)        active cuts
      tour(ii,jj,kk)  possible subtour
      sn(jj,kk)       nodes visited by subtour;
       
  Singleton Set
      cur(stCut)      current cut /c1/;
      
  Parameter
      cc(stCut,ii,jj) cut coefficient
      rhs(stCut);
      
  Equation eq_defdfj(stCut,kk)  'Equations that define cuts. For a subset of nodes S that defines a subtour, the number of arcs between those nodes must be <= card(S)-1';
  
  eq_defdfj(a,k).. sum((i,j)$ijk(i,j,k), cc(a,i,j)*X(i,j,k)) =L= rhs(a);
  
  
  Model vrp_dfj / all /;
  
  vrp_dfj.reslim = %TIMELIMIT%;

  a(stCut)  = no;
  cc(a,i,j) = 0;
  rhs(a)    = 0;
  Scalar foundoptimal / 0 /;
  Singleton Set last(ii,kk);
  
  Repeat
  
      Solve vrp_dfj min TOTDIST use MIP;
      abort$(vrp_dfj.modelStat <> %modelStat.optimal% and vrp_dfj.modelStat <> %modelStat.integerSolution% ) 'problems with MIP solver';
      X.l(i,j,k) = round(X.l(i,j,k));    
   
      !! The following block investigates all tours in the current solution.
      !! If there are illegal subtours, add a cut for the next solve.
      tour(i,j,k) = no;
      sn(i, k) = depot(i); !!start building tour from depot    
      while(card(sn),
          foundoptimal = 1; !! assume we are optimal (but reset to 0 if we detect illegal subtour)
          loop(k$(sum(sn(i,k), 1)),
              last(sn(i,k)) = 1;
              while(sum((i,j)$last(i,k), X.l(i,j,k)),     !! while there are arcs in solution for vehicle k from a node i which is in sn 
                  loop((i,j)$(last(i,k) and X.l(i,j,k)),  !! loop over those arcs 
                      tour(i,j,k)           = yes;        !! add arc to tour of vehicle k
                      X.l(i,j,k)            = 0;          !! remove arc from solution
                      sn(j, k)              = yes;        !! add destination node j of ark to set sn
                      last(j,k)             = yes;
                  );
              );
  
              !! if sn(*,k) does not contain depot, the nodes in sn(*,k) form an illegal subtour. 
              if(not sum(sn(depot,k),1), !! if depot not in sn
                  foundoptimal = 0;
                  $$ifThenE %DEBUG%>=1
                     put_utility$%DEBUG% 'log' / '!!!!!! illegal subtour with ' (sum(sn(j,k), 1)):0 ' nodes for vehicle ' k.tl:0 ' detected';
                     put_utility$%DEBUG% 'log' / '!!!!!! k = ' k.tl:0; 
                     put_utility$%DEBUG% 'log' / '!!!!!! sn(j,k) = ';
                     loop(sn(j,k), put_utility$%DEBUG% 'log' / '!!!!!!           ' j.tl:0;);
                  $$endIf
                  a(cur)                             = yes;                  !! activate current cut
                  rhs(cur)                           = sum(sn(j,k), 1) - 1;  !! compute rhs as card(S)-1
                  cc(cur,i,j)$(sn(i,k) and sn(j, k)) = 1;                    !! set coeefficients for cut to 1
                  cur(stCut+1)                       = cur(stCut);           !! advance current cut to next cut in cut set
                  abort$(card(cur)=0) 'Out of subtour cuts, enlarge set stCut';
              );
  
              sn(j, k) = no;
              loop((i,j)$(sum(sn(i2,k),1) = 0 and x.l(i,j,k)), sn(i, k) = yes);  !! find node in tours of k which has not yet been visited in subtour detection.
          );
      );
  Until foundoptimal or TimeElapsed > %TIMELIMIT%;
  X.l(tour) = 1; !! restore final solution
  display X.l;
$else.TOURELIMINATE
  $$abort Invalid value for --TOURELIMINATE. Allowed are 'mtz' and 'dfj'  
$endIf.TOURELIMINATE