time-varying efficiencies and standing losses for all components

doc/release_notes.rst

@@ -14,6 +14,12 @@ Upcoming Release
 
 * Fixed interference of io routines with linopy optimisation [`#564 <https://github.com/PyPSA/PyPSA/pull/564>`_, `#567 <https://github.com/PyPSA/PyPSA/pull/567>`_]
 
+* Efficiencies and standing losses of stores, storage units and generators can
+  now be specified as time-varying attributes (``efficiency``,
+  ``efficiency_dispatch``, ``efficiency_store``, ``standing_loss``). For
+  example, this allows specifying temperature-dependent generator efficiencies
+  or evaporation in hydro reservoirs.
+
 * The attributes ``lifetime`` and ``build_year`` are now aggregated with a
   capacity-weighted mean when clustering the network. Previously, these
   attributes had to carry identical values for components that were to be

pypsa/component_attrs/generators.csv

@@ -17,7 +17,7 @@ marginal_cost,static or series,currency/MWh,0,Marginal cost of production of 1 M
 build_year,int,year,0,build year,Input (optional)
 lifetime,float,years,inf,lifetime,Input (optional)
 capital_cost,float,currency/MW,0,Capital cost of extending p_nom by 1 MW.,Input (optional)
-efficiency,float,per unit,1,"Ratio between primary energy and electrical energy, e.g. takes value 0.4 MWh_elec/MWh_thermal for gas. This is required for global constraints on primary energy in OPF.",Input (optional)
+efficiency,static or series,per unit,1,"Ratio between primary energy and electrical energy, e.g. takes value 0.4 MWh_elec/MWh_thermal for gas. This is required for global constraints on primary energy in OPF.",Input (optional)
 committable,boolean,n/a,False,Use unit commitment (only possible if p_nom is not extendable).,Input (optional)
 start_up_cost,float,currency,0,Cost to start up the generator. Only read if committable is True.,Input (optional)
 shut_down_cost,float,currency,0,Cost to shut down the generator. Only read if committable is True.,Input (optional)

pypsa/component_attrs/storage_units.csv

@@ -23,9 +23,9 @@ state_of_charge_set,static or series,MWh,NaN,State of charge set points for snap
 cyclic_state_of_charge,boolean,n/a,False,"Switch: if True, then state_of_charge_initial is ignored and the initial state of charge is set to the final state of charge for the group of snapshots in the OPF (soc[-1] = soc[len(snapshots)-1]).",Input (optional)
 cyclic_state_of_charge_per_period,boolean,n/a,True,"Switch: if True, then the cyclic constraints are applied to each period (first snapshot level if multiindexed) separately.",Input (optional)
 max_hours,float,hours,1,Maximum state of charge capacity in terms of hours at full output capacity p_nom,Input (optional)
-efficiency_store,float,per unit,1,Efficiency of storage on the way into the storage.,Input (optional)
-efficiency_dispatch,float,per unit,1,Efficiency of storage on the way out of the storage.,Input (optional)
-standing_loss,float,per unit,0,Losses per hour to state of charge.,Input (optional)
+efficiency_store,static or series,per unit,1,Efficiency of storage on the way into the storage.,Input (optional)
+efficiency_dispatch,static or series,per unit,1,Efficiency of storage on the way out of the storage.,Input (optional)
+standing_loss,static or series,per unit,0,Losses per hour to state of charge.,Input (optional)
 inflow,static or series,MW,0,"Inflow to the state of charge, e.g. due to river inflow in hydro reservoir.",Input (optional)
 p,series,MW,0,active power at bus (positive if net generation),Output
 p_dispatch,series,MW,0,active power dispatch at bus,Output

pypsa/component_attrs/stores.csv

@@ -18,7 +18,7 @@ q_set,static or series,MVar,0,reactive power set point (for PF),Input (optional)
 sign,float,n/a,1,power sign,Input (optional)
 marginal_cost,static or series,currency/MWh,0,Marginal cost of production of 1 MWh.,Input (optional)
 capital_cost,float,currency/MWh,0,Capital cost of extending e_nom by 1 MWh.,Input (optional)
