Задача о сетевом потоке нескольких видов товаров

В этом руководстве рассматривается разновидность задачи из раздела The multi-commodity flow problem, в которой граф является не двудольным, а сетевым.

Цель данного руководства — продемонстрировать стиль моделирования с использованием алгебры отношений.

Требуемые пакеты

В этом руководстве используются следующие пакеты.

using JuMP
import DataFrames
import HiGHS
import SQLite
import SQLite.DBInterface
import Test

Формулировка

Задача о сетевом потоке нескольких видов товаров является расширением задачи из раздела The multi-commodity flow problem, где вместо двудольного графа узлов предложения и спроса граф может содержать набор узлов , каждый из которых имеет (возможно, нулевой) потенциал предложения и (возможно, нулевой) спрос для каждого вида товара . Узлы соединены набором ребер , которые имеют стоимость транспортировки и общую пропускную способность .

Наша задача — выбрать оптимальный объем поставки для каждого узла , а также оптимальный вариант перевозки с минимальной общей стоимостью.

Математическая формулировка:

Цель данного руководства — продемонстрировать построение этой модели с использованием алгебры отношений вместо прямого преобразования математических операций суммирования в код.

Данные

Для целей данного руководства в репозитории JuMP есть пример базы данных под названием commodity_nz.db:

filename = joinpath(@__DIR__, "commodity_nz.db");

Для локального выполнения скачайте файл commodity_nz.db и соответствующим образом измените значение filename.

Загрузите базу данных с помощью SQLite.DB:

db = SQLite.DB(filename)

SQLite.DB("/home/docs/julia_docs/JuMP.jl-master/docs/build/tutorials/linear/commodity_nz.db")

Схему базы данных можно быстро просмотреть с помощью SQLite.tables:

SQLite.tables(db)

4-element Vector{SQLite.DBTable}:
 SQLite.DBTable("products", Tables.Schema:
 :product      Union{Missing, String}
 :cost_per_km  Union{Missing, Float64})
 SQLite.DBTable("shipping", Tables.Schema:
 :origin       Union{Missing, String}
 :destination  Union{Missing, String}
 :product      Union{Missing, String}
 :distance_km  Union{Missing, Float64})
 SQLite.DBTable("supply", Tables.Schema:
 :origin    Union{Missing, String}
 :product   Union{Missing, String}
 :capacity  Union{Missing, Float64}
 :cost      Union{Missing, Float64})
 SQLite.DBTable("demand", Tables.Schema:
 :destination  Union{Missing, String}
 :product      Union{Missing, String}
 :demand       Union{Missing, Float64})

Взаимодействие с базой данных осуществляется путем выполнения запросов и загрузки результатов в DataFrame:

function get_table(db, table)
    query = DBInterface.execute(db, "SELECT * FROM $table")
    return DataFrames.DataFrame(query)
end

get_table (generic function with 1 method)

Таблица shipping содержит набор ребер графа и их расстояний:

df_shipping = get_table(db, "shipping")

16×4 DataFrame

Row	origin	destination	product	distance_km
	String	String	String	Float64
1	auckland	waikato	milk	112.0
2	auckland	tauranga	milk	225.0
3	auckland	christchurch	milk	1070.0
4	waikato	auckland	milk	112.0
5	waikato	tauranga	milk	107.0
6	waikato	wellington	milk	392.0
7	tauranga	auckland	milk	225.0
8	tauranga	waikato	milk	107.0
9	christchurch	auckland	milk	1070.0
10	auckland	waikato	kiwifruit	112.0
11	auckland	christchurch	kiwifruit	1070.0
12	waikato	auckland	kiwifruit	112.0
13	waikato	wellington	kiwifruit	392.0
14	tauranga	auckland	kiwifruit	225.0
15	tauranga	waikato	kiwifruit	107.0
16	christchurch	auckland	kiwifruit	1070.0

В таблице products указана стоимость доставки за километр для каждого товара:

df_products = get_table(db, "products")

2×2 DataFrame

Row	product	cost_per_km
	String	Float64
1	milk	0.001
2	kiwifruit	0.01

Таблица supply содержит потенциал предложения каждого узла, а также стоимость:

df_supply = get_table(db, "supply")

4×4 DataFrame

Row	origin	product	capacity	cost
	String	String	Float64	Float64
1	waikato	milk	10.0	0.5
2	tauranga	milk	6.0	1.0
3	tauranga	kiwifruit	26.0	1.0
4	christchurch	milk	10.0	0.6

Таблица demand содержит спрос каждого узла:

df_demand = get_table(db, "demand")

8×3 DataFrame

Row	destination	product	demand
	String	String	Float64
