I'm new to this platform and would appreciate some guidance on modeling subtour elimination constraints for a multi-vehicle Inventory Routing Problem (IRP) using MILP in OPL CPLEX.
I have attached my current code, which runs without errors; however, I am aware that the subtour elimination constraints are not correctly implemented yet.
Specifically, I'm struggling with how to properly formulate constraints (8) and (9) from the mathematical model within the .mod
file.
Any advice or examples on how to accurately incorporate these subtour constraints would be greatly appreciated!

int numCust = 10; //set customers
int V = 2; //set vehicles
int T = 2; //set time
int P = 3; //Set products
{int} warehouse = {0};
{int} arrware = {11};
{int} customeri = {1,2,3,4,5,6,7,8,9,10};
{int} vehicle = {1,2};
range period = 1..T;
range node = 0..numCust+1;
range product = 1..P;
{int} depNode = warehouse union customeri ;
{int} arrNode = customeri union arrware;
{int} fullNode = warehouse union customeri union arrware;
tuple arc1 {
int depart;
int arrive;
};
{arc1} ArcSet = {<dep,arr> | dep in depNode, arr in arrNode : dep != arr };
tuple arc2 {
int dep;
int arr;
int v;
};
{arc2} ArcVeh = {<dep,arr,v>| dep in depNode, arr in arrNode: dep!=arr, v
in vehicle};
//{IJ} ArcSet = { <i,j> | i in depNode, j in arrNode: i != j && j <= i + numD + numC} ;
tuple arc3 {
int i;
int v;
};
{arc3} CV = {<i,v>|i in customeri, v in vehicle};
tuple arc4 {
int p;
int i;
}
{arc4} PC = {<p,i>|p in product, i in customeri};
tuple arc5 {
int p;
int i;
int v;
}
{arc5} PIV = {<p,i,v>|p in product, i in customeri, v in vehicle};
range r1=1..-2+ftoi(pow(2,card(customeri)));
{int} nodes2 [k in r1] = {i | i in customeri: ((k div (ftoi(pow(2,(ord(customeri,i))))) mod 2) == 1)};
float c[fullNode][fullNode] = ...; //transportation cost
float c0[customeri] = [10,
10,
10,
10,
10,
10,
10,
10,
10,
10
];
float alpha = 0.01; //Discount factor
float vehcapQ = 2000; //CW's vehicles capacity
float inv0[product] = [2000,2000,2000]; //
float g[product][period] = [[20,20],[30,30],[40,40]];
//float r[PC][period] = ...;
//float u[PC][period]=...;
//float std[PC][period]=...;
//float beta = ...;
//float beta_quantile[product][period]=...;
//float s[PC][period]=...;
float deli[PC][period]=...;
dvar boolean x[ArcSet][vehicle][period]; // (i,j) visited by vehicle v in period t
dvar boolean y[CV][period];
dvar float inv[product][period];//
dvar float+ w[PC][period];
dvar float+ q[PC][vehicle][period];
dexpr float obj = sum(t in period)(sum(v in vehicle, <dep,arr> in ArcSet)c[dep][arr]*x[<dep,arr>][v][t]
+ alpha*sum(i in customeri)c0[i]*sum(p in product)w[<p,i>][t]);
minimize obj;
subject to {
constraint5:
forall(p in product, i in customeri, t in period){
sum(v in vehicle) q[<p,i>][v][t] + w[<p,i>][t] == deli[<p,i>][t];
}
//
constraint6:
forall( v in vehicle, t in period){
sum(i in customeri, p in product) q[<p,i>][v][t] <= vehcapQ;
}
constraint7:
forall(v in vehicle, t in period, i in customeri){
sum(p in product) q[<p,i>][v][t] <= vehcapQ * y[<i,v>][t];
}
constraint8:
forall(i in customeri, v in vehicle, t in period){
sum(j in depNode: j!=i) x[<j,i>][v][t] + sum(j in arrNode: j!=i) x[<i,j>][v][t] == 2*y[<i,v>][t];
}
constraint9:
forall(v in vehicle,t in period,k in r1, m in nodes2[k] ) {
sum(i in nodes2[k], j in nodes2[k]:j!=i) x[<i,j>][v][t] <= sum(i in nodes2[k])(y[<i,v>][t]) - y[<m,v>][t];
}
constraint10a: // For t=1, use inv0 as the previous inventory
forall(p in product) {
inv[p][1] == inv0[p] - sum(i in customeri, v in vehicle) q[<p,i>][v][1] + sum(i in customeri) w[<p,i>][1] + g[p][1];
}
//
constraint10b:
forall(p in product, t in period:t>1){
inv[p][t] == inv[p][t-1] - sum(i in customeri, v in vehicle) q[<p,i>][v][t]+ sum(i in customeri) w[<p,i>][1] + g[p][t];
}
//+ sum(i in customeri)w[<p,i>][t]
}
Thank you in advance for your help!
------------------------------
Hieu Truong
------------------------------