Electricity Markets¶

This example gradually builds up more and more complicated energy-only electricity markets in PyPSA, starting from a single bidding zone, going up to multiple bidding zones connected with transmission (NTCs) along with variable renewables and storage.

Preliminaries¶

Here libraries are imported and data is defined.

In [2]:

import numpy as np

import pypsa

import numpy as np import pypsa

In [3]:

# marginal costs in EUR/MWh
marginal_costs = {"Wind": 0, "Hydro": 0, "Coal": 30, "Gas": 60, "Oil": 80}

# power plant capacities (nominal powers in MW) in each country (not necessarily realistic)
power_plant_p_nom = {
    "South Africa": {"Coal": 35000, "Wind": 3000, "Gas": 8000, "Oil": 2000},
    "Mozambique": {
        "Hydro": 1200,
    },
    "Eswatini": {
        "Hydro": 600,
    },
}

# transmission capacities in MW (not necessarily realistic)
transmission = {
    "South Africa": {"Mozambique": 500, "Eswatini": 250},
    "Mozambique": {"Eswatini": 100},
}

# country electrical loads in MW (not necessarily realistic)
loads = {"South Africa": 42000, "Mozambique": 650, "Eswatini": 250}

# marginal costs in EUR/MWh marginal_costs = {"Wind": 0, "Hydro": 0, "Coal": 30, "Gas": 60, "Oil": 80} # power plant capacities (nominal powers in MW) in each country (not necessarily realistic) power_plant_p_nom = { "South Africa": {"Coal": 35000, "Wind": 3000, "Gas": 8000, "Oil": 2000}, "Mozambique": { "Hydro": 1200, }, "Eswatini": { "Hydro": 600, }, } # transmission capacities in MW (not necessarily realistic) transmission = { "South Africa": {"Mozambique": 500, "Eswatini": 250}, "Mozambique": {"Eswatini": 100}, } # country electrical loads in MW (not necessarily realistic) loads = {"South Africa": 42000, "Mozambique": 650, "Eswatini": 250}

Single bidding zone with fixed load, one period¶

In this example we consider a single market bidding zone, South Africa.

The inelastic load has essentially infinite marginal utility (or higher than the marginal cost of any generator).

In [4]:

country = "South Africa"

n = pypsa.Network()

n.add("Bus", country)

for tech in power_plant_p_nom[country]:
    n.add(
        "Generator",
        f"{country} {tech}",
        bus=country,
        p_nom=power_plant_p_nom[country][tech],
        marginal_cost=marginal_costs[tech],
    )


n.add("Load", f"{country} load", bus=country, p_set=loads[country]);

country = "South Africa" n = pypsa.Network() n.add("Bus", country) for tech in power_plant_p_nom[country]: n.add( "Generator", f"{country} {tech}", bus=country, p_nom=power_plant_p_nom[country][tech], marginal_cost=marginal_costs[tech], ) n.add("Load", f"{country} load", bus=country, p_set=loads[country]);

Run optimisation to determine market dispatch:

In [5]:

n.optimize()

n.optimize()

/tmp/ipykernel_8167/1261279110.py:1: FutureWarning: The default value of `include_objective_constant` will change from True to False in version 2.0. Set `include_objective_constant` explicitly to suppress this warning. Using False improves LP numerical conditioning by not including the objective constant as a variable.
  n.optimize()
WARNING:pypsa.consistency:The following buses have carriers which are not defined. Run n.sanitize() to add them. Components with undefined carriers:
Index(['South Africa'], dtype='object', name='name')

INFO:linopy.model: Solve problem using Highs solver

INFO:linopy.model:Solver options:
 - log_to_console: False

INFO:linopy.io: Writing time: 0.02s

INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 4 primals, 9 duals
Objective: 1.29e+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 were not assigned to the network.

Out[5]:

('ok', 'optimal')

Print the load active power (P) consumption:

In [6]:

n.loads_t.p

n.loads_t.p

Out[6]:

name	South Africa load
snapshot
now	42000.0

Print the generator active power (P) dispatch:

In [7]:

n.generators_t.p

n.generators_t.p

Out[7]:

name	South Africa Coal	South Africa Wind	South Africa Gas	South Africa Oil
snapshot
now	35000.0	3000.0	4000.0	-0.0

Print the clearing price, which corresponds to gas:

In [8]:

n.buses_t.marginal_price

n.buses_t.marginal_price

Out[8]:

name	South Africa
snapshot
now	60.0

Two bidding zones connected by transmission, one period¶

In this example we have bidirectional lossless transmission capacity between two bidding zones. The power transfer is treated as controllable (like an A/NTC (Available/Net Transfer Capacity) or HVDC line). Note that in the physical grid, power would flow passively according to the network impedances.

In [9]:

n = pypsa.Network()

countries = ["Mozambique", "South Africa"]

for country in countries:
    n.add("Bus", country)

    for tech in power_plant_p_nom[country]:
        n.add(
            "Generator",
            f"{country} {tech}",
            bus=country,
            p_nom=power_plant_p_nom[country][tech],
            marginal_cost=marginal_costs[tech],
        )

    n.add("Load", f"{country} load", bus=country, p_set=loads[country])

    # add transmission as controllable Link
    if country not in transmission:
        continue

    for other_country in countries:
        if other_country not in transmission[country]:
            continue

        # NB: Link is by default unidirectional, so have to set p_min_pu = -1
        # to allow bidirectional (i.e. also negative) flow
        n.add(
            "Link",
            f"{country} - {other_country} link",
            bus0=country,
            bus1=other_country,
            p_nom=transmission[country][other_country],
            p_min_pu=-1,
        )

n = pypsa.Network() countries = ["Mozambique", "South Africa"] for country in countries: n.add("Bus", country) for tech in power_plant_p_nom[country]: n.add( "Generator", f"{country} {tech}", bus=country, p_nom=power_plant_p_nom[country][tech], marginal_cost=marginal_costs[tech], ) n.add("Load", f"{country} load", bus=country, p_set=loads[country]) # add transmission as controllable Link if country not in transmission: continue for other_country in countries: if other_country not in transmission[country]: continue # NB: Link is by default unidirectional, so have to set p_min_pu = -1 # to allow bidirectional (i.e. also negative) flow n.add( "Link", f"{country} - {other_country} link", bus0=country, bus1=other_country, p_nom=transmission[country][other_country], p_min_pu=-1, )

In [10]:

n.optimize()

n.optimize()

/tmp/ipykernel_8167/1261279110.py:1: FutureWarning: The default value of `include_objective_constant` will change from True to False in version 2.0. Set `include_objective_constant` explicitly to suppress this warning. Using False improves LP numerical conditioning by not including the objective constant as a variable.
  n.optimize()
WARNING:pypsa.consistency:The following buses have carriers which are not defined. Run n.sanitize() to add them. Components with undefined carriers:
Index(['Mozambique', 'South Africa'], dtype='object', name='name')

WARNING:pypsa.consistency:The following links have carriers which are not defined. Run n.sanitize() to add them. Components with undefined carriers:
Index(['South Africa - Mozambique link'], dtype='object', name='name')

INFO:linopy.model: Solve problem using Highs solver

INFO:linopy.model:Solver options:
 - log_to_console: False

INFO:linopy.io: Writing time: 0.02s

INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 6 primals, 14 duals
Objective: 1.26e+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, Link-fix-p-lower, Link-fix-p-upper were not assigned to the network.

Out[10]:

('ok', 'optimal')

In [11]:

n.loads_t.p

n.loads_t.p

Out[11]:

name	Mozambique load	South Africa load
snapshot
now	650.0	42000.0

In [12]:

n.generators_t.p

n.generators_t.p

