diff --git a/CaseStudy.py b/CaseStudy.py index 87dac63..8f1bb5f 100644 --- a/CaseStudy.py +++ b/CaseStudy.py @@ -1,14 +1,17 @@ import concurrent.futures import copy import os +import typing import warnings from pathlib import Path -from typing import Optional, Self +from typing import Optional, Self, Literal import numpy as np import pandas as pd +import tsam.timeseriesaggregation as tsam import ExcelReader +from InOutModule import Utilities from printer import Printer printer = Printer.getInstance() @@ -250,7 +253,34 @@ def __init__(self, self.scale_CaseStudy() def copy(self): - return copy.deepcopy(self) + new_self = copy.deepcopy(self) + return new_self + + def equal_to(self, cs: typing.Self) -> bool: + """ + Check if this CaseStudy is equal to another CaseStudy, checking all dataframes for equality. + :param cs: Other CaseStudy to compare to + :return: True if equal, False otherwise + """ + all_equal = True + for df in (self.rpk_dependent_dataframes + self.rp_only_dependent_dataframes + self.k_only_dependent_dataframes + self.non_time_dependent_dataframes + self.non_dependent_dataframes): + if hasattr(self, df) and hasattr(cs, df): + self_df = getattr(self, df) + cs_df = getattr(cs, df) + + if type(self_df) is pd.DataFrame and type(cs_df) is pd.DataFrame: + if not self_df.equals(cs_df): + printer.error(f"DataFrame '{df}' is not equal.") + all_equal = False + else: + if self_df != cs_df: + printer.error(f"Attribute '{df}' is not equal.") + all_equal = False + else: + printer.error(f"Attribute '{df}' is missing in one of the CaseStudies.") + all_equal = False + + return all_equal def scale_CaseStudy(self): self.scale_dPower_Parameters() @@ -466,23 +496,31 @@ def get_connected_buses(connection_matrix, bus: str): connected_buses.sort() return connected_buses - def merge_single_node_buses(self): - # Create a connection matrix - # TODO check - connectionMatrix = pd.DataFrame(index=self.dPower_BusInfo.index, columns=[self.dPower_BusInfo.index], data=False) + def merge_single_node_buses(self, inplace: bool = True) -> typing.Optional[typing.Self]: + """ + Merge all buses that are only connected via single-node connections (i.e., "SN" technical representation) into one bus. + :param inplace: Whether to perform the operation inplace or return a new CaseStudy object. + :return: The modified CaseStudy object if inplace is False, otherwise None. + """ + if inplace: + cs = self + else: + cs = self.copy() - for index, entry in self.dPower_Network.iterrows(): + # Create a connection matrix + connectionMatrix = pd.DataFrame(index=cs.dPower_BusInfo.index, columns=[cs.dPower_BusInfo.index], data=False) + for index, entry in cs.dPower_Network.iterrows(): if entry["pTecRepr"] == "SN": - connectionMatrix.loc[index] = True + connectionMatrix.loc[index[0], index[1]] = True connectionMatrix.loc[index[1], index[0]] = True + # Merge buses based on connection matrix merged_buses = set() # Set of buses that have been merged already - for index, entry in connectionMatrix.iterrows(): - if index in merged_buses or entry[entry == True].empty: # Skip if bus has already been merged or has no connections + if index in merged_buses or not entry.any(): # Skip if bus has already been merged or has no connections continue - connected_buses = self.get_connected_buses(connectionMatrix, str(index)) + connected_buses = cs.get_connected_buses(connectionMatrix, str(index)) for bus in connected_buses: merged_buses.add(bus) @@ -490,8 +528,10 @@ def merge_single_node_buses(self): new_bus_name = "merged-" + "-".join(connected_buses) ### Adapt dPower_BusInfo - dPower_BusInfo_entry = self.dPower_BusInfo.loc[connected_buses] # Entry for the new bus + dPower_BusInfo_entry = cs.dPower_BusInfo.loc[connected_buses] # Entry for the new bus zoneOfInterest = 1 if any(dPower_BusInfo_entry['zoi'] == 1) else 0 + zone_values = sorted(set(dPower_BusInfo_entry['z'].dropna().unique())) + zone_name = '_'.join(str(v) for v in zone_values) aggregation_methods_for_columns = { # 'System': 'max', # 'BaseVolt': 'mean', @@ -503,98 +543,89 @@ def merge_single_node_buses(self): 'YearCom': 'mean', 'YearDecom': 'mean', 'lat': 'mean', - 'long': 'mean' + 'lon': 'mean' } dPower_BusInfo_entry = dPower_BusInfo_entry.agg(aggregation_methods_for_columns) dPower_BusInfo_entry['zoi'] = zoneOfInterest + dPower_BusInfo_entry['z'] = zone_name dPower_BusInfo_entry = dPower_BusInfo_entry.to_frame().T dPower_BusInfo_entry.index = [new_bus_name] - self.dPower_BusInfo = self.dPower_BusInfo.drop(connected_buses) + cs.dPower_BusInfo = cs.dPower_BusInfo.drop(connected_buses) with warnings.catch_warnings(): # Suppressing FutureWarning because some entries might include NaN values warnings.simplefilter(action='ignore', category=FutureWarning) - self.dPower_BusInfo = pd.concat([self.dPower_BusInfo, dPower_BusInfo_entry]) + cs.dPower_BusInfo = pd.concat([cs.dPower_BusInfo, dPower_BusInfo_entry]) ### Adapt dPower_Network - self.dPower_Network = self.dPower_Network.reset_index() + cs.dPower_Network = cs.dPower_Network.reset_index() rows_to_drop = [] - for i, row in self.dPower_Network.iterrows(): + for i, row in cs.dPower_Network.iterrows(): if row['i'] in connected_buses and row['j'] in connected_buses: rows_to_drop.append(i) elif row['i'] in connected_buses: row['i'] = new_bus_name - self.dPower_Network.iloc[i] = row + cs.dPower_Network.iloc[i] = row elif row['j'] in connected_buses: row['j'] = new_bus_name - self.dPower_Network.iloc[i] = row - self.dPower_Network = self.dPower_Network.drop(rows_to_drop) + cs.dPower_Network.iloc[i] = row + cs.dPower_Network = cs.dPower_Network.drop(rows_to_drop) # Always put new_bus_name to 'j' (handles case where e.g. 2->3 and 4->2 would lead to 2->34 and 34->2 (because 3 and 4 are merged)) - for i, row in self.dPower_Network.iterrows(): + for i, row in cs.dPower_Network.iterrows(): if row['i'] == new_bus_name: row['i'] = row['j'] row['j'] = new_bus_name - self.dPower_Network.loc[i] = row + cs.dPower_Network.loc[i] = row # Handle case where e.g. 2->3 and 2->4 would lead to 2->34 and 2->34 (because 3 and 4 are merged); also incl. handling 2->3 and 4->2 - self.dPower_Network['Technical Representation'] = self.dPower_Network.groupby(['i', 'j'])['Technical Representation'].transform(lambda series: 'DC-OPF' if 'DC-OPF' in series.values else series.iloc[0]) + cs.dPower_Network['pTecRepr'] = cs.dPower_Network.groupby(['i', 'j'])['pTecRepr'].transform(lambda series: 'DC-OPF' if 'DC-OPF' in series.values else series.iloc[0]) aggregation_methods_for_columns = { # 'Circuit ID': 'first', # 'InService': 'max', # 'R': 'mean', - 'X': lambda x: x.map(lambda a: 1 / a).sum() ** -1, # Formula: 1/X = sum((i,j), 1/Xij)) (e.g., 1/X = 1/Xij_1 +1/Xij_2 + 1/Xij_3...) + 'pXline': lambda x: x.map(lambda a: 1 / a).sum() ** -1, # Formula: 1/X = sum((i,j), 1/Xij)) (e.g., 1/X = 1/Xij_1 +1/Xij_2 + 1/Xij_3...) # 'Bc': 'mean', # 'TapAngle': 'mean', # 'TapRatio': 'mean', - 'Pmax': lambda x: x.min() * x.count(), # Number of lines times the minimum Pmax for new Pmax of the merged lines TODO: Calculate this based on more complex method (flow is relative to R, talk to Benjamin) + 'pPmax': lambda x: x.min() * x.count(), # Number of lines times the minimum Pmax for new Pmax of the merged lines TODO: Calculate this based on more complex method (flow is relative to R, talk to Benjamin) # 'FixedCost': 'mean', # 'FxChargeRate': 'mean', - 'Technical Representation': 'first', - 'LineID': 'first', + 'pTecRepr': 'first', 'YearCom': 'mean', 'YearDecom': 'mean' } - self.dPower_Network = self.dPower_Network.groupby(['i', 'j']).agg(aggregation_methods_for_columns) + # Add aggregation for any missing columns + for column in cs.dPower_Network.columns: + if column not in aggregation_methods_for_columns and column not in ['i', 'j', 'c']: + aggregation_methods_for_columns[column] = 'first' + + cs.dPower_Network = cs.dPower_Network.groupby(['i', 'j', 'c']).agg(aggregation_methods_for_columns) ### Adapt dPower_ThermalGen - if hasattr(self, "dPower_ThermalGen"): - for i, row in self.dPower_ThermalGen.iterrows(): - if row['i'] in connected_buses: - row['i'] = new_bus_name - self.dPower_ThermalGen.loc[i] = row + if hasattr(cs, "dPower_ThermalGen"): + cs.dPower_ThermalGen.loc[cs.dPower_ThermalGen['i'].isin(connected_buses), 'i'] = new_bus_name # Adapt dPower_VRES - if hasattr(self, "dPower_VRES"): - for i, row in self.dPower_VRES.iterrows(): - if row['i'] in connected_buses: - row['i'] = new_bus_name - self.dPower_VRES.loc[i] = row + if hasattr(cs, "dPower_VRES"): + cs.dPower_VRES.loc[cs.dPower_VRES['i'].isin(connected_buses), 'i'] = new_bus_name # Adapt dPower_Storage - if hasattr(self, "dPower_Storage"): - for i, row in self.dPower_Storage.iterrows(): - if row['i'] in connected_buses: - row['i'] = new_bus_name - self.dPower_Storage.loc[i] = row + if hasattr(cs, "dPower_Storage"): + cs.dPower_Storage.loc[cs.dPower_Storage['i'].isin(connected_buses), 'i'] = new_bus_name # Adapt dPower_Demand - self.dPower_Demand = self.dPower_Demand.reset_index() - for i, row in self.dPower_Demand.iterrows(): - if row['i'] in connected_buses: - row['i'] = new_bus_name - self.dPower_Demand.loc[i] = row - self.dPower_Demand = self.dPower_Demand.groupby(['rp', 'i', 'k']).sum() - - # Adapt dPower_VRESProfiles - if hasattr(self, "dPower_VRESProfiles"): - self.dPower_VRESProfiles = self.dPower_VRESProfiles.reset_index() - for i, row in self.dPower_VRESProfiles.iterrows(): - if row['i'] in connected_buses: - row['i'] = new_bus_name - self.dPower_VRESProfiles.loc[i] = row + cs.dPower_Demand = cs.dPower_Demand.reset_index() + mask = cs.dPower_Demand['i'].isin(connected_buses) # Create mask for rows to be unified + cs.dPower_Demand.loc[mask, 'i'] = new_bus_name # Update bus names + aggregation_methods_power_demand = { + 'value': 'sum', + 'dataPackage': lambda v: f"merged-{'-'.join(v.unique())}", + 'dataSource': lambda v: f"merged-{'-'.join(v.unique())}", + 'scenario': lambda v: '-'.join(v.unique()) # If there are multiple scenarios, this would probably fail later (which is good - then we know, something isn't right!) + } + cs.dPower_Demand = cs.dPower_Demand.groupby(['rp', 'k', 'i']).agg(aggregation_methods_power_demand) - self.dPower_VRESProfiles = self.dPower_VRESProfiles.groupby(['rp', 'i', 'k', 'tec']).mean() # TODO: Aggregate using more complex method (capacity * productionCapacity * ... * / Total Production Capacity) - self.dPower_VRESProfiles.sort_index(inplace=True) + return cs if not inplace else None # Create transition matrix from Hindex def get_rpTransitionMatrices(self, clip_method: str = "none", clip_value: float = 0) -> tuple[pd.DataFrame, pd.DataFrame, pd.DataFrame]: @@ -628,7 +659,7 @@ def get_rpTransitionMatrices(self, clip_method: str = "none", clip_value: float raise ValueError(f"For 'relative_to_highest', clip_value must be between 0 and 1, not {clip_value}.") for rp in rps: threshold = rpTransitionMatrixAbsolute.loc[rp].max() * clip_value - rpTransitionMatrixAbsolute.loc[rp][rpTransitionMatrixAbsolute.loc[rp] < threshold] = 0 + rpTransitionMatrixAbsolute.loc[rp, rpTransitionMatrixAbsolute.loc[rp] < threshold] = 0 case _: raise ValueError(f"clip_method must be either 'none', 'absolute_count' or 'relative_to_highest', not {clip_method}.") @@ -731,16 +762,21 @@ def filter_scenario(self, scenario_name, inplace: bool = False) -> Optional[Self return None if inplace else caseStudy - def filter_timesteps(self, start: str, end: str, inplace: bool = False) -> Optional[Self]: + def filter_timesteps(self, start: str, end: str, inplace: bool = False, no_weight_k_adjustment: bool = False) -> Optional[Self]: """ Filters each (relevant) dataframe in the case study to only include the timesteps between start and end (both inclusive). :param start: Start timestep (inclusive). :param end: End timestep (inclusive). :param inplace: If True, modifies the current instance. If False, returns a new instance. + :param no_weight_k_adjustment: If True, does not adjust the weights in dPower_WeightsK after filtering. Adjustment is done by default. :return: None if inplace is True, otherwise a new CaseStudy instance. """ case_study = self if inplace else self.copy() + # Calculate total weight before filtering + if not no_weight_k_adjustment and hasattr(case_study, "dPower_WeightsK") and case_study.dPower_WeightsK is not None: + total_weight = case_study.dPower_WeightsK['pWeight_k'].sum() + for df_name in CaseStudy.k_dependent_dataframes: if hasattr(case_study, df_name) and getattr(case_study, df_name) is not None: df = getattr(case_study, df_name) @@ -756,6 +792,19 @@ def filter_timesteps(self, start: str, end: str, inplace: bool = False) -> Optio setattr(case_study, df_name, filtered_df) + if no_weight_k_adjustment: + printer.information("Skipped adjustment of weights in 'dPower_WeightsK' after filtering timesteps as per user request.") + elif (not hasattr(case_study, "dPower_WeightsK")) or (case_study.dPower_WeightsK is None): + printer.information("Skipped adjustment of weights in 'dPower_WeightsK' after filtering timesteps because 'dPower_WeightsK' does not exist.") + else: + # Adjust weights after filtering + filtered_total_weight = case_study.dPower_WeightsK['pWeight_k'].sum() + if filtered_total_weight == 0: + raise ValueError("After filtering timesteps, the total weight in 'dPower_WeightsK' is zero. Cannot adjust weights.") + adjustment_factor = total_weight / filtered_total_weight + case_study.dPower_WeightsK['pWeight_k'] *= adjustment_factor + printer.information(f"Adjusted weights in 'dPower_WeightsK' by a factor of {adjustment_factor} after filtering timesteps.") + return None if inplace else case_study def filter_representative_periods(self, rp: str, inplace: bool = False) -> Optional[Self]: @@ -768,7 +817,7 @@ def filter_representative_periods(self, rp: str, inplace: bool = False) -> Optio case_study = self if inplace else self.copy() for df_name in CaseStudy.rp_dependent_dataframes: - if hasattr(case_study, df_name): + if hasattr(case_study, df_name) and getattr(case_study, df_name) is not None: df = getattr(case_study, df_name) index = df.index.names @@ -821,3 +870,73 @@ def shift_ks(self, shift: int, inplace: bool = False) -> Optional[Self]: setattr(case_study, df_name, df) return None if inplace else case_study + + def apply_kmedoids_aggregation(self, number_rps: int, rp_length: int = 24, + cluster_strategy: Literal["aggregated", "disaggregated"] = "aggregated", + capacity_normalization: Literal["installed", "maxInvestment"] = "maxInvestment", + sum_production: bool = False, inplace: bool = True, verbose: bool = False) -> Optional[Self]: + """ + Apply k-medoids temporal aggregation to a CaseStudy object. + Each scenario from dGlobal_Scenarios is processed independently. + + :param self: The CaseStudy object to aggregate + :param number_rps: Number of representative periods to create + :param rp_length: Hours per representative period (e.g., 24, 48) + :param cluster_strategy: "aggregated" (sum across buses) or "disaggregated" (keep buses separate) + :param capacity_normalization: "installed" or "maxInvestment" for VRES capacity factor weighting + :param sum_production: If True, sum all technologies into single production column + :param inplace: If True, modify the original CaseStudy; otherwise, return a new one + :param verbose: If True, print detailed processing information + + :return: + CaseStudy: New clustered CaseStudy object if inplace is False; otherwise, None + """ + + cs = self if inplace else self.copy() + Utilities.apply_kmedoids_aggregation(cs, number_rps, rp_length, cluster_strategy, capacity_normalization, sum_production, inplace=True, verbose=verbose) + if inplace: + return None + else: + return cs + + def get_kmedoids_representative_periods(self, number_rps: int, rp_length: int = 24, + cluster_strategy: Literal["aggregated", "disaggregated"] = "aggregated", + capacity_normalization: Literal["installed", "maxInvestment"] = "maxInvestment", + sum_production: bool = False, verbose: bool = False) -> dict[str, tsam.TimeSeriesAggregation]: + """ + Get the representative periods using k-medoids temporal aggregation. Does not modify the original CaseStudy. + Each scenario from dGlobal_Scenarios is processed independently. + + :param self: The CaseStudy object to aggregate + :param number_rps: Number of representative periods to create + :param rp_length: Hours per representative period (e.g., 24, 48) + :param cluster_strategy: "aggregated" (sum across buses) or "disaggregated" (keep buses separate) + :param capacity_normalization: "installed" or "maxInvestment" for VRES capacity factor weighting + :param sum_production: If True, sum all technologies into single production column + :param verbose: If True, print detailed processing information + + :return: TSAM TimeSeriesAggregation object with representative periods for each scenario + """ + + return Utilities.get_kmedoids_representative_periods(self, number_rps, rp_length, cluster_strategy, capacity_normalization, sum_production=sum_production, verbose=verbose) + + def apply_representative_periods(self, representative_periods: dict[str, tsam.TimeSeriesAggregation], rp_length: int = 24, + inplace: bool = True, verbose: bool = False) -> Optional[Self]: + """ + Apply precomputed representative periods to a CaseStudy object. + Each scenario from dGlobal_Scenarios is processed independently. + + :param self: The CaseStudy object to aggregate + :param representative_periods: Precomputed TimeSeriesAggregation object + :param rp_length: Hours per representative period (e.g., 24, 48) + :param inplace: If True, modify the original CaseStudy; otherwise, return a new one + :param verbose: If True, print detailed processing information + :returns: New clustered CaseStudy object if inplace is False; otherwise, None + """ + + cs = self if inplace else self.copy() + Utilities.apply_representative_periods(cs, representative_periods, rp_length, inplace=True, verbose=verbose) + if inplace: + return None + else: + return cs diff --git a/ExcelReader.py b/ExcelReader.py index aba2a04..6ae4ab9 100644 --- a/ExcelReader.py +++ b/ExcelReader.py @@ -332,6 +332,11 @@ def get_Power_Network(excel_file_path: str, keep_excluded_entries: bool = False, """ dPower_Network = __read_non_pivoted_file(excel_file_path, "v0.1.2", ["i", "j", "c"], True, keep_excluded_entries, fail_on_wrong_version) + # Check that all values in column pEnableInvest are either 0 or 1 + if not dPower_Network['pEnableInvest'].isin([0, 1]).all(): + invalid_values = dPower_Network.loc[~dPower_Network['pEnableInvest'].isin([0, 1]), 'pEnableInvest'] + raise ValueError(f"dPower_Network: Found invalid values in 'pEnableInvest' column. Only 0 and 1 are allowed, but found: {invalid_values}") + return dPower_Network diff --git a/ExcelWriter.py b/ExcelWriter.py index 940d0cd..f1ee83e 100644 --- a/ExcelWriter.py +++ b/ExcelWriter.py @@ -1,8 +1,11 @@ +from __future__ import annotations + import os import time import xml.etree.ElementTree as ET from copy import copy, deepcopy from pathlib import Path +from typing import TYPE_CHECKING import numpy as np import openpyxl @@ -12,7 +15,9 @@ import ExcelReader import TableDefinition -from CaseStudy import CaseStudy + +if TYPE_CHECKING: + from CaseStudy import CaseStudy from TableDefinition import CellStyle, Alignment, Font, Color, Text, Column, NumberFormat, TableDefinition from printer import Printer @@ -261,6 +266,8 @@ def write_caseStudy(self, cs: CaseStudy, folder_path: str | Path) -> None: self.write_Power_VRES(cs.dPower_VRES, folder_path) if hasattr(cs, "dPower_VRESProfiles"): self.write_Power_VRESProfiles(cs.dPower_VRESProfiles, folder_path) + if hasattr(cs, "dPower_ImportExport") and cs.dPower_ImportExport is not None: + self.write_Power_ImportExport(cs.dPower_ImportExport, folder_path) self.write_Power_WeightsK(cs.dPower_WeightsK, folder_path) self.write_Power_WeightsRP(cs.dPower_WeightsRP, folder_path) diff --git a/SQLiteWriter.py b/SQLiteWriter.py index 7816821..7547aeb 100644 --- a/SQLiteWriter.py +++ b/SQLiteWriter.py @@ -13,6 +13,8 @@ def model_to_sqlite(model: pyo.base.Model, filename: str) -> None: """ Save the model to a SQLite database. + Automatically includes objective decomposition and dual values. + :param model: Pyomo model to save :param filename: Path to the SQLite database file :return: None @@ -40,10 +42,247 @@ def model_to_sqlite(model: pyo.base.Model, filename: str) -> None: df = pd.DataFrame([pyo.value(o)], columns=['values']) case pyomo.core.base.constraint.ConstraintList | pyomo.core.base.constraint.IndexedConstraint | pyomo.core.base.expression.IndexedExpression: # Those will not be saved on purpose continue + case pyomo.core.base.suffix.Suffix: + if str(o) in ["_relaxed_integer_vars", "dual"]: + continue # Not saved on purpose + else: + printer.warning(f"Pyomo-Type {type(o)} not implemented, {o.name} will not be saved to SQLite") + continue case _: printer.warning(f"Pyomo-Type {type(o)} not implemented, {o.name} will not be saved to SQLite") continue df.to_sql(o.name, cnx, if_exists='replace') cnx.commit() cnx.close() + + # Automatically add objective decomposition and dual values + add_objective_decomposition_to_sqlite(filename, model) + add_dual_values_to_sqlite(filename, model) pass + + +def add_solver_statistics_to_sqlite(filename: str, results, work_units=None) -> None: + """ + Add solver statistics (like Gurobi work-units) to an existing SQLite database. + :param filename: Path to the SQLite database file + :param results: Pyomo solver results object + :param work_units: Optional work units value (from Gurobi solver) + :return: None + """ + cnx = sqlite3.connect(filename) + + # Extract solver statistics + stats = {} + + try: + # Add work units if provided + if work_units is not None: + stats['work_units'] = float(work_units) + + # Get basic solver info from solver[0] + if hasattr(results, 'solver') and len(results.solver) > 0: + solver_info = results.solver[0] + + # Status and termination + if hasattr(solver_info, 'status'): + stats['solver_status'] = str(solver_info.status) + if hasattr(solver_info, 'termination_condition'): + stats['termination_condition'] = str(solver_info.termination_condition) + if hasattr(solver_info, 'time'): + try: + time_val = solver_info.time + if time_val is not None and str(type(time_val)) != "": + stats['solver_time'] = float(time_val) + except Exception: + pass + + # Get problem statistics + if hasattr(results, 'problem'): + problem = results.problem + for attr in ['lower_bound', 'upper_bound', 'number_of_constraints', + 'number_of_variables', 'number_of_nonzeros']: + if hasattr(problem, attr): + value = getattr(problem, attr) + if value is not None: + stats[attr] = float(value) if isinstance(value, (int, float)) else str(value) + + # Create a DataFrame with solver statistics + if stats: + df = pd.DataFrame([stats]) + df.to_sql('solver_statistics', cnx, if_exists='replace', index=False) + cnx.commit() + work_units_str = f"{stats['work_units']:.2f}" if 'work_units' in stats else 'N/A' + printer.information(f"Added solver statistics to SQLite database (work_units: {work_units_str})") + else: + printer.warning("No solver statistics found in results object") + + except Exception as e: + printer.error(f"Failed to add solver statistics: {e}") + import traceback + traceback.print_exc() + finally: + cnx.close() + + +def add_run_parameters_to_sqlite(filename: str, **parameters) -> None: + """ + Add run parameters to an existing SQLite database. + Creates a table 'run_parameters' with parameter names and values. + + :param filename: Path to the SQLite database file + :param parameters: Keyword arguments containing parameter names and values + :return: None + + Example: + add_run_parameters_to_sqlite('model.sqlite', + zoi='R1', + dc_buffer=2, + tp_buffer=1, + scale_demand=1.3, + scale_pmax=1.0) + """ + cnx = sqlite3.connect(filename) + + try: + # Convert parameters to DataFrame + params = {} + for key, value in parameters.items(): + # Convert None to string 'None' for storage + if value is None: + params[key] = 'None' + elif isinstance(value, (int, float)): + params[key] = float(value) + else: + params[key] = str(value) + + if params: + df = pd.DataFrame([params]) + df.to_sql('run_parameters', cnx, if_exists='replace', index=False) + cnx.commit() + printer.information(f"Added run parameters to SQLite database: {', '.join([f'{k}={v}' for k, v in params.items()])}") + else: + printer.warning("No run parameters provided") + + except Exception as e: + printer.error(f"Failed to add run parameters: {e}") + import traceback + traceback.print_exc() + finally: + cnx.close() + + +def add_objective_decomposition_to_sqlite(filename: str, model: pyo.ConcreteModel) -> None: + """ + Add objective function decomposition to SQLite database. + This enables recalculation of ZOI objectives without the full model. + + The objective is decomposed into: + - objective_constant: Single row with the constant term + - objective_terms: Variable names, indices, and their coefficients + + :param filename: Path to the SQLite database file + :param model: Pyomo model with objective + :return: None + """ + from pyomo.repn import generate_standard_repn + + cnx = sqlite3.connect(filename) + + try: + # Decompose objective into linear representation + repn = generate_standard_repn(model.objective.expr, quadratic=False) + + # Store objective decomposition as separate tables + # 1. Constant term + df_constant = pd.DataFrame([{'constant': repn.constant if repn.constant else 0.0}]) + df_constant.to_sql('objective_constant', cnx, if_exists='replace', index=False) + + # 2. Variable names, indices, and coefficients + var_names = [var.parent_component().name for var in repn.linear_vars] + var_indices = [str(var.index()) for var in repn.linear_vars] + coefs = list(repn.linear_coefs) + df_terms = pd.DataFrame({ + 'var_name': var_names, + 'var_index': var_indices, + 'coefficient': coefs + }) + df_terms.to_sql('objective_terms', cnx, if_exists='replace', index=False) + + cnx.commit() + printer.information(f"Added objective decomposition to SQLite ({len(var_indices)} terms)") + + except Exception as e: + printer.error(f"Failed to add objective decomposition: {e}") + import traceback + traceback.print_exc() + finally: + cnx.close() + + +def add_dual_values_to_sqlite(filename: str, model: pyo.ConcreteModel) -> None: + """ + Add dual values (shadow prices) from model constraints to SQLite database. + + Dual values are stored in tables named 'dual_' with: + - Index columns for the constraint + - 'dual_value' column containing the dual/shadow price + + :param filename: Path to the SQLite database file + :param model: Solved Pyomo model with dual suffix + :return: None + """ + cnx = sqlite3.connect(filename) + + try: + if not hasattr(model, 'dual'): + printer.warning("Model does not have dual suffix - no dual values to save") + return + + total_duals = 0 + total_constraints = 0 + + # Iterate through all constraints and save their dual values + for constraint in model.component_objects(pyo.Constraint, active=True): + constraint_name = constraint.name + dual_data = [] + + # Get dual values for this constraint + for index in constraint: + try: + dual_value = model.dual[constraint[index]] + if dual_value is not None: + # Store index and dual value + if isinstance(index, tuple): + # Multi-indexed constraint + row_data = {str(i): val for i, val in enumerate(index)} + row_data['dual_value'] = float(dual_value) + else: + # Single-indexed or scalar constraint + row_data = {'0': index, 'dual_value': float(dual_value)} + dual_data.append(row_data) + total_duals += 1 + except (KeyError, AttributeError): + # Dual value not available for this constraint + pass + + total_constraints += 1 + + # Save to database if we have dual values for this constraint + if dual_data: + df = pd.DataFrame(dual_data) + table_name = f'dual_{constraint_name}' + df.to_sql(table_name, cnx, if_exists='replace', index=False) + + cnx.commit() + + if total_duals > 0: + printer.information(f"Added dual values to SQLite ({total_duals} duals from {total_constraints} constraints)") + else: + printer.warning(f"No dual values found in model (checked {total_constraints} constraints)") + + except Exception as e: + printer.error(f"Failed to add dual values: {e}") + import traceback + traceback.print_exc() + finally: + cnx.close() diff --git a/Utilities.py b/Utilities.py index 4964338..92d94ea 100644 --- a/Utilities.py +++ b/Utilities.py @@ -1,3 +1,4 @@ +import typing from typing import Literal, Dict import numpy as np @@ -6,6 +7,8 @@ from InOutModule.printer import Printer +from InOutModule.CaseStudy import CaseStudy + printer = Printer.getInstance() @@ -89,87 +92,6 @@ def capacityFactorsToInflows(vresProfiles_df: pd.DataFrame, vres_df: pd.DataFram return df.set_index(['rp', 'k', 'g']).sort_index(level="k") -def apply_kmedoids_aggregation( - case_study, - k: int, - rp_length: int = 24, - cluster_strategy: Literal["aggregated", "disaggregated"] = "aggregated", - capacity_normalization: Literal["installed", "maxInvestment"] = "maxInvestment", - sum_production: bool = False -): - """ - Apply k-medoids temporal aggregation to a CaseStudy object. - Each scenario from dGlobal_Scenarios is processed independently. - - Args: - case_study: The CaseStudy object to aggregate - k: Number of representative periods to create - rp_length: Hours per representative period (e.g., 24, 48) - cluster_strategy: "aggregated" (sum across buses) or "disaggregated" (keep buses separate) - capacity_normalization: "installed" or "maxInvestment" for VRES capacity factor weighting - sum_production: If True, sum all technologies into single production column - - Returns: - CaseStudy: New clustered CaseStudy object - """ - - # Create a deep copy to avoid modifying the original - aggregated_case_study = case_study.copy() - - # Get scenario names - scenario_names = aggregated_case_study.dGlobal_Scenarios.index.values - - # Process each scenario independently - all_processed_data = {} - for scenario in scenario_names: - print(f"\n=== Processing scenario: {scenario} ===") - - print(f" Step 1: Extracting data for scenario {scenario}") - scenario_clustering_data = _extract_scenario_data(case_study, scenario, capacity_normalization) - - if len(scenario_clustering_data) == 0: - raise ValueError(f"No data found for scenario {scenario}") - - print(f" Found {len(scenario_clustering_data)} data points for clustering") - - print(f" \nStep 2: Preparing data using {cluster_strategy} strategy") - if cluster_strategy == "disaggregated": - pivot_df = _prepare_disaggregated_data(scenario_clustering_data, sum_production) - else: - pivot_df = _prepare_aggregated_data(scenario_clustering_data, sum_production) - - print(f" Prepared {len(pivot_df)} time periods for clustering") - - print(f" \nStep 3: Running k-medoids clustering (k={k}, rp_length={rp_length})") - aggregation_result = _run_kmedoids_clustering(pivot_df, k, rp_length) - - print(f" \nStep 4: Building representative period data") - data = _build_representative_periods( - case_study, scenario, aggregation_result, rp_length - ) - - print(f" \nStep 5: Building weights and hour indices") - weights_rp, weights_k, hindex = _build_scenario_weights_and_indices( - aggregation_result, scenario, rp_length - ) - - all_processed_data[scenario] = { - 'Power_Demand': data["Power_Demand"], - 'Power_VRESProfiles': data["Power_VRESProfiles"] if "Power_VRESProfiles" in data else [], - 'Power_Inflows': data["Power_Inflows"] if "Power_Inflows" in data else [], - 'weights_rp': weights_rp, - 'weights_k': weights_k, - 'hindex': hindex - } - print(f"Scenario {scenario} completed successfully") - - # Update CaseStudy with aggregated data - _update_casestudy_with_scenarios(aggregated_case_study, all_processed_data) - - print(f"\nAll scenarios have been processed and combined successfully!") - return aggregated_case_study - - def _extract_scenario_data(case_study, scenario: str, capacity_normalization_strategy: str) -> pd.DataFrame: """Extract and combine demand, VRES, and inflows data for a single scenario.""" @@ -251,7 +173,7 @@ def _pivot_technologies(df, value_column, index_cols=None): # Try to merge with Power_VRES data if (hasattr(case_study, 'dPower_VRES') and case_study.dPower_VRES is not None and - 'vres_df' in locals() and len(vres_df) > 0): + vres_with_profiles is not None and len(vres_df) > 0): inflows_with_vres = pd.merge( inflows_df, vres_df[['g', 'tec', 'i', 'ExisUnits', 'EnableInvest', 'MaxInvest']], @@ -353,8 +275,9 @@ def _sum_technology_columns(df: pd.DataFrame) -> pd.DataFrame: return result_df -def _run_kmedoids_clustering(pivot_df: pd.DataFrame, k: int, rp_length: int): +def _run_kmedoids_clustering(pivot_df: pd.DataFrame, k: int, rp_length: int, solver: str = None, verbose: bool = False): """Run k-medoids clustering using tsam.""" + printer = Printer.getInstance() # Prepare data for tsam pivot_df_sorted = pivot_df.sort_values('k') @@ -365,21 +288,25 @@ def _run_kmedoids_clustering(pivot_df: pd.DataFrame, k: int, rp_length: int): # Drop grouping columns and set datetime index clustering_data = pivot_df_sorted.drop(columns=['scenario', 'rp', 'k']).set_index('datetime') - print(f" Running k-medoids with {k} clusters, {rp_length} hours/period, {len(clustering_data)} total hours") + if verbose: + printer.information(f" Running k-medoids with {k} clusters, {rp_length} hours/period, {len(clustering_data)} total hours") # Run clustering - aggregation = tsam.TimeSeriesAggregation( - clustering_data, + tsam_kwargs = dict( noTypicalPeriods=k, hoursPerPeriod=rp_length, clusterMethod='k_medoids', rescaleClusterPeriods=False, - solver="gurobi" ) + if solver is not None: + tsam_kwargs['solver'] = solver + + aggregation = tsam.TimeSeriesAggregation(clustering_data, **tsam_kwargs) typical_periods = aggregation.createTypicalPeriods() - print(f" Clustering completed. Created {len(typical_periods)} typical periods.") - print(f" Cluster center indices (medoids): {aggregation.clusterCenterIndices}") + if verbose: + printer.information(f" Clustering completed. Created {len(typical_periods)} typical periods.") + printer.information(f" Cluster center indices (medoids): {aggregation.clusterCenterIndices}") return aggregation @@ -467,8 +394,9 @@ def _build_scenario_weights_and_indices(aggregation, scenario: str, rp_length: i return weights_rp, weights_k, hindex -def _update_casestudy_with_scenarios(case_study, all_processed_data: Dict): +def _update_casestudy_with_scenarios(case_study, all_processed_data: Dict, verbose: bool = False): """Update CaseStudy with aggregated data, maintaining original index structures.""" + printer = Printer.getInstance() # Collect all data across scenarios all_demand_data = [] @@ -486,36 +414,202 @@ def _update_casestudy_with_scenarios(case_study, all_processed_data: Dict): all_weights_k_data.extend(scenario_data['weights_k']) all_hindex_data.extend(scenario_data['hindex']) - print(f"Updating CaseStudy with combined data:") + if verbose: + printer.information(f"Updating CaseStudy with combined data:") if all_demand_data: demand_df = pd.DataFrame(all_demand_data) - case_study.dPower_Demand = demand_df.set_index(['rp', 'i', 'k']) - print(f" - Updated demand: {len(all_demand_data)} entries") + case_study.dPower_Demand = demand_df.set_index(['rp', 'k', 'i']) + if verbose: + printer.information(f" - Updated demand: {len(all_demand_data)} entries") if all_vres_data: vres_df = pd.DataFrame(all_vres_data) case_study.dPower_VRESProfiles = vres_df.set_index(['rp', 'k', 'g']) - print(f" - Updated VRES profiles: {len(all_vres_data)} entries") + if verbose: + printer.information(f" - Updated VRES profiles: {len(all_vres_data)} entries") if all_inflows_data: inflows_df = pd.DataFrame(all_inflows_data) case_study.dPower_Inflows = inflows_df.set_index(['rp', 'k', 'g']) - print(f" - Updated inflows: {len(all_inflows_data)} entries") + if verbose: + printer.information(f" - Updated inflows: {len(all_inflows_data)} entries") if all_weights_rp_data: weights_rp_df = pd.DataFrame(all_weights_rp_data) case_study.dPower_WeightsRP = weights_rp_df.set_index(['rp']) - print(f" - Updated RP weights: {len(all_weights_rp_data)} entries") + if verbose: + printer.information(f" - Updated RP weights: {len(all_weights_rp_data)} entries") if all_weights_k_data: weights_k_df = pd.DataFrame(all_weights_k_data) case_study.dPower_WeightsK = weights_k_df.set_index(['k']) - print(f" - Updated K weights: {len(all_weights_k_data)} entries") + if verbose: + printer.information(f" - Updated K weights: {len(all_weights_k_data)} entries") if all_hindex_data: hindex_df = pd.DataFrame(all_hindex_data) case_study.dPower_Hindex = hindex_df.set_index(['p', 'rp', 'k']) - print(f" - Updated Hindex: {len(all_hindex_data)} entries") + if verbose: + printer.information(f" - Updated Hindex: {len(all_hindex_data)} entries") + + if verbose: + printer.information("CaseStudy update completed successfully!") + + +def get_kmedoids_representative_periods(case_study, number_rps: int, rp_length: int = 24, + cluster_strategy: Literal["aggregated", "disaggregated"] = "aggregated", + capacity_normalization: Literal["installed", "maxInvestment"] = "maxInvestment", + sum_production: bool = False, solver: str = "gurobi", + verbose: bool = False) -> dict[str, tsam.TimeSeriesAggregation]: + """ + Get the representative periods using k-medoids temporal aggregation. Does not modify the original CaseStudy. + Each scenario from dGlobal_Scenarios is processed independently. + + :param case_study: The CaseStudy object to aggregate + :param number_rps: Number of representative periods to create + :param rp_length: Hours per representative period (e.g., 24, 48) + :param cluster_strategy: "aggregated" (sum across buses) or "disaggregated" (keep buses separate) + :param capacity_normalization: "installed" or "maxInvestment" for VRES capacity factor weighting + :param sum_production: If True, sum all technologies into single production column + :param solver: Solver to use for k-medoids clustering (e.g. "gurobi", "glpk"). Defaults to "gurobi". + :param verbose: If True, print detailed processing information + + :return: TSAM TimeSeriesAggregation object with representative periods for each scenario + """ + + # Get scenario names + scenario_names = case_study.dGlobal_Scenarios.index.values + + # Process each scenario independently + all_scenario_results = {} + for scenario in scenario_names: + printer.information(f"Scenario: {scenario}") if verbose else None + + printer.information(f"Extracting data for scenario {scenario}") if verbose else None + scenario_clustering_data = _extract_scenario_data(case_study, scenario, capacity_normalization) + + if len(scenario_clustering_data) == 0: + raise ValueError(f"No data found for scenario {scenario}") + + printer.information(f"Found {len(scenario_clustering_data)} data points for clustering") if verbose else None + + printer.information(f"Preparing data using {cluster_strategy} strategy") if verbose else None + if cluster_strategy == "disaggregated": + pivot_df = _prepare_disaggregated_data(scenario_clustering_data, sum_production) + else: + pivot_df = _prepare_aggregated_data(scenario_clustering_data, sum_production) + + printer.information(f"Prepared {len(pivot_df)} time periods for clustering") if verbose else None + + printer.information(f"Running k-medoids clustering (k={number_rps}, rp_length={rp_length})") if verbose else None + aggregation_result = _run_kmedoids_clustering(pivot_df, number_rps, rp_length, solver=solver, verbose=verbose) + + printer.information(f"Aggregation result for scenario {scenario} received after {aggregation_result.clusteringDuration} seconds") if verbose else None + all_scenario_results[scenario] = aggregation_result + return all_scenario_results + + +def apply_representative_periods( + case_study, + aggregation: dict[str, tsam.TimeSeriesAggregation], + rp_length: int = 24, + inplace: bool = False, + verbose: bool = False) -> typing.Optional[CaseStudy]: + """ + Apply precomputed representative periods to a CaseStudy object. + Each scenario from dGlobal_Scenarios is processed independently. + + :param case_study: The CaseStudy object to aggregate + :param aggregation: Precomputed TimeSeriesAggregation object + :param rp_length: Hours per representative period (e.g., 24, 48) + :param inplace: If True, modify the original CaseStudy; otherwise, return a new one + :param verbose: If True, print detailed processing information + :returns: New clustered CaseStudy object if inplace is False; otherwise, None + """ + + # Create a deep copy to avoid modifying the original + aggregated_case_study = case_study.copy() if not inplace else case_study + scenario_names = aggregated_case_study.dGlobal_Scenarios.index.values + + # Process each scenario independently + all_processed_data = {} + for scenario in scenario_names: + printer.information(f"Scenario: {scenario}") if verbose else None + + printer.information(f"Building representative period data") if verbose else None + data = _build_representative_periods( + case_study, scenario, aggregation[scenario], rp_length + ) - print("CaseStudy update completed successfully!") + printer.information(f"Building weights and hour indices") if verbose else None + weights_rp, weights_k, hindex = _build_scenario_weights_and_indices( + aggregation[scenario], scenario, rp_length + ) + + all_processed_data[scenario] = { + 'Power_Demand': data["Power_Demand"], + 'Power_VRESProfiles': data["Power_VRESProfiles"] if "Power_VRESProfiles" in data else [], + 'Power_Inflows': data["Power_Inflows"] if "Power_Inflows" in data else [], + 'weights_rp': weights_rp, + 'weights_k': weights_k, + 'hindex': hindex + } + printer.information(f"Scenario {scenario} completed successfully") if verbose else None + + # Update CaseStudy with aggregated data + _update_casestudy_with_scenarios(aggregated_case_study, all_processed_data, verbose=verbose) + + printer.information(f"\nAll scenarios have been processed and combined successfully!") if verbose else None + if not inplace: + return aggregated_case_study + else: + return None + + +def apply_kmedoids_aggregation( + case_study, + k: int, + rp_length: int = 24, + cluster_strategy: Literal["aggregated", "disaggregated"] = "aggregated", + capacity_normalization: Literal["installed", "maxInvestment"] = "maxInvestment", + sum_production: bool = False, + solver: str = "gurobi", + inplace: bool = False, + verbose: bool = False): + """ + Apply k-medoids temporal aggregation to a CaseStudy object. + Each scenario from dGlobal_Scenarios is processed independently. + + :param case_study: The CaseStudy object to aggregate + :param k: Number of representative periods to create + :param rp_length: Hours per representative period (e.g., 24, 48) + :param cluster_strategy: "aggregated" (sum across buses) or "disaggregated" (keep buses separate) + :param capacity_normalization: "installed" or "maxInvestment" for VRES capacity factor weighting + :param sum_production: If True, sum all technologies into single production column + :param solver: Solver to use for k-medoids clustering (e.g. "gurobi", "glpk"). Defaults to "gurobi". + :param inplace: If True, modify the original CaseStudy; otherwise, return a new one + :param verbose: If True, print detailed processing information + + :return: + CaseStudy: New clustered CaseStudy object if inplace is False; otherwise, None + """ + + aggregation_results = get_kmedoids_representative_periods( + case_study, + number_rps=k, + rp_length=rp_length, + cluster_strategy=cluster_strategy, + capacity_normalization=capacity_normalization, + sum_production=sum_production, + solver=solver, + verbose=verbose + ) + + return apply_representative_periods( + case_study, + aggregation=aggregation_results, + rp_length=rp_length, + inplace=inplace, + verbose=verbose + ) diff --git a/environment.yml b/environment.yml index f1e9f44..c066c5b 100644 --- a/environment.yml +++ b/environment.yml @@ -15,6 +15,5 @@ dependencies: - rich=13.9.4 - rich-argparse=1.7.1 - tsam=2.3.9 - - xlrd=2.0.2 - pip: - pulp==2.9.0 diff --git a/printer.py b/printer.py index 5ddc0d0..1517c9c 100644 --- a/printer.py +++ b/printer.py @@ -1,6 +1,9 @@ +from __future__ import annotations + import datetime from rich.console import Console +from rich.markup import escape class Printer: @@ -117,9 +120,9 @@ def error(self, text: str, prefix: str = "Error: ", hard_wrap_chars: str = None) text = self.handle_hard_wrap_chars(text, prefix, hard_wrap_chars) if len(prefix) > 0: - self.console.print(f"[red]{prefix}[/red]{text}") # Only have prefix in color if it is set + self.console.print(f"[red]{escape(prefix)}[/red]{escape(text)}") # Only have prefix in color if it is set else: - self.console.print(f"[red]{text}[/red]") + self.console.print(f"[red]{escape(text)}[/red]") self._log(f"{prefix}{text}") return None @@ -138,9 +141,9 @@ def warning(self, text: str, prefix: str = "Warning: ", hard_wrap_chars: str = N text = self.handle_hard_wrap_chars(text, prefix, hard_wrap_chars) if len(prefix) > 0: - self.console.print(f"[yellow]{prefix}[/yellow]{text}") # Only have prefix in color if it is set + self.console.print(f"[yellow]{escape(prefix)}[/yellow]{escape(text)}") # Only have prefix in color if it is set else: - self.console.print(f"[yellow]{text}[/yellow]") + self.console.print(f"[yellow]{escape(text)}[/yellow]") self._log(f"{prefix}{text}") return None @@ -159,9 +162,9 @@ def success(self, text: str, prefix: str = "", hard_wrap_chars: str = None): text = self.handle_hard_wrap_chars(text, prefix, hard_wrap_chars) if len(prefix) > 0: - self.console.print(f"[green]{prefix}[/green]{text}") # Only have prefix in color if it is set + self.console.print(f"[green]{escape(prefix)}[/green]{escape(text)}") # Only have prefix in color if it is set else: - self.console.print(f"[green]{text}[/green]") + self.console.print(f"[green]{escape(text)}[/green]") self._log(f"{prefix}{text}") return None @@ -179,10 +182,23 @@ def information(self, text: str, prefix: str = "", hard_wrap_chars: str = None): """ text = self.handle_hard_wrap_chars(text, prefix, hard_wrap_chars) - self.console.print(f"{prefix}{text}") + self.console.print(escape(f"{prefix}{text}")) self._log(f"{prefix}{text}") return None + def linear_expression(self, expr) -> None: + """ + Pretty-prints a linear expression to the console and logs it to the logfile if one is set. + + :param expr: The linear expression to be printed + :return: None + """ + + expr_str = " ".join([f"{"+" if coef >= 0 else ""}{coef:.1f}*{var}" for coef, var in zip(expr.linear_coefs, expr.linear_vars)]) + self.console.print(escape(expr_str)) + self._log(expr_str) + return None + def _log(self, text: str) -> None: """ Logs the provided text message to a logfile (by appending) if one is specified.