# The multi-commodity flow problem

This tutorial was generated using Literate.jl. Download the source as a .jl file.

This tutorial was originally contributed by Louis Luangkesorn.

This tutorial is a JuMP implementation of the multi-commodity transportation model described in AMPL: A Modeling Language for Mathematical Programming, by R. Fourer, D.M. Gay and B.W. Kernighan.

The purpose of this tutorial is to demonstrate creating a JuMP model from an SQLite database.

## Required packages

This tutorial uses the following packages

using JuMP
import DataFrames
import HiGHS
import SQLite
import Tables

const DBInterface = SQLite.DBInterface
DBInterface

## Formulation

The multi-commondity flow problem is a simple extension of The transportation problem to multiple types of products. Briefly, we start with the formulation of the transportation problem:

\begin{aligned} \min && \sum_{i \in O, j \in D} c_{i,j} x_{i,j} \\ s.t. && \sum_{j \in D} x_{i, j} \le s_i && \forall i \in O \\ && \sum_{i \in O} x_{i, j} = d_j && \forall j \in D \\ && x_{i, j} \ge 0 && \forall i \in O, j \in D \end{aligned}

but introduce a set of products $P$, resulting in:

\begin{aligned} \min && \sum_{i \in O, j \in D, k \in P} c_{i,j,k} x_{i,j,k} \\ s.t. && \sum_{j \in D} x_{i, j, k} \le s_{i,k} && \forall i \in O, k \in P \\ && \sum_{i \in O} x_{i, j, k} = d_{j,k} && \forall j \in D, k \in P \\ && x_{i, j,k} \ge 0 && \forall i \in O, j \in D, k \in P \\ && \sum_{k \in P} x_{i, j, k} \le u_{i,j} && \forall i \in O, j \in D \end{aligned}

Note that the last constraint is new; it says that there is a maximum quantity of goods (of any type) that can be transported from origin $i$ to destination $j$.

## Data

For the purpose of this tutorial, the JuMP repository contains an example database called multi.sqlite.

filename = joinpath(@__DIR__, "multi.sqlite");

To run locally, download multi.sqlite and update filename appropriately.

Load the database using SQLite.DB:

db = SQLite.DB(filename)
SQLite.DB("/home/runner/work/JuMP.jl/JuMP.jl/docs/build/tutorials/linear/multi.sqlite")

A quick way to see the schema of the database is via SQLite.tables:

SQLite.tables(db)
5-element Vector{SQLite.DBTable}:
SQLite.DBTable("locations", Tables.Schema:
:location  Union{Missing, String}
:type      Union{Missing, String})
SQLite.DBTable("products", Tables.Schema:
:product  Union{Missing, String})
SQLite.DBTable("supply", Tables.Schema:
:origin   Union{Missing, String}
:product  Union{Missing, String}
:supply   Union{Missing, Float64})
SQLite.DBTable("demand", Tables.Schema:
:destination  Union{Missing, String}
:product      Union{Missing, String}
:demand       Union{Missing, Float64})
SQLite.DBTable("cost", Tables.Schema:
:origin       Union{Missing, String}
:destination  Union{Missing, String}
:product      Union{Missing, String}
:cost         Union{Missing, Float64})

We interact with the database by executing queries, and then piping the results to an appropriate table. One example is a DataFrame:

DBInterface.execute(db, "SELECT * FROM locations") |> DataFrames.DataFrame
10×2 DataFrame
Rowlocationtype
StringString
1GARYorigin
2CLEVorigin
3PITTorigin
5DETdestination
6LANdestination
7WINdestination
8STLdestination
9FREdestination
10LAFdestination

But other table types are supported, such as Tables.rowtable:

DBInterface.execute(db, "SELECT * FROM locations") |> Tables.rowtable
10-element Vector{NamedTuple{(:location, :type), Tuple{String, String}}}:
(location = "GARY", type = "origin")
(location = "CLEV", type = "origin")
(location = "PITT", type = "origin")
(location = "FRA", type = "destination")
(location = "DET", type = "destination")
(location = "LAN", type = "destination")
(location = "WIN", type = "destination")
(location = "STL", type = "destination")
(location = "FRE", type = "destination")
(location = "LAF", type = "destination")

A rowtable is a Vector of NamedTuples.

You can construct more complicated SQL queries:

origins =
DBInterface.execute(
db,
"SELECT location FROM locations WHERE type = \"origin\"",
) |> Tables.rowtable
3-element Vector{NamedTuple{(:location,), Tuple{String}}}:
(location = "GARY",)
(location = "CLEV",)
(location = "PITT",)

But for our purpose, we just want the list of strings:

origins = map(y -> y.location, origins)
3-element Vector{String}:
"GARY"
"CLEV"
"PITT"

We can compose these two operations to get a list of destinations:

destinations =
DBInterface.execute(
db,
"SELECT location FROM locations WHERE type = \"destination\"",
) |>
Tables.rowtable |>
x -> map(y -> y.location, x)
7-element Vector{String}:
"FRA"
"DET"
"LAN"
"WIN"
"STL"
"FRE"
"LAF"

And a list of products from our products table:

products =
DBInterface.execute(db, "SELECT product FROM products") |>
Tables.rowtable |>
x -> map(y -> y.product, x)
3-element Vector{String}:
"bands"
"coils"
"plate"

## JuMP formulation

We start by creating a model and our decision variables:

