Economic Dispatch

Economic Dispatch#

Note

If you have not yet set up Python on your computer, you can execute this tutorial in your browser via Google Colab. Click on the rocket in the top right corner and launch “Colab”. If that doesn’t work download the .ipynb file and import it in Google Colab.

Then install the following packages by executing the following command in a Jupyter cell at the top of the notebook.

!pip install pypsa highspy

From structured data to optimisation#

The design principle of PyPSA is that each component is associated with a set of variables and constraints that will be added to the optimisation model based on the input data stored for the components.

For an hourly electricity market simulation, PyPSA will solve an optimisation problem that looks like this

(1)#\[\begin{equation} \min_{g_{i,s,t}; f_{\ell,t}; g_{i,r,t,\text{charge}}; g_{i,r,t,\text{discharge}}; e_{i,r,t}} \sum_s o_{s} g_{i,s,t} \end{equation}\]

such that

(2)#\[\begin{align} 0 & \leq g_{i,s,t} \leq \hat{g}_{i,s,t} G_{i,s} & \text{generation limits : generator} \\ -F_\ell &\leq f_{\ell,t} \leq F_\ell & \text{transmission limits : line} \\ d_{i,t} &= \sum_s g_{i,s,t} + \sum_r g_{i,r,t,\text{discharge}} - \sum_r g_{i,r,t,\text{charge}} - \sum_\ell K_{i\ell} f_{\ell,t} & \text{KCL : bus} \\ 0 &=\sum_\ell C_{\ell c} x_\ell f_{\ell,t} & \text{KVL : cycles} \\ 0 & \leq g_{i,r,t,\text{discharge}} \leq G_{i,r,\text{discharge}}& \text{discharge limits : storage unit} \\ 0 & \leq g_{i,r,t,\text{charge}} \leq G_{i,r,\text{charge}} & \text{charge limits : storage unit} \\ 0 & \leq e_{i,r,t} \leq E_{i,r} & \text{energy limits : storage unit} \\ e_{i,r,t} &= \eta^0_{i,r,t} e_{i,r,t-1} + \eta^1_{i,r,t}g_{i,r,t,\text{charge}} - \frac{1}{\eta^2_{i,r,t}} g_{i,r,t,\text{discharge}} & \text{consistency : storage unit} \\ e_{i,r,0} & = e_{i,r,|T|-1} & \text{cyclicity : storage unit} \end{align}\]

Decision variables:

\(g_{i,s,t}\) is the generator dispatch at bus \(i\), technology \(s\), time step \(t\),
\(f_{\ell,t}\) is the power flow in line \(\ell\),
\(g_{i,r,t,\text{dis-/charge}}\) denotes the charge and discharge of storage unit \(r\) at bus \(i\) and time step \(t\),
\(e_{i,r,t}\) is the state of charge of storage \(r\) at bus \(i\) and time step \(t\).

Parameters:

\(o_{i,s}\) is the marginal generation cost of technology \(s\) at bus \(i\),
\(x_\ell\) is the reactance of transmission line \(\ell\),
\(K_{i\ell}\) is the incidence matrix,
\(C_{\ell c}\) is the cycle matrix,
\(G_{i,s}\) is the nominal capacity of the generator of technology \(s\) at bus \(i\),
\(F_{\ell}\) is the rating of the transmission line \(\ell\),
\(E_{i,r}\) is the energy capacity of storage \(r\) at bus \(i\),
\(\eta^{0/1/2}_{i,r,t}\) denote the standing (0), charging (1), and discharging (2) efficiencies.

Note

For a full reference to the optimisation problem description, see https://pypsa.readthedocs.io/en/latest/optimal_power_flow.html

Simple electricity market example#

Let’s get acquainted with PyPSA by building simple electricity market models.

We have the following data:

fuel costs in € / MWh\(_{th}\)

fuel_cost = dict(
    coal=8,
    gas=100,
    oil=48,
    hydro=1,
    wind=0,
)

efficiencies of thermal power plants in MWh\(_{el}\) / MWh\(_{th}\)

efficiency = dict(
    coal=0.33,
    gas=0.58,
    oil=0.35,
    hydro=1,
    wind=1,
)

specific emissions in t\(_{CO_2}\) / MWh\(_{th}\)

# t/MWh thermal
emissions = dict(
    coal=0.34,
    gas=0.2,
    oil=0.26,
    hydro=0,
    wind=0,
)

power plant capacities in MW

power_plants = {
    "A": {"coal": 35000, "wind": 3000, "gas": 8000, "oil": 2000},
    "B": {"hydro": 1200},
}

perfectly inelastic electrical load in MW

loads = {
    "A": 42000,
    "B": 650,
}

Building a basic network#

By convention, PyPSA is imported without an alias:

import pypsa

First, we create a new network object which serves as the overall container for all components.

n = pypsa.Network()

The second component we need are buses. Buses are the fundamental nodes of the network, to which all other components like loads, generators and transmission lines attach. They enforce energy conservation for all elements feeding in and out of it (i.e. Kirchhoff’s Current Law).

Components can be added to the network n using the n.add() function. It takes the component name as a first argument, the name of the component as a second argument and possibly further parameters as keyword arguments. Let’s use this function, to add buses for each country to our network:

n.add("Bus", "A", v_nom=380, carrier="AC")
n.add("Bus", "B", v_nom=380, carrier="AC")

Index(['B'], dtype='object')

For each class of components, the data describing the components is stored in a pandas.DataFrame. For example, all static data for buses is stored in n.buses

n.buses

	v_nom	x	y	carrier	v_mag_pu_set	v_mag_pu_min	v_mag_pu_max	control
Bus
A	380.0	0.0	0.0	AC	1.0	0.0	inf	PQ
B	380.0	0.0	0.0	AC	1.0	0.0	inf	PQ

You see there are many more attributes than we specified while adding the buses; many of them are filled with default parameters which were added. You can look up the field description, defaults and status (required input, optional input, output) for buses here https://pypsa.readthedocs.io/en/latest/components.html#bus, and analogous for all other components.

The method n.add() also allows you to add multiple components at once. For instance, multiple carriers for the fuels with information on specific carbon dioxide emissions, a nice name, and colors for plotting. For this, the function takes the component name as the first argument and then a list of component names and then optional arguments for the parameters. Here, scalar values, lists, dictionary or pandas.Series are allowed. The latter two needs keys or indices with the component names.

n.add(
    "Carrier",
    ["coal", "gas", "oil", "hydro", "wind"],
    co2_emissions=emissions,
    nice_name=["Coal", "Gas", "Oil", "Hydro", "Onshore Wind"],
    color=["grey", "indianred", "black", "aquamarine", "dodgerblue"],
)

Index(['coal', 'gas', 'oil', 'hydro', 'wind'], dtype='object')

The n.add() function is very general. It lets you add any component to the network object n. For instance, in the next step we add generators for all the different power plants.

In Region B:

n.add(
    "Generator",
    "B hydro",
    bus="B",
    carrier="hydro",
    p_nom=1200,  # MW
)

Index(['B hydro'], dtype='object')

In Region A (in a loop):

for tech, p_nom in power_plants["A"].items():
    n.add(
        "Generator",
        f"A {tech}",
        bus="A",
        carrier=tech,
        efficiency=efficiency.get(tech, 1),
        p_nom=p_nom,
        marginal_cost=fuel_cost.get(tech, 0) / efficiency.get(tech, 1),
    )

As a result, the n.generators DataFrame looks like this:

n.generators["marginal_cost"]

Generator
B hydro      0.000000
A coal      24.242424
A wind       0.000000
A gas      172.413793
A oil      137.142857
Name: marginal_cost, dtype: float64

Next, we’re going to add the electricity demand.

A positive value for p_set means consumption of power from the bus.

n.add(
    "Load",
    "A electricity demand",
    bus="A",
    p_set=loads["A"],
)

Index(['A electricity demand'], dtype='object')

n.add(
    "Load",
    "B electricity demand",
    bus="B",
    p_set=loads["B"],
)

Index(['B electricity demand'], dtype='object')

n.loads

	bus	p_set	q_set	sign	active
Load
A electricity demand	A	42000.0	0.0	-1.0	True
B electricity demand	B	650.0	0.0	-1.0	True

Finally, we add the connection between Region B and Region A with a 500 MW line:

n.add(
    "Line",
    "A-B",
    bus0="A",
    bus1="B",
    s_nom=500,
    x=1,
    r=1,
)

Index(['A-B'], dtype='object')

n.lines

	bus0	bus1	type	x	r	g	b	s_nom	s_nom_mod	s_nom_extendable	...	v_ang_min	v_ang_max	sub_network	x_pu	r_pu	g_pu	b_pu	x_pu_eff	r_pu_eff	s_nom_opt
Line
A-B	A	B		1.0	1.0	0.0	0.0	500.0	0.0	False	...	-inf	inf		0.0	0.0	0.0	0.0	0.0	0.0	0.0

1 rows × 31 columns

Optimisation#

With all input data transferred into PyPSA’s data structure, we can now build and run the resulting optimisation problem. In PyPSA, building, solving and retrieving results from the optimisation model is contained in a single function call n.optimize(). This function optimizes dispatch and investment decisions for least cost.

The n.optimize() function can take a variety of arguments. The most relevant for the moment is the choice of the solver (e.g. “highs” and “gurobi”). They need to be installed on your computer, to use them here!

n.optimize()

Let’s have a look at the results.

Since the power flow and dispatch are generally time-varying quantities, these are stored in a different location than e.g. n.generators. They are stored in n.generators_t. Thus, to find out the dispatch of the generators, run

n.generators_t.p

Generator	B hydro	A coal	A wind	A gas	A oil
snapshot
now	1150.0	35000.0	3000.0	1500.0	2000.0

or if you prefer it in relation to the generators nominal capacity

n.generators_t.p / n.generators.p_nom

Generator	B hydro	A coal	A wind	A gas	A oil
snapshot
now	0.958333	1.0	1.0	0.1875	1.0

You see that the time index has the value ‘now’. This is the default index when no time series data has been specified and the network only covers a single state (e.g. a particular hour).

Similarly you will find the power flow in transmission lines at

n.lines_t.p0

Line	A-B
snapshot
now	-500.0

n.lines_t.p1

Line	A-B
snapshot
now	500.0

The p0 will tell you the flow from bus0 to bus1. p1 will tell you the flow from bus1 to bus0.

What about the market clearing prices in the electricity market?

n.buses_t.marginal_price

Bus	A	B
snapshot
now	172.413793	-0.0

Modifying networks#

Modifying data of components in an existing PyPSA network can be done by modifying the entries of a pandas.DataFrame. For instance, if we want to reduce the cross-border transmission capacity between Region A and Region B, we’d run:

n.lines.loc["A-B", "s_nom"] = 400

n.lines

	bus0	bus1	type	x	r	g	b	s_nom	s_nom_mod	s_nom_extendable	...	v_ang_max	sub_network	x_pu	r_pu	g_pu	b_pu	x_pu_eff	r_pu_eff	s_nom_opt	v_nom
Line
A-B	A	B		1.0	1.0	0.0	0.0	400.0	0.0	False	...	inf	0	0.000007	0.000007	0.0	0.0	0.000007	0.000007	500.0	380.0

1 rows × 32 columns

n.optimize(solver_name="highs")

You can see that the production of the hydro power plant was reduced and that of the gas power plant increased owing to the reduced transmission capacity.

n.generators_t.p

Generator	B hydro	A coal	A wind	A gas	A oil
snapshot
now	1050.0	35000.0	3000.0	1600.0	2000.0

