feat: Improve numerical conditioning

CHANGELOG.md

@@ -50,6 +50,24 @@ introduce breaking changes to configuration and outputs.
 
 ### Changed
 
+- Reformulated the **L1 deviation penalties** (production, animal-feed, diet
+  stability) from an absolute-value auxiliary variable with two inequality
+  rows per link to an equivalent equality split into non-negative
+  positive/negative deviation parts, and priced the zero-baseline
+  land-conversion penalty directly on link flows. Together with a faster
+  nodal-balance construction in the vendored PyPSA fork, this cuts
+  full-resolution solve times by roughly a third (about 40% fewer
+  constraint rows after presolve) with identical optima up to solver
+  tolerance.
+- Improved the optimisation model's **numerical conditioning** to remove
+  Gurobi's "large matrix coefficient range" warning. The CH₄ and N₂O emission
+  buses are now denominated in kilotonnes (previously tonnes) so their flow
+  coefficients sit within a few orders of the CO₂ bus, and a new `numerics`
+  config block clips physically-negligible coefficients at build time
+  (sub-hectare areas, trace irrigation/carbon fluxes, rounding-level cost
+  corrections). The former `land.filtering` thresholds now live under
+  `numerics`. Emission totals and the objective are unchanged (to within
+  solver tolerance); only reported CH₄/N₂O bus flows change units.
 - The baseline diet is now derived from **FAOSTAT Food Balance Sheets** by
   default (`diet.source: fbs`), computed from per-country food supply energy
   at model-basis densities and corrected for consumer waste. The GDD-IA

config/default.yaml

@@ -76,6 +76,26 @@ netcdf:
     zlib: true
     complevel: 4
 
+# --- section: numerics ---
+# Build-time numerical conditioning. A few upstream datasets carry values
+# many orders of magnitude below anything that affects results (sub-hectare
+# disaggregated areas, trace yields/irrigation/carbon fluxes, rounding-level
+# cost corrections); leaving them in widens the solver's coefficient/bounds/RHS
+# ranges to ~20 orders and triggers numerical-scaling warnings. Every threshold
+# is deliberately conservative -- far below any value that changes the solution.
+# Set any to 0 to disable that clip.
+numerics:
+  # Structural pruning: drop links/rows whose defining input is implausibly
+  # small (these components would be inert or agronomically non-viable).
+  min_land_area_ha: 100                   # drop resource classes below this area
+  min_crop_yield_t_per_ha: 0.01           # drop crop links below this yield (~1% of entries)
+  min_grassland_yield_t_per_ha: 0.05      # drop grassland links below this yield (~6% of entries)
+  # Coefficient clipping: zero negligible coefficients on retained components.
+  min_link_area_mha: 0.000001            # ~1 ha; smaller land areas -> 0
+  min_water_requirement_m3_per_ha: 0.1   # <0.1 mm/ha irrigation -> no water draw
+  min_co2_coefficient_tco2_per_ha: 0.001 # near-zero land carbon flux -> 0
+  min_cost_correction_bnusd: 0.000001    # ~0.001 USD per unit flow -> 0
+
 # --- section: validation ---
 validation:
   use_actual_yields: true
@@ -227,10 +247,6 @@ land:
   conversion_cost_nonforest_usd_per_ha: 2000 # Overnight investment cost for converting non-forested land to agriculture (2024 USD/ha); sources in docs/land_use.rst
   investment_horizon: 30 # Years over which to annualize land conversion investment costs (harmonized with luc.horizon_years and the ~30-year Cook-Patton regrowth window)
   discount_rate: 0.05 # Annual discount rate for annualizing land conversion investment costs
-  filtering:
-    min_crop_yield_t_per_ha: 0.01      # Minimum yield for crop links (t/ha); filters ~1% of entries
-    min_grassland_yield_t_per_ha: 0.05 # Minimum yield for grassland links (t/ha); filters ~6% of entries
-    min_area_ha: 100                    # Minimum land area (ha); filters very small resource classes
 
 # --- section: water ---
 water:

