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

nodeId = nodes['id'].values
depotId = nodes[nodes['type']== 0]['id'].values
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);

\$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``````
GAMS Development Corp.
GAMS Software GmbH

General Information and Sales
U.S. (+1) 202 342-0180
Europe: (+49) 221 949-9170