VRP with working time windows (and not-working dvar "r") (r = number of palettes in truck at that moment)

camions.dat

/*********************************************
 * OPL 12.5 Data
 * Author: crash
 * Creation Date: Mar 6, 2013 at 3:38:56 PM
 *********************************************/

// 17 clients, WGS84
//Clients = [
//<"Poitiers",  		46.589069,0.340576>,
//<"La Rochelle",			46.1760-1.145325>,
//<"Angouleme",			45.66010.140076>,
//<"Perigueux",			45.18978,0.722351>,
//<"Limoges",				45.830713,1.263428>,
//<"Clermont-Ferrand",	45.784764,3.109131>,
//<"Chateauroux",			46.80758,1.672668>,
//<"Bourges",				47.073863,2.384033>,
//<"Blois",				47.585789,1.340332>,
//<"Angers",				47.476376,-0.563049>,
//<"Nantes",				47.215837,-1.554565>,
//<"Cholet",				47.062638,-0.884399>,
//<"La Roche-sur-Yon",	46.675826,-1.422729>,
//<"Saintes",				45.73686,-0.6427>,
//<"Tours",				47.394631,0.686646>,
//<"Orleans",				47.908978,1.895142>,
//<"Le Mans",				48.004625,0.192261>
//];
//

// 17 nodes, first node is depot.
// <Name, Lam93-X, Lam93-Y, TimeWin-Min, TimeWin-Max>
AllNodes = {
<"Poitiers"     		,496421,6613320,360,35000>,
<"La Rochelle"  		,380354,6572413,361,22000>,
<"Angouleme"    		,477343,6510763,362,22000>,
<"Perigueux"    		,521155,6457075,360,22000>,
<"Limoges"      		,565198,6527156,360,22000>,
<"Clermont-Ferrand"     ,708479,6520577,360,20000>,
<"Chateauroux"  		,598786,6635010,360,22000>,
<"Bourges"      		,653259,6663917,360,22000>,
<"Blois"        		,575264,6721912,360,22000>,
<"Angers"       		,431735,6714499,360,22000>,
<"Nantes"       		,355488,6689443,360,22000>,
<"Cholet"       		,405299,6669738,360,22000>,
<"La Roche-sur-Yon"     ,362092,6628992,360,22000>,
<"Saintes"      		,416833,6521784,360,22000>,
<"Tours"        		,525525,6701921,360,22000>,
<"Orleans"      		,617461,6757086,360,22000>,
<"Le Mans"      		,490662,6770859,360,22000>
};

Dist = [ // distance pairs, in meter, /1000 als kostbepaling?
123,104,158,110,231,104,164,134,120,160,107,135,121,93,187,157,
114,182,190,332,227,287,245,151,119,100,59,62,194,300,227,69,89,
231,173,233,232,208,216,174,165,61,197,283,260,82,197,194,245,270,
272,285,242,234,122,244,315,315,143,112,162,195,230,265,214,227,
148,179,235,254,158,153,241,337,391,337,362,291,257,253,331,61,
90,184,249,196,236,214,99,123,173,97,227,298,248,293,275,133,99,
194,143,222,177,232,255,53,54,97,80,51,110,193,94,190,81,53,60,178,
170,270,157,59,148,124,229,132,120,178,285,191,210,309,259,107,77,127
];

// Livraisons. Noeud 1 attend 1 palette de noeud 14
// <ToNode, FromNode, NbPalettes>
Deliveries = {
<1,14,1>,
<1,10,3>,
<2,5,7>,
<2,4,4>
};

// gewicht van paletten blijft hetzelfde

// 17 trucks? initial value?
// initial: 27 palettes, 33 max, 20 ton max, max tijd op route
// Tuple = distinct!!
Camions = {
	<1,33,20,600>,
	<2,33,20,600>,
	<3,33,20,600>,
	<4,33,20,600>,
	<5,33,20,600>,
	<6,33,20,600>,
	<7,33,20,600>,
	<8,33,20,600>,
	<9,33,20,600>,
	<10,33,20,600>,
		
	<9,33,20,600>,
	
	<15,33,20,600>,
	<16,33,20,600>,
	<17,33,20,600>,
	<18,33,20,600>,
	<19,33,20,600>
};

// laadtijd per node voor elke camion: werken met camiontypes!
//LoadTime = [];



/*
Location = {A, B, C, D, E, F, G, H};
TruckTypes = { SmallTruck, BigTruck };

Spokes = {A, B, C, D, E, F }; // 
Hubs = {G, H}; // laadpunten

Spoke = [<360, 1080>,<400,1150>,<380,1200>,
         <340,900>,<420,800>,<370,1070>];

TruckTypeInfos = [<400,10,55>,<700,15,45>];

LoadTime = [[30,55],[35,50]];

Routes = {<A,G,200>,<A,H,50>,<B,G,120>,<B,H,100>,
          <C,H,110>,<D,G,70>,<D,H,100>,<E,G,120>,
          <E,H,100>,<F,H,105>};

Shipments = {<A,B,300>,<A,C,250>,<A,D,350>,<A,E,145>,<A,F,300>,
             <B,A,185>,<B,C,200>,<B,D,221>,<B,E,263>,<B,F,197>,
             <C,A,143>,<C,B,178>,<C,D,258>,<C,E,221>,<C,F,106>,
             <D,A,75>,<D,B,135>,<D,C,245>,<D,E,283>,<D,F,155>,
             <E,A,123>,<E,B,234>,<E,C,143>,<E,D,78>,<E,F,107>,
             <F,A,201>,<F,B,157>,<F,C,169>,<F,D,212>,<F,E,104>};
 */

vrp.mod

/*********************************************
 * OPL 12.5 Model
 * Author: crash
 * Creation Date: Mar 6, 2013 at 3:33:35 PM
 *********************************************/

using CP;

// Nodes
tuple node { string name; int x; int y; int tmin; int tmax; }
{node} AllNodes = ...;
int N = card(AllNodes);
range nodes = 0..N-1;

// Edges
tuple edge {int i; int j;}
setof(edge) AllEdges = { <i,j> | i,j in nodes : i != j }; // not filtered!
setof(edge) OrderedEdges = { <i,j> | ordered i,j in nodes};
int Dist[OrderedEdges] = ...;
 
// Deliveries
tuple delivery {int destination; int from; int numpalettes;}
{delivery} Deliveries = ...;

// Remove superfluous thingies
sorted {int} clients = { d.from | d in Deliveries : d.numpalettes > 0 } union { d.destination | d in Deliveries : d.numpalettes > 0 };
sorted {int} cities = {0} union clients;
int numCities = card(cities);

int demand[t in cities] = sum (d in Deliveries : d.destination == t) d.numpalettes;
int outgoing[t in cities] = sum (d in Deliveries : d.from == t) d.numpalettes;
int tmin[t in cities] = item(AllNodes, t).tmin;
int tmax[t in cities] = item(AllNodes, t).tmax;

setof(edge) Edges = { <i,j> | i,j in cities : i != j }; // filtered! :)
//
//execute {
// writeln(AllEdges);
// writeln(tmin);
// writeln(tmax);
//}  

// Camions
tuple camion {key int id;int capacity; int maxweight; int maxtime;}
{camion} Camions = ...;
 
// Travel time
int vitesse = 65; // km/h
int tempstrajet[<i,j> in OrderedEdges] = ftoi(round(Dist[<i,j>]/vitesse*60));
int tempsservice = 30;
 
// Decision variables
dvar boolean x[Edges,Camions];
dvar boolean y[Camions,cities];
dvar int arrival[Camions][cities] in 0..numCities-2;
dvar int+ tarrival[Camions][cities];
dvar int+ quantity[cities][Camions];
dexpr int r[j in cities][k in Camions] = sum(i in cities) (y[k][i] == 1 ? ((arrival[k][i] <= arrival[k][j]) ? quantity[i][k] : 0) : 0);

//dvar interval timeofdelivery[i in clients][m] in i.tmin..i.tmax size 30;


execute {
//  cp.param.NoOverlapInferenceLevel = "Medium";
//  cp.param.ElementInferenceLevel   = "Low";
  //cp.param.TemporalRelaxation      = "Off";
  cp.param.TimeMode                = "ElapsedTime";
  cp.param.TimeLimit               = 10;
};

