diff --git a/.claude/rules/claude-md-maintenance.md b/.claude/rules/claude-md-maintenance.md new file mode 100644 index 0000000..aa5c6e5 --- /dev/null +++ b/.claude/rules/claude-md-maintenance.md @@ -0,0 +1,23 @@ +--- +description: Rules for keeping CLAUDE.md and README.md files up to date after code changes in InOutModule/ +--- + +After making any code changes in `InOutModule/`, update the relevant documentation files. + +**`CLAUDE.md`** — context for Claude that is hard to derive from reading the code: non-obvious patterns, gotchas, architectural decisions, naming conventions, and implementation constraints. + +**`README.md`** — human-facing documentation: usage instructions, API overview, key concepts, and examples. + +**Rules:** +- Update the file(s) closest to the code that changed. +- If a change also affects `LEGO/` or root-level code, update the corresponding docs there too (see root `.claude/rules/claude-md-maintenance.md` for the full index). +- Do not duplicate content between `README.md` and `CLAUDE.md` — if something is in `README.md`, `CLAUDE.md` should reference it, not repeat it. +- Use prose references (`See README.md for ...`) to point to large human-facing docs, not `@`-imports. +- Delete stale entries rather than commenting them out. + +**Index of documentation files in InOutModule/:** + +| File | Audience | Contents | +|-------------------------|----------|-------------------------------------------------------------------------| +| `InOutModule/CLAUDE.md` | Claude | CaseStudy read order, ExcelReader/Writer patterns, SQLiteWriter, Caller | +| `InOutModule/README.md` | Humans | CaseStudy API, SQLiteWriter usage, Caller, Excel file format | diff --git a/CLAUDE.md b/CLAUDE.md new file mode 100644 index 0000000..e8974d3 --- /dev/null +++ b/CLAUDE.md @@ -0,0 +1,53 @@ +# CLAUDE.md — InOutModule + +This file provides guidance to Claude Code when working with code in this directory. +See `README.md` for usage, key concepts, and data structure. + +## Architecture Notes + +### CaseStudy + +- Sequential reads happen first (`Global_Parameters`, `Global_Scenarios`, `Power_Parameters`), because subsequent file selection depends on `pEnable*` flags. +- All remaining files are read in parallel via `ThreadPoolExecutor` — order of assignment is non-deterministic, so no file read may depend on another parallel read. +- `dPower_WeightsRP` is **computed** from `dPower_Hindex` (counting occurrences per `rp`); if a `Power_WeightsRP.xlsx` file also exists, it is read and compared — a mismatch triggers a warning but uses the **file** value, not the computed one. +- `merge_single_node_buses()` preserves the `z` (zone) column as a sorted unique union string of all merged zones (e.g. `"R1_R2"`). This is documented in root `CLAUDE.md` as well. +- `merge_generators()` collapses all generators sharing the same `(tec, i)` into one representative generator with ID `"{i}_{tec}"`. It must be called after construction (scaling is not required). VRESProfiles and Inflows are merged **before** VRES so that the original per-generator `MaxProd` weights are available for the capacity-factor weighted average. Generators in VRESProfiles/Inflows that have no matching entry in `dPower_VRES` (left-join miss) are grouped under `(tec=NaN, i=NaN)` — filter to a single scenario first to avoid mixing scenarios in the groupby. +- `filter_zone(zone)` keeps only buses whose `z` column exactly matches the given zone name(s). Cascade order: BusInfo → Network (both endpoints must survive) → ThermalGen / VRES / Storage (by `i`) → Demand (by `i` index level) → VRESProfiles (by surviving VRES gen IDs) → Inflows (by union of surviving VRES + Storage gen IDs) → ImportExport (by `i` index level). `z` is an exact-match filter; merged-bus zone strings like `"R1_R2"` are not matched by `"R1"` alone. +- `shift_transition_matrix(positions, inplace=True, seed=42)` cyclically shifts every row of `rpTransitionMatrixAbsolute` right by `positions` (via `np.roll`, negative shifts left), normalizes, then rebuilds Hindex by Markov chain sampling from the shifted row distributions — anchoring the first period of each scenario to its original RP assignment. Recomputes `dPower_WeightsRP` and all three TM attributes from the new Hindex. The RNG seed defaults to 42. +- `perturb_transition_matrix(randomness, inplace=True, seed=42)` interpolates each row between its original distribution and a random draw: `new_prob = (1-randomness)*orig_prob + randomness*random_draw`, where `random_draw` is `n` uniform values normalized to sum to 1. `randomness=0.0` leaves the matrix unchanged; `randomness=1.0` replaces it fully. Rebuilds Hindex identically to `shift_transition_matrix` (Markov chain sampling, first period anchored, missing-RP correction loop). Uses a single RNG instance for both the random draws and the chain sampling. +- `CaseStudy.copy()` is a full `deepcopy` — safe to modify independently. +- `to_full_hourly_model()` processes each enabled scenario from `dGlobal_Scenarios` independently. It is a no-op (returns unchanged) if all scenarios are already hourly (no two p-values share the same 'k' within any scenario). The produced Hindex uses consecutive h0001…hN p-labels and k0001…kN k-labels per scenario, all mapped to rp01 with weight 1. +- Transition matrices (`rpTransitionMatrixAbsolute`, `rpTransitionMatrixRelativeTo`, `rpTransitionMatrixRelativeFrom`) are computed in the constructor and attached as attributes. + +### ExcelReader + +- Excel sheets whose name starts with `~` are silently skipped (used to disable scenarios in a multi-sheet file without deleting them). +- All Excel files have a version specifier in cell `C2` of each sheet. `check_LEGOExcel_version()` warns (or raises, if `fail_on_wrong_version=True`) on mismatch — wrong version can cause silent column misreads. +- The reader uses `calamine` engine (fast), not `openpyxl`. + +### ExcelWriter + +- All cell styles, column definitions, and table layouts are declared in `TableDefinitions.xml`, not in Python code. When adding a new output table, define its columns there first. +- `ExcelWriter.__init__()` parses the XML once and stores resolved objects (`self.columns`, `self.cell_styles`, etc.). Avoid re-instantiating per row. + +### SQLiteWriter + +- `model_to_sqlite()` automatically calls `add_objective_decomposition_to_sqlite()` and `add_dual_values_to_sqlite()` — these do not need to be called separately. +- `add_run_parameters_to_sqlite()` stores all run configuration in the `run_parameters` table. Evaluation scripts should read from this table (more reliable than filename parsing). +- Pyomo component types not handled by the writer emit a `printer.warning()` and are skipped silently — add new `case` branches if new Pyomo types need to be stored. + +### Utilities + +- `inflowsToCapacityFactors()` joins inflows onto `vresProfiles_df` by dividing by `MaxProd`; generators with missing or zero `MaxProd` are dropped with a warning. +- `capacityFactorsToInflows()` is the inverse; the `remove_Inflows_from_VRESProfiles_inplace` flag modifies the input DataFrame in place when set. +- `plot_transition_matrix()` emits **two** figures per call — a row-normalised ("from" perspective) and a column-normalised ("to" perspective) version — via the inner `_render()` helper. When `output` is given, `os.path.splitext` inserts a `-rowNorm` / `-colNorm` suffix before the extension (e.g. `MK-foo.png` → `MK-foo-rowNorm.png`, `MK-foo-colNorm.png`). Colour intensity and the per-cell percentage encode the respective normalisation; absolute counts and the row/col/grand totals are identical between the two. Cell text flips from black to white where the Blue fill is too dark to read on (perceived luminance `0.299·R + 0.587·G + 0.114·B < 0.5`). It passes `bbox=[0, 0, 1, 1]` to `ax.table()` so the table is forced to fill the entire Axes — without it, matplotlib sizes table rows from font metrics rather than the Axes height, leaving large dead space above/below when the figure is taller than the table's natural size. `fig_w`/`fig_h` are tuned per-cell-inch (not just a generic min-size heuristic) since they now directly determine the rendered cell aspect ratio. + +### Printer + +- `Printer` is a singleton — obtain the instance with `Printer.getInstance()`, never call the constructor directly. +- `set_logfile(path)` redirects all subsequent output to a file (appending). Set to `None` to stop logging. + +### Caller + +- `Caller.py` is a parallel job runner reading from a text file. It uses sentinel files (`.finished{n}`, `.error{n}`) for barrier synchronization across parallel workers. +- Lines containing only `---` act as barriers — workers wait until all prior jobs are complete before continuing past the barrier. diff --git a/Caller.py b/Caller.py new file mode 100644 index 0000000..0124040 --- /dev/null +++ b/Caller.py @@ -0,0 +1,187 @@ +import argparse +import ctypes +import datetime +import os +import subprocess +import sys +import time + +if __name__ == "__main__": + sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__)))) + +from InOutModule.printer import Printer + +printer = Printer.getInstance() + + +def _tail_file(filepath, n=20): + """Return the last n lines of a file, or all if fewer.""" + try: + with open(filepath, 'r', encoding='utf-8', errors='replace') as f: + lines = f.readlines() + except OSError: + return "(could not read log file)" + if not lines: + return "(no output)" + if len(lines) <= n: + return ''.join(lines) + return f"... ({len(lines) - n} lines omitted)\n" + ''.join(lines[-n:]) + + +parser = argparse.ArgumentParser(description='Calls the exact lines from the given file, can be called multiple times.') + +parser.add_argument('jobs', type=str, help='Path to the text-file containing the commands to be called.') +parser.add_argument("--spawn", type=int, help='Number of jobs to spawn (if this is specified, it will call itself multiple times)', nargs='?', default=0) +args = parser.parse_args() +printer.information(f"Using jobs from '{args.jobs}'") + +if args.spawn >= 1: + printer.information(f"Spawning {args.spawn} parallel jobs") + for i in range(args.spawn): + subprocess.Popen([ + "cmd", "/c", "start", f"Caller {i}: {args.jobs}", "cmd", "/k", + f"set POST_ACTIVATE_COMMAND=python {os.path.abspath(__file__)} {args.jobs} && call Conda-Activation-Scripts/activate_environment_windows.bat" + ]) + printer.information(f"Spawned {args.spawn} parallel jobs, exiting... ") + exit(0) + +barrier_waited = {} # tracks total seconds waited per barrier line index; resets when barrier passes or a job is picked up + + +def _all_previous_done(jobs_file, lines, barrier_index): + return all( + os.path.exists(f"{jobs_file}.finished{j}") or os.path.exists(f"{jobs_file}.error{j}") + for j in range(barrier_index) if lines[j].strip() not in ("---", "") + ) + + +def _any_previous_unclaimed(jobs_file, lines, barrier_index): + return any( + not os.path.exists(f"{jobs_file}.started{j}") + and not os.path.exists(f"{jobs_file}.finished{j}") + and not os.path.exists(f"{jobs_file}.error{j}") + for j in range(barrier_index) if lines[j].strip() not in ("---", "") + ) + + +while True: + with open(args.jobs, 'r', encoding='utf-8') as f: + lines = f.readlines() + + found_one = False + restart = False + for i, line in enumerate(lines): + started_job_flag = f"{args.jobs}.started{i}" + finished_job_flag = f"{args.jobs}.finished{i}" + error_job_flag = f"{args.jobs}.error{i}" + if line.strip() == "---": + dot_i = 0 + found_unclaimed = False + last_len = 0 + while not _all_previous_done(args.jobs, lines, i): + if _any_previous_unclaimed(args.jobs, lines, i): + found_unclaimed = True + break + barrier_waited[i] = barrier_waited.get(i, 0) + 3 + dots = "." * ((dot_i % 5) + 1) + msg = f"Barrier '---' at line {i}: Checking every 3s, waited {barrier_waited[i]}s already{dots}" + sys.stdout.write(f"\r{msg}{' ' * max(0, last_len - len(msg))}") + sys.stdout.flush() + last_len = len(msg) + dot_i += 1 + time.sleep(3) + if last_len: + sys.stdout.write('\n') + sys.stdout.flush() + if found_unclaimed: + restart = True + break # restart outer while loop to pick up the unclaimed job + barrier_waited.pop(i, None) + continue # barrier cleared, keep scanning for next job + + try: + fd = os.open(started_job_flag, os.O_CREAT | os.O_EXCL | os.O_WRONLY) + except FileExistsError: + continue # Another worker already claimed this job + if os.path.exists(finished_job_flag) or os.path.exists(error_job_flag): + os.close(fd) + continue + + if barrier_waited: + barrier_waited.clear() + start_datetime = datetime.datetime.now().strftime("%Y-%m-%d %H:%M:%S") + with os.fdopen(fd, 'w', encoding='utf-8') as f: + f.write(f"Command: {line.strip()}\n") + f.write(f"Started at: {start_datetime}") + found_one = True + log_file = f"{args.jobs}.log{i}" + try: + printer.information(f"Executing job {i} from '{args.jobs}': {line.strip()}") + ctypes.windll.kernel32.SetConsoleTitleW(f"Job {i} from '{args.jobs}': {line.strip()}") # type: ignore[attr-defined] - it is Windows-only anyway, so we can ignore the fact that this attribute doesn't exist on other platforms + + start_time = time.time() + with open(log_file, 'w', encoding='utf-8') as log_f: + log_f.write(f"Command: {line.strip()}\n") + log_f.write(f"Started at: {start_datetime}\n") + log_f.write(f"{'=' * 60}\n") + log_f.flush() + + proc = subprocess.Popen( + line.strip(), + shell=True, + stdout=subprocess.PIPE, + stderr=subprocess.STDOUT, + ) + try: + for raw_line in proc.stdout: + decoded = raw_line.decode(errors='replace').replace('\r\n', '\n').replace('\r', '\n') + sys.stdout.write(decoded) + sys.stdout.flush() + log_f.write(decoded) + except KeyboardInterrupt: + # CTRL+C also signals Gurobi, which does a graceful shutdown and saves + # results. Drain remaining output so the log is complete, then let the + # subprocess finish — do not re-raise here. + try: + for raw_line in proc.stdout: + decoded = raw_line.decode(errors='replace').replace('\r\n', '\n').replace('\r', '\n') + sys.stdout.write(decoded) + sys.stdout.flush() + log_f.write(decoded) + except KeyboardInterrupt: + remaining = proc.stdout.read() + if remaining: + decoded = remaining.decode(errors='replace').replace('\r\n', '\n').replace('\r', '\n') + sys.stdout.write(decoded) + log_f.write(decoded) + proc.wait() + end_time = time.time() + + if proc.returncode != 0: + raise RuntimeError( + f"Command exited with code {proc.returncode}. " + f"See log: {log_file}\n" + f"Last output:\n{_tail_file(log_file, 20)}" + ) + + with open(finished_job_flag, 'w', encoding='utf-8') as f: + f.write(f"Command: {line.strip()}\n") + f.write(f"Started at: {start_datetime}\n") + f.write(f"Finished at: {datetime.datetime.now().strftime('%Y-%m-%d %H:%M:%S')}\n") + f.write(f"Execution time: {end_time - start_time:.2f} seconds (= {(end_time - start_time) / 60 / 60:.2f} hours)\n") + f.write(f"Log file: {log_file}\n") + + printer.information(f"Finished job {i} from '{args.jobs}' after {end_time - start_time:.2f} seconds (= {(end_time - start_time) / 60 / 60:.2f} hours).") + except Exception as e: + printer.error(f"Error while executing job {i}: {e}") + with open(error_job_flag, 'w', encoding='utf-8') as f: + f.write(f"Command: {line.strip()}\n") + f.write(f"Error while executing job {i} from '{args.jobs}': {e}\n") + f.write(f"Started at: {start_datetime}\n") + f.write(f"Occurred at: {datetime.datetime.now().strftime('%Y-%m-%d %H:%M:%S')}\n") + f.write(f"Log file: {log_file}\n") + break + + if not found_one and not restart: + printer.information(f"No more jobs to execute in '{args.jobs}', exiting.") + break diff --git a/CaseStudy.py b/CaseStudy.py index 614a128..1ab153f 100644 --- a/CaseStudy.py +++ b/CaseStudy.py @@ -44,11 +44,12 @@ def __init__(self, data_folder: str | Path, do_not_scale_units: bool = False, do_not_merge_single_node_buses: bool = False, + do_not_filter_unused_scenarios: bool = False, parallel_read: bool = True, n_jobs: int = 4, - global_parameters_file: str = "Global_Parameters.xlsx", dGlobal_Parameters: pd.DataFrame = None, + global_parameters_file: str = "Global_Parameters.xlsx", dGlobal_Parameters: dict = None, global_scenarios_file: str = "Global_Scenarios.xlsx", dGlobal_Scenarios: pd.DataFrame = None, - power_parameters_file: str = "Power_Parameters.xlsx", dPower_Parameters: pd.DataFrame = None, + power_parameters_file: str = "Power_Parameters.xlsx", dPower_Parameters: dict = None, power_businfo_file: str = "Power_BusInfo.xlsx", dPower_BusInfo: pd.DataFrame = None, power_network_file: str = "Power_Network.xlsx", dPower_Network: pd.DataFrame = None, power_thermalgen_file: str = "Power_ThermalGen.xlsx", dPower_ThermalGen: pd.DataFrame = None, @@ -65,6 +66,7 @@ def __init__(self, self.data_folder = str(data_folder) if str(data_folder).endswith("/") else str(data_folder) + "/" self.do_not_scale_units = do_not_scale_units self.do_not_merge_single_node_buses = do_not_merge_single_node_buses + self.do_not_filter_unused_scenarios = do_not_filter_unused_scenarios # === SEQUENTIAL READS === if dGlobal_Parameters is not None: @@ -203,22 +205,7 @@ def __init__(self, self.dPower_WeightsRP = dPower_WeightsRP else: self.power_weightsrp_file = power_weightsrp_file - # Calculate dPower_WeightsRP from Hindex - dPower_WeightsRPs = [] - for scenario in self.dPower_Hindex['scenario'].unique().tolist(): - # Count occurences of each value in column 'rp' of dPower_Hindex - dPower_WeightsRP_scenario = pd.DataFrame(self.dPower_Hindex[self.dPower_Hindex['scenario'] == scenario].reset_index()['rp'].value_counts().sort_index()) - dPower_WeightsRP_scenario = dPower_WeightsRP_scenario.rename(columns={'count': 'pWeight_rp'}) - dPower_WeightsRP_scenario['scenario'] = scenario # Add scenario ID - - # Add other columns with default values - dPower_WeightsRP_scenario['id'] = np.nan - dPower_WeightsRP_scenario['dataPackage'] = np.nan - dPower_WeightsRP_scenario['dataSource'] = np.nan - - dPower_WeightsRPs.append(dPower_WeightsRP_scenario) - - dPower_WeightsRP = pd.concat(dPower_WeightsRPs, ignore_index=False) + dPower_WeightsRP = self.calculatePowerWeightsRP(db_id=np.nan, dataPackage=np.nan, dataSource=np.nan) if os.path.exists(self.data_folder + self.power_weightsrp_file): # Compare with given file if it exists self.dPower_WeightsRP = ExcelReader.get_Power_WeightsRP(self.data_folder + self.power_weightsrp_file) @@ -226,21 +213,25 @@ def __init__(self, calculated = dPower_WeightsRP.reset_index().set_index(["rp", "scenario"]) fromFile = self.dPower_WeightsRP.reset_index().set_index(["rp", "scenario"]) - # Normalize both to sum to 1 for comparison - calc_norm = calculated['pWeight_rp'] / calculated['pWeight_rp'].sum() - file_norm = fromFile['pWeight_rp'] / fromFile['pWeight_rp'].sum() # Align indices and fill missing with 0 for comparison - combined = pd.concat([calc_norm, file_norm], axis=1, keys=['calculated', 'fromFile']).fillna(0) + combined = pd.concat([calculated["pWeight_rp"], fromFile["pWeight_rp"]], axis=1, keys=['calculated', 'fromFile']).fillna(0) diff_mask = ~np.isclose(combined['calculated'], combined['fromFile']) if diff_mask.any(): - printer.warning(f"Values for 'pWeight_rp' in `{self.data_folder + self.power_weightsrp_file}` do not match the calculated values based on `{self.power_hindex_file}`. Please check if this is intended, using the file `{self.data_folder + self.power_weightsrp_file}` instead of the calculated values.") + printer.warning(f"Values for 'pWeight_rp' in `{self.data_folder + self.power_weightsrp_file}` do not match the calculated values based on `{self.power_hindex_file}`. Please check if this is intended, now using the file `{self.data_folder + self.power_weightsrp_file}` instead of the calculated values.") # Print all differing lines diffs = combined[diff_mask] printer.warning("Differing entries (index -> calculated | fromFile):\n" + diffs.to_string()) else: # Use calculated dPower_WeightsRP otherwise - printer.warning(f"Executing without 'Power_WeightsRP' (since no file was found at '{self.data_folder + self.power_weightsrp_file}').") + printer.warning(f"Executing without 'Power_WeightsRP' (calculating from Power_Hindex, since no file was found at '{self.data_folder + self.power_weightsrp_file}').") self.dPower_WeightsRP = dPower_WeightsRP + if not self.do_not_filter_unused_scenarios: + self.dGlobal_Scenarios = self.dGlobal_Scenarios[self.dGlobal_Scenarios['relativeWeight'] != 0] # Drop rows in dGlobal_Scenarios where relativeWeight is 0 + if len(self.dGlobal_Scenarios) == 0: + raise ValueError("No scenarios are present in the 'Global_Scenarios' table. Please check if the file exists and contains valid data.") + elif len(self.dGlobal_Scenarios) == 1: + self.filter_scenario(self.dGlobal_Scenarios.index[0], inplace=True) # Filter case study to only include actually present scenarios + self.rpTransitionMatrixAbsolute, self.rpTransitionMatrixRelativeTo, self.rpTransitionMatrixRelativeFrom = self.get_rpTransitionMatrices(clip_method=clip_method, clip_value=clip_value) if not do_not_merge_single_node_buses: @@ -668,6 +659,188 @@ def merge_single_node_buses(self, inplace: bool = True) -> typing.Optional[typin return cs if not inplace else None + def merge_generators(self, inplace: bool = False) -> Optional['CaseStudy']: + """ + Merge generators of the same technology at the same bus into one representative generator. + Affects dPower_ThermalGen, dPower_VRES, dPower_VRESProfiles, and dPower_Inflows. + The new generator ID is '{i}_{tec}'. + + :param inplace: If True, modifies the current instance. If False, returns a new instance. + :return: None if inplace is True, otherwise a new CaseStudy instance. + """ + cs = self if inplace else self.copy() + + # Save original VRES mapping before any merges (needed for VRESProfiles and Inflows weighting) + original_vres_info = None + if hasattr(cs, 'dPower_VRES') and cs.dPower_VRES is not None and 'MaxProd' in cs.dPower_VRES.columns: + original_vres_info = cs.dPower_VRES[['tec', 'i', 'MaxProd']].copy() + + ### Merge dPower_ThermalGen + if hasattr(cs, 'dPower_ThermalGen') and cs.dPower_ThermalGen is not None: + df = cs.dPower_ThermalGen.reset_index() + groups = ['tec', 'i'] + + thermal_simple_agg = { + 'ExisUnits': 'max', + 'MaxProd': 'sum', + 'MinProd': 'min', + 'RampUp': 'sum', + 'RampDw': 'sum', + 'MinUpTime': 'min', + 'MinDownTime': 'min', + 'Qmax': 'sum', + 'Qmin': 'sum', + 'EnableInvest': 'max', + 'YearCom': 'min', + 'YearDecom': 'max', + 'lat': 'mean', + 'lon': 'mean', + } + thermal_weighted_cols = ['InertiaConst', 'FuelCost', 'Efficiency', 'CommitConsumption', + 'OMVarCost', 'StartupConsumption', 'EFOR', 'InvestCost', + 'FirmCapCoef', 'CO2Emis'] + + agg_dict = {} + skip_cols = set(groups + ['g'] + thermal_weighted_cols) + for col in df.columns: + if col in skip_cols: + continue + agg_dict[col] = thermal_simple_agg.get(col, 'first') + + merged = df.groupby(groups).agg(agg_dict).reset_index() + + for col in thermal_weighted_cols: + if col not in df.columns: + continue + numer = (df[col] * df['MaxProd']).groupby([df['tec'], df['i']]).sum() + denom = df['MaxProd'].groupby([df['tec'], df['i']]).sum() + wavg = (numer / denom.replace(0, np.nan)).fillna(df.groupby(groups)[col].mean()) + wavg.name = col + merged = merged.merge(wavg.reset_index(), on=groups, how='left') + + merged['g'] = merged['i'].astype(str) + '_' + merged['tec'] + cs.dPower_ThermalGen = merged.set_index('g') + + ### Merge dPower_VRESProfiles (before dPower_VRES so original MaxProd weights are available) + if (hasattr(cs, 'dPower_VRESProfiles') and cs.dPower_VRESProfiles is not None + and original_vres_info is not None): + df = cs.dPower_VRESProfiles.reset_index() + vres_cols = original_vres_info.reset_index()[['g', 'tec', 'i', 'MaxProd']] + df = df.merge(vres_cols, on='g', how='left') + + groups = ['rp', 'k', 'scenario', 'tec', 'i'] + key = [df['rp'], df['k'], df['scenario'], df['tec'], df['i']] + + numer = (df['value'] * df['MaxProd']).groupby(key).sum() + denom = df['MaxProd'].groupby(key).sum() + merged_value = (numer / denom.replace(0, np.nan)).fillna(df.groupby(groups)['value'].mean()) + merged_value.name = 'value' + + meta_cols = [c for c in ['dataPackage', 'dataSource', 'id'] if c in df.columns] + meta = df.groupby(groups)[meta_cols].first().reset_index() + merged = meta.merge(merged_value.reset_index(), on=groups, how='left') + merged['g'] = merged['i'].astype(str) + '_' + merged['tec'] + merged = merged.drop(columns=['tec', 'i']) + cs.dPower_VRESProfiles = merged.set_index(['rp', 'k', 'g']) + + ### Merge dPower_Inflows + if (hasattr(cs, 'dPower_Inflows') and cs.dPower_Inflows is not None + and original_vres_info is not None): + df = cs.dPower_Inflows.reset_index() + vres_cols = original_vres_info.reset_index()[['g', 'tec', 'i']] + df = df.merge(vres_cols, on='g', how='left') + + groups = ['rp', 'k', 'scenario', 'tec', 'i'] + key = [df['rp'], df['k'], df['scenario'], df['tec'], df['i']] + + merged_value = df['value'].groupby(key).sum() + merged_value.name = 'value' + + meta_cols = [c for c in ['dataPackage', 'dataSource', 'id'] if c in df.columns] + meta = df.groupby(groups)[meta_cols].first().reset_index() + merged = meta.merge(merged_value.reset_index(), on=groups, how='left') + merged['g'] = merged['i'].astype(str) + '_' + merged['tec'] + merged = merged.drop(columns=['tec', 'i']) + cs.dPower_Inflows = merged.set_index(['rp', 'k', 'g']) + + ### Merge dPower_VRES (last, after VRESProfiles and Inflows) + if hasattr(cs, 'dPower_VRES') and cs.dPower_VRES is not None: + df = cs.dPower_VRES.reset_index() + groups = ['tec', 'i'] + + vres_simple_agg = { + 'ExisUnits': 'sum', + 'EnableInvest': 'max', + 'Qmax': 'sum', + 'Qmin': 'sum', + 'YearCom': 'min', + 'YearDecom': 'max', + 'lat': 'mean', + 'lon': 'mean', + } + vres_weighted_cols = ['InvestCost', 'OMVarCost', 'FirmCapCoef', 'InertiaConst'] + special_cols = {'MaxProd', 'MaxInvest'} + + agg_dict = {} + skip_cols = set(groups + ['g'] + vres_weighted_cols + list(special_cols)) + for col in df.columns: + if col in skip_cols: + continue + agg_dict[col] = vres_simple_agg.get(col, 'first') + + merged = df.groupby(groups).agg(agg_dict).reset_index() + + # Special: newMaxProd = sum(ExisUnits * MaxProd) / sum(ExisUnits); fallback to sum when all units are greenfield + if 'MaxProd' in df.columns: + total_mw = (df['ExisUnits'] * df['MaxProd']).groupby([df['tec'], df['i']]).sum() + total_units = df['ExisUnits'].groupby([df['tec'], df['i']]).sum() + new_maxprod = (total_mw / total_units.replace(0, np.nan)).fillna( + df['MaxProd'].groupby([df['tec'], df['i']]).sum() + ) + new_maxprod.name = 'MaxProd' + merged = merged.merge(new_maxprod.reset_index(), on=groups, how='left') + + # Special: newMaxInvest = sum(MaxInvest * MaxProd) / newMaxProd + if 'MaxInvest' in df.columns and 'MaxProd' in df.columns: + invest_mw = (df['MaxInvest'] * df['MaxProd']).groupby([df['tec'], df['i']]).sum() + new_maxprod_s = merged.set_index(groups)['MaxProd'] + new_maxinvest = (invest_mw / new_maxprod_s.replace(0, np.nan)).fillna(0) + new_maxinvest.name = 'MaxInvest' + merged = merged.merge(new_maxinvest.reset_index(), on=groups, how='left') + + for col in vres_weighted_cols: + if col not in df.columns: + continue + numer = (df[col] * df['MaxProd']).groupby([df['tec'], df['i']]).sum() + denom = df['MaxProd'].groupby([df['tec'], df['i']]).sum() + wavg = (numer / denom.replace(0, np.nan)).fillna(df.groupby(groups)[col].mean()) + wavg.name = col + merged = merged.merge(wavg.reset_index(), on=groups, how='left') + + merged['g'] = merged['i'].astype(str) + '_' + merged['tec'] + cs.dPower_VRES = merged.set_index('g') + + return None if inplace else cs + + def calculatePowerWeightsRP(self, db_id, dataPackage, dataSource): + dPower_WeightsRPs = [] + for scenario in self.dPower_Hindex['scenario'].unique().tolist(): + # Count occurrences of each value in column 'rp' of dPower_Hindex + dPower_WeightsRP_scenario = pd.DataFrame(self.dPower_Hindex[self.dPower_Hindex['scenario'] == scenario].reset_index()['rp'].value_counts().sort_index()) + dPower_WeightsRP_scenario = dPower_WeightsRP_scenario.rename(columns={'count': 'pWeight_rp'}) + dPower_WeightsRP_scenario['scenario'] = scenario # Add scenario ID + dPower_WeightsRPs.append(dPower_WeightsRP_scenario) + + dPower_WeightsRP = pd.concat(dPower_WeightsRPs, ignore_index=False) + dPower_WeightsRP['id'] = db_id + dPower_WeightsRP['dataPackage'] = dataPackage + dPower_WeightsRP['dataSource'] = dataSource + + scenario_sums = self.dPower_WeightsK.groupby('scenario')['pWeight_k'].sum() + dPower_WeightsRP['pWeight_rp'] = dPower_WeightsRP['pWeight_rp'] / dPower_WeightsRP['scenario'].map(scenario_sums) + + return dPower_WeightsRP + # Create transition matrix from Hindex def get_rpTransitionMatrices(self, clip_method: str = "none", clip_value: float = 0) -> tuple[pd.DataFrame, pd.DataFrame, pd.DataFrame]: rps = sorted(self.dPower_Hindex.index.get_level_values('rp').unique().tolist()) @@ -716,68 +889,103 @@ def to_full_hourly_model(self, inplace: bool) -> Optional['CaseStudy']: or return a new `CaseStudy` instance if `inplace` is `False`. The adjustments align the data to represent hourly indices and corresponding weights. + Each enabled scenario in `dGlobal_Scenarios` is processed independently. + + If the case study is already a full hourly model (no two p-values share the same 'k' within + any scenario), no adjustment is made and the original instance is returned unchanged. + :param inplace: If `True`, modifies the given instance. If `False`, returns a new `CaseStudy` instance. :return: Adjusted `CaseStudy` instance if `inplace` is `False`, otherwise `None`. """ caseStudy = self.copy() if not inplace else self - # First Adjustment of Hindex (important if the case study was filtered before, to get a coherent p-index) - caseStudy.dPower_Hindex = caseStudy.dPower_Hindex.reset_index() - for i in caseStudy.dPower_Hindex.index: - caseStudy.dPower_Hindex.loc[i, "p"] = f"h{i + 1:0>4}" - caseStudy.dPower_Hindex = caseStudy.dPower_Hindex.set_index(["p", "rp", "k"]) - - # Adjust Demand - adjusted_demand = [] - for i in caseStudy.dPower_BusInfo.index: - for h in caseStudy.dPower_Hindex.index: - adjusted_demand.append(["rp01", h[0].replace("h", "k"), i, caseStudy.dPower_Demand.loc[(h[1], h[2], i), "value"], "ScenarioA", None, None, None]) - - caseStudy.dPower_Demand = pd.DataFrame(adjusted_demand, columns=["rp", "k", "i", "value", "scenario", "id", "dataPackage", "dataSource"]) - caseStudy.dPower_Demand = caseStudy.dPower_Demand.set_index(["rp", "k", "i"]) - - # Adjust VRESProfiles - if hasattr(caseStudy, "dPower_VRESProfiles"): - adjusted_vresprofiles = [] - caseStudy.dPower_VRESProfiles.sort_index(inplace=True) - for g in caseStudy.dPower_VRESProfiles.index.get_level_values('g').unique().tolist(): - for h in caseStudy.dPower_Hindex.index: - adjusted_vresprofiles.append(["rp01", h[0].replace("h", "k"), g, caseStudy.dPower_VRESProfiles.loc[(h[1], h[2], g), "value"], "ScenarioA", None, None, None]) - - caseStudy.dPower_VRESProfiles = pd.DataFrame(adjusted_vresprofiles, columns=["rp", "k", "g", "value", "scenario", "id", "dataPackage", "dataSource"]) - caseStudy.dPower_VRESProfiles = caseStudy.dPower_VRESProfiles.set_index(["rp", "k", "g"]) - - # Adjust Inflows - if hasattr(caseStudy, "dPower_Inflows"): - adjusted_inflows = [] - caseStudy.dPower_Inflows.sort_index(inplace=True) - for g in caseStudy.dPower_Inflows.index.get_level_values('g').unique().tolist(): - for h in caseStudy.dPower_Hindex.index: - adjusted_inflows.append(["rp01", h[0].replace("h", "k"), g, caseStudy.dPower_Inflows.loc[(h[1], h[2], g), "value"], "ScenarioA", None, None, None]) - caseStudy.dPower_Inflows = pd.DataFrame(adjusted_inflows, columns=["rp", "k", "g", "value", "scenario", "id", "dataPackage", "dataSource"]) - caseStudy.dPower_Inflows = caseStudy.dPower_Inflows.set_index(["rp", "k", "g"]) - - # Adjust Hindex - caseStudy.dPower_Hindex = caseStudy.dPower_Hindex.reset_index() - for i in caseStudy.dPower_Hindex.index: - caseStudy.dPower_Hindex.loc[i] = f"h{i + 1:0>4}", f"rp01", f"k{i + 1:0>4}", None, None, None, "ScenarioA" - caseStudy.dPower_Hindex = caseStudy.dPower_Hindex.set_index(["p", "rp", "k"]) - - # Adjust WeightsK - caseStudy.dPower_WeightsK = caseStudy.dPower_WeightsK.reset_index() - caseStudy.dPower_WeightsK = caseStudy.dPower_WeightsK.drop(caseStudy.dPower_WeightsK.index) - for i in range(len(caseStudy.dPower_Hindex)): - caseStudy.dPower_WeightsK.loc[i] = f"{caseStudy.dPower_Hindex.index[i][2]}", None, 1, None, None, "ScenarioA" - caseStudy.dPower_WeightsK = caseStudy.dPower_WeightsK.set_index("k") - - # Adjust WeightsRP - caseStudy.dPower_WeightsRP = caseStudy.dPower_WeightsRP.drop(caseStudy.dPower_WeightsRP.index) - caseStudy.dPower_WeightsRP.loc["rp01"] = None, 1, None, None, "ScenarioA" - - if not inplace: - return caseStudy - else: - return None + # Check if already hourly: within each scenario no 'k' appears more than once + hindex_flat = caseStudy.dPower_Hindex.reset_index() + already_hourly = all( + not (grp.groupby(['k']).size() > 1).any() + for _, grp in hindex_flat.groupby('scenario') + ) + if already_hourly: + return None if inplace else caseStudy + + scenario_names = caseStudy.dGlobal_Scenarios.index.tolist() + + all_demand = [] + all_vresprofiles = [] + all_inflows = [] + all_hindex = [] + all_weightsk = [] + all_weightsrp = [] + + demand_flat = caseStudy.dPower_Demand.reset_index() + vres_flat = caseStudy.dPower_VRESProfiles.reset_index() if hasattr(caseStudy, 'dPower_VRESProfiles') and caseStudy.dPower_VRESProfiles is not None else None + inflows_flat = caseStudy.dPower_Inflows.reset_index() if hasattr(caseStudy, 'dPower_Inflows') and caseStudy.dPower_Inflows is not None else None + + for scenario in scenario_names: + scen_hindex = hindex_flat[hindex_flat['scenario'] == scenario].copy().reset_index(drop=True) + num_hours = len(scen_hindex) + scen_hindex['new_k'] = [f"k{i + 1:0>4}" for i in range(num_hours)] + scen_hindex['new_p'] = [f"h{i + 1:0>4}" for i in range(num_hours)] + + # Demand: merge each original (rp, k) with the buses that have that (rp, k) + demand_scen = demand_flat[demand_flat['scenario'] == scenario][['rp', 'k', 'i', 'value']] + m = scen_hindex[['rp', 'k', 'new_k']].merge(demand_scen, on=['rp', 'k']) + all_demand.append(pd.DataFrame({'rp': 'rp01', 'k': m['new_k'], 'i': m['i'], 'value': m['value'], + 'scenario': scenario, 'id': None, 'dataPackage': None, 'dataSource': None})) + + # VRESProfiles + if vres_flat is not None: + vres_scen = vres_flat[vres_flat['scenario'] == scenario][['rp', 'k', 'g', 'value']] + m = scen_hindex[['rp', 'k', 'new_k']].merge(vres_scen, on=['rp', 'k']) + all_vresprofiles.append(pd.DataFrame({'rp': 'rp01', 'k': m['new_k'], 'g': m['g'], 'value': m['value'], + 'scenario': scenario, 'id': None, 'dataPackage': None, 'dataSource': None})) + + # Inflows + if inflows_flat is not None: + inflows_scen = inflows_flat[inflows_flat['scenario'] == scenario][['rp', 'k', 'g', 'value']] + m = scen_hindex[['rp', 'k', 'new_k']].merge(inflows_scen, on=['rp', 'k']) + all_inflows.append(pd.DataFrame({'rp': 'rp01', 'k': m['new_k'], 'g': m['g'], 'value': m['value'], + 'scenario': scenario, 'id': None, 'dataPackage': None, 'dataSource': None})) + + # Hindex + hindex_scen = pd.DataFrame({ + 'p': scen_hindex['new_p'], + 'rp': 'rp01', + 'k': scen_hindex['new_k'], + 'id': None, + 'dataPackage': None, + 'dataSource': None, + 'scenario': scenario, + }) + all_hindex.append(hindex_scen) + + # WeightsK + weightsk_scen = pd.DataFrame({ + 'k': scen_hindex['new_k'], + 'id': None, + 'pWeight_k': 1, + 'dataPackage': None, + 'dataSource': None, + 'scenario': scenario, + }) + all_weightsk.append(weightsk_scen) + + # WeightsRP + all_weightsrp.append({'rp': 'rp01', 'id': None, 'pWeight_rp': 1, 'dataPackage': None, 'dataSource': None, 'scenario': scenario}) + + caseStudy.dPower_Demand = pd.concat(all_demand).set_index(['rp', 'k', 'i']) + + if all_vresprofiles: + caseStudy.dPower_VRESProfiles = pd.concat(all_vresprofiles).set_index(['rp', 'k', 'g']) + + if all_inflows: + caseStudy.dPower_Inflows = pd.concat(all_inflows).set_index(['rp', 'k', 'g']) + + caseStudy.dPower_Hindex = pd.concat(all_hindex).set_index(['p', 'rp', 'k']) + caseStudy.dPower_WeightsK = pd.concat(all_weightsk).set_index('k') + caseStudy.dPower_WeightsRP = pd.DataFrame(all_weightsrp).set_index('rp') + + return None if inplace else caseStudy def filter_scenario(self, scenario_name, inplace: bool = False) -> Optional[Self]: """ @@ -803,6 +1011,89 @@ def filter_scenario(self, scenario_name, inplace: bool = False) -> Optional[Self return None if inplace else caseStudy + def filter_zone(self, zone: str | list[str], inplace: bool = False) -> Optional[Self]: + """ + Filters the case study to only include buses in the given zone(s). All generators, + network lines, demand entries, and time series profiles connected to buses outside the + zone are removed. + :param zone: Zone name (value of the 'z' column in Power_BusInfo) or list of zone names to keep. + :param inplace: If True, modifies the current instance. If False, returns a new instance. + :return: None if inplace is True, otherwise a new CaseStudy instance. + """ + case_study = self if inplace else self.copy() + + zones = [zone] if isinstance(zone, str) else list(zone) + + # Filter BusInfo and derive remaining bus set + case_study.dPower_BusInfo = case_study.dPower_BusInfo[case_study.dPower_BusInfo['z'].astype(str).isin(zones)] + remaining_buses = set(case_study.dPower_BusInfo.index) + + # Filter Network: drop lines where either endpoint is outside the zone + network_reset = case_study.dPower_Network.reset_index() + case_study.dPower_Network = network_reset[ + network_reset['i'].astype(str).isin(remaining_buses) & network_reset['j'].astype(str).isin(remaining_buses) + ] + case_study.dPower_Network['i'] = case_study.dPower_Network['i'].astype(str) + case_study.dPower_Network['j'] = case_study.dPower_Network['j'].astype(str) + case_study.dPower_Network['c'] = case_study.dPower_Network['c'].astype(str) + case_study.dPower_Network.set_index(['i', 'j', 'c'], inplace=True) + + # Filter ThermalGen + if hasattr(case_study, 'dPower_ThermalGen') and case_study.dPower_ThermalGen is not None: + case_study.dPower_ThermalGen = case_study.dPower_ThermalGen[ + case_study.dPower_ThermalGen['i'].astype(str).isin(remaining_buses) + ] + + # Filter VRES; collect remaining VRES generator IDs for VRESProfiles / Inflows + remaining_vres_gens: set = set() + if hasattr(case_study, 'dPower_VRES') and case_study.dPower_VRES is not None: + case_study.dPower_VRES = case_study.dPower_VRES[ + case_study.dPower_VRES['i'].astype(str).isin(remaining_buses) + ] + remaining_vres_gens = set(case_study.dPower_VRES.index) + + # Filter Storage; collect remaining storage generator IDs for Inflows + remaining_storage_gens: set = set() + if hasattr(case_study, 'dPower_Storage') and case_study.dPower_Storage is not None: + case_study.dPower_Storage = case_study.dPower_Storage[ + case_study.dPower_Storage['i'].astype(str).isin(remaining_buses) + ] + remaining_storage_gens = set(case_study.dPower_Storage.index) + + # Filter Demand + demand_reset = case_study.dPower_Demand.reset_index() + case_study.dPower_Demand = demand_reset[ + demand_reset['i'].astype(str).isin(remaining_buses) + ] + case_study.dPower_Demand['i'] = case_study.dPower_Demand['i'].astype(str) + case_study.dPower_Demand.set_index(['rp', 'k', 'i'], inplace=True) + + # Filter VRESProfiles by remaining VRES generator IDs + if hasattr(case_study, 'dPower_VRESProfiles') and case_study.dPower_VRESProfiles is not None: + profiles_reset = case_study.dPower_VRESProfiles.reset_index() + case_study.dPower_VRESProfiles = profiles_reset[ + profiles_reset['g'].isin(remaining_vres_gens) + ].set_index(['rp', 'k', 'g']) + + # Filter Inflows by remaining VRES + Storage generator IDs + if hasattr(case_study, 'dPower_Inflows') and case_study.dPower_Inflows is not None: + remaining_gens = remaining_vres_gens | remaining_storage_gens + inflows_reset = case_study.dPower_Inflows.reset_index() + case_study.dPower_Inflows = inflows_reset[ + inflows_reset['g'].isin(remaining_gens) + ].set_index(['rp', 'k', 'g']) + + # Filter ImportExport by remaining buses + if hasattr(case_study, 'dPower_ImportExport') and case_study.dPower_ImportExport is not None: + ie_reset = case_study.dPower_ImportExport.reset_index() + case_study.dPower_ImportExport = ie_reset[ + ie_reset['i'].astype(str).isin(remaining_buses) + ] + case_study.dPower_ImportExport['i'] = case_study.dPower_ImportExport['i'].astype(str) + case_study.dPower_ImportExport.set_index(['hub', 'i', 'rp', 'k'], inplace=True) + + return None if inplace else case_study + 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). @@ -912,6 +1203,131 @@ def shift_ks(self, shift: int, inplace: bool = False) -> Optional[Self]: return None if inplace else case_study + def _resample_hindex_from_target_probs(self, target_probs: dict[str, np.ndarray], rps: list[str], rng: np.random.Generator) -> None: + hindex_flat = self.dPower_Hindex.reset_index() + new_parts = [] + + for scenario in hindex_flat['scenario'].unique().tolist(): + sc = hindex_flat[hindex_flat['scenario'] == scenario].copy() + sc = sc.sort_values(['p']) + n_ks_per_rp = len(self.dPower_WeightsK['scenario'] == scenario) + + period_labels = sc['p'].tolist() + first_rp = sc['rp'].iloc[0] + n_periods = len(period_labels) + + # Sample new RP sequence, anchoring the first period to the original + new_rp_seq = [first_rp for _ in range(n_ks_per_rp)] + for _ in range(n_periods - 1): + probs = target_probs[new_rp_seq[-1]] + next_rp = str(rng.choice(rps, p=probs)) + new_rp_seq.extend([next_rp for _ in range(n_ks_per_rp)]) + + sc['rp'] = sc['p'].map(dict(zip(period_labels, new_rp_seq))) + + # Ensure every RP appears at least once; if not, replace the most frequent period with the missing one + rp_count = sc.groupby('rp')['p'].count() / n_ks_per_rp + missing_rps = [rp for rp in rps if rp_count.get(rp, 0) == 0] + while missing_rps: + missing_rp = missing_rps.pop(0) + max_rp = rp_count.idxmax() + if rp_count[max_rp] <= 1: + raise ValueError(f"It seems like there are more RPs than periods in the case study. Check your data and settings.") + + first_occurence_of_max_rp = sc.rp.eq(max_rp).idxmax() + sc.loc[sc.index[first_occurence_of_max_rp:first_occurence_of_max_rp + n_ks_per_rp], 'rp'] = missing_rp + + rp_count.loc[max_rp] -= 1 + rp_count.loc[missing_rp] = 1 + + new_parts.append(sc) + + new_hindex = pd.concat(new_parts, ignore_index=False) + self.dPower_Hindex = new_hindex.set_index(['p', 'rp', 'k']) + + def shift_transition_matrix(self, positions: int, inplace: bool = True, seed: int = 42) -> Optional['CaseStudy']: + """ + Adjust the transition matrix by shifting the probabilities by positions, resampling Hindex from the adjusted + target distribution. dPower_WeightsRP and all three transition-matrices are recomputed from the new Hindex. + + :param positions: Number of positions to shift to the right (negative shifts to the left). + :param inplace: If True, modifies the current instance. If False, returns a new instance. + :param seed: Random seed to guarantee reproduceability. + :return: None if inplace is True, otherwise a new CaseStudy instance. + """ + cs = self if inplace else self.copy() + + rps = cs.rpTransitionMatrixAbsolute.index.tolist() + n = len(rps) + + target_probs: dict[str, np.ndarray] = {} + + printer.information(f"Adjusting transition matrix, shifting it by {positions} positions") + for rp in rps: + c = cs.rpTransitionMatrixAbsolute.loc[rp].values.astype(float) + total = c.sum() + if total == 0: + target_probs[rp] = np.ones(n) / n + else: + c_shifted = np.roll(c, positions) + target_probs[rp] = c_shifted / total + + rng = np.random.default_rng(seed) + cs._resample_hindex_from_target_probs(target_probs, rps, rng) + + # Recompute WeightsRP from new Hindex + cs.dPower_WeightsRP = cs.calculatePowerWeightsRP(np.nan, np.nan, np.nan) + + # Recompute actual TM from new Hindex (approximates the target distributions) + cs.rpTransitionMatrixAbsolute, cs.rpTransitionMatrixRelativeTo, cs.rpTransitionMatrixRelativeFrom = cs.get_rpTransitionMatrices() + + return None if inplace else cs + + def perturb_transition_matrix(self, randomness: float, inplace: bool = True, seed: int = 42) -> Optional['CaseStudy']: + """ + Adjust the transition matrix by interpolating each row between its original distribution and a random draw: + new_prob = (1 - randomness) * orig_prob + randomness * random_draw + Resamples Hindex from the adjusted target distribution. dPower_WeightsRP and all three transition-matrices + are recomputed from the new Hindex. + + :param randomness: Interpolation factor [0.0, 1.0]. 0.0 leaves the matrix unchanged; 1.0 replaces it fully with random draws. + :param inplace: If True, modifies the current instance. If False, returns a new instance. + :param seed: Random seed to guarantee reproducibility. + :return: None if inplace is True, otherwise a new CaseStudy instance. + """ + if not 0.0 <= randomness <= 1.0: + raise ValueError(f"randomness must be in [0.0, 1.0], got {randomness}") + + cs = self if inplace else self.copy() + + rps = cs.rpTransitionMatrixAbsolute.index.tolist() + n = len(rps) + + target_probs: dict[str, np.ndarray] = {} + rng = np.random.default_rng(seed) + + printer.information(f"Adjusting transition matrix, perturbing it with randomness={randomness}") + for rp in rps: + c = cs.rpTransitionMatrixAbsolute.loc[rp].values.astype(float) + total = c.sum() + if total == 0: + target_probs[rp] = np.ones(n) / n + else: + orig_prob = c / total + raw = rng.random(n) + random_draw = raw / raw.sum() + target_probs[rp] = (1 - randomness) * orig_prob + randomness * random_draw + + cs._resample_hindex_from_target_probs(target_probs, rps, rng) + + # Recompute WeightsRP from new Hindex + cs.dPower_WeightsRP = cs.calculatePowerWeightsRP(np.nan, np.nan, np.nan) + + # Recompute actual TM from new Hindex (approximates the target distributions) + cs.rpTransitionMatrixAbsolute, cs.rpTransitionMatrixRelativeTo, cs.rpTransitionMatrixRelativeFrom = cs.get_rpTransitionMatrices() + + return None if inplace else cs + def apply_kmedoids_aggregation(self, number_rps: int, rp_length: int = 24, cluster_strategy: Literal["aggregated", "disaggregated"] = "aggregated", capacity_normalization: Literal["installed", "maxInvestment"] = "maxInvestment", diff --git a/ExcelReader.py b/ExcelReader.py index d940c42..021ea1e 100644 --- a/ExcelReader.py +++ b/ExcelReader.py @@ -66,6 +66,8 @@ def __read_non_pivoted_file(excel_file_path: str, version_specifier: str, indice df = df[df["excl"].isnull()] # Only keep rows that are not excluded (i.e., have no value in the "Excl." column) else: df = df.drop(df.columns[0], axis=1) # Drop the first column (which is empty) + + df[indices] = df[indices].astype(str) # Convert index columns to string to prevent issues with integer-indices df = df.set_index(indices) if len(indices) > 0 else df df["scenario"] = scenario @@ -91,6 +93,7 @@ def __read_pivoted_file(excel_file_path: str, version_specifier: str, indices: l df = __read_non_pivoted_file(excel_file_path, version_specifier, [], has_excl_column, keep_excluded_columns, fail_on_wrong_version) df = df.melt(id_vars=melt_indices + ["scenario"], var_name=pivoted_variable_name, value_name="value") + df[indices] = df[indices].astype(str) # Convert index columns to string to prevent issues with integer-indices (this has to be done here, since we don't pass indices to __read_non_pivoted_file) df = df.set_index(indices) return df diff --git a/README.md b/README.md index 6c553be..9670220 100644 --- a/README.md +++ b/README.md @@ -1 +1,104 @@ -# InOutModule \ No newline at end of file +# InOutModule + +Data I/O package for LEGO-Pyomo. Handles reading Excel case study files, writing results to Excel and SQLite, and utility data transformations. + +## Key Files + +| File | Purpose | +|------------------------|-------------------------------------------------------------------------| +| `CaseStudy.py` | Loads all Excel input files into a single object; data manipulation | +| `ExcelReader.py` | Low-level Excel parsing (per-file reader functions, version checking) | +| `ExcelWriter.py` | Writes formatted Excel output files, driven by `TableDefinitions.xml` | +| `SQLiteWriter.py` | Exports Pyomo model results to SQLite; stores solver stats & run params | +| `Utilities.py` | Data transformations (inflows ↔ capacity factors, Printer helpers) | +| `printer.py` | Singleton console/logfile printer with severity levels | +| `Caller.py` | Parallel job runner for batch experiments | +| `TableDefinitions.xml` | Declarative column/style definitions used by `ExcelWriter` | +| `PypsaReader.py` | Imports case studies from PyPSA network objects | +| `nrel118-reader.py` | Converts NREL 118-bus data into LEGO Excel format | + +## CaseStudy + +### Constructor + +```python +CaseStudy( + data_folder: str | Path, + do_not_scale_units: bool = False, + do_not_merge_single_node_buses: bool = False, + parallel_read: bool = True, + n_jobs: int = 4, + # Per-file overrides: pass a DataFrame to skip reading from disk + dPower_ThermalGen: pd.DataFrame = None, + ... +) +``` + +All Excel files in `data_folder` are read automatically. Any `d*` parameter can be passed directly as a DataFrame to bypass file reading (useful for programmatic construction or testing). + +### DataFrame Attributes + +| Attribute | Time dependency | Source file | +|------------------------|-----------------|------------------------------| +| `dGlobal_Parameters` | none | `Global_Parameters.xlsx` | +| `dGlobal_Scenarios` | none | `Global_Scenarios.xlsx` | +| `dPower_Parameters` | none | `Power_Parameters.xlsx` | +| `dPower_BusInfo` | none | `Power_BusInfo.xlsx` | +| `dPower_Network` | none | `Power_Network.xlsx` | +| `dPower_ThermalGen` | none | `Power_ThermalGen.xlsx` | +| `dPower_VRES` | none | `Power_VRES.xlsx` | +| `dPower_Storage` | none | `Power_Storage.xlsx` | +| `dPower_Demand` | rp + k | `Power_Demand.xlsx` | +| `dPower_VRESProfiles` | rp + k | `Power_VRESProfiles.xlsx` | +| `dPower_Inflows` | rp + k | `Power_Inflows.xlsx` | +| `dPower_ImportExport` | rp + k | `Power_ImportExport.xlsx` | +| `dPower_Hindex` | rp + k | `Power_Hindex.xlsx` | +| `dPower_WeightsRP` | rp only | `Power_WeightsRP.xlsx` | +| `dPower_WeightsK` | k only | `Power_WeightsK.xlsx` | + +### Key Methods + +| Method | Description | +|-------------------------------------|----------------------------------------------------------| +| `copy()` | Deep copy — safe to modify independently | +| `equal_to(cs)` | Compare all DataFrames with another CaseStudy | +| `merge_single_node_buses()` | Collapse single-bus zones; preserves `z` as union string | +| `merge_generators()` | Merge generators of the same `(tec, i)` into one | +| `scale_CaseStudy()` | Applies power and cost scaling factors from parameters | +| `get_rpTransitionMatrices()` | Returns absolute and relative transition matrices | +| `filter_scenario(name)` | Keep only rows matching the given scenario | +| `filter_zone(zone)` | Keep only buses (and all connected data) in a zone | +| `filter_timesteps(start, end)` | Keep only timesteps in the given `k` range | +| `filter_representative_periods(rp)` | Keep only one representative period | +| `apply_kmedoids_aggregation()` | Temporally aggregate using k-medoids clustering | + +## SQLiteWriter + +```python +from InOutModule.SQLiteWriter import model_to_sqlite, add_run_parameters_to_sqlite, add_solver_statistics_to_sqlite + +model_to_sqlite(model, "output/results.sqlite") +add_solver_statistics_to_sqlite("output/results.sqlite", lego) +add_run_parameters_to_sqlite("output/results.sqlite", zoi="R1", dc_buffer=2) +``` + +- `model_to_sqlite()` exports all Pyomo variables, parameters, and sets; automatically appends objective decomposition and dual values. +- `add_run_parameters_to_sqlite()` creates a `run_parameters` table — evaluation scripts should read from this table rather than parsing filenames. + +## Caller (Parallel Job Runner) + +```bash +python InOutModule/Caller.py jobs.txt +python InOutModule/Caller.py jobs.txt --spawn 4 # open 4 parallel terminal windows +``` + +`jobs.txt` is a plain-text file with one shell command per line. Lines containing only `---` act as barriers: all workers wait until every job above the barrier is finished before continuing. + +## Excel File Format + +All input Excel files follow a versioned multi-sheet format: +- Each sheet corresponds to one scenario (or `ScenarioA` for deterministic runs). +- Sheets whose name starts with `~` are skipped. +- Cell `C2` on each sheet contains a version specifier (e.g., `v0.1.0`); mismatches produce a warning. + +See `changelog-LEGOExcels.md` for version history of the Excel format. diff --git a/SQLiteWriter.py b/SQLiteWriter.py index b55de38..ff266c6 100644 --- a/SQLiteWriter.py +++ b/SQLiteWriter.py @@ -6,13 +6,14 @@ import pyomo.environ as pyo from InOutModule.printer import Printer +from LEGO.LEGO import LEGO printer = Printer.getInstance() -def model_to_sqlite(model: pyo.base.Model, filename: str) -> None: +def model_to_sqlite(model: pyo.base.ConcreteModel, filename: str) -> None: """ - Save the model to a SQLite database. + Save the model to an SQLite database. Automatically includes objective decomposition and dual values. :param model: Pyomo model to save @@ -61,14 +62,17 @@ def model_to_sqlite(model: pyo.base.Model, filename: str) -> None: pass -def add_solver_statistics_to_sqlite(filename: str, results, work_units=None) -> None: +def add_solver_statistics_to_sqlite(filename: str, lego: LEGO) -> None: """ - Add solver statistics (like Gurobi work-units) to an existing SQLite database. + Add solver statistics to an existing SQLite database, reading from a LEGO instance. :param filename: Path to the SQLite database file - :param results: Pyomo solver results object - :param work_units: Optional work units value (from Gurobi solver) + :param lego: Solved LEGO instance (reads lego.results, lego.work_units, lego.mip_gap) :return: None """ + results = lego.results + work_units = lego.work_units + mip_gap = lego.mip_gap + cnx = sqlite3.connect(filename) # Extract solver statistics @@ -106,13 +110,31 @@ def add_solver_statistics_to_sqlite(filename: str, results, work_units=None) -> if value is not None: stats[attr] = float(value) if isinstance(value, (int, float)) else str(value) + # MIP gap + if mip_gap is not None: + stats['mip_gap'] = float(mip_gap) + else: + # Try to derive from bounds as fallback + lb = stats.get('lower_bound') + ub = stats.get('upper_bound') + if lb is not None and ub is not None and isinstance(lb, float) and isinstance(ub, float) and ub != 0: + derived_gap = abs(ub - lb) / abs(ub) + if derived_gap > 1e-10: + stats['mip_gap'] = derived_gap + printer.information(f"MIP gap derived from bounds: {derived_gap:.6f}") + else: + printer.information("MIP gap not stored: bounds are equal (LP or MIP solved to optimality)") + else: + printer.information("MIP gap not available: not provided and bounds not found in results") + # 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})") + mip_gap_str = f"{stats['mip_gap']:.6f}" if 'mip_gap' in stats else 'N/A' + printer.information(f"Added solver statistics to SQLite database (work_units: {work_units_str}, mip_gap: {mip_gap_str})") else: printer.warning("No solver statistics found in results object") diff --git a/Utilities.py b/Utilities.py index 2585430..712835a 100644 --- a/Utilities.py +++ b/Utilities.py @@ -1,8 +1,10 @@ from __future__ import annotations +import os import typing from typing import TYPE_CHECKING, Literal, Dict +import matplotlib.pyplot as plt import numpy as np import pandas as pd import tsam.timeseriesaggregation as tsam @@ -224,6 +226,7 @@ def _pivot_technologies(df, value_column, index_cols=None): if combined_tech_data is not None: # Use right join to keep ALL demand data (even nodes without technology data) # Replicates demand for nodes with technology, and preserves demand-only nodes + combined_tech_data['i'] = combined_tech_data['i'].astype(str) scenario_df = pd.merge( combined_tech_data, scenario_df, @@ -321,7 +324,9 @@ def _extract_numeric_and_calc_p(df, rp_length): """Extract numeric values from rp/k strings and calculate absolute hour.""" df['rp_num'] = df['rp'].str[2:].astype(int) df['k_num'] = df['k'].str[1:].astype(int) - df['p'] = (df['rp_num'] - 1) * rp_length + df['k_num'] + # Normalize k so the first k value maps to 1, regardless of offset (e.g. --limitK k2161-k4320) + min_k_num = df['k_num'].min() + df['p'] = (df['rp_num'] - 1) * rp_length + (df['k_num'] - min_k_num + 1) return df time_series_tables = [("Power_Demand", case_study.dPower_Demand)] @@ -616,3 +621,125 @@ def apply_kmedoids_aggregation( inplace=inplace, verbose=verbose ) + + +def plot_transition_matrix(tm: pd.DataFrame, title: str | None = None, output: str | None = None): + """Plot a transition matrix as row- and column-normalised Blue-tinted tables. + + Two figures are produced: one normalised by row sum ("from" perspective) and + one normalised by column sum ("to" perspective). In each, cell text shows + absolute counts and the respective normalised percentage; cell colour + intensity encodes that probability (darker = more likely). Row and column + sums are appended in bold. + + :param tm: Transition matrix DataFrame (index and columns are RP labels). + Typically ``cs.rpTransitionMatrixAbsolute``. + :param title: Optional subtitle shown below the main heading. + :param output: If given, save the figures to this path with ``-rowNorm`` / + ``-colNorm`` inserted before the extension, instead of + displaying them. + """ + labels = list(tm.index) + n = len(labels) + data = tm.values.astype(float) + + row_totals = data.sum(axis=1) # shape (n,) + col_totals = data.sum(axis=0) # shape (n,) + grand_total = data.sum() + + row_rel = data / np.where(row_totals == 0, 1, row_totals)[:, np.newaxis] + col_rel = data / np.where(col_totals == 0, 1, col_totals)[np.newaxis, :] + + def _render(rel: np.ndarray, pct_label: str, out_path: str | None): + title_h = 0.55 if title else 0.35 + fig_w = max(4.0, 0.95 * (n + 1) + 0.3) + fig_h = 0.4 * (n + 2) + title_h + fig, ax = plt.subplots(figsize=(fig_w, fig_h)) + ax.set_axis_off() + + cmap = plt.cm.Blues + + cell_text = [ + [f"{int(data[r, c])}\n{rel[r, c]:.0%}" for c in range(n)] + [f"{int(row_totals[r])}"] + for r in range(n) + ] + cell_text.append([f"{int(col_totals[c])}" for c in range(n)] + [f"{int(grand_total)}"]) + + cell_colors = [ + [cmap(rel[r, c]) for c in range(n)] + ["white"] + for r in range(n) + ] + cell_colors.append(["white"] * (n + 1)) + + sum_labels = labels + ["Sum"] + + tbl = ax.table( + cellText=cell_text, + rowLabels=sum_labels, + colLabels=sum_labels, + cellColours=cell_colors, + loc="center", + cellLoc="center", + bbox=[0, 0, 1, 1], + ) + tbl.auto_set_font_size(False) + tbl.set_fontsize(8) + + for r in range(n + 2): + tbl[r, n].get_text().set_fontweight("bold") + for c in range(-1, n + 1): + tbl[n + 1, c].get_text().set_fontweight("bold") + + # Flip to white text where the Blue shade is too dark to read black on. + # (data cells are table rows 1..n, cols 0..n-1; row 0 is the header) + for r in range(n): + for c in range(n): + red, green, blue, _ = cmap(rel[r, c]) + luminance = 0.299 * red + 0.587 * green + 0.114 * blue + if luminance < 0.5: + tbl[r + 1, c].get_text().set_color("white") + + heading = f"Transition matrix (count / {pct_label})" + if title: + heading += f"\n{title}" + ax.set_title(heading, fontsize=10, pad=8) + fig.tight_layout(pad=0.3) + + # Force layout so cell bounding boxes are finalised + fig.canvas.draw() + renderer = fig.canvas.get_renderer() + + def _to_axes(xd, yd): + return ax.transAxes.inverted().transform((xd, yd)) + + bb0 = tbl[n + 1, 0].get_window_extent(renderer) + bb_rc = tbl[n + 1, n].get_window_extent(renderer) + bb_rl = tbl[n + 1, -1].get_window_extent(renderer) + x0_ax, y_ax = _to_axes(bb_rl.x0, bb0.y1) + x1_ax, _ = _to_axes(bb_rc.x1, bb0.y1) + ax.plot([x0_ax, x1_ax], [y_ax, y_ax], transform=ax.transAxes, + color="black", linewidth=2, clip_on=False, zorder=10) + + bb_hdr = tbl[0, n].get_window_extent(renderer) + bb_bot = tbl[n + 1, n].get_window_extent(renderer) + x_ax, y0_ax = _to_axes(bb_hdr.x0, bb_bot.y0) + _, y1_ax = _to_axes(bb_hdr.x0, bb_hdr.y1) + ax.plot([x_ax, x_ax], [y0_ax, y1_ax], transform=ax.transAxes, + color="black", linewidth=2, clip_on=False, zorder=10) + + if out_path: + fig.savefig(out_path, dpi=150, bbox_inches="tight") + print(f"Saved to {out_path}") + else: + plt.show() + plt.close(fig) + + for rel_matrix, norm_suffix, norm_label in ( + (row_rel, "rowNorm", "row %"), + (col_rel, "colNorm", "column %"), + ): + target = None + if output: + base, ext = os.path.splitext(output) + target = f"{base}-{norm_suffix}{ext}" + _render(rel_matrix, norm_label, target) diff --git a/printer.py b/printer.py index 19baf42..e7f5b93 100644 --- a/printer.py +++ b/printer.py @@ -222,6 +222,37 @@ def _log(self, text: str) -> None: f.write(line + "\n") return None + def note(self, text: str, prefix: str = "Note: ", hard_wrap_chars: str = None) -> None: + """ + Handles a note message. Printed in yellow (same as warning but semantically different). + + :param text: Text to be printed + :param prefix: Prefix (default: "Note: ") + :param hard_wrap_chars: Chars to be added at the end of the text if it exceeds the console width, + None if text should not be truncated + :return: None + """ + text = self.handle_hard_wrap_chars(text, prefix, hard_wrap_chars) + if len(prefix) > 0: + self.console.print(f"[yellow]{escape(prefix)}[/yellow]{escape(text)}") + else: + self.console.print(f"[yellow]{escape(text)}[/yellow]") + self._log(f"{prefix}{text}") + return None + + def rule(self, title: str = "", style: str = "white") -> None: + """ + Prints a horizontal rule (line) to the console, optionally with a centered title. + Delegates to rich.console.Console.rule(). + + :param title: Optional title to display in the center of the rule + :param style: Rich style string for the rule (default: "white") + :return: None + """ + self.console.rule(title, style=style) + self._log(f"--- {title} ---" if title else "-" * self.console.width) + return None + def separator(self) -> None: """ Prints a separator line to the console. The line is made up of dashes and @@ -285,3 +316,8 @@ def pprint_zoi_var(var, zoi, index_positions: list = None, decimals: int = 2): # Iterate over all lists and print the values for i in range(len(value_list)): print(f" {key_list[i]:>{key_spacer}} : {lower_list[i]:>{lower_spacer}} : {value_list[i]:>{value_spacer}} : {upper_list[i]:>{upper_spacer}} : {fixed_list[i]:>{fixed_spacer}} : {stale_list[i]:>{stale_spacer}} : {domain_list[i]:>{domain_spacer}}") + + +def timestamp() -> str: + """Returns a timestamp string in the format 'YYYYMMDD_HHMMSS'.""" + return datetime.datetime.now().strftime('%Y%m%d_%H%M%S')