diff --git a/examples/one-producer/run_mts_easy.py b/examples/one-producer/run_mts_easy.py index 396b9f5..1d5997b 100644 --- a/examples/one-producer/run_mts_easy.py +++ b/examples/one-producer/run_mts_easy.py @@ -12,7 +12,6 @@ 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') @@ -20,7 +19,7 @@ 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}') diff --git a/tests/mts_easy.py b/tests/mts_easy.py index f61d11c..7a12db3 100644 --- a/tests/mts_easy.py +++ b/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 diff --git a/topotherm/__init__.py b/topotherm/__init__.py index acc1fbe..08ef63d 100644 --- a/topotherm/__init__.py +++ b/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 diff --git a/topotherm/fileio.py b/topotherm/fileio.py index fe5b1ae..4903cf8 100644 --- a/topotherm/fileio.py +++ b/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. diff --git a/topotherm/model_old.py b/topotherm/model_old.py index 4ebea39..523bbd9 100644 --- a/topotherm/model_old.py +++ b/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( diff --git a/topotherm/multiple_timestep.py b/topotherm/multiple_timestep.py index e69de29..927826d 100644 --- a/topotherm/multiple_timestep.py +++ b/topotherm/multiple_timestep.py @@ -0,0 +1,319 @@ +"""This module contains the optimization models for the multiple-timestep +district heating network design. + +The module contains the following functions: + * annuity: Calculate the annuity factor + * model: Create the optimization model for the multiple time steps operation +""" +import numpy as np +import pyomo.environ as pyo + +from topotherm.settings import Economics + + +def annuity(c_i, n): + """Calculate the annuity factor. + + Args: + c_i (float): Interest rate + n (float): Number of years + + Returns: + float: annuity factor + """ + a = ((1 + c_i) ** n * c_i) / ((1 + c_i) ** n - 1) + return a + + +def model(matrices: dict, + sets: 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 + """Create the optimization model for the thermo-hydraulic coupled with + multiple time step operation. + + Args: + matrices (dict): Dictionary with the matrices of the district heating + network with keys a_i, a_p, a_c, l_i, position, q_c + sets (dict): Dictionary with the sets for the optimization, obtained + from matrices + regression_inst (dict): Dictionary with the regression coefficients + for the thermal capacity + regression_losses (dict): Dictionary with the regression coefficients + for the heat losses + economics (topotherm.settings.Economics): Object with the economic + parameters + optimization_mode (str): Optimization mode, either 'economic' for + economic or 'forced' for forced operation + + Returns: + mdl (pyomo.environ.ConcreteModel): pyomo model + """ + + # Check if the optimization mode is implemented + if optimization_mode not in ['economic', 'forced']: + raise NotImplementedError( + "Optimization mode %s not implemented" % optimization_mode) + + # Initialize model + mdl = pyo.ConcreteModel() + + # Big-M-Constraint for pipes + p_max_pipe_const = float(regression_inst['power_flow_max_kW'].max()) + # Big-M-Constraint for source + p_max_source = matrices['q_c'].sum() * 2 + + # Define index sets + mdl.set_n_i = pyo.Set(initialize=range(sets['a_i_shape'][1]), + doc='Number of pipe connections supply/return line') + mdl.set_n_p = pyo.Set(initialize=range(sets['a_p_shape'][1]), + doc='Number of producers') + mdl.set_n_c = pyo.Set(initialize=range(sets['a_c_shape'][1]), + doc='Number of consumers') + mdl.set_n = pyo.Set(initialize=range(sets['a_i_shape'][0]), + doc='Nodes in supply/return line') + mdl.set_t = pyo.Set(initialize=range(matrices['q_c'].shape[1]), + doc='Time steps') + mdl.dirs = pyo.Set(initialize=['ij', 'ji'], + doc='Set of pipe directions.') + mdl.flow = pyo.Set(initialize=['in', 'out'], + doc='Flow direction in the pipe') + + # Define the combined set for pipes with consumers in both directions + mdl.cons = pyo.Set( + initialize=[('ij', edge) for edge in sets['connection_c_ij']] + + [('ji', edge) for edge in sets['connection_c_ji']], + dimen=2, + doc='Pipes with consumer in both directions') + + mdl.forced_edges = pyo.Set( + initialize=np.concatenate([sets['connection_c_ij'], sets['connection_c_ji']]), + doc='Pipes with a consumer in direction ij or ji' + ) + + mdl.consumer_edges = pyo.Set( + initialize=[(i, np.where((matrices['a_i'][np.where(matrices['a_c'][:, i] == 1)[0].item(), :] == 1) | + (matrices['a_i'][np.where(matrices['a_c'][:, i] == 1)[0].item(), :] == -1) + )[0].item()) for i in range(sets['a_c_shape'][1])], + dimen=2, + doc='Assign to each consumer the corresponding pipe' + ) + + # Define variables + pipe_power = {'bounds': (0, p_max_pipe_const), + 'domain': pyo.NonNegativeReals, + 'initialize': p_max_pipe_const} + mdl.P = pyo.Var( + mdl.dirs, mdl.flow, mdl.set_n_i, mdl.set_t, + doc='Heat power at the pipes', + **pipe_power) + + mdl.P_cap = pyo.Var( + mdl.set_n_i, + doc='Thermal capacity of each pipe', + **pipe_power + ) + + # Binaries for the chosen direction of a pipe + mdl.lambda_ = pyo.Var( + mdl.dirs, mdl.set_n_i, mdl.set_t, + initialize=1, + domain=pyo.Binary, + doc='Binary direction decisions') + + # Binary building decision of a pipe + mdl.lambda_b = pyo.Var( + mdl.set_n_i, + initialize=1, + domain=pyo.Binary, + doc='Binary building decision' + ) + + # Thermal power of the source + source_power = {'bounds': (0, p_max_source), + 'domain': pyo.PositiveReals, + 'initialize': p_max_source} + mdl.P_source = pyo.Var( + mdl.set_n_p, mdl.set_t, + doc='Thermal power of the source', + **source_power) + mdl.P_source_inst = pyo.Var( + mdl.set_n_p, + doc='Thermal capacity of the heat source', + **source_power) + + def heat_source_inst(m, j, t): + """Never exceed the installed capacity of the heat source.""" + return m.P_source[j, t] <= m.P_source_inst[j] + + mdl.cons_heat_source_inst = pyo.Constraint( + mdl.set_n_p, mdl.set_t, + rule=heat_source_inst, + doc='Upper bound for the heat source supply delivery') + + def nodal_power_balance(m, j, t): + """REFERENCE DIRECTION: left to right + P_ji, in P_ji, out + PIPE <------- NODE <------- PIPE + -------> -------> + P_ij, out P_ij, in + + Energy balance system: out - in = 0 + """ + pipe_to_node = sum(m.P['ji', 'in', k, t] + - m.P['ij', 'out', k, t] + for k in sets['a_i_in'][j]) + node_to_pipe = sum(m.P['ij', 'in', k, t] + - m.P['ji', 'out', k, t] + for k in sets['a_i_out'][j]) + sources = sum(- m.P_source[k, t] + for k in sets['a_p_in'][j]) + sink = 0 + if optimization_mode == "forced": + sink = sum(matrices['q_c'][k, t] + for k in sets['a_c_out'][j]) + elif optimization_mode == "economic": + sink = sum(m.lambda_b[k] for k in sets['a_c_out_edge'][j]) \ + * sum(matrices['q_c'][k, t] for k in sets['a_c_out'][j]) + return node_to_pipe + pipe_to_node + sources + sink == 0 + + mdl.cons_nodal_balance = pyo.Constraint( + mdl.set_n, mdl.set_t, + rule=nodal_power_balance, + doc='Nodal power balance for each time step') + + def power_balance_pipe(m, d, j, t): + """Power balance for the pipes. + + P_ji, out P_ji, in + <------- PIPE <------- + -------> -------> + P_ij, in P_ij, out + + """ + # flows into and out of pipe + flows = m.P[d, 'in', j, t] - m.P[d, 'out', j, t] + # thermal losses calculation + variable = (regression_losses['a'] + * regression_inst['power_flow_max_partload'] + * m.P[d, 'in', j, t]) + fix = regression_losses['b'] * m.lambda_[d, j, t] + return flows - (variable + fix) * matrices['l_i'][j] == 0 + + mdl.cons_power_balance_pipe = pyo.Constraint( + mdl.dirs, mdl.set_n_i, mdl.set_t, + rule=power_balance_pipe, + doc='Power balance for each pipe and time step') + + def power_bigm_P(m, d, j, t): + lhs = m.P[d, 'in', j, t] - p_max_pipe_const * m.lambda_[d, j, t] + rhs = 0 + return lhs <= rhs + mdl.cons_power_bigm_P = pyo.Constraint( + mdl.dirs, mdl.set_n_i, mdl.set_t, + rule=power_bigm_P, doc='Big-M constraint for power flow') + + def power_max_p_built_const(m, j): + return m.P_cap[j] - p_max_pipe_const * m.lambda_b[j] <= 0 + mdl.cons_power_max_P_built_const = pyo.Constraint( + mdl.set_n_i, rule=power_max_p_built_const, doc='Maximum Powerflow constant i->j / j->i') + + def power_max_p_built(m, d, j, t): + return m.P[d, 'in', j, t] - m.P_cap[j] <= 0 + mdl.cons_power_max_P_built = pyo.Constraint( + mdl.dirs, mdl.set_n_i, mdl.set_t, rule=power_max_p_built, + doc='Maximum Powerflow according to capacity i->j / j->i') + + def built_usage_mapping(m, d, j, t): + return m.lambda_[d, j, t] - m.lambda_b[j] <= 0 + mdl.cons_built_usage_mapping_help1 = pyo.Constraint(mdl.dirs, mdl.set_n_i, mdl.set_t, + rule=built_usage_mapping, + doc='Map lambda direction according to lambda_built') + + def one_pipe(m, j, t): + return m.lambda_['ij', j, t] + m.lambda_['ji', j, t] <= 1 + + mdl.one_pipe = pyo.Constraint(mdl.set_n_i, mdl.set_t, + rule=one_pipe, + doc='Just one Direction for each pipe') + + def connection_to_consumer_eco(m, d, j, t): + return m.lambda_[d, j, t] <= sets[f'lambda_c_{d}'][j] + + def connection_to_consumer_fcd(m, d, j, t): + return m.lambda_[d, j, t] == sets[f'lambda_c_{d}'][j] + + if optimization_mode == "economic": + msg_ = """Constraint if houses have their own connection-pipe + and set the direction (ij or ji)""" + mdl.cons_connection_to_consumer = pyo.Constraint( + mdl.cons, mdl.set_t, + rule=connection_to_consumer_eco, + doc=msg_) + elif optimization_mode == "forced": + msg_ = """Constraint if houses have their own connection-pipe + and set the direction (ij or ji)""" + mdl.cons_connection_to_consumer = pyo.Constraint( + 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): + return m.lambda_b[j] >= 1 + mdl.cons_connection_forced = pyo.Constraint( + mdl.forced_edges, rule=connection_to_consumer_built_fcd, + doc='The house connection has to be built to satisfy the demand') + + if optimization_mode == "forced": + def total_energy_conservation(m, t): + return sum(m.P_source[k, t] for k in m.set_n_p) \ + - sum(m.P['ij', 'in', k, t] - m.P['ij', 'out', k, t] for k in m.set_n_i) \ + - sum(m.P['ji', 'in', k, t] - m.P['ji', 'out', k, t] for k in m.set_n_i) \ + - sum(matrices['q_c'][k, t] for k in m.set_n_c) == 0 + mdl.cons_total_energy_cons = pyo.Constraint(mdl.set_t, rule=total_energy_conservation, + doc='Total energy conservation') + + def objective_function(m): + fuel = sum( + sum(m.P_source[k, t] + * economics.source_price[k] + * (economics.source_flh[k] / flh_scaling) + for k in m.set_n_p) + for t in mdl.set_t) + + pipes = sum( + ( + (m.P_cap[k] * regression_inst['a'] + + regression_inst['b'] * m.lambda_b[k] + ) * annuity(economics.pipes_c_irr, economics.pipes_lifetime) * matrices['l_i'][k] + ) for k in m.set_n_i + ) + + source = sum(m.P_source_inst[k] + * economics.source_c_inv[k] + * annuity(economics.source_c_irr[k], + 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['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) + else: + revenue = 0 + + return fuel + pipes + source + revenue + + mdl.obj = pyo.Objective(rule=objective_function, doc='Objective function') + + return mdl diff --git a/topotherm/postprocessing.py b/topotherm/postprocessing.py index 2da5bb0..b29ba63 100644 --- a/topotherm/postprocessing.py +++ b/topotherm/postprocessing.py @@ -107,8 +107,8 @@ def sts(model: pyo.ConcreteModel, # values lambda_ij for ji. This is necessary for the postprocessing. elif lambda_ji[q] == 1: matrices['a_i'][:, q] = matrices['a_i'][:, q] * (-1) - lambda_ij[q] = 1 - lambda_ji[q] = 0 + lambda_ij[q] = 1 # @TODO: For sts this can be removed + lambda_ji[q] = 0 # @TODO: For sts this can be removed p_lin = p_ij + p_ji # Power of the pipes @@ -153,8 +153,8 @@ def sts(model: pyo.ConcreteModel, a_c=a_c_opt, q_c=q_c_opt, l_i=l_i_opt, - lambda_ij_opt=lambda_ij, - lambda_ji_opt=lambda_ji, + lambda_ij_opt=lambda_ij, #@TODO: Remove lambda_ij_opt and lambda_ji_opt? -> lambda_ij_opt should only + lambda_ji_opt=lambda_ji, #@TODO: contain 1 and lamdba_ji_opt only 0 (for sts!). d_i_0=d_lin, m_i_0=m_lin, position=pos_opt, @@ -164,6 +164,129 @@ def sts(model: pyo.ConcreteModel, return res +def mts(model: pyo.ConcreteModel, + matrices: dict, + settings: Settings): + """Postprocessing for the multiple time step model. This includes the + calculation of the diameter and velocity of the pipes, the elimination of + unused pipes and nodes. + + Args: + model (pyo.ConcreteModel): solved pyomo model + matrices (dict): dict containing the matrices + settings (tt.settings.Settings): settings for the optimization + + Returns: + _dict: containing the variables and postprocessed data + """ + # Get the values from the model + p_ij = np.reshape(np.array(pyo.value(model.P['ij', 'in', :, :])), (-1, matrices['q_c'].shape[1])) + p_ji = np.reshape(np.array(pyo.value(model.P['ji', 'in', :, :])), (-1, matrices['q_c'].shape[1])) + p_cap = np.array(pyo.value(model.P_cap[:])) + + # flow direction, binary + lambda_ij = np.reshape(np.around(np.array(pyo.value(model.lambda_['ij', :, :])), 0), (-1, matrices['q_c'].shape[1])) + lambda_ji = np.reshape(np.around(np.array(pyo.value(model.lambda_['ji', :, :])), 0), (-1, matrices['q_c'].shape[1])) + lambda_b = np.around(np.array(pyo.value(model.lambda_b[:])), 0) + + q_c_opt = np.zeros([matrices['a_c'].shape[1], len(model.set_t)]) + + # Exclude non-connected consumers in Q_c, only affects the economic case + # Check for consumers connected in direction ij + for d, e in model.cons: + if d == 'ij': + # edge in incidence matrix where pipe exits into node n (==-1) + a_i_idx = np.where(matrices['a_i'][:, e] == -1) + # location where a_i_idx is connected to a_c + a_c_idx = np.where(matrices['a_c'][a_i_idx[0], :][0] == 1) + if len(a_i_idx) != 1 or len(a_c_idx) != 1: + raise ValueError('Error in the incidence matrix!') + # assign the heat demand to the connected consumer if lambda is 1 + q_c_opt[a_c_idx[0], :] = lambda_b[e] * matrices['q_c'][a_c_idx[0], :] + elif d == 'ji': + a_i_idx = np.where(matrices['a_i'][:, e] == 1) + a_c_idx = np.where(matrices['a_c'][a_i_idx[0], :][0] == 1) + if len(a_i_idx) != 1 or len(a_c_idx) != 1: + raise ValueError('Error in the incidence matrix!') + q_c_opt[a_c_idx[0], :] = lambda_b[e] * matrices['q_c'][a_c_idx[0], :] + + # Remove nonzero elements row-wise + q_c_opt = q_c_opt[q_c_opt.any(axis=1)] + + # Adaption of Incidence Matrix for further postprocessing + for q in model.set_n_i: + # if not active, all is 0 + if lambda_b[q] == 0: + matrices['a_i'][:, q] = 0 + matrices['l_i'][q] = 0 + # if opposite direction operational, switch a_i with -1 and switch + # values lambda_ij for ji. This is necessary for the postprocessing. + elif (lambda_b[q] == 1) & (lambda_ji[q, 0] == 1): + matrices['a_i'][:, q] = matrices['a_i'][:, q] * (-1) + lambda_ij[q, np.where(lambda_ij[q, 1:] == 0)[0]] = 1 + lambda_ji[q, np.where(lambda_ji[q, 1:] == 1)[0]] = 0 + + + p_lin = p_cap # Capacity of the pipes + + # drop entries with 0 in the incidence matrix to reduce size + valid_columns = matrices['a_i'].any(axis=0) + valid_rows = matrices['a_i'].any(axis=1) + + p_lin_opt = p_lin[valid_columns] + p_ij_opt = p_ij[valid_columns, :] + p_ji_opt = p_ji[valid_columns, :] + lambda_ij_opt = lambda_ij[valid_columns, :] + lambda_ji_opt = lambda_ji[valid_columns, :] + pos_opt = matrices['position'][valid_rows, :] + a_c_opt = matrices['a_c'][valid_rows, :] + a_p_opt = matrices['a_p'][valid_rows, :] + a_i_opt = matrices['a_i'][valid_rows, :][:, valid_columns] + l_i_opt = matrices['l_i'][valid_columns] + + a_i_shape_opt = np.shape(a_i_opt) # (rows 0, columns 1) + d_lin = np.zeros(a_i_shape_opt[1]) # Initialize linear diameters + v_lin = np.zeros(a_i_shape_opt[1]) # Initialize velocities + # Assign supply and return temperatures + supply_temp_opt = np.ones(a_i_shape_opt[1]) * settings.temperatures.supply + return_temp_opt = np.ones(a_i_shape_opt[1]) * settings.temperatures.return_ + + # Calculate the mass flow for each pipe with m cp deltaT = P + m_lin = (p_lin_opt * 1000 + / (settings.water.heat_capacity_cp + * (supply_temp_opt - return_temp_opt) + )) + + # Calculate the diameter and velocity for each pipe + for h in range(a_i_shape_opt[1]): + mass_lin = m_lin[h] + sol = root(lambda v: calc_diam_and_velocity(v, mass_lin, settings), + (0.5, 0.02), + method='lm') + if sol.success: + v_lin[h], d_lin[h] = sol.x + else: + print(h, 'failed to calculate diameter and velocity!') + + res = dict( + a_i=a_i_opt, + a_p=a_p_opt, + a_c=a_c_opt, + q_c=q_c_opt, + l_i=l_i_opt, + lambda_ij_opt=lambda_ij_opt, + lambda_ji_opt=lambda_ji_opt, + d_i_0=d_lin, + m_i_0=m_lin, + position=pos_opt, + p=p_lin_opt, + p_ij=p_ij_opt, + p_ji=p_ji_opt + ) + + return res + + def to_networkx_graph(matrices): """Export the postprocessed, optimal district as a networkx graph. @@ -214,98 +337,3 @@ def to_networkx_graph(matrices): # drop all edges with p=0 G.remove_edges_from([(u, v) for u, v, d in G.edges(data=True) if d['p'] == 0]) return G - - -def mts(model, matrices, sets, t_supply, t_return): - """returns all matrices and results for further processing. Essentially, - it simplifies the results from the optimization, including pipe diameter - and mass flow, eliminating the unused pipes and nodes. - - Args: - model (pyomo.environ.ConcreteModel): solved pyomo model - matrices (dict): dict containing the incidence matrices - sets (dict): dict containing the sets - temperatures (dict): dict containing the supply and return temperatures - - Returns: - dict: dict containing the updated matrices, including diameter and mass flow - """ - data_dict = {} - - p_cap = np.zeros([sets['a_i_shape'][1]]) - lambda_built = np.zeros([sets['a_i_shape'][1]]) - lambda_dir = np.zeros([sets['a_i_shape'][1], len(model.set_t)]) - p_source_built = np.zeros(sets['a_p_shape'][1]) - - for v in model.component_objects(pyo.Var, active=True): - var_dict = {(v.name, index): pyo.value(v[index]) for index in v} - data_dict.update(var_dict) - if v.name == "lambda_built": - for index in v: - lambda_built[index] = pyo.value(v[index]) - if v.name == "lambda_dir_1": - for index in v: - lambda_dir[index] = pyo.value(v[index]) - if v.name == "P_cap": - for index in v: - p_cap[index] = pyo.value(v[index]) - if v.name == "P_source_cap": - for index in v: - p_source_built[index] = pyo.value(v[index]) - - lambda_built = np.around(lambda_built, 0) - lambda_dir = np.around(lambda_dir, 0) - - # @TODO for all timesteps (with new formulation) - for q, _ in enumerate(lambda_dir_1): - if lambda_built[q] == 0: - matrices['a_i'][:, q] = 0 - matrices['l_i'][q] = 0 - elif (lambda_built[q] == 1) & (lambda_dir_1[q, 0] == 0): - matrices['a_i'][:, q] = matrices['a_i'][:, q] * (-1) - lambda_dir_1[q, np.where(lambda_dir_1[q, 1:] == 0)[0]] = 1 - lambda_dir_2[q, np.where(lambda_dir_2[q, 1:] == 1)[0]] = 0 - - p_cap_opt = np.delete(p_cap, np.where(~matrices['a_i'].any(axis=0))) - pos_opt = np.delete(matrices['position'], np.where(~matrices['a_i'].any(axis=1)), axis=0) - a_c_opt = np.delete(matrices['a_c'], np.where(~matrices['a_i'].any(axis=1)), axis=0) - a_p_opt = np.delete(matrices['a_p'], np.where(~matrices['a_i'].any(axis=1)), axis=0) - a_i_opt = matrices['a_i'] - a_i_opt = np.delete(a_i_opt, np.where(~a_i_opt.any(axis=0)), axis=1) - a_i_opt = np.delete(a_i_opt, np.where(~a_i_opt.any(axis=1)), axis=0) - l_i_opt = matrices['l_i'][matrices['l_i'] != 0] - - a_i_shape_opt = np.shape(a_i_opt) # (rows 0, columns 1) - d_lin2 = np.zeros(a_i_shape_opt[1]) - v_lin2 = np.zeros(a_i_shape_opt[1]) - supply_temp_opt = np.ones(a_i_shape_opt[1]) * t_supply - return_temp_opt = np.ones(a_i_shape_opt[1]) * t_return - - def equations(v): - vel, d = v - reynolds = (settings.Water.density * vel * d) / settings.Water.dynamic_viscosity - f = (-1.8 * np.log10((settings.Piping.roughness / (3.7 * d)) ** 1.11 + 6.9 / reynolds))**-2 - eq1 = vel - np.sqrt((2 * settings.Piping.max_pr_loss * d) / (f * settings.Water.density)) - eq2 = mass_lin - settings.Water.density * vel * (np.pi / 4) * d ** 2 - return [eq1, eq2] - - m_lin = (p_cap_opt*1000)/(settings.Water.heat_capacity_cp * (supply_temp_opt - return_temp_opt)) - - for h in range(a_i_shape_opt[1]): - mass_lin = m_lin[h] - v_lin2[h], d_lin2[h] = fsolve(equations, (0.5, 0.02)) - - res = dict( - a_i=a_i_opt, - a_p=a_p_opt, - a_c=a_c_opt, - q_c=matrices['q_c'], - l_i=l_i_opt, - d_i_0=d_lin2, - m_i_0=m_lin, - position=pos_opt, - - ) - - return res - diff --git a/topotherm/sets.py b/topotherm/sets.py index 4575fde..95275e8 100644 --- a/topotherm/sets.py +++ b/topotherm/sets.py @@ -3,6 +3,7 @@ """ import numpy as np + def create(matrices): """Create sets for the optimization. The sets are used to define the variables and constraints. @@ -25,13 +26,13 @@ def create(matrices): s['connection_c_ij'] = np.where( matrices['a_i'][consumers, :].sum(axis=0) == -1 - )[0] + )[0] s['lambda_c_ij'] = np.zeros(s['a_i_shape'][1]) s['lambda_c_ij'][s['connection_c_ij']] = 1 s['connection_c_ji'] = np.where( matrices['a_i'][consumers, :].sum(axis=0) == 1 - )[0] + )[0] s['lambda_c_ji'] = np.zeros(s['a_i_shape'][1]) s['lambda_c_ji'][s['connection_c_ji']] = 1 @@ -47,6 +48,12 @@ def create(matrices): s['a_p_in'][i] = np.where(matrices['a_p'][i, :] == -1)[0] s['a_c_out'] = {} + s['a_c_out_edge'] = {} for i in range(s['a_c_shape'][0]): s['a_c_out'][i] = np.where(matrices['a_c'][i, :] == 1)[0] + if matrices['a_c'][i, :].any() == 1: + s['a_c_out_edge'][i] = np.where((matrices['a_i'][i, :] == 1) | (matrices['a_i'][i, :] == -1))[0] + else: + s['a_c_out_edge'][i] = [] + return s diff --git a/topotherm/settings.py b/topotherm/settings.py index 2c04830..d89cb77 100644 --- a/topotherm/settings.py +++ b/topotherm/settings.py @@ -105,7 +105,7 @@ class Solver(BaseModel): description="Time limit for the optimization in seconds") log: str = Field("solver.log", description="Log file for the solver") - +# @TODO: remove flh from setting and incorporate into fileio when modelling consumer specific flh class Economics(BaseModel): """Economic properties for the optimization problem. Used for the optimization model.""" diff --git a/topotherm/single_timestep.py b/topotherm/single_timestep.py index 2896d93..58d3896 100644 --- a/topotherm/single_timestep.py +++ b/topotherm/single_timestep.py @@ -51,27 +51,23 @@ def model(matrices: dict, Returns: mdl (pyomo.environ.ConcreteModel): pyomo model """ - # @TODO: Look with Amedeo if q_c can be adapted to dimensionless vector, - # (in theory it is possible to do - # @TODO: a unidirectional flow formulation with multiple time step with - # topotherm sts) - # check if the optimization mode is implemented + # Check if the optimization mode is implemented if optimization_mode not in ['economic', 'forced']: raise NotImplementedError( "Optimization mode %s not implemented" % optimization_mode) - # init model + # Initialize model mdl = pyo.ConcreteModel() # Big-M-Constraint for pipes - p_max_pipe_const = float(regression_inst['power_flow_max_kW'].max())*20 + p_max_pipe_const = float(regression_inst['power_flow_max_kW'].max()) # Big-M-Constraint for source - p_max_source = matrices['q_c'].sum()*20 + p_max_source = matrices['q_c'].sum()*2 # Define index sets mdl.set_n_i = pyo.Set(initialize=range(sets['a_i_shape'][1]), - doc='N° of Pipe connections supply/return line') + doc='Number of pipe connections supply/return line') mdl.set_n_p = pyo.Set(initialize=range(sets['a_p_shape'][1]), doc='Number of producers') mdl.set_n_c = pyo.Set(initialize=range(sets['a_c_shape'][1]), @@ -100,6 +96,7 @@ def model(matrices: dict, mdl.dirs, mdl.flow, mdl.set_n_i, mdl.set_t, doc='Heat power at the pipes', **pipe_power) + # Building decisions of a pipe mdl.lambda_ = pyo.Var( mdl.dirs, mdl.set_n_i, @@ -120,7 +117,6 @@ def model(matrices: dict, doc='Thermal capacity of the heat source', **source_power) - def heat_source_inst(m, j, t): """Never exceed the installed capacity of the heat source.""" return m.P_source[j, t] <= m.P_source_inst[j] @@ -130,7 +126,6 @@ def heat_source_inst(m, j, t): rule=heat_source_inst, doc='Upper bound for the heat source supply delivery') - def nodal_power_balance(m, j, t): """REFERENCE DIRECTION: left to right P_ji, in P_ji, out @@ -170,7 +165,6 @@ def nodal_power_balance(m, j, t): rule=nodal_power_balance, doc='Nodal Power Balance') - def power_balance_pipe(m, d, j, t): """Power balance for the pipes. @@ -194,7 +188,6 @@ def power_balance_pipe(m, d, j, t): rule=power_balance_pipe, doc='Power balance pipe') - def power_bigm_P(m, d, j, t): lhs = m.P[d, 'in', j, t] - p_max_pipe_const * m.lambda_[d, j] rhs = 0 @@ -203,14 +196,12 @@ def power_bigm_P(m, d, j, t): mdl.dirs, mdl.set_n_i, mdl.set_t, rule=power_bigm_P, doc='Big-M constraint for power flow') - def connection_to_consumer_eco(m, d, j): return m.lambda_[d, j] <= sets[f'lambda_c_{d}'][j] def connection_to_consumer_fcd(m, d, j): return m.lambda_[d, j] == sets[f'lambda_c_{d}'][j] - if optimization_mode == "economic": msg_ = """Constraint if houses have their own connection-pipe and set the direction (ij)""" @@ -281,7 +272,7 @@ def pipes_var(k): # @TODO Implement consumer-specific flh in the economic mode if optimization_mode == "economic": - term4 = (sum( + revenue = (sum( sum( sum(m.lambda_['ij', sets['a_i_in'][j].item()] * matrices['q_c'][k, t] @@ -296,9 +287,9 @@ def pipes_var(k): for t in mdl.set_t) * economics.consumers_flh[0] * economics.heat_price * (-1)) else: - term4 = 0 + revenue = 0 - return fuel + pipes + source + term4 + return fuel + pipes + source + revenue mdl.obj = pyo.Objective(rule=objective_function, doc='Objective function')