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

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

Введение

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

Расположение объектов

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

производственных предприятий;
складов хранения продукции;
точек продаж.

Эти объекты расположены в отдельных целочисленных точках сетки от до в направлениях и . Чтобы местоположения объектов не совпадали, необходимо выполнение условия: .

В данном примере примем , , , .

Производство и распределение

Существует P видов продукции, производимых предприятиями. Примем .

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

Ограничения на производственную мощность предприятий и ёмкость складов:

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

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

Затраты

Стоимость перевозки продуктов с завода на склад и со склада в точку продаж зависит от расстояния между объектами и от конкретного продукта. Расстояние в этом примере — это расстояние по координатной сетке. Оно равно сумме абсолютных разниц координат и .

Стоимость производства единицы продукта на заводе составляет: или .

Задача оптимизации

Имея набор местоположений объектов, а также ограничения на производственную мощность, ёмкость складов и спрос, найдём:

Уровень производства каждого вида продукции на каждом предприятии;
Схему движения продукции с производственных предприятий на склады;
Схему движения продукции со складских объектов в точки продаж.

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

Переменные и уравнения для задачи оптимизации

Управляющие переменные, то есть те, которые можно изменять в процессе оптимизации:

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

Целевая функция для минимизации:

Ограничения целевой функции:

(производственная мощность предприятия).
(удовлетворение спроса на продукцию).
(вместимость складских объектов).
(каждая торговая точка связана с одним складским объектом).
(значение произведённой продукции является неотрицательной величиной).
(логическое значение )

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

Генерируем исходные данные для задачи

Сгенерируем случайное местоположение объектов:

using Random, Plots, JuMP, HiGHS
Random.seed!(1)

N = 25
N2 = N ^ 2

f = 0.05  
w = 0.05    
s = 0.1   

# Расчёт количества объектов
F = floor(Int, f * N2)  # количество заводов
W = floor(Int, w * N2)  # количество складов
S = floor(Int, s * N2)  # количество точек продаж

# Уникальные местоположения объектов
координаты = randperm(Random.GLOBAL_RNG, N2)[1:(F+W+S)]

# Преобразование линейных индексов в координаты
x_коорд = similar(координаты, Int)
y_коорд = similar(координаты, Int)

for (i, idx) in enumerate(координаты)
    x_коорд[i] = ((idx - 1) % N) + 1
    y_коорд[i] = ((idx - 1) ÷ N) + 1
end

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

p1 = plot(legend=:outerbottom, xlims=(0, N+0.5), ylims=(0, N+0.5), xticks=0:1:N, yticks=0:1:N, grid=false)
scatter!(x_коорд[1:F], y_коорд[1:F], marker=:square, markerstrokewidth=1, markercolor=:white, markerstrokecolor=:red, label="Производственные предприятия")
scatter!(x_коорд[F+1:F+W], y_коорд[F+1:F+W], marker=:utriangle, markerstrokewidth=1, markercolor=:white, markerstrokecolor=:green, label="Складские объекты")
scatter!(x_коорд[F+W+1:F+W+S], y_коорд[F+W+1:F+W+S], marker=:circle, markerstrokewidth=1, markercolor=:white, markerstrokecolor=:blue, label="Торговые точки")
display(p1)

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

P = 20;  # количество видов продукции

# Производственные затраты
затраты = 80 .* rand(F, P) .+ 20;

# Производственная мощность заводов
мощ = 1000 .* rand(F, P) .+ 500;

#  Производственная ёмкость складов
ёмк = P * 400 .* rand(W) .+ P * 400;

# Коэффициент оборачиваемости продукта
оборот = 2 .* rand(P) .+ 1;

# Стоимость транспортировки продукта
трансп = 5 .* rand(P) .+ 5;

# Спрос на продукцию по торговым точкам
d = 300 .* rand(S, P) .+ 200;

Создадим матрицы распределений расстояний:

# Матрица распределения расстояний между заводами и складами
завод_склад = zeros(F, W)
for ii in 1:F
    for jj in 1:W
        завод_склад[ii, jj] = abs(x_коорд[ii] - x_коорд[F + jj]) + abs(y_коорд[ii] - y_коорд[F + jj])
    end
end

# Матрица распределения расстояний между торговыми точками и складами
склад_торг = zeros(S, W)
for ii in 1:S
    for jj in 1:W
        склад_торг[ii, jj] = abs(x_коорд[F + W + ii] - x_коорд[F + jj]) + abs(y_коорд[F + W + ii] - y_коорд[F + jj]);
    end
end

Создадим задачу оптимизации, переменные и ограничения:

# Создание задачи оптимизации
задача = Model(HiGHS.Optimizer);

# Переменная x: поток продуктов от фабрик к складам
@variable(задача, x[1:P, 1:F, 1:W] >= 0);

# Переменная y: бинарная переменная назначения складов точкам продаж
@variable(задача, 0 <= y[1:S, 1:W] <= 1, Int);

Ограничение производственной мощности:

@constraint(задача, огр_мощ[p in 1:P, f in 1:F], sum(x[p, f, w] for w in 1:W) <= мощ[f, p]);

Ограничение на удовлетворение спроса:

@constraint(задача, огр_спрос[p in 1:P, w in 1:W], sum(x[p, f, w] for f in 1:F) == sum(d[s, p] * y[s, w] for s in 1:S));

Ограничение на вместимость складских объектов:

@constraint(задача, огр_склад_ёмк[w in 1:W], sum((1.0 / оборот[p]) * sum(d[s, p] * y[s, w] for s in 1:S) for p in 1:P) <= ёмк[w]);

Ограничение на назначение складов точкам продаж:

@constraint(задача, огр_торг_склад[s in 1:S], sum(y[s, w] for w in 1:W) == 1);

Первая часть целевой функции — производственные затраты:

ЦФ1 = sum(sum(x, dims=3) .* затраты');

Вторая часть целевой функции — транспортные затраты от фабрик к складам:

ЦФ2 = 0.0
for p in 1:P
    ЦФ2 = ЦФ2 + трансп[p] * sum(sum(x[p, :, :] .* завод_склад));
end

Третья часть целевой функции — сумма транспортных расходов со складов до торговых точек:

r = sum(склад_торг .* y, dims=2);
v = d * трансп;
ЦФ3 = sum(v .* r);

Полный вид целевой функции:

@objective(задача, Min, ЦФ1 + ЦФ2 + ЦФ3);

Настройка параметров решателя:

set_optimizer_attribute(задача, "output_flag", true)  
set_optimizer_attribute(задача, "log_to_console", true)

Решение оптимизационной задачи:

optimize!(задача)
решение = value.(x);
ЦФ = objective_value(задача);
статус = termination_status(задача);
вывод = raw_status(задача);

Running HiGHS 1.11.0 (git hash: 364c83a51e): Copyright (c) 2025 HiGHS under MIT licence terms
MIP  has 1333 rows; 21142 cols; 80724 nonzeros; 1922 integer variables (1922 binary)
Coefficient ranges:
  Matrix [1e+00, 4e+03]
  Cost   [3e+01, 2e+06]
  Bound  [1e+00, 1e+00]
  RHS    [1e+00, 2e+04]
Presolving model
1333 rows, 21142 cols, 80724 nonzeros  0s
1333 rows, 21142 cols, 80724 nonzeros  0s

Solving MIP model with:
   1333 rows
   21142 cols (1922 binary, 0 integer, 0 implied int., 19220 continuous, 0 domain fixed)
   80724 nonzeros

Src: B => Branching; C => Central rounding; F => Feasibility pump; J => Feasibility jump;
     H => Heuristic; L => Sub-MIP; P => Empty MIP; R => Randomized rounding; Z => ZI Round;
     I => Shifting; S => Solve LP; T => Evaluate node; U => Unbounded; X => User solution;
     z => Trivial zero; l => Trivial lower; u => Trivial upper; p => Trivial point

        Nodes      |    B&B Tree     |            Objective Bounds              |  Dynamic Constraints |       Work      
Src  Proc. InQueue |  Leaves   Expl. | BestBound       BestSol              Gap |   Cuts   InLp Confl. | LpIters     Time

         0       0         0   0.00%   0               inf                  inf        0      0      0         0     1.0s
 R       0       0         0   0.00%   46513964.721    50547653.24682     7.98%        0      0      0      2284     1.5s
 C       0       0         0   0.00%   46798746.52812  49053386.94956     4.60%     1202     87      0      2826     2.7s
 L       0       0         0   0.00%   46803901.39706  46804109.73308     0.00%     1491    112      0      2892     8.2s
         1       0         1 100.00%   46803901.39706  46804109.73308     0.00%     1491    112      0      3959     8.3s

Solving report
  Status            Optimal
  Primal bound      46804109.7331
  Dual bound        46803901.3971
  Gap               0.000445% (tolerance: 0.01%)
  P-D integral      0.345085736919
  Solution status   feasible
                    46804109.7331 (objective)
                    0 (bound viol.)
                    0 (int. viol.)
                    0 (row viol.)
  Timing            8.28 (total)
                    0.00 (presolve)
                    0.00 (solve)
                    0.00 (postsolve)
  Max sub-MIP depth 2
  Nodes             1
  Repair LPs        0 (0 feasible; 0 iterations)
  LP iterations     3959 (total)
                    0 (strong br.)
                    608 (separation)
                    1067 (heuristics)

# Если решение не найдено
if !has_values(задача)
    println("Решение не найдено.")
    return 
end

println("Статус: ", статус)

Статус: OPTIMAL

Округлим значение до целого:

решение_y = Int.(round.(value.(y)));

Ответим на вопрос — Сколько торговых точек обслуживает каждый склад? Обратите внимание, что в данном случае некоторым складам соответствует число 0. В найденном оптимальном решении эти склады не задействованы.

торг = sum(решение_y, dims=1);
println("\nТочки продаж: ", торг)

Точки продаж: [0 1 3 3 2 3 3 3 2 1 3 2 2 3 4 0 2 3 2 2 3 2 0 1 0 4 1 4 0 1 2]

Построим схему связи между каждой торговой точкой и соответствующим ей складом:

# Создание графика с точками
p2 = plot(legend=:outerbottom, xlims=(0, N+0.5), ylims=(0, N+0.95), xticks=0:1:N, yticks=0:1:N, grid=false, title = "Оптимальная схема распределения продукции", titlefontsize = 10)
scatter!(x_коорд[1:F], y_коорд[1:F], marker=:square, markerstrokewidth=1, markercolor=:white, markerstrokecolor=:red, label="Производственные предприятия")
scatter!(x_коорд[F+1:F+W], y_коорд[F+1:F+W], marker=:utriangle, markerstrokewidth=1, markercolor=:white, markerstrokecolor=:green, label="Складские объекты")
scatter!(x_коорд[F+W+1:F+W+S], y_коорд[F+W+1:F+W+S], marker=:circle, markerstrokewidth=1, markercolor=:white, markerstrokecolor=:blue, label="Торговые точки")

for ii in 1:S
    jj_list = findall(решение_y[ii, :] .== 1)
    
    for jj in jj_list  
        x_торг = x_коорд[F+W+ii]
        y_торг = y_коорд[F+W+ii]
        x_склад = x_коорд[F+jj]
        y_склад = y_коорд[F+jj]
        
        if rand() < 0.5 
            plot!(p2, [x_торг, x_торг, x_склад], [y_торг, y_склад, y_склад], line=:dash, color=:green, label="")
        else 
            plot!(p2, [x_торг, x_склад, x_склад], [y_торг, y_торг, y_склад], line=:dash, color=:green, label="")
        end
    end
end

display(p2)

Заключение

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

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

Оптимизация распределения продукции