I am trying to solve a Vehicle Routing Problem with multiple pickups and dropoffs with multiple products carried by just one car. After solving this problem I am going to extend to multiple types of cars as well.
One special setting is that it has a starting point and an ending point which needs not be same. I assumed to be different and setted 1 and n to be the dummy nodes for start and end.
I partially used an example TSP code provided by IBM to solve subtour problem, and got help from internet to print out the optimal tour.
Since I need to find a optimal path that go through all the points. This is NP-hard. But as a first time using ILog, I would like to solve this problem using MIP for practice purpose.
I am having trouble with keeping track of picked up products and droped off products in each arc.
I am trying to minimize the cost of transportation which I set to be
// Decision variables
dvar boolean x[Edges]; //car goes through this arc?
dvar boolean y[Edges][Products]; //at each e, currently loaded products in the car
dvar boolean z[Cities][Products]; //at each cities, what products to load or unload?
// Cost function
// at each arc that car goes through (distance + sum(products in the car*their weights))
dexpr float cost = sum (e in Edges) x[e] *
(dist[e] + (sum (p in Products) y[e][p] * weight[p]));
y is a variable that associates each arc with currently loaded products. z accounts for what to load or unload in each nodes. Since there is only one car I don't think z is actually needed, but for a extension with multiple cars, I think this will be a good thing to have.
If some of these dvars are not necessary, please give me some insight! Below are setups.
// Cities
int n = ...;
range Cities = 1..n;
range Cities1 = 2..n-1; //just for start/end restriction
range Pickups = 2..3;
range Dropoffs= 4..n-1;
// Edges -- sparse set
tuple edge {int i; int j;}
setof(edge) Edges = {<i,j> | ordered i,j in Cities};
int dist[Edges] = ...;
// Products
int p = ...;
range Products = 1..p;
float Pickup[Cities][Products] = ...;
float Dropoff[Cities][Products] = ...;
float weight[Products] = ...;
// Products in pickup and dropoff sums up to equal amount
assert
forall(p in Products)
sum(o in Cities) Pickup[o][p] == sum(d in Cities) Dropoff[d][p];
tuple Subtour { int size; int subtour[Cities]; }
{Subtour} subtours = ...;
Any help on below restrictions will be very helpful. Especially on keeping track of loaded products along the route.
// Objective
minimize cost;
subject to {
// Each city is linked with two other cities
forall (j in Cities1)
sum (<i,j> in Edges) x[<i,j>] + sum (<j,k> in Edges) x[<j,k>] == 2;
// Must start at node 1 and end at node n
sum (e in Edges : e.i==1) x[e] == 1;
sum (e in Edges : e.j==n) x[e] == 1;
// no product remains at 1,n (not necessary?)
sum (p in Products, e in Edges : e.i==1) y[e][p] == 0;
sum (p in Products, e in Edges : e.j==n) y[e][p] == 0;
sum (p in Products) z[1][p] == 0;
sum (p in Products) z[n][p] == 0;
// must pickup all
forall (p in Products) {
sum(i in Pickups) z[i][p] == sum(i in Cities) Pickup[i][p];
sum(i in Dropoffs) z[i][p] == sum(i in Cities) Dropoff[i][p];
}
forall (i in Pickups, p in Products)
z[i][p] <= Pickup[i][p];
//tried to keep track of picked ups, but it is not working
forall (i in Pickups, j,k in Cities, p in Products : k < i < j)
y[<i,j>][p] == y[<k,i>][p] + z[i][p];
// forall (j in Cities, p in Products)
// ctDemand:
// sum(<i,j> in Edges) y[<i,j>][p] + sum(<j,i> in Edges) y[<j,i>][p] == z[j][p];
// tried keeping track of dropoffs. It is partially working but not sure of it
forall (i, k in Cities, j in Dropoffs, p in Products : i < j < k) (y[<i,j>][p] == 1)
=> y[<j,k>][p] == y[<i,j>][p] - Dropoff[j][p];
// Subtour elimination constraints.
forall (s in subtours)
sum (i in Cities : s.subtour[i] != 0)
x[<minl(i, s.subtour[i]), maxl(i, s.subtour[i])>] <= s.size-1;
};
And post processing to find the subtours
// POST-PROCESSING to find the subtours
// Solution information
int thisSubtour[Cities];
int newSubtourSize;
int newSubtour[Cities];
// Auxiliar information
int visited[i in Cities] = 0;
setof(int) adj[j in Cities] = {i | <i,j> in Edges : x[<i,j>] == 1} union
{k | <j,k> in Edges : x[<j,k>] == 1};
execute {
newSubtourSize = n;
for (var i in Cities) { // Find an unexplored node
if (visited[i]==1) continue;
var start = i;
var node = i;
var thisSubtourSize = 0;
for (var j in Cities)
thisSubtour[j] = 0;
while (node!=start || thisSubtourSize==0) {
visited[node] = 1;
var succ = start;
for (i in adj[node])
if (visited[i] == 0) {
succ = i;
break;
}
thisSubtour[node] = succ;
node = succ;
++thisSubtourSize;
}
writeln("Found subtour of size : ", thisSubtourSize);
if (thisSubtourSize < newSubtourSize) {
for (i in Cities)
newSubtour[i] = thisSubtour[i];
newSubtourSize = thisSubtourSize;
}
}
if (newSubtourSize != n)
writeln("Best subtour of size ", newSubtourSize);
}
main {
var opl = thisOplModel
var mod = opl.modelDefinition;
var dat = opl.dataElements;
var status = 0;
var it =0;
while (1) {
var cplex1 = new IloCplex();
opl = new IloOplModel(mod,cplex1);
opl.addDataSource(dat);
opl.generate();
it++;
writeln("Iteration ",it, " with ", opl.subtours.size, " subtours.");
if (!cplex1.solve()) {
writeln("ERROR: could not solve");
status = 1;
opl.end();
break;
}
opl.postProcess();
writeln("Current solution : ", cplex1.getObjValue());
if (opl.newSubtourSize == opl.n) {
// This prints the tour as a cycle
var c = 1; // current city
var lastc = -1; // city visited right before C
write(c);
while (true) {
var nextc = -1; // next city to visit
// Find the next city to visit. To this end we
// find the edge that leaves city C and does not
// end in city LASTC. We know that exactly one such
// edge exists, otherwise the solution would be infeasible.
for (var e in opl.Edges) {
if (opl.x[e] > 0.5) {
if (e.i == c && e.j != lastc) {
nextc = e.j;
break;
}
else if (e.j == c && e.i != lastc) {
nextc = e.i;
break;
}
}
}
// Stop if we are back at the origin.
if (nextc == -1) {
break;
}
// Write next city and update current and last city.
write(" -> ", nextc);
lastc = c;
c = nextc;
}
opl.end();
cplex1.end();
break; // not found
}
dat.subtours.add(opl.newSubtourSize, opl.newSubtour);
opl.end();
cplex1.end();
}
status;
}
Here is a sample dataset that I created. I hope my explanation makes sence to everyone! Thank you very much!!
n = 10;
dist = [
633
257
91
412
150
80
134
259
505
390
661
227
488
572
530
555
289
228
169
112
196
154
372
262
383
120
77
105
175
476
267
351
309
338
196
63
34
264
360
29
232
444
249
402
495
];
// Products
p = 8;
Pickup = [
// 1,2,3,4,5,6,7,8 products
[0 0 0 0 0 0 0 0],//city1
[0 1 0 1 0 1 1 0],//city2
[1 0 1 0 1 0 0 1],//city3
[0 0 0 0 0 0 0 0],//city4
[0 0 0 0 0 0 0 0],//city5
[0 0 0 0 0 0 0 0],//city6
[0 0 0 0 0 0 0 0],//city7
[0 0 0 0 0 0 0 0],//city8
[0 0 0 0 0 0 0 0],//city9
[0 0 0 0 0 0 0 0] //city10
];
Dropoff = [
[0 0 0 0 0 0 0 0],
[0 0 0 0 0 0 0 0],
[0 0 0 0 1 0 0 0],
[0 0 0 0 0 1 0 0],
[0 0 0 0 0 0 1 0],
[1 0 0 0 0 0 0 1],//city6
[0 1 0 0 0 0 0 0],//city7
[0 0 1 0 0 0 0 0],//city8
[0 0 0 1 0 0 0 0],//city9
[0 0 0 0 0 0 0 0] //city10
];
weight = [1, 2, 3, 4, 5, 6, 7, 8];
OK, I have been thinking for a long time to solve this problem. And at last, I have solved it with reasonable amount of running time. I will share it for someone who perhaps wants to dispute what i have done or improve running time. Meanwhile, it will be awesome if someone could verify what i have done is indeed correct logically.
Oh, I have solved it an extension that i have planned with extra several constraints. I hope my comments are enough for understanding what i have done.
// Constants
int TotalTimeLimit = ...;
int LoadTime = ...;
int UnloadTime = ...;
int TransitCost = ...;
int MaxTransit = ...;
// Cities
int n = ...;
int pickupN= ...; //4
range Cities = 1..n;
range Cities1 = 2..n-1; //just for start/end restriction
range Pickups = 2..pickupN+1; //2-5
range Dropoffs= pickupN+2..n-1; //7-17
// Edges -- sparse set
tuple edge {int i; int j;}
setof(edge) Edges = {<i,j> | i,j in Cities};
int dist[Edges] = ...;
int time[Edges] = ...;
// Products
int P = ...;
range Products = 1..P;
int Pickup[Cities][Products] = ...;
int Dropoff[Cities][Products] = ...;
int weight[Products] = ...;
int volume[Products] = ...;
// Products in pickup and dropoff sums up to equal amount
assert
forall(p in Products)
sum(o in Cities) Pickup[o][p] == sum(d in Cities) Dropoff[d][p];
//Trucks
{string} Type = ...;
int Max[Type] = ...;
int c = sum(t in Type) Max[t];
range Cars= 1..c;
int VolumeCapacity[Cars] = ...;
int WeightCapacity[Cars] = ...;
int NumberCapacity[Cars] = ...;
int FixedCost[Cars] = ...;
int forbid[Cities][Cars] = ...;
// Decision variables
dvar boolean x[Cars][Edges];
dvar boolean y[Cars][Edges][Products];
dvar boolean z[Cars][Cities][Products];
dvar boolean isCarUsed[Cars];
// Cost function
dexpr float cost = sum (c in Cars) (isCarUsed[c]*FixedCost[c] +
(sum(e in Edges) (x[c][e]*(dist[e] + TransitCost) + (sum (p in Products) y[c][e][p] * weight[p]))))
- 30 - sum(p in Products) weight[p];
// Objective
minimize cost;
subject to {
// Total pickups for each p equal the total sum of each colum of z in pickups
forall (p in Products)
sum(i in Pickups, c in Cars) z[c][i][p] == sum(i in Pickups) Pickup[i][p];
// For each node, each car can pick at most the node's stock
forall(i in Pickups, p in Products, c in Cars)
z[c][i][p] <= Pickup[i][p];
// For each car, picked up products should be smaller than car's capacities
forall (c in Cars, e in Edges) {
sum(p in Products) y[c][e][p] <= NumberCapacity[c];
sum(p in Products) y[c][e][p]*weight[p] <= WeightCapacity[c];
sum(p in Products) y[c][e][p]*volume[p] <= VolumeCapacity[c];
}
// For each car and product, its picked up amount should equal dropped off amount
forall (c in Cars, p in Products)
sum(i in Pickups) z[c][i][p] == sum(i in Dropoffs) z[c][i][p];
// Total dropoffs for each p equal the total sum of each column of z in dropoffs
forall (p in Products)
sum(i in Dropoffs, c in Cars) z[c][i][p] == sum(i in Dropoffs) Dropoff[i][p];
// For each node, each car can drop at most its needs
forall(i in Dropoffs, p in Products, c in Cars)
z[c][i][p] <= Dropoff[i][p];
// For each car, it cannot stay at one place
forall (c in Cars)
sum(i in Cities) x[c][<i,i>] == 0;
// Prevents looping, must start at node 1 and end at node n
forall (c in Cars) {
sum (e in Edges : e.i==n) x[c][e] == 0;
sum (e in Edges : e.j==1) x[c][e] == 0;
}
// For each car, it cannot go back and forth
forall (c in Cars, i,j in Cities)
x[c][<i,j>]+x[c][<j,i>] <= 1;
// Each city is linked with two other cities
forall (j in Cities1, c in Cars)
sum (<i,j> in Edges) x[c][<i,j>] - sum (<j,k> in Edges) x[c][<j,k>] == 0;
// If car is used, must start at node 1 and end at node n
forall(c in Cars) {
100*sum(e in Edges : e.i==1) x[c][e] >= sum (p in Products, i in Cities1) z[c][i][p];
100*sum(e in Edges : e.j==n) x[c][e] >= sum (p in Products, i in Cities1) z[c][i][p];
100*isCarUsed[c] >= sum(p in Products, i in Cities1) z[c][i][p];
}
forall (c in Cars) {
sum(e in Edges : e.i==1) x[c][e] <= 1;
sum(e in Edges : e.j==n) x[c][e] <= 1;
}
// For each car, it needs to cross nodes that picks or drops something
forall (c in Cars, i in Cities)
100*sum (j in Cities) x[c][<i,j>] >= sum (p in Products) z[c][i][p];
// For each car, its transit count <= maxTransit
forall (c in Cars)
sum(e in Edges) x[c][e] <= MaxTransit;
// For each car, and for each node, we need to check whether time took is <= totalTimeLimit
forall(c in Cars)
sum(e in Edges : e.i in Pickups || e.i in Dropoffs) x[c][e]*time[e]
+ LoadTime*sum(i in Pickups, p in Products) z[c][i][p] + UnloadTime*sum(i in Dropoffs, p in Products) z[c][i][p]
<= TotalTimeLimit;
// Certain type of cars cannot go to certain city
forall (c in Cars, i in Cities)
if (forbid[i][c] == 1)
sum(j,k in Cities) (x[c][<i,j>] + x[c][<k,i>]) == 0;
// Keeps culmulated pacakges through each arcs
forall (c in Cars, i in Pickups, p in Products)
sum(k in Cities) y[c][<i,k>][p] == sum(j in Cities) y[c][<j,i>][p] + z[c][i][p];
forall (c in Cars, i in Pickups, j in Cities)
100*x[c][<i,j>] >= sum(p in Products) y[c][<i,j>][p];
// MTZ subtour elimination constraints
forall (c in Cars, j in Dropoffs, i,k in Cities, p in Products : j != i && j != k)
y[c][<i,j>][p] - y[c][<j,k>][p] >= z[c][j][p] - 100*(1-x[c][<i,j>]);
};
Here is an example data. I made the setting such that we first pickup, and dropoff later. I also made 3 dummy variables for start and end, and transit position between pickups and dropoffs.
n = 12;
pickupN = 3;
dist = [
9999999
0
0
0
9999999
9999999
9999999
9999999
9999999
9999999
9999999
9999999
9999999
9999999
390
661
0
9999999
9999999
9999999
9999999
9999999
9999999
9999999
9999999
390
9999999
228
0
9999999
9999999
9999999
9999999
9999999
9999999
9999999
9999999
661
228
9999999
0
9999999
9999999
9999999
9999999
9999999
9999999
9999999
9999999
9999999
9999999
9999999
9999999
63
34
264
360
208
329
9999999
9999999
9999999
9999999
9999999
9999999
9999999
232
444
292
297
47
0
9999999
9999999
9999999
9999999
9999999
232
9999999
29
232
444
292
0
9999999
9999999
9999999
9999999
9999999
444
29
9999999
249
402
250
0
9999999
9999999
9999999
9999999
9999999
292
232
249
9999999
495
352
0
9999999
9999999
9999999
9999999
9999999
297
444
402
495
9999999
154
0
9999999
9999999
9999999
9999999
9999999
47
292
250
352
154
9999999
0
9999999
9999999
9999999
9999999
9999999
9999999
9999999
9999999
9999999
9999999
9999999
9999999
];
time = [
9999999
0
0
0
9999999
9999999
9999999
9999999
9999999
9999999
9999999
9999999
9999999
9999999
20
52
0
9999999
9999999
9999999
9999999
9999999
9999999
9999999
9999999
14
9999999
26
0
9999999
9999999
9999999
9999999
9999999
9999999
9999999
9999999
43
31
9999999
0
9999999
9999999
9999999
9999999
9999999
9999999
9999999
9999999
9999999
9999999
9999999
9999999
12
31
10
17
26
26
9999999
9999999
9999999
9999999
9999999
9999999
9999999
21
19
20
18
7
0
9999999
9999999
9999999
9999999
9999999
35
9999999
6
13
30
28
0
9999999
9999999
9999999
9999999
9999999
39
23
9999999
38
23
38
0
9999999
9999999
9999999
9999999
9999999
23
32
20
9999999
45
17
0
9999999
9999999
9999999
9999999
9999999
28
35
37
24
9999999
24
0
9999999
9999999
9999999
9999999
9999999
2
39
8
37
21
9999999
0
9999999
9999999
9999999
9999999
9999999
9999999
9999999
9999999
9999999
9999999
9999999
9999999
];
TotalTimeLimit = 150;
LoadTime = 10;
UnloadTime = 10;
TransitCost = 10;
//DropoffTransitCost = 10;
MaxTransit = 10;
// Products
P = 12;
Pickup = [
// 1,2,3,4,5,6,7,8 9 101112 products
[0 0 0 0 0 0 0 0 0 0 0 0],//city1 pickstart
[0 0 0 1 0 0 1 0 1 1 0 0],//city2
[1 1 0 0 0 0 0 1 0 0 1 0],//city3
[0 0 1 0 1 1 0 0 0 0 0 1],//city4
[0 0 0 0 0 0 0 0 0 0 0 0],//city5 pickend & dropstart
[0 0 0 0 0 0 0 0 0 0 0 0],//city6
[0 0 0 0 0 0 0 0 0 0 0 0],//city7
[0 0 0 0 0 0 0 0 0 0 0 0],//city8
[0 0 0 0 0 0 0 0 0 0 0 0],//city9
[0 0 0 0 0 0 0 0 0 0 0 0],//city10
[0 0 0 0 0 0 0 0 0 0 0 0],//city11
[0 0 0 0 0 0 0 0 0 0 0 0]//city12
];
Dropoff = [
[0 0 0 0 0 0 0 0 0 0 0 0],
[0 0 0 0 0 0 0 0 0 0 0 0],
[0 0 0 0 0 0 0 0 0 0 0 0],
[0 0 0 0 0 0 0 0 0 0 0 0],
[0 0 0 0 0 0 0 0 0 0 0 0],//city5 pickend & dropstart
[0 0 0 0 1 0 0 0 0 0 0 1],//city6
[1 0 0 0 0 0 0 0 0 0 0 0],//city7
[0 1 0 0 0 0 1 0 0 0 0 0],//city8
[0 0 0 0 0 1 0 0 1 0 0 0],//city9
[0 0 0 1 0 0 0 0 0 1 0 0],//city10
[0 0 1 0 0 0 0 1 0 0 1 0],//city11
[0 0 0 0 0 0 0 0 0 0 0 0]//city12 dropend
];
weight = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12];
volume = [20, 500, 30, 20, 60, 100, 300, 50, 10, 400, 1, 15];
//Cars
Type ={"SmallTruck", "MiddleTruck", "BigTruck"};
Max = #[
SmallTruck : 2,
MiddleTruck : 2,
BigTruck : 2
]#;
VolumeCapacity = [100, 100, 500, 500, 1000, 1000];
WeightCapacity = [10, 10, 30, 30, 50, 50];
NumberCapacity = [4, 4, 6, 6, 12, 12];
FixedCost = [100, 100, 200, 200, 300, 300];
forbid = [
[0 0 0 0 0 0],//city1
[0 0 0 0 0 0],//city2
[0 0 0 0 0 0],//city3
[0 0 0 0 0 0],//city4 pickend & dropstart
[0 0 0 0 0 0],//city5
[0 0 0 0 0 0],//city6
[0 0 0 0 0 0],//city7
[0 0 0 0 0 0],//city8
[0 0 0 0 0 0],//city9
[0 0 0 0 1 1],//city10
[0 0 0 0 0 0],//city11
[0 0 0 0 0 0]//city12
];