1	auckland	milk	16.0
2	auckland	kiwifruit	16.0
3	tauranga	milk	2.0
4	tauranga	kiwifruit	2.0
5	wellington	milk	2.0
6	wellington	kiwifruit	2.0
7	christchurch	milk	4.0
8	christchurch	kiwifruit	4.0

Формулировка на языке JuMP

Начнем с создания модели и переменных решения:

model = Model(HiGHS.Optimizer)
set_silent(model)

Для принятия решений о транспортировке мы создадим в df_shipping новый столбец под названием x_flow с одной неотрицательной переменной решения для каждой строки:

df_shipping.x_flow = @variable(model, x[1:size(df_shipping, 1)] >= 0)
df_shipping

16×5 DataFrame

Row	origin	destination	product	distance_km	x_flow
	String	String	String	Float64	GenericV…
1	auckland	waikato	milk	112.0	x
2	auckland	tauranga	milk	225.0	x
3	auckland	christchurch	milk	1070.0	x
4	waikato	auckland	milk	112.0	x
5	waikato	tauranga	milk	107.0	x
6	waikato	wellington	milk	392.0	x
7	tauranga	auckland	milk	225.0	x
8	tauranga	waikato	milk	107.0	x
9	christchurch	auckland	milk	1070.0	x
10	auckland	waikato	kiwifruit	112.0	x
11	auckland	christchurch	kiwifruit	1070.0	x
12	waikato	auckland	kiwifruit	112.0	x
13	waikato	wellington	kiwifruit	392.0	x
14	tauranga	auckland	kiwifruit	225.0	x
15	tauranga	waikato	kiwifruit	107.0	x
16	christchurch	auckland	kiwifruit	1070.0	x

Для узлов предложения мы добавим переменную в каждую строку, а затем установим верхнюю границу равной потенциалу предложения узла:

df_supply.x_supply = @variable(model, s[1:size(df_supply, 1)] >= 0)
set_upper_bound.(df_supply.x_supply, df_supply.capacity)
df_supply

4×5 DataFrame

Row	origin	product	capacity	cost	x_supply
	String	String	Float64	Float64	GenericV…
1	waikato	milk	10.0	0.5	s
2	tauranga	milk	6.0	1.0	s
3	tauranga	kiwifruit	26.0	1.0	s
4	christchurch	milk	10.0	0.6	s

Наша цель — минимизировать стоимость доставки и стоимость поставки. Чтобы вычислить стоимость потока, необходимо объединить таблицу доставки, которая содержит distance_km, с таблицей товаров, которая содержит cost_per_km:

df_cost = DataFrames.leftjoin(df_shipping, df_products; on = [:product])
df_cost.flow_cost = df_cost.cost_per_km .* df_cost.distance_km
df_cost

16×7 DataFrame

Row	origin	destination	product	distance_km	x_flow	cost_per_km	flow_cost
	String	String	String	Float64	GenericV…	Float64?	Float64
1	auckland	waikato	milk	112.0	x	0.001	0.112
2	auckland	tauranga	milk	225.0	x	0.001	0.225
3	auckland	christchurch	milk	1070.0	x	0.001	1.07
4	waikato	auckland	milk	112.0	x	0.001	0.112
5	waikato	tauranga	milk	107.0	x	0.001	0.107
6	waikato	wellington	milk	392.0	x	0.001	0.392
7	tauranga	auckland	milk	225.0	x	0.001	0.225
8	tauranga	waikato	milk	107.0	x	0.001	0.107
9	christchurch	auckland	milk	1070.0	x	0.001	1.07
10	auckland	waikato	kiwifruit	112.0	x	0.01	1.12
11	auckland	christchurch	kiwifruit	1070.0	x	0.01	10.7
12	waikato	auckland	kiwifruit	112.0	x	0.01	1.12
13	waikato	wellington	kiwifruit	392.0	x	0.01	3.92
14	tauranga	auckland	kiwifruit	225.0	x	0.01	2.25
15	tauranga	waikato	kiwifruit	107.0	x	0.01	1.07
16	christchurch	auckland	kiwifruit	1070.0	x	0.01	10.7

Затем с помощью линейной алгебры можно вычислить внутреннее произведение двух столбцов:

@objective(
    model,
    Min,
    df_cost.flow_cost' * df_shipping.x_flow +
    df_supply.cost' * df_supply.x_supply
);

Для пропускной способности потока по каждому ребру мы используем DataFrames.groupby для разделения переменных потока на основе :origin и :destination, а затем ограничим их сумму так, чтобы она была меньше фиксированной пропускной способности.

capacity = 30
for df in DataFrames.groupby(df_shipping, [:origin, :destination])
    @constraint(model, sum(df.x_flow) <= capacity)
