diff --git a/CHANGES.md b/CHANGES.md index 3879a96..91194f7 100644 --- a/CHANGES.md +++ b/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 diff --git a/docs/content/constraints.md b/docs/content/constraints.md index f3fec71..7c64f5d 100644 --- a/docs/content/constraints.md +++ b/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 diff --git a/docs/content/foreground_modeling.md b/docs/content/foreground_modeling.md index d16fbe7..b6a4730 100644 --- a/docs/content/foreground_modeling.md +++ b/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": [ diff --git a/docs/content/postprocessing_guide.md b/docs/content/postprocessing_guide.md index 56cbd56..77f202b 100644 --- a/docs/content/postprocessing_guide.md +++ b/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}") diff --git a/docs/content/quickstart.md b/docs/content/quickstart.md index f88467e..f97d8ee 100644 --- a/docs/content/quickstart.md +++ b/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 | diff --git a/docs/content/theory_how_it_works.md b/docs/content/theory_how_it_works.md index a8ddf47..d8a89ca 100644 --- a/docs/content/theory_how_it_works.md +++ b/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. diff --git a/notebooks/ethylene_case_study.ipynb b/notebooks/ethylene_case_study.ipynb index 0135d32..616d029 100644 --- a/notebooks/ethylene_case_study.ipynb +++ b/notebooks/ethylene_case_study.ipynb @@ -2180,15 +2180,27 @@ ")\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", @@ -2196,12 +2208,22 @@ " {\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" ] }, diff --git a/notebooks/plots/paper_figures.py b/notebooks/plots/paper_figures.py index 8ba53d6..00f9571 100644 --- a/notebooks/plots/paper_figures.py +++ b/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) diff --git a/src/optimex/converter.py b/src/optimex/converter.py index 8291492..b885b49 100644 --- a/src/optimex/converter.py +++ b/src/optimex/converter.py @@ -101,7 +101,9 @@ class OptimizationModelInputs(BaseModel): foreground_production: Dict[Tuple[str, str, int], float] = Field( ..., description=( - "Maps (process, product, process_time) to produced amount." + "Maps (process, product, process_time) to produced amount. The entry is " + "the output of one running unit in that year of its lifecycle; the sum " + "over the operation window is the lifetime output of one installed unit." ), ) @@ -248,16 +250,31 @@ class OptimizationModelInputs(BaseModel): ), ) process_deployment_limits_max: Optional[Dict[Tuple[str, int], float]] = Field( - None, description="Upper bounds on (process, system_time) deployment." + None, + description=( + "Upper bounds on (process, system_time) deployment, in process units. One " + "unit delivers its production temporal distribution over its whole lifetime." + ), ) process_deployment_limits_min: Optional[Dict[Tuple[str, int], float]] = Field( - None, description="Lower bounds on (process, system_time) deployment." + None, + description=( + "Lower bounds on (process, system_time) deployment, in process units." + ), ) process_operation_limits_max: Optional[Dict[Tuple[str, int], float]] = Field( - None, description="Upper bounds on (process, system_time) operation." + None, + description=( + "Upper bounds on (process, system_time) operation, i.e. on the number of " + "units running in that year, summed over all vintages." + ), ) process_operation_limits_min: Optional[Dict[Tuple[str, int], float]] = Field( - None, description="Lower bounds on (process, system_time) operation." + None, + description=( + "Lower bounds on (process, system_time) operation, i.e. on the number of " + "units running in that year, summed over all vintages." + ), ) cumulative_process_limits_max: Optional[Dict[str, float]] = Field( None, description=("Global upper bound on cumulative deployment for a process.") @@ -277,7 +294,8 @@ class OptimizationModelInputs(BaseModel): None, description=( "Existing (brownfield) capacity installed before the optimization horizon. " - "Maps (process, installation_year) to capacity amount. Installation years " + "Maps (process, installation_year) to a number of process units, the same " + "unit as var_installation. Installation years " "must be before min(SYSTEM_TIME). These capacities contribute to operation " "and production but their installation impacts are excluded (sunk costs)." ), diff --git a/src/optimex/optimizer.py b/src/optimex/optimizer.py index 2e53e10..1ddda40 100644 --- a/src/optimex/optimizer.py +++ b/src/optimex/optimizer.py @@ -4,13 +4,35 @@ This module creates and solves Pyomo optimization models that minimize environmental impacts over time while meeting demand constraints and respecting process limits. +## Unit Convention + +One "unit" of a process is one reference-flow unit of the underlying LCA process: +installing 1 unit delivers exactly the production stated by its production temporal +distribution, i.e. `sum_tau foreground_production[p, r, tau]`, spread over the whole +operation window. The per-tau entry is therefore the unit's output *per year of +operation*, and the sum over the window is its *lifetime* output. + +Consequences: +- `var_installation[p, v]` counts units of vintage v; its installation-dependent + flows (construction, end-of-life) are incurred once per unit. +- `var_operation[p, v, t]` counts how many of those units are running in year t, + so the capacity bound is `var_operation[p, v, t] <= var_installation[p, v]`. +- Annual output at time t of vintage v is + `foreground_production[p, r, t - v] * var_operation[p, v, t]`, never the sum over + the operation window (that would let one unit deliver its lifetime output every + year and under-count installation impacts by the number of operating years). +- To convert installed units into an annual production capacity for reporting or + plotting, multiply by the same per-tau rate — see + `PostProcessor.get_production_capacity()`. + ## Scaling Convention The optimization uses a two-tier scaling system for numerical stability: ### Decision Variables (REAL UNITS) - `var_installation[p, t]`: Number of process units installed (dimensionless) -- `var_operation[p, t]`: Operation level (dimensionless, 0 to capacity) +- `var_operation[p, v, t]`: Units of vintage v running at time t (dimensionless, + 0 to the units available from that vintage) Both decision variables remain in REAL (unscaled) units to: 1. Maintain physical interpretability @@ -52,12 +74,12 @@ Example constraint dimensional analysis: ``` ProductDemandFulfillment: - production [kg SCALED/operation] × var_operation [#] = demand [kg SCALED] ✓ + production[p, r, t-v] [kg SCALED/(unit·year)] × var_operation [# units running] + = demand [kg SCALED] ✓ -OperationLimit (LHS): - var_operation [#] × production [kg SCALED/operation] × fg_scale = [kg REAL] -OperationLimit (RHS): - production [kg SCALED/process] × fg_scale × var_installation [#] = [kg REAL] +OperationCapacity: + var_operation [# units running] <= var_installation [# units installed] ✓ + (both REAL unit counts, so no scaling factor appears) ``` """ @@ -617,37 +639,31 @@ def scale_tensor_by_operation(tensor: pyo.Param, flow_set: str, overrides: dict, Scale tensor by operation, summing flows across active vintages. With 3D var_operation[p, v, t], we sum flow contributions from each - vintage cohort operating at time t. + vintage cohort operating at time t. The flow taken from the tensor is the + one at the vintage's current lifecycle stage tau = t - v, i.e. the flow + *per operating unit and year* — not the sum over the whole operation + window, which is the unit's lifetime total. """ def expr(m, p, x, t): # Only apply operational scaling to flows marked as operational if pyo.value(m.operation_flow[p, x]) == 0: return 0 - op_start = pyo.value(m.process_operation_start[p]) - op_end = pyo.value(m.process_operation_end[p]) - total = 0 # Sum flows across all active vintages at time t for (proc, v, time) in m.ACTIVE_VINTAGE_TIME: if proc != p or time != t: continue - # Get vintage-specific flow rate + tau = t - v # lifecycle stage of this vintage in year t + if tau not in m.PROCESS_TIME: + continue + + # Get vintage-specific flow rate for this lifecycle stage if (p, x) in overrides_index: - # Vintage-aware: sum flow rates across all operating taus for this vintage - flow_rate = sum( - overrides.get((p, x, tau, v), tensor[p, x, tau]) - for tau in m.PROCESS_TIME - if op_start <= tau <= op_end - ) + flow_rate = overrides.get((p, x, tau, v), tensor[p, x, tau]) else: - # No overrides: all vintages have same flow rate - flow_rate = sum( - tensor[p, x, tau] - for tau in m.PROCESS_TIME - if op_start <= tau <= op_end - ) + flow_rate = tensor[p, x, tau] total += flow_rate * m.var_operation[p, v, t] @@ -738,60 +754,32 @@ def has_production_overrides(p, r): """Check if any vintage overrides exist for this process/product.""" return (p, r) in model._production_overrides_index - def operation_capacity_constraint_rule(model, p, v, t, r): + def operation_capacity_constraint_rule(model, p, v, t): """ - Per-vintage capacity constraint: var_operation[p, v, t] ≤ capacity_for_vintage(p, v, t) + Per-vintage capacity constraint: var_operation[p, v, t] ≤ units of vintage v. - This constraint ensures each vintage's operation level cannot exceed its - own production capacity. + Both variables count process UNITS: var_installation[p, v] is how many units + of vintage v exist, var_operation[p, v, t] is how many of them are running in + year t. A unit's production temporal distribution states what one unit yields + over its whole lifetime, so the annual output of a running unit is the + per-tau entry, and the bound here is a plain unit count comparison. - For greenfield (v in SYSTEM_TIME): capacity = production_rate * var_installation[p, v] - For brownfield (v not in SYSTEM_TIME): capacity = production_rate * existing_capacity[p, v] + For greenfield (v in SYSTEM_TIME): bound is var_installation[p, v] + For brownfield (v not in SYSTEM_TIME): bound is existing_capacity[p, v] """ - fg_scale = model.scales["foreground"] - op_start = pyo.value(model.process_operation_start[p]) - op_end = pyo.value(model.process_operation_end[p]) - - # Calculate production rate for this vintage - if has_production_overrides(p, r): - production_per_unit = sum( - pyo.value(get_production_value(p, r, tau_op, v)) - for tau_op in model.PROCESS_TIME - if op_start <= tau_op <= op_end - ) - else: - production_per_unit = sum( - model.foreground_production[p, r, tau_op] - for tau_op in model.PROCESS_TIME - if op_start <= tau_op <= op_end - ) - - if pyo.value(production_per_unit) == 0: - return pyo.Constraint.Skip - - # Determine capacity based on vintage type if v in model.SYSTEM_TIME: - # Greenfield: capacity from var_installation - capacity = production_per_unit * model.var_installation[p, v] * fg_scale - else: - # Brownfield: capacity from existing_capacity dict - existing_cap = model._existing_capacity_dict.get((p, v), 0) - if existing_cap == 0: - return pyo.Constraint.Skip - capacity = pyo.value(production_per_unit) * existing_cap * fg_scale - - return model.var_operation[p, v, t] <= capacity + # Greenfield: units available come from var_installation + return model.var_operation[p, v, t] <= model.var_installation[p, v] - # Build constraint over ACTIVE_VINTAGE_TIME × PRODUCT - def _build_operation_capacity_constraints(m): - """Generate constraint indices for per-vintage capacity bounds.""" - for (p, v, t) in m.ACTIVE_VINTAGE_TIME: - for r in m.PRODUCT: - yield (p, v, t, r) + # Brownfield: units available come from the existing_capacity dict + existing_cap = model._existing_capacity_dict.get((p, v), 0) + if existing_cap == 0: + return pyo.Constraint.Skip + return model.var_operation[p, v, t] <= existing_cap model.OperationCapacity = pyo.Constraint( - _build_operation_capacity_constraints(model), - rule=lambda m, p, v, t, r: operation_capacity_constraint_rule(m, p, v, t, r), + model.ACTIVE_VINTAGE_TIME, + rule=operation_capacity_constraint_rule, ) def product_demand_fulfillment_rule(model, r, t): @@ -799,7 +787,9 @@ def product_demand_fulfillment_rule(model, r, t): Demand constraint: total_production == external_demand + internal_consumption With 3D var_operation[p, v, t], sum production across all active vintages. - Each vintage may have different production rates (if overrides exist). + The rate used is the vintage's annual output at its current lifecycle stage + tau = t - v; the sum over the whole operation window is the unit's lifetime + output, not its annual output. """ total_production = 0 @@ -808,22 +798,15 @@ def product_demand_fulfillment_rule(model, r, t): if time != t: continue - op_start = pyo.value(model.process_operation_start[p]) - op_end = pyo.value(model.process_operation_end[p]) + tau = t - v # lifecycle stage of this vintage in year t + if tau not in model.PROCESS_TIME: + continue - # Get production rate for this vintage + # Get annual production rate for this vintage at this lifecycle stage if has_production_overrides(p, r): - production_rate = sum( - pyo.value(get_production_value(p, r, tau_op, v)) - for tau_op in model.PROCESS_TIME - if op_start <= tau_op <= op_end - ) + production_rate = get_production_value(p, r, tau, v) else: - production_rate = sum( - model.foreground_production[p, r, tau_op] - for tau_op in model.PROCESS_TIME - if op_start <= tau_op <= op_end - ) + production_rate = model.foreground_production[p, r, tau] total_production += production_rate * model.var_operation[p, v, t] @@ -917,29 +900,22 @@ def cumulative_category_impact_limit_rule(model, c): def total_product_flow_rule(model, r, t): """ Calculate total product output at time t, summing across all active vintages. - Consistent with 3D var_operation[p, v, t]. + Uses each vintage's annual output at its lifecycle stage tau = t - v. """ total = 0 for (p, v, time) in model.ACTIVE_VINTAGE_TIME: if time != t: continue - op_start = pyo.value(model.process_operation_start[p]) - op_end = pyo.value(model.process_operation_end[p]) + tau = t - v # lifecycle stage of this vintage in year t + if tau not in model.PROCESS_TIME: + continue - # Get production rate for this vintage + # Get annual production rate for this vintage at this lifecycle stage if has_production_overrides(p, r): - production_rate = sum( - pyo.value(get_production_value(p, r, tau_op, v)) - for tau_op in model.PROCESS_TIME - if op_start <= tau_op <= op_end - ) + production_rate = get_production_value(p, r, tau, v) else: - production_rate = sum( - model.foreground_production[p, r, tau_op] - for tau_op in model.PROCESS_TIME - if op_start <= tau_op <= op_end - ) + production_rate = model.foreground_production[p, r, tau] total += production_rate * model.var_operation[p, v, t] return total @@ -1275,11 +1251,12 @@ def validate_operation_bounds(model: pyo.ConcreteModel, tolerance: float = 1e-6) Validate that operation levels respect capacity constraints. This function performs post-solve validation to ensure that var_operation[p, v, t] - does not exceed the production capacity for each vintage. + does not exceed the units available from that vintage. Both variables count + process units, so the comparison needs no production rate. With 3D operation variables: - - Greenfield: var_operation[p, v, t] <= production_rate * var_installation[p, v] * fg_scale - - Brownfield: var_operation[p, v, t] <= production_rate * existing_capacity[p, v] * fg_scale + - Greenfield: var_operation[p, v, t] <= var_installation[p, v] + - Brownfield: var_operation[p, v, t] <= existing_capacity[p, v] Parameters ---------- @@ -1307,60 +1284,19 @@ def validate_operation_bounds(model: pyo.ConcreteModel, tolerance: float = 1e-6) violations = [] max_violation = 0.0 - fg_scale = model.scales["foreground"] - # Get sparse overrides for production with precomputed index for O(1) lookup - production_overrides = getattr(model, "_production_vintage_overrides", {}) or {} - production_overrides_index = getattr(model, "_production_overrides_index", frozenset()) existing_cap_dict = getattr(model, "_existing_capacity_dict", {}) - def get_prod_value(p, r, tau, vintage): - """Get production value, checking sparse overrides first.""" - key = (p, r, tau, vintage) - if key in production_overrides: - return production_overrides[key] - return pyo.value(model.foreground_production[p, r, tau]) - - def has_prod_overrides(p, r): - """O(1) check if any vintage overrides exist for this process/product.""" - return (p, r) in production_overrides_index - # Validate per-vintage operation bounds for (p, v, t) in model.ACTIVE_VINTAGE_TIME: operation_value = pyo.value(model.var_operation[p, v, t]) - op_start = pyo.value(model.process_operation_start[p]) - op_end = pyo.value(model.process_operation_end[p]) - - max_capacity = 0.0 - for r in model.PRODUCT: - # Get production rate for this vintage - if has_prod_overrides(p, r): - production_rate = sum( - get_prod_value(p, r, tau_op, v) - for tau_op in model.PROCESS_TIME - if op_start <= tau_op <= op_end - ) - else: - production_rate = sum( - pyo.value(model.foreground_production[p, r, tau_op]) - for tau_op in model.PROCESS_TIME - if op_start <= tau_op <= op_end - ) - - if production_rate == 0: - continue - # Calculate capacity for this vintage - if v in model.SYSTEM_TIME: - # Greenfield: capacity from var_installation - installed = pyo.value(model.var_installation[p, v]) - capacity = production_rate * installed * fg_scale - else: - # Brownfield: capacity from existing_capacity dict - existing_cap = existing_cap_dict.get((p, v), 0) - capacity = production_rate * existing_cap * fg_scale - - max_capacity = max(max_capacity, capacity) + # Units available from this vintage (operation and installation are both + # unit counts, so no production rate enters here) + if v in model.SYSTEM_TIME: + max_capacity = pyo.value(model.var_installation[p, v]) + else: + max_capacity = existing_cap_dict.get((p, v), 0) # Check if operation exceeds capacity if operation_value < -tolerance: diff --git a/src/optimex/postprocessing.py b/src/optimex/postprocessing.py index 71291fb..5dc6ac2 100644 --- a/src/optimex/postprocessing.py +++ b/src/optimex/postprocessing.py @@ -617,8 +617,25 @@ def plot_characterized_dynamic_inventory( def get_installation(self) -> pd.DataFrame: """ Extracts the installation data from the model and returns it as a DataFrame. - The DataFrame will have a MultiIndex with 'Time' and 'Process'. - The values are the installed capacities for each process at each time step. + + Values are the decision variable `var_installation[p, v]` itself, unchanged: + the number of process UNITS built in vintage year v. One unit delivers its full + lifetime production (the sum of its production temporal distribution) spread + over its operation window. + + Two things follow, and both matter when reading plots: + + 1. The index is the year of INSTALLATION, not a year of production. A unit + installed in 2030 with an operation window of tau 1-20 produces in 2031-2050. + 2. The values are a lifetime quantity, so they are not an annual capacity and + must not be compared with the per-year values from `get_production()` or + `get_demand()`. Use `get_production_capacity()` for the annual capacity that + installed units and existing stock make available in each year. + + Returns + ------- + pd.DataFrame + Time (vintage year) as index, Process as columns, units as values. """ # var_installation is already in real units, no scaling needed installation_matrix = { @@ -639,8 +656,15 @@ def get_operation(self, aggregate_vintages: bool = True) -> pd.DataFrame: """ Extracts the operation data from the model and returns it as a DataFrame. - With 3D var_operation[p, v, t], this method can either aggregate across - vintages (backward compatible) or return per-vintage data. + Values are the decision variable `var_operation[p, v, t]` itself: the number of + UNITS of vintage v that run in year t, summed over vintages by default. Unlike + `get_installation()`, the index is the year of OPERATION. + + Units running are not a production volume: multiply by the output per unit and + year (the production entry at the vintage's lifecycle stage) to get production, + or simply use `get_production()`. Units running can be compared directly with + `get_installation()` only per vintage, since operation of a vintage is bounded + by the units installed in that vintage. Parameters ---------- @@ -655,6 +679,7 @@ def get_operation(self, aggregate_vintages: bool = True) -> pd.DataFrame: If aggregate_vintages=True: DataFrame with Time as index, Process as columns. If aggregate_vintages=False: DataFrame with Time as index, (Process, Vintage) MultiIndex columns. + Values are counts of running units in both cases. Note: var_operation is not scaled because when both demand and foreground_production are scaled by the same factor, the scaling @@ -719,9 +744,6 @@ def has_production_overrides(p, r): for p in self.m.PROCESS: for f in self.m.PRODUCT: - op_start = pyo.value(self.m.process_operation_start[p]) - op_end = pyo.value(self.m.process_operation_end[p]) - for t in self.m.SYSTEM_TIME: # Sum production across all active vintages at time t total_production = 0 @@ -729,18 +751,16 @@ def has_production_overrides(p, r): if proc != p or time != t: continue - # Get production rate for this vintage + tau = t - v # lifecycle stage of this vintage in year t + if tau not in self.m.PROCESS_TIME: + continue + + # Annual output per running unit at this lifecycle stage if has_production_overrides(p, f): - production_rate = sum( - pyo.value(get_production_value(p, f, tau_op, v)) - for tau_op in self.m.PROCESS_TIME - if op_start <= tau_op <= op_end - ) + production_rate = pyo.value(get_production_value(p, f, tau, v)) else: - production_rate = sum( - pyo.value(self.m.foreground_production[p, f, tau_op]) - for tau_op in self.m.PROCESS_TIME - if op_start <= tau_op <= op_end + production_rate = pyo.value( + self.m.foreground_production[p, f, tau] ) # Production from this vintage @@ -863,6 +883,12 @@ def plot_installation(self, df_installation=None, annotated=True): """ Plot a stacked bar chart for installation data. + Bars show `var_installation`: the number of process units built in each vintage + year, a LIFETIME quantity plotted at the year of installation. This is not an + annual capacity and does not line up with production in the same year - use + `plot_capacity_balance()` to compare production against the annual capacity + those units make available. + Parameters ---------- df_installation : DataFrame, optional @@ -885,9 +911,11 @@ def plot_installation(self, df_installation=None, annotated=True): fig, axes = self._create_clean_axes() ax = axes[0] self._apply_bar_styles( - df_installation, ax, colors, title="Installed Capacity" + df_installation, ax, colors, title="Installed Units by Vintage" + ) + ax.set_ylabel( + "Installed units (lifetime)", fontsize=self._plot_config["label_fontsize"] ) - ax.set_ylabel("Installed Capacity", fontsize=self._plot_config["label_fontsize"]) # Legend at bottom @@ -911,6 +939,10 @@ def plot_operation(self, df_operation=None, annotated=True): """ Plot a stacked bar chart for operation data. + Bars show `var_operation` summed over vintages: how many units run in each year. + Note the different meaning of the x-axis compared to `plot_installation()`, + which is indexed by the year of installation. + Parameters ---------- df_operation : DataFrame, optional @@ -933,9 +965,9 @@ def plot_operation(self, df_operation=None, annotated=True): fig, axes = self._create_clean_axes() ax = axes[0] self._apply_bar_styles( - df_operation, ax, colors, title="Operational Level" + df_operation, ax, colors, title="Units Running" ) - ax.set_ylabel("Operation Level", fontsize=self._plot_config["label_fontsize"]) + ax.set_ylabel("Units running", fontsize=self._plot_config["label_fontsize"]) # Legend at bottom @@ -1008,20 +1040,23 @@ def get_existing_capacity(self) -> pd.DataFrame: def get_production_capacity(self) -> pd.DataFrame: """ - Calculate maximum available production capacity for each product at each time step. - - Capacity is determined by counting installations in their operation phase and - multiplying by their production coefficients. This includes both new installations - (from var_installation) and existing (brownfield) capacity. - - Note: Uses vintage-aware 4D calculation when production overrides exist, - matching the optimizer's capacity constraint calculation. + Calculate maximum available ANNUAL production capacity for each product at + each time step. + + Installed units (`get_installation()`) are counted in process units, and one + unit yields its full lifetime production over the whole operation window. This + method converts those units into the output they can deliver *in a given year*: + for every vintage active at time t, its unit count is multiplied by the + production coefficient at its current lifecycle stage tau = t - v. This is the + quantity to compare against actual production (`get_production()`), which is + also per year. Includes both new installations (from var_installation) and + existing (brownfield) capacity. Returns ------- pd.DataFrame DataFrame with Time as index and Products as columns. - Values represent maximum production capacity (not actual production). + Values represent maximum annual production capacity (not actual production). """ capacity_tensor = {} fg_scale = getattr(self.m, "scales", {}).get("foreground", 1.0) @@ -1042,73 +1077,33 @@ def has_production_overrides(p, r): """Check if any vintage overrides exist for this process/product.""" return (p, r) in production_overrides_index + def units_available(p, v): + """Units of vintage v that can run: greenfield installs or brownfield stock.""" + if v in self.m.SYSTEM_TIME: + return pyo.value(self.m.var_installation[p, v]) + return existing_cap_dict.get((p, v), 0) + for f in self.m.PRODUCT: for t in self.m.SYSTEM_TIME: - # Calculate total capacity across all processes + # Sum annual capacity over all vintages active at time t total_capacity = 0 - for p in self.m.PROCESS: - op_start = pyo.value(self.m.process_operation_start[p]) - op_end = pyo.value(self.m.process_operation_end[p]) + for (p, v, time) in self.m.ACTIVE_VINTAGE_TIME: + if time != t: + continue + + tau = t - v # lifecycle stage of this vintage in year t + if tau not in self.m.PROCESS_TIME: + continue if has_production_overrides(p, f): - # 4D vintage-aware capacity calculation - # Each vintage may have different production rates - process_capacity = 0 - - # New installations: sum capacity by vintage - for tau in self.m.PROCESS_TIME: - vintage = t - tau - if vintage in self.m.SYSTEM_TIME and op_start <= tau <= op_end: - # Production rate for this vintage (sum over all operating taus) - production_per_unit = sum( - get_production_value(p, f, tau_op, vintage) - for tau_op in self.m.PROCESS_TIME - if op_start <= tau_op <= op_end - ) - installation = pyo.value(self.m.var_installation[p, vintage]) - process_capacity += production_per_unit * installation - - # Existing (brownfield) capacity - for (proc, inst_year), capacity in existing_cap_dict.items(): - if proc == p: - tau_existing = t - inst_year - if op_start <= tau_existing <= op_end: - nearest_vintage = min(self.m.SYSTEM_TIME) - production_per_unit = sum( - get_production_value(p, f, tau_op, nearest_vintage) - for tau_op in self.m.PROCESS_TIME - if op_start <= tau_op <= op_end - ) - process_capacity += production_per_unit * capacity - - total_capacity += process_capacity + annual_production_per_unit = get_production_value(p, f, tau, v) else: - # 3D calculation: no overrides, all vintages have same production rate - # Count new installations in operation phase at time t - installations_operating = sum( - pyo.value(self.m.var_installation[p, t - tau]) - for tau in self.m.PROCESS_TIME - if (t - tau in self.m.SYSTEM_TIME) - and op_start <= tau <= op_end + annual_production_per_unit = pyo.value( + self.m.foreground_production[p, f, tau] ) - # Add existing (brownfield) capacity in operation phase - for (proc, inst_year), capacity in existing_cap_dict.items(): - if proc == p: - tau_existing = t - inst_year - if op_start <= tau_existing <= op_end: - installations_operating += capacity - - # Production capacity per installation (sum over operation phase) - production_per_installation = sum( - pyo.value(self.m.foreground_production[p, f, tau]) - for tau in self.m.PROCESS_TIME - if op_start <= tau <= op_end - ) - - # Total capacity for this process - total_capacity += installations_operating * production_per_installation + total_capacity += annual_production_per_unit * units_available(p, v) # Store denormalized capacity capacity_tensor[(f, t)] = total_capacity * fg_scale @@ -1172,6 +1167,9 @@ def _plot_capacity_balance_on_ax( """ Plot production vs capacity lines on a given axis. + Both series are per year and in product units; the capacity line is the annual + capacity from `get_production_capacity()`, not the installed unit counts. + Parameters ---------- ax : matplotlib.axes.Axes @@ -1181,7 +1179,7 @@ def _plot_capacity_balance_on_ax( prod_df : pd.DataFrame Production DataFrame from get_production(). capacity_df : pd.DataFrame - Capacity DataFrame from get_production_capacity(). + Annual capacity DataFrame from get_production_capacity(). annotated : bool, default=True If True, show human-readable names instead of codes. show_legend : bool, default=True @@ -1214,7 +1212,7 @@ def _plot_capacity_balance_on_ax( max_capacity.values, marker='s', linewidth=self._plot_config["line_width"], - label='Max Capacity', + label='Max annual capacity', color='#000000', linestyle='--', zorder=3 @@ -1234,7 +1232,7 @@ def _plot_capacity_balance_on_ax( # Set labels and title self._set_smart_xticks(ax, actual_production.index) - ax.set_ylabel("Quantity", fontsize=self._plot_config["label_fontsize"]) + ax.set_ylabel("Quantity per year", fontsize=self._plot_config["label_fontsize"]) ax.set_axisbelow(True) ax.grid( axis="both", @@ -1268,7 +1266,10 @@ def _compute_capacity_breakdown(self, product): existing_additions_df, existing_removals_df, operation_df. All DataFrames have process columns and time index. - Note: Uses vintage-aware 4D calculation when production overrides exist. + All capacity figures are ANNUAL production capacity: unit counts multiplied by + the production coefficient at the relevant lifecycle stage, so they are + directly comparable with per-year production. Uses vintage-specific rates when + production overrides exist. """ fg_scale = getattr(self.m, "scales", {}).get("foreground", 1.0) existing_cap_dict = getattr(self.m, "_existing_capacity_dict", {}) @@ -1288,6 +1289,14 @@ def has_production_overrides(p, r): """Check if any vintage overrides exist for this process/product.""" return (p, r) in production_overrides_index + def annual_rate(p, r, tau, vintage): + """Annual output per running unit of `vintage` at lifecycle stage `tau`.""" + if tau not in self.m.PROCESS_TIME: + return 0.0 + if has_production_overrides(p, r): + return pyo.value(get_production_value(p, r, tau, vintage)) + return pyo.value(self.m.foreground_production[p, r, tau]) + capacity_additions = {p: {} for p in self.m.PROCESS} capacity_removals = {p: {} for p in self.m.PROCESS} existing_additions = {p: {} for p in self.m.PROCESS} @@ -1299,14 +1308,12 @@ def has_production_overrides(p, r): op_start = pyo.value(self.m.process_operation_start[p]) op_end = pyo.value(self.m.process_operation_end[p]) - # Base production per installation (3D, for processes without overrides) - prod_per_inst_3d = sum( - pyo.value(self.m.foreground_production[p, product, tau]) + # Skip processes that never produce this product + if all( + annual_rate(p, product, tau, min(self.m.SYSTEM_TIME)) == 0 for tau in self.m.PROCESS_TIME if op_start <= tau <= op_end - ) - - if prod_per_inst_3d == 0: + ): capacity_additions[p][t] = 0 capacity_removals[p][t] = 0 existing_additions[p][t] = 0 @@ -1314,115 +1321,65 @@ def has_production_overrides(p, r): operation[p][t] = 0 continue - if has_production_overrides(p, product): - # 4D vintage-aware calculations - - # New capacity entering operation (vintage = t - op_start) - t_entering = t - op_start - if t_entering in self.m.SYSTEM_TIME: - installation_entering = pyo.value(self.m.var_installation[p, t_entering]) - prod_per_inst_vintage = sum( - get_production_value(p, product, tau_op, t_entering) - for tau_op in self.m.PROCESS_TIME - if op_start <= tau_op <= op_end - ) - capacity_additions[p][t] = installation_entering * prod_per_inst_vintage * fg_scale - else: - capacity_additions[p][t] = 0 - - # Capacity exiting operation (vintage = t - op_end - 1) - t_exiting = t - op_end - 1 - if t_exiting in self.m.SYSTEM_TIME: - installation_exiting = pyo.value(self.m.var_installation[p, t_exiting]) - prod_per_inst_vintage = sum( - get_production_value(p, product, tau_op, t_exiting) - for tau_op in self.m.PROCESS_TIME - if op_start <= tau_op <= op_end - ) - capacity_removals[p][t] = installation_exiting * prod_per_inst_vintage * fg_scale - else: - capacity_removals[p][t] = 0 - - # Existing capacity changes (use nearest vintage for rate) - existing_add = 0 - existing_rem = 0 - nearest_vintage = min(self.m.SYSTEM_TIME) - prod_per_inst_existing = sum( - get_production_value(p, product, tau_op, nearest_vintage) - for tau_op in self.m.PROCESS_TIME - if op_start <= tau_op <= op_end + # New capacity entering operation (vintage = t - op_start), valued at + # the annual output it delivers in its first operating year + t_entering = t - op_start + if t_entering in self.m.SYSTEM_TIME: + installation_entering = pyo.value(self.m.var_installation[p, t_entering]) + capacity_additions[p][t] = ( + installation_entering + * annual_rate(p, product, op_start, t_entering) + * fg_scale ) - for (proc, inst_year), capacity in existing_cap_dict.items(): - if proc == p: - tau_existing = t - inst_year - tau_existing_prev = (t - 1) - inst_year - if op_start <= tau_existing <= op_end: - if tau_existing_prev < op_start: - existing_add += capacity * prod_per_inst_existing * fg_scale - if tau_existing > op_end: - if op_start <= tau_existing_prev <= op_end: - existing_rem += capacity * prod_per_inst_existing * fg_scale - existing_additions[p][t] = existing_add - existing_removals[p][t] = existing_rem - - # Operation level - sum production across all active vintages - total_operation = 0 - for (proc, v, time) in self.m.ACTIVE_VINTAGE_TIME: - if proc != p or time != t: - continue - # Get production rate for this vintage - production_rate = sum( - pyo.value(get_production_value(p, product, tau_op, v)) - for tau_op in self.m.PROCESS_TIME - if op_start <= tau_op <= op_end - ) - total_operation += production_rate * pyo.value(self.m.var_operation[p, v, t]) + else: + capacity_additions[p][t] = 0 - operation[p][t] = total_operation * fg_scale + # Capacity exiting operation (vintage = t - op_end - 1), valued at the + # annual output it delivered in its last operating year + t_exiting = t - op_end - 1 + if t_exiting in self.m.SYSTEM_TIME: + installation_exiting = pyo.value(self.m.var_installation[p, t_exiting]) + capacity_removals[p][t] = ( + installation_exiting + * annual_rate(p, product, op_end, t_exiting) + * fg_scale + ) else: - # 3D calculation: no overrides - prod_per_inst = prod_per_inst_3d - - # New capacity entering operation - t_entering = t - op_start - if t_entering in self.m.SYSTEM_TIME: - installation_entering = pyo.value(self.m.var_installation[p, t_entering]) - capacity_additions[p][t] = installation_entering * prod_per_inst * fg_scale - else: - capacity_additions[p][t] = 0 + capacity_removals[p][t] = 0 - # Capacity exiting operation - t_exiting = t - op_end - 1 - if t_exiting in self.m.SYSTEM_TIME: - installation_exiting = pyo.value(self.m.var_installation[p, t_exiting]) - capacity_removals[p][t] = installation_exiting * prod_per_inst * fg_scale - else: - capacity_removals[p][t] = 0 - - # Existing capacity changes - existing_add = 0 - existing_rem = 0 - for (proc, inst_year), capacity in existing_cap_dict.items(): - if proc == p: - tau_existing = t - inst_year - tau_existing_prev = (t - 1) - inst_year - if op_start <= tau_existing <= op_end: - if tau_existing_prev < op_start: - existing_add += capacity * prod_per_inst * fg_scale - if tau_existing > op_end: - if op_start <= tau_existing_prev <= op_end: - existing_rem += capacity * prod_per_inst * fg_scale - existing_additions[p][t] = existing_add - existing_removals[p][t] = existing_rem - - # Operation level - sum production across all active vintages - total_operation = 0 - for (proc, v, time) in self.m.ACTIVE_VINTAGE_TIME: - if proc != p or time != t: - continue - total_operation += prod_per_inst * pyo.value(self.m.var_operation[p, v, t]) + # Existing (brownfield) capacity entering/leaving operation + existing_add = 0 + existing_rem = 0 + for (proc, inst_year), capacity in existing_cap_dict.items(): + if proc != p: + continue + tau_existing = t - inst_year + tau_existing_prev = (t - 1) - inst_year + if op_start <= tau_existing <= op_end and tau_existing_prev < op_start: + existing_add += ( + capacity + * annual_rate(p, product, op_start, inst_year) + * fg_scale + ) + if tau_existing > op_end and op_start <= tau_existing_prev <= op_end: + existing_rem += ( + capacity + * annual_rate(p, product, op_end, inst_year) + * fg_scale + ) + existing_additions[p][t] = existing_add + existing_removals[p][t] = existing_rem + + # Actual production at time t, summed over all active vintages + total_operation = 0 + for (proc, v, time) in self.m.ACTIVE_VINTAGE_TIME: + if proc != p or time != t: + continue + total_operation += annual_rate(p, product, t - v, v) * pyo.value( + self.m.var_operation[p, v, t] + ) - operation[p][t] = total_operation * fg_scale + operation[p][t] = total_operation * fg_scale # Convert to DataFrames capacity_additions_df = pd.DataFrame(capacity_additions) @@ -1474,6 +1431,12 @@ def _plot_capacity_balance_detailed_on_ax( """ Plot detailed capacity balance with grouped bars on a given axis. + Everything is in product units per year. The capacity bars show the annual + capacity entering and leaving operation - unit counts converted with the output + per unit and year, and placed in the year operation starts or ends. They are + therefore neither `var_installation` nor the year it was installed, unlike + `plot_installation()`. + Parameters ---------- ax : matplotlib.axes.Axes @@ -1565,8 +1528,10 @@ def _plot_capacity_balance_detailed_on_ax( process_legend = [Patch(facecolor=self._color_map.get(process_codes[i], 'black'), edgecolor='black', linewidth=0.5, label=col) for i, col in enumerate(capacity_additions_df.columns)] type_legend = [ - Patch(facecolor="white", edgecolor='#30A834', linewidth=2, label='+ Cap'), - Patch(facecolor="white", edgecolor='#CD221F', linewidth=2, label='− Cap'), + Patch(facecolor="white", edgecolor='#30A834', linewidth=2, + label='+ Annual cap.'), + Patch(facecolor="white", edgecolor='#CD221F', linewidth=2, + label='− Annual cap.'), ] # Plot production and capacity lines @@ -1574,7 +1539,7 @@ def _plot_capacity_balance_detailed_on_ax( linewidth=self._plot_config["line_width"], label='Production / Demand', color='#00549F', linestyle='-', zorder=3) ax.plot(x_positions, max_capacity.values, marker='s', - linewidth=self._plot_config["line_width"], label='Max Capacity', + linewidth=self._plot_config["line_width"], label='Max annual capacity', color='#000000', linestyle='--', zorder=3) # Line legend entries @@ -1582,11 +1547,11 @@ def _plot_capacity_balance_detailed_on_ax( Line2D([0], [0], color='#00549F', marker='o', linestyle='-', linewidth=self._plot_config["line_width"], label='Production / Demand'), Line2D([0], [0], color='#000000', marker='s', linestyle='--', - linewidth=self._plot_config["line_width"], label='Max Capacity'), + linewidth=self._plot_config["line_width"], label='Max annual capacity'), ] self._set_smart_xticks(ax, actual_production.index) - ax.set_ylabel("Quantity", fontsize=self._plot_config["label_fontsize"]) + ax.set_ylabel("Quantity per year", fontsize=self._plot_config["label_fontsize"]) ax.set_axisbelow(True) ax.grid( axis="both", @@ -1615,17 +1580,26 @@ def plot_capacity_balance(self, product=None, prod_df=None, capacity_df=None, de """ Plot actual production vs maximum available capacity. + Everything in this plot is in PRODUCT UNITS PER YEAR, which is what makes + production and capacity comparable. The capacity shown is therefore not + `var_installation`: it is `get_production_capacity()`, i.e. the units of every + vintage that is in its operation phase multiplied by the output that vintage + yields per unit and year, and it is indexed by the year the capacity is + available rather than the year it was installed. + When a specific product is given, plots a single chart. When product is None, auto-detects all products with non-zero demand or production and plots a grid of subplots. Shows two lines per product: - Production (demand is assumed equal and overlaid) - - Maximum available capacity (dashed line) + - Maximum available annual capacity (dashed line) When detailed=True, also shows grouped bars per time step: - - Left bar: Capacity changes (additions/removals stacked by process) - - Right bar: Operation level (stacked by process) + - Left bar: Annual capacity entering/leaving operation, stacked by process. + A cohort installed in year v appears here in the year it starts operating + (v + operation start), again converted to output per year. + - Right bar: Production of the running units, stacked by process Parameters ---------- @@ -1786,6 +1760,12 @@ def plot_utilization_heatmap(self, product=None, annotated=True, show_values=Tru This provides a clean, dedicated view of which processes are being operated vs sitting idle at each time step. + Utilization is computed per year as actual production divided by the annual + capacity of the vintages in their operation phase - both in product units per + year. It is not `var_operation / var_installation`: those are unit counts + indexed by different years (operation year vs vintage year), and a unit only + counts towards capacity while it is inside its operation window. + Parameters ---------- product : str, optional @@ -1795,7 +1775,7 @@ def plot_utilization_heatmap(self, product=None, annotated=True, show_values=Tru show_values : bool, default=True If True, show utilization percentages in cells. - Note: Uses vintage-aware 4D calculation when production overrides exist. + Note: Uses vintage-specific production rates when overrides exist. """ # Get demand to determine product demand_df = self.get_demand() @@ -1824,6 +1804,20 @@ def has_production_overrides(p, r): """Check if any vintage overrides exist for this process/product.""" return (p, r) in production_overrides_index + def annual_rate(p, r, tau, vintage): + """Annual output per running unit of `vintage` at lifecycle stage `tau`.""" + if tau not in self.m.PROCESS_TIME: + return 0.0 + if has_production_overrides(p, r): + return pyo.value(get_production_value(p, r, tau, vintage)) + return pyo.value(self.m.foreground_production[p, r, tau]) + + def units_available(p, v): + """Units of vintage v that can run: greenfield installs or brownfield stock.""" + if v in self.m.SYSTEM_TIME: + return pyo.value(self.m.var_installation[p, v]) + return existing_cap_dict.get((p, v), 0) + # Calculate utilization for each process at each time utilization_data = {} capacity_data = {} @@ -1833,92 +1827,32 @@ def has_production_overrides(p, r): op_start = pyo.value(self.m.process_operation_start[p]) op_end = pyo.value(self.m.process_operation_end[p]) - # Check if process produces this product (3D base rate) - prod_per_inst_3d = sum( - pyo.value(self.m.foreground_production[p, product, tau]) + # Skip processes that don't produce this product + if all( + annual_rate(p, product, tau, min(self.m.SYSTEM_TIME)) == 0 for tau in self.m.PROCESS_TIME if op_start <= tau <= op_end - ) - - if prod_per_inst_3d == 0: - continue # Skip processes that don't produce this product + ): + continue utilization_data[p] = {} capacity_data[p] = {} operation_data[p] = {} for t in self.m.SYSTEM_TIME: - if has_production_overrides(p, product): - # 4D vintage-aware calculations - # Capacity from new installations (sum by vintage) - capacity = 0 - for tau in self.m.PROCESS_TIME: - vintage = t - tau - if vintage in self.m.SYSTEM_TIME and op_start <= tau <= op_end: - production_per_unit = sum( - get_production_value(p, product, tau_op, vintage) - for tau_op in self.m.PROCESS_TIME - if op_start <= tau_op <= op_end - ) - installation = pyo.value(self.m.var_installation[p, vintage]) - capacity += production_per_unit * installation - - # Add existing (brownfield) capacity - nearest_vintage = min(self.m.SYSTEM_TIME) - prod_per_inst_existing = sum( - get_production_value(p, product, tau_op, nearest_vintage) - for tau_op in self.m.PROCESS_TIME - if op_start <= tau_op <= op_end - ) - for (proc, inst_year), cap in existing_cap_dict.items(): - if proc == p: - tau_existing = t - inst_year - if op_start <= tau_existing <= op_end: - capacity += prod_per_inst_existing * cap - - capacity *= fg_scale - - # Operation - sum production across all active vintages - operation = 0 - for (proc, v, time) in self.m.ACTIVE_VINTAGE_TIME: - if proc != p or time != t: - continue - # Get production rate for this vintage - production_rate = sum( - pyo.value(get_production_value(p, product, tau_op, v)) - for tau_op in self.m.PROCESS_TIME - if op_start <= tau_op <= op_end - ) - operation += production_rate * pyo.value(self.m.var_operation[p, v, t]) - operation *= fg_scale - else: - # 3D calculation: no overrides - prod_per_inst = prod_per_inst_3d - - # Calculate capacity from new installations - installations_operating = sum( - pyo.value(self.m.var_installation[p, t - tau]) - for tau in self.m.PROCESS_TIME - if (t - tau in self.m.SYSTEM_TIME) - and op_start <= tau <= op_end - ) - - # Add existing (brownfield) capacity in operation phase - for (proc, inst_year), cap in existing_cap_dict.items(): - if proc == p: - tau_existing = t - inst_year - if op_start <= tau_existing <= op_end: - installations_operating += cap - - capacity = installations_operating * prod_per_inst * fg_scale - - # Calculate operation - sum across all active vintages - operation = 0 - for (proc, v, time) in self.m.ACTIVE_VINTAGE_TIME: - if proc != p or time != t: - continue - operation += prod_per_inst * pyo.value(self.m.var_operation[p, v, t]) - operation *= fg_scale + # Both capacity and operation are annual quantities: unit counts times + # the output per unit and year at each vintage's lifecycle stage + capacity = 0 + operation = 0 + for (proc, v, time) in self.m.ACTIVE_VINTAGE_TIME: + if proc != p or time != t: + continue + rate = annual_rate(p, product, t - v, v) + capacity += rate * units_available(p, v) + operation += rate * pyo.value(self.m.var_operation[p, v, t]) + + capacity *= fg_scale + operation *= fg_scale capacity_data[p][t] = capacity operation_data[p][t] = operation diff --git a/tests/test_background_process_product_separation.py b/tests/test_background_process_product_separation.py index 662cd80..c3ae4df 100644 --- a/tests/test_background_process_product_separation.py +++ b/tests/test_background_process_product_separation.py @@ -139,15 +139,19 @@ def test_separated_background_matches_standard_lca(setup_separated_background_sy Test that optimex produces same results as standard LCA when background has separate process and product nodes. - Expected calculation for 100 kg Widget: - - Direct CO2 from operation: 100 * 5 * 1.0 = 500 kg CO2 - - Background electricity at construction: 100 * 10 * 0.5 = 500 kg CO2 - - Total: 1000 kg CO2 + The process has a two-year operation window and yields 0.5 kg per unit and year, + so the demand is placed in both operating years of one cohort: every unit is then + fully used, which is the condition for optimex and standard LCA to agree. + + Expected calculation for 200 kg Widget: + - Direct CO2 from operation: 200 * 5 * 1.0 = 1000 kg CO2 + - Background electricity at construction: 200 * 10 * 0.5 = 1000 kg CO2 + - Total: 2000 kg CO2 """ # Standard LCA calculation widget = bd.get_node(database="foreground", name="Widget") - lca = bc.LCA({widget: 100}, method=("GWP", "example")) + lca = bc.LCA({widget: 200}, method=("GWP", "example")) lca.lci() lca.lcia() expected_gwp = lca.score @@ -158,7 +162,7 @@ def test_separated_background_matches_standard_lca(setup_separated_background_sy years = range(2020, 2030) td_demand = TemporalDistribution( date=np.array([datetime(year, 1, 1).isoformat() for year in years], dtype='datetime64[s]'), - amount=np.asarray([0, 0, 100, 0, 0, 0, 0, 0, 0, 0]), # 100 kg at year 2022 + amount=np.asarray([0, 0, 100, 100, 0, 0, 0, 0, 0, 0]), # 100 kg in 2022 and 2023 ) lca_config = lca_processor.LCAConfig( diff --git a/tests/test_brownfield.py b/tests/test_brownfield.py index 02310cd..d5685e0 100644 --- a/tests/test_brownfield.py +++ b/tests/test_brownfield.py @@ -557,8 +557,9 @@ def test_production_capacity_includes_existing(self): capacity_df = pp.get_production_capacity() # At 2020: existing capacity (tau=1) is in operation - # Capacity should include the existing 10 units * 2 (production per unit) = 20 - assert capacity_df.loc[2020, "product"] >= 20.0, ( + # get_production_capacity reports ANNUAL capacity, so the existing 10 units + # contribute 10 units * 1.0 production per unit and year = 10 + assert capacity_df.loc[2020, "product"] >= 10.0, ( f"Production capacity at 2020 should include existing capacity, " f"got {capacity_df.loc[2020, 'product']}" ) @@ -934,9 +935,10 @@ def test_brownfield_multiple_existing_entries_same_process(self): "characterization": {("GWP", "CO2", t): 1.0 for t in range(2025, 2036)}, # MULTIPLE existing capacity entries for the same process # This should NOT double the production rate + # Together they cover the demand of 1000/year (10 units * 100 per year) "existing_capacity": { - ("Plant", 2005): 0.5, # 0.5 units from 2005 - ("Plant", 2015): 0.5, # 0.5 units from 2015 + ("Plant", 2005): 5.0, # 5 units from 2005 + ("Plant", 2015): 5.0, # 5 units from 2015 }, } @@ -953,15 +955,13 @@ def test_brownfield_multiple_existing_entries_same_process(self): assert results.solver.termination_condition == pyo.TerminationCondition.optimal - # Check that operation is reasonable (should be ~demand/production_rate = 1000/3000 ≈ 0.33) - # If the bug exists, operation would be ~0.16 (half of expected) due to doubled rate + # Each running unit produces 100 per year, so meeting a demand of 1000 + # requires 10 units running - the two existing entries taken together. + # If the production rate were counted once per existing entry, half as many + # units would appear to be enough. operation_2025 = get_total_operation(solved_model, "Plant", 2025) - # Production per unit of operation = sum of production rates for operating taus - # With operation_time_limits (1, 30), there are 30 operating taus - # Each tau has production = 100, so total = 3000 per unit of operation - # To meet demand of 1000, we need operation = 1000/3000 ≈ 0.333 - expected_operation = 1000 / 3000 # ≈ 0.333 + expected_operation = 1000 / 100 # 10 units running assert operation_2025 == pytest.approx(expected_operation, rel=0.01), ( f"Operation at 2025 should be ~{expected_operation:.4f}, got {operation_2025:.4f}. " "If operation is half of expected, the production rate is being doubled incorrectly." diff --git a/tests/test_merit_order.py b/tests/test_merit_order.py index 5f53f5b..4426b8d 100644 --- a/tests/test_merit_order.py +++ b/tests/test_merit_order.py @@ -323,6 +323,11 @@ def test_operation_limits_apply_to_total(self): Scenario: - Multiple vintages operating at same time - Operation limit constrains total (not per-vintage) + + With a demand of 50 and a production of 1.0 per unit and year, 50 units must + run at 2022, spread over the 2020, 2021 and 2022 vintages. A limit of 50 is + therefore exactly binding, while a limit of 40 must make the model infeasible — + which it only does if the limit applies to the sum across vintages. """ model_inputs_dict = { "PROCESS": ["P1"], @@ -349,7 +354,7 @@ def test_operation_limits_apply_to_total(self): "mapping": {("db", t): 1.0 for t in [2020, 2021, 2022]}, "characterization": {("GWP", "CO2", t): 1.0 for t in [2020, 2021, 2022]}, # Set a maximum operation limit at 2022 - "process_operation_limits_max": {("P1", 2022): 30}, + "process_operation_limits_max": {("P1", 2022): 50}, } model_inputs = converter.OptimizationModelInputs(**model_inputs_dict) @@ -373,9 +378,24 @@ def test_operation_limits_apply_to_total(self): ) # Total operation should respect the limit - assert total_op_2022 <= 30 + 1e-6, ( - f"Total operation at 2022 ({total_op_2022:.4f}) should not exceed limit of 30" + assert total_op_2022 <= 50 + 1e-6, ( + f"Total operation at 2022 ({total_op_2022:.4f}) should not exceed limit of 50" + ) + assert pytest.approx(50.0, rel=0.01) == total_op_2022, ( + f"50 units must run at 2022 to meet the demand, got {total_op_2022:.4f}" + ) + + # A limit below the required total must be infeasible: this only holds if the + # limit bounds the sum over vintages, not each vintage separately + model_inputs_dict["process_operation_limits_max"] = {("P1", 2022): 40} + tight_inputs = converter.OptimizationModelInputs(**model_inputs_dict) + tight_model = optimizer.create_model( + inputs=tight_inputs, + objective_category="GWP", + name="operation_limits_test_tight", ) + with pytest.raises(RuntimeError): + optimizer.solve_model(tight_model, solver_name="glpk", tee=False) class TestPerVintageCapacity: diff --git a/tests/test_operation_constraints.py b/tests/test_operation_constraints.py index 0679678..c00878c 100644 --- a/tests/test_operation_constraints.py +++ b/tests/test_operation_constraints.py @@ -24,13 +24,13 @@ def test_operation_capacity_single_installation(): Test that operation is correctly bounded by a single installation. Setup: - - Single process with production capacity of 1.0 kg/operation (scaled) - - Install 10 units at year 2020 - - Demand at years 2021-2024 during operation phase + - Single process producing 0.5 kg per unit and operating year (tau=1 and tau=2), + i.e. 1.0 kg over a unit's lifetime + - Demand of 5 kg in each of the two operating years 2021 and 2022 Expected: - - Operation should be bounded by installed capacity - - var_operation <= 10 for years when installation is in operation phase + - 10 units installed at 2020: 10 units * 0.5 kg = 5 kg per year + - All 10 units run in both years, so operation equals the installed units """ model_inputs_dict = { "PROCESS": ["P1"], @@ -87,21 +87,21 @@ def test_operation_capacity_single_installation(): assert results.solver.status == pyo.SolverStatus.ok assert results.solver.termination_condition == pyo.TerminationCondition.optimal - # Check that 5 units were installed at 2020 to meet demand - # Total demand is 10 kg, each process produces 1.0 kg total across its lifecycle - # So 5 processes are sufficient (they produce at both tau=1 and tau=2) + # Check that 10 units were installed at 2020 to meet demand + # Each unit delivers 0.5 kg per operating year, so 5 kg/year needs 10 units. + # Their combined lifetime output is 10 kg, matching the total demand of 10 kg. installed_2020 = pyo.value(solved_model.var_installation["P1", 2020]) - assert pytest.approx(5.0, rel=0.01) == installed_2020 + assert pytest.approx(10.0, rel=0.01) == installed_2020 # Check that operation levels are bounded correctly - # var_operation represents number of operating units - # At 2021: 5 units operating * 1.0 kg/unit = 5 kg (demand) - # At 2022: 5 units operating * 1.0 kg/unit = 5 kg (demand) + # var_operation counts the units running in that year + # At 2021: 10 units running * 0.5 kg/unit = 5 kg (demand) + # At 2022: 10 units running * 0.5 kg/unit = 5 kg (demand) operation_2021 = get_total_operation(solved_model, "P1", 2021) operation_2022 = get_total_operation(solved_model, "P1", 2022) - assert pytest.approx(5.0, rel=0.01) == operation_2021 - assert pytest.approx(5.0, rel=0.01) == operation_2022 + assert pytest.approx(10.0, rel=0.01) == operation_2021 + assert pytest.approx(10.0, rel=0.01) == operation_2022 def test_operation_capacity_multiple_installations(): @@ -114,10 +114,10 @@ def test_operation_capacity_multiple_installations(): - High demand requiring multiple installations to operate Expected: - - At year 2022, capacity should include: - * Installations from 2021 (at tau=1) - * Installations from 2020 (at tau=2) - - Operation should be bounded by total capacity + - At year 2022, the units able to run are those installed in 2021 (at tau=1) + and in 2020 (at tau=2) + - Each running unit delivers 0.5 kg, so a demand of 15 needs 30 units running, + and therefore at least 30 units installed across 2020 and 2021 """ model_inputs_dict = { "PROCESS": ["P1"], @@ -169,19 +169,19 @@ def test_operation_capacity_multiple_installations(): install_2021 = pyo.value(solved_model.var_installation["P1", 2021]) # At 2022: - # - Installations from 2021 are at tau=1, contribute 0.5 each - # - Installations from 2020 are at tau=2, contribute 0.5 each - # Total capacity = (install_2020 + install_2021) * 0.5 + # - Installations from 2021 are at tau=1, delivering 0.5 kg each + # - Installations from 2020 are at tau=2, delivering 0.5 kg each + # Total annual capacity = (install_2020 + install_2021) * 0.5 # Get operation at 2022 operation_2022 = get_total_operation(solved_model, "P1", 2022) - # Operation should equal demand (15) since we're minimizing emissions - assert pytest.approx(15.0, rel=0.01) == operation_2022 + # 30 units must run to deliver the demand of 15 kg + assert pytest.approx(30.0, rel=0.01) == operation_2022 - # Total installations should be at least 15 (to provide capacity of 15) + # Total installations must cover those 30 running units total_installations = install_2020 + install_2021 - assert total_installations >= 14.9 # Allow small tolerance + assert total_installations >= 29.9 # Allow small tolerance def test_operation_capacity_with_varying_demand(): @@ -197,10 +197,9 @@ def test_operation_capacity_with_varying_demand(): - Operation should be bounded by capacity when demand > capacity With the capacity constraint: - - capacity = total_production × fg_scale × installations_in_operation - - Production = 2.0 (1.0 at tau=1 + 1.0 at tau=2) - - At t=2021: capacity from install_2020 (at tau=1) - - At t=2022: capacity from install_2020 (at tau=2) + install_2021 (at tau=1) + - A running unit delivers 1.0 (production at tau=1 and tau=2 is 1.0 each) + - Units able to run at t=2021: those installed in 2020 (at tau=1) + - Units able to run at t=2022: install_2020 (at tau=2) + install_2021 (at tau=1) """ model_inputs_dict = { "PROCESS": ["P1"], @@ -248,28 +247,26 @@ def test_operation_capacity_with_varying_demand(): assert results.solver.status == pyo.SolverStatus.ok # Check that operations match demands - # With production of 1.0 at tau=1 and 1.0 at tau=2: - # - total_production = 2.0 (sum over operation phase) - # - To produce 3.0, need var_operation = 3.0 / 2.0 = 1.5 - # - To produce 8.0, need var_operation = 8.0 / 2.0 = 4.0 + # With production of 1.0 per unit and operating year: + # - To produce 3.0, 3 units must run + # - To produce 8.0, 8 units must run operation_2021 = get_total_operation(solved_model, "P1", 2021) operation_2022 = get_total_operation(solved_model, "P1", 2022) - assert pytest.approx(1.5, rel=0.01) == operation_2021 # 1.5 * 2.0 = 3.0 - assert pytest.approx(4.0, rel=0.01) == operation_2022 # 4.0 * 2.0 = 8.0 + assert pytest.approx(3.0, rel=0.01) == operation_2021 # 3 units * 1.0 = 3.0 + assert pytest.approx(8.0, rel=0.01) == operation_2022 # 8 units * 1.0 = 8.0 # Verify capacity constraint was respected install_2020 = pyo.value(solved_model.var_installation["P1", 2020]) install_2021 = pyo.value(solved_model.var_installation["P1", 2021]) - # Capacity constraint: var_operation <= total_production × installations - # At 2021: capacity = 2.0 * install_2020, need >= 1.5 -> install_2020 >= 0.75 - # At 2022: capacity = 2.0 * (install_2020 + install_2021), need >= 4.0 - # -> install_2020 + install_2021 >= 2.0 + # Capacity constraint: var_operation <= installed units of that vintage + # At 2021: only the 2020 vintage can run -> install_2020 >= 3.0 + # At 2022: both vintages can run -> install_2020 + install_2021 >= 8.0 total_installations = install_2020 + install_2021 - assert total_installations >= 1.9 # At least 2 units total needed for capacity - assert install_2020 >= 0.7 # At least ~0.75 for 2021 demand + assert total_installations >= 7.9 # At least 8 units total needed + assert install_2020 >= 2.9 # At least 3 units to serve the 2021 demand def test_operation_capacity_constraint_violation_prevented(): @@ -285,12 +282,11 @@ def test_operation_capacity_constraint_violation_prevented(): - Operation should be bounded by total capacity from both installations With the capacity constraint: - - capacity = total_production × fg_scale × installations_in_operation - - Production = 2.0 (1.0 at tau=1 + 1.0 at tau=2) - - At t=2022: capacity from install_2020 (at tau=2) + install_2021 (at tau=1) - - To meet demand of 8: var_operation = 8 / 2.0 = 4.0 - - Need capacity >= 4.0 -> 2.0 * (install_2020 + install_2021) >= 4.0 - - So install_2020 + install_2021 >= 2.0 + - A running unit delivers 1.0 per year (production 1.0 at tau=1 and tau=2) + - At t=2022 the runnable units are install_2020 (at tau=2) + install_2021 (at tau=1) + - To meet a demand of 8, 8 units must run at 2022 + - Per-year deployment is capped at 5, so neither vintage can cover it alone and + the model must install in both years """ model_inputs_dict = { "PROCESS": ["P1"], @@ -323,8 +319,8 @@ def test_operation_capacity_constraint_violation_prevented(): }, # Limit each year's installation to force distribution "process_deployment_limits_max": { - ("P1", 2020): 1, # Can install max 1 at 2020 - ("P1", 2021): 1, # Can install max 1 at 2021 + ("P1", 2020): 5, # Can install max 5 at 2020 + ("P1", 2021): 5, # Can install max 5 at 2021 }, } @@ -348,19 +344,19 @@ def test_operation_capacity_constraint_violation_prevented(): # At 2022: # - Installations from 2021 are at tau=1 # - Installations from 2020 are at tau=2 - # - Total capacity = 2.0 × (install_2020 + install_2021) - # - To produce 8 with production=2.0 per operation: need var_operation = 4.0 - # - So need capacity of at least 4.0 -> need 2 total installations + # - Each running unit delivers 1.0, so 8 units must run to produce 8 operation_2022 = get_total_operation(solved_model, "P1", 2022) - # Operation should equal demand requirement (8 / 2.0 = 4.0) - assert pytest.approx(4.0, rel=0.01) == operation_2022 + assert pytest.approx(8.0, rel=0.01) == operation_2022 - # Total installations should be at least 2.0 - # (since capacity = 2.0 × installations, need capacity >= 4.0) + # 8 units must be installed, and with a cap of 5 per year both vintages are used total_installations = install_2020 + install_2021 - assert total_installations >= 1.9 # At least 2 units needed + assert total_installations >= 7.9 + assert install_2020 > 0 and install_2021 > 0, ( + "The per-year deployment cap of 5 forces installation in both years, " + f"got 2020={install_2020}, 2021={install_2021}" + ) def test_operation_capacity_with_non_constant_production(): diff --git a/tests/test_optimization.py b/tests/test_optimization.py index 63593fd..5d2f59a 100644 --- a/tests/test_optimization.py +++ b/tests/test_optimization.py @@ -92,12 +92,16 @@ def test_model_solution_is_optimal(solved_system_model): ) +# Note on magnitudes: one installed unit produces 0.5 R1 per operating year +# (production 0.5 at tau=1 and tau=2, i.e. 1.0 over its lifetime), so meeting a demand +# of 10 in a year requires 20 units running. Installation-dependent flows are incurred +# per unit, operation-dependent flows per running unit. @pytest.mark.parametrize( "model_type, expected_value", [ # ("fixed", 3.15417e-10), # Expected value for the fixed model - ("flex", 1.9172462799736082e-10), # Expected value for the flexible model - ("constrained", 1.9197314619763485e-10), # Constrained by limiting P1 to 10 installations (0.13% higher) + ("flex", 2.8168480231004236e-10), # Expected value for the flexible model + ("constrained", 2.827868523922197e-10), # Constrained by limiting P1 to 10 installations (0.39% higher) ], ids=["flex_result", "constrained_process_limit"], # "fixed_result", ) @@ -124,19 +128,22 @@ def test_model_scaling_values_within_tolerance(solved_system_model): model.name == "abstract_system_model_fixed" or model.name == "abstract_system_model_flex" ): + # 20 units per cohort: each unit yields 0.5 R1 per operating year, and the + # cohort has to cover the demand of 10 in its first operating year expected_values = { - ("P1", 2025): 10.0, - ("P1", 2027): 10.0, - ("P2", 2021): 10.0, - ("P2", 2023): 10.0, + ("P1", 2025): 20.0, + ("P1", 2027): 20.0, + ("P2", 2021): 20.0, + ("P2", 2023): 20.0, } elif model.name == "abstract_system_model_constrained": - # With P1 limited to 10 total, optimizer shifts significantly more to P2 + # With P1 limited to 10 units total, the optimizer shifts almost everything to P2 expected_values = { ("P1", 2027): 10.0, - ("P2", 2021): 10.0, - ("P2", 2023): 10.0, - ("P2", 2025): 10.0, + ("P2", 2021): 20.0, + ("P2", 2023): 20.0, + ("P2", 2025): 20.0, + ("P2", 2027): 10.0, } else: pytest.skip(f"Unknown model name: {model.name}") @@ -347,20 +354,20 @@ def test_capacity_constraint_with_high_production(): inst = pyo.value(solved_model.var_installation["P1", 2020]) oper = get_total_operation(solved_model, "P1", 2020) - # With production = 5.0 and demand = 10: - # - Need var_operation = 10/5 = 2 to produce 10 units - # - Capacity = 5.0 * installations, so need 5.0 * inst >= 2 - # - Minimum installations = 2/5 = 0.4 - assert pytest.approx(0.4, rel=0.01) == inst, ( - f"With production=5.0, only 0.4 installations needed, got {inst}" + # With production = 5.0 per unit and year (operation window is the single tau=0) + # and demand = 10: + # - 2 units must be running to produce 10 units of product + # - Operation is bounded by installed units, so 2 units must be installed + assert pytest.approx(2.0, rel=0.01) == inst, ( + f"With production=5.0 per unit, 2 units are needed for a demand of 10, got {inst}" ) assert pytest.approx(2.0, rel=0.01) == oper, ( - f"Operation should be 2.0 to produce 10 units, got {oper}" + f"Operation should be 2.0 units running to produce 10 units, got {oper}" ) - # Verify emissions (capacity-dependent) = 0.4 * 10 = 4 kg CO2 - assert pytest.approx(4.0, rel=0.01) == objective, ( - f"Emissions should be 4.0 kg CO2, got {objective}" + # Verify emissions (installation-dependent) = 2 * 10 = 20 kg CO2 + assert pytest.approx(20.0, rel=0.01) == objective, ( + f"Emissions should be 20.0 kg CO2, got {objective}" ) diff --git a/tests/test_single_route_lca_comparison.py b/tests/test_single_route_lca_comparison.py index 32eff90..b4e4c89 100644 --- a/tests/test_single_route_lca_comparison.py +++ b/tests/test_single_route_lca_comparison.py @@ -104,24 +104,41 @@ def setup_single_route_system(): ]) -def test_single_route_matches_standard_lca(setup_single_route_system): +def test_single_year_demand_builds_stranded_capacity(setup_single_route_system): """ - Test that optimex produces same results as standard LCA for a single-route system. + Demand in a single year, for a process whose operation window spans two years. - Expected calculation for 100 kg Widget: - - Direct CO2 from operation: 100 * 5 * 1.0 = 500 kg CO2 - - Background electricity at construction: 100 * 10 * 0.5 = 500 kg CO2 - - Total: 1000 kg CO2 + One installed unit yields 1 kg Widget over its lifetime, delivered as 0.5 kg in + each of its two operating years. Asking for 100 kg in 2022 alone therefore forces + the model to install 200 units and leave the second operating year of each unit + idle: the delivered amount is 100 kg but the built capacity is worth 200 kg. + + This is NOT expected to match `bc.LCA({widget: 100})`, because standard LCA + implicitly assumes every unit is fully used. Expected here: + - Construction electricity: 200 * 10 * 0.5 = 1000 kg CO2 + - Operation CO2 in 2022: 200 * 2.5 * 1.0 = 500 kg CO2 + - Total: 1500 kg CO2 + which is the LCA score of 100 kg (1000) plus the construction impact of the + 100 units' worth of capacity that is never used (500). + + See test_multi_period_demand_matches_standard_lca for the case where the demand + profile allows full utilization and the two methods must agree exactly. """ + import pyomo.environ as pyo - # Standard LCA calculation + # Standard LCA calculation, for reference widget = bd.get_node(database="foreground", name="Widget") lca = bc.LCA({widget: 100}, method=("GWP", "example")) lca.lci() lca.lcia() - expected_gwp = lca.score + lca_gwp = lca.score - print(f"\nStandard LCA GWP: {expected_gwp}") + # Capacity worth 200 kg is built, so the impact is the LCA score of 200 kg + # minus the operation emissions of the 100 kg that is never produced. + expected_gwp = 1500.0 + + print(f"\nStandard LCA GWP (100 kg, full utilization): {lca_gwp}") + print(f"Expected optimex GWP (single-year demand): {expected_gwp}") # optimex calculation years = range(2020, 2030) @@ -155,14 +172,25 @@ def test_single_route_matches_standard_lca(setup_single_route_system): objective_category="climate_change", ) - _, obj_real, results = optimizer.solve_model(model, solver_name="glpk") + solved_model, obj_real, results = optimizer.solve_model(model, solver_name="glpk") print(f"optimex GWP: {obj_real}") - print(f"Difference: {abs(obj_real - expected_gwp)}") - # They should match (within numerical tolerance) + # 200 units must be installed to deliver 100 kg in a single year + total_installed = sum( + pyo.value(solved_model.var_installation[p, t]) + for p in solved_model.PROCESS + for t in solved_model.SYSTEM_TIME + ) + print(f"Total installed units: {total_installed}") + assert pytest.approx(total_installed, rel=1e-3) == 200.0, ( + f"Delivering 100 kg in one year needs 200 units (0.5 kg per unit and year), " + f"got {total_installed}" + ) + assert pytest.approx(obj_real, rel=1e-3) == expected_gwp, ( - f"optimex result ({obj_real}) should match standard LCA ({expected_gwp})" + f"optimex result ({obj_real}) should be {expected_gwp}: the LCA score of " + f"100 kg ({lca_gwp}) plus the construction impact of the stranded capacity" ) # Additional check: verify postprocessing extracts correct unscaled values @@ -175,9 +203,142 @@ def test_single_route_matches_standard_lca(setup_single_route_system): total_cc_from_pp = df_impacts[climate_change_cols].sum().sum() print(f"\nPostprocessing climate_change total: {total_cc_from_pp}") - print(f"Expected (from LCA): {expected_gwp}") - # Postprocessing should also match standard LCA assert pytest.approx(total_cc_from_pp, rel=1e-3) == expected_gwp, ( - f"Postprocessing climate_change sum ({total_cc_from_pp}) should match standard LCA ({expected_gwp})" + f"Postprocessing climate_change sum ({total_cc_from_pp}) should match " + f"the objective ({expected_gwp})" + ) + + +def test_multi_period_demand_matches_standard_lca(setup_single_route_system): + """ + Same single-route system, but demand at MULTIPLE time steps. + + The point of this test is installation-impact accounting. The process has a + 2-year operation window (tau 1-2) and its production temporal distribution sums + to 1 kg over that window (0.5 kg per year), so one installed unit delivers 1 kg + of Widget over its lifetime. Serving 100 kg in 2022 AND 100 kg in 2023 therefore + requires 200 units, all installed in 2021 and fully utilized in both years — the + demand profile matches the production profile exactly, so this is a case where + optimex must reproduce the standard LCA result. + + Note on test design: demand must line up with the operation window for LCA + equivalence to hold. A single isolated demand year would force the model to build + capacity whose second operating year is never used (legitimate stranded capacity), + which is more impact than the standard LCA of the delivered amount. + + Demand: 100 kg in 2022, 100 kg in 2023. + Standard LCA for 200 kg: + - Direct CO2 from operation: 200 * 5 * 1.0 = 1000 kg CO2 + - Background electricity: 200 * 10 * 0.5 = 1000 kg CO2 + - Total: 2000 kg CO2 + """ + import pyomo.environ as pyo + + demand_by_year = {2022: 100.0, 2023: 100.0} + total_demand = sum(demand_by_year.values()) + + # Standard LCA calculation for the same total amount + widget = bd.get_node(database="foreground", name="Widget") + lca = bc.LCA({widget: total_demand}, method=("GWP", "example")) + lca.lci() + lca.lcia() + expected_gwp = lca.score + + print(f"\nStandard LCA GWP ({total_demand} kg): {expected_gwp}") + + years = list(range(2020, 2030)) + td_demand = TemporalDistribution( + date=np.array( + [datetime(year, 1, 1).isoformat() for year in years], dtype="datetime64[s]" + ), + amount=np.asarray([demand_by_year.get(year, 0.0) for year in years]), + ) + + lca_config = lca_processor.LCAConfig( + demand={widget: td_demand}, + temporal={ + "start_date": datetime(2020, 1, 1), + "temporal_resolution": "year", + "time_horizon": 100, + }, + characterization_methods=[ + { + "category_name": "climate_change", + "brightway_method": ("GWP", "example"), + }, + ], + ) + + lca_data_processor = lca_processor.LCADataProcessor(lca_config) + manager = converter.ModelInputManager() + optimization_model_inputs = manager.parse_from_lca_processor(lca_data_processor) + + model = optimizer.create_model( + optimization_model_inputs, + name="test_multi_period_demand", + objective_category="climate_change", + ) + + _, obj_real, _ = optimizer.solve_model(model, solver_name="glpk") + + # Diagnostics: how much was installed, and when + installations = { + (p, t): pyo.value(model.var_installation[p, t]) + for p in model.PROCESS + for t in model.SYSTEM_TIME + if pyo.value(model.var_installation[p, t]) > 1e-9 + } + operations = { + (p, v, t): pyo.value(model.var_operation[p, v, t]) + for (p, v, t) in model.ACTIVE_VINTAGE_TIME + if pyo.value(model.var_operation[p, v, t]) > 1e-9 + } + total_installed = sum(installations.values()) + + print(f"optimex GWP: {obj_real}") + print(f"Installations (process, year) -> units: {installations}") + print(f"Total installed units: {total_installed}") + print(f"Operation (process, vintage, year) -> level: {operations}") + + # One unit of the process delivers 1 kg Widget over its lifetime, so meeting + # 300 kg of demand must require 300 installed units — no more, no less. + assert pytest.approx(total_installed, rel=1e-3) == total_demand, ( + f"Total installed units ({total_installed}) should equal total demand " + f"({total_demand}); a lower value means installation impacts are amortized " + f"over more production than the process actually delivers" + ) + + assert pytest.approx(obj_real, rel=1e-3) == expected_gwp, ( + f"optimex result ({obj_real}) should match standard LCA ({expected_gwp})" + ) + + # Per-year impact breakdown from postprocessing + pp = postprocessing.PostProcessor(model) + df_impacts = pp.get_impacts() + climate_change_cols = [c for c in df_impacts.columns if c[0] == "climate_change"] + impacts_per_year = df_impacts[climate_change_cols].sum(axis=1) + print("\nclimate_change impact per year:") + print(impacts_per_year[impacts_per_year.abs() > 1e-9]) + + total_cc_from_pp = df_impacts[climate_change_cols].sum().sum() + print(f"Postprocessing climate_change total: {total_cc_from_pp}") + + assert pytest.approx(total_cc_from_pp, rel=1e-3) == expected_gwp, ( + f"Postprocessing climate_change sum ({total_cc_from_pp}) should match " + f"standard LCA ({expected_gwp})" + ) + + # Construction electricity is consumed at tau=0, i.e. in the installation year. + # Its impact must therefore show up in exactly the years where units were installed. + install_years = {t for (_, t) in installations} + construction_impact_per_year = { + t: 0.5 * 10 * sum(v for (_, ty), v in installations.items() if ty == t) + for t in install_years + } + print(f"Expected construction impact per year: {construction_impact_per_year}") + for t, expected_construction in construction_impact_per_year.items(): + assert impacts_per_year.loc[t] >= expected_construction - 1e-6, ( + f"Impact in {t} ({impacts_per_year.loc[t]}) is below the construction " + f"impact of the units installed that year ({expected_construction})" ) diff --git a/tests/test_two_level_supply_chain.py b/tests/test_two_level_supply_chain.py index e100828..3a1a8a3 100644 --- a/tests/test_two_level_supply_chain.py +++ b/tests/test_two_level_supply_chain.py @@ -159,22 +159,31 @@ def setup_two_level_system(): def test_two_level_supply_chain_matches_lca(setup_two_level_system): - """Test that optimex produces same results as standard LCA for two-level system.""" + """ + Test that optimex produces same results as standard LCA for two-level system. + + Both processes have a two-year operation window and produce 0.5 kg per unit and + year (1 kg over a unit's lifetime). The demand is placed in both operating years + of one cohort so that every installed unit is fully used - the condition under + which optimex and standard LCA must agree. The same holds one level down: the + Product 2 cohort consumes 0.5 kg Product 1 per unit and year, which a fully + utilized Product 1 cohort supplies. + """ # Standard LCA calculation product_2 = bd.get_node(database="foreground", name="Product 2") - lca = bc.LCA({product_2: 10}, method=("GWP", "example")) + lca = bc.LCA({product_2: 20}, method=("GWP", "example")) lca.lci() lca.lcia() expected_gwp = lca.score print(f"\nStandard LCA GWP: {expected_gwp}") - # optimex calculation + # optimex calculation: 10 kg in each of the two operating years of one cohort years = range(2020, 2030) td_demand = TemporalDistribution( date=np.array([datetime(year, 1, 1).isoformat() for year in years], dtype='datetime64[s]'), - amount=np.asarray([0, 0, 10, 0, 0, 0, 0, 0, 0, 0]), + amount=np.asarray([0, 0, 10, 10, 0, 0, 0, 0, 0, 0]), ) lca_config = lca_processor.LCAConfig( @@ -239,15 +248,21 @@ def test_two_level_supply_chain_multi_temporal_demand(setup_two_level_system): 2. Installation and operation are correctly distributed over time 3. Internal demands (Product 1 for Product 2) are handled correctly 4. PostProcessor extracts correct values for each time period + + Demand comes in pairs of consecutive years with equal amounts, matching the + two-year operation window: each cohort is then fully utilized over its lifetime, + which is what makes the comparison with standard LCA exact. """ product_2 = bd.get_node(database="foreground", name="Product 2") - # Define multi-temporal demand: 10 units in 2022, 5 units in 2024, 10 units in 2026 + # Define multi-temporal demand: two fully utilized cohorts, in 2022/2023 and + # in 2026/2027 demand_schedule = { 2022: 10, - 2024: 5, - 2026: 10, + 2023: 10, + 2026: 5, + 2027: 5, } # Calculate expected total impact from standard LCA @@ -340,14 +355,18 @@ def test_two_level_supply_chain_multi_temporal_demand(setup_two_level_system): if p2_process_id in df_operation.columns: print(df_operation[[p2_process_id]]) - # Verify operation matches demand at each time point + # Verify operation matches demand at each time point. get_operation() counts + # RUNNING UNITS, and each unit yields 0.5 kg per year, so twice as many units + # run as there are kg demanded. + annual_production_per_unit = 0.5 for year, amount in demand_schedule.items(): if year in df_operation.index: operation = df_operation.loc[year, p2_process_id] - print(f"\nYear {year}: Operation={operation:.2f}, Demand={amount}") - # Operation should equal demand for single-route system - assert pytest.approx(operation, rel=1e-2) == amount, ( - f"Operation at {year} ({operation}) should match demand ({amount})" + produced = operation * annual_production_per_unit + print(f"\nYear {year}: Units running={operation:.2f}, " + f"Production={produced:.2f}, Demand={amount}") + assert pytest.approx(produced, rel=1e-2) == amount, ( + f"Production at {year} ({produced}) should match demand ({amount})" ) print("\n" + "="*80)