Resolve the issue of zero capacity factor

RELEASE_NOTES.rst

@@ -11,6 +11,7 @@ All changes
 -----------
 
 - Correct typo in GAMS formulation, :ref:`equation_renewables_equivalence` (:pull:`581`).
+- Handle zero values in ``capacity_factor`` in models with sub-annual time resolution; expand tests (:issue:`515`, :pull:`561`).
 - Improve configurability of :mod:`.macro`; see the :doc:`documentation <macro>` (:pull:`327`).
 - Split :meth:`.Reporter.add_tasks` for use without an underlying :class:.`Scenario` (:pull:`567`).
 - Allow setting the “model_dir” and “solve_options” options for :class:`.GAMSModel` (and subclasses :class:`.MESSAGE`, :class:`.MACRO`, and :class:`.MESSAGE_MACRO`) through the user's ixmp configuration file; expand documentation (:pull:`557`).

message_ix/model/MESSAGE/data_load.gms

@@ -12,6 +12,7 @@ $LOAD emission, land_scenario, land_type, relation
 $LOAD level_resource, level_renewable
 $LOAD lvl_spatial, lvl_temporal, map_spatial_hierarchy, map_temporal_hierarchy
 $LOAD map_node, map_time, map_commodity, map_resource, map_stocks, map_tec, map_tec_time, map_tec_mode
+$LOAD is_capacity_factor
 $LOAD map_land, map_relation
 $LOAD type_tec, cat_tec, type_year, cat_year, type_emission, cat_emission, type_tec_land
 $LOAD inv_tec, renewable_tec
@@ -142,7 +143,7 @@ map_rating(node,inv_tec,commodity,level,rating,year_all)$(
 
 * set the default capacity factor for technologies where no parameter value is provided in the input data
 capacity_factor(node,tec,year_all2,year_all,time)$( map_tec_time(node,tec,year_all,time)
-    AND map_tec_lifetime(node,tec,year_all2,year_all) AND NOT capacity_factor(node,tec,year_all2,year_all,time) ) = 1 ;
+    AND map_tec_lifetime(node,tec,year_all2,year_all) AND NOT is_capacity_factor(node,tec,year_all2,year_all,time) ) = 1 ;
 
 * assign the yearly average capacity factor (used in equation OPERATION_CONSTRAINT)
 capacity_factor(node,tec,year_all2,year_all,'year') =

message_ix/model/MESSAGE/sets_maps_def.gms

@@ -444,6 +444,7 @@ Sets
     is_dynamic_new_capacity_lo(node,tec,year_all)    flag whether lower dynamic constraint exists for new capacity (investment)
     is_dynamic_activity_up(node,tec,year_all,time)   flag whether upper dynamic constraint exists for a technology (activity)
     is_dynamic_activity_lo(node,tec,year_all,time)   flag whether lower dynamic constraint exists for a technology (activity)
+    is_capacity_factor(node,tec,year_all2,year_all,time)     flag whether capacity factor is defined
 
     is_bound_emission(node,type_emission,type_tec,type_year) flag whether emissions bound exists
 

message_ix/models.py

@@ -118,6 +118,7 @@ def item(ix_type, expr):
     "cat_relation": item("set", "type_relation relation"),
     "cat_tec": item("set", "type_tec t"),
     "cat_year": item("set", "type_year y"),
+    "is_capacity_factor": item("set", "nl t yv ya h"),
     "level_renewable": dict(ix_type="set", idx_sets=["level"]),
     "level_resource": dict(ix_type="set", idx_sets=["level"]),
     "level_stocks": dict(ix_type="set", idx_sets=["level"]),
@@ -331,6 +332,13 @@ def initialize(cls, scenario):
         # Initialize the ixmp items
         cls.initialize_items(scenario, MESSAGE_ITEMS)
 
+    @staticmethod
+    def enforce(scenario):
+        """Enforce data consistency in `scenario`."""
+        # Implemented in MESSAGE sub-class, below
+        # TODO make this an optional method of the ixmp.model.base.Model abstract class
+        pass
+
     def __init__(self, name=None, **model_options):
         # Update the default options with any user-provided options
         model_options.setdefault("model_dir", config.get("message model dir"))
