Engee documentation
Notebook

Optimization of generator power distribution

Introduction

This example demonstrates how to optimally schedule two electric generators to maximize revenue minus costs. Although the example is somewhat simplified, it shows how to account for costs that depend on decision-making time.

Setting the task

The electricity market is characterized by prices that change during the day. If you have generating units at your disposal, you can use this variable price to plan their operation for periods with high prices. Let's assume that you control two generators. Each generator has three modes (off, rated and peak). The specific fuel consumption and the amount of electricity generated at each power level are set for each generator. When the generator is turned off, there is no fuel consumption.

You can set the power level for each generator for every half-hour interval during the day (24 hours, which is 48 intervals). Let's assume that we know the revenue per megawatt hour (MWh) received in each time interval. The purpose of the task is to determine a schedule for turning on generators in which the profitability from their use will be maximum.

Adding the necessary libraries:

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

Let's denote the amount of revenue received in each time interval.

In [ ]:
Prices = [
    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
]

Let's plot the dependence of the cost of electricity on time periods.

In [ ]:
p1 = bar(Prices, 
    bar_width=0.5,
    legend=false,
    xlims=(0.5, 48.5),
    xlabel="Periods",
    ylabel="₽/MWh",
    ylim=(0, 299),
    title="The cost of electricity",
    color=:blue,
    alpha=0.7)
display(p1)

There are costs to start the generator after it is idle. In addition, there is a limit on the maximum fuel consumption per day. This limitation is due to the fact that fuel is purchased a day in advance, and only the pre-purchased volume can be used.

Task designations and parameters

The planning problem can be formulated as a binary integer programming problem. We define the indexes i, j, k and the binary planning vector y as follows:

  • периодов = the number of time periods, in this case 48.

  • i = time period, .

  • j = generator index, for this example.

  • y[i,j,k] = 1 if the generator j is operating at the power level during period I. . Let the rated power correspond to , and the peak — . The generator is turned off wheny[i,j,k] = 0.

Let's determine when the generator starts after idle time. To do this, we will introduce an auxiliary binary variable. z[i,j], which indicates whether a fee is charged for turning on the generator j during period I.

  • z[i,j] = 1 if the generator turned off during the , but included in the period . z[i,j] = 0 otherwise. That is, z[i,j] = 1, when y[i,j,k] = 0 and y[i+1,j,k] = 1.

The values must be set automatically based on the values . This is achieved by using the linear constraint given below.

You will also need task parameters: costs, generation levels for each unit, fuel consumption levels, and available fuel volume.

  • Цены[i] — The amount of revenue per MWh in the range ;

  • генератор[j,k] — The power (in MW) generated by the generator at the level of ;

  • топливо[j,k] — The amount of fuel consumed by the generator at the level of ;

  • общееТопливо — Available amount of fuel per day;

  • стоимостьЗапуска — The cost (in dollars) of starting the generator after downtime;

  • ценаТоплива — The cost of a unit of fuel.

Setting other parameters:

In [ ]:
Fuel price = 8
General fuel = 6.5e4
периодов = length(Цены)  # 48 periods
генераторов = 2  # Two generators
генератор = [71 172; 60 150] # Nominal and peak modes of generators
топливо = [400 950; 280 800] # Fuel consumption for generators in different operating modes
стоимостьЗапуска = 1000.0;  # The cost of starting the generator after downtime

Fuel consumption

We study the specific performance of generators in nominal and peak modes.

In [ ]:
performance = generator ./ fuel
p2 = bar([0.875, 1.875], [performance[1,1], performance[1,2]],
        color = :green,
        label = "Generator 1",
        alpha = 0.9,
        bar_width = 0.15,
        xlabel = "Operating mode",
        ylabel = "Power per unit of fuel",
        title = "Specific productivity",
        ylim = (0.1, 0.2),
        legend = :outerright)
bar!(p2, [1.125, 2.125], [performance[2.1], performance[2.2]],
     color = :cyan,
     label = "Generator 2",
     alpha = 0.9,
     bar_width = 0.15)
xticks!([1, 2], ["Nominal", "Peak"])
xlims!(0.4, 2.6)
ylims!(0.0, 0.24)
display(p2)

