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
import Test

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:

but introduce a set of products , resulting in:

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 to destination .

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/docs/julia_docs/JuMP.jl-master/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

Row	location	type
	String	String
1	GARY	origin
2	CLEV	origin
3	PITT	origin
4	FRA	destination
5	DET	destination
6	LAN	destination
7	WIN	destination
8	STL	destination
9	FRE	destination
10	LAF	destination

Row

location

type

String

GARY

origin

CLEV

origin

PITT

origin

FRA

destination

DET

destination

LAN

destination

WIN

destination

STL

destination

FRE

destination

LAF

destination

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

DBInterface.execute(db, "SELECT * FROM locations") |> Tables.rowtable

10-element Vector{@NamedTuple{location::String, type::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::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

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/docs/julia_docs/JuMP.jl-master/docs/build/tutorials/linear/multi.sqlite"), Base.RefValue{Ptr{SQLite.C.sqlite3_stmt}}(Ptr{SQLite.C.sqlite3_stmt} @0x000000002590ab58), 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)
Test.@test is_solved_and_feasible(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)   : 7.43151e-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    .