From 8ceeba9b67e8b6af369efc4ed15150726d51d6ed Mon Sep 17 00:00:00 2001 From: Amos Schledorn Date: Fri, 9 Jan 2026 15:15:03 +0100 Subject: [PATCH 1/4] create first working version --- config/config.default.yaml | 6 ++ rules/collect.smk | 28 +++++ rules/solve_electricity.smk | 39 +++++++ rules/solve_overnight.smk | 40 +++++++ scripts/solve_network.py | 201 ++++++++++++++++++++++++++++++++++-- 5 files changed, 305 insertions(+), 9 deletions(-) diff --git a/config/config.default.yaml b/config/config.default.yaml index c1eaa966ec..8ee941f0fc 100644 --- a/config/config.default.yaml +++ b/config/config.default.yaml @@ -1002,6 +1002,12 @@ adjustments: capital_cost: 1.0 absolute: false +mga: + enable: false + epsilon: [0.05] + sense: [min, max] + runs: + # docs in https://pypsa-eur.readthedocs.io/en/latest/configuration.html#data data: hotmaps_industrial_sites: diff --git a/rules/collect.smk b/rules/collect.smk index f34b4050b3..754f839341 100644 --- a/rules/collect.smk +++ b/rules/collect.smk @@ -9,7 +9,9 @@ localrules: prepare_elec_networks, prepare_sector_networks, solve_elec_networks, + solve_elec_mga_networks, solve_sector_networks, + solve_sector_mga_networks, rule process_costs: @@ -68,6 +70,19 @@ rule solve_elec_networks: ), +rule solve_elec_mga_networks: + input: + expand( + RESULTS + + "networks/mga/base_s_{clusters}_elec_{opts}_{mga_run}_eps{epsilon}_sense{sense}.nc", + **config["scenario"], + mga_run=config["mga"]["runs"], + epsilon=config["mga"]["epsilon"], + sense=config["mga"]["sense"], + run=config["run"]["name"], + ), + + rule solve_sector_networks: input: expand( @@ -78,6 +93,19 @@ rule solve_sector_networks: ), +rule solve_sector_mga_networks: + input: + expand( + RESULTS + + "networks/mga/base_s_{clusters}_{opts}_{sector_opts}_{planning_horizons}_{mga_run}_eps{epsilon}_sense{sense}.nc", + **config["scenario"], + mga_run=config["mga"]["runs"], + epsilon=config["mga"]["epsilon"], + sense=config["mga"]["sense"], + run=config["run"]["name"], + ), + + rule solve_sector_networks_perfect: input: expand( diff --git a/rules/solve_electricity.smk b/rules/solve_electricity.smk index e94db1e904..9088335168 100644 --- a/rules/solve_electricity.smk +++ b/rules/solve_electricity.smk @@ -11,6 +11,7 @@ rule solve_network: "sector", "co2_sequestration_potential", default=200 ), custom_extra_functionality=input_custom_extra_functionality, + mga=config_provider("mga"), input: network=resources("networks/base_s_{clusters}_elec_{opts}.nc"), output: @@ -34,6 +35,44 @@ rule solve_network: "../scripts/solve_network.py" +rule solve_mga_network: + params: + solving=config_provider("solving"), + foresight=config_provider("foresight"), + co2_sequestration_potential=config_provider( + "sector", "co2_sequestration_potential", default=200 + ), + custom_extra_functionality=input_custom_extra_functionality, + mga=config_provider("mga"), + input: + network=RESULTS + "networks/base_s_{clusters}_elec_{opts}.nc", + output: + network=RESULTS + + "networks/mga/base_s_{clusters}_elec_{opts}_{mga_run}_eps{epsilon}_sense{sense}.nc", + log: + solver=normpath( + RESULTS + + "logs/solve_mga_network/base_s_{clusters}_elec_{opts}_{mga_run}_eps{epsilon}_sense{sense}_solver.log" + ), + memory=RESULTS + + "logs/solve_mga_network/base_s_{clusters}_elec_{opts}_{mga_run}_eps{epsilon}_sense{sense}_memory.log", + python=RESULTS + + "logs/solve_mga_network/base_s_{clusters}_elec_{opts}_{mga_run}_eps{epsilon}_sense{sense}_python.log", + benchmark: + ( + RESULTS + + "benchmarks/solve_mga_network/base_s_{clusters}_elec_{opts}_{mga_run}_eps{epsilon}_sense{sense}" + ) + threads: solver_threads + resources: + mem_mb=memory, + runtime=config_provider("solving", "runtime", default="6h"), + shadow: + shadow_config + script: + "../scripts/solve_network.py" + + rule solve_operations_network: params: options=config_provider("solving", "options"), diff --git a/rules/solve_overnight.smk b/rules/solve_overnight.smk index 8a4830e0c4..fca95b4879 100644 --- a/rules/solve_overnight.smk +++ b/rules/solve_overnight.smk @@ -2,6 +2,8 @@ # # SPDX-License-Identifier: MIT +import math + rule solve_sector_network: params: @@ -11,6 +13,7 @@ rule solve_sector_network: "sector", "co2_sequestration_potential", default=200 ), custom_extra_functionality=input_custom_extra_functionality, + mga=config_provider("mga"), input: network=resources( "networks/base_s_{clusters}_{opts}_{sector_opts}_{planning_horizons}.nc" @@ -40,3 +43,40 @@ rule solve_sector_network: ) script: "../scripts/solve_network.py" + + +rule solve_sector_mga_network: + params: + solving=config_provider("solving"), + foresight=config_provider("foresight"), + co2_sequestration_potential=config_provider( + "sector", "co2_sequestration_potential", default=200 + ), + custom_extra_functionality=input_custom_extra_functionality, + mga=config_provider("mga"), + input: + network=RESULTS + + "networks/base_s_{clusters}_{opts}_{sector_opts}_{planning_horizons}.nc", + output: + network=RESULTS + + "networks/mga/base_s_{clusters}_{opts}_{sector_opts}_{planning_horizons}_{mga_run}_eps{epsilon}_sense{sense}.nc", + log: + solver=RESULTS + + "logs/solve_sector_mga_network/base_s_{clusters}_{opts}_{sector_opts}_{planning_horizons}_{mga_run}_eps{epsilon}_sense{sense}_solver.log", + memory=RESULTS + + "logs/solve_sector_mga_network/base_s_{clusters}_{opts}_{sector_opts}_{planning_horizons}_{mga_run}_eps{epsilon}_sense{sense}_memory.log", + python=RESULTS + + "logs/solve_sector_mga_network/base_s_{clusters}_{opts}_{sector_opts}_{planning_horizons}_{mga_run}_eps{epsilon}_sense{sense}_python.log", + threads: solver_threads + resources: + mem_mb=config_provider("solving", "mem_mb"), + runtime=config_provider("solving", "runtime", default="6h"), + benchmark: + ( + RESULTS + + "benchmarks/solve_sector_mga_network/base_s_{clusters}_{opts}_{sector_opts}_{planning_horizons}_{mga_run}_eps{epsilon}_sense{sense}" + ) + shadow: + shadow_config + script: + "../scripts/solve_network.py" diff --git a/scripts/solve_network.py b/scripts/solve_network.py index 9d5fcb4a8b..6464843032 100644 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -1270,6 +1270,159 @@ def check_objective_value(n: pypsa.Network, solving: dict) -> None: ) +def build_mga_weights(n: pypsa.Network, run_config: dict) -> dict: + """ + Build MGA weights dictionary from regex patterns. + + Parameters + ---------- + n : pypsa.Network + The PyPSA network instance + run_config : dict + Single MGA run config with component types and their parameters. + Example: {"Generator": {"regex": ".*wind", "variable": "p_nom", "weights": 1}} + + Returns + ------- + dict + Weights dictionary for n.optimize.optimize_mga() + """ + weights = {} + component_map = { + "Generator": n.generators, + "Link": n.links, + "StorageUnit": n.storage_units, + "Store": n.stores, + "Line": n.lines, + } + + for component_type, cfg in run_config.items(): + if component_type not in component_map: + logger.warning(f"Unknown component type '{component_type}' in MGA config") + continue + + df = component_map[component_type] + regex = cfg.get("regex", "") + variable = cfg.get("variable", "p_nom") + weight = cfg.get("weights", 1.0) + + matching = df.index[df.index.str.contains(regex, regex=True, na=False)] + + if matching.empty: + logger.warning(f"No {component_type}s found matching pattern '{regex}'") + continue + + logger.info( + f"Found {len(matching)} {component_type}s matching pattern '{regex}'" + ) + + weights.setdefault(component_type, {})[variable] = { + idx: weight for idx in matching + } + + return weights + + +def run_mga_optimization( + n: pypsa.Network, + mga_params: dict, + mga_run: str, + epsilon: float, + sense: str, + model_kwargs: dict, + solve_kwargs: dict, + planning_horizons: str | None = None, +) -> tuple[str, str]: + """ + Run MGA (Modelling to Generate Alternatives) optimization. + + Parameters + ---------- + n : pypsa.Network + The PyPSA network instance (already contains cost-optimal solution) + mga_params : dict + MGA configuration from snakemake.params.mga + mga_run : str + Name of the MGA run to execute + epsilon : float + Slack parameter - how much the objective can deviate from optimal (e.g., 0.05 = 5%) + sense : str + Optimization sense for MGA objective: 'min' or 'max' + model_kwargs : dict + Arguments for model creation + solve_kwargs : dict + Arguments for solving + planning_horizons : str, optional + The current planning horizon year + + Returns + ------- + tuple[str, str] + (status, condition) from the optimization + """ + logger.info(f"Starting MGA optimization for run '{mga_run}'") + logger.info(f"MGA parameters: epsilon={epsilon}, sense={sense}") + + # Get run configuration + mga_runs = mga_params.get("runs", {}) + if mga_run not in mga_runs: + raise RuntimeError(f"MGA run '{mga_run}' not found in configuration") + + run_config = mga_runs[mga_run] + + # Build weights dictionary + weights = build_mga_weights(n, run_config) + + if not weights: + raise RuntimeError(f"No valid weights found for MGA run '{mga_run}'") + + logger.info(f"MGA weights: {weights}") + + # Store original objective for reference + # original_objective = n.objective + # logger.info(f"Original (cost-optimal) objective value: {original_objective}") + + # Call PyPSA's native MGA function + status, condition = n.optimize.optimize_mga( + snapshots=None, + multi_investment_periods=model_kwargs.get("multi_investment_periods", False), + weights=weights, + sense=sense, + slack=epsilon, + solver_name=solve_kwargs["solver_name"], + solver_options=solve_kwargs.get("solver_options", {}), + model_kwargs={ + key: val + for key, val in model_kwargs.items() + if key != "multi_investment_periods" + }, + assign_all_duals=solve_kwargs.get("assign_all_duals", False), + io_api=solve_kwargs.get("io_api", None), + keep_files=solve_kwargs.get("keep_files", False), + extra_functionality=partial( + extra_functionality, planning_horizons=planning_horizons + ), + ) + + # logger.info(f"MGA '{mga_run}' completed with status: {status}, condition: {condition}") + # logger.info(f"MGA objective value: {n.objective}") + + # # Store MGA metadata + # if not hasattr(n, "mga_results"): + # n.mga_results = {} + # n.mga_results[mga_run] = { + # "original_objective": original_objective, + # "mga_objective": n.objective, + # "status": status, + # "condition": condition, + # "epsilon": epsilon, + # "sense": sense, + # "weights": weights, + # } + + return status, condition + + def collect_kwargs( config: dict, solving: dict, @@ -1447,6 +1600,10 @@ def create_optimization_model( rolling_horizon = cf_solving.get("rolling_horizon", False) skip_iterations = cf_solving.get("skip_iterations", False) + # Check for MGA mode + mga_run = getattr(snakemake.wildcards, "mga_run", None) + is_mga_run = snakemake.params["mga"]["enable"] and mga_run is not None + if not n.lines.s_nom_extendable.any(): skip_iterations = True logger.info("No expandable lines found. Skipping iterative solving.") @@ -1459,7 +1616,31 @@ def create_optimization_model( with memory_logger( filename=getattr(snakemake.log, "memory", None), interval=logging_frequency ) as mem: - if rolling_horizon and snakemake.rule == "solve_operations_network": + if is_mga_run: + logger.info("Using MGA...") + model_kwargs, solve_kwargs = collect_kwargs( + snakemake.config, + snakemake.params.solving, + planning_horizons, + log_fn=snakemake.log.solver, + mode="single", + ) + + # Get epsilon and sense from wildcards + n.config = snakemake.config + n.params = snakemake.params + + status, condition = run_mga_optimization( + n, + mga_params=snakemake.params["mga"], + mga_run=mga_run, + epsilon=float(snakemake.wildcards.epsilon), + sense=snakemake.wildcards.sense, + model_kwargs=model_kwargs, + solve_kwargs=solve_kwargs, + planning_horizons=planning_horizons, + ) + elif rolling_horizon and snakemake.rule == "solve_operations_network": logger.info("Using rolling horizon optimization...") all_kwargs, _ = collect_kwargs( snakemake.config, @@ -1540,11 +1721,13 @@ def create_optimization_model( n.meta = dict(snakemake.config, **dict(wildcards=dict(snakemake.wildcards))) n.export_to_netcdf(snakemake.output.network) - with open(snakemake.output.config, "w") as file: - yaml.dump( - n.meta, - file, - default_flow_style=False, - allow_unicode=True, - sort_keys=False, - ) + # Only write config for non-MGA runs + if hasattr(snakemake.output, "config"): + with open(snakemake.output.config, "w") as file: + yaml.dump( + n.meta, + file, + default_flow_style=False, + allow_unicode=True, + sort_keys=False, + ) From e3eb62e6d746836510b7a1b2cb02cdbc147f0d34 Mon Sep 17 00:00:00 2001 From: Amos Schledorn Date: Mon, 12 Jan 2026 16:55:18 +0100 Subject: [PATCH 2/4] Remove commented-out code and logging from MGA optimization function --- scripts/solve_network.py | 20 -------------------- 1 file changed, 20 deletions(-) diff --git a/scripts/solve_network.py b/scripts/solve_network.py index 6464843032..5f1079f1da 100644 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -1378,10 +1378,6 @@ def run_mga_optimization( logger.info(f"MGA weights: {weights}") - # Store original objective for reference - # original_objective = n.objective - # logger.info(f"Original (cost-optimal) objective value: {original_objective}") - # Call PyPSA's native MGA function status, condition = n.optimize.optimize_mga( snapshots=None, @@ -1404,22 +1400,6 @@ def run_mga_optimization( ), ) - # logger.info(f"MGA '{mga_run}' completed with status: {status}, condition: {condition}") - # logger.info(f"MGA objective value: {n.objective}") - - # # Store MGA metadata - # if not hasattr(n, "mga_results"): - # n.mga_results = {} - # n.mga_results[mga_run] = { - # "original_objective": original_objective, - # "mga_objective": n.objective, - # "status": status, - # "condition": condition, - # "epsilon": epsilon, - # "sense": sense, - # "weights": weights, - # } - return status, condition From 6744f75f9ecfb6fe24ddbb50e1096fb046f4e33a Mon Sep 17 00:00:00 2001 From: Amos Schledorn Date: Wed, 25 Feb 2026 13:11:39 +0100 Subject: [PATCH 3/4] Use alternative solver options for MGA runs --- config/config.default.yaml | 1 + scripts/solve_network.py | 2 ++ 2 files changed, 3 insertions(+) diff --git a/config/config.default.yaml b/config/config.default.yaml index f79b140423..ed49f1630a 100644 --- a/config/config.default.yaml +++ b/config/config.default.yaml @@ -1275,6 +1275,7 @@ mga: enable: false epsilon: [0.05] sense: [min, max] + solver_options: gurobi-fallback runs: # docs in https://pypsa-eur.readthedocs.io/en/latest/configuration.html#data diff --git a/scripts/solve_network.py b/scripts/solve_network.py index 6300d3cbac..5b3bcd425c 100644 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -1609,6 +1609,8 @@ def create_optimization_model( ) as mem: if is_mga_run: logger.info("Using MGA...") + solving_params = snakemake.params.solving + solving_params["solver"]["options"] = snakemake.params.mga["solver_options"] model_kwargs, solve_kwargs = collect_kwargs( snakemake.config, snakemake.params.solving, From 141f4fd03d182c70ee113d2aebe15df970b5b0af Mon Sep 17 00:00:00 2001 From: "pypsa[bot]" <181215446+pypsa[bot]@users.noreply.github.com> Date: Tue, 14 Apr 2026 14:38:57 +0000 Subject: [PATCH 4/4] [pypsa-bot] run `generate-config` - `pixi run generate-config` --- config/config.default.yaml | 7 ------- 1 file changed, 7 deletions(-) diff --git a/config/config.default.yaml b/config/config.default.yaml index ed49f1630a..cdec450f07 100644 --- a/config/config.default.yaml +++ b/config/config.default.yaml @@ -1271,13 +1271,6 @@ solving: memory_logging_frequency: 5 runtime: 48h -mga: - enable: false - epsilon: [0.05] - sense: [min, max] - solver_options: gurobi-fallback - runs: - # docs in https://pypsa-eur.readthedocs.io/en/latest/configuration.html#data data: hotmaps_industrial_sites: