Merge pull request #501 from EnergySystemsModellingLab/carbon_budget_fix

Small improvements to carbon budget algorithm

Author: tsmbland

docs/inputs/toml.rst

@@ -148,6 +148,7 @@ Example
 
    - `refine_price`: If True, applies an upper limit on the carbon price. Defaults to False.
    - `price_too_high_threshold`: upper limit on the carbon price. Defaults to 10.
+   - `resolution`: Number of decimal places to solve the carbon price to. When using the bisection method, increasing this value may increase the time taken to solve the carbon market. Defaults to 2.
 
 
 ------------------

docs/tutorial-code/carbon-budget/1-carbon-budget/Results/MCAPrices.csv

@@ -47,51 +47,51 @@ electricity,all-week,evening,all-year,3.98810000000,R1,5,MUS$2010/PJ,2025
 gas,all-week,evening,all-year,0.47220000000,R1,5,MUS$2010/PJ,2025
 heat,all-week,evening,all-year,4.91770000000,R1,5,MUS$2010/PJ,2025
 CO2f,all-week,evening,all-year,0.21000000000,R1,5,MUS$2010/kt,2025
-electricity,all-week,night,all-year,9.02490000000,R1,0,MUS$2010/PJ,2030
+electricity,all-week,night,all-year,8.60540000000,R1,0,MUS$2010/PJ,2030
 gas,all-week,night,all-year,0.47220000000,R1,0,MUS$2010/PJ,2030
-heat,all-week,night,all-year,3.85050000000,R1,0,MUS$2010/PJ,2030
-CO2f,all-week,night,all-year,0.18000000000,R1,0,MUS$2010/kt,2030
-electricity,all-week,morning,all-year,9.02490000000,R1,1,MUS$2010/PJ,2030
+heat,all-week,night,all-year,3.68270000000,R1,0,MUS$2010/PJ,2030
+CO2f,all-week,night,all-year,0.17000000000,R1,0,MUS$2010/kt,2030
+electricity,all-week,morning,all-year,8.60540000000,R1,1,MUS$2010/PJ,2030
 gas,all-week,morning,all-year,0.47220000000,R1,1,MUS$2010/PJ,2030
-heat,all-week,morning,all-year,3.85050000000,R1,1,MUS$2010/PJ,2030
-CO2f,all-week,morning,all-year,0.18000000000,R1,1,MUS$2010/kt,2030
-electricity,all-week,afternoon,all-year,9.02490000000,R1,2,MUS$2010/PJ,2030
+heat,all-week,morning,all-year,3.68270000000,R1,1,MUS$2010/PJ,2030
+CO2f,all-week,morning,all-year,0.17000000000,R1,1,MUS$2010/kt,2030
+electricity,all-week,afternoon,all-year,8.60540000000,R1,2,MUS$2010/PJ,2030
 gas,all-week,afternoon,all-year,0.47220000000,R1,2,MUS$2010/PJ,2030
-heat,all-week,afternoon,all-year,3.85050000000,R1,2,MUS$2010/PJ,2030
-CO2f,all-week,afternoon,all-year,0.18000000000,R1,2,MUS$2010/kt,2030
-electricity,all-week,early-peak,all-year,9.02490000000,R1,3,MUS$2010/PJ,2030
+heat,all-week,afternoon,all-year,3.68270000000,R1,2,MUS$2010/PJ,2030
+CO2f,all-week,afternoon,all-year,0.17000000000,R1,2,MUS$2010/kt,2030
+electricity,all-week,early-peak,all-year,8.60540000000,R1,3,MUS$2010/PJ,2030
 gas,all-week,early-peak,all-year,0.47220000000,R1,3,MUS$2010/PJ,2030
-heat,all-week,early-peak,all-year,3.85050000000,R1,3,MUS$2010/PJ,2030
-CO2f,all-week,early-peak,all-year,0.18000000000,R1,3,MUS$2010/kt,2030
-electricity,all-week,late-peak,all-year,9.02490000000,R1,4,MUS$2010/PJ,2030
+heat,all-week,early-peak,all-year,3.68270000000,R1,3,MUS$2010/PJ,2030
+CO2f,all-week,early-peak,all-year,0.17000000000,R1,3,MUS$2010/kt,2030
+electricity,all-week,late-peak,all-year,8.60540000000,R1,4,MUS$2010/PJ,2030
 gas,all-week,late-peak,all-year,0.47220000000,R1,4,MUS$2010/PJ,2030
