StorageUnit as Link and Store Components¶
This example demonstrates how the StorageUnit component can be replaced by the more fundamental Store and Link components, and how their parameters map to each other.
In [2]:
Copied!
import pandas as pd
from numpy.testing import assert_almost_equal, assert_array_almost_equal
import pypsa
import pandas as pd
from numpy.testing import assert_almost_equal, assert_array_almost_equal
import pypsa
We define two functions we use in the following.
In [3]:
Copied!
def replace_storage_unit(n: pypsa.Network, name: str) -> tuple:
"""
Replace the storage unit with `name` 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 optimisation to implement this constraint.
"""
su = n.storage_units.loc[name]
bus_name = f"{su['bus']} {su['carrier']}"
link_1_name = f"{name} converter {su['carrier']} to AC"
link_2_name = f"{name} converter AC to {su['carrier']}"
store_name = f"{name} store {su['carrier']}"
gen_name = f"{name} inflow"
n.add("Bus", bus_name, carrier=su["carrier"])
# dispatch link
n.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
n.add(
"Link",
link_2_name,
bus0=su["bus"],
bus1=bus_name,
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 (
name in n.storage_units_t.state_of_charge_set.columns
and (~pd.isnull(n.storage_units_t.state_of_charge_set[name])).any()
):
e_max_pu = pd.Series(data=1, index=n.snapshots)
e_min_pu = pd.Series(data=0, index=n.snapshots)
non_null = ~pd.isnull(n.storage_units_t.state_of_charge_set[name])
e_max_pu[non_null] = n.storage_units_t.state_of_charge_set[name][non_null]
e_min_pu[non_null] = n.storage_units_t.state_of_charge_set[name][non_null]
else:
e_max_pu = 1
e_min_pu = 0
n.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"],
)
n.add("Carrier", "rain", co2_emissions=0)
# inflow from a variable generator, which can be curtailed (i.e. spilled)
inflow_max = n.storage_units_t.inflow[name].max()
if inflow_max == 0:
inflow_pu = 0
else:
inflow_pu = n.storage_units_t.inflow[name] / inflow_max
n.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(n: pypsa.Network, sns: pd.Index) -> None:
m = n.model
lhs = (
m["Store-e_nom"].at[store_name]
- m["Link-p_nom"].at[link_1_name] * ratio1
)
m.add_constraints(lhs == 0, name="store_fix_1")
lhs = (
m["Store-e_nom"].at[store_name]
- m["Link-p_nom"].at[link_2_name] * ratio2
)
m.add_constraints(lhs == 0, name="store_fix_2")
else:
extra_functionality = None
n.remove("StorageUnit", name)
return bus_name, link_1_name, link_2_name, store_name, gen_name, extra_functionality
def replace_storage_unit(n: pypsa.Network, name: str) -> tuple:
"""
Replace the storage unit with `name` 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 optimisation to implement this constraint.
"""
su = n.storage_units.loc[name]
bus_name = f"{su['bus']} {su['carrier']}"
link_1_name = f"{name} converter {su['carrier']} to AC"
link_2_name = f"{name} converter AC to {su['carrier']}"
store_name = f"{name} store {su['carrier']}"
gen_name = f"{name} inflow"
n.add("Bus", bus_name, carrier=su["carrier"])
# dispatch link
n.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
n.add(
"Link",
link_2_name,
bus0=su["bus"],
bus1=bus_name,
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 (
name in n.storage_units_t.state_of_charge_set.columns
and (~pd.isnull(n.storage_units_t.state_of_charge_set[name])).any()
):
e_max_pu = pd.Series(data=1, index=n.snapshots)
e_min_pu = pd.Series(data=0, index=n.snapshots)
non_null = ~pd.isnull(n.storage_units_t.state_of_charge_set[name])
e_max_pu[non_null] = n.storage_units_t.state_of_charge_set[name][non_null]
e_min_pu[non_null] = n.storage_units_t.state_of_charge_set[name][non_null]
else:
e_max_pu = 1
e_min_pu = 0
n.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"],
)
n.add("Carrier", "rain", co2_emissions=0)
# inflow from a variable generator, which can be curtailed (i.e. spilled)
inflow_max = n.storage_units_t.inflow[name].max()
if inflow_max == 0:
inflow_pu = 0
else:
inflow_pu = n.storage_units_t.inflow[name] / inflow_max
n.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(n: pypsa.Network, sns: pd.Index) -> None:
m = n.model
lhs = (
m["Store-e_nom"].at[store_name]
- m["Link-p_nom"].at[link_1_name] * ratio1
)
m.add_constraints(lhs == 0, name="store_fix_1")
lhs = (
m["Store-e_nom"].at[store_name]
- m["Link-p_nom"].at[link_2_name] * ratio2
)
m.add_constraints(lhs == 0, name="store_fix_2")
else:
extra_functionality = None
n.remove("StorageUnit", name)
return bus_name, link_1_name, link_2_name, store_name, gen_name, extra_functionality
Now, take an existing example which has already been solved as a reference case:
In [4]:
Copied!
n_r = pypsa.examples.storage_hvdc()
n_r.optimize()
n_r = pypsa.examples.storage_hvdc()
n_r.optimize()
/home/runner/work/PyPSA/PyPSA/pypsa/network/io.py:2082: FutureWarning: pandas infers the `str` dtype for string data since its version 3.0. PyPSA still converts it back to numpy object dtype on import, but will keep it from PyPSA 2.0 on. Set `pypsa.options.api.legacy_string_dtype` explicitly to suppress this warning. new_static = _coerce_string_dtypes(new_static)
INFO:pypsa.network.io:Imported network 'Storage-HVDC' has buses, carriers, generators, global_constraints, lines, links, loads, storage_units
/tmp/ipykernel_7768/2597710260.py:2: 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_r.optimize() WARNING:pypsa.consistency:The following buses have carriers which are not defined. Run n.sanitize() to add them. Components with undefined carriers: Index(['0', '1', '2', '3', '4', '5'], 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(['TL 0', 'TL 1'], dtype='object', name='name')
WARNING:pypsa.consistency:The following lines have carriers which are not defined. Run n.sanitize() to add them. Components with undefined carriers: Index(['0', '1', '2', '3', '4', '5'], dtype='object', name='name')
WARNING:pypsa.consistency:The following lines have zero r, which could break the linear load flow: Index(['0', '1', '2', '3', '4', '5'], dtype='object', name='name')
INFO:linopy.model: Solve problem using Gurobi solver
INFO:linopy.model:Solver options: - log_to_console: False
INFO:linopy.io: Writing time: 0.13s
Set parameter WLSAccessID
Set parameter WLSSecret
Set parameter LicenseID to value 2537914
Academic license 2537914 - for non-commercial use only - registered to l.___@tu-berlin.de
Read LP format model from file /tmp/linopy-problem-g8l4jpcq.lp
Reading time = 0.00 seconds
obj: 1104 rows, 490 columns, 2359 nonzeros
Set parameter LogToConsole to value 0
INFO:linopy.constants: Optimization successful: Status: ok Termination condition: optimal Solution: 550 primals, 1174 duals Objective: 1.47e+07 Solver: gurobi Runtime: 0.01s Dual bound: 1.47e+07 Solver model: available Solver message: 2
INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Generator-ext-p-lower, Generator-ext-p-upper, Line-fix-s-lower, Line-fix-s-upper, Line-ext-s-lower, Line-ext-s-upper, Link-fix-p-lower, Link-fix-p-upper, Link-ext-p-lower, Link-ext-p-upper, StorageUnit-fix-p_dispatch-lower, StorageUnit-fix-p_dispatch-upper, StorageUnit-ext-p_dispatch-lower, StorageUnit-ext-p_dispatch-upper, StorageUnit-fix-p_store-lower, StorageUnit-fix-p_store-upper, StorageUnit-ext-p_store-lower, StorageUnit-ext-p_store-upper, StorageUnit-fix-state_of_charge-lower, StorageUnit-fix-state_of_charge-upper, StorageUnit-ext-state_of_charge-lower, StorageUnit-ext-state_of_charge-upper, StorageUnit-state_of_charge_set, Kirchhoff-Voltage-Law, StorageUnit-energy_balance were not assigned to the network.
Out[4]:
('ok', 'optimal')
We load it again, but now we replace the StorageUnit components with Store and Link components.
In [5]:
Copied!
n = pypsa.examples.storage_hvdc()
name = "Storage 0"
(
bus_name,
link_1_name,
link_2_name,
store_name,
gen_name,
extra_functionality,
) = replace_storage_unit(n, name)
n.optimize(extra_functionality=extra_functionality)
n = pypsa.examples.storage_hvdc()
name = "Storage 0"
(
bus_name,
link_1_name,
link_2_name,
store_name,
gen_name,
extra_functionality,
) = replace_storage_unit(n, name)
n.optimize(extra_functionality=extra_functionality)
INFO:pypsa.network.io:Imported network 'Storage-HVDC' has buses, carriers, generators, global_constraints, lines, links, loads, storage_units
/tmp/ipykernel_7768/1974559407.py:13: 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(extra_functionality=extra_functionality) WARNING:pypsa.consistency:The following buses have carriers which are not defined. Run n.sanitize() to add them. Components with undefined carriers: Index(['0', '1', '2', '3', '4', '5'], 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(['TL 0', 'TL 1', 'Storage 0 converter AC to battery'], dtype='object', name='name')
WARNING:pypsa.consistency:The following lines have carriers which are not defined. Run n.sanitize() to add them. Components with undefined carriers: Index(['0', '1', '2', '3', '4', '5'], dtype='object', name='name')
WARNING:pypsa.consistency:The following lines have zero r, which could break the linear load flow: Index(['0', '1', '2', '3', '4', '5'], dtype='object', name='name')
INFO:linopy.model: Solve problem using Gurobi solver
INFO:linopy.model:Solver options: - log_to_console: False
INFO:linopy.io: Writing time: 0.14s
Set parameter WLSAccessID
Set parameter WLSSecret
Set parameter LicenseID to value 2537914
Academic license 2537914 - for non-commercial use only - registered to l.___@tu-berlin.de
Read LP format model from file /tmp/linopy-problem-7uspiorp.lp
Reading time = 0.00 seconds
obj: 1144 rows, 504 columns, 2413 nonzeros
Set parameter LogToConsole to value 0
INFO:linopy.constants: Optimization successful: Status: ok Termination condition: optimal Solution: 504 primals, 1202 duals Objective: 1.47e+07 Solver: gurobi Runtime: 0.01s Dual bound: 1.47e+07 Solver model: available Solver message: 2
INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Generator-ext-p-lower, Generator-ext-p-upper, Line-fix-s-lower, Line-fix-s-upper, Line-ext-s-lower, Line-ext-s-upper, Link-fix-p-lower, Link-fix-p-upper, Link-ext-p-lower, Link-ext-p-upper, Store-ext-e-lower, Store-ext-e-upper, StorageUnit-fix-p_dispatch-lower, StorageUnit-fix-p_dispatch-upper, StorageUnit-ext-p_dispatch-lower, StorageUnit-ext-p_dispatch-upper, StorageUnit-fix-p_store-lower, StorageUnit-fix-p_store-upper, StorageUnit-ext-p_store-lower, StorageUnit-ext-p_store-upper, StorageUnit-fix-state_of_charge-lower, StorageUnit-fix-state_of_charge-upper, StorageUnit-ext-state_of_charge-lower, StorageUnit-ext-state_of_charge-upper, StorageUnit-state_of_charge_set, Kirchhoff-Voltage-Law, StorageUnit-energy_balance, Store-energy_balance, store_fix_1, store_fix_2 were not assigned to the network.
Out[5]:
('ok', 'optimal')
Let's ensure that the conversion was accurate by comparing the results of the original and modified networks.
In [6]:
Copied!
assert_almost_equal(n_r.objective, n.objective, decimal=2)
assert_array_almost_equal(
n_r.storage_units_t.state_of_charge[name],
n.stores_t.e[store_name],
)
assert_array_almost_equal(
n_r.storage_units_t.p[name],
-n.links_t.p1[link_1_name] - n.links_t.p0[link_2_name],
)
assert_array_almost_equal(
n_r.storage_units.at[name, "p_nom_opt"],
n.links.at[link_2_name, "p_nom_opt"],
)
assert_array_almost_equal(
n_r.storage_units.at[name, "p_nom_opt"],
n.links.at[link_1_name, "p_nom_opt"]
* n_r.storage_units.at[name, "efficiency_dispatch"],
)
assert_almost_equal(n_r.objective, n.objective, decimal=2)
assert_array_almost_equal(
n_r.storage_units_t.state_of_charge[name],
n.stores_t.e[store_name],
)
assert_array_almost_equal(
n_r.storage_units_t.p[name],
-n.links_t.p1[link_1_name] - n.links_t.p0[link_2_name],
)
assert_array_almost_equal(
n_r.storage_units.at[name, "p_nom_opt"],
n.links.at[link_2_name, "p_nom_opt"],
)
assert_array_almost_equal(
n_r.storage_units.at[name, "p_nom_opt"],
n.links.at[link_1_name, "p_nom_opt"]
* n_r.storage_units.at[name, "efficiency_dispatch"],
)