Repository: oogeso/oogeso

Issue:
Update CCS model
Improve CCS model using information from CleanOffHub project

Discussion:

Patch:
diff --git a/docs/changelog.md b/docs/changelog.md
index 1b85da4..29a8570 100644
--- a/docs/changelog.md
+++ b/docs/changelog.md
@@ -4,6 +4,10 @@
 The format is based on [Keep a Changelog](http://keepachangelog.com/)
 and this project adheres to [Semantic Versioning](http://semver.org/).
 
+## v1.5.0
+2024-12-05
+- Add steam cycle and modified CCS (allowing improved modelling of Gas turbine / steam-cycle /CCS setup)
+
 ## v1.4.0
 2024-09-30
 - add gas linepack storage device type (simplified model)
diff --git a/docs/input_data_devices.md b/docs/input_data_devices.md
index 5267ea8..3c3670e 100644
--- a/docs/input_data_devices.md
+++ b/docs/input_data_devices.md
@@ -86,6 +86,7 @@ startupCost | float | cost (NOK) for each start-up
 startupDelay    | float | delay (min) from start-up activation to powr output
 hydrogen_blend_max | float | maximum fraction of hydrogen in fuel blend (Sm3 hydrogen vs Sm3 total)
 hydrogen_blend_min | float | minimum fraction of hydrogen in fuel blend (Sm3 hydrogen vs Sm3 total)
+exhaust_gas_recirculation   | float | Exhaust gas recirculation (EGR) rate (default=0)
 
 ### ```source_el```
 
@@ -143,15 +144,19 @@ f_inj           | float | Ratio of gas injection rate (Sm3/s) to oil production
 injectionpressure   | float | Gas injection pressure (MPa)
 separatorpressure   | float | Pressure (from well) into separator (MPa)
 
+### ```steam_cycle```
+parameter | type | description
+----------|------|------------
+exhaust_gas_recirculation | float | Exhaust gas recirculation (EGR) rate, typically in the range 0-0.6
+alpha                     | float | Factorspecifying how much steam is extracted for CCS (e.g. 0.238)
+
 ### ```carbon_capture```
 parameter | type | description
 ----------|------|------------
-carbon_capture_rate | float | Carbon capture rate (CCR), typically 0.9
-exhaust_gas_recirculation   | float | Exhaust gas recirculatio rate (EGR), typically in the range 0-0.6
-compressor_pressure_in      | float | Inlet pressure for compression of captured CO2 (MPa)
-compressor_pressure_out     | float | Outlet pressure for compression of captured CO2 (MPa)
-compressor_eta              | float | Efficiency of CO2 compresstion (0-1)
-compressor_temp_in          | float | Inlet gas temperature (K)
+carbon_capture_rate                 | float | Carbon capture rate (CCR), typically 0.9
+capture_el_demand_MJ_per_kgCO2      | float | Electricity demand per CO2 captured (MJ/kg)
+capture_heat_demand_MJ_per_kgCO2    | float | Heat demand per CO2 captured (MJ/kg)
+compressor_el_demand_MJ_per_kgCO2   | float | CO2 compressor demand per CO2 captured (MJ/kg)
 
 ### ```storage_hydrogen```
 parameter | type | description
@@ -170,4 +175,15 @@ parameter | type | description
 E_init   | float | Storage level initially (Sm3)
 volume_m3 | float | Volume of pipeline
 
-Maximum storage can be specified by setting pressure limits on associated terminals (done via edge parameters)
\ No newline at end of file
+Maximum storage can be specified by setting pressure limits on associated terminals (done via edge parameters)
+
+### ```water_injection```
+parameter | type | description
+----------|------|------------
+flow_avg | float | Baseline injection rate (Sm3/s)
+E_max    | float | Max positive deviation from baseline, positive value (Sm3)
+E_min    | float | Max negative deviation from baselinem, negative value (Sm3)
+E_init   | float | Baseline value (typically 0) (Sm3)
+E_cost   | float | Cost for deviation from baseline
+target_profile | string | (optional) name of profile used for desired injection rate
+
diff --git a/pyproject.toml b/pyproject.toml
index a1f3189..a232150 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,6 +1,6 @@
 [tool.poetry]
 name = "oogeso"
-version = "1.4.0"
+version = "1.5.0"
 description = "Offshore Oil and Gas Field Energy System Operational Optimisation (OOGESO)"
 authors = ["Harald Svendsen <harald.svendsen@sintef.no>"]
 license = "MIT License (http://opensource.org/licenses/MIT)"
diff --git a/src/oogeso/core/devices/__init__.py b/src/oogeso/core/devices/__init__.py
index bdb2942..fe60151 100644
--- a/src/oogeso/core/devices/__init__.py
+++ b/src/oogeso/core/devices/__init__.py
@@ -3,11 +3,11 @@ from .compressor import CompressorEl, CompressorGas
 from .electrolyser import Electrolyser
 from .fuelcell import FuelCell
 from .gasheater import GasHeater
-from .gasturbine import GasTurbine
+from .gasturbine import GasTurbine, SteamCycle
 from .heatpump import HeatPump
 from .pump import PumpOil, PumpWater, PumpWellStream
 from .separator import Separator, Separator2
 from .sink import PowerSink, SinkCarbon, SinkEl, SinkGas, SinkHeat, SinkOil, SinkWater
 from .source import PowerSource, SourceEl, SourceGas, SourceHydrogen, SourceOil, SourceWater
-from .storage import StorageEl, StorageGasLinepack, StorageHydrogen
+from .storage import StorageEl, StorageGasLinepack, StorageHydrogen, WaterInjection
 from .well import WellGasLift, WellProduction
diff --git a/src/oogeso/core/devices/carbon_capture.py b/src/oogeso/core/devices/carbon_capture.py
index 0415c35..7576bf6 100644
--- a/src/oogeso/core/devices/carbon_capture.py
+++ b/src/oogeso/core/devices/carbon_capture.py
@@ -27,79 +27,53 @@ class CarbonCapture(Device):
         self.carrier_data = carrier_data_dict
 
     @classmethod
-    def co2frac(cls, egr):
-        """The %CO2 is a function of the GT load and of the EGR rate selected.
-        In a first approximation we can make it only a function of the EGR with
-        a piecewise linear relationship:"""
-        # TODO: compare this with gas_co2_content - two parameters for the same thing?
-        # 3.3% if no recirculation.
-        # compare with default co2_content = 2.34 # kg_co2/Sm3_gas (SSB 2016 report)
-        if egr < 0.3:
-            y = (4.96 * egr + 3.3358) / 100
-        else:
-            y = (12.88 * egr + 0.7874) / 100
-        return y
-
-    @classmethod
-    def specific_heat_demand(cls, frac_co2):
-        y = 0
-        if frac_co2 < 0.0489:
-            y = -3.57 * frac_co2 + 4.1551
-        else:
-            y = -2.11 * frac_co2 + 4.0837
+    def required_heat(cls, flow_in_co2, ccr, specific_demand_MW_per_kgCO2):
+        # f_co2 = cls.co2frac(egr)
+        # y = flow_in_co2 * (1 - egr) * ccr * cls.specific_heat_demand(f_co2)
+        y = flow_in_co2 * ccr * specific_demand_MW_per_kgCO2
         return y
 
     @classmethod
-    def specific_electricity_demand(cls, frac_co2):
-        y = 0
-        if frac_co2 < 0.0489:
-            y = 0.16 * frac_co2 + 0.0232
-        else:
-            y = 0.009 * frac_co2 + 0.0305
+    def required_electricity(self, flow_in_co2, ccr, specific_demand_MW_per_kgCO2):
+        # f_co2 = cls.co2frac(egr)
+        # y = flow_in_co2 * (1 - egr) * ccr * cls.specific_electricity_demand(f_co2)
+        y = flow_in_co2 * ccr * specific_demand_MW_per_kgCO2
         return y
 
     @classmethod
-    def required_heat(cls, flow_in_co2, egr, ccr):
-        f_co2 = cls.co2frac(egr)
-        y = flow_in_co2 * (1 - egr) * ccr * cls.specific_heat_demand(f_co2)
-        return y
-
-    @classmethod
-    def required_electricity(cls, flow_in_co2, egr, ccr):
-        f_co2 = cls.co2frac(egr)
-        y = flow_in_co2 * (1 - egr) * ccr * cls.specific_electricity_demand(f_co2)
-        return y
-
-    def required_electricity_compressor(self, flow_in_co2):
+    def required_electricity_compressor(cls, flow_in_co2, energy_demand_MJ_per_kg):
         Q_kg_per_s = flow_in_co2  # kg/s
-        energy_demand_MJ_per_kg = self.dev_data.compressor_energy_demand  # MJ/kg
         P_MW = Q_kg_per_s * energy_demand_MJ_per_kg
         return P_MW
 
     def _rules_carbon_capture(self, model: pyo.Model, t: int, i: int) -> Union[pyo.Expression, pyo.Constraint.Skip]:
         dev = self.id
         ccr = self.dev_data.carbon_capture_rate  # 0-1
-        egr = self.dev_data.exhaust_gas_recirculation  # 0-0.6
+        ccs_el_demand = self.dev_data.capture_el_demand_MJ_per_kgCO2
+        ccs_heat_demand = self.dev_data.capture_heat_demand_MJ_per_kgCO2
+        compressor_el_demand = self.dev_data.compressor_el_demand_MJ_per_kgCO2
 
         co2_flow_in = model.varDeviceFlow[dev, "carbon", "in", t]  # kg/s
 
         if i == 1:
             # heat consumption (heat in) is a linear function of carbon flow
             lhs = model.varDeviceFlow[dev, "heat", "in", t]
-            rhs = self.required_heat(co2_flow_in, egr=egr, ccr=ccr)
+            rhs = self.required_heat(co2_flow_in, ccr=ccr, specific_demand_MW_per_kgCO2=ccs_heat_demand)
             return lhs == rhs
         elif i == 2:
             # electricity consumption for carbon capture process and co2 compression
             lhs = model.varDeviceFlow[dev, "el", "in", t]
-            electricity_carbon_capture = self.required_electricity(co2_flow_in, egr=egr, ccr=ccr)
+            electricity_carbon_capture = self.required_electricity(
+                co2_flow_in, ccr=ccr, specific_demand_MW_per_kgCO2=ccs_el_demand
+            )
             # compression:
-            carbon_captured = co2_flow_in * (1 - egr) * ccr
-            electricity_compressor = self.required_electricity_compressor(carbon_captured)
+            carbon_captured = co2_flow_in * ccr
+            electricity_compressor = self.required_electricity_compressor(carbon_captured, compressor_el_demand)
             rhs = electricity_carbon_capture + electricity_compressor
             return lhs == rhs
         elif i == 3:
             # carbon emitted
-            co2_emitted = co2_flow_in * (1 - egr) * (1 - ccr)
+            co2_emitted = co2_flow_in * (1 - ccr)
             lhs = model.varDeviceFlow[dev, "carbon", "out", t]
             rhs = co2_emitted
             return lhs == rhs
diff --git a/src/oogeso/core/devices/gasturbine.py b/src/oogeso/core/devices/gasturbine.py
index 41be67a..89b4ade 100644
--- a/src/oogeso/core/devices/gasturbine.py
+++ b/src/oogeso/core/devices/gasturbine.py
@@ -56,10 +56,20 @@ class GasTurbine(Device):
             rhs = (energy_in - model.varDeviceFlow[dev, "el", "out", t]) * self.dev_data.eta_heat
             return lhs == rhs
         elif i == 3:
-            """carbon emitted given by gas consumption"""
-            gasflow_co2 = gas_data.co2_content  # kg/m3
+            """carbon emitted given by gas consumption and EGR. The value is co2 after recirculation"""
+            gasflow_co2 = gas_data.co2_content  # kgCO2/Sm3gas
+            egr = self.dev_data.exhaust_gas_recirculation
+            if egr == 0:
+                y_egr = 1  # to not break with previous test cases
+            elif 0 <= egr <= 0.3:
+                y_egr = 1.411 * egr + 0.9838
+            elif 0.3 < egr <= 0.6:
+                y_egr = 3.5142 * egr + 0.3089
+            else:
+                raise ValueError("EGR parameter invalid.")
             lhs = model.varDeviceFlow[dev, "carbon", "out", t]
-            rhs = gasflow_co2 * (model.varDeviceFlow[dev, "gas", "in", t])
+            # multiply by (1-egr) to give flow after gas for recycling has been taken out:
+            rhs = gasflow_co2 * y_egr * (1 - egr) * (model.varDeviceFlow[dev, "gas", "in", t])
             return lhs == rhs
         elif (i == 4) and (has_hydrogen):
             """hydrogen blend max"""
@@ -86,3 +96,76 @@ class GasTurbine(Device):
 
     def get_flow_var(self, pyomo_model: pyo.Model, t: int):
         return pyomo_model.varDeviceFlow[self.id, "el", "out", t]
+
+
+class SteamCycle(Device):
+    """Steam cycle generator"""
+
+    carrier_in = ["heat"]
+    carrier_out = ["el", "heat"]
+    serial = ["heat"]
+
+    def __init__(self, *args, **kwargs):
+        super().__init__(*args, **kwargs)
+
+    #    def __init__(
+    #        self,
+    #        dev_data: dto.DeviceSteamCycleData,
+    #        carrier_data_dict: Dict[str, dto.CarrierData],
+    #    ):
+    #        super().__init__(dev_data=dev_data, carrier_data_dict=carrier_data_dict)
+    #        self.dev_data = dev_data
+    #        self.id = dev_data.id
+    #        self.carrier_data = carrier_data_dict
+
+    def _rules_misc(self, model: pyo.Model, t: int, i: int) -> Union[pyo.Expression, pyo.Constraint.Skip]:
+        dev = self.id
+        alpha = self.dev_data.alpha
+        egr = self.dev_data.exhaust_gas_recirculation
+        p_sc_nominal = self.dev_data.flow_max
+        if i == 1:
+            """power out vs heat in, expressed in normalised variables"""
+            # heat consumption is a linear function of el power output
+            # fuel = B + A*power
+            # => efficiency = power/(A+B*power)
+
+            # TODO make this user input
+            (linA, linB) = (1.96, -0.96)
+
+            heat_input_norm = model.varDeviceFlow[dev, "heat", "in", t] / p_sc_nominal
+            heat_extracted = model.varDeviceFlow[dev, "heat", "out", t]
+            el_output = model.varDeviceFlow[dev, "el", "out", t]
+
+            # normalised equation
+            y_heat = linA * heat_input_norm + linB * (model.varDeviceIsOn[dev, t] + model.varDeviceIsPrep[dev, t])
+            # egr dependency (piecewise linear factor):
+            if egr < 0.1:
+                y_egr = 0.1763 * egr + 1.0
+            else:
+                y_egr = 0.0225 * egr + 1.0153
+            power_computed = y_heat * y_egr * p_sc_nominal
+            # power equivalent of steam extracted for CCS (1/alpha >1 since 1 unit heat gives < 1 unit el)
+            power_extracted = heat_extracted / alpha
+            return el_output == power_computed - power_extracted
+        elif i == 2:
+            # TODO: this should not be necessary
+            # make sure that energy out < energy in
+            # heat out <= heat in
+            return pyo.Constraint.Skip
+            max_efficiency = 1.0
+            energy_in = model.varDeviceFlow[dev, "heat", "in", t]
+            energy_out = model.varDeviceFlow[dev, "heat", "out", t] + model.varDeviceFlow[dev, "el", "out", t]
+            return energy_out <= energy_in * max_efficiency
+        else:
+            return pyo.Constraint.Skip
+
+    def define_constraints(self, pyomo_model: pyo.Model):
+        """Specifies the list of constraints for the device"""
+        list_to_reconstruct = super().define_constraints(pyomo_model)
+
+        constr = pyo.Constraint(pyomo_model.setHorizon, pyo.RangeSet(1, 2), rule=self._rules_misc)
+        setattr(pyomo_model, "constr_{}_{}".format(self.id, "misc"), constr)
+        return list_to_reconstruct
+
+    def get_flow_var(self, pyomo_model: pyo.Model, t: int):
+        return pyomo_model.varDeviceFlow[self.id, "el", "out", t]
diff --git a/src/oogeso/core/devices/sink.py b/src/oogeso/core/devices/sink.py
index 3c8f7e6..991ba18 100644
--- a/src/oogeso/core/devices/sink.py
+++ b/src/oogeso/core/devices/sink.py
@@ -122,48 +122,6 @@ class SinkWater(Device):
         self.id = dev_data.id
         self.carrier_data = carrier_data_dict
 
-    def rule_devmodel_sink_water(self, model: pyo.Model, t: int, i: int):
-        dev = self.id
-        dev_data = self.dev_data
-        time_delta_minutes = model.paramTimestepDeltaMinutes
-
-        if dev_data.flow_avg is None:
-            return pyo.Constraint.Skip
-        if dev_data.max_accumulated_deviation is None:
-            return pyo.Constraint.Skip
-        if dev_data.max_accumulated_deviation == 0:
-            return pyo.Constraint.Skip
-        if i == 1:
-            # FLEXIBILITY
-            # (water_in-water_avg)*dt = delta buffer
-            delta_t = time_delta_minutes / 60  # hours
-            lhs = (model.varDeviceFlow[dev, "water", "in", t] - dev_data.flow_avg) * delta_t
-            if t > 0:
-                E_prev = model.varDeviceStorageEnergy[dev, t - 1]
-            else:
-                E_prev = model.paramDeviceEnergyInitially[dev]
-            rhs = model.varDeviceStorageEnergy[dev, t] - E_prev
-            return lhs == rhs
-        elif i == 2:
-            # energy buffer limit
-            E_max = dev_data.max_accumulated_deviation
-            return pyo.inequality(-E_max / 2, model.varDeviceStorageEnergy[dev, t], E_max / 2)
-        else:
-            raise Exception("impossible")
-
-    def define_constraints(self, pyomo_model: pyo.Model):
-        """Specifies the list of constraints for the device"""
-        list_to_reconstruct = super().define_constraints(pyomo_model)
-
-        constr = pyo.Constraint(
-            pyomo_model.setHorizon,
-            pyo.RangeSet(1, 2),
-            rule=self.rule_devmodel_sink_water,
-        )
-        # add constraints to model:
-        setattr(pyomo_model, "constr_{}_{}".format(self.id, "flex"), constr)
-        return list_to_reconstruct
-
     def get_flow_var(self, pyomo_model: pyo.Model, t: int):
         return pyomo_model.varDeviceFlow[self.id, "water", "in", t]
 
diff --git a/src/oogeso/core/devices/storage.py b/src/oogeso/core/devices/storage.py
index 698ddc3..07dd858 100644
--- a/src/oogeso/core/devices/storage.py
+++ b/src/oogeso/core/devices/storage.py
@@ -13,6 +13,21 @@ class StorageDevice(Device):
     def get_flow_var(self, pyomo_model: pyo.Model, t: int):
         pass
 
+    def compute_cost_for_depleted_storage(self, pyomo_model: pyo.Model, timesteps: pyo.Set):
+        # cost rate CURRENCY/s
+        # Cost associated with deviation from max value at end of horizon
+        dev_data = self.dev_data
+        if dev_data.E_cost is None:
+            return 0
+        if dev_data.E_max is None:
+            return 0
+        t_end = timesteps.at(-1)
+        var_E = pyomo_model.varDeviceStorageEnergy[self.id, t_end]
+        storage_cost = dev_data.E_cost * (dev_data.E_max - var_E)
+        # from cost (CURRENCY) to cost rate (CURRENCY/s):
+        # storage_cost = storage_cost / len(timesteps)  # 2024-11-01: No, don't divide - that means E_cost parameter must be changed if horizon length is changed.
+        return storage_cost
+
 
 class StorageEl(StorageDevice):
     """Electric storage (battery)"""
@@ -130,7 +145,7 @@ class StorageEl(StorageDevice):
         # even if Pmax=2 MW)
         return pyomo_model.varDeviceStoragePmax[self.id, t] + pyomo_model.varDeviceFlow[self.id, "el", "in", t]
 
-    def compute_cost_for_depleted_storage(self, pyomo_model: pyo.Model, timesteps: Optional[pyo.Set] = None):
+    def OBSOLETE_compute_cost_for_depleted_storage(self, pyomo_model: pyo.Model, timesteps: Optional[pyo.Set] = None):
         stor_cost = 0
         dev_data = self.dev_data
         if dev_data.E_cost is not None:
@@ -223,10 +238,10 @@ class StorageHydrogen(StorageDevice):
     def get_flow_var(self, pyomo_model: pyo.Model, t: int):
         return pyomo_model.varDeviceFlow[self.id, "hydrogen", "out", t]
 
-    def compute_cost_for_depleted_storage(self, pyomo_model: pyo.Model, timesteps: Optional[pyo.Set] = None):
+    def OBSOLETE_compute_cost_for_depleted_storage(self, pyomo_model: pyo.Model, timesteps: Optional[pyo.Set] = None):
         return self.compute_cost_for_depleted_storage_alt2(pyomo_model, timesteps)
 
-    def compute_cost_for_depleted_storage_alt1(self, pyomo_model: pyo.Model):
+    def OBSOLETE_compute_cost_for_depleted_storage_alt1(self, pyomo_model: pyo.Model):
         # cost if storage level at end of optimisation deviates from
         # target profile (user input based on expectations)
         # absolute value of deviation (filling too much also a cost)
@@ -237,7 +252,7 @@ class StorageHydrogen(StorageDevice):
         stor_cost = dev_data.E_cost * deviation
         return stor_cost
 
-    def compute_cost_for_depleted_storage_alt2(self, pyomo_model: pyo.Model, timesteps: pyo.Set):
+    def OBSOLOETE_compute_cost_for_depleted_storage_alt2(self, pyomo_model: pyo.Model, timesteps: pyo.Set):
         # cost rate kr/s
         # Cost associated with deviation from target value
         # below target = cost, above target = benefit   => gives signal to fill storage
@@ -326,3 +341,69 @@ class StorageGasLinepack(StorageDevice):
 
         else:
             raise ValueError(f"Argument i must be 1 or 2. {i} was given.")
+
+
+class WaterInjection(StorageDevice):
+    """
+    Water injection - flexible water sink
+    """
+
+    carrier_in = ["water"]
+    carrier_out = []
+    serial = []
+
+    def __init__(
+        self,
+        dev_data: dto.DeviceWaterInjectionData,
+        carrier_data_dict: Dict[str, dto.CarrierWaterData],
+    ):
+        super().__init__(dev_data=dev_data, carrier_data_dict=carrier_data_dict)
+        self.dev_data = dev_data
+        self.id = dev_data.id
+        self.carrier_data = carrier_data_dict
+
+    def rule_devmodel_sink_water(self, model: pyo.Model, t: int, i: int):
+        dev = self.id
+        dev_data = self.dev_data
+        time_delta_minutes = model.paramTimestepDeltaMinutes
+
+        if dev_data.flow_avg is None:
+            return pyo.Constraint.Skip
+        if dev_data.E_max is None:
+            return pyo.Constraint.Skip
+        if dev_data.E_max == 0:
+            return pyo.Constraint.Skip
+        if i == 1:
+            # FLEXIBILITY
+            # (water_in-water_avg)*dt = delta buffer
+            delta_t = time_delta_minutes / 60  # hours
+            lhs = (model.varDeviceFlow[dev, "water", "in", t] - dev_data.flow_avg) * delta_t
+            if t > 0:
+                E_prev = model.varDeviceStorageEnergy[dev, t - 1]
+            else:
+                E_prev = model.paramDeviceEnergyInitially[dev]
+            rhs = model.varDeviceStorageEnergy[dev, t] - E_prev
+            return lhs == rhs
+        elif i == 2:
+            # energy buffer limit
+            E_max = dev_data.E_max
+            E_min = dev_data.E_min
+            return pyo.inequality(E_min, model.varDeviceStorageEnergy[dev, t], E_max)
+        else:
+            raise Exception("impossible")
+
+    def define_constraints(self, pyomo_model: pyo.Model):
+        """Specifies the list of constraints for the device"""
+        list_to_reconstruct = super().define_constraints(pyomo_model)
+
+        constr = pyo.Constraint(
+            pyomo_model.setHorizon,
+            pyo.RangeSet(1, 2),
+            rule=self.rule_devmodel_sink_water,
+        )
+        # add constraints to model:
+        setattr(pyomo_model, "constr_{}_{}".format(self.id, "flex"), constr)
+        return list_to_reconstruct
+
+    def get_flow_var(self, pyomo_model: pyo.Model, t: int):
+        return pyomo_model.varDeviceFlow[self.id, "water", "in", t]
diff --git a/src/oogeso/core/devices/well.py b/src/oogeso/core/devices/well.py
index 6ba185d..d57e505 100644
--- a/src/oogeso/core/devices/well.py
+++ b/src/oogeso/core/devices/well.py
@@ -106,11 +106,12 @@ class WellGasLift(Device):
             else:
                 return pyo.Constraint.Skip  # noqa
         elif i == 4:
-            # gas injection pressure is fixed
-            if carrier == "gas":
+            # gas injection pressure must be larger than minimum value
+            # TODO: This constraint is not really needed, possible to specify pressure deviation bound on edge instead - all tests are ok
+            if (carrier == "gas") & (dev_data.injection_pressure is not None):
                 lhs = pyomo_model.varPressure[(node, carrier, "in", t)]
                 rhs = dev_data.injection_pressure
-                return lhs == rhs
+                return lhs >= rhs
             else:
                 return pyo.Constraint.Skip  # noqa
         else:
diff --git a/src/oogeso/core/optimiser.py b/src/oogeso/core/optimiser.py
index 4da14a3..8ad2eff 100644
--- a/src/oogeso/core/optimiser.py
+++ b/src/oogeso/core/optimiser.py
@@ -263,7 +263,7 @@ class OptimisationModel(pyo.ConcreteModel):
         self.varDeviceIsOn = pyo.Var(self.setDevice, self.setHorizon, within=pyo.Binary, initialize=1)
         self.varDeviceStarting = pyo.Var(self.setDevice, self.setHorizon, within=pyo.Binary, initialize=None)
         self.varDeviceStopping = pyo.Var(self.setDevice, self.setHorizon, within=pyo.Binary, initialize=None)
-        self.varDeviceStorageEnergy = pyo.Var(self.setDevice, self.setHorizon, within=pyo.Reals)
+        self.varDeviceStorageEnergy = pyo.Var(self.setDevice, self.setHorizon, within=pyo.Reals, initialize=0)
         # available reserve power from storage (linked to power rating and storage level):
         self.varDeviceStoragePmax = pyo.Var(self.setDevice, self.setHorizon, within=pyo.NonNegativeReals, initialize=0)
         # binary variable related to available powr from storage:
@@ -296,6 +296,7 @@ class OptimisationModel(pyo.ConcreteModel):
             self.setTerminal,
             self.setHorizon,
             within=pyo.Reals,
+            initialize=0,
         )
 
     def _specify_objective(self):
diff --git a/src/oogeso/dto/__init__.py b/src/oogeso/dto/__init__.py
index 7ab9f25..ffe0540 100644
--- a/src/oogeso/dto/__init__.py
+++ b/src/oogeso/dto/__init__.py
@@ -47,9 +47,11 @@ from oogeso.dto.devices import (
     DeviceSourceHydrogenData,
     DeviceSourceOilData,
     DeviceSourceWaterData,
+    DeviceSteamCycleData,
     DeviceStorageElData,
     DeviceStorageGasLinepackData,
     DeviceStorageHydrogenData,
+    DeviceWaterInjectionData,
     DeviceWellGasLiftData,
     DeviceWellProductionData,
 )
diff --git a/src/oogeso/dto/devices.py b/src/oogeso/dto/devices.py
index 05670b2..2f7984a 100644
--- a/src/oogeso/dto/devices.py
+++ b/src/oogeso/dto/devices.py
@@ -2,7 +2,7 @@ from typing import Dict, Optional
 
 from pydantic import Field
 
-from oogeso.dto import DeviceData
+from oogeso.dto import DeviceData, StartStopData
 from oogeso.dto.types import ModelType
 
 
@@ -72,11 +72,18 @@ class DeviceSinkOilData(DeviceData):
 
 class DeviceSinkWaterData(DeviceData):
     price: Dict[str, float] = Field(default_factory=lambda: {})
-    flow_avg: Optional[float] = None  # required average flow
-    max_accumulated_deviation: Optional[float] = None  # buffer size (max accumulated deviation from average)
     model: ModelType = ModelType.SINK_WATER
 
 
+class DeviceWaterInjectionData(DeviceData):
+    flow_avg: Optional[float] = None  # required average flow
+    target_profile: Optional[str] = None  # target profile
+    E_max: Optional[float] = None  # buffer size (max positive accumulated deviation from average)
+    E_min: Optional[float] = None  # buffer size (max negative accumulated deviation from average)
+    E_cost: float = 0  # cost for deviating from baseline
+    model: ModelType = ModelType.WATER_INJECTION
+
+
 class DeviceSinkCarbonData(DeviceData):
     price: Dict[str, float] = Field(default_factory=lambda: {})
     model: ModelType = ModelType.SINK_CARBON
@@ -119,6 +126,7 @@ class DeviceGasTurbineData(DeviceData):
     eta_heat: float = None
     hydrogen_blend_max: float = 0
     hydrogen_blend_min: float = 0
+    exhaust_gas_recirculation: float = 0  # egr
     reserve_factor: float = 1  # not used capacity contributes fully to spinning reserve
     model: ModelType = ModelType.GAS_TURBINE
 
@@ -166,6 +174,7 @@ class DeviceStorageHydrogenData(DeviceData):
 class DeviceStorageGasLinepackData(DeviceData):
     E_init: float = 0  # Sm3
     volume_m3: float = 0
+    E_cost: Optional[float] = None
     model: ModelType = ModelType.STORAGE_GAS_LINEPACK
 
 
@@ -185,6 +194,14 @@ class DeviceWellGasLiftData(DeviceData):
 
 class DeviceCarbonCaptureData(DeviceData):
     carbon_capture_rate: float = None  # ccr
-    exhaust_gas_recirculation: float = None  # egr
-    compressor_energy_demand: float = None  # MJ/kgCO2
+    capture_el_demand_MJ_per_kgCO2: float = None
+    capture_heat_demand_MJ_per_kgCO2: float = None
+    compressor_el_demand_MJ_per_kgCO2: float = None
     model: ModelType = ModelType.CARBON_CAPTURE
+
+
+class DeviceSteamCycleData(DeviceData):
+    exhaust_gas_recirculation: float = 0  # egr
+    alpha: float = None
+    model: ModelType = ModelType.STEAM_CYCLE
+    start_stop: StartStopData = StartStopData()
diff --git a/src/oogeso/dto/mapper.py b/src/oogeso/dto/mapper.py
index 70e481f..7c843ba 100644
--- a/src/oogeso/dto/mapper.py
+++ b/src/oogeso/dto/mapper.py
@@ -49,9 +49,11 @@ def get_device_from_model_name(model_name: str) -> Callable:
         "sourcehydrogen": devices.SourceHydrogen,
         "storagehydrogen": devices.StorageHydrogen,
         "storagegaslinepack": devices.StorageGasLinepack,
+        "waterinjection": devices.WaterInjection,
         "wellgaslift": devices.WellGasLift,
         "wellproduction": devices.WellProduction,
         "carboncapture": devices.CarbonCapture,
+        "steamcycle": devices.SteamCycle,
     }
     if model_name in map_device_name_to_class:
         return map_device_name_to_class[model_name]
