Add support for units in storage constraints

The implementation of pint units is almost done now.

syfop/node.py

@@ -105,8 +105,7 @@ def __init__(
     def _get_size_commodity(self):
         # TODO refactor this into a property
 
-        # there can be only one output commodity because we don't have ouptut_properties here,
-        # this is tested elsewehere already
+        # there can be only one output commodity because we don't have convertion_factors here
         assert len(set(self.output_commodities)) == 1
         return self.output_commodities[0]
 
@@ -278,6 +277,10 @@ class Storage:
     commodity of its node. Storage for nodes with multiple output commodities are not supported at
     the moment.
 
+    The storage for one node does not support multiple different output commodities at the moment.
+    If you need a storage for different output commodities, create a separate nodes for each
+    commodity and attach a separate storage there.
+
     **Examples:** hydrogen storage, CO2 storage, battery.
 
     Attributes
@@ -294,6 +297,12 @@ class Storage:
         but it is not forbidden in any way. However, such a case will not be optimal if
         ``charging_loss>0``.
 
+    **Note:** The units of the variables ``size``, ``level``, ``charge`` and ``discharge`` are
+    given by the unit of the commodity times hours (independently of the interval size between time
+    stamps). This means for a battery, the variables will be given in `MWh` if the unit for
+    'electricity' is set to `MW`. This means that the values in ``charge`` and ``discharge`` depend
+    on the interval of time stamps.
+
     """
 
     # Note: atm this is not implemented as node class. Probably could be done, but might be more
@@ -305,17 +314,19 @@ def __init__(self, costs, max_charging_speed, storage_loss, charging_loss):
         costs : pint.Quantity
             Storage costs per unit of size, e.g. ``1000 * ureg.EUR/ureg.kWh``.
         max_charging_speed : float
-            Maximum charging speed, i.e. the share of the total size that can be charged per time
-            stamp. For example, if the maximum charging speed is 0.5, two time stamps are needed to
-            charge the storage completely.
+            Maximum charging speed, i.e. the share of the total size that can be charged per hour
+            (indepenent of the length of the interval between time stamps). For example, if the
+            maximum charging speed is 0.5, two hours are needed to charge the storage completely.
+            The same limit is applied for discharging speed.
         storage_loss : float
-            Loss of stored commodity per time stamp as share of the stored commodity. For example,
-            if the storage loss for a battery is 0.01 and the battery is half full, 0.5% of the
-            battery capacity is lost in the next time stamp.
+            Loss of stored commodity per hour (indepenent of the length of the interval between
+            time stamps) as share of the stored commodity. For example, if the storage loss for a
+            battery is 0.01 and the battery is half full, 0.5% of the battery capacity is lost in
+            the next hour.
         charging_loss : float
-            Loss of charged commodity per time stamp as share of the charged commodity. For
-            example, if ``charging_loss`` is 0.01 and there is 100kg of excess hydrogen to be
-            stored in a certain timestamp, only 99kg will end up in the storage.
+            Loss of charged commodity as share of the charged commodity. For example, if
+            ``charging_loss`` is 0.01 and there is 100kg of excess hydrogen to be stored in a
+            certain timestamp, only 99kg will end up in the storage.
 
         """
         self.costs = costs  # per size
@@ -326,6 +337,8 @@ def __init__(self, costs, max_charging_speed, storage_loss, charging_loss):
         # a loss which equals to 1 does not make sense, because everything would be lost
         # if charging_loss == 0. solutions might be indeterministic because charging and
         # discharging might be done in the same time stamp
-        assert 0 <= storage_loss < 1, "storage_loss must be smaller than 1"
-        assert 0 <= charging_loss < 1, "charging_loss must be smaller than 1"
-        assert 0 < max_charging_speed <= 1, "max_charging_speed must not be greater than 1"
+        assert 0 <= storage_loss < 1, "storage_loss must be non-negative and smaller than 1"
+        assert 0 <= charging_loss < 1, "charging_loss must be non-negative and smaller than 1"
+        assert (
+            0 < max_charging_speed <= 1
+        ), "max_charging_speed must be positive and not be greater than 1"

syfop/node_base.py

@@ -1,6 +1,6 @@
 import linopy
 
-from syfop.units import default_units, strip_unit, ureg
+from syfop.units import default_units, interval_length, strip_unit, ureg
 from syfop.util import timeseries_variable
 
 
@@ -135,19 +135,32 @@ def _create_proportion_constraints(self, model, proportions, get_flows):
                 name=f"proportion_{self.name}_{commodity}",
             )
 
+    def _interval_length_h(self):
+        # XXX refactor: pass time_coords as parameter similar to create_variables
+        #
+        # this trusts on the check in Network that all time_coords are identical
+        # also assumes that there is at list one input flow (which is always the case, right?)
+        first_input_flow = list(self.input_flows.values())[0]
+        time_coords = first_input_flow.coords["time"]
+        interval_length_h = interval_length(time_coords).to(ureg.h).magnitude
+        return interval_length_h
+
     def _create_storage_constraints(self, model):
         """This method is not supposed to be called if the node does not have a storage."""
         size = self.storage.size
         level = self.storage.level
         charge = self.storage.charge
         discharge = self.storage.discharge
 
