Merge pull request #96 from SPF-OST/solar_ice_example

Solar ice example

data/examples/solar_ice_hp.py

@@ -6,14 +6,8 @@
 from optihood.energy_network import EnergyNetworkGroup as EnergyNetwork
 
 if __name__ == '__main__':
-    """
-    
-    This example is currently under development (use with caution) !
-    
-    """
-
     # set a time period for the optimization problem
-    timePeriod = pd.date_range("2018-01-01 00:00:00", "2018-01-31 23:00:00", freq="60min")
+    timePeriod = pd.date_range("2018-06-01 00:00:00", "2018-12-31 23:00:00", freq="60min")
 
     # define paths for input and result files
     curDir = _pl.Path(__file__).resolve().parent
@@ -29,9 +23,9 @@
     # initialize parameters
     numberOfBuildings = 4
     optimizationType = "costs"  # set as "env" for environmental optimization and "costs" for cost optimization
-    mergeLinkBuses = True
-    mergeBuses = ["electricity", "heat_buses", "heatPumpInputBus"]
-    dispatchMode = False  # Set to True to run the optimization in dispatch mode
+    mergeLinkBuses = False
+    mergeBuses = [] # "electricity", "heatPumpInputBus"
+    dispatchMode = True  # Set to True to run the optimization in dispatch mode
 
     # solver specific command line options
     optimizationOptions = {

optihood/buildings.py

@@ -107,13 +107,13 @@ def addToBusDict(self, busDictBuilding1, mergeType):
 
     def linkBuses(self, busesToMerge):
         for b in busesToMerge:
-            if b=="electricity":
+            if b == "electricity":
                 self.__linkBuses.extend(["electricityBus", "electricityInBus"])
-            elif b=="space_heat":
+            elif b == "space_heat":
                 self.__linkBuses.extend(["spaceHeatingBus", "shDemandBus"])
-            elif b=="domestic_hot_water":
+            elif b == "domestic_hot_water":
                 self.__linkBuses.extend(["domesticHotWaterBus", "dhwDemandBus"])
-            elif b=="heat_buses":
+            elif b == "heat_buses":
                 self.__linkBuses.extend(["heatDemandBus0", "heatDemandBus2"])
             else:
                 self.__linkBuses.append(b)
@@ -885,7 +885,7 @@ def addStorage(self, data, storageParams, ambientTemperature, opt, mergeLinkBuse
                                                        tStorInit=float(storageParams["ice_storage"].at[s["label"],"intitial_temp"]),
                                                        fMax=float(storageParams["ice_storage"].at[s["label"],"max_ice_fraction"]),
                                                        rho=float(storageParams["ice_storage"].at[s["label"],"rho_fluid"]),
-                                                       V=float(s["capacity max"]),
+                                                       V=float(s["capacity max"])/1000, # conversion from L to m3
                                                        hfluid=float(storageParams["ice_storage"].at[s["label"],"h_fluid"]),
                                                        cp=float(storageParams["ice_storage"].at[s["label"],"cp_fluid"]),
                                                        Tamb=ambientTemperature,

optihood/combined_prod.py

@@ -41,40 +41,24 @@ def _create(self, group=None):
         m = self.parent_block()
 
         for n in group:
-            n.inflow = [i for i in list(n.inputs) if "electricity" in i.label or "naturalGasBus" in i.label][0]
+            n.inflow = [i for i in list(n.inputs) if "electricity" in i.label or list(n.inputs).__len__()==1][0]
             if len(list(n.inputs)) > 1:
                 n.inflowQevap = [i for i in list(n.inputs) if "electricity" not in i.label][0]
             else:
                 n.inflowQevap = 0
             flows = [k for k, v in n.efficiency.items()]
-            n.flowT0 = flows[0]
-            n.flowT1 = flows[1]
-            n.outputT0 = [o for o in n.outputs if n.flowT0 == o][0]
-            n.outputT1 = [o for o in n.outputs if n.flowT1 == o][0]
-            n.efficiency_sq = (
-                n.efficiency[n.outputT0],
-                n.efficiency[n.outputT1]
-            )
-            if len(flows)==3:
-                n.flowT2 = flows[2]
-                n.outputT2 = [o for o in n.outputs if n.flowT2 == o][0]
-                n.efficiency_sq = (
-                    n.efficiency[n.outputT0],
-                    n.efficiency[n.outputT1],
-                    n.efficiency[n.outputT2])
+            n.outputs_ordered = []
+            n.efficiency_sq = ()
+            for i in range(len(flows)):
+                n.outputs_ordered.append([o for o in n.outputs if flows[i] == o][0])
+                n.efficiency_sq = n.efficiency_sq + (n.efficiency[n.outputs_ordered[i]],)
 
         def _input_output_relation_rule(block):
             """Connection between input and outputs."""
             for t in m.TIMESTEPS:
                 for g in group:
                     lhs = m.flow[g.inflow, g, t]
-                    if len(g.efficiency_sq)==3:
-                        rhs = (m.flow[g, g.outputT0, t] / g.efficiency_sq[0][t]
-                               + m.flow[g, g.outputT1, t] / g.efficiency_sq[1][t]
-                               + m.flow[g, g.outputT2, t] / g.efficiency_sq[2][t])
-                    else:
-                        rhs = (m.flow[g, g.outputT0, t] / g.efficiency_sq[0][t]
-                               + m.flow[g, g.outputT1, t] / g.efficiency_sq[1][t])
+                    rhs = sum(m.flow[g, g.outputs_ordered[i], t] / g.efficiency_sq[i][t] for i in range(len(g.efficiency_sq)))
                     block.input_output_relation.add((g, t), (lhs == rhs))
 
         self.input_output_relation = Constraint(
@@ -90,7 +74,7 @@ def _second_input_relation_rule(block):
                 for g in group:
                     if len(list(g.inputs)) > 1:
                         lhs = m.flow[g.inflowQevap, g, t]
-                        rhs = (m.flow[g, g.outputT0, t] + m.flow[g, g.outputT1, t] + m.flow[g, g.outputT2, t] - m.flow[g.inflow, g, t])
+                        rhs = (sum(m.flow[g, g.outputs_ordered[i], t] for i in range(len(g.efficiency_sq))) - m.flow[g.inflow, g, t])
                         block.input_relation.add((g, t), (lhs == rhs))
 
         self.input_relation = Constraint(

optihood/constraints.py

@@ -66,29 +66,29 @@ def connectInvestmentRule(om):
     dhwLinkOutputFlows = [(i, o) for (i, o) in om.flows if i.label == "dhwLink"]
 
     if elLinkOutputFlows:
-        first = om.InvestmentFlow.invest[next(iter(elLinkOutputFlows))]
+        first = om.InvestmentFlowBlock.invest[next(iter(elLinkOutputFlows))]
         for (i, o) in elLinkOutputFlows:
-            expr = (first == om.InvestmentFlow.invest[i, o])
+            expr = (first == om.InvestmentFlowBlock.invest[i, o])
             setattr(
                 om,
                 "elLinkConstr_" + o.label,
                 pyo.Constraint(expr=expr),
             )
 
     if shLinkOutputFlows:
-        first = om.InvestmentFlow.invest[next(iter(shLinkOutputFlows))]
+        first = om.InvestmentFlowBlock.invest[next(iter(shLinkOutputFlows))]
         for (i, o) in shLinkOutputFlows:
-            expr = (first == om.InvestmentFlow.invest[i, o])
+            expr = (first == om.InvestmentFlowBlock.invest[i, o])
             setattr(
                 om,
                 "shLinkConstr_" + o.label,
                 pyo.Constraint(expr=expr),
             )
 
     if dhwLinkOutputFlows:
-        first = om.InvestmentFlow.invest[next(iter(dhwLinkOutputFlows))]
+        first = om.InvestmentFlowBlock.invest[next(iter(dhwLinkOutputFlows))]
         for (i, o) in dhwLinkOutputFlows:
-            expr = (first == om.InvestmentFlow.invest[i, o])
+            expr = (first == om.InvestmentFlowBlock.invest[i, o])
             setattr(
                 om,
                 "dhwLinkConstr_" + o.label,
@@ -228,12 +228,12 @@ def electricRodCapacityConstaint(om, numBuildings):
     groundHeatPumpCapacityTotal = 0
 
     for b in range(1,numBuildings+1):
-        elRodCapacity = [om.InvestmentFlow.invest[i, o] for (i, o) in electricRodInputFlows if ((f'__Building{b}') in o.label)]
-        airHeatPumpCapacity = [om.InvestmentFlow.invest[i, o] for (i, o) in airHeatPumpInputFlows if ((f'__Building{b}') in o.label)]
+        elRodCapacity = [om.InvestmentFlowBlock.invest[i, o] for (i, o) in electricRodInputFlows if ((f'__Building{b}') in o.label)]
+        airHeatPumpCapacity = [om.InvestmentFlowBlock.invest[i, o] for (i, o) in airHeatPumpInputFlows if ((f'__Building{b}') in o.label)]
         if groundHeatPumpInputFlows:
-            groundHeatPumpCapacity = [om.InvestmentFlow.invest[i, o] for (i, o) in groundHeatPumpInputFlows if ((f'__Building{b}') in o.label)]
+            groundHeatPumpCapacity = [om.InvestmentFlowBlock.invest[i, o] for (i, o) in groundHeatPumpInputFlows if ((f'__Building{b}') in o.label)]
         else:
-            groundHeatPumpCapacity = [om.InvestmentFlow.invest[i, o] for (i, o) in groundHeatPumpOutFlows if ((f'__Building{b}') in i.label)]
+            groundHeatPumpCapacity = [om.InvestmentFlowBlock.invest[i, o] for (i, o) in groundHeatPumpOutFlows if ((f'__Building{b}') in i.label)]
         if elRodCapacity:
             elRodCapacity = elRodCapacity[0]
             airHeatPumpCapacity = airHeatPumpCapacity[0] if airHeatPumpCapacity else 0
@@ -257,7 +257,7 @@ def totalPVCapacityConstraint(om, numBuildings):
     pvOutFlows = [(i, o) for (i, o) in om.flows if ("pv" in i.label)]
     pvCapacityTotal = 0
     for b in range(1,numBuildings+1):
-        pvCapacity = [om.InvestmentFlow.invest[i, o] for (i, o) in pvOutFlows if ((f'__Building{b}') in o.label)]
+        pvCapacity = [om.InvestmentFlowBlock.invest[i, o] for (i, o) in pvOutFlows if ((f'__Building{b}') in o.label)]
         if pvCapacity:
             pvCapacity = pvCapacity[0]
             pvCapacityTotal = pvCapacityTotal + pvCapacity

optihood/energy_network.py

@@ -1404,9 +1404,13 @@ def _addLinks(self, data, numberOfBuildings, mergeLinkBuses):  # connects buses
                     elif "heat" in l["label"]:
                         busesOut.append(self._busDict["heatBus" + l["label"][-1] + '__Building' + str(b + 1)])
                         busesIn.append(self._busDict["heatDemandBus" + l["label"][-1] + '__Building' + str(b + 1)])
-                    else:
+                    elif "electricity" in l["label"]:
                         busesOut.append(self._busDict["electricityBus" + '__Building' + str(b + 1)])
                         busesIn.append(self._busDict["electricityInBus" + '__Building' + str(b + 1)])
+                    elif "dhLink" in l["label"]:
+                        if "districtHeatingInputBus" + '__Building' + str(b + 1) in self._busDict:
+                            busesOut.append(self._busDict["districtHeatingInputBus"  + '__Building' + str(b + 1)])
+                        busesIn.append(self._busDict["districtHeatingBus" + '__Building' + str(b + 1)])
 
                 self._nodesList.append(Link(
                     label=l["label"],

optihood/storages.py

@@ -297,20 +297,19 @@ class IceStorage(on.Node):
     fMax: maximum allowed ice fraction in the ice storage
     rho: density of water [kg/m3]
     V: Volume of the ice storage [m3]
-    hfluid: []
-    cp: specific heat capacity of water []
+    hfluid: latent heat of fusion of water [kJ/kg]
+    cp: specific heat capacity of water [kJ/(kg.°C)]
     T_amb: ambient temperature time-series [°C]
-    UA_tank: heat loss coefficient of the ice storage tank []
+    UA_tank: heat loss coefficient of the ice storage tank [kW/°C]
     inflow_conversion_factor: efficiency of heat exchanger at the inlet of the ice storage
     outflow_conversion_factor: efficiency of heat exchanger at the outlet of the ice storage
     """
     def __init__(self, label, input, output, tStorInit, fMax, rho, V, hfluid, cp, Tamb, UAtank, inflow_conversion_factor, outflow_conversion_factor):
         if tStorInit <= 0:
             raise ValueError("The initial temperature of ice storage should be greater than 0°C, i.e. without any ice formation.")
         self.tStorInit = tStorInit
-        if fMax < 0.5 or fMax >0.8:
-            raise ValueError(
-                "Maximum ice fraction should be defined within the range 0.5-0.8.")
+        if fMax < 0.5 or fMax > 0.8:
+            raise ValueError("Maximum ice fraction should be defined within the range 0.5-0.8")
         self.fMax = fMax
         self.massWaterMax = rho*V
         self.hf = hfluid
@@ -374,19 +373,22 @@ def _create(self, group=None):
         o = {n: [o for o in n.outputs][0] for n in group}
 
         #  ************* SET OF ICE STORAGES *****************************
-        self.ICESTORAGES = Set(initizalize=[n for n in group])
+        self.icestorages = Set(initialize=[n for n in group])
 
         #  ************* DECISION VARIABLES *****************************
         # temperature of ice storage
-        self.tStor = Var(self.ICESTORAGES, m.TIMESTEPS, within=NonNegativeReals, bounds=(0,200))
-        self.tStor_prev = Var(self.ICESTORAGES, m.TIMESTEPS, within=NonNegativeReals, bounds=(0,200))
+        self.tStor = Var(self.icestorages, m.TIMESTEPS, within=NonNegativeReals, bounds=(0,200))
+        self.tStor_prev = Var(self.icestorages, m.TIMESTEPS, within=NonNegativeReals, bounds=(0,200))
         # mass of ice in storage
-        self.mIceStor = Var(self.ICESTORAGES, m.TIMESTEPS, within=NonNegativeReals, bounds=(0,1000))
-        self.mIceStor_prev = Var(self.ICESTORAGES, m.TIMESTEPS, within=NonNegativeReals, bounds=(0,1000))
+        self.mIceStor = Var(self.icestorages, m.TIMESTEPS, within=NonNegativeReals, bounds=(0,100000))
+        self.mIceStor_prev = Var(self.icestorages, m.TIMESTEPS, within=NonNegativeReals, bounds=(0,100000))
         # ice fraction
-        self.fIce = Var(self.ICESTORAGES, m.TIMESTEPS, within=NonNegativeReals, bounds=(0,self.fMax))
+        self.fIce = Var(self.icestorages, m.TIMESTEPS, within=NonNegativeReals, bounds=(0,1))
         # binary variable defining the status of ice formation in the storage
-        self.iceStatus = Var(self.ICESTORAGES, m.TIMESTEPS, within=Binary)
+        self.iceStatus = Var(self.icestorages, m.TIMESTEPS, within=Binary)
+
+        # for linearization of non-linear constraints with big M method
+        M = 2000
 
         #  ************* CONSTRAINTS *****************************
 
@@ -423,12 +425,23 @@ def _prev_temperature_rule(block):
         self.prev_temperature = Constraint(group, m.TIMESTEPS, noruleinit=True)
         self.prev_temperature_build = BuildAction(rule=_prev_temperature_rule)
 
+        def _max_ice_fraction_rule(block):
+            """set the value of ice fraction"""
+            for g in group:
+                for t in m.TIMESTEPS:
+                    lhs = self.fIce[g,t]
+                    rhs = g.fMax
+                    block.max_ice_fraction.add((g, t), (lhs <= rhs))
+
+        self.max_ice_fraction = Constraint(group, m.TIMESTEPS, noruleinit=True)
+        self.max_ice_fraction_build = BuildAction(rule=_max_ice_fraction_rule)
+
         def _ice_fraction_rule(block):
             """set the value of ice fraction"""
             for g in group:
                 for t in m.TIMESTEPS:
                     lhs = self.fIce[g,t]
-                    rhs = self.mIceStor[g,t]/g.massWaterMax
+                    rhs = self.mIceStor[g,t]/g.massWaterMax     # division is most likely not allowed!!
                     block.ice_fraction.add((g, t), (lhs == rhs))
 
         self.ice_fraction = Constraint(group, m.TIMESTEPS, noruleinit=True)
@@ -439,11 +452,14 @@ def _storage_balance_rule(block):
             for g in group:
                 for t in m.TIMESTEPS:
                     expr = 0
-                    expr += (g.rho*g.V*g.cp*(self.tStor[g,t] - self.tStor_prev[g,t]))
-                    expr += (g.UAtank*(self.tStor_prev[g,t] - g.Tamb[t])*m.timeincrement[t])
-                    expr += (-g.hf*(self.mIceStor[g,t] - self.mIceStor_prev[g,t]))
-                    expr += (-m.flow[i[g], g, t]*g.inflow_conversion_factor[t]*m.timeincrement[t])
-                    expr += ((m.flow[g, o[g], t]/g.outflow_conversion_factor[t])*m.timeincrement[t])
+                    expr += g.rho*g.V*g.cp*self.tStor[g,t]
+                    expr += - g.rho * g.V * g.cp * self.tStor_prev[g, t]
+                    expr += g.UAtank*self.tStor_prev[g,t]*m.timeincrement[t]
+                    expr += - g.UAtank * g.Tamb[t] * m.timeincrement[t]
+                    expr += - g.hf*self.mIceStor[g,t]
+                    expr += g.hf * self.mIceStor_prev[g, t]
+                    expr += (-m.flow[i[g], g, t]*g.inflow_conversion_factor*m.timeincrement[t])
+                    expr += ((m.flow[g, o[g], t]/g.outflow_conversion_factor)*m.timeincrement[t])
                     block.storage_balance.add((g, t), (expr == 0))
 
         self.storage_balance = Constraint(group, m.TIMESTEPS, noruleinit=True)
@@ -455,30 +471,41 @@ def _mass_ice_rule(block):
                 for t in m.TIMESTEPS:
                     lhs = self.mIceStor[g, t]
                     rhs = self.iceStatus[g, t]*(self.mIceStor_prev[g, t] +
-                                                (((m.flow[g, o[g], t]/g.outflow_conversion_factor[t])*m.timeincrement[t]
-                                                  - m.flow[i[g], g, t]*g.inflow_conversion_factor[t]*m.timeincrement[t]
+                                                (((m.flow[g, o[g], t]/g.outflow_conversion_factor)*m.timeincrement[t]
+                                                  - m.flow[i[g], g, t]*g.inflow_conversion_factor*m.timeincrement[t]
                                                   + g.UAtank*(self.tStor_prev[g, t] - g.Tamb[t]) * m.timeincrement[t]
                                                   - g.rho*g.V*g.cp*self.tStor_prev[g, t])/g.hf))
                     block.mass_ice.add((g, t), (lhs == rhs))
 
         self.mass_ice = Constraint(group, m.TIMESTEPS, noruleinit=True)
         self.mass_ice_build = BuildAction(rule=_mass_ice_rule)
 
-        def _ice_state_rule(block):
+        def _ice_state_rule_1(block):
+            """rule for calculating the mass of ice in each timestep"""
+            for g in group:
+                for t in m.TIMESTEPS:
+                    lhs = self.iceStatus[g,t]
+                    rhs = self.tStor[g,t]
+                    block.ice_state_1.add((g, t), (lhs <= rhs))
+
+        self.ice_state_1 = Constraint(group, m.TIMESTEPS, noruleinit=True)
+        self.ice_state_1_build = BuildAction(rule=_ice_state_rule_1)
+
+        def _ice_state_rule_2(block):
             """rule for calculating the mass of ice in each timestep"""
             for g in group:
                 for t in m.TIMESTEPS:
                     lhs = self.iceStatus[g,t]
-                    rhs = (self.tStor[g,t] == 0)
-                    block.ice_state.add((g, t), (lhs == rhs))
+                    rhs = self.tStor[g,t] - M*(1-self.tStor[g,t])
+                    block.ice_state_2.add((g, t), (lhs >= rhs))
 
-        self.ice_state = Constraint(group, m.TIMESTEPS, noruleinit=True)
-        self.ice_state_build = BuildAction(rule=_ice_state_rule)
+        self.ice_state_2 = Constraint(group, m.TIMESTEPS, noruleinit=True)
+        self.ice_state_2_build = BuildAction(rule=_ice_state_rule_2)
 
     def _objective_expression(self):
         """objective expression for storages with no investment"""
         m = self.parent_block()
-        if not hasattr(self, "ICESTORAGES"):
+        if not hasattr(self, "icestorages"):
             return 0
         costs = 0   # no cost from using the storage with a fixed capacity
         self.costs = Expression(expr=costs)