diff --git a/rules/supply_curves.smk b/rules/supply_curves.smk index 194711e..3ee7c8d 100644 --- a/rules/supply_curves.smk +++ b/rules/supply_curves.smk @@ -35,13 +35,32 @@ def _product_uses_renewables(product): return True +def _demand_levels_for_product(product): + """Demand levels to sweep when solving LCoX for a product. + + Renewable-fed routes (e.g. hbi/dri) get costlier at scale as cheaper + resource classes are exhausted, so the full sweep is needed to trace + the curve. Grid-connected routes (e.g. steel/eaf) draw on grid + electricity at a fixed, uncapped marginal cost in this stage, so their + LCoX is flat with respect to volume (verified: identical across all + demand levels for every region). Solving all levels for those products + is redundant; one point at the top of the range is enough and keeps + capacity headroom for the trade model. + """ + demand_levels = config.get("steel_demand_levels") + if _product_uses_renewables(product): + return demand_levels + return [max(demand_levels)] + + def _all_supply_curve_targets(): targets = [] + cost_year = config["trade_chains"].get("cost_year", 2050) for region in config["regions"]: wacc = config["trade_chains"].get("wacc", "uniform") for product in SUPPLY_CURVE_PRODUCTS: targets.append( - f"resources/supply_curves/cost_year~2050/wacc~{wacc}/{region}_marginal_cost_{product}.csv" + f"resources/supply_curves/cost_year~{cost_year}/wacc~{wacc}/{region}_marginal_cost_{product}.csv" ) return targets @@ -161,12 +180,12 @@ if config["enable"].get("run_supply_curve", True): scenario=config.get("supply_curve", {}).get( "default_scenario", "allocated_share" ), - product_demand_mt=config.get("steel_demand_levels"), + product_demand_mt=_demand_levels_for_product(wildcards.product), ), lco_unreserved=lambda wildcards: ( expand( f"resources/lco-{wildcards.product}/cost_year~{wildcards.cost_year}/wacc~{wildcards.wacc}/{wildcards.region}_unreserved/results_{{product_demand_mt}}.csv", - product_demand_mt=config.get("steel_demand_levels"), + product_demand_mt=_demand_levels_for_product(wildcards.product), ) if config.get("supply_curve", {}).get("generate_unreserved", False) and _product_uses_renewables(wildcards.product) diff --git a/workflow/notebooks/prepare-wacc.ipynb b/workflow/notebooks/prepare-wacc.ipynb index 8663705..3dff537 100644 --- a/workflow/notebooks/prepare-wacc.ipynb +++ b/workflow/notebooks/prepare-wacc.ipynb @@ -8,7 +8,6 @@ "outputs": [], "source": [ "import pandas as pd\n", - "import pycountry\n", "import wbdata" ] }, @@ -97,36 +96,6 @@ "regions = config['regions']" ] }, - { - "cell_type": "code", - "execution_count": null, - "id": "a532a77c", - "metadata": {}, - "outputs": [], - "source": [ - "country_name_corrections = {\n", - " \"Democratic Republic of the Congo\": \"Congo, The Democratic Republic of the\",\n", - " \"Republic of the Congo\": \"Republic of the Congo\",\n", - " \"Kosovo\": \"Republic of Kosovo\", # pycountry not supported\n", - " \"Russia\": \"Russian Federation\",\n", - " \"Turkey\": \"Türkiye\",\n", - " \"Venezuela\": \"Venezuela, Bolivarian Republic of\",\n", - " \"Tanzania\": \"United Republic of Tanzania\",\n", - " \"Bolivia\": \"Plurinational State of Bolivia\",\n", - " \"Vietnam\": \"Viet Nam\",\n", - " \"South Korea\": \"Korea, Republic of\",\n", - " \"North Korea\": \"Korea, Democratic People's Republic of\",\n", - " \"Taiwan\": \"Taiwan, Province of China\",\n", - " \"Laos\": \"Lao People's Democratic Republic\",\n", - " \"Brunei\": \"Brunei Darussalam\",\n", - " \"Equatorial French Guiana\": \"French Guiana\",\n", - " \"Syria\": \"Syrian Arab Republic\", \n", - " \"Palestine\": \"Palestine, State of\",\n", - " \"Moldova\": \"Republic of Moldova\",\n", - " \"North Korea\": \"Korea, Democratic People's Republic of\",\n", - "}" - ] - }, { "cell_type": "code", "execution_count": null, @@ -134,13 +103,13 @@ "metadata": {}, "outputs": [], "source": [ - "# Build a mapping from country name to ISO-3 code (still needed for GDP World Bank lookup)\n", - "country_name_to_iso = {}\n", - "for country in pycountry.countries:\n", - " country_name_to_iso[country.name] = country.alpha_3\n", - " # Add common names\n", - " if hasattr(country, 'official_name'):\n", - " country_name_to_iso[country.official_name] = country.alpha_3\n", + "# Build a mapping from World Bank's own country name to ISO-3 code.\n", + "# Using World Bank's own metadata (rather than pycountry) avoids name-matching\n", + "# mismatches for countries whose WB display name differs from pycountry's\n", + "# (e.g. WB \"Turkiye\" vs pycountry \"Türkiye\", WB \"Korea, Rep.\" vs pycountry\n", + "# \"Korea, Republic of\"), which previously caused those countries' GDP rows to be\n", + "# dropped and, for single-country regions, produced NaN GDP-weighted WACC.\n", + "country_name_to_iso = {c[\"name\"]: c[\"id\"] for c in wbdata.get_countries()}\n", "\n", "# Build a mapping from ISO-3 code to region (config regions now use ISO-3 codes directly)\n", "iso_to_region = {}\n", diff --git a/workflow/scripts/calculate_lcox.py b/workflow/scripts/calculate_lcox.py index 47f4eb2..ec55f81 100644 --- a/workflow/scripts/calculate_lcox.py +++ b/workflow/scripts/calculate_lcox.py @@ -408,7 +408,7 @@ def solve_network(network, config): logger.error(f"Solver exception: {e}") raise - return network + return network, status_ok # ============================================================================ @@ -416,7 +416,7 @@ def solve_network(network, config): # ============================================================================ -def extract_lcox(network, product, demands): +def extract_lcox(network, product, demands, status_ok=True): """Extract LCOX from an optimized network. Parameters @@ -428,6 +428,13 @@ def extract_lcox(network, product, demands): demands : dict Must include key `'product_demand_mt'` (float, Mt/year) used to compute annual production and per-unit LCOX. + status_ok : bool + Whether the solver reported a genuinely optimal termination status. + A non-optimal solve (e.g. numerical trouble during the barrier + method) can leave a stale but finite `network.objective` behind even + though the result isn't trustworthy; checking only `None`/`NaN` lets + that garbage value through as if it were a valid LCOX. Requiring + `status_ok` closes that gap. Returns ------- @@ -455,6 +462,9 @@ def extract_lcox(network, product, demands): ) try: + if not status_ok: + raise ValueError("Optimization did not report an optimal status") + obj_value = network.objective if obj_value is None or np.isnan(obj_value): raise ValueError("Optimization failed to return valid objective") @@ -598,16 +608,19 @@ def extract_lcox(network, product, demands): if snakemake.config.get("debug_network_inspection", False): inspect_network(network, product) # Debug inspection try: - solve_network(network, snakemake.config) + network, status_ok = solve_network(network, snakemake.config) optimization_status = ( "optimal" - if network.objective is not None and not np.isnan(network.objective) + if status_ok + and network.objective is not None + and not np.isnan(network.objective) else "infeasible" ) except Exception as e: logger.warning( f"Solver error for product demand {product_demand_mt} Mt/year: {e}" ) + status_ok = False optimization_status = "error" if optimization_status != "optimal": @@ -620,6 +633,7 @@ def extract_lcox(network, product, demands): results_df = extract_lcox( network=network, product=product, + status_ok=status_ok, demands=scaled_demands, )