Merge pull request #90 from SPF-OST/heat_pump_pvt_source_example

updates to PVT class and examples

data/examples/pvt_example.py

@@ -0,0 +1,71 @@
+import pandas as pd
+import os
+import pathlib as _pl
+
+from optihood.IO.writers import ScenarioFileWriterExcel
+from optihood.energy_network import EnergyNetworkGroup as EnergyNetwork
+
+if __name__ == '__main__':
+    # 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")
+
+    # define paths for input and result files
+    curDir = _pl.Path(__file__).resolve().parent
+    inputFilePath = curDir / ".." / "excels" / "pvt_example"
+    inputfileName = "scenario.xls"
+
+    resultFilePath = curDir / ".." / "results"
+    resultFileName = "results_pvt_example.xlsx"
+
+    print("Scenario file path: " + os.path.join(inputFilePath, inputfileName))
+    print("Result file path: " + os.path.join(resultFilePath, resultFileName))
+
+    # initialize parameters
+    numberOfBuildings = 4
+    optimizationType = "costs"  # set as "env" for environmental optimization and "costs" for cost optimization
+    mergeLinkBuses = True
+    mergeBuses = ["electricity", "heat_buses"]
+    dispatchMode = False  # Set to True to run the optimization in dispatch mode
+
+    # solver specific command line options
+    optimizationOptions = {
+        "gurobi": {
+            "BarConvTol": 0.5,
+            # The barrier solver terminates when the relative difference between the primal and dual objective values
+            # is less than the specified tolerance (with a GRB_OPTIMAL status)
+            "OptimalityTol": 1e-4,
+            # Reduced costs must all be smaller than OptimalityTol in the improving direction in order for a model to
+            # be declared optimal
+            "MIPGap": 1e-2,
+            # Relative Tolerance between the best integer objective and the objective of the best node remaining
+            "MIPFocus": 2
+            # 1 feasible solution quickly. 2 proving optimality. 3 if the best objective bound is moving very
+            # slowly/focus on the bound "Cutoff": #Indicates that you aren't interested in solutions whose objective
+            # values are worse than the specified value., could be dynamically be used in moo
+        }
+        # add the command line options of the solver here, For example to use CBC add a new item to the dictionary
+        # "cbc": {"tee": False}
+    }
+
+    # create an energy network and set the network parameters from an excel file
+    network = EnergyNetwork(timePeriod, temperatureLevels=True)
+    network.setFromExcel(os.path.join(inputFilePath, inputfileName), numberOfBuildings, opt=optimizationType,
+                         mergeLinkBuses=mergeLinkBuses, mergeBuses=mergeBuses, dispatchMode=dispatchMode)
+
+    # optimize the energy network
+    limit, capacitiesTransformers, capacitiesStorages = network.optimize(solver='gurobi',
+                                                                         numberOfBuildings=numberOfBuildings,
+                                                                         options=optimizationOptions,
+                                                                         mergeLinkBuses=mergeLinkBuses)
+
+    # print optimization outputs i.e. costs, environmental impact and capacities selected for different components (
+    # with investment optimization)
+    network.printInvestedCapacities(capacitiesTransformers, capacitiesStorages)
+    network.printCosts()
+    network.printEnvImpacts()
+    network.printMetaresults()
+
+    # save results
+    if not os.path.exists(resultFilePath):
+        os.makedirs(resultFilePath)
+    network.exportToExcel(os.path.join(resultFilePath, resultFileName), mergeLinkBuses=mergeLinkBuses)

data/examples/pvt_source_heat_pump.py

