Unit Commitment for Electrical Power Generation

unit_commitment.ipynb Open In Colab Kaggle Gradient Open In SageMaker Studio Lab Hits

Description: This notebook illustrates the power generation problem using AMPL. The original version featured the Gurobi solver. By default, this notebook uses the HiGHS and CBC solvers.

Major electric power companies around the world utilize mathematical optimization to manage the flow of energy across their electrical grids. In this example, you’ll discover the power of mathematical optimization in addressing a common energy industry problem: unit commitment for electrical power generation. We’ll show you how to figure out the optimal set of power stations to turn on in order to satisfy anticipated power demand over a 24-hour time horizon.

This model is example 15 from the fifth edition of Model Building in Mathematical Programming by H. Paul Williams on pages 270 – 271 and 325 – 326.

This example is at the intermediate level, where we assume that you know Python and the AMPL’s Python API and that you have some knowledge of building mathematical optimization models.

Tags: amplpy, energy, power-generation, unit-commitment

Notebook author: Gyorgy Matyasfalvi <gyorgy@ampl.com>

References:

  1. Electrical power generation 1 notebook

  2. H. Paul Williams, Model Building in Mathematical Programming, fifth edition.

# Install dependencies
%pip install -q amplpy pandas
# Google Colab & Kaggle integration
from amplpy import AMPL, ampl_notebook

ampl = ampl_notebook(
    modules=["highs", "cbc"],  # modules to install
    license_uuid="default",  # license to use
)  # instantiate AMPL object and register magics

Problem Description

In this problem, power generation units are grouped into three distinct types, with different characteristics for each type (power output, cost per megawatt hour, startup cost, etc.). A unit can be on or off, with a startup cost associated with transitioning from off to on, and power output that can fall anywhere between a specified minimum and maximum value when the unit is on. A 24-hour time horizon is divided into 5 discrete time periods, each with an expected total power demand. The model decides which units to turn on, and when, in order to satisfy demand for each time period. The model also captures a reserve requirement, where the selected power plants must be capable of increasing their output, while still respecting their maximum output, in order to cope with the situation where actual demand exceeds predicted demand.

A set of generators is available to satisfy power demand for the following day. Anticipated demand is as follows:

Time Period

Demand (megawatts)

12 am to 6 am

15000

6 am to 9 am

30000

9 am to 3 pm

25000

3 pm to 6 pm

40000

6 pm to 12 am

27000

Generators are grouped into three types, with the following minimum and maximum output for each type (when they are on):

Type

Number available

Minimum output (MW)

Maximum output (MW)

wind

12

850

2000

gas

10

1250

1750

hydro

5

1500

4000

There are costs associated with using a generator: a cost per hour when the generator is on (and generating its minimum output), a cost per megawatt hour above its minimum, and a startup cost for turning a generator on:

Type

Cost per hour (when on)

Cost per MWh above minimum

Startup cost

wind

$$1000$

$$2.00$

$$2000$

gas

$$2600$

$$1.30$

$$1000$

hydro

$$3000$

$$3.00$

$$500$

Generators must meet predicted demand, but they must also have sufficient reserve capacity to be able to cope with the situation where actual demand exceeds predicted demand. For this model, the set of selected generators must be able to produce as much as 115% of predicted demand.

Which generators should be committed to meeting anticipated demand in order to minimize total cost?


Model Formulation

Sets and Indices

$t \in \text{Types}={ \text{wind},; \text{gas},; \text{hydro} }$: Types of generators.

$p \in \text{Periods}={ \text{‘12 am to 6 am’},; \text{‘6 am to 9 am’},; \text{‘9 am to 3 pm’},; \text{‘3 pm to 6 pm’},; \text{‘6 pm to 12 am’} }$: Time periods.

Parameters

$\text{period_hours}_p \in \mathbb{N}^+$: Number of hours in each time period.

$\text{generators}_t \in \mathbb{N}^+$: Number of generators of type $t$.

$\text{demand}_p \in \mathbb{R}^+$: Total power demand for time period $p$.

$\text{maxstart0} \in \mathbb{N}^+$: Number of generators that are on at the beginning of the time horizon (and available in time period 0 without paying a startup cost).

$\text{min_output}_t \in \mathbb{R}^+$: Minimum output for generator type $t$ (when on).

$\text{max_output}_t \in \mathbb{R}^+$: Maximum output for generator type $t$.

$\text{base_cost}_t \in \mathbb{R}^+$: Minimum operating cost (per hour) for a generator of type $t$.

$\text{per_mw_cost}_t \in \mathbb{R}^+$: Cost to generate one additional MW (per hour) for a generator of type $t$.

$\text{startup_cost}_t \in \mathbb{R}^+$: Startup cost for generator of type $t$.

Decision Variables

$\text{ngen}_{t,p} \in \mathbb{N}^+$: Number of generators of type $t$ that are on in time period $p$.

$\text{output}_{t,p} \in \mathbb{R}^+$: Total power output from generators of type $t$ in time period $p$.