@@ -89,9 +91,11 @@ def get_device_data_class_from_str(model_name: str) -> Callable:
         "sourcehydrogen": dto.DeviceSourceHydrogenData,
         "storagehydrogen": dto.DeviceStorageHydrogenData,
         "storagegaslinepack": dto.DeviceStorageGasLinepackData,
+        "waterinjection": dto.DeviceWaterInjectionData,
         "wellgaslift": dto.DeviceWellGasLiftData,
         "wellproduction": dto.DeviceWellProductionData,
         "carboncapture": dto.DeviceCarbonCaptureData,
+        "steamcycle": dto.DeviceSteamCycleData,
     }
     if model_name in map_device_name_to_class:
         return map_device_name_to_class[model_name]
diff --git a/src/oogeso/dto/types.py b/src/oogeso/dto/types.py
index 1b4eb47..2c7b5ac 100644
--- a/src/oogeso/dto/types.py
+++ b/src/oogeso/dto/types.py
@@ -29,6 +29,7 @@ class ModelType(str, Enum):
     SINK_OIL = "sinkoil"
     SINK_WATER = "sinkwater"
     SINK_CARBON = "sinkcarbon"
+    WATER_INJECTION = "waterinjection"
     COMPRESSOR_EL = "compressorel"
     COMPRESSOR_GAS = "compressorgas"
     ELECTROLYSER = "electrolyser"
