Optimizer¶

def create_model(
    inputs: OptimizationModelInputs,
    name: str,
    objective_category: str,
    debug_path: str = None,
) -> pyo.ConcreteModel:
    """
    Build a Pyomo ConcreteModel for the optimization problem based on the provided
    inputs.

    This function constructs a fully defined Pyomo model using data from a `OptimizationModelInputs`
    instance. It uses flexible operation mode where processes can operate between 0 and
    their maximum installed capacity.

    Parameters
    ----------
    inputs : OptimizationModelInputs
        Structured input data containing all flows, mappings, and constraints
        required for model construction.
    name : str
        Name of the Pyomo model instance.
    objective_category : str
        The category of impact to be minimized in the optimization problem.
    debug_path : str, optional
        If provided, specifies the directory path where intermediate model data (such as
        the LP formulation) or diagnostics may be stored.

    Returns
    -------
    pyo.ConcreteModel
        A fully constructed Pyomo model ready for optimization.
    """

    model = pyo.ConcreteModel(name=name)
    model._objective_category = objective_category
    scaled_inputs, scales = inputs.get_scaled_copy()
    model.scales = scales  # Store scales for denormalization later

    logger.info("Creating sets")
    # Sets
    model.PROCESS = pyo.Set(
        doc="Set of processes (or activities), indexed by p",
        initialize=scaled_inputs.PROCESS,
    )
    model.PRODUCT = pyo.Set(
        doc="Set of foreground products, indexed by r",
        initialize=scaled_inputs.PRODUCT,
    )
    model.INTERMEDIATE_FLOW = pyo.Set(
        doc="Set of background products (intermediate flows), indexed by i",
        initialize=scaled_inputs.INTERMEDIATE_FLOW,
    )
    model.ELEMENTARY_FLOW = pyo.Set(
        doc="Set of elementary flows, indexed by e",
        initialize=scaled_inputs.ELEMENTARY_FLOW,
    )
    model.FLOW = pyo.Set(
        initialize=lambda m: m.PRODUCT
        | m.INTERMEDIATE_FLOW
        | m.ELEMENTARY_FLOW,
        doc="Set of all flows, indexed by f",
    )
    model.CATEGORY = pyo.Set(
        doc="Set of impact categories, indexed by c", initialize=scaled_inputs.CATEGORY
    )

    model.BACKGROUND_ID = pyo.Set(
        doc="Set of identifiers of the prospective background databases, indexed by b",
        initialize=scaled_inputs.BACKGROUND_ID,
    )
    model.PROCESS_TIME = pyo.Set(
        doc="Set of process time points, indexed by tau",
        initialize=scaled_inputs.PROCESS_TIME,
    )
    model.SYSTEM_TIME = pyo.Set(
        doc="Set of system time points, indexed by t",
        initialize=scaled_inputs.SYSTEM_TIME,
    )

    # Parameters
    logger.info("Creating parameters")
    model.process_names = pyo.Param(
        model.PROCESS,
        within=pyo.Any,
        doc="Names of the processes",
        default=None,
        initialize=scaled_inputs.process_names,
    )
    model.demand = pyo.Param(
        model.PRODUCT,
        model.SYSTEM_TIME,
        within=pyo.Reals,
        doc="time-explicit external demand vector d",
        default=0,
        initialize=scaled_inputs.demand,
    )
    # Always use 3D base tensors - vintage overrides are applied via sparse lookup
    model.foreground_technosphere = pyo.Param(
        model.PROCESS,
        model.INTERMEDIATE_FLOW,
        model.PROCESS_TIME,
        within=pyo.Reals,
        doc="time-explicit foreground technosphere tensor A (background flows)",
        default=0,
        initialize=scaled_inputs.foreground_technosphere,
    )
    model.foreground_biosphere = pyo.Param(
        model.PROCESS,
        model.ELEMENTARY_FLOW,
        model.PROCESS_TIME,
        within=pyo.Reals,
        doc="time-explicit foreground biosphere tensor B",
        default=0,
        initialize=scaled_inputs.foreground_biosphere,
    )
    model.foreground_production = pyo.Param(
        model.PROCESS,
        model.PRODUCT,
        model.PROCESS_TIME,
        within=pyo.Reals,
        doc="time-explicit foreground production tensor F",
        default=0,
        initialize=scaled_inputs.foreground_production,
    )

    # Store sparse vintage overrides as Python dicts (not Pyomo params)
    # These are looked up at expression construction time
    model._technosphere_vintage_overrides = getattr(
        scaled_inputs, 'foreground_technosphere_vintage_overrides', None
    ) or {}
    model._biosphere_vintage_overrides = getattr(
        scaled_inputs, 'foreground_biosphere_vintage_overrides', None
    ) or {}
    model._production_vintage_overrides = getattr(
        scaled_inputs, 'foreground_production_vintage_overrides', None
    ) or {}

    # Precompute sets of (process, flow) pairs that have overrides for O(1) lookup
    model._technosphere_overrides_index = frozenset(
        (k[0], k[1]) for k in model._technosphere_vintage_overrides
    )
    model._biosphere_overrides_index = frozenset(
        (k[0], k[1]) for k in model._biosphere_vintage_overrides
    )
    model._production_overrides_index = frozenset(
        (k[0], k[1]) for k in model._production_vintage_overrides
    )

    model.internal_demand_technosphere = pyo.Param(
        model.PROCESS,
        model.PRODUCT,
        model.PROCESS_TIME,
        within=pyo.Reals,
        doc="time-explicit internal demand tensor A^{internal}",
        default=0,
        initialize=scaled_inputs.internal_demand_technosphere,
    )
    model.background_inventory = pyo.Param(
        model.BACKGROUND_ID,
        model.INTERMEDIATE_FLOW,
        model.ELEMENTARY_FLOW,
        within=pyo.Reals,
        doc="prospective background inventory tensor G",
        default=0,
        initialize=scaled_inputs.background_inventory,
    )
    model.mapping = pyo.Param(
        model.BACKGROUND_ID,
        model.SYSTEM_TIME,
        within=pyo.Reals,
        doc="time-explicit background mapping tensor M",
        default=0,
        initialize=scaled_inputs.mapping,
    )
    model.characterization = pyo.Param(
        model.CATEGORY,
        model.ELEMENTARY_FLOW,
        model.SYSTEM_TIME,
        within=pyo.Reals,
        doc="time-explicit characterization tensor Q",
        default=0,
        initialize=scaled_inputs.characterization,
    )
    model.operation_flow = pyo.Param(
        model.PROCESS,
        model.FLOW,
        within=pyo.Binary,
        doc="operation flow matrix",
        default=0,
        initialize=scaled_inputs.operation_flow,
    )
    model.process_operation_start = pyo.Param(
        model.PROCESS,
        within=pyo.NonNegativeIntegers,
        doc="start time of process operation",
        default=0,
        initialize={k: v[0] for k, v in scaled_inputs.operation_time_limits.items()},
    )
    model.process_operation_end = pyo.Param(
        model.PROCESS,
        within=pyo.NonNegativeIntegers,
        doc="end time of process operation",
        default=0,
        initialize={k: v[1] for k, v in scaled_inputs.operation_time_limits.items()},
    )
    model.process_deployment_limits_max = pyo.Param(
        model.PROCESS,
        model.SYSTEM_TIME,
        within=pyo.Reals,
        doc="maximum time specific process deployment limit S_max",
        default=scaled_inputs.process_deployment_limits_max_default,
        initialize=(
            scaled_inputs.process_deployment_limits_max
            if scaled_inputs.process_deployment_limits_max is not None
            else {}
        ),
    )
    model.process_deployment_limits_min = pyo.Param(
        model.PROCESS,
        model.SYSTEM_TIME,
        within=pyo.Reals,
        doc="minimum time specific process deployment limit S_min",
        default=scaled_inputs.process_deployment_limits_min_default,
        initialize=(
            scaled_inputs.process_deployment_limits_min
            if scaled_inputs.process_deployment_limits_min is not None
            else {}
        ),
    )
    model.process_operation_limits_max = pyo.Param(
        model.PROCESS,
        model.SYSTEM_TIME,
        within=pyo.Reals,
        doc="maximum time specific process operation limit O_max",
        default=scaled_inputs.process_operation_limits_max_default,
        initialize=(
            scaled_inputs.process_operation_limits_max
            if scaled_inputs.process_operation_limits_max is not None
            else {}
        ),
    )
    model.process_operation_limits_min = pyo.Param(
        model.PROCESS,
        model.SYSTEM_TIME,
        within=pyo.Reals,
        doc="minimum time specific process operation limit O_min",
        default=scaled_inputs.process_operation_limits_min_default,
        initialize=(
            scaled_inputs.process_operation_limits_min
            if scaled_inputs.process_operation_limits_min is not None
            else {}
        ),
    )
    model.cumulative_process_limits_max = pyo.Param(
        model.PROCESS,
        within=pyo.Reals,
        doc="maximum cumulatative process limit S_max,cum",
        default=scaled_inputs.cumulative_process_limits_max_default,
        initialize=(
            scaled_inputs.cumulative_process_limits_max
            if scaled_inputs.cumulative_process_limits_max is not None
            else {}
        ),
    )
    model.cumulative_process_limits_min = pyo.Param(
        model.PROCESS,
        within=pyo.Reals,
        doc="minimum cumulatative process limit S_min,cum",
        default=scaled_inputs.cumulative_process_limits_min_default,
        initialize=(
            scaled_inputs.cumulative_process_limits_min
            if scaled_inputs.cumulative_process_limits_min is not None
            else {}
        ),
    )
    model.process_coupling = pyo.Param(
        model.PROCESS,
        model.PROCESS,
        within=pyo.NonNegativeReals,
        doc="coupling matrix",
        initialize=(
            scaled_inputs.process_coupling
            if scaled_inputs.process_coupling is not None
            else {}
        ),
        default=0,  # Set default coupling value to 0 if not defined
    )

    # Existing (brownfield) capacity: capacity installed before SYSTEM_TIME
    # These contribute to operation capacity but NOT to installation-phase impacts
    model.existing_capacity = pyo.Param(
        model.PROCESS,
        pyo.Any,  # Installation year (can be any year before SYSTEM_TIME)
        within=pyo.NonNegativeReals,
        doc="Existing capacity (process, installation_year) -> amount",
        initialize=(
            scaled_inputs.existing_capacity
            if scaled_inputs.existing_capacity is not None
            else {}
        ),
        default=0,
    )
    # Store the existing capacity dict for iteration
    model._existing_capacity_dict = (
        scaled_inputs.existing_capacity
        if scaled_inputs.existing_capacity is not None
        else {}
    )

    # Build ACTIVE_VINTAGE_TIME index: valid (process, vintage, time) tuples
    # where vintage is active (in operation phase) at time t
    def _build_active_vintage_index(m):
        """
        Build set of valid (process, vintage, time) tuples where vintage is active at time t.

        A vintage v is active at time t if:
        - tau = t - v (process-relative time)
        - op_start <= tau <= op_end (within operation phase)

        This includes both:
        - Greenfield vintages: v in SYSTEM_TIME (new installations)
        - Brownfield vintages: v from existing_capacity (pre-existing installations)
        """
        indices = []
        for p in m.PROCESS:
            op_start = pyo.value(m.process_operation_start[p])
            op_end = pyo.value(m.process_operation_end[p])

            # Greenfield vintages (new installations within system time horizon)
            for v in m.SYSTEM_TIME:
                for t in m.SYSTEM_TIME:
                    tau = t - v  # Process-relative time
                    if v <= t and op_start <= tau <= op_end:
                        indices.append((p, v, t))

            # Brownfield vintages (existing capacity installed before system time)
            for (proc, inst_year), _ in model._existing_capacity_dict.items():
                if proc == p:
                    for t in m.SYSTEM_TIME:
                        tau = t - inst_year
                        if op_start <= tau <= op_end:
                            indices.append((p, inst_year, t))
        return indices

    model.ACTIVE_VINTAGE_TIME = pyo.Set(
        dimen=3,
        initialize=_build_active_vintage_index,
        doc="Set of (process, vintage, time) tuples where vintage is in operation phase at time t"
    )

    # Store category impact limit data for constraint generation
    model._category_impact_limits = (
        scaled_inputs.category_impact_limits
        if scaled_inputs.category_impact_limits is not None
        else {}
    )

    model.cumulative_category_impact_limits = pyo.Param(
        model.CATEGORY,
        within=pyo.Reals,
        doc="cumulative maximum impact limit per category",
        default=float("inf"),
        initialize=(
            scaled_inputs.cumulative_category_impact_limits
            if scaled_inputs.cumulative_category_impact_limits is not None
            else {}
        ),
    )

    # Variables
    logger.info("Creating variables")
    model.var_installation = pyo.Var(
        model.PROCESS,
        model.SYSTEM_TIME,
        within=pyo.NonNegativeReals,
        doc="Installation of the process",
    )

    # Deployment limits
    model.ProcessDeploymentLimitMax = pyo.Constraint(
        model.PROCESS,
        model.SYSTEM_TIME,
        rule=lambda m, p, t: m.var_installation[p, t] <= m.process_deployment_limits_max[p, t],
    )

    model.ProcessDeploymentLimitMin = pyo.Constraint(
        model.PROCESS,
        model.SYSTEM_TIME,
        rule=lambda m, p, t: m.var_installation[p, t] >= m.process_deployment_limits_min[p, t],
    )

    model.CumulativeProcessLimitMax = pyo.Constraint(
        model.PROCESS,
        rule=lambda m, p: sum(m.var_installation[p, t] for t in m.SYSTEM_TIME)
        <= m.cumulative_process_limits_max[p],
    )
    model.CumulativeProcessLimitMin = pyo.Constraint(
        model.PROCESS,
        rule=lambda m, p: sum(m.var_installation[p, t] for t in m.SYSTEM_TIME)
        >= m.cumulative_process_limits_min[p],
    )

    # Process coupling
    def process_coupling_rule(model, p1, p2, t):
        if (
            model.process_coupling[p1, p2] > 0
        ):  # only create constraint for non-zero coupling
            return (
                model.var_installation[p1, t]
                == model.process_coupling[p1, p2] * model.var_installation[p2, t]
            )
        else:
            return pyo.Constraint.Skip

    model.ProcessCouplingConstraint = pyo.Constraint(
        model.PROCESS, model.PROCESS, model.SYSTEM_TIME, rule=process_coupling_rule
    )

    def in_operation_phase(p, tau):
        return model.process_operation_start[p] <= tau <= model.process_operation_end[p]

    # 3D operation variable: (process, vintage, time)
    # Each vintage cohort has its own operation variable for merit-order dispatch
    model.var_operation = pyo.Var(
        model.ACTIVE_VINTAGE_TIME,
        within=pyo.NonNegativeReals,
        doc="Operational activity level per vintage (process, vintage, time)",
    )

    # Operation limits apply to TOTAL operation across all vintages
    def process_operation_limit_max_rule(m, p, t):
        active_vintages = [
            (proc, v, time)
            for (proc, v, time) in m.ACTIVE_VINTAGE_TIME
            if proc == p and time == t
        ]
        if not active_vintages:
            return pyo.Constraint.Skip
        total_op = sum(m.var_operation[proc, v, time] for (proc, v, time) in active_vintages)
        return total_op <= m.process_operation_limits_max[p, t]

    model.ProcessOperationLimitMax = pyo.Constraint(
        model.PROCESS,
        model.SYSTEM_TIME,
        rule=process_operation_limit_max_rule,
    )

    def process_operation_limit_min_rule(m, p, t):
        active_vintages = [
            (proc, v, time)
            for (proc, v, time) in m.ACTIVE_VINTAGE_TIME
            if proc == p and time == t
        ]
        if not active_vintages:
            return pyo.Constraint.Skip
        total_op = sum(m.var_operation[proc, v, time] for (proc, v, time) in active_vintages)
        return total_op >= m.process_operation_limits_min[p, t]

    model.ProcessOperationLimitMin = pyo.Constraint(
        model.PROCESS,
        model.SYSTEM_TIME,
        rule=process_operation_limit_min_rule,
    )

    # Expression builders using sparse vintage override lookup
    # Base tensor is always 3D; overrides are checked first for vintage-specific values
    def scale_tensor_by_installation(tensor: pyo.Param, flow_set: str, overrides: dict, overrides_index: frozenset):
        def expr(m, p, x, t):
            result = 0
            for tau in m.PROCESS_TIME:
                vintage = t - tau
                if (vintage in m.SYSTEM_TIME) and (
                    not in_operation_phase(p, tau) or not m.operation_flow[p, x]
                ):
                    # Check sparse override first, fall back to base 3D tensor
                    key = (p, x, tau, vintage)
                    if key in overrides:
                        flow_value = overrides[key]
                    else:
                        flow_value = tensor[p, x, tau]
                    result += flow_value * m.var_installation[p, vintage]
            return result

        return pyo.Expression(
            model.PROCESS, getattr(model, flow_set), model.SYSTEM_TIME, rule=expr
        )

    def scale_tensor_by_operation(tensor: pyo.Param, flow_set: str, overrides: dict, overrides_index: frozenset):
        """
        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.
        """
        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
                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
                    )
                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
                    )

                total += flow_rate * m.var_operation[p, v, t]

            return total

        return pyo.Expression(
            model.PROCESS, getattr(model, flow_set), model.SYSTEM_TIME, rule=expr
        )

    # Create expressions with sparse vintage override lookup
    model.scaled_technosphere_dependent_on_installation = (
        scale_tensor_by_installation(
            model.foreground_technosphere,
            "INTERMEDIATE_FLOW",
            model._technosphere_vintage_overrides,
            model._technosphere_overrides_index,
        )
    )
    model.scaled_biosphere_dependent_on_installation = scale_tensor_by_installation(
        model.foreground_biosphere,
        "ELEMENTARY_FLOW",
        model._biosphere_vintage_overrides,
        model._biosphere_overrides_index,
    )
    model.scaled_technosphere_dependent_on_operation = scale_tensor_by_operation(
        model.foreground_technosphere,
        "INTERMEDIATE_FLOW",
        model._technosphere_vintage_overrides,
        model._technosphere_overrides_index,
    )
    model.scaled_biosphere_dependent_on_operation = scale_tensor_by_operation(
        model.foreground_biosphere,
        "ELEMENTARY_FLOW",
        model._biosphere_vintage_overrides,
        model._biosphere_overrides_index,
    )
    # Internal demand has no vintage overrides
    _empty_index = frozenset()
    model.scaled_internal_demand_dependent_on_installation = (
        scale_tensor_by_installation(
            model.internal_demand_technosphere, "PRODUCT", {}, _empty_index
        )
    )
    model.scaled_internal_demand_dependent_on_operation = (
        scale_tensor_by_operation(
            model.internal_demand_technosphere, "PRODUCT", {}, _empty_index
        )
    )

    def scaled_inventory_tensor(model, p, e, t):
        """
        Returns a Pyomo expression for the total inventory impact for a given
        process p, elementary flow e, and time step t.
        """

        return sum(
            (
                model.scaled_technosphere_dependent_on_installation[p, i, t]
                + model.scaled_technosphere_dependent_on_operation[p, i, t]
            )
            * sum(
                model.background_inventory[bkg, i, e] * model.mapping[bkg, t]
                for bkg in model.BACKGROUND_ID
            )
            for i in model.INTERMEDIATE_FLOW
        ) + (
            model.scaled_biosphere_dependent_on_installation[p, e, t]
            + model.scaled_biosphere_dependent_on_operation[p, e, t]
        )

    model.scaled_inventory = pyo.Expression(
        model.PROCESS,
        model.ELEMENTARY_FLOW,
        model.SYSTEM_TIME,
        rule=scaled_inventory_tensor,
    )

    # Helper to get production value with sparse override lookup
    def get_production_value(p, r, tau, vintage):
        """Get production value, checking sparse overrides first."""
        key = (p, r, tau, vintage)
        if key in model._production_vintage_overrides:
            return model._production_vintage_overrides[key]
        return model.foreground_production[p, r, tau]

    # O(1) check if production overrides exist for a given process/product
    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):
        """
        Per-vintage capacity constraint: var_operation[p, v, t] ≤ capacity_for_vintage(p, v, t)

        This constraint ensures each vintage's operation level cannot exceed its
        own production capacity.

        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]
        """
        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

    # 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)

    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),
    )

    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).
        """
        total_production = 0

        # Sum production across all active vintages at time t
        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])

            # Get production rate for this vintage
            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
                )
            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
                )

            total_production += production_rate * model.var_operation[p, v, t]

        external_demand = model.demand[r, t]
        internal_consumption = sum(
            model.scaled_internal_demand_dependent_on_installation[p, r, t]
            + model.scaled_internal_demand_dependent_on_operation[p, r, t]
            for p in model.PROCESS
        )
        return total_production == external_demand + internal_consumption

    model.ProductDemandFulfillment = pyo.Constraint(
        model.PRODUCT, model.SYSTEM_TIME, rule=product_demand_fulfillment_rule
    )

    def category_process_time_specific_impact(model, c, p, t):
        return sum(
            model.characterization[c, e, t]
            * (model.scaled_inventory[p, e, t])  # Total inventory impact
            for e in model.ELEMENTARY_FLOW
        )

    # impact of process p at time t in category c
    model.specific_impact = pyo.Expression(
        model.CATEGORY,
        model.PROCESS,
        model.SYSTEM_TIME,
        rule=category_process_time_specific_impact,
    )

    # Total impact
    def total_impact_in_category(model, c):
        return sum(
            model.specific_impact[c, p, t]
            for p in model.PROCESS
            for t in model.SYSTEM_TIME
        )

    model.total_impact = pyo.Expression(model.CATEGORY, rule=total_impact_in_category)

    # Time-specific impact (impact at a specific time across all processes)
    def time_specific_impact_rule(model, c, t):
        return sum(model.specific_impact[c, p, t] for p in model.PROCESS)

    model.time_specific_impact = pyo.Expression(
        model.CATEGORY, model.SYSTEM_TIME, rule=time_specific_impact_rule
    )

    # Time-specific category impact limits
    def category_impact_limits_rule(model, c, t):
        if (c, t) in model._category_impact_limits:
            return model.time_specific_impact[c, t] <= model._category_impact_limits[(c, t)]
        return pyo.Constraint.Skip

    model.CategoryImpactLimits = pyo.Constraint(
        model.CATEGORY, model.SYSTEM_TIME, rule=category_impact_limits_rule
    )

    # Cumulative category impact limit
    def cumulative_category_impact_limit_rule(model, c):
        return model.total_impact[c] <= model.cumulative_category_impact_limits[c]

    model.CumulativeCategoryImpactLimits = pyo.Constraint(
        model.CATEGORY, rule=cumulative_category_impact_limit_rule
    )

    # Flow limits
    # Store flow limits data for constraint generation
    model._flow_limits_max = (
        scaled_inputs.flow_limits_max
        if scaled_inputs.flow_limits_max is not None
        else {}
    )
    model._flow_limits_min = (
        scaled_inputs.flow_limits_min
        if scaled_inputs.flow_limits_min is not None
        else {}
    )
    model._cumulative_flow_limits_max = (
        scaled_inputs.cumulative_flow_limits_max
        if scaled_inputs.cumulative_flow_limits_max is not None
        else {}
    )
    model._cumulative_flow_limits_min = (
        scaled_inputs.cumulative_flow_limits_min
        if scaled_inputs.cumulative_flow_limits_min is not None
        else {}
    )

    # Expression for total product output at time t (in SCALED units)
    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].
        """
        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])

            # Get production rate for this vintage
            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
                )
            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
                )

            total += production_rate * model.var_operation[p, v, t]
        return total

    model.total_product_flow = pyo.Expression(
        model.PRODUCT, model.SYSTEM_TIME, rule=total_product_flow_rule
    )

    # Expression for total intermediate flow consumed at time t (in SCALED units)
    def total_intermediate_flow_rule(model, i, t):
        return sum(
            model.scaled_technosphere_dependent_on_installation[p, i, t]
            + model.scaled_technosphere_dependent_on_operation[p, i, t]
            for p in model.PROCESS
        )

    model.total_intermediate_flow = pyo.Expression(
        model.INTERMEDIATE_FLOW, model.SYSTEM_TIME, rule=total_intermediate_flow_rule
    )

    # Expression for total elementary flow at time t (in SCALED units)
    # This includes both foreground biosphere flows AND background inventory flows
    # (flows from intermediate flows going through background databases)
    def total_elementary_flow_rule(model, e, t):
        # Foreground biosphere contribution
        foreground_flow = sum(
            model.scaled_biosphere_dependent_on_installation[p, e, t]
            + model.scaled_biosphere_dependent_on_operation[p, e, t]
            for p in model.PROCESS
        )
        # Background inventory contribution (via intermediate flows)
        background_flow = sum(
            (
                model.scaled_technosphere_dependent_on_installation[p, i, t]
                + model.scaled_technosphere_dependent_on_operation[p, i, t]
            )
            * sum(
                model.background_inventory[bkg, i, e] * model.mapping[bkg, t]
                for bkg in model.BACKGROUND_ID
            )
            for i in model.INTERMEDIATE_FLOW
            for p in model.PROCESS
        )
        return foreground_flow + background_flow

    model.total_elementary_flow = pyo.Expression(
        model.ELEMENTARY_FLOW, model.SYSTEM_TIME, rule=total_elementary_flow_rule
    )

    # Helper function to get total flow for any flow type
    def get_total_flow(model, f, t):
        if f in model.PRODUCT:
            return model.total_product_flow[f, t]
        elif f in model.INTERMEDIATE_FLOW:
            return model.total_intermediate_flow[f, t]
        elif f in model.ELEMENTARY_FLOW:
            return model.total_elementary_flow[f, t]
        else:
            return 0

    # Time-specific flow limit constraints (max)
    def flow_limit_max_rule(model, f, t):
        if (f, t) not in model._flow_limits_max:
            return pyo.Constraint.Skip
        fg_scale = model.scales["foreground"]
        limit = model._flow_limits_max[(f, t)]
        total_flow = get_total_flow(model, f, t)
        # total_flow is SCALED, convert limit to scaled units
        return total_flow <= limit / fg_scale

    model.FlowLimitMax = pyo.Constraint(
        model.FLOW, model.SYSTEM_TIME, rule=flow_limit_max_rule
    )

    # Time-specific flow limit constraints (min)
    def flow_limit_min_rule(model, f, t):
        if (f, t) not in model._flow_limits_min:
            return pyo.Constraint.Skip
        fg_scale = model.scales["foreground"]
        limit = model._flow_limits_min[(f, t)]
        total_flow = get_total_flow(model, f, t)
        # total_flow is SCALED, convert limit to scaled units
        return total_flow >= limit / fg_scale

    model.FlowLimitMin = pyo.Constraint(
        model.FLOW, model.SYSTEM_TIME, rule=flow_limit_min_rule
    )

    # Cumulative flow limit constraints (max)
    def cumulative_flow_limit_max_rule(model, f):
        if f not in model._cumulative_flow_limits_max:
            return pyo.Constraint.Skip
        fg_scale = model.scales["foreground"]
        limit = model._cumulative_flow_limits_max[f]
        total_flow = sum(get_total_flow(model, f, t) for t in model.SYSTEM_TIME)
        # total_flow is SCALED, convert limit to scaled units
        return total_flow <= limit / fg_scale

    model.CumulativeFlowLimitMax = pyo.Constraint(
        model.FLOW, rule=cumulative_flow_limit_max_rule
    )

    # Cumulative flow limit constraints (min)
    def cumulative_flow_limit_min_rule(model, f):
        if f not in model._cumulative_flow_limits_min:
            return pyo.Constraint.Skip
        fg_scale = model.scales["foreground"]
        limit = model._cumulative_flow_limits_min[f]
        total_flow = sum(get_total_flow(model, f, t) for t in model.SYSTEM_TIME)
        # total_flow is SCALED, convert limit to scaled units
        return total_flow >= limit / fg_scale

    model.CumulativeFlowLimitMin = pyo.Constraint(
        model.FLOW, rule=cumulative_flow_limit_min_rule
    )

    def objective_function(model):
        return model.total_impact[model._objective_category]

    model.OBJ = pyo.Objective(sense=pyo.minimize, rule=objective_function)

    if debug_path is not None:
        model.write(
            debug_path,
            io_options={"symbolic_solver_labels": True},
            format=ProblemFormat.cpxlp,
        )
    return model