Replace components with fundamental Links and Stores
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.22.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
INFO:pypsa.linopf:Prepare linear problem
INFO:pypsa.linopf:Total preparation time: 0.3s
INFO:pypsa.linopf:Solve linear problem using Glpk solver
INFO:pypsa.linopf:Optimization successful. Objective value: 1.47e+07
[3]:
('ok', '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.22.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
INFO:pypsa.linopf:Prepare linear problem
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
Cell In[4], line 13
3 su_to_replace = "Storage 0"
5 (
6 bus_name,
7 link_1_name,
(...)
11 extra_functionality,
12 ) = replace_su(network, su_to_replace)
---> 13 network.lopf(
14 network.snapshots, extra_functionality=extra_functionality, formulation="kirchhoff"
15 )
17 assert_almost_equal(network_r.objective, network.objective, decimal=2)
18 assert_array_almost_equal(
19 network_r.storage_units_t.state_of_charge[su_to_replace],
20 network.stores_t.e[store_name],
21 )
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.22.0/lib/python3.11/site-packages/pypsa/components.py:767, in Network.lopf(self, snapshots, pyomo, solver_name, solver_options, solver_logfile, formulation, keep_files, extra_functionality, multi_investment_periods, **kwargs)
765 if pyomo:
766 return network_lopf(self, **args)
--> 767 return network_lopf_lowmem(self, **args)
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.22.0/lib/python3.11/site-packages/pypsa/linopf.py:1460, in network_lopf(n, snapshots, solver_name, solver_logfile, extra_functionality, multi_investment_periods, skip_objective, skip_pre, extra_postprocessing, formulation, keep_references, keep_files, keep_shadowprices, solver_options, warmstart, store_basis, solver_dir)
1457 n.determine_network_topology()
1459 logger.info("Prepare linear problem")
-> 1460 fdp, problem_fn = prepare_lopf(
1461 n, snapshots, keep_files, skip_objective, extra_functionality, solver_dir
1462 )
1463 fds, solution_fn = mkstemp(prefix="pypsa-solve", suffix=".sol", dir=solver_dir)
1465 if warmstart == True:
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.22.0/lib/python3.11/site-packages/pypsa/linopf.py:1129, in prepare_lopf(n, snapshots, keep_files, skip_objective, extra_functionality, solver_dir)
1126 define_objective(n, snapshots)
1128 if extra_functionality is not None:
-> 1129 extra_functionality(n, snapshots)
1131 n.binaries_f.write("end\n")
1133 # explicit closing with file descriptor is necessary for windows machines
Cell In[2], line 160, in replace_su.<locals>.extra_functionality(network, snapshots)
159 def extra_functionality(network, snapshots):
--> 160 model = network.model
161 model.store_fix_1 = Constraint(
162 rule=lambda model: model.store_e_nom[store_name]
163 == model.link_p_nom[link_1_name] * ratio1
164 )
165 model.store_fix_2 = Constraint(
166 rule=lambda model: model.store_e_nom[store_name]
167 == model.link_p_nom[link_2_name] * ratio2
168 )
AttributeError: 'Network' object has no attribute 'model'
[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.22.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
INFO:pypsa.linopf:Prepare linear problem
INFO:pypsa.linopf:Total preparation time: 0.41s
INFO:pypsa.linopf:Solve linear problem using Glpk solver
INFO:pypsa.linopf:Optimization successful. Objective value: 1.47e+07
---------------------------------------------------------------------------
AssertionError Traceback (most recent call last)
Cell In[5], line 8
5 bus_name, link_name, store_name = replace_gen(network, gen_to_replace)
6 network.lopf(network.snapshots)
----> 8 assert_almost_equal(network_r.objective, network.objective, decimal=2)
10 # check dispatch
11 assert_array_almost_equal(
12 -network.links_t.p1[link_name], network_r.generators_t.p[gen_to_replace]
13 )
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.22.0/lib/python3.11/contextlib.py:81, in ContextDecorator.__call__.<locals>.inner(*args, **kwds)
78 @wraps(func)
79 def inner(*args, **kwds):
80 with self._recreate_cm():
---> 81 return func(*args, **kwds)
File ~/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.22.0/lib/python3.11/site-packages/numpy/testing/_private/utils.py:604, in assert_almost_equal(actual, desired, decimal, err_msg, verbose)
602 pass
603 if abs(desired - actual) >= np.float64(1.5 * 10.0**(-decimal)):
--> 604 raise AssertionError(_build_err_msg())
AssertionError:
Arrays are not almost equal to 2 decimals
ACTUAL: 14670510.15
DESIRED: 14670510.2