Merge PR #61 and migrate ethylene brownfield capacity

CHANGES.md

@@ -5,10 +5,12 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
-## [Unreleased]
-* Added foreground_db_name argument to LCAConfig
-
-## [0.4.2] - 2026-04-27
+## [Unreleased]
+* Added foreground_db_name argument to LCAConfig
+* Fixed under-counting of installation impacts: one installed unit now delivers its production over its whole lifetime, not in every operating year (changes results for processes with multi-year operation windows)
+* Postprocessing capacity outputs (`get_production_capacity`, `plot_capacity_balance`, `plot_utilization_heatmap`) now report annual capacity, comparable with production
+
+## [0.4.2] - 2026-04-27
 * Fix conda publish workflow
 
 ## [0.4.1] - 2026-04-27

docs/content/constraints.md

@@ -35,7 +35,10 @@ solved_model, objective, results = optimizer.solve_model(model)
 
 ## Process Deployment Limits
 
-Control how much capacity can be installed for each process.
+Control how many units of a process can be installed.
+
+!!! note "Limits are in process units"
+    Deployment limits bound `var_installation`, which counts process **units**. One unit delivers the production of its production exchange over its whole lifetime, so a limit of 50 allows a lifetime output of 50 reference-flow units, not an annual output of 50. See [Foreground Modeling](foreground_modeling.md#what-one-installed-unit-means).
 
 ### Time-Specific Deployment Limits
 
@@ -100,6 +103,8 @@ Control how much a process can operate in each time period.
 !!! note "Operation vs Installation"
     Operation limits constrain how much of the installed capacity is actually used in each period. A process can only operate up to its installed capacity, but operation limits can further restrict this.
 
+    Both variables are counted in process units: `var_operation[p, v, t]` is how many units of vintage `v` run in year `t`, bounded by `var_installation[p, v]`. An operation limit therefore caps the units running in that year, summed over all vintages — multiply by the process's output per unit and year to express it as a production limit (or use [flow limits](#flow-limits) to bound production directly).
+
 **Example: Limit coal plant operation during phase-out**
 ```python
 model_inputs.process_operation_limits_max = {
@@ -266,7 +271,9 @@ Model systems with pre-existing infrastructure that was installed before the opt
 
 | Field | Type | Description |
 |-------|------|-------------|
-| `existing_capacity` | `Dict[Tuple[str, int], float]` | Capacity at `(process, installation_year)` |
+| `existing_capacity` | `Dict[Tuple[str, int], float]` | Units at `(process, installation_year)` |
+
+Existing capacity is given in the same process units as `var_installation`, and the entries behave like vintages of that installation year: they can run while their lifecycle stage falls inside the operation window, and they carry no installation impacts.
 
 !!! info "Brownfield vs Greenfield"
     - **Greenfield**: Optimization starts from scratch with no existing capacity

docs/content/foreground_modeling.md

@@ -132,9 +132,37 @@ TemporalDistribution(
 
 ---
 
+## What One Installed Unit Means
+
+This convention decides how installation impacts are amortized, so it is worth stating explicitly:
+
+**One unit of a process delivers the production stated by its production exchange — over its entire lifetime, not per year.**
+
+For the example above (`amount=1`, spread as `[0, 0.5, 0.5, 0]` over a two-year operation window):
+
+| Quantity | Value |
+|----------|-------|
+| Lifetime output of one unit | 1 kg (the sum of the temporal distribution) |
+| Output of one running unit per year | 0.5 kg (the per-τ entry) |
+| Units needed to deliver 10 kg in a single year | 20 |
+| Units needed to deliver 10 kg in each of two consecutive years | 20 (one cohort, fully used) |
+
+Consequences for the optimization:
+
+- `var_installation[p, v]` counts **units**, and installation-dependent exchanges (construction, end-of-life) are incurred once per unit.
+- `var_operation[p, v, t]` counts how many of those units **run in year t**, bounded by the units installed: `var_operation[p, v, t] ≤ var_installation[p, v]`.
+- Annual output is `production[τ = t - v] × var_operation[p, v, t]`.
+
+This is what makes optimex agree with a standard LCA of the delivered amount whenever every unit is fully utilized. If the demand profile does not allow full utilization — for instance a single demand year for a process with a multi-year operation window — the model builds capacity that is partly idle, and the resulting impact is legitimately *higher* than the standard LCA score of the delivered amount.
+
+!!! tip "Reporting an annual capacity"
+    Because installed units cover a whole lifetime, `PostProcessor.get_installation()` is not comparable with per-year production. Use `PostProcessor.get_production_capacity()`, which multiplies unit counts by the output per unit and year.
+
+---
+
 ## The Operation Flag
 
-Exchanges marked with `"operation": True` **scale with the operational level** $\mathbf{o}_{t,v}$. This enables flexible operation where a process can run below its installed capacity $\mathbf{s}_v$ (see [Theory: Flexible Operation](theory.md#step-5-flexible-operation)).
+Exchanges marked with `"operation": True` **scale with the operational level** $\mathbf{o}_{t,v}$, the number of units running. This enables flexible operation where a process can run below its installed capacity $\mathbf{s}_v$ (see [Theory: Flexible Operation](theory.md#step-5-flexible-operation)).
 
 ```python
 "exchanges": [

docs/content/postprocessing_guide.md

@@ -75,7 +75,10 @@ df_installation = pp.get_installation()
 Returns a DataFrame with:
 - **Index**: System time (years)
 - **Columns**: Process IDs
-- **Values**: Capacity installed at that time
+- **Values**: Number of process **units** installed at that time
+
+!!! warning "Installed units are a lifetime quantity"
+    One installed unit delivers the production of its production exchange over its **whole lifetime** (the sum of its temporal distribution), spread across its operation window. These numbers are therefore not an annual capacity and must not be compared directly with production or demand. For the per-year quantity, use [`get_production_capacity()`](#production-capacity).
 
 ```
 Process     ProcessA    ProcessB
@@ -88,10 +91,10 @@ Time
 
 **Common analyses:**
 ```python
-# Total capacity installed per process
-total_capacity = df_installation.sum()
+# Total units installed per process
+total_units = df_installation.sum()
 
-# Cumulative capacity over time
+# Cumulative units installed over time
 cumulative = df_installation.cumsum()
 
 # When each process was first deployed
@@ -102,7 +105,7 @@ first_deployment = df_installation[df_installation > 0].idxmax()
 
 ### Operation Levels
 
-Get how much each process operates per time period:
+Get how many units of each process run per time period:
 
 ```python
 df_operation = pp.get_operation()
@@ -111,7 +114,7 @@ df_operation = pp.get_operation()
 Returns a DataFrame with:
 - **Index**: System time (years)
 - **Columns**: Process IDs
-- **Values**: Operation level (units operating)
+- **Values**: Number of units running in that year
 
 ```
 Process     ProcessA    ProcessB
@@ -124,11 +127,31 @@ Time
 
 **Capacity utilization:**
 ```python
-# Calculate utilization rate
-cumulative_capacity = df_installation.cumsum()
-utilization = df_operation / cumulative_capacity
+# Units running vs. units available: compare production against annual capacity
+df_capacity = pp.get_production_capacity()   # per year
+df_production = pp.get_production()          # per year
+utilization = df_production.groupby(level="Product", axis=1).sum() / df_capacity
 ```
 
+Do not divide `df_operation` by the cumulative sum of `df_installation`: units only run within their operation window, so a cumulative sum counts units that were never or are no longer available. `plot_utilization_heatmap()` does this correctly.
+
+---
+
+### Production Capacity
+
+Get the maximum **annual** production available in each year:
+
+```python
+df_capacity = pp.get_production_capacity()
+```
+
+Returns a DataFrame with:
+- **Index**: System time (years)
+- **Columns**: Products
+- **Values**: Output the installed and existing units could deliver in that year
+
+This multiplies the units of every vintage that is in its operation phase by the output that vintage yields per unit and year, so it is directly comparable with `get_production()` and `get_demand()`.
+
 ---
 
 ### Production
@@ -203,7 +226,7 @@ fig = pp.plot_installation()
 
 ### Operation Plot
 
-Shows operation levels over time:
+Shows how many units run over time:
 
 ```python
 fig = pp.plot_operation()
@@ -213,7 +236,7 @@ fig = pp.plot_operation()
 
 ### Capacity Balance Plot
 
-Compares actual production with maximum available capacity:
+Compares actual production with the maximum annual capacity available:
 
 ```python
 fig = pp.plot_capacity_balance()
@@ -232,12 +255,16 @@ For custom analysis, access the solved model directly:
 ```python
 import pyomo.environ as pyo
 
-# Get decision variable values
+# Get decision variable values: units installed per vintage, and units of each
+# vintage running in a given year
 for p in solved_model.PROCESS:
     for t in solved_model.SYSTEM_TIME:
         installation = pyo.value(solved_model.var_installation[p, t])
-        operation = pyo.value(solved_model.var_operation[p, t])
-        print(f"{p}, {t}: install={installation:.2f}, operate={operation:.2f}")
+        print(f"{p}, installed in {t}: {installation:.2f} units")
+
+for (p, v, t) in solved_model.ACTIVE_VINTAGE_TIME:
+    operation = pyo.value(solved_model.var_operation[p, v, t])
+    print(f"{p}, vintage {v}, year {t}: {operation:.2f} units running")
 
 # Access expressions
 for c in solved_model.CATEGORY:
@@ -315,7 +342,7 @@ for cat in df_impacts.columns.get_level_values(0).unique():
     total = df_impacts[cat].sum().sum()
     print(f"  {cat}: {total:.2f}")
 
-print(f"\nTotal capacity installed by process:")
+print(f"\nTotal units installed by process:")
 for proc in df_installation.columns:
     total = df_installation[proc].sum()
     print(f"  {proc}: {total:.2f}")

docs/content/quickstart.md

@@ -190,8 +190,9 @@ model_inputs.existing_capacity = {("old_plant", 2010): 500.0}
 | Method | Returns | Description |
 |--------|---------|-------------|
 | `get_impacts()` | DataFrame | Impact by process, category, time |
-| `get_installation()` | DataFrame | Capacity installed per time |
-| `get_operation()` | DataFrame | Operation level per time |
+| `get_installation()` | DataFrame | Process units installed per time (a lifetime quantity, not an annual capacity) |
+| `get_production_capacity()` | DataFrame | Maximum annual production available per time |
+| `get_operation()` | DataFrame | Units running per time |
 | `get_production()` | DataFrame | Production by process and product |
 | `get_demand()` | DataFrame | Demand fulfillment over time |
 | `plot_impacts()` | Figure | Stacked area plot of impacts |

docs/content/theory_how_it_works.md

@@ -45,8 +45,10 @@ Both the foreground and background systems can evolve over time:
 
 Traditional LCO assumes processes always run at full capacity. `optimex` separates the decision into two components:
 
-- **Capacity installation**: How much capacity is built at a given time (vintage)
-- **Operational level**: How much of that capacity is actually used at each point in time
+- **Capacity installation**: How many process units are built at a given time (vintage)
+- **Operational level**: How many of those units are actually running at each point in time
+
+One installed unit corresponds to one reference-flow unit of the underlying LCA process: it delivers the production stated by its production exchange over its **whole lifetime**, so the per-time-step entry of that exchange is what one running unit yields **per year**. Installation-dependent exchanges are charged once per unit, which means the construction burden of a unit is amortized over its lifetime output — no more and no less. A demand profile that leaves part of a unit's lifetime unused therefore shows up as idle capacity and a correspondingly higher impact.
 
 This separation is important because it enables **vintage-specific dispatch**: when multiple cohorts of the same technology coexist, the optimizer can preferentially utilize cleaner vintages — creating an emissions-aware merit order. It also allows the model to identify strategic overcapacities, where early investment in clean technologies offsets the stranded cost of idled fossil infrastructure.
 

notebooks/ethylene_case_study.ipynb

@@ -2180,28 +2180,50 @@
     ")\n",
     "functional_demand = {ethylene: ethylene_demand}\n",
     "\n",
+    "# PR #61 semantics: existing_capacity stores process units. One steam-cracker\n",
+    "# unit delivers the 1 kg production exchange over its 50-year lifetime.\n",
+    "steam_cracking_lifetime_years = int(LIFETIMES_YEARS[\"steam_cracking\"])\n",
+    "steam_cracking_annual_output_per_unit_kg = 1.0 / steam_cracking_lifetime_years\n",
+    "brownfield_annual_capacity_per_vintage_kg = float(\n",
+    "    assumptions.loc[\"brownfield.capacity_per_vintage\", \"value\"]\n",
+    ")\n",
+    "brownfield_units_per_vintage = (\n",
+    "    brownfield_annual_capacity_per_vintage_kg\n",
+    "    / steam_cracking_annual_output_per_unit_kg\n",
+    ")\n",
+    "\n",
     "existing_capacities = {\n",
     "    (\n",
     "        \"steam_cracking\",\n",
     "        int(assumptions.loc[\"brownfield.vintage_1\", \"value\"]),\n",
-    "    ): float(assumptions.loc[\"brownfield.capacity_per_vintage\", \"value\"]),\n",
+    "    ): brownfield_units_per_vintage,\n",
     "    (\n",
     "        \"steam_cracking\",\n",
     "        int(assumptions.loc[\"brownfield.vintage_2\", \"value\"]),\n",
-    "    ): float(assumptions.loc[\"brownfield.capacity_per_vintage\", \"value\"]),\n",
+    "    ): brownfield_units_per_vintage,\n",
     "}\n",
     "\n",
     "brownfield_table = pd.DataFrame(\n",
     "    [\n",
     "        {\n",
     "            \"process\": PROCESS_NAMES[process],\n",
     "            \"installation_year\": year,\n",
-    "            \"capacity_kg_per_year\": capacity,\n",
+    "            \"installed_process_units\": units,\n",
+    "            \"annual_capacity_kg_per_year\": (\n",
+    "                units * steam_cracking_annual_output_per_unit_kg\n",
+    "            ),\n",
     "        }\n",
-    "        for (process, year), capacity in existing_capacities.items()\n",
+    "        for (process, year), units in existing_capacities.items()\n",
     "    ]\n",
     ")\n",
-    "assert brownfield_table[\"capacity_kg_per_year\"].sum() == 1e9\n",
+    "assert np.allclose(\n",
+    "    brownfield_table[\"annual_capacity_kg_per_year\"],\n",
+    "    brownfield_annual_capacity_per_vintage_kg,\n",
+    ")\n",
+    "assert np.isclose(\n",
+    "    brownfield_table[\"annual_capacity_kg_per_year\"].sum(),\n",
+    "    2 * brownfield_annual_capacity_per_vintage_kg,\n",
+    ")\n",
     "brownfield_table\n"
    ]
   },

notebooks/plots/paper_figures.py

@@ -109,7 +109,15 @@ def load_scenario_data(scenario: str) -> dict:
 
 
 def clean_capacity_df(df: pd.DataFrame) -> pd.DataFrame:
-    """Clean capacity DataFrame and rename columns."""
+    """
+    Clean capacity DataFrame and rename columns.
+
+    The capacity files are exported from `PostProcessor.get_production_capacity()`,
+    which reports ANNUAL production capacity (installed units times the output per
+    unit and year). That is the quantity to plot against annual production - the raw
+    installation numbers from `get_installation()` are unit counts covering a unit's
+    whole lifetime and are not comparable with per-year production.
+    """
     df = df.copy()
     df.index = df.index.astype(int)
     df.columns = [PRODUCT_NAMES.get(c, c) for c in df.columns]
@@ -405,7 +413,7 @@ def create_combined_results_figure(scenarios_data: dict):
                 capacity_values = cap_df[product_col].values[cap_mask] / 1e6
                 if np.any(capacity_values > 0.001):
                     ax.plot(x_positions, capacity_values, color="#000000", linestyle="--",
-                            linewidth=1, marker="", label="Capacity", zorder=4)
+                            linewidth=1, marker="", label="Annual capacity", zorder=4)
 
             # if not is_intermediate:
             #     ax.axhline(y=1, color="red", linestyle="--", linewidth=1)
@@ -447,7 +455,7 @@ def create_combined_results_figure(scenarios_data: dict):
         all_handles.append(Patch(facecolor=color, edgecolor="white", linewidth=0.5))
         all_labels.append(process)
     all_handles.append(plt.Line2D([0], [0], color="#000000", linestyle="--", linewidth=1))
-    all_labels.append("Available capacity")
+    all_labels.append("Available annual capacity")
     # all_handles.append(plt.Line2D([0], [0], color="red", linestyle="--", linewidth=1))
     # all_labels.append("Demand")
     # all_handles.append(Patch(facecolor="#BDCD00", edgecolor="#41811C", linewidth=1, hatch="///"))
@@ -1213,7 +1221,7 @@ def create_combined_results_and_impacts_figure(scenarios_data: dict):
                 capacity_values = cap_df[product_col].values[cap_mask] / 1e6
                 if np.any(capacity_values > 0.001):
                     ax.plot(x_positions, capacity_values, color="#000000", linestyle="--",
-                            linewidth=1, marker="", label="Capacity", zorder=4)
+                            linewidth=1, marker="", label="Annual capacity", zorder=4)
 
             ax.axhline(y=0, color="gray", linewidth=0.5, zorder=0)
 
@@ -1449,7 +1457,7 @@ def add_row_group_label(start_row, end_row, label_text):
 
     # Add capacity line
     all_handles.append(plt.Line2D([0], [0], color="#000000", linestyle="--", linewidth=1))
-    all_labels_legend.append("Available capacity")
+    all_labels_legend.append("Available annual capacity")
 
     fig.legend(all_handles, all_labels_legend, loc="lower center", bbox_to_anchor=(0.5, 0.07),
                ncol=4, frameon=False, fontsize=9)