Out[12]:

name	Mozambique Hydro	South Africa Coal	South Africa Wind	South Africa Gas	South Africa Oil
snapshot
now	1150.0	35000.0	3000.0	3500.0	-0.0

The clearing price corresponds to hydro in Mozambique and gas in South Africa.

In [13]:

n.buses_t.marginal_price

n.buses_t.marginal_price

Out[13]:

name	Mozambique	South Africa
snapshot
now	-0.0	60.0

In [14]:

n.links_t.p0

n.links_t.p0

Out[14]:

name	South Africa - Mozambique link
snapshot
now	-500.0

The shadow prices of the links measure the inframarginal rent of the link, i.e. the difference between the marginal price of the two bidding zones.

In [15]:

n.links_t.mu_lower

n.links_t.mu_lower

Out[15]:

name
snapshot
now

Three bidding zones connected by transmission, one period¶

In this example we have bidirectional lossless transmission capacity between three bidding zones. The power transfer is treated as controllable (like an A/NTC (Available/Net Transfer Capacity) or HVDC line). Note that in the physical grid, power would flow passively according to the network impedances.

In [16]:

n = pypsa.Network()

countries = ["Eswatini", "Mozambique", "South Africa"]

for country in countries:
    n.add("Bus", country)

    for tech in power_plant_p_nom[country]:
        n.add(
            "Generator",
            f"{country} {tech}",
            bus=country,
            p_nom=power_plant_p_nom[country][tech],
            marginal_cost=marginal_costs[tech],
        )

    n.add("Load", f"{country} load", bus=country, p_set=loads[country])

    if country not in transmission:
        continue

    for other_country in countries:
        if other_country not in transmission[country]:
            continue

        n.add(
            "Link",
            f"{country} - {other_country} link",
            bus0=country,
            bus1=other_country,
            p_nom=transmission[country][other_country],
            p_min_pu=-1,
        )

n = pypsa.Network() countries = ["Eswatini", "Mozambique", "South Africa"] for country in countries: n.add("Bus", country) for tech in power_plant_p_nom[country]: n.add( "Generator", f"{country} {tech}", bus=country, p_nom=power_plant_p_nom[country][tech], marginal_cost=marginal_costs[tech], ) n.add("Load", f"{country} load", bus=country, p_set=loads[country]) if country not in transmission: continue for other_country in countries: if other_country not in transmission[country]: continue n.add( "Link", f"{country} - {other_country} link", bus0=country, bus1=other_country, p_nom=transmission[country][other_country], p_min_pu=-1, )

In [17]:

n.optimize()

n.optimize()

/tmp/ipykernel_8167/1261279110.py:1: FutureWarning: The default value of `include_objective_constant` will change from True to False in version 2.0. Set `include_objective_constant` explicitly to suppress this warning. Using False improves LP numerical conditioning by not including the objective constant as a variable.
  n.optimize()
WARNING:pypsa.consistency:The following buses have carriers which are not defined. Run n.sanitize() to add them. Components with undefined carriers:
Index(['Eswatini', 'Mozambique', 'South Africa'], dtype='object', name='name')

WARNING:pypsa.consistency:The following links have carriers which are not defined. Run n.sanitize() to add them. Components with undefined carriers:
Index(['Mozambique - Eswatini link', 'South Africa - Eswatini link',
       'South Africa - Mozambique link'],
      dtype='object', name='name')

INFO:linopy.model: Solve problem using Highs solver

INFO:linopy.model:Solver options:
 - log_to_console: False

INFO:linopy.io: Writing time: 0.02s

INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 9 primals, 21 duals
Objective: 1.24e+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, Link-fix-p-lower, Link-fix-p-upper were not assigned to the network.

Out[17]:

('ok', 'optimal')

In [18]:

n.loads_t.p

n.loads_t.p

Out[18]:

name	Eswatini load	Mozambique load	South Africa load
snapshot
now	250.0	650.0	42000.0

