$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