config/schemas/config.schema.yaml

@@ -16,6 +16,7 @@ required:
   - paths
   - calibration
   - netcdf
+  - numerics
   - validation
   - consumer_values
   - food_incentives
@@ -160,6 +161,48 @@ properties:
           - type: "null"
         description: "Compression settings passed to xarray.Dataset.to_netcdf; null disables compression"
 
+  numerics:
+    type: object
+    description: "Build-time numerical-conditioning thresholds (structural pruning + coefficient clipping)"
+    required:
+      - min_land_area_ha
+      - min_crop_yield_t_per_ha
+      - min_grassland_yield_t_per_ha
+      - min_link_area_mha
+      - min_water_requirement_m3_per_ha
+      - min_co2_coefficient_tco2_per_ha
+      - min_cost_correction_bnusd
+    additionalProperties: false
+    properties:
+      min_land_area_ha:
+        type: number
+        minimum: 0
+        description: "Resource classes below this land area (ha) are dropped"
+      min_crop_yield_t_per_ha:
+        type: number
+        minimum: 0
+        description: "Crop production links below this yield (t/ha) are dropped"
+      min_grassland_yield_t_per_ha:
+        type: number
+        minimum: 0
+        description: "Grassland feed links below this yield (t/ha) are dropped"
+      min_link_area_mha:
+        type: number
+        minimum: 0
+        description: "Land areas (p_nom_max, baseline) below this many Mha are zeroed"
+      min_water_requirement_m3_per_ha:
+        type: number
+        minimum: 0
+        description: "Irrigation water coefficients below this many m3/ha are zeroed"
+      min_co2_coefficient_tco2_per_ha:
+        type: number
+        minimum: 0
+        description: "Land carbon-flux coefficients below this many tCO2/ha are zeroed"
+      min_cost_correction_bnusd:
+        type: number
+        minimum: 0
+        description: "Cost-calibration corrections below this many bnUSD per unit flow are zeroed"
+
   validation:
     type: object
     required:
@@ -306,7 +349,7 @@ properties:
 
   land:
     type: object
-    required: [regional_limit, reforestation_cap, land_use_cost_usd_per_ha, conversion_cost_forest_usd_per_ha, conversion_cost_nonforest_usd_per_ha, investment_horizon, discount_rate, filtering]
+    required: [regional_limit, reforestation_cap, land_use_cost_usd_per_ha, conversion_cost_forest_usd_per_ha, conversion_cost_nonforest_usd_per_ha, investment_horizon, discount_rate]
     additionalProperties: false
     properties:
       regional_limit:
@@ -361,23 +404,6 @@ properties:
         minimum: 0
         maximum: 1
         description: "Annual discount rate for annualizing land conversion investment costs"
-      filtering:
-        type: object
-        required: [min_crop_yield_t_per_ha, min_grassland_yield_t_per_ha, min_area_ha]
-        additionalProperties: false
-        properties:
-          min_crop_yield_t_per_ha:
-            type: number
-            minimum: 0
-            description: "Minimum yield threshold for crop production links (t/ha)"
-          min_grassland_yield_t_per_ha:
-            type: number
-            minimum: 0
-            description: "Minimum yield threshold for grassland feed links (t/ha)"
-          min_area_ha:
-            type: number
-            minimum: 0
-            description: "Minimum land area threshold (ha)"
 
   water:
     type: object

data/curated/calibration/default/provenance.yaml

@@ -4,8 +4,8 @@
 # tools/calibrate; do not edit by hand.
 # Licensing: see the annotation in REUSE.toml.
 base_config: config/default.yaml
-generated_at: '2026-07-02T21:35:52+00:00'
-git_commit: 08a68eccbce283e7a4b01940e002fd266f000cef
+generated_at: '2026-07-03T17:40:28+00:00'
+git_commit: 22281c512133c39d1fd438f77f076a0b931c3d33
 source: default
 structural_config:
   aggregation.irrigated_area_source: current
