diff --git a/rules/collect.smk b/rules/collect.smk index 743407ee92..4ca6d06505 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: @@ -76,6 +78,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: message: "Collecting solved sector-coupled network files" @@ -88,6 +103,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: message: "Collecting solved sector-coupled network files with perfect foresight" diff --git a/rules/solve_electricity.smk b/rules/solve_electricity.smk index 1259943104..07d186c906 100644 --- a/rules/solve_electricity.smk +++ b/rules/solve_electricity.smk @@ -13,6 +13,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: @@ -41,6 +42,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: message: "Solving electricity network operations optimization for {wildcards.clusters} clusters and {wildcards.opts} electric options" diff --git a/rules/solve_overnight.smk b/rules/solve_overnight.smk index a56dfb2809..99620e7061 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: message: @@ -13,6 +15,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" @@ -48,3 +51,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 0fa8e90fc9..5b3bcd425c 100644 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -1280,6 +1280,139 @@ 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}") + + # 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 + ), + ) + + return status, condition + + def collect_kwargs( config: dict, solving: dict, @@ -1458,6 +1591,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.") @@ -1470,7 +1607,33 @@ def create_optimization_model( with memory_logger( filename=getattr(snakemake.log, "memory", None), interval=logging_frequency ) as mem: - if rolling_horizon: + 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, + 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: logger.info("Using rolling horizon optimization...") all_kwargs, _ = collect_kwargs( snakemake.config, @@ -1554,11 +1717,13 @@ def create_optimization_model( if snakemake.output.get("model"): n.model.to_netcdf(snakemake.output.model) - 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, + )