Оптимизация распределения мощности генераторов

Введение

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

Постановка задачи

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

Вы можете назначать уровень мощности для каждого генератора на каждые получасовые интервалы в течение суток (24 часа, что составляет 48 интервалов). Предположим, что нам известна выручка за мегаватт-час (МВт×ч), получаемая в каждом временном интервале. Цель задачи — определить такое расписание включений генераторов, при котором доходность от их использования будет максимальной.

Присоединим необходимые библиотеки:

#import Pkg
#Pkg.add(["JuMP", "HiGHS"])
using JuMP, HiGHS

Обозначим величины выручки получаемые в каждом временном интервале.

Цены = [
    58.0, 56.5, 55.8, 54.0, 52.1, 52.5, 58.2, 58.0, 56.0, 48.5,
    49.5, 54.5, 53.0, 52.5, 55.0, 55.0, 62.0, 62.0, 62.0, 62.0,
    62.0, 62.0, 62.0, 62.5, 62.5, 62.5, 69.5, 87.0, 97.0, 157.0,
    280.0, 280.0, 280.0, 103.0, 97.0, 69.0, 61.0, 67.0, 82.0, 64.0,
    61.5, 69.0, 62.5, 62.0, 58.5, 62.0, 62.0, 57.5
]

Построим график зависимости стоимости электроэнергии от временных периодов.

p1 = bar(Цены, 
    bar_width=0.5,
    legend=false,
    xlims=(0.5, 48.5),
    xlabel="Периоды",
    ylabel="₽/МВт*ч",
    ylim=(0, 299),
    title="Стоимость электроэнергии",
    color=:blue,
    alpha=0.7)
display(p1)

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

Обозначения и параметры задачи

Задачу планирования можно сформулировать как задачу бинарного целочисленного программирования. Определим индексы i, j, k и бинарный вектор планирования y следующим образом:

периодов = количество временных периодов, в данном случае 48.
i = временной период, .
j = индекс генератора, для данного примера.
y[i,j,k] = 1, если в период i генератор j работает на уровне мощности . Пусть номинальная мощность соответствует , а пиковая — . Генератор выключен, когдаy[i,j,k] = 0.

Определим, когда генератор запускается после простоя. Для этого введем вспомогательную бинарную переменную z[i,j], которая указывает, взимается ли плата за включение генератора j в период i.

z[i,j] = 1, если генератор выключен в период , но включен в период . z[i,j] = 0 в противном случае. То есть, z[i,j] = 1, когда y[i,j,k] = 0 и y[i+1,j,k] = 1.

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

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

Цены[i] — Величины выручки за МВт·ч в интервале ;
генератор[j,k] — Мощность (в МВт), вырабатываемая генератором на уровне ;
топливо[j,k] — Величина топлива, потребляемого генератором на уровне ;
общееТопливо — Доступный объём топлива за сутки;
стоимостьЗапуска — Затраты (в долларах) на запуск генератора после простоя;
ценаТоплива — Стоимость единицы топлива.

Зададим другие параметры:

ценаТоплива = 8
общееТопливо = 6.5e4
периодов = length(Цены)  # 48 периодов
генераторов = 2  # Два генератора
генератор = [71 172; 60 150] # Номинальные и пиковые режимы генераторов
топливо = [400 950; 280 800] # Расход топлива для генераторов в разных режимах работы
стоимостьЗапуска = 1000.0;  # Стоимость запуска генератора после простоя

Расход топлива

Исследуем удельную производительность генераторов в номинальном и пиковом режимах.

производительность = генератор ./ топливо
p2 = bar([0.875, 1.875], [производительность[1,1], производительность[1,2]],
        color = :green,
        label = "Генератор 1",
        alpha = 0.9,
        bar_width = 0.15,
        xlabel = "Режим работы",
        ylabel = "Мощность на единицу топлива",
        title = "Удельная производительность",
        ylim = (0.1, 0.2),
        legend = :outerright)
bar!(p2, [1.125, 2.125], [производительность[2,1], производительность[2,2]],
     color = :cyan,
     label = "Генератор 2",
     alpha = 0.9,
     bar_width = 0.15)
xticks!([1, 2], ["Номинальный", "Пиковый"])
xlims!(0.4, 2.6)
ylims!(0.0, 0.24)
display(p2)

В номинальном режиме удельная производительность Генератора 2 несколько превышает удельную производительность Генератора 1. В пиковом режиме удельная производительность Генератора 1 чуть выше, чем в номинальном, а Генератора 2 несколько ниже, чем в номинальном. При этом в пиковом режиме удельная производительность Генератора 2 чуть выше, чем Генератора 1.

Переменные для решения