@@ -801,9 +801,6 @@ structural_config:
   land.conversion_cost_forest_usd_per_ha: 8000
   land.conversion_cost_nonforest_usd_per_ha: 2000
   land.discount_rate: 0.05
-  land.filtering.min_area_ha: 100
-  land.filtering.min_crop_yield_t_per_ha: 0.01
-  land.filtering.min_grassland_yield_t_per_ha: 0.05
   land.investment_horizon: 30
   land.land_use_cost_usd_per_ha: 0.0
   luc.cropland_source: gaez
@@ -833,6 +830,13 @@ structural_config:
   - alfalfa
   - silage-maize
   - biomass-sorghum
+  numerics.min_co2_coefficient_tco2_per_ha: 0.001
+  numerics.min_cost_correction_bnusd: 1.0e-06
+  numerics.min_crop_yield_t_per_ha: 0.01
+  numerics.min_grassland_yield_t_per_ha: 0.05
+  numerics.min_land_area_ha: 100
+  numerics.min_link_area_mha: 1.0e-06
+  numerics.min_water_requirement_m3_per_ha: 0.1
   optimal_taxes.enabled: false
   residues.max_feed_fraction: 0.3
   residues.max_feed_fraction_by_region.Asia: 0.7

data/curated/calibration/gbd-anchored/provenance.yaml

@@ -4,8 +4,8 @@
 # tools/calibrate; do not edit by hand.
 # Licensing: see the annotation in REUSE.toml.
 base_config: config/gsa.yaml
-generated_at: '2026-07-02T22:28:37+00:00'
-git_commit: 08a68eccbce283e7a4b01940e002fd266f000cef
+generated_at: '2026-07-03T17:40:31+00:00'
+git_commit: 22281c512133c39d1fd438f77f076a0b931c3d33
 source: gbd-anchored
 structural_config:
   aggregation.irrigated_area_source: current
@@ -799,9 +799,6 @@ structural_config:
   land.conversion_cost_forest_usd_per_ha: 8000
   land.conversion_cost_nonforest_usd_per_ha: 2000
   land.discount_rate: 0.05
-  land.filtering.min_area_ha: 100
-  land.filtering.min_crop_yield_t_per_ha: 0.01
-  land.filtering.min_grassland_yield_t_per_ha: 0.05
   land.investment_horizon: 30
   land.land_use_cost_usd_per_ha: 0.0
   luc.cropland_source: gaez
@@ -831,6 +828,13 @@ structural_config:
   - alfalfa
   - silage-maize
   - biomass-sorghum
+  numerics.min_co2_coefficient_tco2_per_ha: 0.001
+  numerics.min_cost_correction_bnusd: 1.0e-06
+  numerics.min_crop_yield_t_per_ha: 0.01
+  numerics.min_grassland_yield_t_per_ha: 0.05
+  numerics.min_land_area_ha: 100
+  numerics.min_link_area_mha: 1.0e-06
+  numerics.min_water_requirement_m3_per_ha: 0.1
   optimal_taxes.enabled: false
   residues.max_feed_fraction: 0.3
   residues.max_feed_fraction_by_region.Asia: 0.7

data/curated/calibration/gdd-ia/provenance.yaml

@@ -4,8 +4,8 @@
 # tools/calibrate; do not edit by hand.
 # Licensing: see the annotation in REUSE.toml.
 base_config: config/central.yaml
-generated_at: '2026-07-02T21:56:54+00:00'
-git_commit: 08a68eccbce283e7a4b01940e002fd266f000cef
+generated_at: '2026-07-03T17:40:30+00:00'
+git_commit: 22281c512133c39d1fd438f77f076a0b931c3d33
 source: gdd-ia
 structural_config:
   aggregation.irrigated_area_source: current
@@ -799,9 +799,6 @@ structural_config:
   land.conversion_cost_forest_usd_per_ha: 8000
   land.conversion_cost_nonforest_usd_per_ha: 2000
   land.discount_rate: 0.05