end

Для каждого узла графа необходимо вычислить ограничение массового баланса, которое гласит, что для каждого товара сумма предложения и потока в узел минус поток из узла равна спросу.

Выражение для потока из каждого узла можно вычислить, используя DataFrames.groupby применительно к столбцам origin и product таблицы df_shipping:

df_flow_out = DataFrames.DataFrame(
    (node = i.origin, product = i.product, x_flow_out = sum(df.x_flow)) for
    (i, df) in pairs(DataFrames.groupby(df_shipping, [:origin, :product]))
)

8×3 DataFrame

Row	node	product	x_flow_out
	String	String	AffExpr
1	auckland	milk	x[1] + x[2] + x[3]
2	waikato	milk	x[4] + x[5] + x[6]
3	tauranga	milk	x[7] + x[8]
4	christchurch	milk	x[9]
5	auckland	kiwifruit	x[10] + x[11]
6	waikato	kiwifruit	x[12] + x[13]
7	tauranga	kiwifruit	x[14] + x[15]
8	christchurch	kiwifruit	x[16]

Row

node

product

x_flow_out

String

AffExpr

auckland

milk

x[1] + x[2] + x[3]

waikato

milk

x[4] + x[5] + x[6]

tauranga

milk

x[7] + x[8]

christchurch

milk

x[9]

auckland

kiwifruit

x[10] + x[11]

waikato

kiwifruit

x[12] + x[13]

tauranga

kiwifruit

x[14] + x[15]

christchurch

kiwifruit

x[16]

Выражение для потока в каждый узел можно вычислить, используя DataFrames.groupby применительно к столбцам destination и product таблицы df_shipping:

df_flow_in = DataFrames.DataFrame(
    (node = i.destination, product = i.product, x_flow_in = sum(df.x_flow))
    for (i, df) in
    pairs(DataFrames.groupby(df_shipping, [:destination, :product]))
)

9×3 DataFrame

Row	node	product	x_flow_in
	String	String	AffExpr
1	waikato	milk	x[1] + x[8]
2	tauranga	milk	x[2] + x[5]
3	christchurch	milk	x[3]
4	auckland	milk	x[4] + x[7] + x[9]
5	wellington	milk	x[6]
6	waikato	kiwifruit	x[10] + x[15]
7	christchurch	kiwifruit	x[11]
8	auckland	kiwifruit	x[12] + x[14] + x[16]
9	wellington	kiwifruit	x[13]

Row

node

product

x_flow_in

String

AffExpr

waikato

milk

x[1] + x[8]

tauranga

milk

x[2] + x[5]

christchurch

milk

x[3]

auckland

milk

x[4] + x[7] + x[9]

wellington

milk

x[6]

waikato

kiwifruit

x[10] + x[15]

christchurch

kiwifruit

x[11]

auckland

kiwifruit

x[12] + x[14] + x[16]

wellington

kiwifruit

x[13]

Объединить таблицы можно с помощью DataFrames.outerjoin. Здесь приходится использовать outerjoin, так как некоторые строки могут быть пропущены.

df = DataFrames.outerjoin(df_flow_in, df_flow_out; on = [:node, :product])

10×4 DataFrame

Row	node	product	x_flow_in	x_flow_out
	String	String	AffExpr?	AffExpr?
1	waikato	milk	x[1] + x[8]	x[4] + x[5] + x[6]
2	tauranga	milk	x[2] + x[5]	x[7] + x[8]
3	christchurch	milk	x[3]	x[9]
4	auckland	milk	x[4] + x[7] + x[9]	x[1] + x[2] + x[3]
5	waikato	kiwifruit	x[10] + x[15]	x[12] + x[13]
6	christchurch	kiwifruit	x[11]	x[16]
7	auckland	kiwifruit	x[12] + x[14] + x[16]	x[10] + x[11]
8	wellington	milk	x[6]	missing
9	wellington	kiwifruit	x[13]	missing
10	tauranga	kiwifruit	missing	x[14] + x[15]

Row

node

product

x_flow_in

x_flow_out

String

AffExpr?

waikato

milk

x[1] + x[8]

x[4] + x[5] + x[6]

tauranga

milk

x[2] + x[5]

x[7] + x[8]

christchurch

milk

x[3]

x[9]

auckland

milk

x[4] + x[7] + x[9]

x[1] + x[2] + x[3]

waikato

kiwifruit

x[10] + x[15]

x[12] + x[13]

christchurch

kiwifruit

x[11]

x[16]

auckland

kiwifruit

x[12] + x[14] + x[16]

x[10] + x[11]

wellington

milk

x[6]

missing

wellington

kiwifruit

x[13]

missing

tauranga

kiwifruit