+        max_charging_per_timestamp = (
+            size * self._interval_length_h() * self.storage.max_charging_speed
+        )
         model.add_constraints(
-            charge - size * self.storage.max_charging_speed <= 0,
+            charge - max_charging_per_timestamp <= 0,
             name=f"max_charging_speed_{self.name}",
         )
         model.add_constraints(
-            discharge - size * self.storage.max_charging_speed <= 0,
+            discharge - max_charging_per_timestamp <= 0,
             name=f"max_discharging_speed_{self.name}",
         )
         model.add_constraints(level - size <= 0, name=f"storage_max_level_{self.name}")
@@ -303,7 +316,10 @@ def _create_constraint_inout_flow_balance_commodity(
             rhs = 0
 
         if self.storage is not None:
-            lhs = lhs + self.storage.charge - self.storage.discharge
+            # no need for unit conversion, Example: charge and discharge are MWh and lhs in MW
+            lhs = lhs + 1 / self._interval_length_h() * (
+                self.storage.charge - self.storage.discharge
+            )
 
         model.add_constraints(lhs == rhs, name=f"inout_flow_balance_{self.name}")
 
@@ -333,8 +349,8 @@ def create_constraints(self, model):
 
     def storage_cost_magnitude(self, currency_unit):
         assert hasattr(self, "storage") and self.storage is not None, "node has no storage"
-        storage_unit = default_units[self._get_size_commodity(self)]
-        return self.storage.costs.to(currency_unit / storage_unit).magnitude
+        storage_unit = default_units[self._get_size_commodity()]
+        return self.storage.costs.to(currency_unit / (storage_unit * ureg.h)).magnitude
 
 
 class NodeScalableBase(NodeBase):

tests/test_network.py

@@ -92,7 +92,7 @@ def test_simple_co2_storage(storage_type):
     """
 
     wind_flow = const_time_series(0.5)
-    co2_flow = const_time_series(0.5)
+    co2_flow = const_time_series(0.5) * ureg.t / ureg.h
     co2_storage = None
     electricity_storage = None
     hydrogen_storage = None
@@ -104,7 +104,7 @@ def test_simple_co2_storage(storage_type):
 
     if storage_type == "co2_storage":
         co2_flow = 2 * co2_flow
-        co2_flow[1::2] = 0
+        co2_flow[1::2] = np.array(0.0) * ureg.t / ureg.h
         co2_storage = Storage(
             costs=1000 * ureg.EUR / (ureg.t / ureg.h),  # price not relevant, see comment above
             max_charging_speed=1.0,
@@ -115,7 +115,7 @@ def test_simple_co2_storage(storage_type):
         expected_storage_costs = 1000 * 0.5
     elif storage_type == "electricity_storage":
         electricity_storage = Storage(
-            costs=100,  # price not relevant, see comment above
+            costs=100 * ureg.EUR / ureg.MWh,  # price not relevant, see comment above
             max_charging_speed=1.0,
             storage_loss=0.0,
             charging_loss=0.0,
@@ -126,7 +126,7 @@ def test_simple_co2_storage(storage_type):
         expected_storage_costs = 100 * 3.0 * 0.5
     elif storage_type == "hydrogen_storage":
         hydrogen_storage = Storage(
-            costs=30,  # price not relevant, see comment above
+            costs=30 * ureg.EUR / ureg.t,  # price not relevant, see comment above
             max_charging_speed=1.0,
             storage_loss=0.0,
             charging_loss=0.0,
@@ -145,30 +145,38 @@ def test_simple_co2_storage(storage_type):
     wind = NodeScalableInput(
         name="wind",
         input_profile=wind_flow,
-        costs=1.3,
+        costs=1.3 * ureg.EUR / ureg.MW,
         output_unit="MW",
         storage=electricity_storage,
     )
     hydrogen = Node(
         name="hydrogen",
         inputs=[wind],
         input_commodities="electricity",
-        costs=3,
+        costs=3 * ureg.EUR / (ureg.t / ureg.h),
+        convert_factor=ureg.t / ureg.h / ureg.MW,
         output_unit="t",
         storage=hydrogen_storage,
     )
     co2 = NodeFixInput(
-        name="co2", input_flow=co2_flow, storage=co2_storage, costs=0, output_unit="t"
+        name="co2",
+        input_flow=co2_flow,
+        storage=co2_storage,
+        costs=0,
+        output_unit="t",
     )
 
     methanol_synthesis = Node(
         name="methanol_synthesis",
         inputs=[co2, hydrogen],
         input_commodities=["co2", "hydrogen"],
-        costs=1.2,
+        costs=1.2 * ureg.EUR / (ureg.t / ureg.h),
         output_unit="t",
         size_commodity="methanol",
-        input_proportions={"co2": 0.25, "hydrogen": 0.75},
+        input_proportions={"co2": 1 * ureg.t / ureg.h, "hydrogen": 3 * ureg.t / ureg.h},
+        convert_factors={
+            "methanol": ("hydrogen", 1 / 0.75),
+        },
     )
 
     network = Network([wind, hydrogen, co2, methanol_synthesis])