Emissions and temperature

def _get_emissions_constraints(m: AbstractModel) -> Sequence[GeneralConstraint]:
    """
    ## Baseline emissions

    In MIMOSA, emissions are represented by CO<sub>2</sub> emissions only[^1]. The emissions are calculated
    relative to baseline emissions: emissions that would occur in the absence of any climate policy.

    There are two ways to calculate baseline emissions: either directly imported exogenously from
    the SSP scenarios, or calculated from the baseline carbon intensity from the SSPs. The latter
    accounts for the fact that in the absence of climate policy, baseline emissions would go down if
    the GDP goes down, and vice versa. This can be chosen with the parameter [`use_baseline_carbon_intensity`](../parameters.md#emissions.baseline carbon intensity).
    If this parameter is set to true, baseline emissions are calculated as:

    $$
    \\text{baseline emissions}_{t,r} = \\text{baseline carbon intensity}_{t,r} \\cdot \\text{GDP}_{\\text{net}, t,r}
    $$

    ??? info "Baseline emissions values"

        <div style="overflow: scroll;" markdown>
        ``` plotly
        {"file_path": "./assets/plots/baseline_emissions.json"}
        ```
        </div>

    ## Emission reductions
    To calculate the actual emissions per region, the baseline emissions are reduced by a relative
    abatement factor, which in term is determined by the Marginal Abatement Cost curve and the carbon price
    (see [Mitigation](mitigation.md)). The regional emissions are then calculated as:

    $$
    \\text{regional emissions}_{t,r} = (1 - \\text{relative abatement}_{t,r}) \\cdot \\text{baseline emissions}_{t,r}
    $$

    In the first period, reductions are assumed to be zero:

    $$
    \\text{regional emissions}_{0,r} = \\text{baseline emissions}_{0,r}.
    $$

    ## Global and cumulative emissions

    The regional emissions are aggregated to global emissions:

    $$
    \\text{global emissions}_{t} = \\sum_r \\text{regional emissions}_{t,r},
    $$

    which are used to calculate the cumulative emissions. There are two ways to calculate them:
    using trapezoidal integration or by simply adding up all the years. The former is more accurate,
    while the latter is numerically more stable. This is chosen with the parameter [`cumulative_emissions_trapz`](../parameters.md#emissions.cumulative_emissions_trapz).

    === "Trapezoidal integration `default`"

        $$
        \\begin{aligned}
        \\text{cumulative emissions}_{t} & = \\text{cumulative emissions}_{t-1} \\\\ 
        & + \\frac{\\Delta t}{2} \\cdot (\\text{global emissions}_{t} + \\text{global emissions}_{t-1}) ,\\text{for } t > 0.
        \\end{aligned}
        $$

    === "Simple sum"

        $$
        \\text{cumulative emissions}_{t} = \\text{cumulative emissions}_{t-1} + \\Delta t \\cdot \\text{global emissions}_{t}
        $$


    ## Parameters defined in this module
    - use_carbon_intensity_for_baseline: If true, use baseline carbon intensity to calculate baseline emissions. If false, the SSP baseline emissions are used, regardless of lower GDP.. Type: bool. Default: True.
    - cumulative_emissions_trapz: If true, calculate cumulative emissions using trapezoidal interpolation. If false, cum. emissions are simply cum_emissions[t-1] + dt * cum_emissions[t]. This is less accurate, but better for numerical stability. For small dt the  approximation is usually still acceptable.. Type: bool. Default: True.

    [^1]: The effect of other greenhouse gases is implicitly accounted for in the TCRE which
        translates cumulative CO<sub>2</sub> emissions into temperature change. This assumes a linear
        relation between CO<sub>2</sub> emissions and other greenhouse gases.
    """

    constraints = []

    m.regional_emissions = Var(m.t, m.regions, units=quant.unit("emissionsrate_unit"))
    m.baseline = Var(
        m.t,
        m.regions,
        initialize=lambda m, t, r: m.baseline_emissions[t, r],
        units=quant.unit("emissionsrate_unit"),
    )
    m.use_carbon_intensity_for_baseline = Param(
        doc="::emissions.baseline carbon intensity"
    )

    m.relative_abatement = Var(
        m.t,
        m.regions,
        initialize=0,
        bounds=(0, 2.5),
        units=quant.unit("fraction_of_baseline_emissions"),
    )
    m.regional_emission_reduction = Var(
        m.t, m.regions, units=quant.unit("emissionsrate_unit")
    )
    m.cumulative_emissions = Var(m.t, units=quant.unit("emissions_unit"))
    m.global_emissions = Var(m.t, units=quant.unit("emissionsrate_unit"))

    m.cumulative_emissions_trapz = Param(doc="::emissions.cumulative_emissions_trapz")

    constraints.extend(
        [
            # Baseline emissions based on emissions or carbon intensity
            RegionalConstraint(
                lambda m, t, r: (
                    (
                        m.baseline[t, r]
                        == m.baseline_carbon_intensity[t, r] * m.GDP_net[t, r]
                    )
                    if value(m.use_carbon_intensity_for_baseline)
                    else (m.baseline[t, r] == m.baseline_emissions[t, r])
                ),
                name="baseline_emissions",
            ),
            # Regional emissions from baseline and relative abatement
            RegionalConstraint(
                lambda m, t, r: (
                    m.regional_emissions[t, r]
                    == (1 - m.relative_abatement[t, r])
                    * (
                        m.baseline[t, r]
                        if value(m.use_carbon_intensity_for_baseline)
                        else m.baseline_emissions[t, r]
                        # Note: this should simply be m.baseline[t,r], but this is numerically less stable
                        # than m.baseline_emissions[t, r] whenever baseline intensity
                        # is used instead of baseline emissions. In fact, m.baseline_emissions[t, r]
                        # is just a fixed number, whereas m.baseline[t,r] is a variable depending on
                        # GDP.
                    )
                    if t > 0
                    else Constraint.Skip
                ),
                "regional_abatement",
            ),
            RegionalInitConstraint(
                lambda m, r: m.regional_emissions[0, r] == m.baseline_emissions[0, r]
            ),
            RegionalConstraint(
                lambda m, t, r: m.regional_emission_reduction[t, r]
                == m.baseline[t, r] - m.regional_emissions[t, r],
                "regional_emission_reduction",
            ),
            # Global emissions (sum from regional emissions)
            GlobalConstraint(
                lambda m, t: (
                    m.global_emissions[t]
                    == sum(m.regional_emissions[t, r] for r in m.regions)
                    if t > 0
                    else Constraint.Skip
                ),
                "global_emissions",
            ),
            GlobalInitConstraint(
                lambda m: m.global_emissions[0]
                == sum(m.baseline_emissions[0, r] for r in m.regions),
                "global_emissions_init",
            ),
            # Cumulative global emissions
            GlobalConstraint(
                lambda m, t: (
                    m.cumulative_emissions[t]
                    == m.cumulative_emissions[t - 1]
                    + (
                        (m.dt * (m.global_emissions[t] + m.global_emissions[t - 1]) / 2)
                        if value(m.cumulative_emissions_trapz)
                        else (m.dt * m.global_emissions[t])
                    )
                    if t > 0
                    else Constraint.Skip
                ),
                "cumulative_emissions",
            ),
            GlobalInitConstraint(lambda m: m.cumulative_emissions[0] == 0),
        ]
    )

    m.emission_relative_cumulative = Var(m.t, initialize=1)
    constraints.extend(
        [
            GlobalConstraint(
                lambda m, t: (
                    (
                        m.emission_relative_cumulative[t]
                        == m.cumulative_emissions[t]
                        / m.cumulative_global_baseline_emissions[t]
                    )
                    if t > 0
                    else Constraint.Skip
                ),
                name="relative_cumulative_emissions",
            ),
            GlobalInitConstraint(lambda m: m.emission_relative_cumulative[0] == 1),
        ]
    )

    return constraints

