Merge pull request #20 from jylambert/implement-consumer-specific-flh

Implement consumer specific flh

examples/one-producer/data/config.yaml → examples/one-producer/data_mts/config.yaml

@@ -25,12 +25,10 @@ solver:
   time_limit: 10000
 
 economics:
-  source_flh: [2500]  # h/y
   source_c_inv: [0]  # Investment costs for each source
   source_c_irr: [0.08]  # Interest rate for each source
   source_lifetime: [40.]  # years for source investments
   source_price: [80e-3]  # Price for heat production at supply in €/kW
-  consumer_flh: 2500  # h/y
   heat_price: 120e-3  # Selling Price for heat in €/kW
   pipes_c_irr: 0.08  # Interest rate
   pipes_lifetime: 40.  # years for piping investments

examples/one-producer/data_sts/config.yaml

@@ -0,0 +1,34 @@
+water:
+  dynamic_viscosity: 4.041e-4
+  density: 977.76
+  heat_capacity_cp: 4.187e3
+
+ground:
+  thermal_conductivity: 2.4
+
+temperatures:
+  ambient: -20
+  supply: 70
+  return: 55
+
+piping:
+  diameter: [0.0216, 0.0285, 0.0372, 0.0431, 0.0545, 0.0703, 0.0825, 0.1071, 0.1325, 0.1603, 0.2101, 0.263, 0.3127, 0.3444, 0.3938]
+  outer_diameter: [0.09, 0.09, 0.11, 0.11, 0.125, 0.14, 0.16, 0.2, 0.225, 0.25, 0.315, 0.4, 0.45, 0.5, 0.56]
+  cost: [390, 400, 430, 464, 498, 537, 602, 670, 754, 886, 1171, 1184, 1197, 1401, 1755]
+  number_diameter: 15
+  max_pr_loss: 250
+  roughness: 0.05e-3
+  thermal_conductivity: 0.024
+
+solver:
+  mip_gap: 1e-4
+  time_limit: 10000
+
+economics:
+  source_c_inv: [0]  # Investment costs for each source
+  source_c_irr: [0.08]  # Interest rate for each source
+  source_lifetime: [40.]  # years for source investments
+  source_price: [80e-3]  # Price for heat production at supply in €/kW
+  heat_price: 120e-3  # Selling Price for heat in €/kW
+  pipes_c_irr: 0.08  # Interest rate
+  pipes_lifetime: 40.  # years for piping investments

examples/one-producer/data_sts/regression.csv

@@ -0,0 +1,2 @@
+id,heat_loss_a,heat_loss_b,capacity_a,capacity_b,power_flow_max_kW,power_flow_max_partload
+0,4.348e-07,0.02189,0.018377,567.335,6.9e+04,1.

examples/one-producer/run_mts_easy.py

@@ -14,7 +14,7 @@
 import topotherm as tt
 
 
-DATAPATH = os.path.join(os.path.dirname(os.path.abspath(__file__)), 'data')
+DATAPATH = os.path.join(os.path.dirname(os.path.abspath(__file__)), 'data_mts')
 OUTPUTPATH = os.path.join(os.path.dirname(os.path.abspath(__file__)), 'results/mts_eco_easy/')
 REGRESSION = 'regression.csv'  # regression coefficients for thermal capacity and heat losses
 TIMESERIES = 'timeseries.csv'  # timeseries for heat scaling
@@ -29,7 +29,7 @@ def read_regression(path, i):
     # read df and force floats
     df = pd.read_csv(path, sep=',', index_col=0, header=0, dtype=float)
     r_thermal_cap = {
-        "power_flow_max_kW" : df.loc[i, "power_flow_max_kW"],
+        "power_flow_max_kW": df.loc[i, "power_flow_max_kW"],
         "a": df.loc[i, "capacity_a"],
         "b": df.loc[i, "capacity_b"],
         "power_flow_max_partload": df.loc[i, "power_flow_max_partload"]
@@ -48,12 +48,7 @@ def main(filepath, outputpath, plots=True, solver='gurobi', mode='forced'):
 
     # Load the district
     mat = tt.fileio.load(filepath)
-    # read in demand profile
-    timeseries = pd.read_csv(os.path.join(filepath, TIMESERIES),
-                             sep=';', index_col=0, header=0).iloc[7:9, :].values.squeeze() #4:9
-    # create dummy profile, q_c should already contain the timeseries of all consumer demands
-    mat['q_c'] = mat['q_c'] * timeseries  # convert to timeseries
-
+
     if plots:
         f = tt.plotting.district(mat, isnot_init=False)  # Save initial District
         f.savefig(os.path.join(outputpath, 'district_initial.svg'), bbox_inches='tight')
@@ -65,8 +60,6 @@ def main(filepath, outputpath, plots=True, solver='gurobi', mode='forced'):
     settings = tt.settings.load(os.path.join(filepath, 'config.yaml'))
     # modify either in code or in the config file
     settings.economics.source_c_inv = [0.]  # no investment costs for sources
-    settings.economics.source_flh = [2500.]  # full load hours
-    settings.economics.consumers_flh = [2500.]  # full load hours
     settings.economics.pipes_lifetime = 40
     settings.economics.source_lifetime = [40]
     settings.temperatures.supply = 90
@@ -79,11 +72,7 @@ def main(filepath, outputpath, plots=True, solver='gurobi', mode='forced'):
         regression_inst=r_thermal_cap,
         regression_losses=r_heat_loss,
         economics=settings.economics,
-        optimization_mode=mode,
-        flh_scaling=timeseries.sum())
-    """model = tt.model_old.mts_easy_orig(
-        mat, model_sets, r_thermal_cap, r_heat_loss,
-        settings.economics, "eco", flh_scaling=1.9254)"""
+        optimization_mode=mode)
 
 
     # Optimization initialization