Exploring the optimization model#

n.model

Linopy LP model
===============

Variables:
----------
 * Generator-p (snapshot, Generator)
 * Line-s (snapshot, Line)

Constraints:
------------
 * Generator-fix-p-lower (snapshot, Generator-fix)
 * Generator-fix-p-upper (snapshot, Generator-fix)
 * Line-fix-s-lower (snapshot, Line-fix)
 * Line-fix-s-upper (snapshot, Line-fix)
 * Bus-nodal_balance (Bus, snapshot)

Status:
-------
ok

n.model.constraints["Generator-fix-p-upper"]

Constraint `Generator-fix-p-upper` [snapshot: 1, Generator-fix: 5]:
-------------------------------------------------------------------
[now, B hydro]: +1 Generator-p[now, B hydro] ≤ 1200.0
[now, A coal]: +1 Generator-p[now, A coal]   ≤ 35000.0
[now, A wind]: +1 Generator-p[now, A wind]   ≤ 3000.0
[now, A gas]: +1 Generator-p[now, A gas]     ≤ 8000.0
[now, A oil]: +1 Generator-p[now, A oil]     ≤ 2000.0

n.model.constraints["Bus-nodal_balance"]

Constraint `Bus-nodal_balance` [Bus: 2, snapshot: 1]:
-----------------------------------------------------
[A, now]: +1 Generator-p[now, A coal] + 1 Generator-p[now, A wind] + 1 Generator-p[now, A gas] + 1 Generator-p[now, A oil] - 1 Line-s[now, A-B] = 42000.0
[B, now]: +1 Generator-p[now, B hydro] + 1 Line-s[now, A-B]                                                                                     = 650.0

n.model.objective

Objective:
----------
LinearExpression: +24.24 Generator-p[now, A coal] + 172.4 Generator-p[now, A gas] + 137.1 Generator-p[now, A oil]
Sense: min
Value: 1398632.6317348

Exercises#

Task 1: Model an outage of the transmission line by setting it’s capacity to zero. How does the model compensate for the lack of transmission?

Show code cell content Hide code cell content

n.lines["s_nom"] = 0.
n.optimize()
n.generators_t.p

WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['A', 'B'], dtype='object', name='Bus')

WARNING:pypsa.consistency:The following sub_networks have carriers which are not defined:
Index(['0'], dtype='object', name='SubNetwork')

WARNING:pypsa.consistency:The following lines have carriers which are not defined:
Index(['A-B'], dtype='object', name='Line')

INFO:linopy.model: Solve problem using Highs solver

INFO:linopy.io: Writing time: 0.02s

INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 6 primals, 14 duals
Objective: 1.47e+06
Solver model: available
Solver message: optimal

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Line-fix-s-lower, Line-fix-s-upper were not assigned to the network.

Running HiGHS 1.9.0 (git hash: fa40bdf): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e+00, 1e+00]
  Cost   [2e+01, 2e+02]
  Bound  [0e+00, 0e+00]
  RHS    [6e+02, 4e+04]
Presolving model
1 rows, 2 cols, 2 nonzeros  0s
0 rows, 0 cols, 0 nonzeros  0s
Presolve : Reductions: rows 0(-14); columns 0(-6); elements 0(-19) - Reduced to empty
Solving the original LP from the solution after postsolve
Model name          : linopy-problem-zdokxnoi
Model status        : Optimal
Objective value     :  1.4675981490e+06
Relative P-D gap    :  0.0000000000e+00
HiGHS run time      :          0.00
Writing the solution to /tmp/linopy-solve-oqnz21ey.sol

Generator	B hydro	A coal	A wind	A gas	A oil
snapshot
now	650.0	35000.0	3000.0	2000.0	2000.0

Task 2: Double the wind capacity. How is the electricity price affected? What generator is price-setting?

Show code cell content Hide code cell content

