Fixing materials flows in myopic mode (take 2)

temoa/_internal/table_data_puller.py

@@ -36,12 +36,12 @@ def _marks(num: int) -> str:
     return marks
 
 
-def ritvo(fi: FI) -> tuple[Region, Commodity, Technology, Vintage, Commodity]:
+def ritvo(fi: FI) -> tuple[Region, Commodity | None, Technology, Vintage, Commodity | None]:
     """convert FI to ritvo index"""
     return fi.r, fi.i, fi.t, fi.v, fi.o
 
 
-def rpetv(fi: FI, e: Commodity) -> tuple[Region, Period, Commodity, Technology, Vintage]:
+def rpetv(fi: FI, e: Commodity) -> tuple[Region, Period, Commodity | None, Technology, Vintage]:
     """convert FI and emission to rpetv index"""
     return fi.r, fi.p, e, fi.t, fi.v
 
@@ -58,7 +58,18 @@ def poll_capacity_results(model: TemoaModel, epsilon: float = 1e-5) -> CapData:
     for r, t, v in model.v_new_capacity.keys():
         if v in model.time_optimize:
             val = value(model.v_new_capacity[r, t, v])
-            if abs(val) < epsilon:
+            if val < -epsilon:
+                logger.warning(
+                    'Negative built capacity for %s, %s, %s: %s. '
+                    'This should not be possible. Could be a result of '
+                    'numerical instability or a code problem.',
+                    r,
+                    t,
+                    v,
+                    val,
+                )
+                continue
+            if val < epsilon:
                 continue
             new_cap = (r, t, v, val)
             built.append(new_cap)
@@ -67,7 +78,19 @@ def poll_capacity_results(model: TemoaModel, epsilon: float = 1e-5) -> CapData:
     net = []
     for r, p, t, v in model.v_capacity.keys():
         val = value(model.v_capacity[r, p, t, v])
-        if abs(val) < epsilon:
+        if val < -epsilon:
+            logger.warning(
+                'Negative net capacity for %s, %s, %s, %s: %s. '
+                'This should not be possible. Could be a result of '
+                'numerical instability or a code problem.',
+                r,
+                p,
+                t,
+                v,
+                val,
+            )
+            continue
+        if val < epsilon:
             continue
         new_net_cap = (r, p, t, v, val)
         net.append(new_net_cap)
@@ -84,8 +107,19 @@ def poll_capacity_results(model: TemoaModel, epsilon: float = 1e-5) -> CapData:
             if t in model.tech_retirement and v < p <= v + lifetime - value(model.period_length[p]):
                 early = value(model.v_retired_capacity[r, p, t, v])
                 eol -= early
-            early = 0 if abs(early) < epsilon else early
-            eol = 0 if abs(eol) < epsilon else eol
+            if early < -epsilon or eol < -epsilon:
+                logger.warning(
+                    'Negative retirement components for %s, %s, %s, %s: cap_eol=%s, cap_early=%s',
+                    r,
+                    p,
+                    t,
+                    v,
+                    eol,
+                    early,
+                )
+                continue
+            early = 0 if early < epsilon else early
+            eol = 0 if eol < epsilon else eol
             if early == 0 and eol == 0:
                 continue
             new_retired_cap = (r, p, t, v, eol, early)
@@ -190,7 +224,7 @@ def poll_flow_results(model: TemoaModel, epsilon: float = 1e-5) -> dict[FI, dict
         )
         for s in model.time_season[v]:
             for d in model.time_of_day:
-                fi = FI(r, v, s, d, i, t, v, cast('Commodity', 'construction_input'))
+                fi = FI(r, v, s, d, i, t, v, cast('Commodity', None))
                 flow = annual * value(model.segment_fraction[v, s, d])
                 if abs(flow) < epsilon:
                     continue
@@ -206,7 +240,7 @@ def poll_flow_results(model: TemoaModel, epsilon: float = 1e-5) -> dict[FI, dict
             )
             for s in model.time_season[p]:
                 for d in model.time_of_day:
-                    fi = FI(r, p, s, d, cast('Commodity', 'end_of_life_output'), t, v, o)
+                    fi = FI(r, p, s, d, cast('Commodity', None), t, v, o)
                     flow = annual * value(model.segment_fraction[p, s, d])
                     if abs(flow) < epsilon:
                         continue

temoa/_internal/table_writer.py

