Пример планирования работы фабрик

Это руководство было изначально предоставлено @Crghilardi.

Это руководство представляет собой перевод части 5 из статьи Введение в линейное программирование на Python на язык Julia.

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

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

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

using JuMP
import CSV
import DataFrames
import HiGHS
import StatsPlots

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

Задача планирования работы фабрик предполагает оптимизацию производства товара на фабриках в течение 12 месяцев .

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

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

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

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

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

Данные

Данные, которые понадобятся для этого руководства, содержатся в двух текстовых файлах в репозитории JuMP на GitHub.

В первом файле содержится набор данных по фабрикам A и B с помесячными уровнями производства и затрат. Для документации файл находится по следующему пути:

factories_filename = joinpath(@__DIR__, "factory_schedule_factories.txt");

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

Содержимое файла выглядит так:

print(read(factories_filename, String))

factory month min_production max_production fixed_cost variable_cost
A       1     20000          100000         500        10
A       2     20000          110000         500        11
A       3     20000          120000         500        12
A       4     20000          145000         500        9
A       5     20000          160000         500        8
A       6     20000          140000         500        8
A       7     20000          155000         500        5
A       8     20000          200000         500        7
A       9     20000          210000         500        9
A       10    20000          197000         500        10
A       11    20000          80000          500        8
A       12    20000          150000         500        8
B       1     20000          50000          600        5
B       2     20000          55000          600        4
B       3     20000          60000          600        3
B       4     20000          100000         600        5
B       5     0              0              0          0
B       6     20000          70000          600        6
B       7     20000          60000          600        4
B       8     20000          100000         600        6
B       9     20000          100000         600        8
B       10    20000          100000         600        11
B       11    20000          120000         600        10
B       12    20000          150000         600        12

Для его считывания в Julia мы используем пакеты CSV и DataFrames:

factory_df = CSV.read(
    factories_filename,
    DataFrames.DataFrame;
    delim = ' ',
    ignorerepeated = true,
)

24×6 DataFrame

Row	factory	month	min_production	max_production	fixed_cost	variable_cost
	String1	Int64	Int64	Int64	Int64	Int64
1	A	1	20000	100000	500	10
2	A	2	20000	110000	500	11
3	A	3	20000	120000	500	12
4	A	4	20000	145000	500	9
5	A	5	20000	160000	500	8
6	A	6	20000	140000	500	8
7	A	7	20000	155000	500	5
8	A	8	20000	200000	500	7
9	A	9	20000	210000	500	9
10	A	10	20000	197000	500	10
11	A	11	20000	80000	500	8
12	A	12	20000	150000	500	8
13	B	1	20000	50000	600	5
14	B	2	20000	55000	600	4
15	B	3	20000	60000	600	3
16	B	4	20000	100000	600	5
17	B	5	0	0	0	0
18	B	6	20000	70000	600	6
19	B	7	20000	60000	600	4
20	B	8	20000	100000	600	6
21	B	9	20000	100000	600	8
22	B	10	20000	100000	600	11
23	B	11	20000	120000	600	10
24	B	12	20000	150000	600	12

Второй файл содержит помесячные данные по спросу:

demand_filename = joinpath(@__DIR__, "factory_schedule_demand.txt");

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

demand_df = CSV.read(
    demand_filename,
    DataFrames.DataFrame;
    delim = ' ',
    ignorerepeated = true,
)

12×2 DataFrame

Row	month	demand
	Int64	Int64
1	1	120000
2	2	100000
3	3	130000
4	4	130000
5	5	140000
6	6	130000
7	7	150000
8	8	170000
9	9	200000
10	10	190000
11	11	140000
12	12	100000

Проверка данных

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

function valiate_data(
    demand_df::DataFrames.DataFrame,
    factory_df::DataFrames.DataFrame,
)
    # Минимальный объем производства не должен превышать максимальный.
    @assert all(factory_df.min_production .<= factory_df.max_production)
    # Спрос, минимальный объем производства, постоянные и переменные издержки должны быть
    # неотрицательными.
    @assert all(demand_df.demand .>= 0)
    @assert all(factory_df.min_production .>= 0)
    @assert all(factory_df.fixed_cost .>= 0)
    @assert all(factory_df.variable_cost .>= 0)
    return
end

valiate_data(demand_df, factory_df)

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

Далее необходимо изложить формулировку JuMP в виде кода. Как показано в разделе Design patterns for larger models, всегда лучше программировать модель в виде функции, которая принимает четко определенные входные данные и возвращает четко определенный результат.