-  land.filtering.min_area_ha: 100
-  land.filtering.min_crop_yield_t_per_ha: 0.01
-  land.filtering.min_grassland_yield_t_per_ha: 0.05
   land.investment_horizon: 30
   land.land_use_cost_usd_per_ha: 0.0
   luc.cropland_source: gaez
@@ -831,6 +828,13 @@ structural_config:
   - alfalfa
   - silage-maize
   - biomass-sorghum
+  numerics.min_co2_coefficient_tco2_per_ha: 0.001
+  numerics.min_cost_correction_bnusd: 1.0e-06
+  numerics.min_crop_yield_t_per_ha: 0.01
+  numerics.min_grassland_yield_t_per_ha: 0.05
+  numerics.min_land_area_ha: 100
+  numerics.min_link_area_mha: 1.0e-06
+  numerics.min_water_requirement_m3_per_ha: 0.1
   optimal_taxes.enabled: false
   residues.max_feed_fraction: 0.3
   residues.max_feed_fraction_by_region.Asia: 0.7

docs/model_framework.rst

@@ -177,9 +177,9 @@ The model uses consistent units throughout:
   * Nutritional energy (calories): kcal/person/day → PJ/year
 
 **Emissions**
-  * Methane (CH₄) and nitrous oxide (N₂O): tonnes of gas, converted downstream to Mt CO₂-eq
-  * Aggregate greenhouse gases: tCO₂-eq (tonnes CO₂-equivalent)
-  * Conversion factors: CH₄ (28 GWP100), N₂O (265 GWP100)
+  * Methane (CH₄) and nitrous oxide (N₂O): kilotonnes of gas (kept close in magnitude to the CO₂ bus for numerical conditioning), converted downstream to Mt CO₂-eq
+  * Carbon dioxide (CO₂) and aggregate greenhouse gases: Mt CO₂-eq
+  * Conversion factors: CH₄ (27 GWP100), N₂O (273 GWP100)
 
 **Water**
   * Water use: Mm³ (million cubic meters)

pixi.toml

@@ -66,7 +66,7 @@ glade-workflow = { path = ".", editable = true }
 snakemake-logger-plugin-compact = { path = "tools/snakemake-logger-plugin-compact" }
 
 # Optimization and modeling; own branches of pypsa and linopy with custom features needed by GLADE
-pypsa = { git = "https://github.com/koen-vg/PyPSA.git", tag = "v1.2.0+glade2" }
+pypsa = { git = "https://github.com/koen-vg/PyPSA.git", tag = "v1.2.0+glade3" }
 linopy = { git = "https://github.com/koen-vg/linopy.git", tag = "v0.8.0+glade2" }
 
 # Geospatial processing

workflow/rules/model.smk

@@ -217,6 +217,7 @@ rule build_model:
         baseline_year=config["baseline_year"],
         validation=config["validation"],
         deviation_penalty=config["deviation_penalty"],
+        numerics=config["numerics"],
         netcdf=config["netcdf"],
         # Add health-cluster stores when health is enabled in the base config or
         # any scenario (the build is shared across scenarios).

workflow/scripts/analysis/extract_ghg_attribution.py

@@ -44,7 +44,7 @@
 import pandas as pd
 import pypsa
 
-from workflow.scripts.constants import TONNE_TO_MEGATONNE
+from workflow.scripts.constants import KILOTONNE_TO_MEGATONNE
 
 _BUS_COL_PATTERN = re.compile(r"^bus(\d+)$")
 
@@ -136,8 +136,8 @@ def build_ghg_links_dataframe(
     gwp = pd.Series(
         {
             "emission:co2": 1.0,
-            "emission:ch4": ch4_gwp * TONNE_TO_MEGATONNE,
-            "emission:n2o": n2o_gwp * TONNE_TO_MEGATONNE,
+            "emission:ch4": ch4_gwp * KILOTONNE_TO_MEGATONNE,
+            "emission:n2o": n2o_gwp * KILOTONNE_TO_MEGATONNE,
         }
     )