feat: add mga support

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"

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"

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"

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,
+            )