Для настройки задачи необходимо представить все данные и ограничения в форме, требуемой для решения. Переменные y[i,j,k] представляют решение задачи, а вспомогательные переменные z[i,j] указывают, взимается ли плата за запуск генератора. Все переменные являются бинарными.

модель = Model(HiGHS.Optimizer)
@variable(модель, y[1:периодов, 1:генераторов, 1:2], Bin) 
@variable(модель, z[1:периодов, 1:генераторов], Bin)

Линейные ограничения

Чтобы гарантировать, что уровень мощности имеет не более одной компоненты, равной 1, зададим линейное ограничение-неравенство.

ограничениеМощности = @constraint(модель, [i=1:периодов, j=1:генераторов], y[i,j,1] + y[i,j,2] <= 1)

Эксплуатационные затраты за период — это стоимость топлива за этот период. Для генератора j, работающего на уровне k, затраты составляют ценаТоплива * Топливо[j,k].

Создадим выражение, которое учитывает всё использованный топливо.

использованноеТопливо = @expression(модель, sum(y[i,j,k] * топливо[j,k] for i in 1:периодов, j in 1:генераторов, k in 1:2))

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

ограничениеТоплива = @constraint(модель, использованноеТопливо <= общееТопливо)

Переменные индикатора запуска генератора

Как обеспечить, чтобы решатель автоматически устанавливал переменные z в соответствии с активными/неактивными периодами переменных ? Напомним, условие, которое должно выполняться: z(i,j) = 1 точно тогда, когда y[i,j,k] = 0 и y[i+1,j,k] = 1.

Заметим, что ( - y[i,j,k] + y[i+1,j,k] ) > 0 как раз в том случае, когда требуется z[i,j] = 1.

Следовательно, включим в формулировку задачи следующие линейные ограничения-неравенства:

( - y[i,j,k] + y[i+1,j,k] ) - z[i,j] <= 0.

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

Создадим вспомогательную переменную , которая представляет y(i+1,j,k) - y(i,j,k). Представим ограничение на запуск генератора через .

w = @expression(модель, [i=1:периодов, j=1:генераторов], 
    if i < периодов
        y[i+1,j,1] - y[i,j,1] + y[i+1,j,2] - y[i,j,2]
    else
        y[1,j,1] - y[i,j,1] + y[1,j,2] - y[i,j,2]
    end)
ограничениеПереключения = @constraint(модель, [i=1:периодов, j=1:генераторов], w[i,j] - z[i,j] <= 0)

Определение целевой функции

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

уровеньГенератора = zeros(периодов, генераторов, 2)
уровеньГенератора[:, 1, 1] .= генератор[1, 1]  
уровеньГенератора[:, 1, 2] .= генератор[1, 2]  
уровеньГенератора[:, 2, 1] .= генератор[2, 1]  
уровеньГенератора[:, 2, 2] .= генератор[2, 2]

уровеньГенератора

Определим выражение для величины выручки.

выручка = @expression(модель, [i=1:периодов, j=1:генераторов, k=1:2],
    Цены[i] * y[i,j,k] * уровеньГенератора[i,j,k])

Определим выражение для величины затрат на топливо.

затратыНаТопливо = использованноеТопливо * ценаТоплива

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

затратыНаЗапуск = z*стоимостьЗапуска

Определим выражение для величины прибыли.

прибыль = @expression(модель, 
    sum(Цены[i] * y[i,j,k] * уровеньГенератора[i,j,k] 
        for i in 1:периодов, j in 1:генераторов, k in 1:2)  - sum(y[i,j,k] * топливо[j,k] * ценаТоплива 
        for i in 1:периодов, j in 1:генераторов, k in 1:2) - sum(z[i,j] 
        for i in 1:периодов, j in 1:генераторов) * стоимостьЗапуска)

Решение задачи

Создадим задачу оптимизации, включим в неё целевую функцию и ограничения.

@objective(модель, Max, прибыль)

Решим задачу.

optimize!(модель)
решение = (y = value.(y), z = value.(z))  
значение = objective_value(модель)

Running HiGHS 1.12.0 (git hash: 755a8e027a): Copyright (c) 2025 HiGHS under MIT licence terms
MIP has 193 rows; 288 cols; 864 nonzeros; 288 integer variables (288 binary)
Coefficient ranges:
  Matrix  [1e+00, 1e+03]
  Cost    [2e+02, 4e+04]
  Bound   [1e+00, 1e+00]
  RHS     [1e+00, 6e+04]
Presolving model
193 rows, 288 cols, 864 nonzeros  0s
193 rows, 288 cols, 864 nonzeros  0s
Presolve reductions: rows 193(-0); columns 288(-0); nonzeros 864(-0) - Not reduced
Objective function is integral with scale 10