@@ -597,7 +597,9 @@ def write_flow_tables(self, iteration: int | None = None) -> None:
 
         self.connection.commit()
 
-    def _get_flow_units(self, flow_type: FlowType, input_comm: str, output_comm: str) -> str | None:
+    def _get_flow_units(
+        self, flow_type: FlowType, input_comm: str | None, output_comm: str | None
+    ) -> str | None:
         """
         Get units for flow based on flow type.
 
@@ -616,11 +618,12 @@ def _get_flow_units(self, flow_type: FlowType, input_comm: str, output_comm: str
         if not unit_prop:
             return None
 
-        if flow_type == FlowType.IN:
+        if flow_type == FlowType.IN and input_comm is not None:
             return unit_prop.get_flow_in_units(input_comm)
-        else:
+        elif output_comm is not None:
             # OUT, CURTAIL, FLEX all use output commodity units
             return unit_prop.get_flow_out_units(output_comm)
+        return None
 
     def write_summary_flow(self, model: TemoaModel, iteration: int | None = None) -> None:
         flow_data = self.calculate_flows(model=model)
@@ -636,9 +639,9 @@ def _insert_summary_flow_results(
         self.flow_register = flow_data
 
         # Aggregate flows (sum across seasons/time of day)
-        output_flows: defaultdict[tuple[str, Period, str, Technology, Vintage, str], float] = (
-            defaultdict(float)
-        )
+        output_flows: defaultdict[
+            tuple[str, Period, str | None, Technology, Vintage, str | None], float
+        ] = defaultdict(float)
 
         for fi, flows in self.flow_register.items():
             val = flows.get(FlowType.OUT)
@@ -676,7 +679,7 @@ def check_flow_balance(self, model: TemoaModel) -> bool:
         for fi, flow_vals in flows.items():
             if fi.t in model.tech_storage:
                 continue
-            if fi.i == 'end_of_life_output' or fi.o == 'construction_input':
+            if fi.i is None or fi.o is None:
                 continue
 
             fin = flow_vals.get(FlowType.IN, 0)

temoa/components/capacity.py

@@ -604,29 +604,54 @@ def create_capacity_and_retirement_sets(model: TemoaModel) -> None:
 
     logger.debug('Creating capacity, retirement, and construction/EOL sets.')
     # Calculate retirement periods based on lifetime and survival curves
-    for r, _i, t, v, _o in model.efficiency.sparse_iterkeys():
+    unique_rtv = {(r, t, v) for r, _i, t, v, _o in model.efficiency.sparse_iterkeys()} | set(
+        model.existing_capacity.sparse_iterkeys()
+    )
+    for r, t, v in unique_rtv:
+        if t in model.tech_uncap:
+            # No capacity to retire
+            continue
+        if t not in model.tech_all:
+            # Not an active technology so wont have a lifetime
+            # If it has an EOLoutput it will be in tech_all
+            continue
         lifetime = value(model.lifetime_process[r, t, v])
         for p in model.time_optimize:
+            # retires bang on start of horizon or survives into planning periods
+            is_p0_eol = (
+                (p == model.time_optimize.first())
+                and (v + lifetime == p)
+                and value(model.existing_capacity[r, t, v]) > 0
+            )
+            is_living = (r, t, v) in model.process_periods
+            if not (is_p0_eol or is_living):
+                continue
+
             is_natural_eol = p <= v + lifetime < p + value(model.period_length[p])
             is_early_retire = t in model.tech_retirement and v < p <= v + lifetime - value(
                 model.period_length[p]
             )
             is_survival_curve = model.is_survival_curve_process[r, t, v] and v <= p <= v + lifetime
 
-            if t not in model.tech_uncap and any(
-                (is_natural_eol, is_early_retire, is_survival_curve)
-            ):
+            if any((is_natural_eol, is_early_retire, is_survival_curve)):
                 model.retirement_periods.setdefault((r, t, v), set()).add(p)
 
     # Link construction materials to technologies
     for r, i, t, v in model.construction_input.sparse_iterkeys():
         model.capacity_consumption_techs.setdefault((r, v, i), set()).add(t)
+        model.used_techs.add(t)
 
     # Link end-of-life materials to retiring technologies
     for r, t, v, o in model.end_of_life_output.sparse_iterkeys():
         if (r, t, v) in model.retirement_periods:
             for p in model.retirement_periods[r, t, v]:
                 model.retirement_production_processes.setdefault((r, p, o), set()).add((t, v))
+                model.used_techs.add(t)
+
+    # Link end-of-life emissions to retiring technologies
+    for r, _e, t, v in model.emission_end_of_life.sparse_iterkeys():
+        if (r, t, v) in model.retirement_periods:
+            model.used_techs.add(t)
 
     # Create active capacity index sets from the now-populated process_vintages
     model.new_capacity_rtv = {

temoa/components/flows.py

@@ -110,22 +110,24 @@ def create_commodity_balance_and_flow_sets(model: TemoaModel) -> None:
     """
     logger.debug('Creating commodity balance and active flow index sets.')
     # 1. Commodity Balance