@@ -350,6 +358,10 @@ def run(self, scenario):
            simultaneously; but using *different* CPLEX options in each process may
            produce unexpected results.
         """
+        # Ensure the data in `scenario` is consistent with the MESSAGE formulation
+        # TODO move this call to ixmp.model.base.Model.run(); remove here
+        self.enforce(scenario)
+
         # If two runs are kicked off simulatenously with the same self.model_dir, then
         # they will try to write the same optfile, and may write different contents.
         #
@@ -380,6 +392,32 @@ class MESSAGE(GAMSModel):
 
     name = "MESSAGE"
 
+    @staticmethod
+    def enforce(scenario):
+        """Enforce data consistency in `scenario`."""
+        # Check masks ("mapping sets") that indicate which elements of corresponding
+        # parameters are active/non-zero. Note that there are other masks currently
+        # handled in JDBCBackend. For the moment, this code does not backstop that
+        # behaviour.
+        # TODO Extend to handle all masks, e.g. for new backends.
+        for par_name in ("capacity_factor",):
+            # Name of the corresponding set
+            set_name = f"is_{par_name}"
+
+            # Existing and expected contents
+            existing = scenario.set(set_name)
+            expected = scenario.par(par_name).drop(columns=["value", "unit"])
+
+            if existing.equals(expected):
+                continue  # Contents are as expected; do nothing
+            # else:
+            #     scenario.add_set(set_name, expected)
+
+            # Not consistent; empty and then re-populate the set
+            with scenario.transact(f"Enforce consistency of {set_name} and {par_name}"):
+                scenario.remove_set(set_name, existing)
+                scenario.add_set(set_name, expected)
+
 
 class MACRO(GAMSModel):
     """Model class for MACRO."""

message_ix/testing/__init__.py

@@ -587,3 +587,143 @@ def make_westeros(mp, emissions=False, solve=False, quiet=True):
         scen.solve(quiet=quiet)
 
     return scen
+
+
+def make_subannual(
+    request,
+    tec_dict,
+    time_steps,
+    demand,
+    time_relative=[],
+    com_dict={"gas_ppl": {"input": "fuel", "output": "electr"}},
+    capacity={"gas_ppl": {"inv_cost": 0.1, "technical_lifetime": 5}},
+    capacity_factor={},
+    var_cost={},
+):
+    """Return an :class:`message_ix.Scenario` with subannual time resolution.
+
+    The scenario contains a simple model with two technologies, and a number of time
+    slices.
+
+    Parameters
+    ----------
+    request :
+        The pytest ``request`` fixture.
+    tec_dict : dict
+        A dictionary for a technology and required info for time-related parameters.
+        (e.g., ``tec_dict = {"gas_ppl": {"time_origin": ["summer"], "time": ["summer"],
+        "time_dest": ["summer"]}``)
+    time_steps : list of tuples
+        Information about each time slice, packed in a tuple with four elements,
+        including: time slice name, duration relative to "year", "temporal_lvl",
+        and parent time slice (e.g., ``time_steps = [("summer", 1, "season", "year")]``)
+    demand : dict
+        A dictionary for information of "demand" in each time slice.
+        (e.g., 11demand = {"summer": 2.5}``)
+    time_relative: list of str, optional
+        List of parent "time" slices, for which a relative duration time is maintained.
+        This will be used to specify parameter "duration_time_rel" for these "time"s.
+    com_dict : dict, optional
+        A dictionary for specifying "input" and "output" commodities.
+        (e.g., ``com_dict = {"gas_ppl": {"input": "fuel", "output": "electr"}}``)
+    capacity : dict, optional
+        Data for "inv_cost" and "technical_lifetime" per technology.
+    capacity_factor : dict, optional
+        "capacity_factor" with technology as key and "time"/"value" pairs as value.
+    var_cost : dict, optional
+        "var_cost" with technology as key and "time"/"value" pairs as value.
+    """
+    # Get the `test_mp` fixture for the requesting test function
+    mp = request.getfixturevalue("test_mp")
+
+    # Build an empty scenario
+    scen = Scenario(mp, request.node.name, scenario="test", version="new")
+
+    # Add required sets
+    scen.add_set("node", "node")
+    for c in com_dict.values():
+        scen.add_set("commodity", [x for x in list(c.values()) if x])
+
+    # Fixed values
+    y = 2020
+    unit = "GWa"
+
+    scen.add_set("level", "level")
+    scen.add_set("year", y)
+    scen.add_set("type_year", y)
+    scen.add_set("mode", "mode")
+
+    scen.add_set("technology", list(tec_dict.keys()))
+
+    # Add "time" and "duration_time" to the model
+    for (h, dur, tmp_lvl, parent) in time_steps:
+        scen.add_set("time", h)
+        scen.add_set("time", parent)
+        scen.add_set("lvl_temporal", tmp_lvl)
+        scen.add_set("map_temporal_hierarchy", [tmp_lvl, h, parent])
+        scen.add_par("duration_time", [h], dur, "-")
+
+    scen.add_set("time_relative", time_relative)
+
+    # Common dimensions for parameter data
+    common = dict(
+        node="node",
+        node_loc="node",
+        mode="mode",
+        level="level",
+        year=y,
+        year_vtg=y,
+        year_act=y,
+    )
+
+    # Define demand; unpack (key, value) pairs into individual pd.DataFrame rows
+    df = make_df(
+        "demand",
+        **common,
+        commodity="electr",
+        time=demand.keys(),
+        value=demand.values(),
+        unit=unit,
+    )
+    scen.add_par("demand", df)
+
+    # Add "input" and "output" parameters of technologies
+    common.update(value=1.0, unit="-")
+    base_output = make_df("output", **common, node_dest="node")
+    base_input = make_df("input", **common, node_origin="node")
+    for tec, times in tec_dict.items():
+        c = com_dict[tec]
+        for h1, h2 in zip(times["time"], times.get("time_dest", [])):
+            scen.add_par(
+                "output",
+                base_output.assign(
+                    technology=tec, commodity=c["output"], time=h1, time_dest=h2
+                ),
+            )
+        for h1, h2 in zip(times["time"], times.get("time_origin", [])):
+            scen.add_par(
+                "input",
+                base_input.assign(
+                    technology=tec, commodity=c["input"], time=h1, time_origin=h2
+                ),
+            )
+
+    # Add capacity related parameters
+    for year, tec in product([y], capacity.keys()):
+        for parname, val in capacity[tec].items():
+            scen.add_par(parname, ["node", tec, year], val, "-")
+
+    common.pop("value")
+
+    # Add capacity factor and variable cost data, both optional
+    for name, arg in [("capacity_factor", capacity_factor), ("var_cost", var_cost)]:
+        for tec, data in arg.items():
+            df = make_df(
+                name, **common, technology=tec, time=data.keys(), value=data.values()
+            )
+            scen.add_par(name, df)
+
+    scen.commit(f"Scenario with subannual time resolution for {request.node.name}")
+    scen.solve()
+
+    return scen

message_ix/tests/test_feature_capacity_factor.py

@@ -0,0 +1,156 @@
+"""Test ``capacity_factor`` effects, mainly for models with sub-annual resolution."""
+import pytest
+
+from message_ix import ModelError, Scenario
+from message_ix.testing import make_subannual
+
+
+def check_solution(scen: Scenario) -> None:
+    """Perform several assertions about the solution of `scen`."""
+    # Reading "ACT" and "CAP" from the solution
+    act = scen.var("ACT").set_index(["technology", "time"])
+    cap = scen.var("CAP").set_index(["technology"])
+
+    # 1) ACT is zero when capacity factor is zero
+    cf = scen.par("capacity_factor").set_index(["technology", "time"])
+    cf_zero = cf.loc[cf["value"] == 0]
+    for i in cf_zero.index:
+        assert act.loc[i, "lvl"] == 0
+
+    # 2) CAP is correctly calculated based on ACT and capacity_factor
+    for i in act.loc[act["lvl"] > 0].index:
+        # Correct ACT based on duration of each time slice
+        duration = float(scen.par("duration_time", {"time": i[1]})["value"])
+        act.loc[i, "duration-corrected"] = act.loc[i, "lvl"] / duration
+        # Divide by (non-zero) capacity factor
+        act.loc[i, "cf-corrected"] = act.loc[i, "duration-corrected"] / float(
+            cf.loc[i, "value"]
+        )
+    act = act.fillna(0).reset_index().set_index(["technology"])
+    # CAP = max("ACT" / "duration_time" / "capcity_factor")
+    for i in cap.index:
+        assert max(act.loc[i, "cf-corrected"]) == float(cap.loc[i, "lvl"])
+
+
+TD_0 = {
+    "gas_ppl": {
+        "time_origin": [],
+        "time": ["summer", "winter"],
+        "time_dest": ["summer", "winter"],
+    },
+}
+
+
+def test_capacity_factor_time(request):
+    """``capacity_factor`` is calculated correctly when it varies by time slice."""
+    # Build model and solve
+    scen = make_subannual(
+        request,
+        TD_0,
+        time_steps=[
+            ("summer", 0.5, "season", "year"),
+            ("winter", 0.5, "season", "year"),
+        ],
+        demand={"summer": 2, "winter": 1},
+        capacity_factor={"gas_ppl": {"summer": 0.8, "winter": 0.6}},
+        var_cost={"gas_ppl": {"summer": 0.2, "winter": 0.2}},
+    )
+    check_solution(scen)
+
+
+def test_capacity_factor_unequal_time(request):
+    """``capacity_factor`` is calculated correctly when ``duration_time`` is uneven."""
+    # Build model and solve
+    scen = make_subannual(
+        request,
+        TD_0,
+        time_steps=[
+            ("summer", 0.3, "season", "year"),
+            ("winter", 0.7, "season", "year"),
+        ],
+        demand={"summer": 2, "winter": 1},
+        capacity_factor={"gas_ppl": {"summer": 0.8, "winter": 0.8}},
+        var_cost={"gas_ppl": {"summer": 0.2, "winter": 0.2}},
+    )
+    check_solution(scen)
+
+
+TS_0 = [
+    ("day", 0.5, "subannual", "year"),
+    ("night", 0.5, "subannual", "year"),
+]
+
+
+def test_capacity_factor_zero(request):
+    """Test zero capacity factor (CF) in a time slice.
+
+    "solar_pv_ppl" is active in "day" and NOT at "night" (CF = 0). It is expected that
+    the model will be infeasible, because "demand" at night cannot be met.
+    """
+    # Technology input/output
+    tec_dict = {
+        "solar_pv_ppl": {
+            "time_origin": [],
+            "time": ["day", "night"],
+            "time_dest": ["day", "night"],
+        },
+    }
+
+    # Build model and solve (should raise GAMS error)
+    with pytest.raises(ModelError):
+        make_subannual(
+            request,
+            tec_dict,
+            com_dict={"solar_pv_ppl": {"input": "fuel", "output": "electr"}},
+            time_steps=TS_0,
+            demand={"day": 2, "night": 1},
+            capacity={"solar_pv_ppl": {"inv_cost": 0.2, "technical_lifetime": 5}},
+            capacity_factor={"solar_pv_ppl": {"day": 0.8, "night": 0}},
+        )
+
+
+def test_capacity_factor_zero_two(request):
+    """Test zero capacity factor (CF) in a time slice.
+
+    "solar_pv_ppl" is active in "day" and NOT at "night" (CF = 0). The model output
+    should show no activity of "solar_pv_ppl" at "night". So, "gas_ppl" is active at
+    "night", even though a more expensive technology.
+    """
+    # Technology input/output
+    tec_dict = {
+        "solar_pv_ppl": {
+            "time_origin": [],
+            "time": ["day", "night"],
+            "time_dest": ["day", "night"],
+        },
+        "gas_ppl": {
+            "time_origin": [],
+            "time": ["day", "night"],
+            "time_dest": ["day", "night"],
+        },
+    }
+
+    # Build model and solve
+    scen = make_subannual(
+        request,
+        tec_dict,
+        com_dict={
+            "solar_pv_ppl": {"input": "fuel", "output": "electr"},
+            "gas_ppl": {"input": "fuel", "output": "electr"},
+        },
+        time_steps=TS_0,
+        demand={"day": 2, "night": 1},
+        capacity={
+            "solar_pv_ppl": {"inv_cost": 0.1, "technical_lifetime": 5},
+            "gas_ppl": {"inv_cost": 0.1, "technical_lifetime": 5},
+        },
+        capacity_factor={
+            "solar_pv_ppl": {"day": 0.8, "night": 0},
+            "gas_ppl": {"day": 0.8, "night": 0.8},
+        },
+        var_cost={
+            "solar_pv_ppl": {"day": 0, "night": 0},
+            "gas_ppl": {"day": 0.2, "night": 0.2},
+        },
+    )
+    check_solution(scen)