Генерация столбцов

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

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

using JuMP
import DataFrames
import HiGHS
import Plots
import SparseArrays

Вводная информация

Задача раскроя заключается в разрезании больших рулонов бумаги на более мелкие части.

Обозначим множество возможных размеров кусков, на которые можно разрезать рулон, как . У каждого куска есть ширина и спрос . Ширина большого рулона равна .

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

В руководстве будут использоваться следующие данные:

struct Piece
    w::Float64
    d::Int
end

struct Data
    pieces::Vector{Piece}
    W::Float64
end

function Base.show(io::IO, d::Data)
    println(io, "Data for the cutting stock problem:")
    println(io, "  W = $(d.W)")
    println(io, "with pieces:")
    println(io, "   i   w_i d_i")
    println(io, "  ------------")
    for (i, p) in enumerate(d.pieces)
        println(io, lpad(i, 4), " ", lpad(p.w, 5), " ", lpad(p.d, 3))
    end
    return
end

function get_data()
    data = [
        75.0 38
        75.0 44
        75.0 30
        75.0 41
        75.0 36
        53.8 33
        53.0 36
        51.0 41
        50.2 35
        32.2 37
        30.8 44
        29.8 49
        20.1 37
        16.2 36
        14.5 42
        11.0 33
        8.6 47
        8.2 35
        6.6 49
        5.1 42
    ]
    return Data([Piece(data[i, 1], data[i, 2]) for i in axes(data, 1)], 100.0)
end

data = get_data()

Data for the cutting stock problem:
  W = 100.0
with pieces:
   i   w_i d_i
  ------------
   1  75.0  38
   2  75.0  44
   3  75.0  30
   4  75.0  41
   5  75.0  36
   6  53.8  33
   7  53.0  36
   8  51.0  41
   9  50.2  35
  10  32.2  37
  11  30.8  44
  12  29.8  49
  13  20.1  37
  14  16.2  36
  15  14.5  42
  16  11.0  33
  17   8.6  47
  18   8.2  35
  19   6.6  49
  20   5.1  42

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

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

— это двоичная переменная, которая указывает, используется ли рулон , а в подсчитывается, сколько кусков размера было отрезано от рулона .

Таким образом, наша частично целочисленная линейная программа имеет следующий вид:

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

Формулировка этой модели на JuMP:

I = length(data.pieces)
J = 1_000  # Некоторое большое число
model = Model(HiGHS.Optimizer)
set_silent(model)
@variable(model, x[1:I, 1:J] >= 0, Int)
@variable(model, y[1:J], Bin)
@objective(model, Min, sum(y))
@constraint(model, [i in 1:I], sum(x[i, :]) >= data.pieces[i].d)
@constraint(
    model,
    [j in 1:J],
    sum(data.pieces[i].w * x[i, j] for i in 1:I) <= data.W * y[j],
);

К сожалению, решить эту формулировку для реалистичных случаев невозможно, так как это займет очень много времени. (Попробуйте снять ограничение по времени.)

set_time_limit_sec(model, 5.0)
optimize!(model)
solution_summary(model)

* Solver : HiGHS

* Status
  Result count       : 1
  Termination status : TIME_LIMIT
  Message from the solver:
  "kHighsModelStatusTimeLimit"

