Note
Go to the end to download the full example code.
Multicommodity transportation problem#
Please refer to chapters 4 & 20 in the AMPL Book for a detailed description of the problem. We will illustrate how to implement this model with DOcplex using the docplex-extensions library.
Implementation reference: ampl/colab.ampl.com — Copyright (c) 2022-2022, AMPL Optimization inc. (licensed under the MIT License)
Ensure DOcplex and its required dependecies are installed correctly#
from docplex.mp.check_list import run_docplex_check_list
run_docplex_check_list()
* system is: Linux 64bit
* Python version 3.10.17, located at: /home/docs/checkouts/readthedocs.org/user_builds/docplex-extensions/envs/stable/bin/python
* docplex is present, version is 2.29.245
* CPLEX library is present, version is 22.1.2.0, located at: /home/docs/checkouts/readthedocs.org/user_builds/docplex-extensions/envs/stable/lib/python3.10/site-packages
* pandas is present, version is 2.2.3
* Your cplex version 22.1.2.0 is the latest available
! Cplex promotional version, limited to 1000 variables, 1000 constraints
* diagnostics: 2
-- Module cloudpickle is missing, run: pip install cloudpickle
-- Your local CPLEX edition is limited. Consider purchasing a full license.
Import DOcplex, docplex-extensions, and pandas#
import pandas as pd
from docplex.mp.model import Model
import docplex_extensions as dex
Input data#
df_supply = (
pd.DataFrame(
{
'CLEV': {'bands': 700, 'coils': 1600, 'plate': 300},
'GARY': {'bands': 400, 'coils': 800, 'plate': 200},
'PITT': {'bands': 800, 'coils': 1800, 'plate': 300},
}
)
.rename_axis('PROD', axis=0)
.rename_axis('ORI', axis=1)
.transpose()
)
print(df_supply)
PROD bands coils plate
ORI
CLEV 700 1600 300
GARY 400 800 200
PITT 800 1800 300
df_demand = (
pd.DataFrame(
{
'DET': {'bands': 300, 'coils': 750, 'plate': 100},
'FRA': {'bands': 300, 'coils': 500, 'plate': 100},
'FRE': {'bands': 225, 'coils': 850, 'plate': 100},
'LAF': {'bands': 250, 'coils': 500, 'plate': 250},
'LAN': {'bands': 100, 'coils': 400, 'plate': 0},
'STL': {'bands': 650, 'coils': 950, 'plate': 200},
'WIN': {'bands': 75, 'coils': 250, 'plate': 50},
}
)
.rename_axis('PROD', axis=0)
.rename_axis('DES', axis=1)
.transpose()
)
print(df_demand)
PROD bands coils plate
DES
DET 300 750 100
FRA 300 500 100
FRE 225 850 100
LAF 250 500 250
LAN 100 400 0
STL 650 950 200
WIN 75 250 50
df_vcost = (
pd.DataFrame(
{
('CLEV', 'DET'): {'bands': 7, 'coils': 9, 'plate': 9},
('CLEV', 'FRA'): {'bands': 22, 'coils': 27, 'plate': 29},
('CLEV', 'FRE'): {'bands': 82, 'coils': 95, 'plate': 99},
('CLEV', 'LAF'): {'bands': 13, 'coils': 17, 'plate': 18},
('CLEV', 'LAN'): {'bands': 10, 'coils': 12, 'plate': 13},
('CLEV', 'STL'): {'bands': 21, 'coils': 26, 'plate': 28},
('CLEV', 'WIN'): {'bands': 7, 'coils': 9, 'plate': 9},
('GARY', 'DET'): {'bands': 10, 'coils': 14, 'plate': 15},
('GARY', 'FRA'): {'bands': 30, 'coils': 39, 'plate': 41},
('GARY', 'FRE'): {'bands': 71, 'coils': 82, 'plate': 86},
('GARY', 'LAF'): {'bands': 6, 'coils': 8, 'plate': 8},
('GARY', 'LAN'): {'bands': 8, 'coils': 11, 'plate': 12},
('GARY', 'STL'): {'bands': 11, 'coils': 16, 'plate': 17},
('GARY', 'WIN'): {'bands': 10, 'coils': 14, 'plate': 16},
('PITT', 'DET'): {'bands': 11, 'coils': 14, 'plate': 14},
('PITT', 'FRA'): {'bands': 19, 'coils': 24, 'plate': 26},
('PITT', 'FRE'): {'bands': 83, 'coils': 99, 'plate': 104},
('PITT', 'LAF'): {'bands': 15, 'coils': 20, 'plate': 20},
('PITT', 'LAN'): {'bands': 12, 'coils': 17, 'plate': 17},
('PITT', 'STL'): {'bands': 25, 'coils': 28, 'plate': 31},
('PITT', 'WIN'): {'bands': 10, 'coils': 13, 'plate': 13},
}
)
.rename_axis('PROD', axis=0)
.rename_axis(['ORI', 'DES'], axis=1)
.transpose()
)
print(df_vcost)
PROD bands coils plate
ORI DES
CLEV DET 7 9 9
FRA 22 27 29
FRE 82 95 99
LAF 13 17 18
LAN 10 12 13
STL 21 26 28
WIN 7 9 9
GARY DET 10 14 15
FRA 30 39 41
FRE 71 82 86
LAF 6 8 8
LAN 8 11 12
STL 11 16 17
WIN 10 14 16
PITT DET 11 14 14
FRA 19 24 26
FRE 83 99 104
LAF 15 20 20
LAN 12 17 17
STL 25 28 31
WIN 10 13 13
df_fcost = (
pd.DataFrame(
{
'DET': {'CLEV': 1000, 'GARY': 1200, 'PITT': 1200},
'FRA': {'CLEV': 2000, 'GARY': 3000, 'PITT': 2000},
'FRE': {'CLEV': 3000, 'GARY': 3500, 'PITT': 3500},
'LAF': {'CLEV': 2200, 'GARY': 2500, 'PITT': 2200},
'LAN': {'CLEV': 1500, 'GARY': 1200, 'PITT': 1500},
'STL': {'CLEV': 2500, 'GARY': 2500, 'PITT': 2500},
'WIN': {'CLEV': 1200, 'GARY': 1200, 'PITT': 1500},
}
)
.rename_axis('ORI', axis=0)
.rename_axis('DEST', axis=1)
)
print(df_fcost)
DEST DET FRA FRE LAF LAN STL WIN
ORI
CLEV 1000 2000 3000 2200 1500 2500 1200
GARY 1200 3000 3500 2500 1200 2500 1200
PITT 1200 2000 3500 2200 1500 2500 1500
Set up problem data#
Index-sets#
\[\begin{split}\begin{align*}
i &\in ORIG: && \text{Origins to ship products from} \\
j &\in DEST: && \text{Destinations to ship products to} \\
p &\in PROD: && \text{Products to be shipped} \\
\end{align*}\end{split}\]
ORIG = df_supply.index.dex.to_indexset()
DEST = df_demand.index.dex.to_indexset()
PROD = df_supply.columns.dex.to_indexset()
Parameters#
\[\begin{split}\begin{align*}
supply_{i, p} &\in \mathbb{R}_{0}^{+}: && \text{Units of prodcut $p$ available at origin $i$},
\; \forall \; i \in ORIG \; \text{and} \; p \in PROD \\
demand_{j, p} &\in \mathbb{R}_{0}^{+}: && \text{Units of prodcut $p$ required at destination
$j$}, \; \forall \; j \in DEST \; \text{and} \; p \in PROD \\
vcost_{i, j, p} &\in \mathbb{R}_{0}^{+}: && \text{Variable cost of shipping product $p$ from
origin $i$ to destination $j$}, \; \forall \; i \in ORIG \; \text{and} \; j \in DEST \;
\text{and} \; p \in PROD \\
fcost_{i, j} &\in \mathbb{R}_{0}^{+}: && \text{Fixed cost of shipping from origin $i$ to
destination $j$}, \; \forall \; i \in ORIG \; \text{and} \; j \in DEST \\
limit &\in \mathbb{R}_{0}^{+}: && \text{Limit on shipping total units of products from any
origin to any destination} \\
\end{align*}\end{split}\]
supply = df_supply.stack().dex.to_paramdict()
demand = df_demand.stack().dex.to_paramdict()
vcost = df_vcost.stack().dex.to_paramdict()
fcost = df_fcost.stack().dex.to_paramdict()
limit = 625
Set up model#
Instantiate#
model = Model(name='multicommodity')
Add decision variables#
\[\begin{split}\begin{align*}
trans_{i, j, p} &\in \mathbb{R}_{0}^{+}: && \text{Units of prodcut $p$ to ship from origin $i$
to destination $j$}, \; \forall \; i \in ORIG \; \text{and} \; j \in DEST \; \text{and} \;
p \in PROD \\
use_{i, j} &\in \mathbb{B}: && \text{Whether to ship from origin $i$ to destination $j$}, \;
\forall \; i \in ORIG \; \text{and} \; j \in DEST \\
\end{align*}\end{split}\]
trans = dex.add_variables(
model,
indexset=dex.IndexSetND(ORIG, DEST, PROD, names=('ORIG', 'DEST', 'PROD')),
var_type='continuous',
name='NUM-UNITS',
)
use = dex.add_variables(
model,
indexset=dex.IndexSetND(ORIG, DEST, names=('ORIG', 'DEST')),
var_type='binary',
name='USE-ROUTE',
)
Set objective#
\[\begin{split}\begin{align*}
&\text{Minimize the total shipping cost:} \\
&\text{min.} \; \sum_{i \in ORIG} \sum_{j \in DEST} \sum_{p \in PROD} vcost_{i, j, p} \times
trans_{i, j, p} + \sum_{i \in ORIG} \sum_{j \in DEST} fcost_{i, j} \times use_{i, j} \\
\end{align*}\end{split}\]
# fmt: off
model.minimize(
model.sum(
vcost[i, j, p] * trans[i, j, p]
for i in ORIG
for j in DEST
for p in PROD
)
+ model.sum(
fcost[i, j] * use[i, j]
for i in ORIG
for j in DEST
)
)
# fmt: on
Add constraints#
\[\begin{split}\begin{align*}
&\text{Total units of prodcut $p$ shipped from origin $i$ should be equal to its supply:} \\
& \sum_{j \in DEST} trans_{i, j, p} = supply_{i, p}, \; \forall \; i \in ORIG \; \text{and} \;
p \in PROD \\
&\text{Total units of prodcut $p$ shipped to destination $j$ should be equal to its demand:} \\
& \sum_{i \in ORIG} trans_{i, j, p} = demand_{j, p}, \; \forall \; j \in DEST \; \text{and} \;
p \in PROD \\
&\text{Total units of prodcuts shipped from origin $i$ to destination $j$ is limited:} \\
& \sum_{p \in PROD} trans_{i, j, p} \leq limit, \; \forall \; i \in ORIG \; \text{and} \; j
\in DEST \\
\end{align*}\end{split}\]
# fmt: off
supply_cons = model.add_constraints(
(trans.sum(i, '*', p) == supply[i, p], f'supply_{i}_{p}')
for i in ORIG
for p in PROD
)
demand_cons = model.add_constraints(
(trans.sum('*', j, p) == demand[j, p], f'demand_{j}_{p}')
for j in DEST
for p in PROD
)
multi_cons = model.add_constraints(
(trans.sum(i, j, '*') <= limit * use[i,j], f'multi_{i}_{j}')
for i in ORIG
for j in DEST
)
# fmt: on
Solve#
sol = dex.solve(model, log_output=True)
----------------------------- MILP problem statistics -----------------------------
Problem name : multicommodity
Objective sense : Minimize
Variables : 84 [Nneg: 63, Binary: 21]
Objective nonzeros : 84
Linear constraints : 51 [Less: 21, Equal: 30]
Nonzeros : 210
RHS nonzeros : 29
Variables : Min LB: 0.000000 Max UB: 1.000000
Objective nonzeros : Min : 6.000000 Max : 3500.000
Linear constraints :
Nonzeros : Min : 1.000000 Max : 625.0000
RHS nonzeros : Min : 50.00000 Max : 1800.000
------------------------------- CPLEX optimizer log -------------------------------
Version identifier: 22.1.2.0 | 2024-12-10 | f4cec290b
CPXPARAM_Read_DataCheck 1
Tried aggregator 1 time.
MIP Presolve eliminated 1 rows and 3 columns.
MIP Presolve modified 6 coefficients.
Reduced MIP has 50 rows, 81 columns, and 201 nonzeros.
Reduced MIP has 21 binaries, 0 generals, 0 SOSs, and 0 indicators.
Presolve time = 0.00 sec. (0.10 ticks)
Found incumbent of value 249325.000000 after 0.00 sec. (0.31 ticks)
Probing fixed 3 vars, tightened 0 bounds.
Probing time = 0.00 sec. (0.02 ticks)
Tried aggregator 1 time.
Detecting symmetries...
MIP Presolve eliminated 0 rows and 3 columns.
Reduced MIP has 50 rows, 78 columns, and 198 nonzeros.
Reduced MIP has 18 binaries, 0 generals, 0 SOSs, and 0 indicators.
Presolve time = 0.00 sec. (0.11 ticks)
Probing time = 0.00 sec. (0.01 ticks)
Clique table members: 3.
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.16 ticks)
Nodes Cuts/
Node Left Objective IInf Best Integer Best Bound ItCnt Gap
* 0+ 0 249325.0000 7500.0000 96.99%
* 0+ 0 235800.0000 7500.0000 96.82%
0 0 224640.0000 10 235800.0000 224640.0000 31 4.73%
0 0 228881.6667 7 235800.0000 Cuts: 58 57 2.93%
0 0 229455.0909 7 235800.0000 Cuts: 21 71 2.69%
* 0+ 0 229850.0000 229455.0909 0.17%
0 0 229574.1667 4 229850.0000 Cuts: 19 76 0.12%
0 0 229592.3329 11 229850.0000 Cuts: 5 90 0.11%
0 0 229593.3428 11 229850.0000 Cuts: 4 93 0.11%
Detecting symmetries...
0 0 229593.5242 11 229850.0000 Flowcuts: 2 95 0.11%
0 0 cutoff 229850.0000 229850.0000 95 0.00%
Elapsed time = 0.01 sec. (6.96 ticks, tree = 0.01 MB, solutions = 3)
Implied bound cuts applied: 10
Flow cuts applied: 6
Mixed integer rounding cuts applied: 9
Multi commodity flow cuts applied: 5
Gomory fractional cuts applied: 5
Root node processing (before b&c):
Real time = 0.01 sec. (6.97 ticks)
Parallel b&c, 2 threads:
Real time = 0.00 sec. (0.00 ticks)
Sync time (average) = 0.00 sec.
Wait time (average) = 0.00 sec.
------------
Total (root+branch&cut) = 0.01 sec. (6.97 ticks)
--------------------------- Solution quality statistics ---------------------------
Incumbent solution:
MILP objective 2.2985000000e+05
MILP solution norm |x| (Total, Max) 6.91400e+03 6.25000e+02
MILP solution error (Ax=b) (Total, Max) 2.98428e-13 1.13687e-13
MILP x bound error (Total, Max) 6.25278e-13 2.27374e-13
MILP x integrality error (Total, Max) 0.00000e+00 0.00000e+00
MILP slack bound error (Total, Max) 1.13687e-12 3.97904e-13
-------------------------------------------------------------------------------------
Solution#
sol.display(print_zeros=False)
solution for: multicommodity
objective: 229850.000
status: OPTIMAL_SOLUTION(2)
NUM-UNITS_CLEV_DET_coils = 525.000
NUM-UNITS_CLEV_DET_plate = 100.000
NUM-UNITS_CLEV_FRA_bands = 225.000
NUM-UNITS_CLEV_FRA_plate = 50.000
NUM-UNITS_CLEV_FRE_coils = 0.000
NUM-UNITS_CLEV_LAF_bands = 50.000
NUM-UNITS_CLEV_LAF_coils = 375.000
NUM-UNITS_CLEV_LAN_coils = 350.000
NUM-UNITS_CLEV_STL_bands = 350.000
NUM-UNITS_CLEV_STL_coils = 100.000
NUM-UNITS_CLEV_STL_plate = 100.000
NUM-UNITS_CLEV_WIN_bands = 75.000
NUM-UNITS_CLEV_WIN_coils = 250.000
NUM-UNITS_CLEV_WIN_plate = 50.000
NUM-UNITS_GARY_FRE_coils = 525.000
NUM-UNITS_GARY_FRE_plate = 100.000
NUM-UNITS_GARY_LAF_coils = -0.000
NUM-UNITS_GARY_LAN_bands = 100.000
NUM-UNITS_GARY_LAN_coils = 50.000
NUM-UNITS_GARY_STL_bands = 300.000
NUM-UNITS_GARY_STL_coils = 225.000
NUM-UNITS_GARY_STL_plate = 100.000
NUM-UNITS_PITT_DET_bands = 300.000
NUM-UNITS_PITT_DET_coils = 225.000
NUM-UNITS_PITT_FRA_bands = 75.000
NUM-UNITS_PITT_FRA_coils = 500.000
NUM-UNITS_PITT_FRA_plate = 50.000
NUM-UNITS_PITT_FRE_bands = 225.000
NUM-UNITS_PITT_FRE_coils = 325.000
NUM-UNITS_PITT_FRE_plate = -0.000
NUM-UNITS_PITT_LAF_bands = 200.000
NUM-UNITS_PITT_LAF_coils = 125.000
NUM-UNITS_PITT_LAF_plate = 250.000
NUM-UNITS_PITT_STL_coils = 625.000
NUM-UNITS_PITT_STL_plate = -0.000
NUM-UNITS_PITT_WIN_coils = -0.000
USE-ROUTE_CLEV_DET = 1
USE-ROUTE_CLEV_FRA = 1
USE-ROUTE_CLEV_LAF = 1
USE-ROUTE_CLEV_LAN = 1
USE-ROUTE_CLEV_STL = 1
USE-ROUTE_CLEV_WIN = 1
USE-ROUTE_GARY_FRE = 1
USE-ROUTE_GARY_LAN = 1
USE-ROUTE_GARY_STL = 1
USE-ROUTE_PITT_DET = 1
USE-ROUTE_PITT_FRA = 1
USE-ROUTE_PITT_FRE = 1
USE-ROUTE_PITT_LAF = 1
USE-ROUTE_PITT_STL = 1
Total running time of the script: (0 minutes 0.044 seconds)