-heat,all-week,late-peak,all-year,3.85050000000,R1,4,MUS$2010/PJ,2030
-CO2f,all-week,late-peak,all-year,0.18000000000,R1,4,MUS$2010/kt,2030
-electricity,all-week,evening,all-year,9.02490000000,R1,5,MUS$2010/PJ,2030
+heat,all-week,late-peak,all-year,3.68270000000,R1,4,MUS$2010/PJ,2030
+CO2f,all-week,late-peak,all-year,0.17000000000,R1,4,MUS$2010/kt,2030
+electricity,all-week,evening,all-year,8.60540000000,R1,5,MUS$2010/PJ,2030
 gas,all-week,evening,all-year,0.47220000000,R1,5,MUS$2010/PJ,2030
-heat,all-week,evening,all-year,3.85050000000,R1,5,MUS$2010/PJ,2030
-CO2f,all-week,evening,all-year,0.18000000000,R1,5,MUS$2010/kt,2030
-electricity,all-week,night,all-year,11.61590000000,R1,0,MUS$2010/PJ,2035
+heat,all-week,evening,all-year,3.68270000000,R1,5,MUS$2010/PJ,2030
+CO2f,all-week,evening,all-year,0.17000000000,R1,5,MUS$2010/kt,2030
+electricity,all-week,night,all-year,10.91870000000,R1,0,MUS$2010/PJ,2035
 gas,all-week,night,all-year,0.47220000000,R1,0,MUS$2010/PJ,2035
-heat,all-week,night,all-year,4.88690000000,R1,0,MUS$2010/PJ,2035
-CO2f,all-week,night,all-year,0.29000000000,R1,0,MUS$2010/kt,2035
-electricity,all-week,morning,all-year,11.61590000000,R1,1,MUS$2010/PJ,2035
+heat,all-week,night,all-year,4.60800000000,R1,0,MUS$2010/PJ,2035
+CO2f,all-week,night,all-year,0.27000000000,R1,0,MUS$2010/kt,2035
+electricity,all-week,morning,all-year,10.91870000000,R1,1,MUS$2010/PJ,2035
 gas,all-week,morning,all-year,0.47220000000,R1,1,MUS$2010/PJ,2035
-heat,all-week,morning,all-year,4.88690000000,R1,1,MUS$2010/PJ,2035
-CO2f,all-week,morning,all-year,0.29000000000,R1,1,MUS$2010/kt,2035
-electricity,all-week,afternoon,all-year,11.61590000000,R1,2,MUS$2010/PJ,2035
+heat,all-week,morning,all-year,4.60800000000,R1,1,MUS$2010/PJ,2035
+CO2f,all-week,morning,all-year,0.27000000000,R1,1,MUS$2010/kt,2035
+electricity,all-week,afternoon,all-year,10.91870000000,R1,2,MUS$2010/PJ,2035
 gas,all-week,afternoon,all-year,0.47220000000,R1,2,MUS$2010/PJ,2035
-heat,all-week,afternoon,all-year,4.88690000000,R1,2,MUS$2010/PJ,2035
-CO2f,all-week,afternoon,all-year,0.29000000000,R1,2,MUS$2010/kt,2035
-electricity,all-week,early-peak,all-year,11.61590000000,R1,3,MUS$2010/PJ,2035
+heat,all-week,afternoon,all-year,4.60800000000,R1,2,MUS$2010/PJ,2035
+CO2f,all-week,afternoon,all-year,0.27000000000,R1,2,MUS$2010/kt,2035
+electricity,all-week,early-peak,all-year,10.91870000000,R1,3,MUS$2010/PJ,2035
 gas,all-week,early-peak,all-year,0.47220000000,R1,3,MUS$2010/PJ,2035
-heat,all-week,early-peak,all-year,4.88690000000,R1,3,MUS$2010/PJ,2035
-CO2f,all-week,early-peak,all-year,0.29000000000,R1,3,MUS$2010/kt,2035
-electricity,all-week,late-peak,all-year,11.61590000000,R1,4,MUS$2010/PJ,2035
+heat,all-week,early-peak,all-year,4.60800000000,R1,3,MUS$2010/PJ,2035
+CO2f,all-week,early-peak,all-year,0.27000000000,R1,3,MUS$2010/kt,2035
+electricity,all-week,late-peak,all-year,10.91870000000,R1,4,MUS$2010/PJ,2035
 gas,all-week,late-peak,all-year,0.47220000000,R1,4,MUS$2010/PJ,2035
-heat,all-week,late-peak,all-year,4.88690000000,R1,4,MUS$2010/PJ,2035
-CO2f,all-week,late-peak,all-year,0.29000000000,R1,4,MUS$2010/kt,2035
-electricity,all-week,evening,all-year,11.61590000000,R1,5,MUS$2010/PJ,2035
+heat,all-week,late-peak,all-year,4.60800000000,R1,4,MUS$2010/PJ,2035
+CO2f,all-week,late-peak,all-year,0.27000000000,R1,4,MUS$2010/kt,2035
+electricity,all-week,evening,all-year,10.91870000000,R1,5,MUS$2010/PJ,2035
 gas,all-week,evening,all-year,0.47220000000,R1,5,MUS$2010/PJ,2035