In [19]:

n.generators_t.p

n.generators_t.p

Out[19]:

name	Eswatini Hydro	Mozambique Hydro	South Africa Coal	South Africa Wind	South Africa Gas	South Africa Oil
snapshot
now	600.0	1050.0	35000.0	3000.0	3250.0	-0.0

The clearing prices correspond to hydro in Eswatini and Mozambique, and gas in South Africa.

In [20]:

n.buses_t.marginal_price

n.buses_t.marginal_price

Out[20]:

name	Eswatini	Mozambique	South Africa
snapshot
now	-0.0	-0.0	60.0

In [21]:

n.links_t.p0

n.links_t.p0

Out[21]:

name	Mozambique - Eswatini link	South Africa - Eswatini link	South Africa - Mozambique link
snapshot
now	-100.0	-250.0	-500.0

In [22]:

n.links_t.mu_lower

n.links_t.mu_lower

Out[22]:

name
snapshot
now

Single bidding zone with price-sensitive industrial load, one period¶

In this example we consider a single market bidding zone, South Africa.

Now there is a large industrial load with a marginal utility which is low enough to interact with the generation marginal cost. See also the demand elasticity example.

In [23]:

country = "South Africa"

n = pypsa.Network()

n.add("Bus", country)

for tech in power_plant_p_nom[country]:
    n.add(
        "Generator",
        f"{country} {tech}",
        bus=country,
        p_nom=power_plant_p_nom[country][tech],
        marginal_cost=marginal_costs[tech],
    )

# standard high marginal utility consumers
n.add("Load", f"{country} load", bus=country, p_set=loads[country])

# add an industrial load as a negative-dispatch generator with marginal utility of 70 EUR/MWh for 8000 MW
n.add(
    "Generator",
    f"{country} industrial load",
    bus=country,
    p_max_pu=0,
    p_min_pu=-1,
    p_nom=8000,
    marginal_cost=70,
);

country = "South Africa" n = pypsa.Network() n.add("Bus", country) for tech in power_plant_p_nom[country]: n.add( "Generator", f"{country} {tech}", bus=country, p_nom=power_plant_p_nom[country][tech], marginal_cost=marginal_costs[tech], ) # standard high marginal utility consumers n.add("Load", f"{country} load", bus=country, p_set=loads[country]) # add an industrial load as a negative-dispatch generator with marginal utility of 70 EUR/MWh for 8000 MW n.add( "Generator", f"{country} industrial load", bus=country, p_max_pu=0, p_min_pu=-1, p_nom=8000, marginal_cost=70, );

In [24]:

n.optimize()

n.optimize()

/tmp/ipykernel_8167/1261279110.py:1: FutureWarning: The default value of `include_objective_constant` will change from True to False in version 2.0. Set `include_objective_constant` explicitly to suppress this warning. Using False improves LP numerical conditioning by not including the objective constant as a variable.
  n.optimize()
WARNING:pypsa.consistency:The following buses have carriers which are not defined. Run n.sanitize() to add them. Components with undefined carriers:
Index(['South Africa'], dtype='object', name='name')

INFO:linopy.model: Solve problem using Highs solver

INFO:linopy.model:Solver options:
 - log_to_console: False

INFO:linopy.io: Writing time: 0.01s

INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 5 primals, 11 duals
Objective: 1.25e+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 were not assigned to the network.

Out[24]:

('ok', 'optimal')

In [25]:

n.loads_t.p

n.loads_t.p

Out[25]:

name	South Africa load
snapshot
now	42000.0

A look at the generator dispatch shows that only half of industrial load is served, because this maxes out gas; oil is too expensive with a marginal cost of 80 EUR/MWh compared to the industrial load marginal utility of 70 EUR/MWh.

In [26]:

n.generators_t.p

n.generators_t.p

Out[26]:

name	South Africa Coal	South Africa Wind	South Africa Gas	South Africa Oil	South Africa industrial load
snapshot
now	35000.0	3000.0	8000.0	-0.0	-4000.0

In [27]:

n.buses_t.marginal_price

n.buses_t.marginal_price

Out[27]:

name	South Africa
snapshot
now	70.0

Single bidding zone with fixed load, several periods¶

In this example we consider a single market bidding zone, South Africa.

We consider multiple time periods (labelled [0,1,2,3]) to represent variable wind generation.

In [28]:

country = "South Africa"

n = pypsa.Network()

n.set_snapshots(range(4))

n.add("Bus", country)

# availability (p_max_pu) is variable for wind
for tech in power_plant_p_nom[country]:
    n.add(
        "Generator",
        f"{country} {tech}",
        bus=country,
        p_nom=power_plant_p_nom[country][tech],
        marginal_cost=marginal_costs[tech],
        p_max_pu=([0.3, 0.6, 0.4, 0.5] if tech == "Wind" else 1),
    )

# load which varies over the snapshots
n.add(
    "Load",
    f"{country} load",
    bus=country,
    p_set=loads[country] + np.array([0, 1000, 3000, 4000]),
);

country = "South Africa" n = pypsa.Network() n.set_snapshots(range(4)) n.add("Bus", country) # availability (p_max_pu) is variable for wind for tech in power_plant_p_nom[country]: n.add( "Generator", f"{country} {tech}", bus=country, p_nom=power_plant_p_nom[country][tech], marginal_cost=marginal_costs[tech], p_max_pu=([0.3, 0.6, 0.4, 0.5] if tech == "Wind" else 1), ) # load which varies over the snapshots n.add( "Load", f"{country} load", bus=country, p_set=loads[country] + np.array([0, 1000, 3000, 4000]), );

In [29]:

n.optimize()

n.optimize()

/tmp/ipykernel_8167/1261279110.py:1: FutureWarning: The default value of `include_objective_constant` will change from True to False in version 2.0. Set `include_objective_constant` explicitly to suppress this warning. Using False improves LP numerical conditioning by not including the objective constant as a variable.
  n.optimize()
WARNING:pypsa.consistency:The following buses have carriers which are not defined. Run n.sanitize() to add them. Components with undefined carriers:
Index(['South Africa'], dtype='object', name='name')

INFO:linopy.model: Solve problem using Highs solver

INFO:linopy.model:Solver options:
 - log_to_console: False

INFO:linopy.io: Writing time: 0.01s

INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 16 primals, 36 duals
Objective: 6.08e+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 were not assigned to the network.

Out[29]:

('ok', 'optimal')

In [30]:

n.loads_t.p

n.loads_t.p

Out[30]:

name	South Africa load
snapshot
0	42000.0
1	43000.0
2	45000.0
3	46000.0

In [31]:

n.generators_t.p

n.generators_t.p

Out[31]:

name	South Africa Coal	South Africa Wind	South Africa Gas	South Africa Oil
snapshot
0	35000.0	900.0	6100.0	-0.0
1	35000.0	1800.0	6200.0	-0.0
2	35000.0	1200.0	8000.0	800.0
3	35000.0	1500.0	8000.0	1500.0

In [32]:

n.buses_t.marginal_price

n.buses_t.marginal_price

Out[32]:

name	South Africa
snapshot
0	60.0
1	60.0
2	80.0
3	80.0

Single bidding zone with fixed load and storage, several periods¶

In this example we consider a single market bidding zone, South Africa.

We consider multiple time periods (labelled [0,1,2,3]) to represent variable wind generation. Storage is allowed to do price arbitrage to reduce oil consumption.

In [33]:

country = "South Africa"

n = pypsa.Network()

# snapshots labelled by [0,1,2,3]
n.set_snapshots(range(4))

n.add("Bus", country)

