Merge pull request #19 from jylambert/rename_mts

Rename mts

examples/one-producer/run_mts_easy.py

@@ -12,15 +12,14 @@
 import pyomo.environ as pyo
 
 import topotherm as tt
-from topotherm.settings import Optimization
 
 
 DATAPATH = os.path.join(os.path.dirname(os.path.abspath(__file__)), 'data')
 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
 PLOTS = True  # save plots of the district
-SOLVER = 'gurobi' # 'gurobi' or 'cbc'
+SOLVER = 'gurobi'  # 'gurobi' or 'cbc'
 
 
 def read_regression(path, i):
@@ -41,7 +40,8 @@ def read_regression(path, i):
     }
     return r_thermal_cap, r_heat_loss
 
-def main(filepath, outputpath, plots=True, solver='gurobi'):
+
+def main(filepath, outputpath, plots=True, solver='gurobi', mode='forced'):
     """Main function to run the optimization"""
     # Create output directory if it does not exist
     tt.utils.create_dir(outputpath)
@@ -50,26 +50,46 @@ def main(filepath, outputpath, plots=True, solver='gurobi'):
     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
+                             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 = tt.plotting.district(mat, isnot_init=False)  # Save initial District
         f.savefig(os.path.join(outputpath, 'district_initial.svg'), bbox_inches='tight')
 
     # regression
     r_thermal_cap, r_heat_loss = read_regression(os.path.join(filepath, REGRESSION), 0)
 