-standing_loss,float,per unit,0,Losses per hour to energy.,Input (optional)
+standing_loss,static or series,per unit,0,Losses per hour to energy.,Input (optional)
 build_year,int,year,0,build year,Input (optional)
 lifetime,float,years,inf,lifetime,Input (optional)
 p,series,MW,0,active power at bus (positive if net generation),Output

pypsa/linopf.py

@@ -29,8 +29,9 @@
     get_bounds_pu,
     get_extendable_i,
     get_non_extendable_i,
-    nominal_attrs,
 )
+from pypsa.descriptors import get_switchable_as_dense as get_as_dense
+from pypsa.descriptors import nominal_attrs
 from pypsa.linopt import (
     align_with_static_component,
     define_binaries,
@@ -52,7 +53,6 @@
     write_objective,
 )
 from pypsa.pf import _as_snapshots
-from pypsa.pf import get_switchable_as_dense as get_as_dense
 
 agg_group_kwargs = (
     dict(numeric_only=False) if parse(pd.__version__) >= Version("1.3") else {}
@@ -543,9 +543,9 @@ def define_storage_unit_constraints(n, sns):
     # elapsed hours
     eh = expand_series(n.snapshot_weightings.stores[sns], sus_i)
     # efficiencies
-    eff_stand = expand_series(1 - n.df(c).standing_loss, sns).T.pow(eh)
-    eff_dispatch = expand_series(n.df(c).efficiency_dispatch, sns).T
-    eff_store = expand_series(n.df(c).efficiency_store, sns).T
+    eff_stand = (1 - get_as_dense(n, c, "standing_loss", sns)).pow(eh)
+    eff_dispatch = get_as_dense(n, c, "efficiency_dispatch", sns)
+    eff_store = get_as_dense(n, c, "efficiency_store", sns)
 
     soc = get_var(n, c, "state_of_charge")
 
@@ -655,7 +655,7 @@ def define_store_constraints(n, sns):
 
     # elapsed hours
     eh = expand_series(n.snapshot_weightings.stores[sns], stores_i)  # elapsed hours
-    eff_stand = expand_series(1 - n.df(c).standing_loss, sns).T.pow(eh)
+    eff_stand = (1 - get_as_dense(n, c, "standing_loss", sns)).pow(eh)
 
     e = get_var(n, c, "e")
 
@@ -821,11 +821,9 @@ def get_period(n, glc, sns):
         # generators
         gens = n.generators.query("carrier in @emissions.index")
         if not gens.empty:
-            em_pu = gens.carrier.map(emissions) / gens.efficiency
-            em_pu = (
-                weightings["generators"].to_frame("weightings")
-                @ em_pu.to_frame("weightings").T
-            ).loc[period]
+            efficiency = get_as_dense(n, "Generator", "efficiency", inds=gens.index)
+            em_pu = gens.carrier.map(emissions) / efficiency
+            em_pu = em_pu.multiply(weightings.generators, axis=0).loc[period]
             p = get_var(n, "Generator", "p").loc[sns, gens.index].loc[period]
 
             vals = linexpr((em_pu, p), as_pandas=False)

pypsa/opf.py

@@ -39,12 +39,9 @@ class PersistentSolver:
 logger = logging.getLogger(__name__)
 
 
-from pypsa.descriptors import (
-    allocate_series_dataframes,
-    get_switchable_as_dense,
-    get_switchable_as_iter,
-    zsum,
-)
+from pypsa.descriptors import allocate_series_dataframes
+from pypsa.descriptors import get_switchable_as_dense as get_as_dense
+from pypsa.descriptors import get_switchable_as_iter, zsum
 from pypsa.opt import (
     LConstraint,
     LExpression,
@@ -109,8 +106,8 @@ def define_generator_variables_constraints(network, snapshots):
 
     start_i = network.snapshots.get_loc(snapshots[0])
 
-    p_min_pu = get_switchable_as_dense(network, "Generator", "p_min_pu", snapshots)
-    p_max_pu = get_switchable_as_dense(network, "Generator", "p_max_pu", snapshots)
+    p_min_pu = get_as_dense(network, "Generator", "p_min_pu", snapshots)
+    p_max_pu = get_as_dense(network, "Generator", "p_max_pu", snapshots)
 
     ## Define generator dispatch variables ##
 
@@ -689,8 +686,8 @@ def define_storage_variables_constraints(network, snapshots):
 
     ## Define storage dispatch variables ##
 
-    p_max_pu = get_switchable_as_dense(network, "StorageUnit", "p_max_pu", snapshots)
-    p_min_pu = get_switchable_as_dense(network, "StorageUnit", "p_min_pu", snapshots)
+    p_max_pu = get_as_dense(network, "StorageUnit", "p_max_pu", snapshots)
+    p_min_pu = get_as_dense(network, "StorageUnit", "p_min_pu", snapshots)
 
     bounds = {(su, sn): (0, None) for su in ext_sus_i for sn in snapshots}
     bounds.update(
@@ -733,7 +730,7 @@ def su_p_store_bounds(model, su_name, snapshot):
     free_pyomo_initializers(network.model.storage_p_store)
 
     ## Define spillage variables only for hours with inflow>0. ##
-    inflow = get_switchable_as_dense(network, "StorageUnit", "inflow", snapshots)
+    inflow = get_as_dense(network, "StorageUnit", "inflow", snapshots)
     spill_sus_i = sus.index[inflow.max() > 0]  # skip storage units without any inflow
     inflow_gt0_b = inflow > 0
     spill_bounds = {
@@ -834,10 +831,14 @@ def su_p_lower(model, su_name, snapshot):
     # store the combinations with a fixed soc
     fixed_soc = {}
 
-    state_of_charge_set = get_switchable_as_dense(
+    state_of_charge_set = get_as_dense(
         network, "StorageUnit", "state_of_charge_set", snapshots
     )
 
+    standing_loss = get_as_dense(network, "StorageUnit", "standing_loss").T
+    eff_dispatch = get_as_dense(network, "StorageUnit", "efficiency_dispatch").T
+    eff_store = get_as_dense(network, "StorageUnit", "efficiency_store").T
+
     for su in sus.index:
         for i, sn in enumerate(snapshots):
             soc[su, sn] = [[], "==", 0.0]
@@ -847,13 +848,13 @@ def su_p_lower(model, su_name, snapshot):
             if i == 0 and not sus.at[su, "cyclic_state_of_charge"]:
                 previous_state_of_charge = sus.at[su, "state_of_charge_initial"]
                 soc[su, sn][2] -= (
-                    1 - sus.at[su, "standing_loss"]
+                    1 - standing_loss.at[su, sn]
                 ) ** elapsed_hours * previous_state_of_charge
             else:
                 previous_state_of_charge = model.state_of_charge[su, snapshots[i - 1]]
                 soc[su, sn][0].append(
                     (
-                        (1 - sus.at[su, "standing_loss"]) ** elapsed_hours,
+                        (1 - standing_loss.at[su, sn]) ** elapsed_hours,
                         previous_state_of_charge,
                     )
                 )
@@ -873,13 +874,13 @@ def su_p_lower(model, su_name, snapshot):
 
             soc[su, sn][0].append(
                 (
-                    sus.at[su, "efficiency_store"] * elapsed_hours,
+                    eff_store.at[su, sn] * elapsed_hours,
                     model.storage_p_store[su, sn],
                 )
             )
             soc[su, sn][0].append(
                 (
-                    -(1 / sus.at[su, "efficiency_dispatch"]) * elapsed_hours,
+                    -(1 / eff_dispatch.at[su, sn]) * elapsed_hours,
                     model.storage_p_dispatch[su, sn],
                 )
             )
@@ -908,8 +909,8 @@ def define_store_variables_constraints(network, snapshots):
     ext_stores = stores.index[stores.e_nom_extendable]
     fix_stores = stores.index[~stores.e_nom_extendable]
 
-    e_max_pu = get_switchable_as_dense(network, "Store", "e_max_pu", snapshots)
-    e_min_pu = get_switchable_as_dense(network, "Store", "e_min_pu", snapshots)
+    e_max_pu = get_as_dense(network, "Store", "e_max_pu", snapshots)
+    e_min_pu = get_as_dense(network, "Store", "e_min_pu", snapshots)
 
     model = network.model
 
@@ -978,6 +979,8 @@ def store_e_lower(model, store, snapshot):
 
     e = {}
 
+    standing_loss = get_as_dense(network, "Store", "standing_loss").T
+
     for store in stores.index:
         for i, sn in enumerate(snapshots):
             e[store, sn] = LConstraint(sense="==")
@@ -989,13 +992,13 @@ def store_e_lower(model, store, snapshot):
             if i == 0 and not stores.at[store, "e_cyclic"]:
                 previous_e = stores.at[store, "e_initial"]
                 e[store, sn].lhs.constant += (
-                    1 - stores.at[store, "standing_loss"]
+                    1 - standing_loss.at[store, sn]
                 ) ** elapsed_hours * previous_e
             else:
                 previous_e = model.store_e[store, snapshots[i - 1]]
                 e[store, sn].lhs.variables.append(
                     (
-                        (1 - stores.at[store, "standing_loss"]) ** elapsed_hours,
+                        (1 - standing_loss.at[store, sn]) ** elapsed_hours,
                         previous_e,
                     )
                 )
@@ -1053,8 +1056,8 @@ def define_link_flows(network, snapshots):
 
     fixed_links_i = network.links.index[~network.links.p_nom_extendable]
 
-    p_max_pu = get_switchable_as_dense(network, "Link", "p_max_pu", snapshots)
-    p_min_pu = get_switchable_as_dense(network, "Link", "p_min_pu", snapshots)
+    p_max_pu = get_as_dense(network, "Link", "p_max_pu", snapshots)
+    p_min_pu = get_as_dense(network, "Link", "p_min_pu", snapshots)
 
     fixed_lower = p_min_pu.loc[:, fixed_links_i].multiply(
         network.links.loc[fixed_links_i, "p_nom"]
@@ -1429,7 +1432,7 @@ def define_passive_branch_constraints(network, snapshots):
 
     s_max_pu = pd.concat(
         {
-            c: get_switchable_as_dense(network, c, "s_max_pu", snapshots)
+            c: get_as_dense(network, c, "s_max_pu", snapshots)
             for c in network.passive_branch_components
         },
         axis=1,
@@ -1512,9 +1515,9 @@ def define_nodal_balances(network, snapshots):
         (bus, sn): LExpression() for bus in network.buses.index for sn in snapshots
     }
 
-    efficiency = get_switchable_as_dense(network, "Link", "efficiency", snapshots)
+    efficiency = get_as_dense(network, "Link", "efficiency", snapshots)
 
-    filter = (get_switchable_as_dense(network, "Link", "p_min_pu", snapshots) < 0) & (
+    filter = (get_as_dense(network, "Link", "p_min_pu", snapshots) < 0) & (
         efficiency < 1
     )
     links = filter[filter].dropna(how="all", axis=1)
@@ -1543,9 +1546,7 @@ def define_nodal_balances(network, snapshots):
         for col in network.links.columns
         if col[:3] == "bus" and col not in ["bus0", "bus1"]
     ]:
-        efficiency = get_switchable_as_dense(
-            network, "Link", "efficiency{}".format(i), snapshots
-        )
+        efficiency = get_as_dense(network, "Link", "efficiency{}".format(i), snapshots)
         for cb in network.links.index[network.links["bus{}".format(i)] != ""]:
             bus = network.links.at[cb, "bus{}".format(i)]
             for sn in snapshots:
@@ -1561,7 +1562,7 @@ def define_nodal_balances(network, snapshots):
                 (sign, network.model.generator_p[gen, sn])
             )
 
-    load_p_set = get_switchable_as_dense(network, "Load", "p_set", snapshots)
+    load_p_set = get_as_dense(network, "Load", "p_set", snapshots)
     for load in network.loads.index:
         bus = network.loads.at[load, "bus"]
         sign = network.loads.at[load, "sign"]
@@ -1659,11 +1660,12 @@ def define_global_constraints(network, snapshots):
                     continue
                 # for generators, use the prime mover carrier
                 gens = network.generators.index[network.generators.carrier == carrier]
+                efficiency = get_as_dense(network, "Generator", "efficiency").T
                 c.lhs.variables.extend(
                     [
                         (
                             attribute
-                            * (1 / network.generators.at[gen, "efficiency"])
+                            * (1 / efficiency.at[gen, sn])
                             * network.snapshot_weightings.generators[sn],
                             network.model.generator_p[gen, sn],
                         )
@@ -1927,7 +1929,7 @@ def get_shadows(constraint, multiind=True):
         set_from_series(network.stores_t.e, get_values(model.store_e))
 
     if len(network.loads):
-        load_p_set = get_switchable_as_dense(network, "Load", "p_set", snapshots)
+        load_p_set = get_as_dense(network, "Load", "p_set", snapshots)
         network.loads_t["p"].loc[snapshots] = load_p_set.loc[snapshots]
 
     if len(network.buses):
@@ -1965,7 +1967,7 @@ def get_shadows(constraint, multiind=True):
     if len(network.links):
         set_from_series(network.links_t.p0, get_values(model.link_p))
 
-        efficiency = get_switchable_as_dense(network, "Link", "efficiency", snapshots)
+        efficiency = get_as_dense(network, "Link", "efficiency", snapshots)
 
         network.links_t.p1.loc[snapshots] = (
             -network.links_t.p0.loc[snapshots] * efficiency.loc[snapshots]
@@ -1991,7 +1993,7 @@ def get_shadows(constraint, multiind=True):
             for col in network.links.columns
             if col[:3] == "bus" and col not in ["bus0", "bus1"]
         ]:
-            efficiency = get_switchable_as_dense(
+            efficiency = get_as_dense(
                 network, "Link", "efficiency{}".format(i), snapshots
             )
             p_name = "p{}".format(i)

pypsa/optimization/constraints.py

@@ -652,9 +652,9 @@ def define_storage_unit_constraints(n, sns):
     # elapsed hours
     eh = expand_series(n.snapshot_weightings.stores[sns], assets.index)
     # efficiencies
-    eff_stand = expand_series(1 - assets.standing_loss, sns).T.pow(eh)
-    eff_dispatch = expand_series(assets.efficiency_dispatch, sns).T
-    eff_store = expand_series(assets.efficiency_store, sns).T
+    eff_stand = (1 - get_as_dense(n, c, "standing_loss", sns)).pow(eh)
+    eff_dispatch = get_as_dense(n, c, "efficiency_dispatch", sns)
+    eff_store = get_as_dense(n, c, "efficiency_store", sns)
 
     soc = m[f"{c}-state_of_charge"]
 
@@ -737,7 +737,7 @@ def define_store_constraints(n, sns):
     # elapsed hours
     eh = expand_series(n.snapshot_weightings.stores[sns], assets.index)
     # efficiencies
-    eff_stand = expand_series(1 - assets.standing_loss, sns).T.pow(eh)
+    eff_stand = (1 - get_as_dense(n, c, "standing_loss", sns)).pow(eh)
 
     e = m[c + "-e"]
     p = m[c + "-p"]

pypsa/optimization/global_constraints.py

@@ -12,6 +12,7 @@
 from numpy import isnan, nan
 from xarray import DataArray
 
+from pypsa.descriptors import get_switchable_as_dense as get_as_dense
 from pypsa.descriptors import nominal_attrs
 
 logger = logging.getLogger(__name__)
@@ -268,9 +269,9 @@ def define_primary_energy_limit(n, sns):
         # generators
         gens = n.generators.query("carrier in @emissions.index")
         if not gens.empty:
-            w = weightings["generators"].to_frame("weight")
-            em_pu = (gens.carrier.map(emissions) / gens.efficiency).to_frame("weight")
-            em_pu = w @ em_pu.T
+            efficiency = get_as_dense(n, "Generator", "efficiency", inds=gens.index)
+            em_pu = gens.carrier.map(emissions) / efficiency
+            em_pu = em_pu.multiply(weightings.generators, axis=0)
             p = m["Generator-p"].loc[snapshots, gens.index]
             expr = (p * em_pu).sum()
             lhs.append(expr)