In the nominal mode, the specific power of Generator 2 is slightly higher than the specific power of Generator 1. In peak mode, the specific productivity of Generator 1 is slightly higher than in nominal mode, and Generator 2 is slightly lower than in nominal mode. At the same time, in peak mode, the specific productivity of Generator 2 is slightly higher than Generator 1.

Variables to solve

To set up a task, you must present all the data and constraints in the form required for the solution. Variables y[i,j,k] represent the solution to the problem, and auxiliary variables z[i,j] indicate whether there is a charge for starting the generator. All variables are binary.

In [ ]:
model = Model(HiGHS.Optimizer)
@variable(model, y[1:periods, 1:generators, 1:2], Bin) 
@variable(model, z[1:periods, 1:generators], Bin) 

Linear constraints

To ensure that the power level has no more than one component equal to 1, we define a linear constraint-inequality.

In [ ]:
Space constraints = @constraint(model, [i=1:periods, j=1:generators], y[i,j,1] + y[i,j,2] <= 1)

The operating cost for a period is the cost of fuel for that period. For a generator j operating at level k, the costs are ценаТоплива * Топливо[j,k].

Let's create an expression that takes into account all the fuel used.

In [ ]:
Used fuel = @expression(model, sum(y[i,j,k] * fuel[j,k] for i in 1:periods, j in 1:generators, k in 1:2))

The limitation is that the fuel used must not exceed the available volume.

In [ ]:
Fuel Restriction = @constraint(model used by fuel <= General fuel)

Generator start indicator variables

How to ensure that the solver automatically sets the z variables according to the active/inactive periods of the variables ? Recall the condition that must be fulfilled: z(i,j) = 1 exactly when y[i,j,k] = 0 and y[i+1,j,k] = 1.

Note that ( - y[i,j,k] + y[i+1,j,k] ) > 0 just in the case when it is required z[i,j] = 1.

Therefore, we include the following linear constraints in the formulation of the problem-inequalities:

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

In addition, we will include variables to the target function. Because the variables They are included in the objective function, the algorithm seeks to minimize their values, that is, it tries to set them all equal. .

Creating an auxiliary variable , which represents y(i+1,j,k) - y(i,j,k). Let's imagine a restriction on starting the generator via .

In [ ]:
w = @expression(model, [i=1:periods, j=1:generators], 
    if i < periods
        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)
Inclusion constraints = @constraint(model, [i=1:periods, j=1:generators], w[i,j] - z[i,j] <= 0)

Defining the target function

The objective function includes the cost of fuel for generators, the profitability of their operation, and the cost of starting them.

In [ ]:
Generator level = zeros(periods, generators, 2)
The generator level is [:, 1, 1] .= generator[1, 1]  
The generator level is [:, 1, 2] .= generator[1, 2]  
Generator level[:, 2, 1] .= generator[2, 1]  
The generator level is [:, 2, 2] .= generator[2, 2]

Generator Level

Let's define an expression for the amount of revenue.

In [ ]:
revenue = @expression(model, [i=1:periods, j=1:generators, k=1:2],
    Prices[i] * y[i,j,k] * Generator level[i,j,k])

Let's define an expression for the amount of fuel costs.

In [ ]:
Fuel cost = Used fuel * Fuel price

Let's define an expression for the cost of starting the generator.

In [ ]:
Startup Cost = z*Startup cost

Let's define an expression for the amount of profit.

In [ ]:
profit = @expression(model, 
    sum(Prices[i] * y[i,j,k] * Generator level[i,j,k] 
        for i in 1:periods, j in 1:generators, k in 1:2) - sum(y[i,j,k] * fuel[j,k] * Fuel price 
        for i in 1:periods, j in 1:generators, k in 1:2) - sum(z[i,j] 
        for i in 1:periods, j in 1:generators) * Startup cost)

Problem solving

Let's create an optimization task and include the objective function and constraints in it.

In [ ]:
@objective(Model, Max, profit)

Let's solve the problem.

In [ ]:
optimize!(model)
solution = (y = value.(y), z = value.(z))  
value = objective_value(model)                               
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)
Out[0]:
514015.60000000003

Solution analysis