-    commodity_upstream_rpo = set(
+    commodity_upstream_rpo = {
         (r, p, o)
         for r, p, o in (
             model.commodity_up_stream_process
-            | model.retirement_production_processes | model.import_regions
+            | model.retirement_production_processes
+            | model.import_regions
         )
-        if o not in model.commodity_sink # only balanced if input to another process (caught below)
-    )
-    commodity_downstream_rpi = set(
+        if o not in model.commodity_sink  # only balanced if input to another process (caught below)
+    }
+    commodity_downstream_rpi = {
         (r, p, i)
         for r, p, i in (
             model.commodity_down_stream_process
-            | model.capacity_consumption_techs | model.export_regions
+            | model.capacity_consumption_techs
+            | model.export_regions
         )
-        if i not in model.commodity_source # sources are never balanced (infinite source)
-    )
+        if i not in model.commodity_source  # sources are never balanced (infinite source)
+    }
     model.commodity_balance_rpc = commodity_upstream_rpo.union(commodity_downstream_rpi)
 
     # 2. Active Flow Indices (Time-Sliced)

temoa/components/technology.py

@@ -61,6 +61,7 @@ def lifetime_process_indices(model: TemoaModel) -> set[tuple[Region, Technology,
     process indices that may be specified in the lifetime_process parameter.
     """
     indices = {(r, t, v) for r, i, t, v, o in model.efficiency.sparse_iterkeys()}
+    indices = indices | set(model.existing_capacity.sparse_iterkeys())
 
     return indices
 
@@ -159,13 +160,13 @@ def populate_core_dictionaries(model: TemoaModel) -> None:
                     process,
                 )
                 continue
-            if t not in model.tech_uncap and model.existing_capacity[process] == 0:
-                logger.warning(
-                    'Notice: Unnecessary specification of existing_capacity for %s. '
-                    'Declaring a capacity of zero may be omitted.',
-                    process,
+            if t not in model.tech_uncap and value(model.existing_capacity[process]) <= 0:
+                msg = (
+                    f'Notice: Non-positive existing capacity declaration {process}. '
+                    'This is not supported for processes surviving into future periods.'
                 )
-                continue
+                logger.error(msg)
+                raise ValueError(msg)
             if v + lifetime <= first_period:
                 logger.info(
                     '%s specified as existing_capacity, but its '
@@ -227,6 +228,8 @@ def create_survival_curve(model: TemoaModel) -> None:
 
     for r, _, t, v, _ in model.efficiency.sparse_iterkeys():
         model.is_survival_curve_process[r, t, v] = False  # by default
+    for r, t, v in model.existing_capacity.sparse_iterkeys():
+        model.is_survival_curve_process[r, t, v] = False  # by default
 
     # Collect rptv indices into (r, t, v): p dictionary
     for r, p, t, v in model.lifetime_survival_curve.sparse_iterkeys():
@@ -352,7 +355,7 @@ def check_efficiency_indices(model: TemoaModel) -> None:
         f_msg = msg.format(', '.join(diff_str))
         logger.error(f_msg)
         raise ValueError(f_msg)
-    
+
     c_inputs = {i for r, i, t, v, o in model.efficiency.sparse_iterkeys()}
     c_inputs = c_inputs | {i for r, i, t, v in model.construction_input.sparse_iterkeys()}
     c_carrier = c_inputs | c_outputs
@@ -368,8 +371,11 @@ def check_efficiency_indices(model: TemoaModel) -> None:
         f_msg = msg.format(', '.join(symdiff_str))
         logger.error(f_msg)
         raise ValueError(f_msg)
-    
+
     techs = {t for r, i, t, v, o in model.efficiency.sparse_iterkeys()}
+    techs = techs | {t for r, t, v, o in model.end_of_life_output.sparse_iterkeys()}
+    techs = techs | {t for r, i, t, v in model.construction_input.sparse_iterkeys()}
+    techs = techs | {t for r, e, t, v in model.emission_end_of_life.sparse_iterkeys()}
 
     symdiff = techs.symmetric_difference(model.tech_production)
     if symdiff:
@@ -431,3 +437,29 @@ def check_efficiency_variable(model: TemoaModel) -> None:
                     num_seg,
                     (r, p, i, t, v, o),
                 )
+
+
+def check_existing_capacity(model: TemoaModel) -> None:
+    """
+    Check that all existing capacities are properly accounted for in the model.
+    """
+    for r, t, v in model.existing_capacity.sparse_iterkeys():
+        cap = value(model.existing_capacity[r, t, v])
+        if cap <= 0:
+            msg = (
+                f'Existing capacity {r, t, v} has non-positive capacity {cap}. '
+                'This entry will be ignored.'
+            )
+            logger.warning(msg)
+            continue
+        if t not in model.tech_all:
+            continue
+        life = value(model.lifetime_process[r, t, v])
+        if (r, t, v) not in model.process_periods and v + life > model.time_optimize.first():
+            msg = (
+                f'Existing capacity {r, t, v} with lifetime {life} and capacity {cap} '
+                'should extend into future periods but it is not in process periods. '
+                'Was it included in the Efficiency table?'
+            )
+            logger.error(msg)
+            raise ValueError(msg)

temoa/core/model.py

@@ -83,7 +83,7 @@ def create_sparse_dicts(model: 'TemoaModel') -> None:
         for tech in sorted(unused_techs):
             logger.warning(
                 "Notice: '%s' is specified as a technology but is not "
-                'utilized in the efficiency parameter.',
+                'utilized in the process network.',
                 tech,
             )
 
@@ -388,6 +388,9 @@ def __init__(self, *args: object, **kwargs: object) -> None:
         self.end_of_life_output = Param(
             self.regions, self.tech_with_capacity, self.vintage_all, self.commodity_carrier
         )
+        self.emission_end_of_life = Param(
+            self.regions, self.commodity_emissions, self.tech_with_capacity, self.vintage_all
+        )
 
         self.efficiency = Param(
             self.regional_indices,
@@ -511,6 +514,7 @@ def __init__(self, *args: object, **kwargs: object) -> None:
         # equations below.
         self.create_sparse_dicts = BuildAction(rule=create_sparse_dicts)
         self.initialize_demands = BuildAction(rule=commodities.create_demands)
+        self.validate_existing_capacity = BuildAction(rule=technology.check_existing_capacity)
 
         self.capacity_factor_rpsdt = Set(dimen=5, initialize=capacity.capacity_factor_tech_indices)
         self.capacity_factor_tech = Param(
@@ -725,9 +729,6 @@ def __init__(self, *args: object, **kwargs: object) -> None:
             self.tech_with_capacity,
             self.vintage_optimize,
         )
-        self.emission_end_of_life = Param(
-            self.regions, self.commodity_emissions, self.tech_with_capacity, self.vintage_all
-        )
 
         self.myopic_discounting_year = Param(default=0)
 

temoa/data_io/component_manifest.py

@@ -556,6 +556,7 @@ def build_manifest(model: TemoaModel) -> list[LoadItem]:
             columns=['region', 'emis_comm', 'tech', 'vintage', 'value'],
             validator_name='viable_rtv',
             validation_map=(0, 2, 3),
+            custom_loader_name='_load_emission_embodied',
             is_period_filtered=False,
             is_table_required=False,
         ),
@@ -574,14 +575,15 @@ def build_manifest(model: TemoaModel) -> list[LoadItem]:
             columns=['region', 'input_comm', 'tech', 'vintage', 'value'],
             validator_name='viable_rtv',
             validation_map=(0, 2, 3),
+            custom_loader_name='_load_construction_input',
             is_period_filtered=False,
             is_table_required=False,
         ),
         LoadItem(
             component=model.end_of_life_output,
             table='end_of_life_output',
             columns=['region', 'tech', 'vintage', 'output_comm', 'value'],
-            validator_name='viable_rtv',
+            validator_name='viable_rtv_eol',
             validation_map=(0, 1, 2),
             is_period_filtered=False,
             is_table_required=False,