$\text{nstart}_{t,p} \in \mathbb{N}^+$: Number of generators of type $t$ to start in time period $p$.

Objective Function

  • Cost: Minimize the cost (in USD) to satisfy the predicted electricity demand.

\begin{equation} \text{Minimize} \quad Z_{on} + Z_{extra} + Z_{startup} \end{equation}

\begin{equation} Z_{on} = \sum_{(t,p) \in \text{Types} \times \text{Periods}}{\text{base_cost}t*\text{ngen}{t,p}} \end{equation} \begin{equation} Z_{extra} = \sum_{(t,p) \in \text{Types} \times \text{Periods}}{\text{per_mw_cost}t*(\text{output}{t,p} - \text{min_load}t}) \end{equation} \begin{equation} Z{startup} = \sum_{(t,p) \in \text{Types} \times \text{Periods}}{\text{startup_cost}t*\text{nstart}{t,p}} \end{equation}

Constraints

  • Available generators: Number of generators used must be less than or equal to the number available.

\begin{equation} \text{ngen}{t,p} \leq \text{generators}{t} \quad \forall (t,p) \in \text{Types} \times \text{Periods} \end{equation}

  • Demand: Total power generated across all generator types must meet anticipated demand for each time period $p$.

\begin{equation} \sum_{t \in \text{Types}}{\text{output}_{t,p}} \geq \text{demand}_p \quad \forall p \in \text{Periods} \end{equation}

  • Min/max generation: Power generation must respect generator min/max values.

\begin{equation} \text{output}_{t,p} \geq \text{min_output}t*\text{ngen}{t,p} \quad \forall (t,p) \in \text{Types} \times \text{Periods} \end{equation}

\begin{equation} \text{output}_{t,p} \leq \text{max_output}t*\text{ngen}{t,p} \quad \forall (t,p) \in \text{Types} \times \text{Periods} \end{equation}

  • Reserve: Selected generators must be able to satisfy demand that is as much as 15% above the prediction.

\begin{equation} \sum_{t \in \text{Types}}{\text{max_output}t*\text{ngen}{t,p}} \geq 1.15 * \text{demand}_p \quad \forall p \in \text{Periods} \end{equation}

  • Startup: Establish relationship between number of active generators and number of startups (use $maxstart0$ for period before the time horizon starts)

\begin{equation} \text{ngen}{t,p} \leq \text{ngen}{t,p-1} + \text{startup}_{t,p} \quad \forall (t,p) \in \text{Types} \times \text{Periods} \end{equation}


Load necessary modules

import pandas as pd

Reset AMPL

Needed to allow for repeated runs

ampl.reset()

Model Development

We first create the sets, parameters, and the variables in AMPL. For each time period, we have: an integer decision variable to capture the number of active generators of each type (ngen), an integer variable to capture the number of generators of each type we must start (nstart), and a continuous decision variable to capture the total power output for each generator type (output).

%%ampl_eval

# Declare indexing sets used for delcaration of model entities
set TYPES;
set PERIODS ordered;

# Declare parameters
param maxstart0;
param reserve_capacity > 1.0; # Ensure that reserve capcity is always above 100%
param generators {TYPES};
param period_hours {PERIODS};
param demand {PERIODS};
param min_load {TYPES};
param max_load {TYPES};
param base_cost {TYPES};
param per_mw_cost {TYPES};
param startup_cost {TYPES};

#Declare variables
var ngen{TYPES, PERIODS} >= 0, integer;
var nstart{TYPES, PERIODS} >= 0, integer;
var output{TYPES, PERIODS} >= 0;

Load data directly from Python data structures using amplpy

We define all the input data of the model and send it to AMPL. It’s good practice to first load the data of all indexing sets into AMPL as it allows us to send data for indexed entities more readily in later stages.

### First load indexing sets ###
# Generator types used for indexing
TYPES = ["gas", "hydro", "wind"]
ampl.set["TYPES"] = TYPES

# Time Periods used for indexing
PERIODS = [
    "12 am to 6 am",
    "6 am to 9 am",
    "9 am to 3 pm",
    "3 pm to 6 pm",
    "6 pm to 12 am",
]
ampl.set["PERIODS"] = PERIODS

### Second load parameter data ###
# Load reserve capacity and maxtart0
reserve_capcity = 1.15
ampl.param["reserve_capacity"] = reserve_capcity

# Load maxstart0
maxstart0 = 5
ampl.param["maxstart0"] = maxstart0

# Data for generators
generators = [12, 10, 5]
period_hours = [6, 3, 6, 3, 6]
demand = [15000, 30000, 25000, 40000, 27000]
min_load = [850, 1250, 1500]
max_load = [2000, 1750, 4000]
base_cost = [1000, 2600, 3000]
per_mw_cost = [2, 1.3, 3]
startup_cost = [2000, 1000, 500]

