Задача о диете

Цель данного руководства — продемонстрировать, как объекты DataFrame включаются в модель JuMP. В качестве примера мы используем классическую задачу о диете Стиглера.

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

Для этого руководства требуются следующие пакеты.

using JuMP
import CSV
import DataFrames
import HiGHS
import Test

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

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

У каждого продукта есть своя стоимость , а также профиль макронутриентов для каждого макронутриента .

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

Кроме того, поскольку нас интересует экономия, мы ищем решение с минимальной стоимостью.

Приложив немного усилий, мы можем сформулировать нашу задачу о приготовлении обеда в виде следующей линейной программы:

В оставшейся части этого руководства мы построим и решим эту задачу в JuMP, чтобы узнать, что приготовить на обед.

Данные

Во-первых, нам нужны данные для задачи. В этом руководстве мы запишем CSV-файлы во временный каталог из Julia. Если у вас уже есть нужные файлы, вы можете указать их имена.

dir = mktempdir()

"/tmp/jl_XK6tL2"

Первый файл — это список продуктов с указанием их макронутриентного профиля:

food_csv_filename = joinpath(dir, "diet_foods.csv")
open(food_csv_filename, "w") do io
    write(
        io,
        """
        name,cost,calories,protein,fat,sodium
        hamburger,2.49,410,24,26,730
        chicken,2.89,420,32,10,1190
        hot dog,1.50,560,20,32,1800
        fries,1.89,380,4,19,270
        macaroni,2.09,320,12,10,930
        pizza,1.99,320,15,12,820
        salad,2.49,320,31,12,1230
        milk,0.89,100,8,2.5,125
        ice cream,1.59,330,8,10,180
        """,
    )
    return
end
foods = CSV.read(food_csv_filename, DataFrames.DataFrame)

9×6 DataFrame

Row	name	cost	calories	protein	fat	sodium
	String15	Float64	Int64	Int64	Float64	Int64
1	hamburger	2.49	410	24	26.0	730
2	chicken	2.89	420	32	10.0	1190
3	hot dog	1.5	560	20	32.0	1800
4	fries	1.89	380	4	19.0	270
5	macaroni	2.09	320	12	10.0	930
6	pizza	1.99	320	15	12.0	820
7	salad	2.49	320	31	12.0	1230
8	milk	0.89	100	8	2.5	125
9	ice cream	1.59	330	8	10.0	180

Здесь $F$ — это foods.name, а $c_f$ — foods.cost. (Мы немного вольно обращаемся с термином «макронутриент», включая в него калории и натрий.)

Нам также нужны минимальный и максимальный пределы:

nutrient_csv_filename = joinpath(dir, "diet_nutrient.csv")
open(nutrient_csv_filename, "w") do io
    write(
        io,
        """
        nutrient,min,max
        calories,1800,2200
        protein,91,
        fat,0,65
        sodium,0,1779
        """,
    )
    return
end
limits = CSV.read(nutrient_csv_filename, DataFrames.DataFrame)

4×3 DataFrame

Row	nutrient	min	max
	String15	Int64	Int64?
1	calories	1800	2200
2	protein	91	missing
3	fat	0	65
4	sodium	0	1779

Для белков максимальный предел отсутствует. Исправим это с помощью coalesce:

limits.max = coalesce.(limits.max, Inf)
limits

4×3 DataFrame

Row	nutrient	min	max
	String15	Int64	Real
1	calories	1800	2200
2	protein	91	Inf
3	fat	0	65
4	sodium	0	1779

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

Теперь мы готовы преобразовать математическую формулировку в модель JuMP.

Сначала создадим модель JuMP. Так как программа линейная, мы используем оптимизатор HiGHS:

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

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

@variable(model, x[foods.name] >= 0)

1-dimensional DenseAxisArray{VariableRef,1,...} with index sets:
    Dimension 1, String15["hamburger", "chicken", "hot dog", "fries", "macaroni", "pizza", "salad", "milk", "ice cream"]
And data, a 9-element Vector{VariableRef}:
 x[hamburger]
 x[chicken]
 x[hot dog]
 x[fries]
 x[macaroni]
 x[pizza]
 x[salad]
 x[milk]
 x[ice cream]

Чтобы упростить себе задачу в дальнейшем, мы сохраним вектор как новый столбец x в DataFrame foods. Так как x — это массив DenseAxisArray, сначала его нужно преобразовать в Array:

foods.x = Array(x)

9-element Vector{VariableRef}:
 x[hamburger]
 x[chicken]
 x[hot dog]
 x[fries]
 x[macaroni]
 x[pizza]
 x[salad]
 x[milk]
 x[ice cream]

Наша цель — минимизировать общую стоимость продуктов:

@objective(model, Min, sum(foods.cost .* foods.x));

Далее необходимо добавить ограничение, чтобы общее потребление каждого компонента находилось в пределах, указанных в DataFrame limits:

@constraint(
    model,
    [row in eachrow(limits)],
    row.min <= sum(foods[!, row.nutrient] .* foods.x) <= row.max,
);

Как выглядит наша модель?

print(model)

Min 2.49 x[hamburger] + 2.89 x[chicken] + 1.5 x[hot dog] + 1.89 x[fries] + 2.09 x[macaroni] + 1.99 x[pizza] + 2.49 x[salad] + 0.89 x[milk] + 1.59 x[ice cream]
Subject to
 410 x[hamburger] + 420 x[chicken] + 560 x[hot dog] + 380 x[fries] + 320 x[macaroni] + 320 x[pizza] + 320 x[salad] + 100 x[milk] + 330 x[ice cream] ∈ [1800, 2200]
 24 x[hamburger] + 32 x[chicken] + 20 x[hot dog] + 4 x[fries] + 12 x[macaroni] + 15 x[pizza] + 31 x[salad] + 8 x[milk] + 8 x[ice cream] ∈ [91, Inf]
 26 x[hamburger] + 10 x[chicken] + 32 x[hot dog] + 19 x[fries] + 10 x[macaroni] + 12 x[pizza] + 12 x[salad] + 2.5 x[milk] + 10 x[ice cream] ∈ [0, 65]
 730 x[hamburger] + 1190 x[chicken] + 1800 x[hot dog] + 270 x[fries] + 930 x[macaroni] + 820 x[pizza] + 1230 x[salad] + 125 x[milk] + 180 x[ice cream] ∈ [0, 1779]
 x[hamburger] ≥ 0
 x[chicken] ≥ 0
 x[hot dog] ≥ 0
 x[fries] ≥ 0
 x[macaroni] ≥ 0
 x[pizza] ≥ 0
 x[salad] ≥ 0
 x[milk] ≥ 0
 x[ice cream] ≥ 0

Решение

Давайте выполним оптимизацию и посмотрим на решение:

optimize!(model)
@assert 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.18289e+01
  Objective bound    : 0.00000e+00
  Relative gap       : Inf
  Dual objective value : 1.18289e+01

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

Найдено оптимальное решение. Давайте посмотрим, как оно выглядит:

for row in eachrow(foods)
    println(row.name, " = ", value(row.x))
end

hamburger = 0.6045138888888871
chicken = 0.0
hot dog = 0.0
fries = 0.0
macaroni = 0.0
pizza = 0.0
salad = 0.0
milk = 6.9701388888888935
ice cream = 2.5913194444444447

Молока и мороженого очень много, а гамбургера, к сожалению, всего 0.6.

С помощью функции Containers.rowtable можно также легко преобразовать результат в DataFrame:

table = Containers.rowtable(value, x; header = [:food, :quantity])
solution = DataFrames.DataFrame(table)

9×2 DataFrame

Row	food	quantity
	String15	Float64
1	hamburger	0.604514
2	chicken	0.0
3	hot dog	0.0
4	fries	0.0
5	macaroni	0.0
6	pizza	0.0
7	salad	0.0
8	milk	6.97014
9	ice cream	2.59132

Это упрощает анализ решения:

filter!(row -> row.quantity > 0.0, solution)

3×2 DataFrame

Row	food	quantity
	String15	Float64
1	hamburger	0.604514
2	milk	6.97014
3	ice cream	2.59132

Изменение задачи

JuMP позволяет легко изменить существующую модель, добавив дополнительные ограничения. Давайте посмотрим, что произойдет, если добавить ограничение, согласно которому можно купить в сумме не более 6 единиц молока и мороженого.

dairy_foods = ["milk", "ice cream"]
is_dairy = map(name -> name in dairy_foods, foods.name)
dairy_constraint = @constraint(model, sum(foods[is_dairy, :x]) <= 6)
optimize!(model)
Test.@test !is_solved_and_feasible(model)
Test.@test termination_status(model) == INFEASIBLE
Test.@test primal_status(model) == NO_SOLUTION
solution_summary(model)

* Solver : HiGHS

* Status
  Result count       : 1
  Termination status : INFEASIBLE
  Message from the solver:
  "kHighsModelStatusInfeasible"

* Candidate solution (result #1)
  Primal status      : NO_SOLUTION
  Dual status        : INFEASIBILITY_CERTIFICATE
  Objective value    : 1.18289e+01
  Objective bound    : 0.00000e+00
  Relative gap       : Inf
  Dual objective value : 3.56146e+00

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

Допустимого решения задачи не существует. Похоже, на обед придется есть одно мороженое.

Дальнейшие действия

Можно удалить ограничение с помощью delete(model, dairy_constraint). Можно ли добавить другое ограничение, чтобы в рационе было меньше молочных продуктов?
Некоторые продукты (например, гамбургеры) являются дискретными. Вы можете использовать set_integer, чтобы переменная принимала только целочисленные значения. Что произойдет с решением в этом случае?