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 demonstrate how to implement this model with DOcplex, gurobipy, Xpress, and highspy along with the functionality from opti-extensions.
Implementation reference: ampl/colab.ampl.com — Copyright (c) 2022-2022, AMPL Optimization inc. (licensed under the MIT License)
Mathematical formulation#
Index-sets:
- Origins to ship products from:\(i \in ORIG\)
- Destinations to ship products to:\(j \in DEST\)
- Products to be shipped:\(p \in PROD\)
Parameters:
- Units of product \(p\) available at origin \(i\):\(supply_{i, p} \in \mathbb{R}_{0}^{+} \quad \forall \; i \in ORIG \; \& \; p \in PROD\)
- Units of product \(p\) required at destination \(j\):\(demand_{j, p} \in \mathbb{R}_{0}^{+} \quad \forall \; j \in DEST \; \& \; p \in PROD\)
- Variable cost of shipping product \(p\) from origin \(i\) to destination \(j\):\(vcost_{i, j, p} \in \mathbb{R}_{0}^{+} \quad \forall \; i \in ORIG \; \& \; j \in DEST \; \& \; p \in PROD\)
- Fixed cost of shipping from origin \(i\) to destination \(j\):\(fcost_{i, j} \in \mathbb{R}_{0}^{+} \quad \forall \; i \in ORIG \; \& \; j \in DEST\)
- Limit on shipping total units of products from any origin to any destination:\(limit \in \mathbb{R}_{0}^{+}\)
Variables:
- Units of product \(p\) to ship from origin \(i\) to destination \(j\):\(trans_{i, j, p} \in \mathbb{R}_{0}^{+} \quad \forall \; i \in ORIG \; \& \; j \in DEST \; \& \; p \in PROD\)
- Whether to ship from origin \(i\) to destination \(j\):\(use_{i, j} \in \mathbb{B} \quad \forall \; i \in ORIG \; \& \; j \in DEST\)
Objective:
- Minimize the total shipping cost:\(\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}\)
Constraints:
- Total units of product \(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 \; \& \; p \in PROD\)
- Total units of product \(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 \; \& \; p \in PROD\)
- Total units of products shipped from origin \(i\) to destination \(j\) is limited:\(\sum_{p \in PROD} trans_{i, j, p} \leq limit, \; \forall \; i \in ORIG \; \& \; j \in DEST\)
Input data in form of pandas dataframes#
import pandas as pd
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)
# fmt: off
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#
# Need to import the package, in some form, to directly cast pandas objects to our data structures
Index-sets#
ORIG = df_supply.index.opti.to_indexset()
DEST = df_demand.index.opti.to_indexset()
PROD = df_supply.columns.opti.to_indexset()
Parameters#
supply = df_supply.stack().opti.to_paramdict()
demand = df_demand.stack().opti.to_paramdict()
vcost = df_vcost.stack().opti.to_paramdict()
fcost = df_fcost.stack().opti.to_paramdict()
limit = 625
Implement with DOcplex#
from docplex.mp.model import Model
from opti_extensions import IndexSetND
from opti_extensions.docplex import add_variables, solve
# Instantiate model
model = Model(name='multicommodity')
# Add variables
trans = add_variables(
model,
indexset=IndexSetND(ORIG, DEST, PROD, names=('ORIG', 'DEST', 'PROD')),
vartype='continuous',
name='NUM-UNITS',
)
# Instead of:
# trans = model.continuous_var_dict(
# [(i, j, p) for i in ORIG for j in DEST for p in PROD],
# name='NUM-UNITS',
# )
use = add_variables(
model,
indexset=IndexSetND(ORIG, DEST, names=('ORIG', 'DEST')),
vartype='binary',
name='USE-ROUTE',
)
# Instead of:
# use = model.binary_var_dict(
# [(i, j) for i in ORIG for j in DEST],
# name='USE-ROUTE',
# )
# Set objective
model.minimize(
vcost @ trans + fcost @ use
# Instead of:
# 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
# )
)
# Add constraints
model.add_constraints_(
trans.sum(i, '*', p) == supply[i, p]
# Instead of:
# model.sum(trans[i, j, p] for j in DEST) == supply[i, p]
for i in ORIG for p in PROD
)
model.add_constraints_(
trans.sum('*', j, p) == demand[j, p]
# Instead of:
# model.sum(trans[i, j, p] for i in ORIG) == demand[j, p]
for j in DEST for p in PROD
)
model.add_constraints_(
trans.sum(i, j, '*') <= limit * use[i, j]
# Instead of:
# model.sum(trans[i, j, p] for p in PROD) <= limit * use[i,j]
for i in ORIG for j in DEST
)
# Solve with additional logging output for problem and solution statistics
sol = solve(model, log_output=True)
# Instead of:
# sol = model.solve(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 | 2026-03-02 | af0ce9b93
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.95 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.95 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.95 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
-------------------------------------------------------------------------------------
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
Implement with gurobipy#
from gurobipy import GRB, Model
from opti_extensions.gurobipy import addVars
# Instantiate model
model = Model(name='multicommodity')
# Add variables
trans = addVars(
model,
IndexSetND(ORIG, DEST, PROD, names=('ORIG', 'DEST', 'PROD')),
vtype=GRB.CONTINUOUS,
name='NUM-UNITS',
)
# Instead of:
# trans = model.addVars(
# ORIG,
# DEST,
# PROD,
# vtype=GRB.CONTINUOUS,
# name='NUM-UNITS',
# )
use = addVars(
model,
IndexSetND(ORIG, DEST, names=('ORIG', 'DEST')),
vtype=GRB.BINARY,
name='USE-ROUTE',
)
# Instead of:
# use = model.addVars(
# ORIG,
# DEST,
# vtype=GRB.BINARY,
# name='USE-ROUTE',
# )
# Set objective
model.setObjective(
vcost @ trans + fcost @ use,
# Instead of:
# quicksum(
# vcost[i, j, p] * trans[i, j, p]
# for i in ORIG for j in DEST for p in PROD
# )
# + quicksum(
# fcost[i, j] * use[i, j]
# for i in ORIG for j in DEST
# )
sense=GRB.MINIMIZE,
)
# Add constraints
model.addConstrs(
trans.sum(i, '*', p) == supply[i, p]
# Instead of:
# quicksum(trans[i, j, p] for j in DEST) == supply[i, p]
for i in ORIG for p in PROD
)
model.addConstrs(
trans.sum('*', j, p) == demand[j, p]
# Instead of:
# quicksum(trans[i, j, p] for i in ORIG) == demand[j, p]
for j in DEST for p in PROD
)
model.addConstrs(
trans.sum(i, j, '*') <= limit * use[i, j]
# Instead of:
# quicksum(trans[i, j, p] for p in PROD) <= limit * use[i,j]
for i in ORIG for j in DEST
)
# Solve
model.optimize()
Gurobi Optimizer version 13.0.1 build v13.0.1rc0 (linux64 - "Ubuntu 22.04.3 LTS")
CPU model: AMD EPYC 7R13 Processor, instruction set [SSE2|AVX|AVX2]
Thread count: 1 physical cores, 2 logical processors, using up to 2 threads
Optimize a model with 51 rows, 84 columns and 210 nonzeros (Min)
Model fingerprint: 0x76f09fde
Model has 84 linear objective coefficients
Variable types: 63 continuous, 21 integer (21 binary)
Coefficient statistics:
Matrix range [1e+00, 6e+02]
Objective range [6e+00, 4e+03]
Bounds range [1e+00, 1e+00]
RHS range [5e+01, 2e+03]
Found heuristic solution: objective 249325.00000
Presolve removed 1 rows and 6 columns
Presolve time: 0.00s
Presolved: 50 rows, 78 columns, 198 nonzeros
Variable types: 60 continuous, 18 integer (18 binary)
Root relaxation: objective 2.246400e+05, 39 iterations, 0.00 seconds (0.00 work units)
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
0 0 224640.000 0 8 249325.000 224640.000 9.90% - 0s
H 0 0 234300.00000 224640.000 4.12% - 0s
H 0 0 233550.00000 224640.000 3.82% - 0s
0 0 226484.385 0 11 233550.000 226484.385 3.03% - 0s
0 0 226641.538 0 11 233550.000 226641.538 2.96% - 0s
0 0 226641.538 0 10 233550.000 226641.538 2.96% - 0s
0 0 226641.538 0 10 233550.000 226641.538 2.96% - 0s
0 0 227725.444 0 6 233550.000 227725.444 2.49% - 0s
H 0 0 233285.00000 227891.444 2.31% - 0s
0 0 227891.444 0 6 233285.000 227891.444 2.31% - 0s
0 0 228183.000 0 5 233285.000 228183.000 2.19% - 0s
H 0 0 233207.00000 228183.000 2.15% - 0s
0 0 228183.000 0 6 233207.000 228183.000 2.15% - 0s
0 0 228199.605 0 6 233207.000 228199.605 2.15% - 0s
H 0 0 233188.44186 228199.605 2.14% - 0s
0 0 228199.605 0 6 233188.442 228199.605 2.14% - 0s
0 0 228199.605 0 6 233188.442 228199.605 2.14% - 0s
H 0 0 230206.25000 228199.605 0.87% - 0s
H 0 0 230166.66667 228199.605 0.85% - 0s
0 2 228231.000 0 6 230166.667 228231.000 0.84% - 0s
* 6 3 5 230050.00000 228506.000 0.67% 7.5 0s
H 8 1 229850.00000 228506.000 0.58% 6.1 0s
Cutting planes:
Gomory: 2
Cover: 4
Implied bound: 8
MIR: 7
Flow cover: 6
Explored 15 nodes (204 simplex iterations) in 0.02 seconds (0.01 work units)
Thread count was 2 (of 2 available processors)
Solution count 6: 229850 230050 233188 ... 249325
Optimal solution found (tolerance 1.00e-04)
Best objective 2.298500000000e+05, best bound 2.298500000000e+05, gap 0.0000%
model.printAttr('X')
Variable X
-------------------------
NUM-UNITS[CLEV,DET,coils] 525
NUM-UNITS[CLEV,DET,plate] 100
NUM-UNITS[CLEV,FRA,bands] 275
NUM-UNITS[CLEV,LAF,coils] 375
NUM-UNITS[CLEV,LAF,plate] 50
NUM-UNITS[CLEV,LAN,coils] 350
NUM-UNITS[CLEV,STL,bands] 350
NUM-UNITS[CLEV,STL,coils] 100
NUM-UNITS[CLEV,STL,plate] 100
NUM-UNITS[CLEV,WIN,bands] 75
NUM-UNITS[CLEV,WIN,coils] 250
NUM-UNITS[CLEV,WIN,plate] 50
NUM-UNITS[GARY,FRE,coils] 525
NUM-UNITS[GARY,FRE,plate] 100
NUM-UNITS[GARY,LAN,bands] 100
NUM-UNITS[GARY,LAN,coils] 50
NUM-UNITS[GARY,STL,bands] 300
NUM-UNITS[GARY,STL,coils] 225
NUM-UNITS[GARY,STL,plate] 100
NUM-UNITS[PITT,DET,bands] 300
NUM-UNITS[PITT,DET,coils] 225
NUM-UNITS[PITT,FRA,bands] 25
NUM-UNITS[PITT,FRA,coils] 500
NUM-UNITS[PITT,FRA,plate] 100
NUM-UNITS[PITT,FRE,bands] 225
NUM-UNITS[PITT,FRE,coils] 325
NUM-UNITS[PITT,LAF,bands] 250
NUM-UNITS[PITT,LAF,coils] 125
NUM-UNITS[PITT,LAF,plate] 200
NUM-UNITS[PITT,STL,coils] 625
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
Implement with Xpress#
import xpress as xp
from opti_extensions.xpress import addVariables
# Instantiate problem
prob = xp.problem(name='multicommodity')
# Add variables
trans = addVariables(
prob,
IndexSetND(ORIG, DEST, PROD, names=('ORIG', 'DEST', 'PROD')),
vartype=xp.continuous,
name='NUM-UNITS',
)
# Instead of:
# trans = prob.addVariables(
# [(i, j, p) for i in ORIG for j in DEST for p in PROD],
# vartype=xp.continuous,
# name='NUM-UNITS',
# )
use = addVariables(
prob,
IndexSetND(ORIG, DEST, names=('ORIG', 'DEST')),
vartype=xp.binary,
name='USE-ROUTE',
)
# Instead of:
# use = prob.addVariables(
# [(i, j) for i in ORIG for j in DEST],
# vartype=xp.binary,
# name='USE-ROUTE',
# )
# Set objective
prob.setObjective(
vcost @ trans + fcost @ use,
# Instead of:
# xp.Sum(
# vcost[i, j, p] * trans[i, j, p]
# for i in ORIG for j in DEST for p in PROD
# )
# + xp.Sum(
# fcost[i, j] * use[i, j]
# for i in ORIG for j in DEST
# )
sense=xp.minimize,
)
# Add constraints
prob.addConstraint(
trans.sum(i, '*', p) == supply[i, p]
# Instead of:
# xp.Sum(trans[i, j, p] for j in DEST) == supply[i, p]
for i in ORIG for p in PROD
)
prob.addConstraint(
trans.sum('*', j, p) == demand[j, p]
# Instead of:
# xp.Sum(trans[i, j, p] for i in ORIG) == demand[j, p]
for j in DEST for p in PROD
)
prob.addConstraint(
trans.sum(i, j, '*') <= limit * use[i, j]
# Instead of:
# xp.Sum(trans[i, j, p] for p in PROD) <= limit * use[i,j]
for i in ORIG for j in DEST
)
# Solve
prob.optimize()
FICO Xpress v9.8.1, Community, solve started 19:57:17, Apr 11, 2026
Heap usage: 465KB (peak 465KB, 78KB system)
Minimizing MILP multicommodity using up to 2 threads and up to 7734MB memory, with these control settings:
OUTPUTLOG = 1
NLPPOSTSOLVE = 1
XSLP_DELETIONCONTROL = 0
XSLP_OBJSENSE = 1
Original problem has:
51 rows 84 cols 210 elements 21 entities
Presolved problem has:
50 rows 78 cols 198 elements 18 entities
LP relaxation tightened
Presolve finished in 0 seconds
Heap usage: 1641KB (peak 1641KB, 78KB system)
Coefficient range original solved
Coefficients [min,max] : [ 1.00e+00, 6.25e+02] / [ 1.95e-03, 1.95e+00]
RHS and bounds [min,max] : [ 1.00e+00, 1.80e+03] / [ 1.00e+00, 1.80e+03]
Objective [min,max] : [ 6.00e+00, 3.50e+03] / [ 6.00e+00, 3.50e+03]
Autoscaling applied standard scaling
Will try to keep branch and bound tree memory usage below 6.7GB
Starting concurrent solve with dual (1 thread)
Concurrent-Solve, 0s
Dual
objective dual inf
D 224640.00 .0000000
------- optimal --------
Concurrent statistics:
Dual: 73 simplex iterations, 0.00s
Optimal solution found
Its Obj Value S Ninf Nneg Sum Dual Inf Time
73 224640.0000 D 0 0 .000000 0
Dual solved problem
73 simplex iterations in 0.00 seconds at time 0
Final objective : 2.246400000000000e+05
Max primal violation (abs/rel) : 0.0 / 0.0
Max dual violation (abs/rel) : 4.441e-16 / 4.441e-16
Max complementarity viol. (abs/rel) : 0.0 / 0.0
Starting root cutting & heuristics
Its Type BestSoln BestBound Sols Add Del Gap GInf Time
k 235800.0000 224640.0000 1 4.73% 0 0
1 K 235800.0000 228137.5000 1 50 0 3.25% 10 0
2 K 235800.0000 228207.3780 1 28 37 3.22% 8 0
3 K 235800.0000 228236.0000 1 20 27 3.21% 9 0
4 K 235800.0000 228238.8235 1 10 18 3.21% 8 0
5 K 235800.0000 228260.0000 1 13 9 3.20% 8 0
6 K 235800.0000 228260.0000 1 7 13 3.20% 8 0
7 K 235800.0000 228260.0000 1 29 8 3.20% 8 0
8 K 235800.0000 228260.0000 1 10 29 3.20% 8 0
9 K 235800.0000 228265.5556 1 27 10 3.20% 8 0
10 K 235800.0000 228290.0000 1 10 28 3.18% 7 0
11 K 235800.0000 228297.5000 1 15 9 3.18% 8 0
12 K 235800.0000 228301.3693 1 14 15 3.18% 7 0
13 K 235800.0000 228302.9775 1 20 12 3.18% 10 0
14 K 235800.0000 228305.0000 1 5 20 3.18% 8 0
15 K 235800.0000 228308.0469 1 5 6 3.18% 7 0
16 K 235800.0000 228308.0469 1 0 5 3.18% 7 0
17 G 235800.0000 228618.8866 1 7 0 3.05% 7 0
18 G 235800.0000 228720.6522 1 55 50 3.00% 9 0
d 230550.0000 228720.6522 2 0.79% 0 0
Heuristic search 'R' started
Heuristic search 'R' stopped
M 230050.0000 228720.6522 3 0.58% 0 0
Cuts in the matrix : 29
Cut elements in the matrix : 208
Starting tree search.
Deterministic mode with up to 2 running threads and up to 4 tasks.
Heap usage: 4659KB (peak 7042KB, 101KB system)
Node BestSoln BestBound Sols Active Depth Gap GInf Time
1 230050.0000 228920.4789 3 2 1 0.49% 10 0
2 230050.0000 228920.4789 3 2 3 0.49% 3 0
3 230050.0000 228920.4789 3 2 4 0.49% 2 0
4 230050.0000 228920.4789 3 2 4 0.49% 2 0
5 230050.0000 229281.0000 3 3 3 0.33% 3 0
8 230050.0000 229395.4444 3 2 5 0.28% 1 0
* 8 229850.0000 229445.6250 4 1 2 0.18% 0 0
9 229850.0000 229445.6250 4 1 2 0.18% 1 0
*** Search completed ***
Uncrunching matrix
Final MIP objective : 2.298500000000000e+05
Final MIP bound : 2.298497701500000e+05
Solution time / primaldual integral : 0.05s/ 7.221455%
Work / work units per second : 0.06 / 1.18
Number of solutions found / nodes : 4 / 9
Max primal violation (abs/rel) : 4.050e-13 / 4.050e-13
Max integer violation (abs ) : 8.882e-16
(<SolveStatus.COMPLETED: 3>, <SolStatus.OPTIMAL: 1>)
print(f'var: {"value" :>20}')
print('-' * 25)
for vd in (trans, use):
for k, v in prob.getSolution(vd).items():
if abs(v) > 1e-6:
name = f'{vd.value_name}[{k}]:'
print(f'{name:<15} {v :>8.4f}')
var: value
-------------------------
NUM-UNITS[('CLEV', 'DET', 'coils')]: 525.0000
NUM-UNITS[('CLEV', 'DET', 'plate')]: 100.0000
NUM-UNITS[('CLEV', 'FRA', 'bands')]: 275.0000
NUM-UNITS[('CLEV', 'LAF', 'coils')]: 375.0000
NUM-UNITS[('CLEV', 'LAF', 'plate')]: 50.0000
NUM-UNITS[('CLEV', 'LAN', 'coils')]: 350.0000
NUM-UNITS[('CLEV', 'STL', 'bands')]: 350.0000
NUM-UNITS[('CLEV', 'STL', 'coils')]: 100.0000
NUM-UNITS[('CLEV', 'STL', 'plate')]: 100.0000
NUM-UNITS[('CLEV', 'WIN', 'bands')]: 75.0000
NUM-UNITS[('CLEV', 'WIN', 'coils')]: 250.0000
NUM-UNITS[('CLEV', 'WIN', 'plate')]: 50.0000
NUM-UNITS[('GARY', 'FRE', 'coils')]: 525.0000
NUM-UNITS[('GARY', 'FRE', 'plate')]: 100.0000
NUM-UNITS[('GARY', 'LAN', 'bands')]: 100.0000
NUM-UNITS[('GARY', 'LAN', 'coils')]: 50.0000
NUM-UNITS[('GARY', 'STL', 'bands')]: 300.0000
NUM-UNITS[('GARY', 'STL', 'coils')]: 225.0000
NUM-UNITS[('GARY', 'STL', 'plate')]: 100.0000
NUM-UNITS[('PITT', 'DET', 'bands')]: 300.0000
NUM-UNITS[('PITT', 'DET', 'coils')]: 225.0000
NUM-UNITS[('PITT', 'FRA', 'bands')]: 25.0000
NUM-UNITS[('PITT', 'FRA', 'coils')]: 500.0000
NUM-UNITS[('PITT', 'FRA', 'plate')]: 100.0000
NUM-UNITS[('PITT', 'FRE', 'bands')]: 225.0000
NUM-UNITS[('PITT', 'FRE', 'coils')]: 325.0000
NUM-UNITS[('PITT', 'LAF', 'bands')]: 250.0000
NUM-UNITS[('PITT', 'LAF', 'coils')]: 125.0000
NUM-UNITS[('PITT', 'LAF', 'plate')]: 200.0000
NUM-UNITS[('PITT', 'STL', 'coils')]: 625.0000
USE-ROUTE[('CLEV', 'DET')]: 1.0000
USE-ROUTE[('CLEV', 'FRA')]: 1.0000
USE-ROUTE[('CLEV', 'LAF')]: 1.0000
USE-ROUTE[('CLEV', 'LAN')]: 1.0000
USE-ROUTE[('CLEV', 'STL')]: 1.0000
USE-ROUTE[('CLEV', 'WIN')]: 1.0000
USE-ROUTE[('GARY', 'FRE')]: 1.0000
USE-ROUTE[('GARY', 'LAN')]: 1.0000
USE-ROUTE[('GARY', 'STL')]: 1.0000
USE-ROUTE[('PITT', 'DET')]: 1.0000
USE-ROUTE[('PITT', 'FRA')]: 1.0000
USE-ROUTE[('PITT', 'FRE')]: 1.0000
USE-ROUTE[('PITT', 'LAF')]: 1.0000
USE-ROUTE[('PITT', 'STL')]: 1.0000
Implement with highspy#
from highspy import Highs, HighsVarType, ObjSense
from opti_extensions.highspy import addVariables as addVariablesHighs
# Instantiate model
model = Highs()
model.silent()
# Add variables
trans = addVariablesHighs(
model,
IndexSetND(ORIG, DEST, PROD, names=('ORIG', 'DEST', 'PROD')),
type=HighsVarType.kContinuous,
name_prefix='NUM-UNITS',
)
# Instead of:
# trans = model.addVariables(
# [(i, j, p) for i in ORIG for j in DEST for p in PROD],
# type=HighsVarType.kContinuous,
# name_prefix='NUM-UNITS',
# )
use = addVariablesHighs(
model,
IndexSetND(ORIG, DEST, names=('ORIG', 'DEST')),
type=HighsVarType.kInteger,
lb=0,
ub=1,
name_prefix='USE-ROUTE',
)
# Instead of:
# use = model.addVariables(
# [(i, j) for i in ORIG for j in DEST],
# type=HighsVarType.kInteger,
# lb=0,
# ub=1,
# name_prefix='USE-ROUTE',
# )
# Set objective
model.setObjective(
vcost @ trans + fcost @ use,
# Instead of:
# Highs.qsum(
# vcost[i, j, p] * trans[i, j, p]
# for i in ORIG for j in DEST for p in PROD
# )
# + Highs.qsum(
# fcost[i, j] * use[i, j]
# for i in ORIG for j in DEST
# )
ObjSense.kMinimize,
)
# Add constraints
model.addConstrs(
trans.sum(i, '*', p) == supply[i, p]
# Instead of:
# Highs.qsum(trans[i, j, p] for j in DEST) == supply[i, p]
for i in ORIG for p in PROD
)
model.addConstrs(
trans.sum('*', j, p) == demand[j, p]
# Instead of:
# Highs.qsum(trans[i, j, p] for i in ORIG) == demand[j, p]
for j in DEST for p in PROD
)
model.addConstrs(
trans.sum(i, j, '*') <= limit * use[i, j]
# Instead of:
# Highs.qsum(trans[i, j, p] for p in PROD) <= limit * use[i, j]
for i in ORIG for j in DEST
)
# Solve
model.solve()
<HighsStatus.kOk: 0>
print(f'{"var:":<15} {"value":>8}')
print('-' * 25)
for vd in (trans, use):
for idx, var in vd.items():
val = model.val(var)
if abs(val) > 1e-6:
name = f'{vd.value_name}[{idx}]:'
print(f'{name:<15} {val:>8.4f}')
var: value
-------------------------
NUM-UNITS[('CLEV', 'DET', 'bands')]: 50.0000
NUM-UNITS[('CLEV', 'DET', 'coils')]: 475.0000
NUM-UNITS[('CLEV', 'DET', 'plate')]: 100.0000
NUM-UNITS[('CLEV', 'FRA', 'bands')]: 225.0000
NUM-UNITS[('CLEV', 'FRA', 'plate')]: 50.0000
NUM-UNITS[('CLEV', 'LAF', 'coils')]: 425.0000
NUM-UNITS[('CLEV', 'LAN', 'coils')]: 350.0000
NUM-UNITS[('CLEV', 'STL', 'bands')]: 350.0000
NUM-UNITS[('CLEV', 'STL', 'coils')]: 100.0000
NUM-UNITS[('CLEV', 'STL', 'plate')]: 100.0000
NUM-UNITS[('CLEV', 'WIN', 'bands')]: 75.0000
NUM-UNITS[('CLEV', 'WIN', 'coils')]: 250.0000
NUM-UNITS[('CLEV', 'WIN', 'plate')]: 50.0000
NUM-UNITS[('GARY', 'FRE', 'coils')]: 525.0000
NUM-UNITS[('GARY', 'FRE', 'plate')]: 100.0000
NUM-UNITS[('GARY', 'LAN', 'bands')]: 100.0000
NUM-UNITS[('GARY', 'LAN', 'coils')]: 50.0000
NUM-UNITS[('GARY', 'STL', 'bands')]: 300.0000
NUM-UNITS[('GARY', 'STL', 'coils')]: 225.0000
NUM-UNITS[('GARY', 'STL', 'plate')]: 100.0000
NUM-UNITS[('PITT', 'DET', 'bands')]: 250.0000
NUM-UNITS[('PITT', 'DET', 'coils')]: 275.0000
NUM-UNITS[('PITT', 'FRA', 'bands')]: 75.0000
NUM-UNITS[('PITT', 'FRA', 'coils')]: 500.0000
NUM-UNITS[('PITT', 'FRA', 'plate')]: 50.0000
NUM-UNITS[('PITT', 'FRE', 'bands')]: 225.0000
NUM-UNITS[('PITT', 'FRE', 'coils')]: 325.0000
NUM-UNITS[('PITT', 'LAF', 'bands')]: 250.0000
NUM-UNITS[('PITT', 'LAF', 'coils')]: 75.0000
NUM-UNITS[('PITT', 'LAF', 'plate')]: 250.0000
NUM-UNITS[('PITT', 'STL', 'coils')]: 625.0000
USE-ROUTE[('CLEV', 'DET')]: 1.0000
USE-ROUTE[('CLEV', 'FRA')]: 1.0000
USE-ROUTE[('CLEV', 'LAF')]: 1.0000
USE-ROUTE[('CLEV', 'LAN')]: 1.0000
USE-ROUTE[('CLEV', 'STL')]: 1.0000
USE-ROUTE[('CLEV', 'WIN')]: 1.0000
USE-ROUTE[('GARY', 'FRE')]: 1.0000
USE-ROUTE[('GARY', 'LAN')]: 1.0000
USE-ROUTE[('GARY', 'STL')]: 1.0000
USE-ROUTE[('PITT', 'DET')]: 1.0000
USE-ROUTE[('PITT', 'FRA')]: 1.0000
USE-ROUTE[('PITT', 'FRE')]: 1.0000
USE-ROUTE[('PITT', 'LAF')]: 1.0000
USE-ROUTE[('PITT', 'STL')]: 1.0000
Total running time of the script: (0 minutes 0.327 seconds)