examples/one-producer/run_sts.py

@@ -15,7 +15,7 @@
 import topotherm as tt
 
 
-DATAPATH = os.path.join(os.path.dirname(os.path.abspath(__file__)), 'data')
+DATAPATH = os.path.join(os.path.dirname(os.path.abspath(__file__)), 'data_sts')
 OUTPUTPATH = os.path.join(os.path.dirname(os.path.abspath(__file__)),
                           'results', 'sts_forced')
 
@@ -64,11 +64,9 @@ def main(filepath, outputpath, plots=True, solver='gurobi', mode='economic'):
 
     # import settings
     settings = tt.settings.load(os.path.join(filepath, 'config.yaml'))
-    print(settings)
+
     # modify either in code or in the config file
     settings.economics.source_c_inv = [0.]  # no investment costs for sources
-    settings.economics.source_flh = [2500.]  # full load hours
-    settings.economics.consumers_flh = [2500.]  # full load hours
 
     model_sets = tt.sets.create(mat)
     model = tt.single_timestep.model(
@@ -84,7 +82,6 @@ def main(filepath, outputpath, plots=True, solver='gurobi', mode='economic'):
     opt.options['mipgap'] = settings.solver.mip_gap
     opt.options['timelimit'] = settings.solver.time_limit
     opt.options['logfile'] = os.path.join(outputpath, 'optimization.log')
-    #opt.options['Seed'] = 56324978
 
     result = opt.solve(model, tee=True)
 
@@ -115,5 +112,5 @@ def main(filepath, outputpath, plots=True, solver='gurobi', mode='economic'):
 
 if __name__ == '__main__':
     main(filepath=os.path.join(DATAPATH), outputpath=os.path.join(OUTPUTPATH),
-         plots=PLOTS, solver=SOLVER, mode='forced')
+         plots=PLOTS, solver=SOLVER, mode='economic')
     print(f'Finished {OUTPUTPATH}')

topotherm/fileio.py

@@ -44,10 +44,13 @@ def duplicate_columns(data, minoccur=2):
     q_c = (pd.read_parquet(
         os.path.join(path, 'Q_c.parquet')
         ).to_numpy().astype(float)) / 1000  # Data is in W, optimization in kW
+    flh_consumer = pd.read_parquet(os.path.join(path, 'flh_consumer.parquet')).to_numpy()
+    flh_source = pd.read_parquet(os.path.join(path, 'flh_source.parquet')).to_numpy()
     position = pd.read_parquet(
         os.path.join(path, 'rel_positions.parquet')
         ).iloc[:,  [-2,-1]].to_numpy().astype(float)
 
+    # @TODO Implement real warnings/errors and implement checks for full load hours
     if (a_i.sum(axis=0).sum() != 0) | (np.abs(a_i).sum(axis=0).sum()/2 != np.shape(a_i)[1]):
         print("Warning: The structure of A_i is not correct!")
     elif (-a_p.sum(axis=0).sum() != np.shape(a_p)[1]) | (np.abs(a_p).sum(axis=0).sum() != np.shape(a_p)[1]):
@@ -77,6 +80,8 @@ def duplicate_columns(data, minoccur=2):
     r['a_p'] = a_p
     r['a_c'] = a_c
     r['q_c'] = q_c
+    r['flh_consumer'] = flh_consumer
+    r['flh_source'] = flh_source
     r['l_i'] = length
     r['position'] = position
     return r

topotherm/multiple_timestep.py

