The multi-commodity flow problem

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:

db = SQLite.DB(joinpath(@__DIR__, "multi.sqlite"))

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

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

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} @0x0000000028120f48), 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.46863e-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    .