-    model_sets = tt.model.create_sets(mat)
-    settings = Optimization()
-    model = tt.model.mts_easy(mat, model_sets, r_thermal_cap, r_heat_loss,
-                              economics=settings.economics, opt_mode="eco", flh_scaling=timeseries.sum())
+    # import settings
+    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
+    settings.economics.heat_price = 120e-3
+
+    model_sets = tt.sets.create(mat)
+    model = tt.multiple_timestep.model(
+        matrices=mat,
+        sets=model_sets,
+        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 initialization
     opt = pyo.SolverFactory(solver)
-    opt.options['mipgap'] = settings.opt_settings.mip_gap
-    opt.options['timelimit'] = settings.opt_settings.time_limit
+    opt.options['mipgap'] = settings.solver.mip_gap
+    opt.options['timelimit'] = settings.solver.time_limit
     opt.options['logfile'] = os.path.join(outputpath, 'optimization.log')
 
     result = opt.solve(model, tee=True)
@@ -82,10 +102,9 @@ def main(filepath, outputpath, plots=True, solver='gurobi'):
     dfsol = tt.utils.solver_to_df(result, model, solver=solver)
     dfsol.to_csv(os.path.join(outputpath, 'solver.csv'), sep=';')
 
-    opt_mats = tt.postprocessing.postprocess(model, mat, model_sets,
-                                             "mts",
-                                             t_return=settings.temperatures.return_,
-                                             t_supply=settings.temperatures.supply)
+    opt_mats = tt.postprocessing.mts(model=model,
+                                     matrices=mat,
+                                     settings=settings)
 
     # iterate over opt_mats and save each matrix as parquet file
     for key, value in opt_mats.items():
@@ -96,6 +115,8 @@ def main(filepath, outputpath, plots=True, solver='gurobi'):
         f = tt.plotting.district(opt_mats, diameter=opt_mats['d_i_0'], isnot_init=True)
         f.savefig(os.path.join(outputpath, 'district_optimal.svg'), bbox_inches='tight')
 
+
 if __name__ == '__main__':
-    main(filepath=DATAPATH, outputpath=OUTPUTPATH, plots=True)
+    main(filepath=os.path.join(DATAPATH), outputpath=os.path.join(OUTPUTPATH),
+         plots=PLOTS, solver=SOLVER, mode='economic')
     print(f'Finished {OUTPUTPATH}')

tests/mts_easy.py

@@ -4,8 +4,6 @@
 import pandas as pd
 
 import topotherm as tt
-from topotherm.settings import Optimization
-
 
 def read_regression(path, i):
     """Read the regression coefficients for the thermal capacity and heat
@@ -26,23 +24,42 @@ def read_regression(path, i):
     return r_thermal_cap, r_heat_loss
 
 
-def test_mtseasy_forced():
+def test_mtseasy_forced(request):
     """Main function to run the optimization"""
+    solver_name = request.config.getoption("--solver")
+    assert solver_name in ["scip", "gurobi", "cbc"], f"Unsupported solver: {solver_name}"
     # Load the district
     current_path = os.path.dirname(os.path.abspath(__file__))
     mat = tt.fileio.load(os.path.join(current_path, 'data'))
 
+    # read in demand profile
+    timeseries = pd.read_csv(os.path.join(current_path, 'data/timeseries.csv'),
+                             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
     # regression
     r_thermal_cap, r_heat_loss = read_regression(
         os.path.join(current_path, 'data', 'regression.csv'), 0)
 
-    model_sets = tt.model.create_sets(mat)
-    settings = Optimization()
-    model = tt.model.mts_easy(mat, model_sets, r_thermal_cap, r_heat_loss,
-                              economics=settings.economics, opt_mode="forced", flh_scaling=1.9254)
+    # import settings
+    settings = tt.settings.Settings()
+    settings.economics.source_flh = [2500.]  # full load hours
+    settings.economics.consumers_flh = [2500.]  # full load hours
+
+    model_sets = tt.sets.create(mat)
+    model = tt.multiple_timestep.model(
+        matrices=mat,
+        sets=model_sets,
+        regression_inst=r_thermal_cap,
+        regression_losses=r_heat_loss,
+        economics=settings.economics,
+        optimization_mode="forced",
+        flh_scaling=1.9254)
+
 
     # Optimization initialization
-    opt = pyo.SolverFactory('scip')
+    opt = pyo.SolverFactory(solver_name)
     opt.options['mipgap'] = 0.01
     opt.options['timelimit'] = 3600
 
@@ -53,23 +70,41 @@ def test_mtseasy_forced():
     assert (abs(pyo.value(model.obj)) - 4.6259e+06) < 0.02 * 4.6259e+06
 
 
-def test_mtseasy_eco():
+def test_mtseasy_eco(request):
     """Main function to run the optimization"""
+    solver_name = request.config.getoption("--solver")
+    assert solver_name in ["scip", "gurobi", "cbc"], f"Unsupported solver: {solver_name}"
     # Load the district
     current_path = os.path.dirname(os.path.abspath(__file__))
     mat = tt.fileio.load(os.path.join(current_path, 'data'))
 
+    # read in demand profile
+    timeseries = pd.read_csv(os.path.join(current_path, 'data/timeseries.csv'),
+                             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
     # regression
     r_thermal_cap, r_heat_loss = read_regression(
         os.path.join(current_path, 'data', 'regression.csv'), 0)
 
-    model_sets = tt.model.create_sets(mat)
-    settings = Optimization()
-    model = tt.model.mts_easy(mat, model_sets, r_thermal_cap, r_heat_loss, 
-                              economics=settings.economics, opt_mode="eco", flh_scaling=1.9254)
+    # import settings
+    settings = tt.settings.Settings()
+    settings.economics.source_flh = [2500.]  # full load hours
+    settings.economics.consumers_flh = [2500.]  # full load hours
+
+    model_sets = tt.sets.create(mat)
+    model = tt.multiple_timestep.model(
+        matrices=mat,
+        sets=model_sets,
+        regression_inst=r_thermal_cap,
+        regression_losses=r_heat_loss,
+        economics=tt.settings.Settings().economics,
+        optimization_mode="economic",
+        flh_scaling=1.9254)
 
     # Optimization initialization
-    opt = pyo.SolverFactory('scip')
+    opt = pyo.SolverFactory(solver_name)
     opt.options['mipgap'] = 0.01
     opt.options['timelimit'] = 3600
 

topotherm/__init__.py

@@ -5,4 +5,6 @@
 from . import postprocessing
 from . import settings
 from . import single_timestep
+from . import multiple_timestep
 from . import sets
+from . import model_old

topotherm/fileio.py

@@ -5,6 +5,7 @@
 import pandas as pd
 
 
+# @TODO: Insert into the load function consumer specific flh
 def load(path):
     """Read the input data from the given path and return the matrices A_i,
     A_p, A_c, Heat Demand, Length of edges, and positions.

topotherm/model_old.py

@@ -227,6 +227,7 @@ def one_pipe(m, j):
         return m.lambda_dir_1[j] + m.lambda_dir_2[j] <= 1
     model.one_pipe = pyo.Constraint(model.set_n_i, rule=one_pipe,
                                     doc='Just one Direction for each pipe')
+
     # @TODO: Develop total energy conservation equation for the eco mode (testing needed if beneficial)
     if opt_mode == "forced":
         def total_energy_cons(m, t):
@@ -272,7 +273,7 @@ def objective_function(m):
 # @TODO: discuss with jerry simplification strategies, since both models share a lot of equations.
 # @TODO: change the flh_scaling somehow
 # @TODO: implement existing pipes and sources
-def mts_easy(matrices, sets, regression_caps, regression_losses,
+def mts_easy_orig(matrices, sets, regression_caps, regression_losses,
              economics: Economics, opt_mode: str, flh_scaling: float):
     """Create the optimization model for the thermo-hydraulic coupled with multiple time 
     step operation. The model is based on the STS model and implements a simplified themal
@@ -294,7 +295,7 @@ def mts_easy(matrices, sets, regression_caps, regression_losses,
 
     p_max_pipe_const = float(regression_caps['power_flow_max_kW'])  # Big-M-Constraint for pipes
     p_max_source = matrices['q_c'].sum()*2  # Big-M-Constraint for source
-    model.flh = economics.flh / flh_scaling
+    model.flh = economics.consumers_flh[0] / flh_scaling
 
     # Define index sets
     model.set_n_i = pyo.Set(initialize=range(sets['a_i_shape'][1]),
@@ -479,18 +480,18 @@ def built_usage_mapping_help2(m, j, t):
 
     def objective_function(m):
         term1 = sum(
-            sum(m.P_source[k, t] * economics.source_price * model.flh for k in m.set_n_p)
+            sum(m.P_source[k, t] * economics.source_price[k] * model.flh for k in m.set_n_p)
             for t in model.set_t
                 )
         term2 = sum(
             (
                 (m.P_cap[k] * regression_caps['a']
                     + regression_caps['b'] * m.lambda_built[k]
-                 ) * annuity(economics.c_irr, economics.life_time) * matrices['l_i'][k]
+                 ) * annuity(economics.pipes_c_irr, economics.pipes_lifetime) * matrices['l_i'][k]
             ) for k in m.set_n_i
         )
-        term3 = sum(m.P_source_cap[k] * economics.c_inv_source[k]
-                    * annuity(economics.c_irr, economics.life_time) for k in m.set_n_p)
+        term3 = sum(m.P_source_cap[k] * economics.source_c_inv[k]
+                    * annuity(economics.source_c_irr[k], economics.source_lifetime[k]) for k in m.set_n_p)
 
         if opt_mode == "eco":
             term4 = sum(sum(