Solving MIP model with:
   193 rows
   288 cols (288 binary, 0 integer, 0 implied int., 0 continuous, 0 domain fixed)
   864 nonzeros

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

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

 z       0       0         0   0.00%   inf             -0                 Large        0      0      0         0     0.1s
 J       0       0         0   0.00%   inf             560                Large        0      0      0         0     0.1s
 R       0       0         0   0.00%   514387.540741   504531.6           1.95%        0      0      0       115     0.1s
 L       0       0         0   0.00%   514255.480826   514015.6           0.05%     1045     21     11       241     0.6s

82.3% inactive integer columns, restarting
Model after restart has 24 rows, 37 cols (35 bin., 2 int., 0 impl., 0 cont., 0 dom.fix.), and 104 nonzeros

         0       0         0   0.00%   514255.326734   514015.6           0.05%        7      0      0       581     0.7s
         0       0         0   0.00%   514255.326734   514015.6           0.05%        7      7      2       602     0.7s

2.7% inactive integer columns, restarting
Model after restart has 22 rows, 36 cols (35 bin., 1 int., 0 impl., 0 cont., 0 dom.fix.), and 97 nonzeros

         0       0         0   0.00%   514238.65312    514015.6           0.04%        7      0      0       778     0.8s
         0       0         0   0.00%   514238.65312    514015.6           0.04%        7      7      2       804     0.8s

2.8% inactive integer columns, restarting
Model after restart has 21 rows, 35 cols (34 bin., 1 int., 0 impl., 0 cont., 0 dom.fix.), and 93 nonzeros

         0       0         0   0.00%   514232.533911   514015.6           0.04%        7      0      0       889     0.9s
         0       0         0   0.00%   514232.533911   514015.6           0.04%        7      7      4       905     0.9s
        49       0        20 100.00%   514067.0025     514015.6           0.01%      656      7    354      1824     1.2s

Solving report
  Status            Optimal
  Primal bound      514015.6
  Dual bound        514067
  Gap               0.01% (tolerance: 0.01%)
  P-D integral      4.02219234811
  Solution status   feasible
                    514015.6 (objective)
                    0 (bound viol.)
                    6.17542280566e-15 (int. viol.)
                    0 (row viol.)
  Timing            1.20
  Max sub-MIP depth 4
  Nodes             49
  Repair LPs        0
  LP iterations     1824
                    397 (strong br.)
                    347 (separation)
                    679 (heuristics)

514015.60000000003

Анализ решения

Построим графики оптимального расписания работы генераторов как функцию времени.

p3 = plot(layout = (3,1), size = (800, 600))
bar!(p3[1], 
     решение.y[:,1,1] * генератор[1,1] + решение.y[:,1,2] * генератор[1,2],
     bar_width = 0.5,
     color = :green,
     xlims = (0.5, 48.5),
     ylims = (0, 179),
     ylabel = "МВт*ч",
     title = "Генератор 1 — Оптимальное расписание",
     titlefontweight = :bold,
     legend = false)
bar!(p3[2],
     решение.y[:,2,1] * генератор[2,1] + решение.y[:,2,2] * генератор[2,2],
     bar_width = 0.5,
     color = :cyan,
     ylims = (0, 159),
     xlims = (0.5, 48.5),
     ylabel = "МВт*ч",
     title = "Генератор 2 — Оптимальное расписание",
     titlefontweight = :bold,
     legend = false)
bar!(p3[3],
     Цены,
     bar_width = 0.5,
     ylims = (0, 299),
     xlims = (0.5, 48.5),
     xlabel = "Период",
     ylabel = "₽/МВт*ч",
     title = "Цены на электроэнергию",
     titlefontweight = :bold,
     legend = false
)
display(p3)

В соответствии с оптимальным расписанием, Генератор 1 работает в пиковом режиме, отключается на несколько периодов, затем снова включается и работает в пиковом режиме. Генератор 2 работает работает в течение всех 48 периодов, 2 раза переключаясь с пикового на номинальный режим.

Сравнение с наименьшими затратами на запуск

Зададим меньшее значение стоимости запуска, и снова запустим задачу оптимизации.

стоимостьЗапуска = 100
затратыНаЗапуск = sum(z) * стоимостьЗапуска 
прибыль = @expression(модель, 
    sum(выручка[i,j,k] for i in 1:периодов, j in 1:генераторов, k in 1:2) 
    - затратыНаТопливо 
    - затратыНаЗапуск)
@objective(модель, Max, прибыль)
optimize!(модель)
новоеРешение = (y = value.(y), z = value.(z))
новоеЗначение = objective_value(модель)