missing

x[14] + x[15]

Далее необходимо присоединить столбец предложения:

df = DataFrames.leftjoin(
    df,
    DataFrames.select(df_supply, [:origin, :product, :x_supply]);
    on = [:node => :origin, :product],
)

10×5 DataFrame

Row	node	product	x_flow_in	x_flow_out	x_supply
	String	String	AffExpr?	AffExpr?	GenericV…?
1	waikato	milk	x[1] + x[8]	x[4] + x[5] + x[6]	s[1]
2	tauranga	milk	x[2] + x[5]	x[7] + x[8]	s[2]
3	christchurch	milk	x[3]	x[9]	s[4]
4	tauranga	kiwifruit	missing	x[14] + x[15]	s[3]
5	auckland	milk	x[4] + x[7] + x[9]	x[1] + x[2] + x[3]	missing
6	waikato	kiwifruit	x[10] + x[15]	x[12] + x[13]	missing
7	christchurch	kiwifruit	x[11]	x[16]	missing
8	auckland	kiwifruit	x[12] + x[14] + x[16]	x[10] + x[11]	missing
9	wellington	milk	x[6]	missing	missing
10	wellington	kiwifruit	x[13]	missing	missing

Row

node

product

x_flow_in

x_flow_out

x_supply

String

AffExpr?

GenericV…?

waikato

milk

x[1] + x[8]

x[4] + x[5] + x[6]

s[1]

tauranga

milk

x[2] + x[5]

x[7] + x[8]

s[2]

christchurch

milk

x[3]

x[9]

s[4]

tauranga

kiwifruit

missing

x[14] + x[15]

s[3]

auckland

milk

x[4] + x[7] + x[9]

x[1] + x[2] + x[3]

missing

waikato

kiwifruit

x[10] + x[15]

x[12] + x[13]

missing

christchurch

kiwifruit

x[11]

x[16]

missing

auckland

kiwifruit

x[12] + x[14] + x[16]

x[10] + x[11]

missing

wellington

milk

x[6]

missing

wellington

kiwifruit

x[13]

missing

и столбец спроса:

df = DataFrames.leftjoin(
    df,
    DataFrames.select(df_demand, [:destination, :product, :demand]);
    on = [:node => :destination, :product],
)

10×6 DataFrame

Row	node	product	x_flow_in	x_flow_out	x_supply	demand
	String	String	AffExpr?	AffExpr?	GenericV…?	Float64?
1	tauranga	milk	x[2] + x	x[7] + x	s	2.0
2	christchurch	milk	x	x	s	4.0
3	tauranga	kiwifruit	missing	x[14] + x	s	2.0
4	auckland	milk	x[4] + x[7] + x	x[1] + x[2] + x	missing	16.0
5	christchurch	kiwifruit	x	x	missing	4.0
6	auckland	kiwifruit	x[12] + x[14] + x	x[10] + x	missing	16.0
7	wellington	milk	x	missing	missing	2.0
8	wellington	kiwifruit	x	missing	missing	2.0
9	waikato	milk	x[1] + x	x[4] + x[5] + x	s	missing
10	waikato	kiwifruit	x[10] + x	x[12] + x	missing	missing

Теперь мы готовы добавить ограничение массового баланса. Так как в некоторых строках есть значения missing, необходимо использовать coalesce для преобразования missing в числовое значение:

@constraint(
    model,
    [r in eachrow(df)],
    coalesce(r.x_supply, 0.0) + coalesce(r.x_flow_in, 0.0) -
    coalesce(r.x_flow_out, 0.0) == coalesce(r.demand, 0.0),
);

Решение

Наконец, можно оптимизировать модель:

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    : 1.43228e+02
  Objective bound    : 0.00000e+00
  Relative gap       : Inf
  Dual objective value : 1.43228e+02

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

Обновляем решение в объектах DataFrame:

df_shipping.x_flow = value.(df_shipping.x_flow)
df_supply.x_supply = value.(df_supply.x_supply);

И выводим оптимальное решение для потоков:

DataFrames.select(
    filter!(row -> row.x_flow > 0.0, df_shipping),
    [:origin, :destination, :product, :x_flow],
)

9×4 DataFrame

Row	origin	destination	product	x_flow
	String	String	String	Float64
1	waikato	auckland	milk	10.0
2	waikato	wellington	milk	2.0
3	tauranga	auckland	milk	2.0
4	tauranga	waikato	milk	2.0
5	christchurch	auckland	milk	4.0
6	auckland	christchurch	kiwifruit	4.0
7	waikato	auckland	kiwifruit	20.0
8	waikato	wellington	kiwifruit	2.0
9	tauranga	waikato	kiwifruit	22.0