@@ -46,6 +47,7 @@ class ModelType(str, Enum):
     WELL_PRODUCTION = "wellproduction"
     WELL_GAS_LIFT = "wellgaslift"
     CARBON_CAPTURE = "carboncapture"
+    STEAM_CYCLE = "steamcycle"
 
 
 class PressureMethodType(str, Enum):
diff --git a/src/oogeso/plots/plots.py b/src/oogeso/plots/plots.py
index 1ada435..f1e2c52 100644
--- a/src/oogeso/plots/plots.py
+++ b/src/oogeso/plots/plots.py
@@ -51,6 +51,7 @@ def plot_device_profile(
     include_prep=False,
     devs_shareload=None,
     plotter="plotly",
+    time_slice=None,
 ):
     """plot forecast and actual profile (available power), and device output
 
@@ -76,11 +77,13 @@ def plot_device_profile(
     """
     res = sim_result
     optimiser = optimisation_model
+    if time_slice is None:
+        time_slice = slice(None)
     if type(devs) is not list:
         devs = [devs]
     if include_forecasts & (len(devs) > 1):
         raise ValueError("Can only plot a single device when showing forecasts")
-    df = res.device_flow.unstack(["carrier", "terminal"])[("el", "out")].unstack("device")
+    df = res.device_flow.unstack(["carrier", "terminal"])[("el", "out")].unstack("device").loc[time_slice]
     if devs_shareload is None:
         # gas turbines:
         devs_shareload = [d for d, d_obj in optimiser.all_devices.items() if d_obj.dev_data.model == "gasturbine"]
