From 19247d119cb34a9b81b22323fa6316346b57a2d4 Mon Sep 17 00:00:00 2001 From: jylambert <153106800+jylambert@users.noreply.github.com> Date: Mon, 9 Sep 2024 13:57:25 +0200 Subject: [PATCH 1/3] Removed old mts function for postprocessing --- topotherm/postprocessing.py | 90 ------------------------------------- 1 file changed, 90 deletions(-) diff --git a/topotherm/postprocessing.py b/topotherm/postprocessing.py index 86ab159..865242a 100644 --- a/topotherm/postprocessing.py +++ b/topotherm/postprocessing.py @@ -128,93 +128,3 @@ def equations(v): return res - -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 ununsed - 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) - - # Restart, Adaption of Incidence Matrix for the thermo-hydraulic coupled optimization - for q, _ in enumerate(lambda_built): - if lambda_built[q] == 0: - matrices['a_i'][:, q] = 0 - matrices['l_i'][q] = 0 - elif (lambda_built[q] == 1) & (lambda_dir[q, 0] == 0): - matrices['a_i'][:, q] = matrices['a_i'][:, q] * (-1) - - 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 From 682bafcb6d48aa273c06270d9bc6b4413d65647f Mon Sep 17 00:00:00 2001 From: jylambert <153106800+jylambert@users.noreply.github.com> Date: Mon, 9 Sep 2024 14:06:07 +0200 Subject: [PATCH 2/3] Inserted new method to calculate linear diameter --- topotherm/postprocessing.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/topotherm/postprocessing.py b/topotherm/postprocessing.py index 865242a..294b825 100644 --- a/topotherm/postprocessing.py +++ b/topotherm/postprocessing.py @@ -8,7 +8,7 @@ import numpy as np import pyomo.environ as pyo -from scipy.optimize import fsolve +from scipy.optimize import root from topotherm import settings @@ -113,7 +113,11 @@ def equations(v): 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)) + sol = root(equations, (0.5, 0.02), method='lm') + if sol.success: + v_lin2[h], d_lin2[h] = sol.x + else: + print(h, 'failed to calculate diameter and velocity!') res = dict( a_i=a_i_opt, From 7e9c0fd450a44a9f29b2b709f7c54528a0a6c9bc Mon Sep 17 00:00:00 2001 From: jylambert <153106800+jylambert@users.noreply.github.com> Date: Mon, 9 Sep 2024 17:25:15 +0200 Subject: [PATCH 3/3] Added functionalities for postprocessing: Adapted Q_c if a consumer is not connected. Lamda_dir is now incorporated into the result dictonnary. Added some comments --- topotherm/postprocessing.py | 46 +++++++++++++++++++++++++++++-------- 1 file changed, 36 insertions(+), 10 deletions(-) diff --git a/topotherm/postprocessing.py b/topotherm/postprocessing.py index 294b825..0904bfc 100644 --- a/topotherm/postprocessing.py +++ b/topotherm/postprocessing.py @@ -63,10 +63,26 @@ def postprocess(model, matrices, sets, mode, t_supply, t_return): for index in v: p_cap[index] = pyo.value(v[index]) + # Round lamdda_dir and lambda_built to make sure that hey are integer lambda_dir_1 = np.around(lambda_dir_1, 0) lambda_dir_2 = np.around(lambda_dir_2, 0) lambda_built = np.around(lambda_built, 0) + q_c_real = np.zeros([sets['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 i in sets['connection_c_ij']: + q_c_real[np.where(matrices['a_c'][np.where(matrices['a_i'][:, i] == -1)[0], :][0] == 1)[0], :] = \ + lambda_dir_1[i, 0] * matrices['q_c'][np.where(matrices['a_c'][np.where(matrices['a_i'][:, i] == -1)[0], :][0] == 1)[0], :] + + # Check for consumers connected in direction ji + for i in sets['connection_c_ji']: + q_c_real[np.where(matrices['a_c'][np.where(matrices['a_i'][:, i] == 1)[0], :][0] == 1)[0], :] = \ + lambda_dir_2[i, 0] * matrices['q_c'][np.where(matrices['a_c'][np.where(matrices['a_i'][:, i] == 1)[0], :][0] == 1)[0], :] + + q_c_real = q_c_real[q_c_real.any(axis=1)] # Remove nonzero elements row-wise + # Restart, Adaption of Incidence Matrix for the thermo-hydraulic coupled optimization if mode == "sts": for q, _ in enumerate(lambda_dir_1): @@ -82,6 +98,11 @@ def postprocess(model, matrices, sets, mode, t_supply, t_return): 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 + + lambda_dir_1 = lambda_dir_1[lambda_dir_1.any(axis=1)] + lambda_dir_2 = lambda_dir_2[lambda_dir_2.any(axis=1)] if mode == "sts": p_lin = p_11 + p_21 @@ -96,17 +117,19 @@ def postprocess(model, matrices, sets, mode, t_supply, t_return): 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 + # Prepare variables for the calculation of linear diameters + 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 + supply_temp_opt = np.ones(a_i_shape_opt[1]) * t_supply # Assign constant supply temperature + return_temp_opt = np.ones(a_i_shape_opt[1]) * t_return # Assign constant return temperature + 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 + f = (-1.8 * np.log10((settings.Piping.roughness / (3.7 * d)) ** 1.11 + 6.9 / reynolds))**-2 # friction factor + eq1 = vel - np.sqrt((2 * settings.Piping.max_pr_loss * d) / (f * settings.Water.density)) # eq. for diameter + eq2 = mass_lin - settings.Water.density * vel * (np.pi / 4) * d ** 2 # eq. for velocity return [eq1, eq2] m_lin = (p_lin_opt*1000)/(settings.Water.heat_capacity_cp * (supply_temp_opt - return_temp_opt)) @@ -115,7 +138,7 @@ def equations(v): mass_lin = m_lin[h] sol = root(equations, (0.5, 0.02), method='lm') if sol.success: - v_lin2[h], d_lin2[h] = sol.x + v_lin[h], d_lin[h] = sol.x else: print(h, 'failed to calculate diameter and velocity!') @@ -124,9 +147,12 @@ def equations(v): a_p=a_p_opt, a_c=a_c_opt, q_c=matrices['q_c'], + q_c_con=q_c_real, l_i=l_i_opt, - d_i_0=d_lin2, + d_i_0=d_lin, m_i_0=m_lin, + lambda_dir_1=lambda_dir_1, + lambda_dir_2=lambda_dir_2, position=pos_opt, )