n.generators.loc["A wind", "p_nom"] *= 2
n.optimize()
n.buses_t.marginal_price

WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['A', 'B'], dtype='object', name='Bus')

WARNING:pypsa.consistency:The following sub_networks have carriers which are not defined:
Index(['0'], dtype='object', name='SubNetwork')

WARNING:pypsa.consistency:The following lines have carriers which are not defined:
Index(['A-B'], dtype='object', name='Line')

INFO:linopy.model: Solve problem using Highs solver

INFO:linopy.io: Writing time: 0.02s

INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 6 primals, 14 duals
Objective: 9.86e+05
Solver model: available
Solver message: optimal

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Line-fix-s-lower, Line-fix-s-upper were not assigned to the network.

Running HiGHS 1.9.0 (git hash: fa40bdf): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e+00, 1e+00]
  Cost   [2e+01, 2e+02]
  Bound  [0e+00, 0e+00]
  RHS    [6e+02, 4e+04]
Presolving model
1 rows, 3 cols, 3 nonzeros  0s
1 rows, 3 cols, 3 nonzeros  0s
Presolve : Reductions: rows 1(-13); columns 3(-3); elements 3(-16)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     0.0000000000e+00 Pr: 1(36000) 0s
          1     9.8562770563e+05 Pr: 0(0) 0s
Solving the original LP from the solution after postsolve
Model name          : linopy-problem-ry37dvk9
Model status        : Optimal
Simplex   iterations: 1
Objective value     :  9.8562770563e+05
Relative P-D gap    :  9.4490299867e-16
HiGHS run time      :          0.00
Writing the solution to /tmp/linopy-solve-zrtamz1_.sol

Bus	A	B
snapshot
now	137.142857	-0.0

Task 3: Region A brings a new 1000 MW nuclear power plant into operation with an estimated marginal electricity generation cost of 20 €/MWh. Will that affect the electricity price?

Show code cell content Hide code cell content

n.add(
    "Generator",
    "A nuclear",
    bus="A",
    carrier="nuclear",
    p_nom=1000,
    marginal_cost=20,
)
n.optimize()
n.buses_t.marginal_price

WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['A', 'B'], dtype='object', name='Bus')

WARNING:pypsa.consistency:The following sub_networks have carriers which are not defined:
Index(['0'], dtype='object', name='SubNetwork')

WARNING:pypsa.consistency:The following lines have carriers which are not defined:
Index(['A-B'], dtype='object', name='Line')

WARNING:pypsa.consistency:The following generators have carriers which are not defined:
Index(['A nuclear'], dtype='object', name='Generator')

INFO:linopy.model: Solve problem using Highs solver

INFO:linopy.io: Writing time: 0.02s

INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 7 primals, 16 duals
Objective: 8.68e+05
Solver model: available
Solver message: optimal

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Line-fix-s-lower, Line-fix-s-upper were not assigned to the network.

Running HiGHS 1.9.0 (git hash: fa40bdf): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e+00, 1e+00]
  Cost   [2e+01, 2e+02]
  Bound  [0e+00, 0e+00]
  RHS    [6e+02, 4e+04]
Presolving model
1 rows, 3 cols, 3 nonzeros  0s
0 rows, 0 cols, 0 nonzeros  0s
Presolve : Reductions: rows 0(-16); columns 0(-7); elements 0(-22) - Reduced to empty
Solving the original LP from the solution after postsolve
Model name          : linopy-problem-q36tzych
Model status        : Optimal
Objective value     :  8.6848484848e+05
Relative P-D gap    :  0.0000000000e+00
HiGHS run time      :          0.00
Writing the solution to /tmp/linopy-solve-p6o7up58.sol

Bus	A	B
snapshot
now	24.242424	-0.0

Economic Dispatch

Contents

Economic Dispatch#

From structured data to optimisation#

Simple electricity market example#

Building a basic network#

Optimisation#

Modifying networks#

Exploring the optimization model#

Exercises#