-heat,all-week,evening,all-year,4.88690000000,R1,5,MUS$2010/PJ,2035
-CO2f,all-week,evening,all-year,0.29000000000,R1,5,MUS$2010/kt,2035
+heat,all-week,evening,all-year,4.60800000000,R1,5,MUS$2010/PJ,2035
+CO2f,all-week,evening,all-year,0.27000000000,R1,5,MUS$2010/kt,2035

src/muse/carbon_budget.py

@@ -77,6 +77,7 @@ def fitting(
     price_too_high_threshold: float = 10,
     sample_size: int = 5,
     fitter: str = "linear",
+    resolution: int = 2,
 ) -> float:
     """Used to solve the carbon market.
 
@@ -95,6 +96,7 @@ def fitting(
         price_too_high_threshold: threshold on carbon price
         sample_size: sample size for fitting
         fitter: method to fit emissions with carbon price
+        resolution: Number of decimal places to solve the carbon price to
 
     Returns:
         new_price: adjusted carbon price to meet budget
@@ -125,12 +127,12 @@ def fitting(
         threshold,  # type: ignore
     )
 
-    # Cap price between 0.01 and price_too_high_threshold
+    # Cap price between 0.0 and price_too_high_threshold
     if refine_price:
         new_price = min(new_price, price_too_high_threshold)
-    new_price = max(new_price, 0.01)
+    new_price = max(new_price, 0.0)
 
-    return new_price
+    return round(new_price, resolution)
 
 
 def linear_fun(x, a, b):
@@ -297,6 +299,7 @@ def bisection(
     max_iterations: int = 5,
     tolerance: float = 0.1,
     early_termination_count: int = 5,
+    resolution: int = 2,
 ) -> float:
     """Applies bisection algorithm to escalate carbon price and meet the budget.
 
