Link to Capacitated Vehicle Routing problem CVRP == https://github.com/mayashenoi/DO/blob/main/cvrp-cplex(2).ipynb
# Capcitaed Vehicle Routing problem.
# Minimize the number of trips by a vehicle based on the capacity of the vehicle delivering quantities to customers
import numpy as nprnd = np.random
rnd.seed(0)n = 10 #number of customers/locations
Q = 20 #vehicle capacity
N = [i for i in range(1, n+1)] #customer set
V = [0] + N # set of nodes, customer + 1 depot
q = {i: rnd.randint(1, 10) for i in N} #amount to be delivered to this client random number from 1 to 10loc_x = rnd.rand(len(V))*200
loc_y = rnd.rand(len(V))*100import matplotlib.pyplot as pltplt.scatter(loc_x[1:], loc_y[1:], c='b')
for i in N:
plt.annotate('$q_%d=%d$' % (i, q[i]), (loc_x[i]+2, loc_y[i]))
plt.plot(loc_x[0], loc_y[0], c='r', marker='s')
plt.axis('equal')(-5.390764142541891, 202.16699573081289, 7.525723005074169, 102.16355380509556)
A = [(i, j) for i in V for j in V if i != j] #A is set of arcs
c = {(i, j): np.hypot(loc_x[i]-loc_x[j], loc_y[i]-loc_y[j]) for i, j in A} # cost to travel over the arcsfrom docplex.mp.model import Modelmdl = Model('CVRP')x = mdl.binary_var_dict(A, name='x') # dictionary of arcs
u = mdl.continuous_var_dict(N, ub=Q, name='u') #dictionary of packages/quantities to be delivered at each customermdl.minimize(mdl.sum(c[i, j]*x[i, j] for i, j in A)) # minimize distance travelled
mdl.add_constraints(mdl.sum(x[i, j] for j in V if j != i) == 1 for i in N) # ensure each customer is visited only once
mdl.add_constraints(mdl.sum(x[i, j] for i in V if i != j) == 1 for j in N) # ensure vehicle leaves the customer location
mdl.add_indicator_constraints(mdl.indicator_constraint(x[i, j], u[i]+q[j] == u[j]) for i, j in A if i != 0 and j != 0) #route starts from depot and does not go in subtour(??)
mdl.add_constraints(u[i] >= q[i] for i in N) #truck's capacity is greater than customer's quantity to be delivered
mdl.parameters.timelimit = 15
solution = mdl.solve(log_output=True)WARNING: Number of workers has been reduced to 2 to comply with platform limitations.
Version identifier: 22.1.1.0 | 2023-06-15 | d64d5bd77
CPXPARAM_Read_DataCheck 1
CPXPARAM_Threads 2
CPXPARAM_TimeLimit 15
Found incumbent of value 1784.255531 after 0.00 sec. (0.01 ticks)
Tried aggregator 2 times.
MIP Presolve eliminated 10 rows and 0 columns.
MIP Presolve modified 45 coefficients.
Aggregator did 45 substitutions.
Reduced MIP has 65 rows, 165 columns, and 335 nonzeros.
Reduced MIP has 110 binaries, 0 generals, 0 SOSs, and 90 indicators.
Presolve time = 0.01 sec. (0.35 ticks)
Probing time = 0.00 sec. (0.32 ticks)
Tried aggregator 1 time.
Detecting symmetries...
Reduced MIP has 65 rows, 165 columns, and 335 nonzeros.
Reduced MIP has 110 binaries, 0 generals, 0 SOSs, and 90 indicators.
Presolve time = 0.00 sec. (0.23 ticks)
Probing time = 0.00 sec. (0.31 ticks)
Clique table members: 65.
MIP emphasis: balance optimality and feasibility.
MIP search method: dynamic search.
Parallel mode: deterministic, using up to 2 threads.
Root relaxation solution time = 0.00 sec. (0.13 ticks)
Nodes Cuts/
Node Left Objective IInf Best Integer Best Bound ItCnt Gap
* 0+ 0 1784.2555 0.0000 100.00%
0 0 456.3031 10 1784.2555 456.3031 15 74.43%
0 0 477.8467 8 1784.2555 Cuts: 23 26 73.22%
0 0 477.8467 14 1784.2555 Cuts: 18 38 73.22%
0 0 477.8467 12 1784.2555 Cuts: 3 39 73.22%
* 0+ 0 875.2215 477.8467 45.40%
* 0+ 0 851.1001 477.8467 43.86%
Detecting symmetries...
0 2 477.8467 12 851.1001 477.8467 39 43.86%
Elapsed time = 0.05 sec. (3.90 ticks, tree = 0.02 MB, solutions = 3)
* 63 24 integral 0 792.8484 481.2282 155 39.30%
* 82+ 49 766.8162 481.2282 37.24%
* 103 55 integral 0 732.0147 481.2282 290 34.26%
* 299+ 127 726.2486 483.2994 33.45%
Performing restart 1
Repeating presolve.
Tried aggregator 1 time.
Reduced MIP has 65 rows, 165 columns, and 335 nonzeros.
Reduced MIP has 110 binaries, 0 generals, 0 SOSs, and 90 indicators.
Presolve time = 0.00 sec. (0.16 ticks)
Tried aggregator 1 time.
Reduced MIP has 65 rows, 165 columns, and 335 nonzeros.
Reduced MIP has 110 binaries, 0 generals, 0 SOSs, and 90 indicators.
Presolve time = 0.00 sec. (0.24 ticks)
Represolve time = 0.01 sec. (1.30 ticks)
3702 0 484.0331 19 726.2486 Cuts: 13 20701 16.38%
3702 0 490.4492 13 726.2486 Cuts: 3 20704 16.38%
3702 0 491.3074 7 726.2486 Cuts: 36 20717 16.38%
3702 0 500.2467 17 726.2486 Cuts: 31 20723 16.38%
3702 0 503.0517 10 726.2486 Cuts: 20 20734 16.38%
3702 0 503.0517 15 726.2486 Cuts: 22 20740 16.38%
3702 0 503.0517 12 726.2486 Impl Bds: 3 20742 16.38%
4553 282 cutoff 726.2486 607.2656 23920 16.38%
7445 274 667.3517 15 726.2486 648.1099 46416 10.76%
Clique cuts applied: 1
Cover cuts applied: 66
Implied bound cuts applied: 94
Mixed integer rounding cuts applied: 1
Zero-half cuts applied: 2
Gomory fractional cuts applied: 7
Root node processing (before b&c):
Real time = 0.06 sec. (3.87 ticks)
Parallel b&c, 2 threads:
Real time = 0.92 sec. (522.21 ticks)
Sync time (average) = 0.06 sec.
Wait time (average) = 0.00 sec.
------------
Total (root+branch&cut) = 0.98 sec. (526.08 ticks)
print(solution)solution for: CVRP
objective: 726.249
status: OPTIMAL_SOLUTION(2)
x_0_2=1
x_0_3=1
x_0_10=1
x_1_0=1
x_2_9=1
x_3_5=1
x_4_0=1
x_5_0=1
x_6_1=1
x_7_8=1
x_8_4=1
x_9_7=1
x_10_6=1
u_1=18.000
u_2=2.000
u_3=12.000
u_4=20.000
u_5=20.000
u_6=12.000
u_7=13.000
u_8=16.000
u_9=7.000
u_10=8.000
solution.solve_status<JobSolveStatus.OPTIMAL_SOLUTION: 2>active_arcs = [a for a in A if x[a].solution_value > 0.9]plt.scatter(loc_x[1:], loc_y[1:], c='b')
for i in N:
plt.annotate('$q_%d=%d$' % (i, q[i]), (loc_x[i]+2, loc_y[i]))
for i, j in active_arcs:
plt.plot([loc_x[i], loc_x[j]], [loc_y[i], loc_y[j]], c='g', alpha=0.3)
plt.plot(loc_x[0], loc_y[0], c='r', marker='s')
plt.axis('equal')(-5.390764142541891, 202.16699573081289, 7.525723005074169, 102.16355380509556)
Maya Shenoi's Tech Notes 