* Candidate solution (result #1)
  Primal status      : FEASIBLE_POINT
  Dual status        : NO_SOLUTION
  Objective value    : 4.11000e+02
  Objective bound    : 2.93000e+02
  Relative gap       : 2.87105e-01

* Work counters
  Solve time (sec)   : 5.05693e+00
  Simplex iterations : 17086
  Barrier iterations : -1
  Node count         : 0

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

Теория генерации столбцов

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

Например, если рассматривать только рулон , то допустимым решением будет:

(1 штука куска № 1)
(1 штука куска № 13)
Для всех остальных

Еще одно решение:

(19 штук куска № 20)
Для всех остальных

Такие шаблоны раскроя, как и , недопустимы, так как совокупная длина больше .

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

К сожалению, таких шаблонов будет очень много, поэтому часто бывает практически невозможно перечислить все столбцы .

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

Выбор исходного множества шаблонов

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

patterns = map(1:I) do i
    n_pieces = floor(Int, data.W / data.pieces[i].w)
    return SparseArrays.sparsevec([i], [n_pieces], I)
end

20-element Vector{SparseVector{Int64, Int64}}:
   [1]  =  1
   [2]  =  1
   [3]  =  1
   [4]  =  1
   [5]  =  1
   [6]  =  1
   [7]  =  1
   [8]  =  1
   [9]  =  1
   [10]  =  3
   [11]  =  3
   [12]  =  3
   [13]  =  4
   [14]  =  6
   [15]  =  6
   [16]  =  9
   [17]  =  11
   [18]  =  12
   [19]  =  15
   [20]  =  19

Визуализировать шаблоны можно так:

"""
    cutting_locations(data::Data, pattern::SparseArrays.SparseVector)


A function which returns a vector of the locations along the roll at which to
cut in order to produce pattern `pattern`.
"""

function cutting_locations(data::Data, pattern::SparseArrays.SparseVector)
    locations = Float64[]
    offset = 0.0
    for (i, c) in zip(SparseArrays.findnz(pattern)...)
        for _ in 1:c
            offset += data.pieces[i].w
            push!(locations, offset)
        end
    end
    return locations
end

function plot_patterns(data::Data, patterns)
    plot = Plots.bar(;
        xlims = (0, length(patterns) + 1),
        ylims = (0, data.W),
        xlabel = "Pattern",
        ylabel = "Roll length",
    )
    for (i, p) in enumerate(patterns)
        locations = cutting_locations(data, p)
        Plots.bar!(
            plot,
            fill(i, length(locations)),
            reverse(locations);
            bar_width = 0.6,
            label = false,
            color = "#90caf9",
        )
    end
    return plot
end

plot_patterns(data, patterns)

Базовая задача

Используя начальное множество шаблонов, можно создать и оптимизировать базовую модель:

model = Model(HiGHS.Optimizer)
set_silent(model)
@variable(model, x[1:length(patterns)] >= 0, Int)
@objective(model, Min, sum(x))
@constraint(model, demand[i in 1:I], patterns[i]' * x >= data.pieces[i].d)
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        : NO_SOLUTION
  Objective value    : 4.21000e+02
  Objective bound    : 4.21000e+02
  Relative gap       : 0.00000e+00

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

В этом решении требуется 421 рулон. Это решение неоптимальное, так как модель включает в себя не все множество возможных шаблонов.

Как найти новый столбец, позволяющий улучшить решение?

Выбор новых столбцов

При генерации столбцов новый столбец выбирается путем ослабления ограничения целочисленности для и анализа двойственной переменной , связанной с ограничением спроса .

Например, двойственное решение для demand[13] будет таким:

unset_integer.(x)
optimize!(model)
@assert is_solved_and_feasible(model; dual = true)
π_13 = dual(demand[13])

0.25

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

Если эта задача, называемая задачей ценообразования, имеет целевое значение больше , то мы оцениваем, что добавление y в качестве коэффициентов нового столбца уменьшит целевое значение больше, чем на стоимость дополнительного рулона.

Вот код для решения задачи ценообразования:

function solve_pricing(data::Data, π::Vector{Float64})
    I = length(π)
    model = Model(HiGHS.Optimizer)
    set_silent(model)
    @variable(model, y[1:I] >= 0, Int)
    @constraint(model, sum(data.pieces[i].w * y[i] for i in 1:I) <= data.W)
    @objective(model, Max, sum(π[i] * y[i] for i in 1:I))
    optimize!(model)
    @assert is_solved_and_feasible(model)
    number_of_rolls_saved = objective_value(model)
    if number_of_rolls_saved > 1 + 1e-8
        # Выгода от шаблона больше, чем стоимость нового рулона плюс некоторый
        # допуск
        return SparseArrays.sparse(round.(Int, value.(y)))
    end
    return nothing
end

solve_pricing (generic function with 1 method)

Если решить задачу ценообразования с искусственным двойственным вектором:

solve_pricing(data, [1.0 / i for i in 1:I])

20-element SparseVector{Int64, Int64} with 3 stored entries:
  [1 ]  =  1
  [17]  =  1
  [20]  =  3

решение будет рулон с 1 штукой куска № 1, 1 штукой куска № 17 и 3 штуками куска № 20.

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

solve_pricing(data, zeros(I))

Итеративный алгоритм

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

while true
    # Решаем линейное ослабление
    optimize!(model)
    @assert is_solved_and_feasible(model; dual = true)
    # Получаем новый двойственный вектор
    π = dual.(demand)
    # Решаем задачу ценообразования
    new_pattern = solve_pricing(data, π)
    # Останавливаем итерации, если нет новых шаблонов
    if new_pattern === nothing
        @info "No new patterns, terminating the algorithm."
        break
    end
    push!(patterns, new_pattern)
    # Создаем столбец
    push!(x, @variable(model, lower_bound = 0))
    # Обновляем коэффициент целевой функции для нового столбца
    set_objective_coefficient(model, x[end], 1.0)
    # Обновляем ненулевые элементы в матрице коэффициентов
    for (i, count) in zip(SparseArrays.findnz(new_pattern)...)
        set_normalized_coefficient(demand[i], x[end], count)
    end
    println("Found new pattern. Total patterns = $(length(patterns))")
end

Found new pattern. Total patterns = 21
Found new pattern. Total patterns = 22
Found new pattern. Total patterns = 23
Found new pattern. Total patterns = 24
Found new pattern. Total patterns = 25
Found new pattern. Total patterns = 26
Found new pattern. Total patterns = 27
Found new pattern. Total patterns = 28
Found new pattern. Total patterns = 29
Found new pattern. Total patterns = 30
Found new pattern. Total patterns = 31
Found new pattern. Total patterns = 32
Found new pattern. Total patterns = 33
Found new pattern. Total patterns = 34
Found new pattern. Total patterns = 35
[ Info: No new patterns, terminating the algorithm.

Найдено множество новых шаблонов. Вот шаблон 21:

patterns[21]

20-element SparseVector{Int64, Int64} with 3 stored entries:
  [9 ]  =  1
  [13]  =  2
  [17]  =  1

Давайте теперь посмотрим на шаблоны:

plot_patterns(data, patterns)

Изучение решения

Посмотрим, сколько штук каждого столбца требуется:

solution = DataFrames.DataFrame([
    (pattern = p, rolls = value(x_p)) for (p, x_p) in enumerate(x)
])
filter!(row -> row.rolls > 0, solution)

16×2 DataFrame

Row	pattern	rolls
	Int64	Float64
1	1	38.0
2	2	44.0
3	3	30.0
4	21	0.5
5	22	10.2
6	23	14.65
7	24	23.1
8	25	11.25
9	26	21.35
10	28	4.3
11	29	19.55
12	30	11.25
13	31	17.45
14	33	36.0
15	34	11.4
16	35	41.0

Так как мы решили линейную программу, некоторые столбцы имеют дробные решения. Мы можем создать целочисленное допустимое решение, округлив количества. Потребуется 341 рулон:

sum(ceil.(Int, solution.rolls))

В качестве альтернативы можно повторно ввести ограничения целочисленности и заново решить задачу:

set_integer.(x)
optimize!(model)
@assert is_solved_and_feasible(model)
solution = DataFrames.DataFrame([
    (pattern = p, rolls = value(x_p)) for (p, x_p) in enumerate(x)
])
filter!(row -> row.rolls > 0, solution)

16×2 DataFrame

Row	pattern	rolls
	Int64	Float64
1	1	38.0
2	2	44.0
3	3	30.0
4	21	1.0
5	22	9.0
6	23	19.0
7	24	19.0
8	25	13.0
9	26	17.0
10	28	2.0
11	29	19.0
12	30	13.0
13	31	18.0
14	33	36.0
15	34	15.0
16	35	41.0

Теперь требуется 334 рулона:

sum(solution.rolls)

333.99999999999994

Обратите внимание, что, возможно, это не глобальный минимум, так как мы не добавляем новые столбцы во время решения частично целочисленной задачи model (получился бы алгоритм, известный как ветвь и цена). Тем не менее алгоритм генерации столбцов обычно находит хорошие целочисленные допустимые решения для трудноразрешимых иными способами задач оптимизации.

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

Наша целевая функция служит для минимизации общего количества рулонов. Какова общая величина отходов? Как она соотносится с общим спросом?
Написание алгоритма оптимизации — это только часть проблемы. Как лучше всего донести решение до заинтересованных лиц?