# p_max_pu is variable for wind
for tech in power_plant_p_nom[country]:
    n.add(
        "Generator",
        f"{country} {tech}",
        bus=country,
        p_nom=power_plant_p_nom[country][tech],
        marginal_cost=marginal_costs[tech],
        p_max_pu=([0.3, 0.6, 0.4, 0.5] if tech == "Wind" else 1),
    )

# load which varies over the snapshots
n.add(
    "Load",
    f"{country} load",
    bus=country,
    p_set=loads[country] + np.array([0, 1000, 3000, 4000]),
)

# storage unit to do price arbitrage
n.add(
    "StorageUnit",
    f"{country} pumped hydro",
    bus=country,
    p_nom=1000,
    max_hours=6,  # energy storage in terms of hours at full power
)

country = "South Africa" n = pypsa.Network() # snapshots labelled by [0,1,2,3] n.set_snapshots(range(4)) n.add("Bus", country) # p_max_pu is variable for wind for tech in power_plant_p_nom[country]: n.add( "Generator", f"{country} {tech}", bus=country, p_nom=power_plant_p_nom[country][tech], marginal_cost=marginal_costs[tech], p_max_pu=([0.3, 0.6, 0.4, 0.5] if tech == "Wind" else 1), ) # load which varies over the snapshots n.add( "Load", f"{country} load", bus=country, p_set=loads[country] + np.array([0, 1000, 3000, 4000]), ) # storage unit to do price arbitrage n.add( "StorageUnit", f"{country} pumped hydro", bus=country, p_nom=1000, max_hours=6, # energy storage in terms of hours at full power )

In [34]:

n.optimize()

n.optimize()

/tmp/ipykernel_8167/1261279110.py:1: FutureWarning: The default value of `include_objective_constant` will change from True to False in version 2.0. Set `include_objective_constant` explicitly to suppress this warning. Using False improves LP numerical conditioning by not including the objective constant as a variable.
  n.optimize()
WARNING:pypsa.consistency:The following buses have carriers which are not defined. Run n.sanitize() to add them. Components with undefined carriers:
Index(['South Africa'], dtype='object', name='name')

INFO:linopy.model: Solve problem using Highs solver

INFO:linopy.model:Solver options:
 - log_to_console: False

INFO:linopy.io: Writing time: 0.04s

INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 28 primals, 64 duals
Objective: 6.05e+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, StorageUnit-fix-p_dispatch-lower, StorageUnit-fix-p_dispatch-upper, StorageUnit-fix-p_store-lower, StorageUnit-fix-p_store-upper, StorageUnit-fix-state_of_charge-lower, StorageUnit-fix-state_of_charge-upper, StorageUnit-energy_balance were not assigned to the network.

Out[34]:

('ok', 'optimal')

In [35]:

n.loads_t.p

n.loads_t.p

Out[35]:

name	South Africa load
snapshot
0	42000.0
1	43000.0
2	45000.0
3	46000.0

In [36]:

n.generators_t.p

n.generators_t.p

Out[36]:

name	South Africa Coal	South Africa Wind	South Africa Gas	South Africa Oil
snapshot
0	35000.0	900.0	6900.0	-0.0
1	35000.0	1800.0	7200.0	-0.0
2	35000.0	1200.0	8000.0	-0.0
3	35000.0	1500.0	8000.0	500.0

In [37]:

n.storage_units_t.p

n.storage_units_t.p

Out[37]:

name	South Africa pumped hydro
snapshot
0	-800.0
1	-1000.0
2	800.0
3	1000.0

In [38]:

n.storage_units_t.state_of_charge

n.storage_units_t.state_of_charge

Out[38]:

name	South Africa pumped hydro
snapshot
0	800.0
1	1800.0
2	1000.0
3	-0.0

In [39]:

n.buses_t.marginal_price

n.buses_t.marginal_price

Out[39]:

name	South Africa
snapshot
0	60.0
1	60.0
2	60.0
3	80.0