@@ -93,14 +96,18 @@ def plot_device_profile(
             df.loc[mask, c] = devs_mean[mask]
     df.columns.name = "devices"
     df = df[devs]
+    if df.min().min() < 0:
+        print(f"Setting negative values to zero. (min={df.min().min():g})")
+        df = df.clip(lower=0)
     nrows = 1
     if include_on_off:
         nrows = nrows + 1
     if include_prep:
         nrows = nrows + 1
-    df2 = res.device_is_on.unstack("device")[devs]
-    df_prep = res.device_is_prep.unstack("device")[devs]
-    timerange = (res.device_flow.index.get_level_values("time")).unique()
+    df2 = res.device_is_on.unstack("device")[devs].loc[time_slice]
+    df_prep = res.device_is_prep.unstack("device")[devs].loc[time_slice]
+    # timerange = (res.device_flow.index.get_level_values("time")).unique()
+    timerange = df.index  # time_slice
     if plotter == "plotly":
         fig = plotly.subplots.make_subplots(rows=nrows, cols=1, shared_xaxes=True)
         colour = pio.templates[pio.templates.default].layout.colorway
@@ -208,7 +215,7 @@ def plot_device_profile(
                     df_forecast[dev + " forecast"] = (
                         res.profiles_forecast.loc[timerange, dev_data.profile] * dev_data.flow_max
                     )
-            df_forecast.plot.area(zorder=0, linewidth=0, ax=axs[0])
+            df_forecast.plot(zorder=0, ax=axs[0], linestyle="-")
         tmin = df.index.get_level_values("time").min()
         tmax = df.index.get_level_values("time").max() + 1
         ax.set_xlim(tmin, tmax)
@@ -701,6 +708,8 @@ def plot_reserve(
     include_sum=True,
     devs_shareload=None,
     plotter="plotly",
+    time_slice=None,
+    return_devs=False,
 ):
     """Plot unused online capacity by all el devices
     devs_shareload : list ([]=ignore, None=do it for gas turbines)
@@ -711,9 +720,11 @@ def plot_reserve(
         this is included in the plots.
     """
     df_devs = pd.DataFrame()
+    if time_slice is None:
+        time_slice = slice(None)
     res = sim_result
     optimiser = optimisation_model
-    timerange = list(res.el_reserve.index)
+    timerange = list(res.el_reserve.index[time_slice])
     margin_incr = pd.DataFrame(0, index=timerange, columns=["margin"])
     if devs_shareload is None:
         # split load evently amongst gas turbines:
@@ -760,6 +771,10 @@ def plot_reserve(
             reserv = cap_avail - p_generating
             df_devs[d] = reserv
     df_devs.columns.name = "device"
+    if df_devs.min().min() < 0:
+        print(f"Setting negative values to zero. (min={df_devs.min().min():g})")
+        df_devs = df_devs.clip(lower=0)
+    df_devs = df_devs.loc[time_slice]
     if plotter == "plotly":
         fig = px.area(df_devs, line_shape="hv")
         if include_sum:
@@ -799,7 +814,10 @@ def plot_reserve(
         fig = plt.gcf()
     else:
         raise ValueError(f"Plotter: {plotter} has not been implemented for plot reserve.")
-    return fig
+    if return_devs:
+        return fig, df_devs
+    else:
+        return fig
 
 
 def plot_el_backup(sim_result, filename=None, show_margin=False, return_margin=False, plotter="plotly"):


Test patch:
diff --git a/tests/test_carboncapture.py b/tests/test_carboncapture.py
index ff99c9f..b163f55 100644
--- a/tests/test_carboncapture.py
+++ b/tests/test_carboncapture.py
@@ -19,30 +19,28 @@ def testcase_carboncapture_data() -> oogeso.dto.base.EnergySystemData:
 def test_carboncapture_model():
     # dev_obj = CarbonCapture()
 
-    egr = 0.4
+    # egr = 0.4
     ccr = 0.9
     flow_in_co2 = 1
-    frac_co2 = CarbonCapture.co2frac(egr=egr)
-    specific_heat_demand = CarbonCapture.specific_heat_demand(frac_co2=frac_co2)
-    specific_el_demand = CarbonCapture.specific_electricity_demand(frac_co2=frac_co2)
-    heat_demand = CarbonCapture.required_heat(flow_in_co2=1, egr=egr, ccr=ccr)
-    el_demand = CarbonCapture.required_electricity(flow_in_co2=flow_in_co2, egr=egr, ccr=ccr)
-    co2_captured = flow_in_co2 * (1 - egr) * ccr
-
-    assert frac_co2 == (12.88 * egr + 0.7874) / 100  # = 0.059394
-    if frac_co2 < 0.0489:
-        assert specific_heat_demand == -3.57 * frac_co2 + 4.1551
-        assert specific_el_demand == 0.16 * frac_co2 + 0.0232
-    else:
-        assert specific_heat_demand == -2.11 * frac_co2 + 4.0837
-        assert specific_el_demand == 0.009 * frac_co2 + 0.0305
-    assert heat_demand == co2_captured * specific_heat_demand
-    assert el_demand == co2_captured * specific_el_demand
+    # frac_co2 = CarbonCapture.co2frac(egr=egr)
+    # specific_heat_demand = CarbonCapture.specific_heat_demand(frac_co2=frac_co2)
+    # specific_el_demand = CarbonCapture.specific_electricity_demand(frac_co2=frac_co2)
+    heat_demand = CarbonCapture.required_heat(flow_in_co2=1, ccr=ccr, specific_demand_MW_per_kgCO2=0.4)
+    el_demand = CarbonCapture.required_electricity(
+        flow_in_co2=flow_in_co2, ccr=ccr, specific_demand_MW_per_kgCO2=0.032
+    ) + CarbonCapture.required_electricity_compressor(flow_in_co2=flow_in_co2, energy_demand_MJ_per_kg=0.3)
+    capture_heat_demand_MJ_per_kgCO2 = 0.4  # MJ/kgCO2
+    capture_el_demand_MJ_per_kgCO2 = 0.032  # MJ/kgCO2
+    compressor_el_demand_MJ_per_kgCO2 = 0.30  # MJ/kgCO2
+    # co2_captured = flow_in_co2 * (1 - egr) * ccr
+
+    assert el_demand == (capture_el_demand_MJ_per_kgCO2 * ccr + compressor_el_demand_MJ_per_kgCO2) * flow_in_co2
+    assert heat_demand == capture_heat_demand_MJ_per_kgCO2 * flow_in_co2 * ccr
 
 
 @pytest.mark.skipif(not pyo.SolverFactory("cbc").available(), reason="Skipping test because CBC is not available.")
 def test_carboncapture_simple_simulation(testcase_carboncapture_data: dto.EnergySystemData):
     my_sim = oogeso.Simulator(testcase_carboncapture_data)
     res = my_sim.run_simulation("cbc", time_range=(0, 4))
-    assert (res.device_flow["ccs", "carbon", "in"] == 2.34).all()
-    assert (res.device_flow["ccs", "carbon", "out"] == 0.1404).all()
+    assert res.device_flow["ccs", "carbon", "in", 0] == pytest.approx(24.04744, rel=1e-4)
+    assert res.device_flow["ccs", "carbon", "out", 0] == pytest.approx(2.40474, rel=1e-4)
diff --git a/tests/test_data/testdata_carboncapture.yaml b/tests/test_data/testdata_carboncapture.yaml
index f37287b..95972ab 100644
--- a/tests/test_data/testdata_carboncapture.yaml
+++ b/tests/test_data/testdata_carboncapture.yaml
@@ -3,11 +3,11 @@ parameters:
   planning_horizon: 2 #8 # timesteps in each rolling optimisation
   optimisation_timesteps: 1 #4 # timesteps between each optimisation
   forecast_timesteps: 1 #4 # timesteps beyond which forecast (instead of nowcast) profile is used
-  objective: penalty # select objective function to use
+  objective: costs # select objective function to use
 
 carriers:
   - id: el
-    el_reserve_margin: 2 # MWm -1=no limit
+    el_reserve_margin: -1 #2 # MWm -1=no limit
     reserve_storage_minutes: 30 # minutes, how long stored energy must be sustained to count as reserve
     powerflow_method: "transport"
   - id: gas
@@ -24,6 +24,7 @@ carriers:
 
 nodes:
   - id: gt
+  - id: sc
   - id: ccs
   - id: sink
 
@@ -40,17 +41,29 @@ edges:
     node_from: ccs
     node_to: sink
 
+  - id: el3
+    include: 1
+    carrier: el
+    node_from: sc
+    node_to: ccs
+
   - id: heat1
     include: 1
     carrier: heat
     node_from: gt
-    node_to: ccs
+    node_to: sink
 
   - id: heat2
     include: 1
     carrier: heat
-    node_from: ccs
-    node_to: sink
+    node_from: sc
+    node_to: ccs
+
+  - id: heat3
+    include: 1
+    carrier: heat
+    node_from: gt
+    node_to: sc
 
   - id: c1
     include: 1
@@ -78,26 +91,38 @@ devices:
     node_id: gt
     model: gasturbine
     name: GTG
-    fuel_A: 2.35
-    fuel_B: 0.531
-    eta_heat: 0.5    
-    flow_max: 20
+    fuel_A: 3.6
+    fuel_B: 0.3
+    eta_heat: 1    
+    flow_max: 209
     flow_min: 0
+    exhaust_gas_recirculation: 0.0
+
+  - id: sc
+    include: 1
+    node_id: sc
+    model: steamcycle
+    name: "steam cycle"
+    alpha: 0.238
+    flow_max: 93
+    exhaust_gas_recirculation: 0.0
 
   - id: ccs
     include: 1
     node_id: ccs
     model: carboncapture
     name: CCS
-    exhaust_gas_recirculation: 0.4
     carbon_capture_rate: 0.9
-    compressor_energy_demand: 0.30 # MJ/kgCO2
+    capture_heat_demand_MJ_per_kgCO2: 0.4 # MJ/kgCO2
+    capture_el_demand_MJ_per_kgCO2: 0.032 # MJ/kgCO2
+    compressor_el_demand_MJ_per_kgCO2: 0.30 # MJ/kgCO2
 
   - id: el_dump
     include: True
     node_id: sink
     model: sink_el
     name: "el sink"
+    flow_min: 200
 
   - id: heat_dump
     include: True
@@ -110,3 +135,4 @@ devices:
     node_id: sink
     model: sink_carbon
     name: "CO2 emitted"
+    op_cost: 5
diff --git a/tests/test_devices_constraints.py b/tests/test_devices_constraints.py
index bbd7df6..d31decf 100644
--- a/tests/test_devices_constraints.py
+++ b/tests/test_devices_constraints.py
@@ -334,9 +334,11 @@ def test_well_production_constraints():
 def test_carbon_capture_constraints():
     dev_data = dto.DeviceCarbonCaptureData(
         **dev_data_generic,
-        exhaust_gas_recirculation=0.3,
+        # exhaust_gas_recirculation=0.3,
         carbon_capture_rate=0.5,
-        compressor_energy_demand=0.30,
+        compressor_el_demand_MJ_per_kgCO2=0.30,
+        capture_el_demand_MJ_per_kgCO2=0,
+        capture_heat_demand_MJ_per_kgCO2=0
     )
     optimisation_model = _build_lp_problem_with_single_dev(dev_data)
     assert isinstance(optimisation_model, OptimisationModel)
diff --git a/tests/test_devices_create.py b/tests/test_devices_create.py
index e639ca3..d0915c5 100644
--- a/tests/test_devices_create.py
+++ b/tests/test_devices_create.py
@@ -156,13 +156,22 @@ def test_sink_oil():
 
 
 def test_sink_water():
-    dev_data = dto.DeviceSinkWaterData(**dev_data_generic, flow_avg=None, max_accumulated_deviation=None)
+    dev_data = dto.DeviceSinkWaterData(**dev_data_generic)
     carrier_data_dict = {}
     obj = devices.SinkWater(dev_data, carrier_data_dict)
     assert isinstance(obj, devices.SinkWater)
     assert obj.dev_data.flow_max == 20
 
 
+def test_water_injection():
+    dev_data = dto.DeviceWaterInjectionData(**dev_data_generic, flow_avg=None, E_max=10, E_min=30)
+    carrier_data_dict = {}
+    obj = devices.WaterInjection(dev_data, carrier_data_dict)
+    assert isinstance(obj, devices.WaterInjection)
+    assert obj.dev_data.flow_max == 20
+    assert obj.dev_data.E_min == 30
+
+
 def test_sink_heat():
     dev_data = dto.DeviceSinkHeatData(**dev_data_generic)
     carrier_data_dict = {}
@@ -301,13 +310,21 @@ def test_carbon_capture():
     dev_data = dto.DeviceCarbonCaptureData(
         **dev_data_generic,
         carbon_capture_rate=0.9,
-        exhaust_gas_recirculation=0.3,
-        compressor_energy_demand=0.3,
+        capture_heat_demand_MJ_per_kgCO2=0.4,  # MJ/kgCO2
+        capture_el_demand_MJ_per_kgCO2=0.032,  # MJ/kgCO2
+        compressor_el_demand_MJ_per_kgCO2=0.30,  # MJ/kgCO2
     )
     carrier_data_dict = {}
     obj = devices.CarbonCapture(dev_data, carrier_data_dict)
     assert isinstance(obj, devices.CarbonCapture)
     assert obj.dev_data.flow_max == 20
-    assert obj.dev_data.compressor_energy_demand == 0.3
+    assert obj.dev_data.compressor_el_demand_MJ_per_kgCO2 == 0.3
     assert obj.carrier_in == ["carbon", "el", "heat"]
     assert obj.carrier_out == ["carbon"]
+
+
+def test_steamcycle():
+    dev_data = dto.DeviceSteamCycleData(**dev_data_generic, exhaust_gas_recirculation=0.3, alpha=0.5)
+    obj = devices.SteamCycle(dev_data, {})
+    assert isinstance(obj, devices.SteamCycle)
+    assert obj.dev_data.exhaust_gas_recirculation == 0.3