@@ -30,8 +30,7 @@ def model(matrices: dict,
           regression_inst: dict,
           regression_losses: dict,
           economics: Economics,
-          optimization_mode: str,
-          flh_scaling: float): # @TODO flh_scaling will be removed when consumer specific flh are implemented
+          optimization_mode: str):
     """Create the optimization model for the thermo-hydraulic coupled with
     multiple time step operation.
 
@@ -260,9 +259,6 @@ def connection_to_consumer_fcd(m, d, j, t):
             mdl.cons, mdl.set_t,
             rule=connection_to_consumer_fcd,
             doc=msg_)
-    else:
-        raise NotImplementedError(
-            "Optimization mode %s not implemented" % optimization_mode)
 
     if optimization_mode == "forced":
         def connection_to_consumer_built_fcd(m, j):
@@ -284,7 +280,7 @@ def objective_function(m):
         fuel = sum(
             sum(m.P_source[k, t]
                 * economics.source_price[k]
-                * (economics.source_flh[k] / flh_scaling)
+                * matrices['flh_source'][k, t]
                 for k in m.set_n_p)
             for t in mdl.set_t)
 
@@ -302,13 +298,13 @@ def objective_function(m):
                                economics.source_lifetime[k])
                      for k in m.set_n_p)
 
-        # @TODO Implement consumer-specific flh in the economic mode
         if optimization_mode == "economic":
             revenue = sum(
                     sum(m.lambda_b[j]
+                        * matrices['flh_consumer'][k, t]
                         * matrices['q_c'][k, t]
                         for k, j in mdl.consumer_edges)
-                    for t in mdl.set_t) * (economics.consumers_flh[0] / flh_scaling) * economics.heat_price * (-1)
+                    for t in mdl.set_t) * economics.heat_price * (-1)
         else:
             revenue = 0
 

topotherm/settings.py

@@ -109,9 +109,9 @@ class Solver(BaseModel):
 class Economics(BaseModel):
     """Economic properties for the optimization problem. Used for the
     optimization model."""
-    source_flh: List[float] = Field(
-        default_factory=lambda: [1500.],
-        description="Full load hours of sources in h/y")
+    #source_flh: List[float] = Field(
+    #    default_factory=lambda: [1500.],
+    #    description="Full load hours of sources in h/y")
     source_price: List[float] = Field(
         default_factory=lambda: [80e-3],
         description="Price for heat production at supply in €/kW")
@@ -127,39 +127,39 @@ class Economics(BaseModel):
 
     pipes_c_irr: float = Field(0.08,
                                description="Interest rate for pipes")
-    consumers_flh: List[float] = Field(
-        default_factory=lambda: [1500.],
-        description="Full load hours of consumers in h/y")
+    #consumers_flh: List[float] = Field(
+    #    default_factory=lambda: [1500.],
+    #    description="Full load hours of consumers in h/y")
     heat_price: float = Field(120e-3,
                               description="Selling price for heat in €/kW")
     pipes_lifetime: float = Field(
         50.0,
         description="Lifetime for piping investments in years")
 
-    @model_validator(mode='after')
-    def check_len(self) -> Self:
-        """Check if the length of source_flh, source_price, source_c_inv,
-        source_c_irr and source_lifetime is consistent with the defined
-        number of sources.
-        """
-        if len(self.source_flh) != len(self.source_price):
-            raise ValueError(
-                f"""Length of source_flh {len(self.source_flh)} is not equal
-                to source_price {len(self.source_price)}""")
-        if len(self.source_flh) != len(self.source_c_inv):
-            raise ValueError(
-                f"""Length of source_flh {len(self.source_flh)} is not equal
-                to c_inv_source {len(self.source_flh)}""")
-        if len(self.source_flh) != len(self.source_lifetime):
-            raise ValueError(
-                f"""Length of source_flh {len(self.source_flh)} is not equal
-                to source_lifetime {len(self.source_lifetime)}""")
-        if len(self.source_flh) != len(self.source_c_irr):
-            raise ValueError(
-                f"""Length of source_flh {len(self.source_flh)} is not equal
-                to source_c_irr {len(self.source_c_irr)}""")
-
-        return self
+    # @model_validator(mode='after')
+    # def check_len(self) -> Self:
+    #    """Check if the length of source_flh, source_price, source_c_inv,
+    #    source_c_irr and source_lifetime is consistent with the defined
+    #    number of sources.
+    #    """
+    #    if len(self.source_flh) != len(self.source_price):
+    #        raise ValueError(
+    #            f"""Length of source_flh {len(self.source_flh)} is not equal
+    #            to source_price {len(self.source_price)}""")
+    #    if len(self.source_flh) != len(self.source_c_inv):
+    #        raise ValueError(
+    #            f"""Length of source_flh {len(self.source_flh)} is not equal
+    #            to c_inv_source {len(self.source_flh)}""")
+    #    if len(self.source_flh) != len(self.source_lifetime):
+    #        raise ValueError(
+    #            f"""Length of source_flh {len(self.source_flh)} is not equal
+    #            to source_lifetime {len(self.source_lifetime)}""")
+    #    if len(self.source_flh) != len(self.source_c_irr):
+    #        raise ValueError(
+    #            f"""Length of source_flh {len(self.source_flh)} is not equal
+    #            to source_c_irr {len(self.source_c_irr)}""")
+    #
+    #    return self
 
 
 class Settings(BaseModel):