// Data section
int n = ...; // number of nodes from node.dat
range N = 1..n;
range C = N; // set of customers (non-hub nodes)
range H = N; // set of potential hubs
// Data reading from external files
tuple Coordinate {
float x;
float y;
}
Coordinate coordinates[N] = ...;
float d[N][N] = ...; // demand matrix
float F[N] = ...; // fixed costs for hubs
float Qh[N] = ...; // hub capacities
int K = ...; // number of vehicles
float t[N][N] = ...; // travel times
tuple TimeWindow {
float start;
float end;
}
TimeWindow timewindows[N] = ...;
// Constants
float Qv = 1100; // vehicle capacity
float M = 1000000; // big M
float alpha = 0.004; // unit inter-hub transportation cost
float beta = 1.0; // unit local delivery cost
float zeta = 1.0; // early arrival penalty coefficient
float gamma = 2.0; // late arrival penalty coefficient
// Calculate distances between nodes
float c[i in N][j in N] = sqrt(pow(coordinates[i].x - coordinates[j].x, 2) +
pow(coordinates[i].y - coordinates[j].y, 2));
// Decision Variables
dvar boolean x[N][N][H]; // routing variable
dvar boolean y[H]; // hub location variable
dvar float+ z[C][C][H][H]; // flow variable
dvar float+ u[N][H]; // pickup amount
dvar float+ v[N][H]; // delivery amount
dvar float+ s[N]; // service start time
dvar float+ P[N]; // time window penalty
dvar float+ visitOrder[i in N]; // exclude subtour
// Rest of the model remains the same...
// [Previous constraints and objective function code goes here]
// Optimization Model
minimize
beta * sum(i in N, j in N, m in H) c[i][j] * x[i][j][m] +
sum(m in H) F[m] * y[m] +
alpha * sum(m in H, n in H, i in C, j in C) d[i][j] * z[i][j][m][n] * c[m][n] +
sum(i in N) P[i];
subject to {
//self-looping protection
forall(i in N, m in H)
ct0: x[i][i][m] == 0;
//Return to the Hub
forall(m in H)
ct1: sum(i in N: i != m) x[i][m][m] == sum(j in N: j != m) x[m][j][m];
// Constraint (2): All demand must be routed through hubs
forall(i in C, j in C)
ct2:sum(m in H, n in H) z[i][j][m][n] == 1;
// Constraint (3): Hub capacity constraints
forall(m in H)
ct3:sum(i in C, j in C, n in H) d[i][j] * z[i][j][m][n] +
sum(i in C, j in C, n in H: n != m) d[i][j] * z[i][j][n][m] <= Qh[m] * y[m];
// Constraint (4): Vehicles can only depart from established hubs
forall(m in H)
ct4:sum(j in N) x[m][j][m] <= M * y[m];
// Constraint (5): Each customer must be visited exactly once
forall(i in C)
ct5:sum(j in N: j != i, m in H) x[i][j][m] == 1;
// Constraint (6): Flow conservation
forall(i in N, m in H)
ct6:sum(j in N) x[i][j][m] == sum(j in N) x[j][i][m];
// Constraint (7): No inter-hub routes in local delivery
forall(m in H)
ct7:sum(j in N, n in H: n != m) x[n][j][m] == 0;
// Constraint (8): Link between flow and routing variables
forall(i in C, m in H)
ct8:M * sum(j in N) x[i][j][m] >=
sum(j in C, n in H) z[i][j][m][n] + sum(j in C, n in H) z[j][i][n][m];
// Constraint (9): Pickup load continuity
forall(i in N, j in C, m in H: i != j)
ct9:u[i][m] + sum(t in C, n in H) d[j][t] * z[j][t][m][n] -
M * (1 - x[i][j][m]) <= u[j][m];
// Constraint (10): Delivery load continuity
forall(i in C, j in N, m in H: i != j)
ct10:v[i][m] - sum(t in C, n in H) d[t][i] * z[t][i][n][m] +
M * (1 - x[i][j][m]) >= v[j][m];
// Constraint (11): Vehicle capacity constraint for delivery
forall(i in C, m in H)
ct11:v[i][m] <= Qv;
// Constraint (12): Vehicle capacity constraint considering loading/unloading
forall(i in C, m in H)
ct12:u[i][m] + v[i][m] - sum(t in C, n in H) d[t][i] * z[t][i][n][m] <= Qv;
// Constraint (13): Number of vehicles limit
ct13:sum(m in H) sum(j in N) x[m][j][m] <= K;
// Constraint (14): Time consistency
forall(i in N, j in N, m in H)
ct14:s[i] + t[i][j] - s[j] <= M * (1 - x[i][j][m]);
// Constraint (15): Time window penalties
forall(i in C) {
ct15:P[i] >= zeta * (timewindows[i].start - s[i]);
P[i] >= gamma * (s[i] - timewindows[i].end);
P[i] >= 0;};
//go around in turn
forall(i in N, j in N, m in H: i != j && i != m && j != m)
ct16: visitOrder[i] - visitOrder[j] + n * x[i][j][m] <= n - 1;
};
// Solver configuration
execute {
cp.tilim = 10800; // Time limit: 10800 seconds
}
↓yyy.dat~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
n = 4;
coordinates = [
<12944.330389, 19522.690462>,
<23487.769950, 14749.226168>,
<31004.306839, 21458.349354>,
<49728.492698, 17946.318136>
];
d = [
[5.705460, 9.801600, 5.028110, 3.105340],
[21.242070, 38.706030, 18.462480, 10.396320],
[5.763190, 9.909430, 9.328980, 8.665890],
[3.820190, 6.300870, 9.807810, 10.572880]
];
F = [
28766.736921,
28376.761527,
29774.238965,
24301.334212
];
Qh = [
2198.532635,
6106.711881,
2226.593104,
4515.466789
];
K = 4;
timewindows = [
<0,10000>,
<0,10000>,
<0,10000>,
<0,10000>
];
t = [
[0.00 233.68 19.41 529.92],
[233.68 0.00 436.25 352.77],
[19.41 436.25 0.00 339.71],
[529.92 352.77 339.71 0.00]
];