def _get_temperature_constraints(m: AbstractModel) -> Sequence[GeneralConstraint]:
    """

    The global temperature change is calculated as a linear function of cumulative emissions, with a slope
    given by the Transient Climate Response to CO<sub>2</sub> Emissions (TCRE), following [Dietz et al.] who
    showed that this performed better than the default two-box DICE climate module.

    $$
    \\text{temperature}_{t} = T_0 + \\text{TCRE} \\cdot \\text{cumulative emissions}_{t},
    $$

    where [$T_0$](../parameters.md#temperature.initial) is the initial temperature at the start of the run (by default in 2020),
    and the [TCRE](../parameters.md#temperature.TCRE) is the Transient Climate Response to CO<sub>2</sub> Emissions.

    The initial temperature is set to 1.16°C in 2020 by default, following [Visser et al.]. The TCRE
    is calibrated on the IPCC AR5 or AR6 reports (the median value of the TCRE is the same in the AR5
    and AR6 calibration), but the distribution is different.

    ??? info "Calibration of the TCRE"

        The TCRE is calibrated on the IPCC AR5 and AR6 reports, both using the linear relation shown in the SPM figure of the WG1 report.
        The default value of the TCRE is 0.62°C per 1000 GtCO<sub>2</sub> (median value of the AR5 and AR6 calibration). The 5th and 95th percentiles
        differ between the two calibrations, with the AR5 calibration having a wider distribution:

        | Percentile | AR5 TCRE <br><span style="font-weight: normal">(°C per TtCO<sub>2</sub>)</span> | AR6 TCRE <br><span style="font-weight: normal">(°C per TtCO<sub>2</sub>)</span> |
        |------------|----------|----------|
        | 5th        | 0.42     | 0.42     |
        | 50th       | 0.62     | 0.62 `default` |
        | 95th       | 0.82     | 0.75     |

        === "Calibrated on IPCC AR5"

            <div style="overflow: scroll;" markdown>
            ``` plotly
            {"file_path": "./assets/plots/ar5_tcre.json"}
            ```
            </div>
            Source: [IPCC AR5 WG1 Figure SPM.10](https://www.ipcc.ch/report/ar5/wg1/summary-for-policymakers/figspm-10/)

        === "Calibrated on IPCC AR6"

            <div style="overflow: scroll;" markdown>
            ``` plotly
            {"file_path": "./assets/plots/ar6_tcre.json"}
            ```
            </div>

            Source: [IPCC AR6 WG1 Figure SPM.10](https://www.ipcc.ch/report/ar6/wg1/figures/summary-for-policymakers/figure-spm-10)

    ## Temperature target
    The usual way to specify a policy target in MIMOSA is using carbon budgets. However, it is also possible to
    specify a temperature target. This is done by setting the parameter [`temperature_target`](../parameters.md#temperature.target).
    This is an upper bound on the temperature: if the cost-optimal temperature is below this value, this constraint is not binding.

    If this parameter is set, the following constraint is active:

    $$
    \\text{temperature}_{t} \\leq \\text{temperature target}
    $$

    ## Parameters defined in this module
    - T0: Temperature in initial year of model run (2020 by default).. Type: quantity. Default: 1.16 delta_degC. Unit: temperature_unit.
    - TCRE: Transient Climate Response to CO2 Emissions. Type: quantity. Default: 0.62 delta_degC/(TtCO2). Unit: (temperature_unit)/(emissions_unit).
    - temperature_target: Temperature target in 2100 (and beyond). Can also be false, then no temperature target is imposed. Type: quantity. Default: False. Can also be false. Unit: temperature_unit.
    """

    constraints = []

    m.T0 = Param(units=quant.unit("degC_above_PI"), doc="::temperature.initial")
    m.temperature = Var(
        m.t, initialize=lambda m, t: m.T0, units=quant.unit("degC_above_PI")
    )
    m.TCRE = Param(doc="::temperature.TCRE")
    m.temperature_target = Param(doc="::temperature.target")
    constraints.extend(
        [
            GlobalConstraint(
                lambda m, t: m.temperature[t]
                == m.T0 + m.TCRE * m.cumulative_emissions[t],
                "temperature",
            ),
            GlobalInitConstraint(lambda m: m.temperature[0] == m.T0),
            GlobalConstraint(
                lambda m, t: (
                    m.temperature[t] <= m.temperature_target
                    if (m.year(t) >= 2100 and value(m.temperature_target) is not False)
                    else Constraint.Skip
                ),
                name="temperature_target",
            ),
        ]
    )

    m.perc_reversible_damages = Param(doc="::economics.damages.percentage reversible")

    # m.overshoot = Var(m.t, initialize=0)
    # m.overshootdot = DerivativeVar(m.overshoot, wrt=m.t)
    # m.netnegative_emissions = Var(m.t)
    # global_constraints.extend(
    #     [
    #         lambda m, t: m.netnegative_emissions[t]
    #         == m.global_emissions[t] * (1 - tanh(m.global_emissions[t] * 10)) / 2
    #         if value(m.perc_reversible_damages) < 1
    #         else Constraint.Skip,
    #         lambda m, t: m.overshootdot[t]
    #         == (m.netnegative_emissions[t] if t <= value(m.year2100) and t > 0 else 0)
    #         if value(m.perc_reversible_damages) < 1
    #         else Constraint.Skip,
    #     ]
    # )

    # global_constraints_init.extend([lambda m: m.overshoot[0] == 0])

    return constraints