# Set all params indexed over TYPES at once
# First create Pandas DataFrame
params_indexed_over_TYPES = pd.DataFrame(
    {
        "generators": generators,
        "min_load": min_load,
        "max_load": max_load,
        "base_cost": base_cost,
        "per_mw_cost": per_mw_cost,
        "startup_cost": startup_cost,
    },
    index=TYPES,
)
# Set data using set_data() method
# Since we already populated the indexing sets above,
# here it's enough to invoke 'set_data()' without its second argument (the indexing set argument).
ampl.set_data(params_indexed_over_TYPES)

# Set all params indexed over PERIODS at once
params_indexed_over_PERIODS = pd.DataFrame(
    {"period_hours": period_hours, "demand": demand}, index=PERIODS
)
# Again since we already populated the indexing sets above,
# here it's enough to invoke 'set_data()' without its second argument (the indexing set argument).
ampl.set_data(params_indexed_over_PERIODS)

Next we insert the constraints:

The number of active generators can’t exceed the number of generators.

%%ampl_eval

# Generator count
s.t. numgen {type in TYPES, period in PERIODS}:
             ngen[type, period] <= generators[type];

Total power output for a generator type depends on the number of generators of that type that are active.

%%ampl_eval

# Respect minimum and maximum output per generator type
s.t. min_output {type in TYPES, period in PERIODS}:
                 output[type, period] >= min_load[type] * ngen[type, period];
          
s.t. max_output {type in TYPES, period in PERIODS}:
                 output[type, period] <= max_load[type] * ngen[type, period];

Total output for each time period must meet predicted demand.

%%ampl_eval

# Meet demand
s.t. meet_demand {period in PERIODS}:
                  sum {type in TYPES} output[type, period] >= demand[period];

Selected generators must be able to cope with an excess of demand.

%%ampl_eval

# Provide sufficient reserve capacity
s.t. reserve {period in PERIODS}:
             sum {type in TYPES} max_load[type] * ngen[type, period] >= reserve_capacity * demand[period];

Connect the decision variables that capture active generators with the decision variables that count startups.

%%ampl_eval

# Startup constraints
s.t. startup0 {type in TYPES}:
              ngen[type, first(PERIODS)] <= maxstart0 + nstart[type, first(PERIODS)];
          
s.t. startup {type in TYPES, period in PERIODS: period != first(PERIODS)}:
             ngen[type, period] <= ngen[type, prev(period)] + nstart[type, period];

Objective: minimize total cost. Cost consists of three components: the cost for running active generation units, the cost to generate power beyond the minimum for each unit, and the cost to start up generation units.

%%ampl_eval

# Objective: minimize total cost
minimize total_cost:
         sum {type in TYPES, period in PERIODS} 
             ( base_cost[type] * period_hours[period] * ngen[type, period] 
             +
             per_mw_cost[type] * period_hours[period] * ( output[type, period] - min_load[type] * ngen[type, period] )
             + 
             startup_cost[type] * nstart[type, period] );

Set option solver and solve

  1. Solve with HiGHS:

ampl.set_option("solver", "highs")
ampl.solve()
HiGHS 1.5.1: HiGHS 1.5.1: optimal solution; objective 1002540
38 simplex iterations
1 branching nodes
  1. Solve with CBC:

ampl.set_option("solver", "cbc")
ampl.solve()
cbc 2.10.7: cbc 2.10.7: optimal solution; objective 1002540
7 simplex iterations
7 barrier iterations

Analysis

The anticipated demand for electricity over the 24-hour time window can be met for a total cost of $$1,002,540$. The detailed plan for each time period is as follows.

Unit Commitments

The following table shows the number of generators of each type that are active in each time period in the optimal solution using AMPL’s display command:

%%ampl_eval
option display_1col 0, display_transpose -10, display_width 100;
display ngen;
ngen [*,*]
:     '12 am to 6 am' '6 am to 9 am' '9 am to 3 pm' '3 pm to 6 pm' '6 pm to 12 am'    :=
gas          12              12             12             12             12
hydro         3               8              8              9              9
wind          0               0              0              2              0
;

The following shows the number of generators of each type that must be started in each time period to achieve this plan (recall that the model assumes that up to 5 generators of each type are available at the beginning of the time horizon). Here we use Pandas DataFrame and IPython’s display() function:

start_units = ampl.get_data("nstart").to_pandas()
# Capture the original order b/c unstack will mess-up the original ordering
original_order = start_units.index.get_level_values(1).unique()
start_units = start_units.unstack(level=1)
# Convert multi-index to simple index
start_units.columns = start_units.columns.get_level_values(1)
# Reindex using the original order
start_units = start_units.reindex(columns=original_order)
display(start_units)
12 am to 6 am 6 am to 9 am 9 am to 3 pm 3 pm to 6 pm 6 pm to 12 am
gas 7.0 0.0 0.0 0.0 0.0
hydro 0.0 5.0 0.0 1.0 0.0
wind 0.0 0.0 0.0 2.0 0.0