Let's plot the optimal schedule of generators as a function of time.

In [ ]:
p3 = plot(layout = (3,1), size = (800, 600))
bar!(p3[1], 
     decision.y[:,1,1] * generator[1,1] + solution.y[:,1,2] * generator[1,2],
     bar_width = 0.5,
     color = :green,
     xlims = (0.5, 48.5),
     ylims = (0, 179),
     ylabel = "MWh",
     title = "Generator 1 — Optimal schedule",
     titlefontweight = :bold,
     legend = false)
bar!(p3[2],
     decision.y[:,2,1] * generator[2,1] + solution.y[:,2,2] * generator[2,2],
     bar_width = 0.5,
     color = :cyan,
     ylims = (0, 159),
     xlims = (0.5, 48.5),
     ylabel = "MWh",
     title = "Generator 2 — Optimal schedule",
     titlefontweight = :bold,
     legend = false)
bar!(p3[3],
     Prices,
     bar_width = 0.5,
     ylims = (0, 299),
     xlims = (0.5, 48.5),
     xlabel = "Period",
     ylabel = "₽/MWh",
     title = "Electricity prices",
     titlefontweight = :bold,
     legend = false
)
display(p3)

According to the optimal schedule, Generator 1 operates in peak mode, turns off for several periods, then turns on again and operates in peak mode. Generator 2 works for all 48 periods, switching from peak to nominal mode 2 times.

Comparison with the lowest startup cost

Set a lower value for the startup cost, and run the optimization task again.

In [ ]:
Startup Cost = 100
Startup cost = sum(z) * Startup cost 
profit = @expression(model, 
    sum(revenue[i,j,k] for i in 1:periods, j in 1:generators, k in 1:2) 
    - Fuel consumption 
    - Startup cost)
@objective(Model, Max, profit)
optimize!(model)
New generation = (y = value.(y), z = value.(z))
New value = objective_value(model)
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)
Out[0]:
515785.2

Let's plot the optimal operating schedule of generators with reduced startup costs.

In [ ]:
p4 = plot(layout = (3,1), size = (800, 600))
bar!(p4[1], 
     [New generation.y[i,1,1] * generator[1,1] + New generation.y[i,1,2] * generator[1,2] for i in 1:periods],
     bar_width = 0.5, color = :green, xlims = (0.5, 48.5), ylabel = "MWh", ylims = (0, 179),
     title = "Generator 1 — Optimal schedule", titlefontweight = :bold, legend = false)  
bar!(p4[2],
     [New generation.y[i,2,1] * generator[2,1] + New generation.y[i,2,2] * generator[2,2] for i in 1:periods],
     bar_width = 0.5, color = :cyan, xlims = (0.5, 48.5), ylabel = "MWh", ylims = (0, 159),
     title = "Generator 2 — Optimal schedule", titlefontweight = :bold, legend = false)
bar!(p4[3], Prices,
     bar_width = 0.5, xlims = (0.5, 48.5), xlabel = "Period", ylabel = "₽/MWh", ylims = (0, 299),
     title = "Electricity prices", titlefontweight = :bold, legend = false)
display(p4)

The graphs show the optimal operating schedule of the generators, taking into account the reduced cost of starting. Generator 1 mainly operates in peak mode and has two interruptions in operation. Generator 2 initially operates in nominal mode, then switches to peak mode, and only at the end switches to nominal mode twice, and operates continuously during all periods.

Conclusion

The study confirmed the effectiveness of using binary integer programming to optimize the daily schedule of electric generators. The developed model successfully takes into account key factors: changing market prices, discrete power levels, fuel consumption and startup costs.

The results show that the cost structure critically influences the strategy. The high startup cost leads to long continuous work cycles, while its reduction makes it economically advantageous to work more flexibly with multiple inclusions to accurately follow price peaks. The model also quantified the choice of generator modes, taking into account not only their specific efficiency, but also the time-varying cost of energy.

The computational experiment demonstrates the practical applicability of the approach: a medium-scale problem is solved to a close to optimal solution in an acceptable time. This opens up opportunities for the implementation of such models in dispatch support systems, ensuring a balance between maximizing revenue and complying with technological constraints.