//tuple Subtour { 
//  int ssize;
//	camion truck;
//	int subtour[cities];
//};
//{Subtour} subtours = ...;
 
minimize sum (<i,j> in Edges, m in Camions) Dist[<minl(i,j), maxl(i,j)>]*x[<i,j>,m];

subject to {
	forall (k in Camions, j in cities)
	  	(sum(i in cities : j!=i) x[<i,j>][k]) == y[k][j]; 
	
	forall (k in Camions, j in cities)
	  	(sum(i in cities : j!=i) x[<j,i>][k]) == y[k][j];

	forall (i in clients, k in Camions)
		(quantity[i][k] >= 1) => y[k][i] == 1;

 	forall (i in cities)
 	  	((demand[i] + outgoing[i]) > 0) => sum(k in Camions) y[k][i] >= 1;

	// continuité  	
  	forall (p in cities, k in Camions)
  	  	((sum(i in cities : i!=p) x[<i,p>][k]) - (sum(j in cities : j!=p) x[<p,j>][k])) == 0;
	  	  	
  	// chaque véhicule entre et sort qu'un seul fois du dépôt
//	forall (k in Camions)
//	  	(sum(j in clients) x[<0,j>][k]) == 1;
//	  	  	
//  	forall (k in Camions)
//  	  	(sum(i in clients) x[<i,0>][k]) == 1;
  	  	
  	forall (k in Camions, <i,j> in Edges : i!=0 && j != 0)
  	  	(x[<i,j>][k] == 1) => (arrival[k][i] + 1 == arrival[k][j]);
	  	  	
	// For each client, the aggregated quantity of all trucks must match the demand
	forall (i in cities)
	  	sum(k in Camions) quantity[i][k] == demand[i];	
  	  		  	
  	// The amount of palettes the truck is going to deliver must be inferior|equal to his capacity
  	// (perhaps include dvar (c12)
  	forall (m in Camions)
  	  	sum(i in cities) quantity[i][m] <= m.capacity;
  	  	
//	// respect timewindows
	forall(j in clients, k in Camions)
	  (x[<0,j>][k] == 1) => (tarrival[k][j] >= (tempstrajet[<0,j>] + tmin[0]));

  	forall (<i,j> in Edges : i!=0, k in Camions)
  		(x[<i,j>][k] == 1) => ((tarrival[k][i] + tempstrajet[<minl(i,j), maxl(i,j)>] + (15 + 3*(quantity[i][k] + r[j][k] - r[i][k]))) <= tarrival[k][j]);
//  	
  	// min max timeframe for entry at client
  	forall(i in cities, k in Camions)
  		(y[k][i] ==1) => ((tmin[i] <= tarrival[k][i]) && (tarrival[k][i] <= tmax[i]));

  	// truck can't be longer than e.g. 12h on route
  	forall (m in Camions)
		m.maxtime >= sum (<i,j> in Edges) x[<i,j>,m]*(tempstrajet[<minl(i,j),maxl(i,j)>]+tempsservice);
//	
	// Delivery at node 2 must be later than delivery at node 1
	forall (<i,j> in Edges, d in Deliveries : d.destination == j && d.from == i, m in Camions)
		(y[m][d.from] == 1 && y[m][d.destination] == 1) => (arrival[m][d.from] < arrival[m][d.destination]);

};

//execute {
//  	function recurse_results(k, i) {
//  	  	for (var t in y[k]) {
//  	  		if (y[k][t] == 1) {
//  	  		  	var temp = new Array(i,t);
//  	  		  	if (x[temp][k] == 1) {
//  	  		  	  	recurse_results(k, t);
//  	  		  	  	writeln(t + " -> ");
//  	  		  		break;
//          		}  	  		  	  
//  	  			
//       		}  	  		    
//      	}
//   	}  	  
//  
//  
//	for (var k in arrival) {
//	  	var str = "Tour for camion "+k.id+": 0 -> ";
//	  	recurse_results(k, 0);
//		writeln(str + "0");
//	}  	  
//};