MIP has 193 rows; 288 cols; 864 nonzeros; 288 integer variables (288 binary)
Coefficient ranges:
  Matrix  [1e+00, 1e+03]
  Cost    [1e+02, 4e+04]
  Bound   [1e+00, 1e+00]
  RHS     [1e+00, 6e+04]
Presolving model
193 rows, 288 cols, 864 nonzeros  0s
193 rows, 288 cols, 864 nonzeros  0s
Presolve reductions: rows 193(-0); columns 288(-0); nonzeros 864(-0) - Not reduced
Objective function is integral with scale 10

Solving MIP model with:
   193 rows
   288 cols (288 binary, 0 integer, 0 implied int., 0 continuous, 0 domain fixed)
   864 nonzeros

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

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

 z       0       0         0   0.00%   inf             -0                 Large        0      0      0         0     0.0s
 J       0       0         0   0.00%   inf             36560              Large        0      0      0         0     0.0s
 R       0       0         0   0.00%   516033.193684   508237             1.53%        0      0      0       111     0.0s
 C       0       0         0   0.00%   516010.462282   508323.4           1.51%       67      5      0       119     0.0s
 L       0       0         0   0.00%   515939.624512   515785.2           0.03%      662     17      0       196     0.4s

55.9% inactive integer columns, restarting
Model after restart has 36 rows, 53 cols (53 bin., 0 int., 0 impl., 0 cont., 0 dom.fix.), and 136 nonzeros

         0       0         0   0.00%   515939.515044   515785.2           0.03%        7      0      0       506     0.5s
         0       0         0   0.00%   515939.515044   515785.2           0.03%        7      7      2       520     0.5s

5.7% inactive integer columns, restarting
Model after restart has 34 rows, 50 cols (50 bin., 0 int., 0 impl., 0 cont., 0 dom.fix.), and 129 nonzeros

         0       0         0   0.00%   515916.892919   515785.2           0.03%        8      0      0       922     0.9s
         0       0         0   0.00%   515916.892919   515785.2           0.03%        8      8      4       940     0.9s
        44       0        12 100.00%   515835.4        515785.2           0.01%     1547     12    414      2218     1.5s

Solving report
  Status            Optimal
  Primal bound      515785.2
  Dual bound        515835.4
  Gap               0.00973% (tolerance: 0.01%)
  P-D integral      0.0494135065927
  Solution status   feasible
                    515785.2 (objective)
                    0 (bound viol.)
                    8.881784197e-16 (int. viol.)
                    0 (row viol.)
  Timing            1.46
  Max sub-MIP depth 4
  Nodes             44
  Repair LPs        0
  LP iterations     2218
                    538 (strong br.)
                    394 (separation)
                    921 (heuristics)

515785.2

Построим графики оптимального расписания работы генераторов с уменьшенными затратами на запуск.

p4 = plot(layout = (3,1), size = (800, 600))
bar!(p4[1], 
     [новоеРешение.y[i,1,1] * генератор[1,1] + новоеРешение.y[i,1,2] * генератор[1,2] for i in 1:периодов],
     bar_width = 0.5, color = :green, xlims = (0.5, 48.5), ylabel = "МВт*ч", ylims = (0, 179),
     title = "Генератор 1 — Оптимальное расписание", titlefontweight = :bold, legend = false)  
bar!(p4[2],
     [новоеРешение.y[i,2,1] * генератор[2,1] + новоеРешение.y[i,2,2] * генератор[2,2] for i in 1:периодов],
     bar_width = 0.5, color = :cyan, xlims = (0.5, 48.5), ylabel = "МВт*ч", ylims = (0, 159),
     title = "Генератор 2 — Оптимальное расписание", titlefontweight = :bold, legend = false)
bar!(p4[3], Цены,
     bar_width = 0.5, xlims = (0.5, 48.5), xlabel = "Период", ylabel = "₽/МВт*ч", ylims = (0, 299),
     title = "Цены на электроэнергию", titlefontweight = :bold, legend = false)
display(p4)

На графиках представлено оптимальное расписание работы генераторов с учётом сниженной стоимости запуска. Генератор 1 в основном работает в пиковом режиме и имеет два перерыва а работе. Генератор 2 вначале работает в номинальном режиме, затем переключается на пиковый, и только в конце два раза переключается на номинальный режим, и работает непрерывно в течение всех периодов.

Заключение

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

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

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

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

Оптимизация распределения мощности генераторов

Введение

Постановка задачи

Обозначения и параметры задачи

Расход топлива

Переменные для решения

Линейные ограничения

Переменные индикатора запуска генератора

Определение целевой функции

Решение задачи

Анализ решения

Сравнение с наименьшими затратами на запуск

Заключение

Тип

Краткое описание

Связанные материалы

Теги

Категории

Языки

Форматы

Уровень

Статус

Опубликовано

Обновлено

Источник