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