@@ -0,0 +1,71 @@
+import pandas as pd
+import os
+import pathlib as _pl
+
+from optihood.IO.writers import ScenarioFileWriterExcel
+from optihood.energy_network import EnergyNetworkGroup as EnergyNetwork
+
+if __name__ == '__main__':
+    # 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")
+
+    # define paths for input and result files
+    curDir = _pl.Path(__file__).resolve().parent
+    inputFilePath = curDir / ".." / "excels" / "pvt_source_heat_pump"
+    inputfileName = "scenario.xls"
+
+    resultFilePath = curDir / ".." / "results"
+    resultFileName = "results_pvt_source_heat_pump.xlsx"
+
+    print("Scenario file path: " + os.path.join(inputFilePath, inputfileName))
+    print("Result file path: " + os.path.join(resultFilePath, resultFileName))
+
+    # initialize parameters
+    numberOfBuildings = 4
+    optimizationType = "costs"  # set as "env" for environmental optimization and "costs" for cost optimization
+    mergeLinkBuses = True
+    mergeBuses = ["electricity", "heat_buses"]
+    dispatchMode = False  # Set to True to run the optimization in dispatch mode
+
+    # solver specific command line options
+    optimizationOptions = {
+        "gurobi": {
+            "BarConvTol": 0.5,
+            # The barrier solver terminates when the relative difference between the primal and dual objective values
+            # is less than the specified tolerance (with a GRB_OPTIMAL status)
+            "OptimalityTol": 1e-4,
+            # Reduced costs must all be smaller than OptimalityTol in the improving direction in order for a model to
+            # be declared optimal
+            "MIPGap": 1e-2,
+            # Relative Tolerance between the best integer objective and the objective of the best node remaining
+            "MIPFocus": 2
+            # 1 feasible solution quickly. 2 proving optimality. 3 if the best objective bound is moving very
+            # slowly/focus on the bound "Cutoff": #Indicates that you aren't interested in solutions whose objective
+            # values are worse than the specified value., could be dynamically be used in moo
+        }
+        # add the command line options of the solver here, For example to use CBC add a new item to the dictionary
+        # "cbc": {"tee": False}
+    }
+
+    # create an energy network and set the network parameters from an excel file
+    network = EnergyNetwork(timePeriod, temperatureLevels=True)
+    network.setFromExcel(os.path.join(inputFilePath, inputfileName), numberOfBuildings, opt=optimizationType,
+                         mergeLinkBuses=mergeLinkBuses, mergeBuses=mergeBuses, dispatchMode=dispatchMode)
+
+    # optimize the energy network
+    limit, capacitiesTransformers, capacitiesStorages = network.optimize(solver='gurobi',
+                                                                         numberOfBuildings=numberOfBuildings,
+                                                                         options=optimizationOptions,
+                                                                         mergeLinkBuses=mergeLinkBuses)
+
+    # print optimization outputs i.e. costs, environmental impact and capacities selected for different components (
+    # with investment optimization)
+    network.printInvestedCapacities(capacitiesTransformers, capacitiesStorages)
+    network.printCosts()
+    network.printEnvImpacts()
+    network.printMetaresults()
+
+    # save results
+    if not os.path.exists(resultFilePath):
+        os.makedirs(resultFilePath)
+    network.exportToExcel(os.path.join(resultFilePath, resultFileName), mergeLinkBuses=mergeLinkBuses)

optihood/buildings.py

@@ -83,15 +83,16 @@ def addBus(self, data, opt, mergeLinkBuses, mergeHeatSourceSink, electricityImpa
                         # add the excess production cost to self.__costParam
                         self.__costParam["excess"+label] = float(b["excess costs"])
                     if "shortage" in b:
-                        self.__nodesList.append(
-                            solph.Source(
-                                label="shortage"+label,
-                                outputs={
-                                    self.__busDict[label]: solph.Flow(
-                                        variable_costs=float(b["shortage costs"])*(opt == "costs")  # if opt = "env" variable costs should be zero
-                                    )}))
-                        # add the excess production cost to self.__costParam
-                        self.__costParam["shortage"+label] = float(b["shortage costs"])
+                        if b["shortage"]:
+                            self.__nodesList.append(
+                                solph.components.Source(
+                                    label="shortage"+label,
+                                    outputs={
+                                        self.__busDict[label]: solph.Flow(
+                                            variable_costs=float(b["shortage costs"])*(opt == "costs")  # if opt = "env" variable costs should be zero
+                                        )}))
+                            # add the excess production cost to self.__costParam
+                            self.__costParam["shortage"+label] = float(b["shortage costs"])
 
         if (mergeLinkBuses or mergeHeatSourceSink) and self.__buildingLabel=='Building1':
             return self.__busDict
@@ -232,12 +233,12 @@ def addPVT(self, data, data_timeseries, opt, mergeLinkBuses, dispatchMode):
                     inputBusLabel = s["from"]
                 else:
                     inputBusLabel = s["from"] + '__' + self.__buildingLabel
-
-                outputBuses = [self.__busDict[s["to"].split(",")[0] + '__' + self.__buildingLabel],
-                               self.__busDict[s["to"].split(",")[1] + '__' + self.__buildingLabel],
-                               self.__busDict[s["to"].split(",")[2] + '__' + self.__buildingLabel]]
-                connectBuses = [self.__busDict[s["connect"].split(",")[0] + '__' + self.__buildingLabel],
-                               self.__busDict[s["connect"].split(",")[1] + '__' + self.__buildingLabel]]
+                outputBuses = [self.__busDict[o + '__' + self.__buildingLabel] for o in s["to"].split(",")]
+                connectBuses = [self.__busDict[c + '__' + self.__buildingLabel] for c in s["connect"].split(",")]
+                if isinstance(s["delta_temp_n"], float) or isinstance(s["delta_temp_n"], int):
+                    delta_temp_n = [float(s["delta_temp_n"])]
+                else:
+                    delta_temp_n = [float(t) for t in s["delta_temp_n"].split(",")]
                 if opt == "costs":
                     epc = self._calculateInvest(s)[0]
                     base = self._calculateInvest(s)[1]
@@ -270,26 +271,38 @@ def addPVT(self, data, data_timeseries, opt, mergeLinkBuses, dispatchMode):
                                                        float(s["longitude"]), float(s["tilt"]), s["roof_area"],
                                                        float(s["zenith_angle"]), float(s["azimuth"]),
                                                        float(s["eta_0"]), float(s["a_1"]), float(s["a_2"]), float(s["temp_collector_inlet"]),data_timeseries['tre200h0'],
-                                                       float(s["delta_temp_n"]), data_timeseries['gls'], data_timeseries['str.diffus'], float(s["capacity_min"]), float(s["capacity_max"]),
-                                                       epc, base, env_capa, env_flow, varc, dispatchMode,
-                                                        float(s["pv_efficiency"]))
+                                                       delta_temp_n, data_timeseries['gls'], data_timeseries['str.diffus'], float(s["capacity_min"]), float(s["capacity_max"]),
+                                                       epc, base, env_capa, env_flow, varc, dispatchMode)
                 nodes = [pvtcollector.getPVT("el_source")]
                 for t in ["heat_source", "heat_transformer", "excess_heat_sink"]:
-                    nodes.append(pvtcollector.getPVT(t))
+                    sh, T2, dhw = pvtcollector.getPVT(t)
+                    if outputBuses.__len__() == 3:
+                        nodes.extend([sh, dhw])
+                    elif outputBuses.__len__() == 4:
+                        nodes.extend([sh, T2, dhw])
+                    else:
+                        nodes.extend([sh])
                 for x in nodes:
                     self.__nodesList.append(x)
 
-
-                self.__envParam['heatSourceSH_' + s['label'] + '__' + self.__buildingLabel] = envParam
-                self.__envParam['heatSourceDHW_' + s['label'] + '__' + self.__buildingLabel] = [env_flow, 0, 0]
-                self.__costParam['heatSourceSH_' + s['label'] + '__' + self.__buildingLabel] = [self._calculateInvest(s)[0],
-                                                                                            self._calculateInvest(s)[1]]
-                self.__costParam['heatSourceDHW_' + s['label'] + '__' + self.__buildingLabel] = [0, 0]
-
+                self.__envParam['heatSource_SH' + s['label'] + '__' + self.__buildingLabel] = envParam
+                self.__costParam['heatSource_SH' + s['label'] + '__' + self.__buildingLabel] = [self._calculateInvest(s)[0],
+                                                                                                self._calculateInvest(s)[1]]
                 self.__technologies.append(
                     [outputBuses[0], s["label"] + 'SH__' + self.__buildingLabel])
-                self.__technologies.append(
-                    [outputBuses[1], s["label"] + 'DHW__' + self.__buildingLabel])
+
+                if outputBuses.__len__() >= 3:
+                    self.__envParam['heatSource_DHW' + s['label'] + '__' + self.__buildingLabel] = [env_flow, 0, 0]
+                    self.__costParam['heatSource_DHW' + s['label'] + '__' + self.__buildingLabel] = [0, 0]
+                    self.__technologies.append(
+                        [outputBuses[-2], s["label"] + 'DHW__' + self.__buildingLabel])
+
+                if outputBuses.__len__() == 4:
+                    self.__envParam['heatSource_T2' + s['label'] + '__' + self.__buildingLabel] = [env_flow, 0, 0]
+                    self.__costParam['heatSource_T2' + s['label'] + '__' + self.__buildingLabel] = [0, 0]
+                    self.__technologies.append(
+                        [outputBuses[1], s["label"] + 'T2__' + self.__buildingLabel])
+
 
     def addGridSeparation(self, dataGridSeparation, mergeLinkBuses):
         if not dataGridSeparation.empty:
@@ -454,6 +467,10 @@ def _addHeatPump(self, data, operationTempertures, temperatureAmb, opt, mergeLin
         hpSHLabel = data["label"] + '__' + self.__buildingLabel
         if mergeLinkBuses and data["from"] in self.__linkBuses:
             inputBusLabel = [data["from"]]
+        elif mergeLinkBuses and any([b in self.__linkBuses for b in data["from"].split(',')]):
+            inputBusLabel = [i + '__' + self.__buildingLabel for i in data["from"].split(",") if
+                             i not in self.__linkBuses]
+            inputBusLabel.extend([i for i in data["from"].split(",") if i in self.__linkBuses])
         elif mergeHeatSourceSink and any([b in self.__heatSourceSinkBuses for b in data["from"].split(',')]):
             inputBusLabel = [i + '__' + self.__buildingLabel for i in data["from"].split(",") if
                              i not in self.__heatSourceSinkBuses]

optihood/combined_prod.py

@@ -89,8 +89,8 @@ def _second_input_relation_rule(block):
             for t in m.TIMESTEPS:
                 for g in group:
                     if len(list(g.inputs)) > 1:
-                        lhs = (len(list(g.inputs)) > 1) * m.flow[g.inflowQevap, g, t]
-                        rhs = (len(list(g.inputs)) > 1) * (m.flow[g, g.outputSH, t] + m.flow[g, g.outputDHW, t] - m.flow[g.inflow, g, t])
+                        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])
                         block.input_relation.add((g, t), (lhs == rhs))
 
         self.input_relation = Constraint(

optihood/constraints.py

@@ -273,28 +273,39 @@ def PVTElectricalThermalCapacityConstraint(om, numBuildings):
     pvtElOutFlows = [(i, o) for (i, o) in om.flows if ("elSource_pvt" in i.label)]
     pvtShOutFlows = [(i, o) for (i, o) in om.flows if ("heatSource_SHpvt" in i.label)]
     pvtDhwOutFlows = [(i, o) for (i, o) in om.flows if ("heatSource_DHWpvt" in i.label)]
+    pvtT2OutFlows = [(i, o) for (i, o) in om.flows if ("heatSource_T2pvt" in i.label)]
     for b in range(1, numBuildings + 1):
-        elCapacity = [om.InvestmentFlow.invest[i, o] for (i, o) in pvtElOutFlows if ((f'__Building{b}') in o.label)]
-        shCapacity = [om.InvestmentFlow.invest[i, o] for (i, o) in pvtShOutFlows if ((f'__Building{b}') in o.label)]
-        dhwCapacity = [om.InvestmentFlow.invest[i, o] for (i, o) in pvtDhwOutFlows if ((f'__Building{b}') in o.label)]
+        elCapacity = [om.InvestmentFlowBlock.invest[i, o] for (i, o) in pvtElOutFlows if ((f'__Building{b}') in o.label)]
+        shCapacity = [om.InvestmentFlowBlock.invest[i, o] for (i, o) in pvtShOutFlows if ((f'__Building{b}') in o.label)]
+        dhwCapacity = [om.InvestmentFlowBlock.invest[i, o] for (i, o) in pvtDhwOutFlows if ((f'__Building{b}') in o.label)]
+        T2Capacity = [om.InvestmentFlowBlock.invest[i, o] for (i, o) in pvtT2OutFlows if ((f'__Building{b}') in o.label)]
         areaUnitCapEl = [getattr(om.flows[i, o].investment, 'space_el') for (i, o) in pvtElOutFlows if ((f'__Building{b}') in o.label)]
         areaUnitCapSh = [getattr(om.flows[i, o].investment, 'space') for (i, o) in pvtShOutFlows if ((f'__Building{b}') in o.label)]
         if elCapacity or shCapacity:
             elCapacity = elCapacity[0]
             shCapacity = shCapacity[0]
-            dhwCapacity = dhwCapacity[0]
             areaUnitCapEl = areaUnitCapEl[0]
             areaUnitCapSh = areaUnitCapSh[0]
-            expr = (elCapacity*areaUnitCapEl == shCapacity*areaUnitCapSh)
+            expr1 = (elCapacity * areaUnitCapEl == shCapacity * areaUnitCapSh)
             setattr(
                 om,
-                "PVTSizeConstrElTh_B"+str(b),
-                pyo.Constraint(expr=expr),
-            )
-            expr = (dhwCapacity == shCapacity)
-            setattr(
-                om,
-                "PVTSizeConstrDhwSh_B" + str(b),
-                pyo.Constraint(expr=expr),
+                "PVTSizeConstrElTh_B" + str(b),
+                pyo.Constraint(expr=expr1),
             )
+            if dhwCapacity:
+                dhwCapacity = dhwCapacity[0]
+                expr2 = (dhwCapacity == shCapacity)
+                setattr(
+                    om,
+                    "PVTSizeConstrDhwSh_B" + str(b),
+                    pyo.Constraint(expr=expr2),
+                )
+            if T2Capacity:
+                T2Capacity = T2Capacity[0]
+                expr3 = (T2Capacity == shCapacity)
+                setattr(
+                    om,
+                    "PVTSizeConstrT2Sh_B" + str(b),
+                    pyo.Constraint(expr=expr3),
+                )
     return om