@@ -319,36 +322,45 @@ def bisection(
         tolerance: Maximum permitted deviation of emissions from the budget
         early_termination_count: Will terminate the loop early if the last n solutions
             are the same
+        resolution: Number of decimal places to solve the carbon price to
 
     Returns:
         New value of global carbon price
     """
     from logging import getLogger
 
+    # Create cache for emissions at different price points
+    emissions_cache = EmissionsCache(market, equilibrium, commodities)
+
     # Carbon price and emissions threshold in the forecast year
     future = market.year[-1]
     target = carbon_budget.sel(year=future).values.item()
     price = market.prices.sel(year=future, commodity=commodities).mean().values.item()
 
+    # Test if emissions are already below the budget without imposing a carbon price
+    if emissions_cache[0.0] < target:
+        message = (
+            f"Emissions for the year {int(future)} are already below the carbon budget "
+            "without imposing a carbon price. The carbon price has been set to zero."
+        )
+        getLogger(__name__).warning(message)
+        return 0.0
+
     # Initial lower and upper bounds on carbon price for the bisection algorithm
     current = market.year[0]
     time_exp = int(future - current)
-    lb_price = 0.01
-    ub_price = (
-        max(price, 0.01) * 1.1**time_exp
+    lb_price = 0.0
+    epsilon = 10**-resolution  # smallest nonzero price
+    ub_price = round(
+        (max(price, epsilon) * 1.1**time_exp), resolution
     )  # i.e. 10% yearly increase on current price
 
     # Bisection loop
-    emissions_cache = EmissionsCache(market, equilibrium, commodities)
     for _ in range(max_iterations):  # maximum number of iterations before terminating
-        # Cap prices between 0.01 and price_too_high_threshold
+        # Cap prices between 0.0 and price_too_high_threshold
         if refine_price:
             ub_price = min(ub_price, price_too_high_threshold)
-        lb_price = max(lb_price, 0.01)
-
-        # Round prices to 2dp
-        lb_price = round(lb_price, 2)
-        ub_price = round(ub_price, 2)
+        lb_price = max(lb_price, 0.0)
 
         # Calculate/retrieve carbon emissions at new bounds
         lb_price_emissions = emissions_cache[lb_price]
@@ -372,16 +384,35 @@ def bisection(
             ub_price,
             emissions_cache,
             target,
+            resolution,
         )
 
     # If convergence isn't reached, new price is that with emissions closest to
     # threshold. If multiple prices are equally close, it returns the lowest price
-    message = (
-        f"Carbon budget could not be matched for year {int(future)}. "
-        "Try increasing the tolerance or sample size."
+    new_price = min(
+        emissions_cache, key=lambda k: (abs(emissions_cache[k] - target), k)
     )
+
+    # Raise warning message
+    if all(emissions_cache[k] > target for k in emissions_cache):
+        message = (
+            f"Carbon budget could not be met for the year {int(future)} "
+            f"(budget: {target}, emissions: {emissions_cache[new_price]}). "
+            "This may be because there are no processes available that can meet the "
+            "budget, or because emissions from capacity installed earlier in the time "
+            "horizon is preventing the budget from being met. "
+            "The CO2 price in this year should be interpreted with caution."
+        )
+    else:
+        message = (
+            f"Carbon budget could not be matched for the year {int(future)} to within "
+            "the specified tolerance. "
+            "This is sometimes unavoidable due to a discontinuous emissions landscape "
+            "which can make the budget unreachable, but can sometimes be "
+            "fixed by increasing max_iterations, early_termination_count or resolution."
+        )
     getLogger(__name__).warning(message)
-    return min(emissions_cache, key=lambda k: (abs(emissions_cache[k] - target), k))
+    return new_price
 
 
 class EmissionsCache(dict):
@@ -408,6 +439,7 @@ def adjust_bounds(
     ub_price: float,
     emissions_cache: dict[float, float],
     target: float,
+    resolution: int = 2,
 ) -> tuple[float, float]:
     """Adjust the bounds of the carbon price for the bisection algorithm.
 
@@ -420,6 +452,7 @@ def adjust_bounds(
         ub_price: Value of carbon price at upper bound
         emissions_cache: Dictionary of emissions at different price points
         target: Carbon budget
+        resolution: Number of decimal places to solve the carbon price to
 
     Returns:
         New lower and upper bounds for the carbon price.
@@ -454,32 +487,41 @@ def adjust_bounds(
         ub_price,
         emissions_cache,
         target,
+        resolution,
     )
 
 
 def decrease_bounds(
-    lb_price: float, ub_price: float, emissions_cache: dict[float, float], target: float
+    lb_price: float,
+    ub_price: float,
+    emissions_cache: dict[float, float],
+    target: float,
+    resolution: int = 2,
 ) -> tuple[float, float]:
     """Decreases the lb of the carbon price, and sets the ub to the previous lb."""
     denominator = max(target, 1e-3)
     lb_price_emissions = emissions_cache[lb_price]
     exponent = (lb_price_emissions - target) / abs(denominator)  # will be negative
     exponent = max(exponent, -1)  # cap exponent at -1
     ub_price = lb_price
-    lb_price = lb_price * np.exp(exponent)
+    lb_price = round(lb_price * np.exp(exponent), resolution)
     return lb_price, ub_price
 
 
 def increase_bounds(
-    lb_price: float, ub_price: float, emissions_cache: dict[float, float], target: float
+    lb_price: float,
+    ub_price: float,
+    emissions_cache: dict[float, float],
+    target: float,
+    resolution: int = 2,
 ) -> tuple[float, float]:
     """Increases the ub of the carbon price, and sets the lb to the previous ub."""
     denominator = max(target, 1e-3)
     ub_price_emissions = emissions_cache[ub_price]
     exponent = (ub_price_emissions - target) / abs(denominator)  # will be positive
     exponent = min(exponent, 1)  # cap exponent at 1
     lb_price = ub_price
-    ub_price = ub_price * np.exp(exponent)
+    ub_price = round(ub_price * np.exp(exponent), resolution)
     return lb_price, ub_price
 
 
@@ -488,9 +530,10 @@ def bisect_bounds(
     ub_price: float,
     emissions_cache: dict[float, float],
     target: float,
+    resolution: int = 2,
 ) -> tuple[float, float]:
     """Bisects the bounds of the carbon price."""
-    midpoint = round((lb_price + ub_price) / 2.0, 2)
+    midpoint = round((lb_price + ub_price) / 2.0, resolution)
     midpoint_emissions = emissions_cache[midpoint]
     if midpoint_emissions < target:
         ub_price = midpoint
@@ -504,9 +547,10 @@ def bisect_bounds_inverted(
     ub_price: float,
     emissions_cache: dict[float, float],
     target: float,
+    resolution: int = 2,
 ) -> tuple[float, float]:
     """Bisects the bounds of the carbon price, in the case of inverted bounds."""
-    midpoint = round((lb_price + ub_price) / 2.0, 2)
+    midpoint = round((lb_price + ub_price) / 2.0, resolution)
     midpoint_emissions = emissions_cache[midpoint]
     if midpoint_emissions > target:
         ub_price = midpoint