Многоцелевая задача о рюкзаке

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

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

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

using JuMP
import HiGHS
import MultiObjectiveAlgorithms as MOA
import Plots

MultiObjectiveAlgorithms.jl — это пакет, реализующий различные алгоритмы для решения задач многоцелевой оптимизации. Так как у пакета длинное имя, мы импортируем его как MOA.

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

Задача о рюкзаке — это классическая задача частично целочисленного программирования. Для данного набора предметов , у каждого из которых есть определенный вес и ценность , в задаче о рюкзаке необходимо найти подмножество предметов с максимальной ценностью, которое можно упаковать в рюкзак так, чтобы общий вес был меньше вместимости . Математическая формулировка:

где равно , если предмет помещается в рюкзак, и в противном случае.

В этом руководстве одноцелевая задача о рюкзаке расширяется путем добавления еще одной цели: учитывая оценку желательности , нужно максимизировать общую желательность предметов в рюкзаке. В результате математическая формулировка теперь будет такой:

Данные

Данные для этого примера были взяты из репозитория vOptGeneric и изначально созданы пользователем @xgandibleux.

profit = [77, 94, 71, 63, 96, 82, 85, 75, 72, 91, 99, 63, 84, 87, 79, 94, 90]
desire = [65, 90, 90, 77, 95, 84, 70, 94, 66, 92, 74, 97, 60, 60, 65, 97, 93]
weight = [80, 87, 68, 72, 66, 77, 99, 85, 70, 93, 98, 72, 100, 89, 67, 86, 91]
capacity = 900
N = length(profit)

Сравнение вместимости с общим весом всех предметов:

capacity / sum(weight)

0.6428571428571429

показывает, что можно взять примерно 64 % предметов.

Построив график, можно увидеть, что предметы различаются по ценности и желательности. Некоторые предметы имеют высокую ценность и желательность, другие — низкую ценность и высокую желательность (и наоборот).

Plots.scatter(
    profit,
    desire;
    xlabel = "Profit",
    ylabel = "Desire",
    legend = false,
)

Цель двухцелевой задачи о рюкзаке — выбрать подмножество с максимизацией обеих целевых функций.

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

Формулировка в JuMP является прямым переводом математической формулировки:

model = Model()
@variable(model, x[1:N], Bin)
@constraint(model, sum(weight[i] * x[i] for i in 1:N) <= capacity)
@expression(model, profit_expr, sum(profit[i] * x[i] for i in 1:N))
@expression(model, desire_expr, sum(desire[i] * x[i] for i in 1:N))
@objective(model, Max, [profit_expr, desire_expr])

2-element Vector{AffExpr}:
 77 x[1] + 94 x[2] + 71 x[3] + 63 x[4] + 96 x[5] + 82 x[6] + 85 x[7] + 75 x[8] + 72 x[9] + 91 x[10] + 99 x[11] + 63 x[12] + 84 x[13] + 87 x[14] + 79 x[15] + 94 x[16] + 90 x[17]
 65 x[1] + 90 x[2] + 90 x[3] + 77 x[4] + 95 x[5] + 84 x[6] + 70 x[7] + 94 x[8] + 66 x[9] + 92 x[10] + 74 x[11] + 97 x[12] + 60 x[13] + 60 x[14] + 65 x[15] + 97 x[16] + 93 x[17]

Обратите внимание, что многоцелевая программа формируется путем передачи вектора скалярных целевых функций.

Решение

Для решения этой модели нужен оптимизатор, поддерживающий многоцелевые линейные программы. Один из вариантов — пакет MultiObjectiveAlgorithms.jl.

set_optimizer(model, () -> MOA.Optimizer(HiGHS.Optimizer))
set_silent(model)

MultiObjectiveAlgorithms.jl поддерживает множество различных алгоритмов для решения задач многоцелевой оптимизации. Одним из вариантов является метод эпсилон-ограничений:

set_attribute(model, MOA.Algorithm(), MOA.EpsilonConstraint())

Давайте решим задачу и посмотрим на решение:

optimize!(model)
@assert termination_status(model) == OPTIMAL
solution_summary(model)

* Solver : MOA[algorithm=MultiObjectiveAlgorithms.EpsilonConstraint, optimizer=HiGHS]

* Status
  Result count       : 9
  Termination status : OPTIMAL
  Message from the solver:
  "Solve complete. Found 9 solution(s)"

* Candidate solution (result #1)
  Primal status      : FEASIBLE_POINT
  Dual status        : NO_SOLUTION
  Objective value    : [9.18000e+02,9.83000e+02]
  Objective bound    : [9.55000e+02,9.83000e+02]
  Relative gap       : 0.00000e+00

* Work counters
  Solve time (sec)   : 2.20082e-01
  Simplex iterations : 0
  Barrier iterations : 0
  Node count         : 0

Доступно 9 решений. Чтобы узнать, сколько решений доступно, можно также воспользоваться методом result_count:

result_count(model)

Доступ к нескольким решениям

Для доступа к девяти различным решениям модели используйте именованный аргумент result методов solution_summary, value и objective_value:

solution_summary(model; result = 5)

* Solver : MOA[algorithm=MultiObjectiveAlgorithms.EpsilonConstraint, optimizer=HiGHS]

* Status
  Result count       : 9
  Termination status : OPTIMAL

* Candidate solution (result #5)
  Primal status      : FEASIBLE_POINT
  Dual status        : NO_SOLUTION
  Objective value    : [9.36000e+02,9.42000e+02]

@assert primal_status(model; result = 5) == FEASIBLE_POINT

@assert is_solved_and_feasible(model; result = 5)

objective_value(model; result = 5)

2-element Vector{Float64}:
 936.0
 942.0

Обратите внимание: поскольку задан вектор из двух целевых функций, целевое значение представляет собой вектор из двух элементов. Можно также запросить значение каждой целевой функции по отдельности:

value(profit_expr; result = 5)

936.0

Визуализация целевого пространства

В отличие от задач одноцелевой оптимизации, задачи многоцелевой оптимизации не имеют единственного оптимального решения. Вместо этого возвращаемые решения представляют собой возможные компромиссы между двумя целями, между которыми пользователь должен сделать выбор. Распространенный способ визуализации — построение графика значений целевых функций каждого из решений:

plot = Plots.scatter(
    [value(profit_expr; result = i) for i in 1:result_count(model)],
    [value(desire_expr; result = i) for i in 1:result_count(model)];
    xlabel = "Profit",
    ylabel = "Desire",
    title = "Objective space",
    label = "",
    xlims = (915, 960),
)
for i in 1:result_count(model)
    y = objective_value(model; result = i)
    Plots.annotate!(y[1] - 1, y[2], (i, 10))
end
ideal_point = objective_bound(model)
Plots.scatter!([ideal_point[1]], [ideal_point[2]]; label = "Ideal point")

Визуализация целевого пространства позволяет пользователю выбрать решение исходя из личных предпочтений. Например, результат #7 близок к максимальному значению ценности, но обеспечивает значительно большую желательность по сравнению с решениями #8 и #9.

В решении #7 выбираются следующие предметы:

items_chosen = [i for i in 1:N if value(x[i]; result = 7) > 0.9]

11-element Vector{Int64}:
  1
  2
  3
  5
  6
  8
 10
 11
 15
 16
 17

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

В MultiObjectiveAlgorithms.jl реализован и ряд других алгоритмов. Попробуйте решить ту же задачу с помощью MOA.Dichotomy(). Будет ли найдено такое же решение?