Merge branch 'pik-piam:main' into main

Author: irr-github

config/default_config.yaml

@@ -1,10 +1,10 @@
-# SPDX-FileCopyrightText: : 2022 The PyPSA-China Authors
+# SPDX-FileCopyrightText: : 2025 The PyPSA-China-PIK Authors
 #
 # SPDX-License-Identifier: CC0-1.0
 # DO NOT FORGET TO LOOK AT THE TECHNOLOGY CONFIGURATION FILES
 
 run:
-  name: "unamed_default_run"
+  name: "unnamed_run"
 foresight: "overnight"
 
 paths:
@@ -73,7 +73,7 @@ enable:
   retrieve_raster: false
 
 atlite:
-  cutout_name: "China-2020c"
+  cutout_name: "China-2020c_bare"
   freq: "h"
   nprocesses: 1
   show_progress: true
@@ -169,16 +169,23 @@ bus_carrier: {
 Techs:
   ## tech selection
   #TODO: rework vre_techs to split with clean_techs. Merge vre with non_dispatchable
-  vre_techs: [ "onwind","offwind","solar","solar thermal","hydroelectricity", "nuclear","biomass","beccs","heat pump","resistive heater","Sabatier","H2 CHP", "fuel cell"]
+  vre_techs: [ "onwind","offwind","solar","solar thermal","hydroelectricity", "nuclear","biomass","beccs","heat pump","resistive heater","Sabatier", "fuel cell", "H2 CHP"]
   non_dispatchable: ['Offshore Wind', 'Onshore Wind', 'Solar', 'Solar Residential', "Solar Thermal"]
-  conv_techs: ["OCGT", "CCGT", "CHP gas", "gas boiler","coal boiler","coal power plant","CHP coal"]
+  conv_techs: ["CCGT", "CHP gas", "CHP OCGT gas", "gas boiler","coal boiler","coal power plant","CHP coal", "OCGT"]
   store_techs: ["H2","battery","water tanks","PHS"]
   non_dispatchable: ["onwind", "offwind", "solar", "solar thermal"]
   hydrogen_lines: true
-  coal_ccs_retrofit: true # currently myopic pathway only.  CC = co2 cap 
+  coal_ccs_retrofit: false # currently myopic pathway only.  CC = co2 cap 
+
+## Nuclear capacity growth limit
+nuclear_reactors:
+  enable_growth_limit: true  # Switch: true=enable growth limit, false=unrestricted growth
+  max_annual_capacity_addition: 5000  # Fixed MW to add per year
+  base_year: 2020  # Base year for capacity reference (should be <= first planning year)
+  # base_capacity: 50000  # (optional) Base year total capacity in MW. If not set, auto-detect from base_year network
 
 ## add components (overwrites vre tech choice)
-heat_coupling: true
+heat_coupling: false
 add_biomass: True
 add_hydro: True
 add_H2: True
@@ -224,6 +231,12 @@ solving:
     tunnel_port: 27000
     timeout_s: 15
     login_node: "03"
+    gurobi_home: "/p/projects/rd3mod/gurobi1103/linux64"
+    license_path: "/p/projects/rd3mod/gurobi_rc/gurobi.lic"
+    verbose: "1"
+    timeout: "10"
+    ssl_cert: "/p/projects/rd3mod/ssl/ca-bundle.pem_2022-02-08"
+    grb_cafile: "/p/projects/rd3mod/ssl/ca-bundle.pem_2022-02-08"
 
   solver_options:
     highs-default:
@@ -294,13 +307,13 @@ lines:
     base_year: 2020 # base step year, the exp limit is relative to this year
 
 security:
-  line_margin: 70 # max percent of line capacity
+  line_security_margin: 70 # max percent of line capacity
 
 # brownfield options
 existing_capacities:
   add: True # whether to add brownfield capacities
   grouping_years: [1985, 1990, 1995, 2000, 2005, 2010, 2015, 2020, 2025, 2030, 2035, 2040, 2045, 2050, 2055, 2060]
-  collapse_years: False # Treat as a single unit when preparing & solving network 
+  collapse_years: True # Treat as a single unit when preparing & solving network 
   threshold_capacity: 80 # MW
   techs: ['coal','CHP coal', 'CHP gas', 'OCGT', 'CCGT', 'solar', 'onwind', 'offwind', 'nuclear', "PHS"]
   node_assignment: simple # simple | gps
@@ -346,3 +359,13 @@ fetch_regions:
   simplify_tol: 
     eez: 0.5
     land: 0.05
+
+subsidies:
+  # New format: subsidies.fuel_type -> province -> value (EUR/MWh) only allowed for negative values (subsidies)
+  gas:
+    Guangdong: -8
+    Jiangsu:  -8.5
+    Zhejiang: -9
+    Beijing: -10  
+    Tianjin: -10
+    Shanghai: -10

docs/.nav.yml

@@ -20,6 +20,6 @@ nav:
   - Tutorials: tutorials
   - Reference: reference
 
-append_unmatched: true 
+append_unmatched: false 
 
 

workflow/notebooks/land_availability.ipynb

@@ -496,17 +496,26 @@ def setup_gurobi_tunnel_and_env(
         logger.error("SSH tunnel communication timed out.")
 
     os.environ["https_proxy"] = f"socks5://127.0.0.1:{port}"
-    os.environ["SSL_CERT_FILE"] = "/p/projects/rd3mod/ssl/ca-bundle.pem_2022-02-08"
-    os.environ["GRB_CAFILE"] = "/p/projects/rd3mod/ssl/ca-bundle.pem_2022-02-08"
+    os.environ["SSL_CERT_FILE"] = tunnel_config.get(
+        "ssl_cert", "/p/projects/rd3mod/ssl/ca-bundle.pem_2022-02-08"
+    )
+    os.environ["GRB_CAFILE"] = tunnel_config.get(
+        "grb_cafile", "/p/projects/rd3mod/ssl/ca-bundle.pem_2022-02-08"
+    )
 
     # Set up Gurobi environment variables
-    os.environ["GUROBI_HOME"] = "/p/projects/rd3mod/gurobi1103/linux64"
+    # TODO soft code
+    os.environ["GUROBI_HOME"] = tunnel_config.get(
+        "gurobi_home", "/p/projects/rd3mod/gurobi1103/linux64"
+    )
     os.environ["PATH"] += f":{os.environ['GUROBI_HOME']}/bin"
     if "LD_LIBRARY_PATH" in os.environ:
         os.environ["LD_LIBRARY_PATH"] += f":{os.environ['GUROBI_HOME']}/lib"
-    os.environ["GRB_LICENSE_FILE"] = "/p/projects/rd3mod/gurobi_rc/gurobi.lic"
-    os.environ["GRB_CURLVERBOSE"] = "1"
-    os.environ["GRB_SERVER_TIMEOUT"] = "10"
+    os.environ["GRB_LICENSE_FILE"] = tunnel_config.get(
+        "license_path", "/p/projects/rd3mod/gurobi_rc/gurobi.lic"
+    )
+    os.environ["GRB_CURLVERBOSE"] = tunnel_config.get("verbose", "1")
+    os.environ["GRB_SERVER_TIMEOUT"] = tunnel_config.get("timeout", "10")
 
     return socks_proc
 

workflow/scripts/_helpers.py

@@ -246,6 +246,86 @@ def add_co2_capture_support(
     )
 
 
+def add_fuel_subsidies(n: pypsa.Network, subsidy_config: dict):
+    """Apply fuel subsidies to generators as a post-processing step.
+    
+    Subsidies are applied to generators based on their carrier and location.
+    The subsidy values (in EUR/MWh fuel) are divided by efficiency to convert
+    to electricity basis (EUR/MWhel). Links are not modified as they get their
+    fuel from generators.
+    
+    Args:
+        n (pypsa.Network): The network object to modify.
+        subsidy_config (dict): Subsidy configuration dictionary with keys like
+            "coal" or "gas", each containing a dict mapping provinces to subsidy values.
+    """
+    if not subsidy_config:
+        return
+    
+    carriers = subsidy_config.keys()
+    
+    for carrier in carriers:
+        subs_dict = subsidy_config.get(carrier, {})
+        if not subs_dict:
+            continue
+        
+        # Convert subsidy dict to Series indexed by province
+        subs = pd.Series(subs_dict, dtype=float)
+        
+        # Check that subsidies are non-positive (negative = subsidy, positive would be a reward)
+        if (subs > 0).any():
+            raise ValueError(
+                f"Positive subsidy values found for carrier '{carrier}': "
+                f"{subs[subs > 0].to_dict()}. Only zero or negative values are allowed "
+                f"(negative reduces marginal cost, positive would increase it)."
+            )
+        
+        # Check if location column exists
+        if "location" not in n.generators.columns:
+            logger.warning(
+                f"Location column not found in generators. "
+                f"Cannot apply subsidies for carrier '{carrier}'."
+            )
+            continue
+        
+        # Query generators with matching carrier and location in subsidy provinces
+        mask = n.generators.query(
+            f"carrier == @carrier and location in @subs.index"
+        ).index
+        
+        if mask.empty:
+            logger.warning(
+                f"No generators found with carrier '{carrier}' and locations "
+                f"in {list(subs.index)}. Skipping subsidy application."
+            )
+            continue
+        
+        # Merge subsidies with generators by location
+        gen_locs = n.generators.loc[mask, "location"]
+        subs_to_apply = gen_locs.map(subs).fillna(0.0)
+        
+        # Check if all provinces were found
+        missing_provs = set(subs.index) - set(gen_locs.unique())
+        if missing_provs:
+            logger.warning(
+                f"Subsidies specified for provinces {missing_provs} but no "
+                f"generators found with carrier '{carrier}' in these provinces."
+            )
+        
+        # Apply subsidies: divide by efficiency to convert from fuel to electricity basis
+        # Handle cases where efficiency might be NaN or missing
+        efficiencies = n.generators.loc[mask, "efficiency"].fillna(1.0)
+        subs_electricity = subs_to_apply / efficiencies
+        
+        # Subtract subsidy from marginal cost (negative subsidy = cost reduction)
+        n.generators.loc[mask, "marginal_cost"] -= subs_electricity
+        
+        logger.info(
+            f"Applied subsidies for carrier '{carrier}' to {len(mask)} generators "
+            f"in provinces {sorted(gen_locs.unique())}"
+        )
+
+
 def add_conventional_generators(
     network: pypsa.Network,
     nodes: pd.Index,
@@ -660,6 +740,7 @@ def add_voltage_links(network: pypsa.Network, config: dict):
             p_nom=edges["p_nom"].values,
             p_nom_min=edges["p_nom"].values,
             p_min_pu=0,
+            p_max_pu=line_margin,
             efficiency=config["transmission_efficiency"]["DC"]["efficiency_static"]
             * config["transmission_efficiency"]["DC"]["efficiency_per_1000km"] ** (lengths / 1000),
             length=0,
@@ -1640,6 +1721,12 @@ def prepare_network(
         add_voltage_links(network, config)
 
     assign_locations(network)
+    
+    # Apply fuel subsidies as post-processing step
+    subsidy_config = config.get("subsidies", {})
+    if subsidy_config:
+        add_fuel_subsidies(network, subsidy_config)
+    
     return network
 
 

workflow/scripts/prepare_network.py

@@ -29,6 +29,70 @@
 logger = logging.getLogger(__name__)
 
 
+def calc_nuclear_expansion_limit(
+    n: pypsa.Network,
+    config: dict,
+    planning_year: int,
+    network_path: str,
+) -> None:
+    """
+    Calculate and apply the nuclear expansion limit from configuration.
+    
+    Args:
+        n (pypsa.Network): the network object
+        config (dict): full configuration dictionary (mutated in place)
+        planning_year (int): target planning horizon year
+        network_path (str): path to the current network file, used to locate base year
+    """
+    nuclear_cfg = config.setdefault("nuclear_reactors", {})
+    if not nuclear_cfg.get("enable_growth_limit"):
+        return
+    
+    annual_addition = nuclear_cfg.get("max_annual_capacity_addition")
+    if not annual_addition:
+        logger.warning("Nuclear growth limit enabled but max_annual_capacity_addition missing")
+        return
+    
+    base_year = nuclear_cfg.get("base_year", 2020)
+    n_years = planning_year - base_year
+    if n_years <= 0:
+        logger.info(
+            "Planning year %s is not after base year %s; skipping nuclear expansion limit",
+            planning_year,
+            base_year,
+        )
+        return
+    
+    base_capacity = nuclear_cfg.get("base_capacity")
+    if base_capacity is None:
+        base_path = network_path.replace(f"ntwk_{planning_year}.nc", f"ntwk_{base_year}.nc")
+        if os.path.exists(base_path):
+            n_base = pypsa.Network(base_path)
+            base_capacity = n_base.generators[n_base.generators.carrier == "nuclear"]["p_nom"].sum()
+        else:
+            base_capacity = n.generators[n.generators.carrier == "nuclear"]["p_nom"].sum()
+    
+    max_capacity = base_capacity + annual_addition * n_years
+    logger.info(
+        f"Adding nuclear expansion limit for {planning_year}: {max_capacity:.0f} MW "
+        f"[{base_capacity:.0f} + {annual_addition:.0f} × {n_years} years]"
+    )
+
+    nuclear_gens_ext = n.generators[
+        (n.generators.carrier == "nuclear") & (n.generators.p_nom_extendable == True)
+    ].index
+    
+    if len(nuclear_gens_ext) == 0:
+        logger.warning("No extendable nuclear generators found")
+        return
+    
+    n.generators.loc[nuclear_gens_ext, "p_nom_max"] = max_capacity
+    nuclear_cfg["expansion_limit"] = max_capacity
+    logger.info(
+        f"Nuclear expansion limit set: {max_capacity:.0f} MW for {len(nuclear_gens_ext)} generators"
+    )
+
+
 def set_transmission_limit(n: pypsa.Network, kind: str, factor: float, n_years=1):
     """
     Set global transimission limit constraints - adapted from pypsa-eur
@@ -264,6 +328,30 @@ def prepare_network(
     return n
 
 
+def add_nuclear_expansion_constraints(n: pypsa.Network):
+    """
+    Add nuclear expansion limit constraint if configured.
+    
+    Args:
+        n (pypsa.Network): the network object
+    """
+    limit = n.config.get("nuclear_reactors", {}).get("expansion_limit")
+    if limit is None:
+        return
+    
+    nuclear_gens_ext = n.generators[
+        (n.generators.carrier == "nuclear") & (n.generators.p_nom_extendable == True)
+    ].index
+    
+    if len(nuclear_gens_ext) == 0:
+        return
+    
+    # Add global constraint: sum of all nuclear p_nom <= limit
+    lhs = n.model["Generator-p_nom"].loc[nuclear_gens_ext].sum()
+    n.model.add_constraints(lhs <= limit, name="nuclear_expansion_limit")
+    logger.info(f"Applied global nuclear constraint: sum(p_nom) <= {limit:.0f} MW")
+
+
 def add_battery_constraints(n: pypsa.Network):
     """
     Add constraint ensuring that charger = discharger, i.e.
@@ -743,6 +831,8 @@ def extra_functionality(n: pypsa.Network, _) -> None:
     config = n.config
     add_battery_constraints(n)
     add_transmission_constraints(n)
+    add_nuclear_expansion_constraints(n)
+    
     if config["heat_coupling"]:
         add_water_tank_charger_constraints(n, config)
         add_chp_constraints(n)
@@ -831,6 +921,7 @@ def solve_network(
             configfiles="resources/tmp/pseudo-coupled.yaml",
         )
     configure_logging(snakemake)
+    config = snakemake.config
 
     opts = snakemake.wildcards.get("opts", "")
     if "sector_opts" in snakemake.wildcards.keys():
@@ -868,6 +959,13 @@ def solve_network(
     # # TODO: remove ugly hack
     # n.storage_units.p_nom_max = n.storage_units.p_nom * 1.05**exp_years
 
+    calc_nuclear_expansion_limit(
+        n=n,
+        config=config,
+        planning_year=int(snakemake.wildcards.planning_horizons),
+        network_path=snakemake.input.network_name,
+    )
+
     if tunnel:
         logger.info(f"tunnel process alive? {tunnel.poll()}")
 
@@ -880,7 +978,7 @@ def solve_network(
 
         n = solve_network(
             n,
-            config=snakemake.config,
+            config=config,
             solving=snakemake.params.solving,
             opts=opts,
             log_fn=snakemake.log.solver,