$title VRP Model using BCH Facility * enable commenting $eolcom // $onEcho > VRPCut_v2.inc * Read set required input data for sets and parameters for further processing * //$IF not exist testdaten.gdx $CALL GDXXRW.EXE Input=Testdaten.xlsx MaxDupeErrors=150 @Dateneingabe.txt $IF not exist Testdaten_small.gdx $CALL GDXXRW.EXE Input=Testdaten_small.xlsx MaxDupeErrors=150 @Dateneingabe.txt ** Required Sets ** Sets * Generic Set Definitions V set of vertices * Important Subsets Depot (V) set of depots N (V) set of customers M set of vehicles ; * Receive important sets //$GDXIN testdaten.gdx $GDXIN Testdaten_small.gdx $LOAD V, N, M, Depot $GDXIN * Important "artifical"/"indice" sets alias(V,i); alias(V,j); SET Arc(V, V) set of arcs between all vertices V / set.V.set.V /; SET Edge(V, V) set of edges between all vertices V / Depot.#V, Customer1.#V, Customer2.#V, Customer3.#V, Customer4.#V/ ; SET Delta_i_Plus (i,j) set of positive incident arcs / set.V.set.V /; Delta_i_Plus(i,j)$(sameas(i,j)) = no; SET Delta_i_Minus (i,j) set of negative incident arcs / set.V.set.V /; Delta_i_Minus(i,j)$(sameas(i,j)) = no; ** Required Parameters ** Parameters c_ij(V,V) costs per arc q(v) demand at node i GZ big number ; * Receive important parameters $GDXIN testdaten_small.gdx //$GDXIN testdaten.gdx $LOAD c_ij, q $GDXIN GZ = 1000; * See all input information in a GDX file. Execute_Unload 'input_vrpcut.gdx', i, j, q, Delta_i_Plus, Delta_i_Minus, N, V, Edge; Variables x_ij(i,j) path in graph taken by vehicle GOAL goal value; Binary Variable x_ij; Equations objective cost function constraint_1(V) leave constraint constraint_2(V) visit constraint constraint_3 leave depot constraint ; objective .. GOAL =e= sum((i,j), x_ij(i,j) * c_ij(i,j)); constraint_1(i)$(N(i)) .. sum(j$(Delta_i_Plus(i,j)), x_ij(i,j)) =e= 1; constraint_2(j)$(N(j)) .. sum(i$(Delta_i_Minus(i,j)), x_ij(i,j)) =e= 1; constraint_3.. sum(j$(Delta_i_Plus('Depot',j)), x_ij('Depot',j)) =l= card(M); Model vrpcut / all / ; $offEcho $include VRPCut_v2.inc Option Limrow=1000; Option mip=cplex; Option optcr=0; Option SysOut = On; vrpcut.optfile=1; Solve vrpcut using mip minimizing GOAL ; //abort$(vrpcut.modelStat <> 1) 'wrong solution'; execute_unload "result_vrpcut.gdx", x_ij, GOAL; * Cplex Option file to enable the cut generation call at every new integer point $onecho > cplex.opt *usercutcall vrpcuts lf=vrpcuts.log o=vrpcuts.lst mip=cplex UserIncbCall vrpcuts lf=vrpcuts.log o=vrpcuts.lst mip=cplex usercutcall subtourelim.gms $offecho $onEchoV > vrpcuts.gms $title Generate Cutting Planes for VRP * enable commenting $eolcom // ** Required Sets ** Sets V set of vertices // Generic Set Definitions Depot (V) set of depots N (V) set of customers ; * Receive important sets //$GDXIN testdaten.gdx $GDXIN Testdaten_small.gdx $LOAD V, Depot, N $GDXIN * Important "artifical"/"indice" sets alias(V,i,j,ix,ix2); Variables x_ij(i,j) path in graph taken by vehicle; Binary Variable x_ij; set cc 'cuts' /1*1000/; alias(cc,ccc); // allow 1000 cuts max singleton set curcut(cc) current cut always one element; set allcuts(cc) total cuts; * In order to communicate with the MIP solver we need certain conventions * 1. Cut matrix interface requirement with fixed names parameters rhs_c(cc) right hand side sense_c(cc) sense of cut 1=l 2=e 3=g numcuts 'number of cuts passed back' * 2. Nonzeros in cut matrix using the original variable name plus _c x_ij_c(cc, i, j); * read bch content $IF not exist bchout_i.gdx $goto nocuts $GDXIN bchout_i.gdx // opening file $LOAD x_ij set t tours /t1*t10/; singleton set tt(t) contains always one element: the current subtour; sets tour(i,j,t) tours visited(i) flag whether a city is already visited next(j) contains always one element: the to city subtours(t) subtours ; singleton sets from(i) contains always one element: the from city ; * initialize from(i)$(ord(i)=1) = yes; // turn first element on tt(t)$(ord(t)=1) = yes; // turn first element on curcut(cc)$(ord(cc)=1) = yes; // turn first element on // determine subtours scalar iteration; scalar rhs_value; iteration = 0; // init * initialize from(i) = no; tt(t) = no; subtours(t) = no; from(i)$(ord(i)=2) = yes; // turn first element on : First entry after depot tt(t)$(ord(t)=1) = yes; // turn first element on tour(i,j,t)=no; visited(i)=no; Display from; loop(i$(N(i)), // only for customers next(j)=no; // clearance // find city visited after city 'from' loop(j, next(j)$(x_ij.l(from,j)>0.5) = yes); // note: check x.l(from,j)=1 would be dangerous Display next; tour(from,next,tt) = yes; // store sequence in table visited(from) = yes; // mark city 'from' as visited from(j) = next(j); Display tour; // subtour: only check for sequences like c1 - c2 - c1 or c1 - c2 - c3 - c1. Depot checks not included. if (sum(visited(next)$(not Depot(next)),1)>0, // if already visited. Logic is as follows: if next (not yet marked as visited) is already in visited -> subtour (if not depot)! Display "Subtour identified -> enter loop"; subtours(tt) = yes; tt(t) = tt(t-1); // subtour identified. Swap to new tour ("Singleton Set"); Display "New tour added", tt; loop((ix,ix2)$(not visited(ix) and not visited(ix2) and x_ij.l(ix, ix2)>0.5), // find starting point of new subtour from(ix) = yes; Display from; break; // exit loop once finished ) ) ); display tour; display subtours; numcuts = sum(t, max(0,smax(subtours(t),1))); display numcuts; // else: introduce cut loop(t$(ord(t) <= numcuts), Display curcut; rhs_c(curcut)=-1; rhs_value = sum((ix,ix2)$(tour(ix,ix2,t) and (x_ij.l(ix,ix2) > 0.5)), 1) - 1; rhs_c(curcut) = rhs_value; sense_c(curcut) = 1; // must be less Display rhs_c; loop(tour(i,j,t), iteration = iteration+ 1; x_ij_c(curcut, i, j)$(x_ij.l(i,j) > 0.5) = 1; ); Display iteration; Display x_ij_c, rhs_c, sense_c; allcuts(curcut) = yes; // include this cut in set curcut(cc) = curcut(cc-1); ); display numcuts, x_ij_c, rhs_c, sense_c; abort$(not numcuts) 'no cuts found means solution should be accepted'; *unload found cuts for subtour elimination execute_unload 'nextcut', cc = cut, i, numcuts, x_ij_c, rhs_c, sense_c; $exit $label nocuts numcuts = 0; // return no cuts found information $offEcho $onEchoV > subtourelim.gms $title Supply the Cut found in Subtourcheck to Solver Set cut 'cuts' i 'vertices'; Parameter rhs_c(cut) 'cut rhs' sense_c(cut) 'the sense of the cuts 1 =l=, 2 =e=, 3 =g=' numcuts 'number of cuts passed back' x_ij_c(cut,i,i) 'coefficient of the x variables in the cut'; $if not exist nextcut.gdx $goto nocuts $gdxIn nextcut $load cut i rhs_c sense_c x_ij_c numcuts $gdxIn $call rm -f nextcut.gdx $exit $label nocuts numcuts = 0; $offEcho