function solve_factory_scheduling(
    demand_df::DataFrames.DataFrame,
    factory_df::DataFrames.DataFrame,
)
    # Несмотря на то, что данные были проверены ранее, рекомендуется сделать это и
    # здесь.
    valiate_data(demand_df, factory_df)
    months, factories = unique(factory_df.month), unique(factory_df.factory)
    model = Model(HiGHS.Optimizer)
    set_silent(model)
    @variable(model, status[months, factories], Bin)
    @variable(model, production[months, factories], Int)
    @variable(model, unmet_demand[months] >= 0)
    # Мы используем `eachrow` для перебора строк объекта DataFrame и добавления
    # соответствующих ограничений.
    for r in eachrow(factory_df)
        m, f = r.month, r.factory
        @constraint(model, production[m, f] <= r.max_production * status[m, f])
        @constraint(model, production[m, f] >= r.min_production * status[m, f])
    end
    @constraint(
        model,
        [r in eachrow(demand_df)],
        sum(production[r.month, :]) + unmet_demand[r.month] == r.demand,
    )
    @objective(
        model,
        Min,
        10_000 * sum(unmet_demand) + sum(
            r.fixed_cost * status[r.month, r.factory] +
            r.variable_cost * production[r.month, r.factory] for
            r in eachrow(factory_df)
        )
    )
    optimize!(model)
    @assert is_solved_and_feasible(model)
    schedules = Dict{Symbol,Vector{Float64}}(
        Symbol(f) => value.(production[:, f]) for f in factories
    )
    schedules[:unmet_demand] = value.(unmet_demand)
    return (
        termination_status = termination_status(model),
        cost = objective_value(model),
        # Этот оператор `select` изменяет порядок столбцов в DataFrame.
        schedules = DataFrames.select(
            DataFrames.DataFrame(schedules),
            [:unmet_demand, :A, :B],
        ),
    )
end

solve_factory_scheduling (generic function with 1 method)

Решение

Теперь можно вызвать функцию solve_factory_scheduling со считанными ранее данными.

solution = solve_factory_scheduling(demand_df, factory_df);

Посмотрим, что содержится в solution:

solution.termination_status

OPTIMAL::TerminationStatusCode = 1

solution.cost

1.29064e7

solution.schedules

12×3 DataFrame

Row	unmet_demand	A	B
	Float64	Float64	Float64
1	0.0	70000.0	50000.0
2	0.0	45000.0	55000.0
3	0.0	70000.0	60000.0
4	0.0	30000.0	100000.0
5	0.0	140000.0	-0.0
6	0.0	60000.0	70000.0
7	0.0	90000.0	60000.0
8	0.0	70000.0	100000.0
9	0.0	100000.0	100000.0
10	0.0	190000.0	-0.0
11	0.0	80000.0	60000.0
12	0.0	100000.0	-0.0

Эти планы лучше визуализировать в виде диаграммы:

StatsPlots.groupedbar(
    Matrix(solution.schedules);
    bar_position = :stack,
    labels = ["unmet demand" "A" "B"],
    xlabel = "Month",
    ylabel = "Production",
    legend = :topleft,
    color = ["#20326c" "#4063d8" "#a0b1ec"],
)

Обратите внимание, что неудовлетворенного спроса нет.

Что произойдет при увеличении спроса?

Давайте проведем эксперимент, увеличив спрос на 50 % во все временные периоды:

demand_df.demand .*= 1.5

12-element Vector{Float64}:
 180000.0
 150000.0
 195000.0
 195000.0
 210000.0
 195000.0
 225000.0
 255000.0
 300000.0
 285000.0
 210000.0
 150000.0

Теперь решим задачу:

high_demand_solution = solve_factory_scheduling(demand_df, factory_df);

и визуализируем решение:

StatsPlots.groupedbar(
    Matrix(high_demand_solution.schedules);
    bar_position = :stack,
    labels = ["unmet demand" "A" "B"],
    xlabel = "Month",
    ylabel = "Production",
    legend = :topleft,
    color = ["#20326c" "#4063d8" "#a0b1ec"],
)

Увы, удовлетворить спрос полностью не получается.

Насколько чувствительно решение к изменению переменных издержек?

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

Но сначала вернем спрос на исходный уровень:

demand_df.demand ./= 1.5;

Для нашего эксперимента мы будем масштабировать переменные издержки обеих фабрик на набор значений от 0.0 до 1.5:

scale_factors = 0:0.1:1.5

0.0:0.1:1.5

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

cost = zeros(length(scale_factors), length(scale_factors));

Так как factory_df изменяется на месте, исходные переменные издержки необходимо сохранить в новом столбце:

factory_df[!, :old_variable_cost] = copy(factory_df.variable_cost);

Затем потребуется функция для масштабирования столбца :variable_cost для определенной фабрики на значение scale:

function scale_variable_cost(df, factory, scale)
    rows = df.factory .== factory
    df[rows, :variable_cost] .=
        round.(Int, df[rows, :old_variable_cost] .* scale)
    return
end

scale_variable_cost (generic function with 1 method)

Наш эксперимент — это просто вложенный цикл for, изменяющий A и B и сохраняющий издержки:

for (j, a) in enumerate(scale_factors)
    scale_variable_cost(factory_df, "A", a)
    for (i, b) in enumerate(scale_factors)
        scale_variable_cost(factory_df, "B", b)
        cost[i, j] = solve_factory_scheduling(demand_df, factory_df).cost
    end
end

Визуализируем матрицу издержек:

StatsPlots.contour(
    scale_factors,
    scale_factors,
    cost;
    xlabel = "Scale of factory A",
    ylabel = "Scale of factory B",
)

Какие выводы можно сделать из решения?

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