model = Model(HiGHS.Optimizer)
set_silent(model)
@variable(model, x[origins, destinations, products] >= 0)
3-dimensional DenseAxisArray{VariableRef,3,...} with index sets:
Dimension 1, ["GARY", "CLEV", "PITT"]
Dimension 2, ["FRA", "DET", "LAN", "WIN", "STL", "FRE", "LAF"]
Dimension 3, ["bands", "coils", "plate"]
And data, a 3×7×3 Array{VariableRef, 3}:
[:, :, "bands"] =
x[GARY,FRA,bands]  x[GARY,DET,bands]  …  x[GARY,LAF,bands]
x[CLEV,FRA,bands]  x[CLEV,DET,bands]     x[CLEV,LAF,bands]
x[PITT,FRA,bands]  x[PITT,DET,bands]     x[PITT,LAF,bands]

[:, :, "coils"] =
x[GARY,FRA,coils]  x[GARY,DET,coils]  …  x[GARY,LAF,coils]
x[CLEV,FRA,coils]  x[CLEV,DET,coils]     x[CLEV,LAF,coils]
x[PITT,FRA,coils]  x[PITT,DET,coils]     x[PITT,LAF,coils]

[:, :, "plate"] =
x[GARY,FRA,plate]  x[GARY,DET,plate]  …  x[GARY,LAF,plate]
x[CLEV,FRA,plate]  x[CLEV,DET,plate]     x[CLEV,LAF,plate]
x[PITT,FRA,plate]  x[PITT,DET,plate]     x[PITT,LAF,plate]

One approach when working with databases is to extract all of the data into a Julia datastructure. For example, let's pull the cost table into a DataFrame and then construct our objective by iterating over the rows of the DataFrame:

cost = DBInterface.execute(db, "SELECT * FROM cost") |> DataFrames.DataFrame
@objective(
model,
Max,
sum(r.cost * x[r.origin, r.destination, r.product] for r in eachrow(cost)),
);

If we don't want to use a DataFrame, we can use a Tables.rowtable instead:

supply = DBInterface.execute(db, "SELECT * FROM supply") |> Tables.rowtable
for r in supply
@constraint(model, sum(x[r.origin, :, r.product]) <= r.supply)
end

Another approach is to execute the query, and then to iterate through the rows of the query using Tables.rows:

demand = DBInterface.execute(db, "SELECT * FROM demand")
for r in Tables.rows(demand)
@constraint(model, sum(x[:, r.destination, r.product]) == r.demand)
end
Warning

Iterating through the rows of a query result works by incrementing a cursor inside the database. As a consequence, you cannot call Tables.rows twice on the same query result.

The SQLite queries can be arbitrarily complex. For example, here's a query which builds every possible origin-destination pair:

od_pairs = DBInterface.execute(
db,
"""
SELECT a.location as 'origin',
b.location as 'destination'
FROM locations a
INNER JOIN locations b
ON a.type = 'origin' AND b.type = 'destination'
""",
)
SQLite.Query{false}(SQLite.Stmt(SQLite.DB("/home/runner/work/JuMP.jl/JuMP.jl/docs/build/tutorials/linear/multi.sqlite"), Base.RefValue{Ptr{SQLite.C.sqlite3_stmt}}(Ptr{SQLite.C.sqlite3_stmt} @0x0000000040032908), Dict{Int64, Any}()), Base.RefValue{Int32}(100), [:origin, :destination], Type[Union{Missing, String}, Union{Missing, String}], Dict(:origin => 1, :destination => 2), Base.RefValue{Int64}(0))

With a constraint that we cannot send more than 625 units between each pair:

for r in Tables.rows(od_pairs)
@constraint(model, sum(x[r.origin, r.destination, :]) <= 625)
end

## Solution

Finally, we can optimize the model:

optimize!(model)
solution_summary(model)
* Solver : HiGHS

* Status
Result count       : 1
Termination status : OPTIMAL
Message from the solver:
"kHighsModelStatusOptimal"

* Candidate solution (result #1)
Primal status      : FEASIBLE_POINT
Dual status        : FEASIBLE_POINT
Objective value    : 2.25700e+05
Objective bound    : 2.25700e+05
Relative gap       : Inf
Dual objective value : 2.25700e+05

* Work counters
Solve time (sec)   : 8.22783e-04
Simplex iterations : 54
Barrier iterations : 0
Node count         : -1


and print the solution:

begin
println("         ", join(products, ' '))
for o in origins, d in destinations
v = lpad.([round(Int, value(x[o, d, p])) for p in products], 5)
println(o, " ", d, " ", join(replace.(v, "   0" => "  . "), " "))
end
end
         bands coils plate
GARY FRA    25   500   100
GARY DET   125    .     50
GARY LAN    .     .     .
GARY WIN    .     .     50
GARY STL   250   300    .
GARY FRE    .     .     .
GARY LAF    .     .     .
CLEV FRA   275    .     .
CLEV DET   100   200    50
CLEV LAN   100    .     .
CLEV WIN    .     .     .
CLEV STL    .    625    .
CLEV FRE   225   400    .
CLEV LAF    .    375   250
PITT FRA    .     .     .
PITT DET    75   550    .
PITT LAN    .    400    .
PITT WIN    75   250    .
PITT STL   400    25   200
PITT FRE    .    450   100
PITT LAF   250   125    .