def _get_inertia_and_budget_constraints(
    m: AbstractModel,
) -> Sequence[GeneralConstraint]:
    """
    MIMOSA allows several types of constraints on emissions: a global carbon budget, inertia constraints,
    limits on net negative emissions and constraints on emissions after 2100.

    ## Carbon budget

    By default, the carbon budget is set to `False`, which means that MIMOSA runs in CBA mode without carbon budget.
    If a carbon budget is defined using the parameter [`carbonbudget`](../parameters.md#emissions.carbonbudget), the following constraint is active:

    $$
    \\text{cumulative emissions}_{t} - \\text{budget} \\leq 0,
    $$

    for $t \\geq 2100$. Note that the [`carbonbudget`](../parameters.md#emissions.carbonbudget) parameter unit should be
    in GtCO<sub>2</sub> (or TtCO<sub>2</sub> or MtCO<sub>2</sub>).

    ## Inertia

    In MIMOSA, there is a limit on how fast emissions can be reduced. This constraint is called inertia, and represents the
    technological, social, political and economic difficulties of quickly reducing emissions. It is implemented as a hard limit
    on the speed of emission reductions. The parameters [`inertia_global`](../parameters.md#emissions.inertia.global) and
    [`inertia_regional`](../parameters.md#emissions.inertia.regional) are defined as percentages of baseline emissions in the
    initial year. If baseline emissions in a region are 10 GtCO<sub>2</sub>/yr in the initial year, and the global inertia is 5%,
    then the global emissions cannot be reduced by more than 0.5 GtCO<sub>2</sub>/yr throughout the whole run. This way, the parameter
    is independent of the size of the region.

    === "Global inertia"


        $$
        \\text{global emissions}_{t} - \\text{global emissions}_{t-1} \\geq \\Delta t \\cdot \\text{inertia_global} \\cdot \\text{scaling factor},
        $$

        with the scaling factor the global baseline emissions in the starting year:

        $$
        \\text{scaling factor} = \\left( \\sum_r \\text{baseline emissions}_{t=0,r}\\right)
        $$

        By default, the [`inertia_global`](../parameters.md#emissions.inertia.global) is set to `False`, such that there is no global inertia
        constraint. The only inertia constraint happens regionally by default.

    === "Regional inertia"

        $$
        \\text{regional emissions}_{t,r} - \\text{regional emissions}_{t-1,r} \\geq \\Delta t \\cdot \\text{inertia_regional} \\cdot \\text{baseline emissions}_{t=0,r}.
        $$

        By default, the [`inertia_regional`](../parameters.md#emissions.inertia.regional) parameter is active, such that regional emissions cannot be reduced by more than 5% of the initial baseline emissions per year. This is
        based on maximum reduction speeds of the scenarios in the scenario explorer for 1.5°C pathways underpinning the IPCC Special Report on Global Warming
        of 1.5°C (<https://data.ene.iiasa.ac.at/iamc-1.5c-explorer>) (ref https://www.frontiersin.org/articles/10.3389/fclim.2021.785577/full)


    ## Minimum emission levels (limits to net-negative emissions)

    Since the Marginal Abatement Cost curve is a continuous function with no upper bound, reductions can theoretically be without bound.
    To still address the difficulties associated with CDR technologies, a limit on (net) negative emissions can be imposed. These difficulties
    can be of economic (e.g., land becoming increasingly scarce or increased dependence on very expensive storage sites) and socio-political
    (concerns about biodiversity and food security) nature (TODO ref Hotelling-Hof-Wijst paper and Fuss et al.). By the default, the global emissions
    are limited to -20 GtCO<sub>2</sub>/yr, and regionally to -10 GtCO<sub>2</sub>:

    $$
    \\text{global emissions}_t \\geq \\text{global min level},
    $$

    $$
    \\text{regional emissions}_t \\geq \\text{regional min level}.
    $$

    Note that currently, the minimum regional emission level is the same for every region, but this could be easily changed if needed (see
    [Extending MIMOSA](../extending/index.md)).

    ## Constraints on emissions after 2100

    In CBA mode, it can be useful to run MIMOSA until 2150 instead of 2100, to avoid end-of-horizon effects. One way to mitigate
    this effect is to impose a constraint on the growth of emissions after 2100. By default, emissions are not allowed to grow after 2100:

    $$
    \\text{regional emissions}_{t,r} - \\text{regional emissions}_{t-1,r} \\leq 0, \\text{if } t - 1 > 2100,
    $$

    if the parameter [`non_increasing_emissions_after_2100`](../parameters.md#emissions.non increasing emissions after 2100) is set to `True`.

    When a carbon budget is imposed, it is also possible to impose a net-zero emissions constraint after the budget year (2100):

    $$
    \\text{global emissions}_{t} \\leq 0, \\text{if } t \\geq 2100,
    $$
    if the parameter [`no_pos_emissions_after_budget_year`](../parameters.md#emissions.not positive after budget year) is set to `True`,
    **and** if a carbon budget is specified. Otherwise, this constraint is ignored.


    ## Parameters defined in this module
    - budget: Value of the carbon budget. Example: "800 GtCO2" (the unit is important). If set to False, no carbon budget is imposed: this is cost-benefit mode. Default: False.. Type: quantity. Default: False. Can also be false. Unit: emissions_unit.
    - inertia_global: Maximum reduction speed, in % of initial emissions (should be negative) Can also be false, then no inertia limit is imposed. Type: float. Default: False. Can also be false. Min: -inf. Max: 0.
    - inertia_regional: Maximum reduction speed, in % of initial emissions (should be negative) Can also be false, then no inertia limit is imposed. Type: float. Default: -0.05. Can also be false. Min: -inf. Max: 0.
    - global_min_level: Limit on the emission level (globally), mostly used for negative emissions. Can also be false, then no limit is imposed. Type: quantity. Default: -20 GtCO2/yr. Can also be false. Unit: emissionsrate_unit.
    - regional_min_level: Limit on the emission level (per region), mostly used for negative emissions. Can also be false, then no limit is imposed. Type: quantity. Default: -10 GtCO2/yr. Can also be false. Unit: emissionsrate_unit.
    - no_pos_emissions_after_budget_year: If true, impose net-zero emissions after budget year (2100). Type: bool. Default: True.
    - non_increasing_emissions_after_2100: If true, the regional emissions after 2100 are not allowed to climb.. Type: bool. Default: True.

    """

    constraints = []

    m.budget = Param(doc="::emissions.carbonbudget")
    m.inertia_global = Param(doc="::emissions.inertia.global")
    m.inertia_regional = Param(doc="::emissions.inertia.regional")
    m.global_min_level = Param(doc="::emissions.global min level")
    m.regional_min_level = Param(doc="::emissions.regional min level")
    m.non_increasing_emissions_after_2100 = Param(
        doc="::emissions.non increasing emissions after 2100"
    )
    m.no_pos_emissions_after_budget_year = Param(
        doc="::emissions.not positive after budget year"
    )
    constraints.extend(
        [
            # Carbon budget constraints:
            GlobalConstraint(
                lambda m, t: (
                    m.cumulative_emissions[t]
                    - (
                        m.budget
                        + (
                            m.overshoot[t] * (1 - m.perc_reversible_damages)
                            if value(m.perc_reversible_damages) < 1
                            else 0
                        )
                    )
                    <= 0
                    if (m.year(t) >= 2100 and value(m.budget) is not False)
                    else Constraint.Skip
                ),
                name="carbon_budget",
            ),
            GlobalConstraint(lambda m, t: m.cumulative_emissions[t] >= 0),
            # Global and regional inertia constraints:
            GlobalConstraint(
                lambda m, t: (
                    m.global_emissions[t] - m.global_emissions[t - 1]
                    >= m.dt
                    * m.inertia_global
                    * sum(m.baseline_emissions[0, r] for r in m.regions)
                    if value(m.inertia_global) is not False and t > 0
                    else Constraint.Skip
                ),
                name="global_inertia",
            ),
            RegionalConstraint(
                lambda m, t, r: (
                    m.regional_emissions[t, r] - m.regional_emissions[t - 1, r]
                    >= m.dt * m.inertia_regional * m.baseline_emissions[0, r]
                    if value(m.inertia_regional) is not False and t > 0
                    else Constraint.Skip
                ),
                name="regional_inertia",
            ),
            GlobalConstraint(
                lambda m, t: (
                    m.global_emissions[t] >= m.global_min_level
                    if value(m.global_min_level) is not False
                    else Constraint.Skip
                ),
                "global_min_level",
            ),
            RegionalConstraint(
                lambda m, t, r: (
                    m.regional_emissions[t, r] >= m.regional_min_level
                    if value(m.regional_min_level) is not False
                    else Constraint.Skip
                ),
                "regional_min_level",
            ),
            RegionalConstraint(
                lambda m, t, r: (
                    m.regional_emissions[t, r] - m.regional_emissions[t - 1, r] <= 0
                    if m.year(t - 1) > 2100
                    and value(m.non_increasing_emissions_after_2100)
                    else Constraint.Skip
                ),
                name="non_increasing_emissions_after_2100",
            ),
            GlobalConstraint(
                lambda m, t: (
                    m.global_emissions[t] <= 0
                    if (
                        m.year(t) >= 2100
                        and value(m.no_pos_emissions_after_budget_year) is True
                        and value(m.budget) is not False
                    )
                    else Constraint.Skip
                ),
                name="net_zero_after_2100",
            ),
        ]
    )

    return constraints