Note
You can download this example as a Jupyter notebook or start it in interactive mode.
Replace components with fundamental Links and Stores
This notebook demonstrates how generators and storage units can be replaced by more fundamental components, and how their parameters map to each other.
[1]:
import pypsa, os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from numpy.testing import assert_almost_equal, assert_array_almost_equal
from pyomo.environ import Constraint
We define two functions we use in the following.
[2]:
def replace_gen(network, gen_to_replace):
"""Replace the generator gen_to_replace with a bus for the energy
carrier, a link for the conversion from the energy carrier to electricity
and a store to keep track of the depletion of the energy carrier and its
CO2 emissions."""
gen = network.generators.loc[gen_to_replace]
bus_name = "{} {}".format(gen["bus"], gen["carrier"])
link_name = "{} converter {} to AC".format(gen_to_replace, gen["carrier"])
store_name = "{} store {}".format(gen_to_replace, gen["carrier"])
network.add("Bus", bus_name, carrier=gen["carrier"])
network.add(
"Link",
link_name,
bus0=bus_name,
bus1=gen["bus"],
capital_cost=gen["capital_cost"] * gen["efficiency"],
p_nom=gen["p_nom"] / gen["efficiency"],
p_nom_extendable=gen["p_nom_extendable"],
p_nom_max=gen["p_nom_max"] / gen["efficiency"],
p_nom_min=gen["p_nom_min"] / gen["efficiency"],
p_max_pu=network.generators_t.p_max_pu.loc[:, gen_to_replace]
if gen_to_replace in network.generators_t.p_max_pu.columns
else gen["p_max_pu"],
p_min_pu=network.generators_t.p_min_pu.loc[:, gen_to_replace]
if gen_to_replace in network.generators_t.p_min_pu.columns
else gen["p_min_pu"],
marginal_cost=gen["marginal_cost"] * gen["efficiency"],
efficiency=gen["efficiency"],
)
network.add(
"Store",
store_name,
bus=bus_name,
e_nom_min=-float("inf"),
e_nom_max=0,
e_nom_extendable=True,
e_min_pu=1.0,
e_max_pu=0.0,
)
network.remove("Generator", gen_to_replace)
return bus_name, link_name, store_name
def replace_su(network, su_to_replace):
"""Replace the storage unit su_to_replace with a bus for the energy
carrier, two links for the conversion of the energy carrier to and from electricity,
a store to keep track of the depletion of the energy carrier and its
CO2 emissions, and a variable generator for the storage inflow.
Because the energy size and power size are linked in the storage unit by the max_hours,
extra functionality must be added to the LOPF to implement this constraint."""
su = network.storage_units.loc[su_to_replace]
bus_name = "{} {}".format(su["bus"], su["carrier"])
link_1_name = "{} converter {} to AC".format(su_to_replace, su["carrier"])
link_2_name = "{} converter AC to {}".format(su_to_replace, su["carrier"])
store_name = "{} store {}".format(su_to_replace, su["carrier"])
gen_name = "{} inflow".format(su_to_replace)
network.add("Bus", bus_name, carrier=su["carrier"])
# dispatch link
network.add(
"Link",
link_1_name,
bus0=bus_name,
bus1=su["bus"],
capital_cost=su["capital_cost"] * su["efficiency_dispatch"],
p_nom=su["p_nom"] / su["efficiency_dispatch"],
p_nom_extendable=su["p_nom_extendable"],
p_nom_max=su["p_nom_max"] / su["efficiency_dispatch"],
p_nom_min=su["p_nom_min"] / su["efficiency_dispatch"],
p_max_pu=su["p_max_pu"],
marginal_cost=su["marginal_cost"] * su["efficiency_dispatch"],
efficiency=su["efficiency_dispatch"],
)
# store link
network.add(
"Link",
link_2_name,
bus1=bus_name,
bus0=su["bus"],
p_nom=su["p_nom"],
p_nom_extendable=su["p_nom_extendable"],
p_nom_max=su["p_nom_max"],
p_nom_min=su["p_nom_min"],
p_max_pu=-su["p_min_pu"],
efficiency=su["efficiency_store"],
)
if (
su_to_replace in network.storage_units_t.state_of_charge_set.columns
and (
~pd.isnull(network.storage_units_t.state_of_charge_set[su_to_replace])
).any()
):
e_max_pu = pd.Series(data=1.0, index=network.snapshots)
e_min_pu = pd.Series(data=0.0, index=network.snapshots)
non_null = ~pd.isnull(
network.storage_units_t.state_of_charge_set[su_to_replace]
)
e_max_pu[non_null] = network.storage_units_t.state_of_charge_set[su_to_replace][
non_null
]
e_min_pu[non_null] = network.storage_units_t.state_of_charge_set[su_to_replace][
non_null
]
else:
e_max_pu = 1.0
e_min_pu = 0.0
network.add(
"Store",
store_name,
bus=bus_name,
e_nom=su["p_nom"] * su["max_hours"],
e_nom_min=su["p_nom_min"] / su["efficiency_dispatch"] * su["max_hours"],
e_nom_max=su["p_nom_max"] / su["efficiency_dispatch"] * su["max_hours"],
e_nom_extendable=su["p_nom_extendable"],
e_max_pu=e_max_pu,
e_min_pu=e_min_pu,
standing_loss=su["standing_loss"],
e_cyclic=su["cyclic_state_of_charge"],
e_initial=su["state_of_charge_initial"],
)
network.add("Carrier", "rain", co2_emissions=0.0)
# inflow from a variable generator, which can be curtailed (i.e. spilled)
inflow_max = network.storage_units_t.inflow[su_to_replace].max()
if inflow_max == 0.0:
inflow_pu = 0.0
else:
inflow_pu = network.storage_units_t.inflow[su_to_replace] / inflow_max
network.add(
"Generator",
gen_name,
bus=bus_name,
carrier="rain",
p_nom=inflow_max,
p_max_pu=inflow_pu,
)
if su["p_nom_extendable"]:
ratio2 = su["max_hours"]
ratio1 = ratio2 * su["efficiency_dispatch"]
def extra_functionality(network, snapshots):
model = network.model
model.store_fix_1 = Constraint(
rule=lambda model: model.store_e_nom[store_name]
== model.link_p_nom[link_1_name] * ratio1
)
model.store_fix_2 = Constraint(
rule=lambda model: model.store_e_nom[store_name]
== model.link_p_nom[link_2_name] * ratio2
)
else:
extra_functionality = None
network.remove("StorageUnit", su_to_replace)
return bus_name, link_1_name, link_2_name, store_name, gen_name, extra_functionality
Now, take an example from the git repo which has already been solved
[3]:
network_r = pypsa.examples.storage_hvdc(from_master=True)
network_r.lopf()
WARNING:pypsa.io:Importing network from PyPSA version v0.17.1 while current version is v0.21.0. Read the release notes at https://pypsa.readthedocs.io/en/latest/release_notes.html to prepare your network for import.
INFO:pypsa.io:Imported network storage-hvdc.nc has buses, carriers, generators, global_constraints, lines, links, loads, storage_units
WARNING:pypsa.components:Solving optimisation problem with pyomo.In PyPSA version 0.21 the default will change to ``n.lopf(pyomo=False)``.Explicitly set ``n.lopf(pyomo=True)`` to retain current behaviour.
INFO:pypsa.opf:Performed preliminary steps
INFO:pypsa.opf:Building pyomo model using `kirchhoff` formulation
INFO:pypsa.opf:Solving model using glpk
INFO:pypsa.opf:Optimization successful
# ==========================================================
# = Solver Results =
# ==========================================================
# ----------------------------------------------------------
# Problem Information
# ----------------------------------------------------------
Problem:
- Name: unknown
Lower bound: 14670508.8253592
Upper bound: 14670508.8253592
Number of objectives: 1
Number of constraints: 796
Number of variables: 490
Number of nonzeros: 2049
Sense: minimize
# ----------------------------------------------------------
# Solver Information
# ----------------------------------------------------------
Solver:
- Status: ok
Termination condition: optimal
Statistics:
Branch and bound:
Number of bounded subproblems: 0
Number of created subproblems: 0
Error rc: 0
Time: 0.03682851791381836
# ----------------------------------------------------------
# Solution Information
# ----------------------------------------------------------
Solution:
- number of solutions: 0
number of solutions displayed: 0
[3]:
(<SolverStatus.ok: 'ok'>, <TerminationCondition.optimal: 'optimal'>)
[4]:
network = pypsa.examples.storage_hvdc(from_master=True)
su_to_replace = "Storage 0"
(
bus_name,
link_1_name,
link_2_name,
store_name,
gen_name,
extra_functionality,
) = replace_su(network, su_to_replace)
network.lopf(
network.snapshots, extra_functionality=extra_functionality, formulation="kirchhoff"
)
assert_almost_equal(network_r.objective, network.objective, decimal=2)
assert_array_almost_equal(
network_r.storage_units_t.state_of_charge[su_to_replace],
network.stores_t.e[store_name],
)
assert_array_almost_equal(
network_r.storage_units_t.p[su_to_replace],
-network.links_t.p1[link_1_name] - network.links_t.p0[link_2_name],
)
# check optimised size
assert_array_almost_equal(
network_r.storage_units.at[su_to_replace, "p_nom_opt"],
network.links.at[link_2_name, "p_nom_opt"],
)
assert_array_almost_equal(
network_r.storage_units.at[su_to_replace, "p_nom_opt"],
network.links.at[link_1_name, "p_nom_opt"]
* network_r.storage_units.at[su_to_replace, "efficiency_dispatch"],
)
WARNING:pypsa.io:Importing network from PyPSA version v0.17.1 while current version is v0.21.0. Read the release notes at https://pypsa.readthedocs.io/en/latest/release_notes.html to prepare your network for import.
INFO:pypsa.io:Imported network storage-hvdc.nc has buses, carriers, generators, global_constraints, lines, links, loads, storage_units
WARNING:pypsa.components:Solving optimisation problem with pyomo.In PyPSA version 0.21 the default will change to ``n.lopf(pyomo=False)``.Explicitly set ``n.lopf(pyomo=True)`` to retain current behaviour.
INFO:pypsa.opf:Performed preliminary steps
INFO:pypsa.opf:Building pyomo model using `kirchhoff` formulation
INFO:pypsa.opf:Solving model using glpk
INFO:pypsa.opf:Optimization successful
# ==========================================================
# = Solver Results =
# ==========================================================
# ----------------------------------------------------------
# Problem Information
# ----------------------------------------------------------
Problem:
- Name: unknown
Lower bound: 14670508.8253592
Upper bound: 14670508.8253592
Number of objectives: 1
Number of constraints: 846
Number of variables: 504
Number of nonzeros: 2113
Sense: minimize
# ----------------------------------------------------------
# Solver Information
# ----------------------------------------------------------
Solver:
- Status: ok
Termination condition: optimal
Statistics:
Branch and bound:
Number of bounded subproblems: 0
Number of created subproblems: 0
Error rc: 0
Time: 0.029987812042236328
# ----------------------------------------------------------
# Solution Information
# ----------------------------------------------------------
Solution:
- number of solutions: 0
number of solutions displayed: 0
[5]:
network = pypsa.examples.storage_hvdc(from_master=True)
gen_to_replace = "Gas 0"
bus_name, link_name, store_name = replace_gen(network, gen_to_replace)
network.lopf(network.snapshots)
assert_almost_equal(network_r.objective, network.objective, decimal=2)
# check dispatch
assert_array_almost_equal(
-network.links_t.p1[link_name], network_r.generators_t.p[gen_to_replace]
)
# check optimised size
assert_array_almost_equal(
network_r.generators.at[gen_to_replace, "p_nom_opt"],
network.links.at[link_name, "p_nom_opt"]
* network.links.at[link_name, "efficiency"],
)
WARNING:pypsa.io:Importing network from PyPSA version v0.17.1 while current version is v0.21.0. Read the release notes at https://pypsa.readthedocs.io/en/latest/release_notes.html to prepare your network for import.
INFO:pypsa.io:Imported network storage-hvdc.nc has buses, carriers, generators, global_constraints, lines, links, loads, storage_units
WARNING:pypsa.components:Solving optimisation problem with pyomo.In PyPSA version 0.21 the default will change to ``n.lopf(pyomo=False)``.Explicitly set ``n.lopf(pyomo=True)`` to retain current behaviour.
INFO:pypsa.opf:Performed preliminary steps
INFO:pypsa.opf:Building pyomo model using `kirchhoff` formulation
INFO:pypsa.opf:Solving model using glpk
INFO:pypsa.opf:Optimization successful
# ==========================================================
# = Solver Results =
# ==========================================================
# ----------------------------------------------------------
# Problem Information
# ----------------------------------------------------------
Problem:
- Name: unknown
Lower bound: 14670508.8253592
Upper bound: 14670508.8253592
Number of objectives: 1
Number of constraints: 868
Number of variables: 515
Number of nonzeros: 2157
Sense: minimize
# ----------------------------------------------------------
# Solver Information
# ----------------------------------------------------------
Solver:
- Status: ok
Termination condition: optimal
Statistics:
Branch and bound:
Number of bounded subproblems: 0
Number of created subproblems: 0
Error rc: 0
Time: 0.03361344337463379
# ----------------------------------------------------------
# Solution Information
# ----------------------------------------------------------
Solution:
- number of solutions: 0
number of solutions displayed: 0