diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md
index 80925bf7..695b21eb 100644
--- a/docs/CHANGELOG.md
+++ b/docs/CHANGELOG.md
@@ -13,6 +13,13 @@ All notable changes to GemsPy are documented here.
   compute tighter variable bounds from the first solve. Each model declares
   what a heuristic reads/writes via `models[].heuristics` in
   `optim-config.yml`.
+- **`lower_bound(variable_name)`** and **`upper_bound(variable_name)`** operators in the
+  expression language, usable in `extra-outputs` and port-field-definitions. Both take a bare
+  variable identifier and return its *current* lower/upper bound post-solve — in particular
+  reflecting mutations made by thermal heuristics (see "Integer strategy and thermal
+  heuristics" above), which previously had no way to be surfaced in results. Validated at
+  model-build time; using them inside constraints, binding-constraints, objective
+  contributions, or variable bounds raises a `ValueError`.
 
 ## [0.1.3] - 2026-07-24
 
diff --git a/docs/agents/testing.md b/docs/agents/testing.md
index 8db9e246..a7db98b9 100644
--- a/docs/agents/testing.md
+++ b/docs/agents/testing.md
@@ -23,7 +23,7 @@ nature (they build a study with `gems_craft` and solve it with `gems_runner` in
 | Unit — gems_craft / system | `tests/unittests/gems_craft/system/` | Model, network, and port object behaviour |
 | Unit — gems_craft / system parsing | `tests/unittests/gems_craft/system_parsing/` | System YAML parsing |
 | Unit — gems_craft_hybrid | `tests/unittests/gems_craft_hybrid/` | Hybrid GEMS/Antares Simulator schema parsing |
-| Unit — gems_runner / expression | `tests/unittests/gems_runner/expression/` | Solver-output expression evaluation (`dual()`, `reduced_cost()`, `variable()`) |
+| Unit — gems_runner / expression | `tests/unittests/gems_runner/expression/` | Solver-output expression evaluation (`dual()`, `reduced_cost()`, `lower_bound()`, `upper_bound()`, `variable()`) |
 | Unit — gems_runner / simulation | `tests/unittests/gems_runner/simulation/` | Full problem build + solve on small networks |
 | End-to-end — functional | `tests/e2e/functional/` | Cross-cutting tests: library/system combinations, stochastic, investment, scenario builder |
 | End-to-end — models | `tests/e2e/models/` | Model-level tests (andromede-v1 models, operator tests, proof-of-concept models) |
diff --git a/src/gems_craft/expression/copy.py b/src/gems_craft/expression/copy.py
index 90e93ff3..4e352eb5 100644
--- a/src/gems_craft/expression/copy.py
+++ b/src/gems_craft/expression/copy.py
@@ -23,6 +23,7 @@
     ExpressionNode,
     FloorNode,
     LiteralNode,
+    LowerBoundNode,
     MaxNode,
     MinNode,
     ParameterNode,
@@ -34,6 +35,7 @@
     TimeEvalNode,
     TimeShiftNode,
     TimeSumNode,
+    UpperBoundNode,
     VariableNode,
 )
 from .visitor import ExpressionVisitorOperations, visit
@@ -111,6 +113,12 @@ def dual(self, node: DualNode) -> ExpressionNode:
     def reduced_cost(self, node: ReducedCostNode) -> ExpressionNode:
         return ReducedCostNode(node.variable_id)
 
+    def lower_bound(self, node: LowerBoundNode) -> ExpressionNode:
+        return LowerBoundNode(node.variable_id)
+
+    def upper_bound(self, node: UpperBoundNode) -> ExpressionNode:
+        return UpperBoundNode(node.variable_id)
+
 
 def copy_expression(expression: ExpressionNode) -> ExpressionNode:
     return visit(expression, CopyVisitor())
diff --git a/src/gems_craft/expression/degree.py b/src/gems_craft/expression/degree.py
index 97c007fc..4fe13ef0 100644
--- a/src/gems_craft/expression/degree.py
+++ b/src/gems_craft/expression/degree.py
@@ -20,6 +20,7 @@
     CeilNode,
     DualNode,
     FloorNode,
+    LowerBoundNode,
     MaxNode,
     MinNode,
     PortFieldAggregatorNode,
@@ -29,6 +30,7 @@
     TimeEvalNode,
     TimeShiftNode,
     TimeSumNode,
+    UpperBoundNode,
 )
 
 from .expression import (
@@ -132,6 +134,12 @@ def dual(self, node: DualNode) -> int | float:
     def reduced_cost(self, node: ReducedCostNode) -> int | float:
         return math.inf
 
+    def lower_bound(self, node: LowerBoundNode) -> int | float:
+        return math.inf
+
+    def upper_bound(self, node: UpperBoundNode) -> int | float:
+        return math.inf
+
 
 def compute_degree(expression: ExpressionNode) -> int | float:
     return visit(expression, ExpressionDegreeVisitor())
diff --git a/src/gems_craft/expression/equality.py b/src/gems_craft/expression/equality.py
index 1e5b356d..ede03651 100644
--- a/src/gems_craft/expression/equality.py
+++ b/src/gems_craft/expression/equality.py
@@ -32,6 +32,7 @@
     CeilNode,
     DualNode,
     FloorNode,
+    LowerBoundNode,
     MaxNode,
     MinNode,
     PortFieldAggregatorNode,
@@ -42,6 +43,7 @@
     TimeEvalNode,
     TimeShiftNode,
     TimeSumNode,
+    UpperBoundNode,
 )
 
 
@@ -115,6 +117,10 @@ def visit(self, left: ExpressionNode, right: ExpressionNode) -> bool:
             return self.dual(left, right)
         if isinstance(left, ReducedCostNode) and isinstance(right, ReducedCostNode):
             return self.reduced_cost(left, right)
+        if isinstance(left, LowerBoundNode) and isinstance(right, LowerBoundNode):
+            return self.lower_bound(left, right)
+        if isinstance(left, UpperBoundNode) and isinstance(right, UpperBoundNode):
+            return self.upper_bound(left, right)
         raise NotImplementedError(f"Equality not implemented for {left.__class__}")
 
     def literal(self, left: LiteralNode, right: LiteralNode) -> bool:
@@ -217,6 +223,12 @@ def dual(self, left: DualNode, right: DualNode) -> bool:
     def reduced_cost(self, left: ReducedCostNode, right: ReducedCostNode) -> bool:
         return left.variable_id == right.variable_id
 
+    def lower_bound(self, left: LowerBoundNode, right: LowerBoundNode) -> bool:
+        return left.variable_id == right.variable_id
+
+    def upper_bound(self, left: UpperBoundNode, right: UpperBoundNode) -> bool:
+        return left.variable_id == right.variable_id
+
 
 def expressions_equal(
     left: ExpressionNode, right: ExpressionNode, abs_tol: float = 0, rel_tol: float = 0
diff --git a/src/gems_craft/expression/expression.py b/src/gems_craft/expression/expression.py
index 00fdc9eb..45d65da2 100644
--- a/src/gems_craft/expression/expression.py
+++ b/src/gems_craft/expression/expression.py
@@ -357,6 +357,16 @@ class ReducedCostNode(ExpressionNode):
     variable_id: str
 
 
+@dataclass(frozen=True, eq=False)
+class LowerBoundNode(ExpressionNode):
+    variable_id: str
+
+
+@dataclass(frozen=True, eq=False)
+class UpperBoundNode(ExpressionNode):
+    variable_id: str
+
+
 def sum_expressions(expressions: Sequence[ExpressionNode]) -> ExpressionNode:
     if len(expressions) == 0:
         return LiteralNode(0)
diff --git a/src/gems_craft/expression/indexing.py b/src/gems_craft/expression/indexing.py
index a74d6a7c..3168b15d 100644
--- a/src/gems_craft/expression/indexing.py
+++ b/src/gems_craft/expression/indexing.py
@@ -27,6 +27,7 @@
     ExpressionNode,
     FloorNode,
     LiteralNode,
+    LowerBoundNode,
     MaxNode,
     MinNode,
     MultiplicationNode,
@@ -40,6 +41,7 @@
     TimeEvalNode,
     TimeShiftNode,
     TimeSumNode,
+    UpperBoundNode,
     VariableNode,
 )
 from .visitor import ExpressionVisitor, T, visit
@@ -156,6 +158,12 @@ def dual(self, node: DualNode) -> IndexingStructure:
     def reduced_cost(self, node: ReducedCostNode) -> IndexingStructure:
         return self.context.get_variable_structure(node.variable_id)
 
+    def lower_bound(self, node: LowerBoundNode) -> IndexingStructure:
+        return self.context.get_variable_structure(node.variable_id)
+
+    def upper_bound(self, node: UpperBoundNode) -> IndexingStructure:
+        return self.context.get_variable_structure(node.variable_id)
+
 
 def compute_indexation(
     expression: ExpressionNode, provider: IndexingStructureProvider
diff --git a/src/gems_craft/expression/parsing/parse_expression.py b/src/gems_craft/expression/parsing/parse_expression.py
index 661ab69e..cddcad6d 100644
--- a/src/gems_craft/expression/parsing/parse_expression.py
+++ b/src/gems_craft/expression/parsing/parse_expression.py
@@ -21,9 +21,11 @@
     Comparator,
     ComparisonNode,
     DualNode,
+    LowerBoundNode,
     PortFieldAggregatorNode,
     PortFieldNode,
     ReducedCostNode,
+    UpperBoundNode,
     maximum,
     minimum,
 )
@@ -194,6 +196,22 @@ def _visit_reduced_cost(self, arg_exprs: list) -> ExpressionNode:
             raise ValueError(f"'{vid}' is not a variable of the model.")
         return ReducedCostNode(vid)
 
+    def _visit_lower_bound(self, arg_exprs: list) -> ExpressionNode:
+        if len(arg_exprs) != 1:
+            raise ValueError("lower_bound() requires exactly 1 argument.")
+        vid = arg_exprs[0].getText()  # type: ignore
+        if vid not in self.identifiers.variables:
+            raise ValueError(f"'{vid}' is not a variable of the model.")
+        return LowerBoundNode(vid)
+
+    def _visit_upper_bound(self, arg_exprs: list) -> ExpressionNode:
+        if len(arg_exprs) != 1:
+            raise ValueError("upper_bound() requires exactly 1 argument.")
+        vid = arg_exprs[0].getText()  # type: ignore
+        if vid not in self.identifiers.variables:
+            raise ValueError(f"'{vid}' is not a variable of the model.")
+        return UpperBoundNode(vid)
+
     # Visit a parse tree produced by ExprParser#function.
     def visitFunction(self, ctx: ExprParser.FunctionContext) -> ExpressionNode:
         function_name: str = ctx.IDENTIFIER().getText()  # type: ignore
@@ -204,6 +222,10 @@ def visitFunction(self, ctx: ExprParser.FunctionContext) -> ExpressionNode:
             return self._visit_dual(arg_exprs)
         if function_name == "reduced_cost":
             return self._visit_reduced_cost(arg_exprs)
+        if function_name == "lower_bound":
+            return self._visit_lower_bound(arg_exprs)
+        if function_name == "upper_bound":
+            return self._visit_upper_bound(arg_exprs)
 
         args: list[ExpressionNode] = (
             [expr.accept(self) for expr in arg_exprs]  # type: ignore
diff --git a/src/gems_craft/expression/print.py b/src/gems_craft/expression/print.py
index d90ef3c6..79ece5ce 100644
--- a/src/gems_craft/expression/print.py
+++ b/src/gems_craft/expression/print.py
@@ -20,6 +20,7 @@
     DualNode,
     ExpressionNode,
     FloorNode,
+    LowerBoundNode,
     MaxNode,
     MinNode,
     PortFieldAggregatorNode,
@@ -29,6 +30,7 @@
     TimeEvalNode,
     TimeShiftNode,
     TimeSumNode,
+    UpperBoundNode,
 )
 
 from .expression import (
@@ -144,6 +146,12 @@ def dual(self, node: DualNode) -> str:
     def reduced_cost(self, node: ReducedCostNode) -> str:
         return f"reduced_cost({node.variable_id})"
 
+    def lower_bound(self, node: LowerBoundNode) -> str:
+        return f"lower_bound({node.variable_id})"
+
+    def upper_bound(self, node: UpperBoundNode) -> str:
+        return f"upper_bound({node.variable_id})"
+
 
 def print_expr(expression: ExpressionNode) -> str:
     return visit(expression, PrinterVisitor())
diff --git a/src/gems_craft/expression/uses_sum_connections_on.py b/src/gems_craft/expression/uses_sum_connections_on.py
index 7824da28..0da90abe 100644
--- a/src/gems_craft/expression/uses_sum_connections_on.py
+++ b/src/gems_craft/expression/uses_sum_connections_on.py
@@ -21,6 +21,7 @@
     ExpressionNode,
     FloorNode,
     LiteralNode,
+    LowerBoundNode,
     MaxNode,
     MinNode,
     MultiplicationNode,
@@ -34,6 +35,7 @@
     TimeEvalNode,
     TimeShiftNode,
     TimeSumNode,
+    UpperBoundNode,
     VariableNode,
 )
 
@@ -122,6 +124,12 @@ def dual(self, node: DualNode) -> bool:
     def reduced_cost(self, node: ReducedCostNode) -> bool:
         return False
 
+    def lower_bound(self, node: LowerBoundNode) -> bool:
+        return False
+
+    def upper_bound(self, node: UpperBoundNode) -> bool:
+        return False
+
 
 def uses_sum_connections_on(
     expr: ExpressionNode, port_name: str, field_name: str
diff --git a/src/gems_craft/expression/visitor.py b/src/gems_craft/expression/visitor.py
index a7a15cc8..3db8c6d7 100644
--- a/src/gems_craft/expression/visitor.py
+++ b/src/gems_craft/expression/visitor.py
@@ -29,6 +29,7 @@
     ExpressionNode,
     FloorNode,
     LiteralNode,
+    LowerBoundNode,
     MaxNode,
     MinNode,
     MultiplicationNode,
@@ -42,6 +43,7 @@
     TimeEvalNode,
     TimeShiftNode,
     TimeSumNode,
+    UpperBoundNode,
     VariableNode,
 )
 
@@ -126,6 +128,12 @@ def dual(self, node: DualNode) -> T: ...
     @abstractmethod
     def reduced_cost(self, node: ReducedCostNode) -> T: ...
 
+    @abstractmethod
+    def lower_bound(self, node: LowerBoundNode) -> T: ...
+
+    @abstractmethod
+    def upper_bound(self, node: UpperBoundNode) -> T: ...
+
 
 def visit(root: ExpressionNode, visitor: ExpressionVisitor[T]) -> T:
     """
@@ -177,6 +185,10 @@ def visit(root: ExpressionNode, visitor: ExpressionVisitor[T]) -> T:
         return visitor.dual(root)
     elif isinstance(root, ReducedCostNode):
         return visitor.reduced_cost(root)
+    elif isinstance(root, LowerBoundNode):
+        return visitor.lower_bound(root)
+    elif isinstance(root, UpperBoundNode):
+        return visitor.upper_bound(root)
     raise ValueError(f"Unknown expression node type {root.__class__}")
 
 
diff --git a/src/gems_craft/model/port.py b/src/gems_craft/model/port.py
index bf6435d4..e46bf599 100644
--- a/src/gems_craft/model/port.py
+++ b/src/gems_craft/model/port.py
@@ -32,6 +32,7 @@
     CeilNode,
     DualNode,
     FloorNode,
+    LowerBoundNode,
     MaxNode,
     MinNode,
     PortFieldAggregatorNode,
@@ -42,6 +43,7 @@
     TimeEvalNode,
     TimeShiftNode,
     TimeSumNode,
+    UpperBoundNode,
 )
 from gems_craft.expression.visitor import visit
 
@@ -178,6 +180,12 @@ def dual(self, node: DualNode) -> None:
     def reduced_cost(self, node: ReducedCostNode) -> None:
         pass  # reduced_cost() is permitted in port-field definitions
 
+    def lower_bound(self, node: LowerBoundNode) -> None:
+        pass  # lower_bound() is permitted in port-field definitions
+
+    def upper_bound(self, node: UpperBoundNode) -> None:
+        pass  # upper_bound() is permitted in port-field definitions
+
 
 def _validate_port_field_expression(definition: PortFieldDefinition) -> None:
     visit(definition.definition, _PortFieldExpressionChecker())
diff --git a/src/gems_craft/model/resolve_library.py b/src/gems_craft/model/resolve_library.py
index 8603f484..f3ec74f1 100644
--- a/src/gems_craft/model/resolve_library.py
+++ b/src/gems_craft/model/resolve_library.py
@@ -23,6 +23,7 @@
     DualNode,
     FloorNode,
     LiteralNode,
+    LowerBoundNode,
     MaxNode,
     MinNode,
     MultiplicationNode,
@@ -36,6 +37,7 @@
     TimeEvalNode,
     TimeShiftNode,
     TimeSumNode,
+    UpperBoundNode,
     VariableNode,
 )
 from gems_craft.expression.indexing_structure import IndexingStructure
@@ -291,6 +293,12 @@ def dual(self, node: DualNode) -> None:
     def reduced_cost(self, node: ReducedCostNode) -> None:
         pass
 
+    def lower_bound(self, node: LowerBoundNode) -> None:
+        pass
+
+    def upper_bound(self, node: UpperBoundNode) -> None:
+        pass
+
 
 def _forbid_bare_port_field(expr: ExpressionNode, context: str) -> None:
     visit(expr, _ForbidBarePortFieldVisitor(context))
diff --git a/src/gems_runner/expression/evaluate.py b/src/gems_runner/expression/evaluate.py
index c6cf5b97..0219db19 100644
--- a/src/gems_runner/expression/evaluate.py
+++ b/src/gems_runner/expression/evaluate.py
@@ -24,6 +24,7 @@
     ExpressionNode,
     FloorNode,
     LiteralNode,
+    LowerBoundNode,
     MaxNode,
     MinNode,
     ParameterNode,
@@ -35,6 +36,7 @@
     TimeEvalNode,
     TimeShiftNode,
     TimeSumNode,
+    UpperBoundNode,
     VariableNode,
 )
 from gems_craft.expression.indexing import IndexingStructureProvider
@@ -139,6 +141,12 @@ def dual(self, node: DualNode) -> float:
     def reduced_cost(self, node: ReducedCostNode) -> float:
         raise NotImplementedError("reduced_cost() is not statically evaluable.")
 
+    def lower_bound(self, node: LowerBoundNode) -> float:
+        raise NotImplementedError("lower_bound() is not statically evaluable.")
+
+    def upper_bound(self, node: UpperBoundNode) -> float:
+        raise NotImplementedError("upper_bound() is not statically evaluable.")
+
 
 def evaluate(expression: ExpressionNode, value_provider: ValueProvider) -> float:
     return visit(expression, EvaluationVisitor(value_provider))
diff --git a/src/gems_runner/simulation/extra_output.py b/src/gems_runner/simulation/extra_output.py
index 7bc5e988..d0deb5ec 100644
--- a/src/gems_runner/simulation/extra_output.py
+++ b/src/gems_runner/simulation/extra_output.py
@@ -35,7 +35,9 @@
     Comparator,
     ComparisonNode,
     DualNode,
+    LowerBoundNode,
     ReducedCostNode,
+    UpperBoundNode,
     VariableNode,
 )
 from gems_craft.expression.visitor import visit
@@ -110,6 +112,14 @@ class VectorizedExtraOutputBuilder(VectorizedBuilderBase[xr.DataArray]):
     var_reduced_cost_arrays:
         Mapping from (model_id, var_name) to a DataArray of reduced cost values,
         with dims in {component, time, scenario} (or a subset).
+    var_lower_bound_arrays:
+        Mapping from (model_id, var_name) to a DataArray of the variable's
+        current lower bound, with dims in {component, time, scenario} (or a
+        subset).
+    var_upper_bound_arrays:
+        Mapping from (model_id, var_name) to a DataArray of the variable's
+        current upper bound, with dims in {component, time, scenario} (or a
+        subset).
     port_arrays:
         Pre-computed xr.DataArray for each PortFieldId of this model.
         Keyed by PortFieldId(port_name, field_name).
@@ -124,6 +134,12 @@ class VectorizedExtraOutputBuilder(VectorizedBuilderBase[xr.DataArray]):
     var_reduced_cost_arrays: Dict[Tuple[str, str], xr.DataArray] = field(
         default_factory=dict
     )
+    var_lower_bound_arrays: Dict[Tuple[str, str], xr.DataArray] = field(
+        default_factory=dict
+    )
+    var_upper_bound_arrays: Dict[Tuple[str, str], xr.DataArray] = field(
+        default_factory=dict
+    )
 
     def variable(self, node: VariableNode) -> xr.DataArray:
         key = (self.model_id, node.name)
@@ -152,6 +168,24 @@ def reduced_cost(self, node: ReducedCostNode) -> xr.DataArray:
             )
         return self.var_reduced_cost_arrays[key]
 
+    def lower_bound(self, node: LowerBoundNode) -> xr.DataArray:
+        key = (self.model_id, node.variable_id)
+        if key not in self.var_lower_bound_arrays:
+            raise KeyError(
+                f"Lower bound of variable '{node.variable_id}' not found for model "
+                f"{self.model_id!r}."
+            )
+        return self.var_lower_bound_arrays[key]
+
+    def upper_bound(self, node: UpperBoundNode) -> xr.DataArray:
+        key = (self.model_id, node.variable_id)
+        if key not in self.var_upper_bound_arrays:
+            raise KeyError(
+                f"Upper bound of variable '{node.variable_id}' not found for model "
+                f"{self.model_id!r}."
+            )
+        return self.var_upper_bound_arrays[key]
+
     def comparison(self, node: ComparisonNode) -> xr.DataArray:
         """Evaluate a comparison post-solve as a float indicator DataArray.
 
diff --git a/src/gems_runner/simulation/optimization.py b/src/gems_runner/simulation/optimization.py
index 6eb4ff61..cc1905ba 100644
--- a/src/gems_runner/simulation/optimization.py
+++ b/src/gems_runner/simulation/optimization.py
@@ -433,18 +433,22 @@ def get_variable_labels(
         lv = self._linopy_vars.get((model_id, var_name))
         return lv.labels if lv is not None else None
 
-    def get_variable_solution(
-        self, model_id: str, var_name: str
+    def _reassemble_variable_attr(
+        self, model_id: str, var_name: str, attr: str
     ) -> Optional[xr.DataArray]:
-        """Return solved values for *var_name* across all its components.
-
-        Unlike ``linopy_model.solution[<name>]``, this is correct even when the
-        variable was split across relaxed/exact strategy groups: the merged
-        ``_linopy_vars`` copy keeps the ``.name`` of only one of the two really
-        -registered group Variables, so indexing the solver's solution Dataset
-        by that name silently drops the other group's components. This instead
-        reads ``.solution`` directly off the real per-component Variables
-        (``_linopy_vars_by_component``) and reassembles them.
+        """Reassemble a per-instance ``linopy.Variable`` attribute (``solution``,
+        ``lower``, ``upper``) for *var_name* across all its components.
+
+        Unlike reading the attribute off ``self._linopy_vars[(model_id, var_name)]``,
+        this is correct even when the variable was split across relaxed/exact
+        strategy groups: that merged, detached ``_MergedGroupVariable`` copy
+        keeps the ``.name`` of only one of the two really-registered group
+        Variables (so indexing the solver's solution Dataset by that name
+        silently drops the other group's components), and it is rebuilt once at
+        problem-build time so it never reflects later bound mutations (e.g. from
+        heuristics). This instead reads the attribute directly off the real
+        per-component Variables (``_linopy_vars_by_component``) and reassembles
+        them.
         """
         by_name: Dict[str, linopy.Variable] = {}
         for (m, vn, _c), variable in self._linopy_vars_by_component.items():
@@ -456,11 +460,46 @@ def get_variable_solution(
         if not group_vars:
             return None
         if len(group_vars) == 1:
-            return group_vars[0].solution
+            return cast(xr.DataArray, getattr(group_vars[0], attr))
         return cast(
-            xr.DataArray, xr.concat([v.solution for v in group_vars], dim="component")
+            xr.DataArray,
+            xr.concat([getattr(v, attr) for v in group_vars], dim="component"),
         )
 
+    def get_variable_solution(
+        self, model_id: str, var_name: str
+    ) -> Optional[xr.DataArray]:
+        """Return solved values for *var_name* across all its components.
+
+        See :meth:`_reassemble_variable_attr` for why this must bypass the
+        merged ``_linopy_vars`` copy.
+        """
+        return self._reassemble_variable_attr(model_id, var_name, "solution")
+
+    def get_variable_lower_bound(
+        self, model_id: str, var_name: str
+    ) -> Optional[xr.DataArray]:
+        """Return the current lower bound for *var_name* across all its components.
+
+        Reflects any bound mutation applied after problem construction (e.g. by
+        thermal heuristics via :meth:`get_component_variable`). See
+        :meth:`_reassemble_variable_attr` for why this must bypass the merged
+        ``_linopy_vars`` copy.
+        """
+        return self._reassemble_variable_attr(model_id, var_name, "lower")
+
+    def get_variable_upper_bound(
+        self, model_id: str, var_name: str
+    ) -> Optional[xr.DataArray]:
+        """Return the current upper bound for *var_name* across all its components.
+
+        Reflects any bound mutation applied after problem construction (e.g. by
+        thermal heuristics via :meth:`get_component_variable`). See
+        :meth:`_reassemble_variable_attr` for why this must bypass the merged
+        ``_linopy_vars`` copy.
+        """
+        return self._reassemble_variable_attr(model_id, var_name, "upper")
+
 
 # ---------------------------------------------------------------------------
 # Internal builder
diff --git a/src/gems_runner/simulation/simulation_table.py b/src/gems_runner/simulation/simulation_table.py
index ca993644..e3dd1559 100644
--- a/src/gems_runner/simulation/simulation_table.py
+++ b/src/gems_runner/simulation/simulation_table.py
@@ -279,6 +279,8 @@ def _collect_extra_outputs(
 
         constraint_dual_arrays = self._collect_constraint_duals(problem)
         var_reduced_cost_arrays = self._collect_reduced_costs(problem)
+        var_lower_bound_arrays = self._collect_lower_bounds(problem)
+        var_upper_bound_arrays = self._collect_upper_bounds(problem)
 
         for mk, components in problem.study.model_components.items():
             model = problem.study.models[mk]
@@ -295,6 +297,8 @@ def _collect_extra_outputs(
                     var_solution_arrays=var_solution_arrays,
                     constraint_dual_arrays=constraint_dual_arrays,
                     var_reduced_cost_arrays=var_reduced_cost_arrays,
+                    var_lower_bound_arrays=var_lower_bound_arrays,
+                    var_upper_bound_arrays=var_upper_bound_arrays,
                     port_arrays={},
                     block_length=problem.block_length,
                 ),
@@ -307,6 +311,8 @@ def _collect_extra_outputs(
                     var_solution_arrays=var_solution_arrays,
                     constraint_dual_arrays=constraint_dual_arrays,
                     var_reduced_cost_arrays=var_reduced_cost_arrays,
+                    var_lower_bound_arrays=var_lower_bound_arrays,
+                    var_upper_bound_arrays=var_upper_bound_arrays,
                     port_arrays=port_arrays,
                     block_length=problem.block_length,
                 )
@@ -411,6 +417,30 @@ def _collect_reduced_costs(
         except Exception:
             return {}
 
+    @staticmethod
+    def _collect_lower_bounds(
+        problem: OptimizationProblem,
+    ) -> Dict[Tuple[str, str], xr.DataArray]:
+        """Return current variable lower bounds keyed by (model_key, var_name)."""
+        result: Dict[Tuple[str, str], xr.DataArray] = {}
+        for mk, vname in problem._linopy_vars:
+            lb = problem.get_variable_lower_bound(mk, vname)
+            if lb is not None:
+                result[(mk, vname)] = lb
+        return result
+
+    @staticmethod
+    def _collect_upper_bounds(
+        problem: OptimizationProblem,
+    ) -> Dict[Tuple[str, str], xr.DataArray]:
+        """Return current variable upper bounds keyed by (model_key, var_name)."""
+        result: Dict[Tuple[str, str], xr.DataArray] = {}
+        for mk, vname in problem._linopy_vars:
+            ub = problem.get_variable_upper_bound(mk, vname)
+            if ub is not None:
+                result[(mk, vname)] = ub
+        return result
+
     # -------------------------------------------------------------------------
     # Objective value
     # -------------------------------------------------------------------------
diff --git a/src/gems_runner/simulation/vectorized_builder.py b/src/gems_runner/simulation/vectorized_builder.py
index 137bcc47..f6a831f6 100644
--- a/src/gems_runner/simulation/vectorized_builder.py
+++ b/src/gems_runner/simulation/vectorized_builder.py
@@ -51,6 +51,7 @@
     ExpressionNode,
     FloorNode,
     LiteralNode,
+    LowerBoundNode,
     MaxNode,
     MinNode,
     MultiplicationNode,
@@ -64,6 +65,7 @@
     TimeEvalNode,
     TimeShiftNode,
     TimeSumNode,
+    UpperBoundNode,
     VariableNode,
 )
 from gems_craft.expression.visitor import (
@@ -406,6 +408,18 @@ def reduced_cost(self, node: ReducedCostNode) -> VectorizedExpr:
             f"not in {type(self).__name__}."
         )
 
+    def lower_bound(self, node: LowerBoundNode) -> VectorizedExpr:
+        raise NotImplementedError(
+            f"lower_bound() is only available in the extra-output builder, "
+            f"not in {type(self).__name__}."
+        )
+
+    def upper_bound(self, node: UpperBoundNode) -> VectorizedExpr:
+        raise NotImplementedError(
+            f"upper_bound() is only available in the extra-output builder, "
+            f"not in {type(self).__name__}."
+        )
+
     # ------------------------------------------------------------------ #
     # Private helpers                                                       #
     # ------------------------------------------------------------------ #
@@ -573,6 +587,12 @@ def dual(self, node: DualNode) -> xr.DataArray:
     def reduced_cost(self, node: ReducedCostNode) -> xr.DataArray:
         raise NotImplementedError
 
+    def lower_bound(self, node: LowerBoundNode) -> xr.DataArray:
+        raise NotImplementedError
+
+    def upper_bound(self, node: UpperBoundNode) -> xr.DataArray:
+        raise NotImplementedError
+
 
 def _and_mask(
     a: Optional[xr.DataArray], b: Optional[xr.DataArray]
@@ -741,3 +761,9 @@ def dual(self, node: DualNode) -> Optional[xr.DataArray]:
 
     def reduced_cost(self, node: ReducedCostNode) -> Optional[xr.DataArray]:
         return None
+
+    def lower_bound(self, node: LowerBoundNode) -> Optional[xr.DataArray]:
+        return None
+
+    def upper_bound(self, node: UpperBoundNode) -> Optional[xr.DataArray]:
+        return None
diff --git a/tests/e2e/functional/libs/thermal_variants_for_heuristic.yml b/tests/e2e/functional/libs/thermal_variants_for_heuristic.yml
index 95e5d8de..bd7011c2 100644
--- a/tests/e2e/functional/libs/thermal_variants_for_heuristic.yml
+++ b/tests/e2e/functional/libs/thermal_variants_for_heuristic.yml
@@ -154,6 +154,6 @@ library:
           expression: expec(sum(market_bid_cost * generation_power))
       extra-outputs:
         - id: num_units_on
-          expression: ceil(generation_power / max_power_per_unit)
+          expression: max(floor(lower_bound(generation_power) / min_power_per_unit), ceil(generation_power / max_power_per_unit / (1-spinning/100)))
         - id: non_prop_cost
-          expression: startup_cost * max(0,ceil(generation_power / max_power_per_unit)-(ceil(generation_power / max_power_per_unit))[t-1])+ fixed_cost * ceil(generation_power / max_power_per_unit)
+          expression: startup_cost * max(0,max(floor(lower_bound(generation_power) / min_power_per_unit), ceil(generation_power / max_power_per_unit / (1-spinning/100)))-(max(floor(lower_bound(generation_power) / min_power_per_unit), ceil(generation_power / max_power_per_unit / (1-spinning/100))))[t-1])+ fixed_cost * max(floor(lower_bound(generation_power) / min_power_per_unit), ceil(generation_power / max_power_per_unit / (1-spinning/100)))
diff --git a/tests/e2e/functional/test_thermal_heuristic_four_clusters_mixed_strategies.py b/tests/e2e/functional/test_thermal_heuristic_four_clusters_mixed_strategies.py
index 00ceebcb..b8c5a97c 100644
--- a/tests/e2e/functional/test_thermal_heuristic_four_clusters_mixed_strategies.py
+++ b/tests/e2e/functional/test_thermal_heuristic_four_clusters_mixed_strategies.py
@@ -55,7 +55,7 @@
 _G3_GEN = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 100.0, 328.0, 503.0, 563.0, 600.0, 563.0, 503.0, 328.0, 100.0, 0.0, 0.0, 0.0, 0.0]
 _G4_GEN = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 50.0, 50.0, 50.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
 
-_G2_NODU = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 2.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 2.0, 1.0, 1.0]
+_G2_NODU = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 2.0, 2.0, 2.0]
 _G3_NODU = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 2.0, 3.0, 3.0, 3.0, 3.0, 3.0, 2.0, 1.0, 0.0, 0.0, 0.0, 0.0]
 _G4_NODU = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
 # fmt: on
@@ -90,5 +90,5 @@ def test_mixed_integer_strategies() -> None:
     check_output(st, "N", "spilled_energy", [0.0] * 24)
 
     assert total_output_sum(st, THERMAL_COMPONENTS, "non_prop_cost") == pytest.approx(
-        13662
+        13665
     )
diff --git a/tests/e2e/functional/test_thermal_heuristic_two_clusters_low_load.py b/tests/e2e/functional/test_thermal_heuristic_two_clusters_low_load.py
index 6ae6c433..dd641726 100644
--- a/tests/e2e/functional/test_thermal_heuristic_two_clusters_low_load.py
+++ b/tests/e2e/functional/test_thermal_heuristic_two_clusters_low_load.py
@@ -53,7 +53,7 @@
 
 _G1_GEN_FAST = [24600.0, 24600.0, 24600.0, 24600.0, 24600.0, 24600.0, 29152.0, 28840.0, 30246.0, 29912.0, 27437.0, 31543.0, 38677.0, 36685.0, 31115.0, 22643.0, 17874.0, 18589.0, 19041.0, 16725.0, 15600.0, 15600.0, 15600.0, 6724.0, 4800.0, 4800.0, 4800.0, 4800.0, 4800.0, 4800.0, 4800.0, 9000.0, 9000.0, 9000.0, 9000.0, 9000.0, 13147.0, 12872.0, 9425.0, 10200.0, 10200.0, 10200.0, 11457.0, 11380.0, 12093.0, 12033.0, 14753.0, 12000.0, 12000.0, 12827.0, 12000.0, 12000.0, 12000.0, 12172.0, 17738.0, 19634.0, 19837.0, 17800.0, 13800.0, 16393.0, 19334.0, 17565.0, 14124.0, 14137.0, 13316.0, 16819.0, 19358.0, 15619.0, 13200.0, 13200.0, 13200.0, 6600.0, 6600.0, 6600.0, 6600.0, 6600.0, 6600.0, 6600.0, 9673.0, 15000.0, 15000.0, 15000.0, 15000.0, 16806.0, 21925.0, 20016.0, 16679.0, 14606.0, 13176.0, 15619.0, 17136.0, 14818.0, 12000.0, 12000.0, 12000.0, 5400.0, 5400.0, 5400.0, 5400.0, 5400.0, 5400.0, 5400.0, 7893.0, 18000.0, 18000.0, 18000.0, 18000.0, 21118.0, 26935.0, 26345.0, 24368.0, 22581.0, 24429.0, 31131.0, 31414.0, 28160.0, 24430.0, 22302.0, 23479.0, 20638.0, 14400.0, 14400.0, 14400.0, 14400.0, 14400.0, 14400.0, 21222.0, 28128.0, 29724.0, 29975.0, 24316.0, 27662.0, 32183.0, 30524.0, 30265.0, 31729.0, 35121.0, 40762.0, 44185.0, 45000.0, 41891.0, 41048.0, 41509.0, 39384.0, 34533.0, 34661.0, 33049.0, 31777.0, 31616.0, 31693.0, 38073.0, 39023.0, 38405.0, 36189.0, 30393.0, 35535.0, 42843.0, 43891.0, 42826.0, 43974.0, 45000.0, 45000.0, 45000.0, 45000.0, 42494.0, 38634.0, 39132.0, 36502.0]
 _G2_GEN_FAST = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 67.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1937.0, 5358.0, 5539.0, 3045.0, 0.0, 0.0, 0.0, 0.0]
-_G1_NODU_FAST = [41.0, 41.0, 41.0, 41.0, 41.0, 41.0, 41.0, 43.0, 43.0, 43.0, 43.0, 43.0, 43.0, 43.0, 43.0, 26.0, 26.0, 26.0, 26.0, 26.0, 26.0, 26.0, 26.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 15.0, 15.0, 15.0, 15.0, 15.0, 15.0, 15.0, 15.0, 17.0, 17.0, 17.0, 17.0, 17.0, 17.0, 17.0, 17.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 23.0, 23.0, 23.0, 23.0, 23.0, 23.0, 23.0, 23.0, 22.0, 22.0, 22.0, 22.0, 22.0, 22.0, 22.0, 22.0, 11.0, 11.0, 11.0, 11.0, 11.0, 11.0, 11.0, 11.0, 25.0, 25.0, 25.0, 25.0, 25.0, 25.0, 25.0, 25.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 9.0, 9.0, 9.0, 9.0, 9.0, 9.0, 9.0, 9.0, 30.0, 30.0, 30.0, 30.0, 30.0, 30.0, 30.0, 30.0, 35.0, 35.0, 35.0, 35.0, 35.0, 35.0, 35.0, 35.0, 24.0, 24.0, 24.0, 24.0, 24.0, 24.0, 24.0, 24.0, 36.0, 36.0, 36.0, 36.0, 36.0, 36.0, 36.0, 36.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 41.0]
+_G1_NODU_FAST = [41.0, 41.0, 41.0, 41.0, 41.0, 41.0, 41.0, 43.0, 43.0, 43.0, 43.0, 43.0, 43.0, 43.0, 43.0, 26.0, 26.0, 26.0, 26.0, 26.0, 26.0, 26.0, 26.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 15.0, 15.0, 15.0, 15.0, 15.0, 15.0, 15.0, 15.0, 17.0, 17.0, 17.0, 17.0, 17.0, 17.0, 17.0, 17.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 23.0, 23.0, 23.0, 23.0, 23.0, 23.0, 23.0, 23.0, 22.0, 22.0, 22.0, 22.0, 22.0, 22.0, 22.0, 22.0, 11.0, 11.0, 11.0, 11.0, 11.0, 11.0, 11.0, 11.0, 25.0, 25.0, 25.0, 25.0, 25.0, 25.0, 25.0, 25.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 20.0, 9.0, 9.0, 9.0, 9.0, 9.0, 9.0, 9.0, 9.0, 30.0, 30.0, 30.0, 30.0, 30.0, 30.0, 30.0, 30.0, 35.0, 35.0, 35.0, 35.0, 35.0, 35.0, 35.0, 35.0, 24.0, 24.0, 24.0, 24.0, 24.0, 24.0, 24.0, 24.0, 36.0, 36.0, 36.0, 36.0, 36.0, 36.0, 36.0, 36.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 44.0, 44.0, 44.0, 44.0, 44.0, 44.0, 44.0, 44.0, 49.0, 49.0, 49.0, 49.0, 49.0, 49.0, 49.0, 49.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 41.0]
 _G2_NODU_FAST = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.0, 6.0, 6.0, 4.0, 0.0, 0.0, 0.0, 0.0]
 _SPIL_FAST = [4811.0, 4163.0, 4503.0, 2473.0, 697.0, 800.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3102.0, 4107.0, 4344.0, 0.0, 4491.0, 4673.0, 4798.0, 4706.0, 4339.0, 4283.0, 4079.0, 8019.0, 7483.0, 6327.0, 5720.0, 958.0, 0.0, 0.0, 0.0, 3930.0, 4173.0, 400.0, 0.0, 0.0, 0.0, 0.0, 0.0, 939.0, 333.0, 0.0, 1071.0, 1791.0, 2358.0, 0.0, 0.0, 0.0, 0.0, 0.0, 397.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3562.0, 7636.0, 7958.0, 4913.0, 5642.0, 6156.0, 6502.0, 6368.0, 6169.0, 4857.0, 0.0, 1405.0, 1373.0, 1983.0, 3635.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 826.0, 3465.0, 4101.0, 1176.0, 3191.0, 4589.0, 5003.0, 4510.0, 4611.0, 4778.0, 0.0, 4926.0, 4174.0, 3533.0, 5550.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3122.0, 4855.0, 7280.0, 7215.0, 5066.0, 1475.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
 _UNSP_FAST = [0.0] * 168
@@ -106,7 +106,6 @@ def test_accurate_heuristic() -> None:
         2641690
     )
 
-@pytest.mark.xfail(reason="Computation of num_units_on with fast heuristic not implemented.")
 def test_fast_heuristic() -> None:
     """
     Solve the same problem with the fast heuristic.
@@ -127,5 +126,5 @@ def test_fast_heuristic() -> None:
     check_output(st, "N", "spilled_energy", _SPIL_FAST)
     assert sum(_SPIL_FAST) == pytest.approx(255873)
     assert total_output_sum(st, THERMAL_COMPONENTS, "non_prop_cost") == pytest.approx(
-        4640150
+        2273100
     )
diff --git a/tests/unittests/gems_craft/expressions/parsing/test_expression_parsing.py b/tests/unittests/gems_craft/expressions/parsing/test_expression_parsing.py
index ee430edf..a8c253ec 100644
--- a/tests/unittests/gems_craft/expressions/parsing/test_expression_parsing.py
+++ b/tests/unittests/gems_craft/expressions/parsing/test_expression_parsing.py
@@ -17,7 +17,9 @@
 from gems_craft.expression.equality import expressions_equal
 from gems_craft.expression.expression import (
     DualNode,
+    LowerBoundNode,
     ReducedCostNode,
+    UpperBoundNode,
     maximum,
     minimum,
     port_field,
@@ -251,6 +253,8 @@ def test_parsing_visitor(
     [
         (set(), set(), {"balance"}, "dual(balance)", DualNode("balance")),
         ({"p"}, set(), set(), "reduced_cost(p)", ReducedCostNode("p")),
+        ({"p"}, set(), set(), "lower_bound(p)", LowerBoundNode("p")),
+        ({"p"}, set(), set(), "upper_bound(p)", UpperBoundNode("p")),
     ],
 )
 def test_parsing_dual_and_reduced_cost(
@@ -277,6 +281,18 @@ def test_parse_reduced_cost_unknown_variable_raises() -> None:
         parse_expression("reduced_cost(p)", identifiers)
 
 
+def test_parse_lower_bound_unknown_variable_raises() -> None:
+    identifiers = ModelIdentifiers({"x"}, set(), set())
+    with pytest.raises(ParsingException, match="not a variable"):
+        parse_expression("lower_bound(p)", identifiers)
+
+
+def test_parse_upper_bound_unknown_variable_raises() -> None:
+    identifiers = ModelIdentifiers({"x"}, set(), set())
+    with pytest.raises(ParsingException, match="not a variable"):
+        parse_expression("upper_bound(p)", identifiers)
+
+
 @pytest.mark.parametrize(
     "expression_str",
     [
diff --git a/tests/unittests/gems_craft/expressions/visitor/test_degree.py b/tests/unittests/gems_craft/expressions/visitor/test_degree.py
index 0bceccea..0ceb958c 100644
--- a/tests/unittests/gems_craft/expressions/visitor/test_degree.py
+++ b/tests/unittests/gems_craft/expressions/visitor/test_degree.py
@@ -28,8 +28,10 @@
     CeilNode,
     DualNode,
     FloorNode,
+    LowerBoundNode,
     ReducedCostNode,
     RoundNode,
+    UpperBoundNode,
 )
 
 
@@ -87,6 +89,11 @@ def test_dual_reduced_cost_degree() -> None:
     assert visit(ReducedCostNode("p"), ExpressionDegreeVisitor()) == math.inf
 
 
+def test_lower_upper_bound_degree() -> None:
+    assert visit(LowerBoundNode("x"), ExpressionDegreeVisitor()) == math.inf
+    assert visit(UpperBoundNode("x"), ExpressionDegreeVisitor()) == math.inf
+
+
 @pytest.mark.xfail(reason="Degree simplification not implemented")
 def test_degree_computation_should_take_into_account_simplifications() -> None:
     x = var("x")
diff --git a/tests/unittests/gems_craft/expressions/visitor/test_equality.py b/tests/unittests/gems_craft/expressions/visitor/test_equality.py
index 0c831a73..9fdae1d4 100644
--- a/tests/unittests/gems_craft/expressions/visitor/test_equality.py
+++ b/tests/unittests/gems_craft/expressions/visitor/test_equality.py
@@ -14,7 +14,14 @@
 
 from gems_craft.expression import ExpressionNode, copy_expression, literal, param, var
 from gems_craft.expression.equality import expressions_equal
-from gems_craft.expression.expression import DualNode, ReducedCostNode, maximum, minimum
+from gems_craft.expression.expression import (
+    DualNode,
+    LowerBoundNode,
+    ReducedCostNode,
+    UpperBoundNode,
+    maximum,
+    minimum,
+)
 
 
 @pytest.mark.parametrize(
@@ -94,3 +101,11 @@ def test_dual_reduced_cost_equality() -> None:
     )
     assert not expressions_equal(ReducedCostNode("p"), ReducedCostNode("q"))
     assert not expressions_equal(DualNode("p"), ReducedCostNode("p"))
+
+
+def test_lower_upper_bound_equality() -> None:
+    assert expressions_equal(LowerBoundNode("x"), copy_expression(LowerBoundNode("x")))
+    assert not expressions_equal(LowerBoundNode("x"), LowerBoundNode("y"))
+    assert expressions_equal(UpperBoundNode("x"), copy_expression(UpperBoundNode("x")))
+    assert not expressions_equal(UpperBoundNode("x"), UpperBoundNode("y"))
+    assert not expressions_equal(LowerBoundNode("x"), UpperBoundNode("x"))
diff --git a/tests/unittests/gems_craft/expressions/visitor/test_indexing.py b/tests/unittests/gems_craft/expressions/visitor/test_indexing.py
index bf6aa7fa..84e94943 100644
--- a/tests/unittests/gems_craft/expressions/visitor/test_indexing.py
+++ b/tests/unittests/gems_craft/expressions/visitor/test_indexing.py
@@ -12,7 +12,12 @@
 
 
 from gems_craft.expression import param, var
-from gems_craft.expression.expression import DualNode, ReducedCostNode
+from gems_craft.expression.expression import (
+    DualNode,
+    LowerBoundNode,
+    ReducedCostNode,
+    UpperBoundNode,
+)
 from gems_craft.expression.indexing import IndexingStructureProvider, compute_indexation
 from gems_craft.expression.indexing_structure import IndexingStructure
 
@@ -122,3 +127,13 @@ def test_dual_reduced_cost_indexing() -> None:
     assert compute_indexation(ReducedCostNode("p"), provider) == IndexingStructure(
         True, True
     )
+
+
+def test_lower_upper_bound_indexing() -> None:
+    provider = StructureProvider()
+    assert compute_indexation(LowerBoundNode("x"), provider) == IndexingStructure(
+        True, True
+    )
+    assert compute_indexation(UpperBoundNode("x"), provider) == IndexingStructure(
+        True, True
+    )
diff --git a/tests/unittests/gems_craft/expressions/visitor/test_printer.py b/tests/unittests/gems_craft/expressions/visitor/test_printer.py
index 2b6d71de..37d406e7 100644
--- a/tests/unittests/gems_craft/expressions/visitor/test_printer.py
+++ b/tests/unittests/gems_craft/expressions/visitor/test_printer.py
@@ -11,7 +11,12 @@
 # This file is part of the Antares project.
 
 from gems_craft.expression import ExpressionNode, PrinterVisitor, param, var, visit
-from gems_craft.expression.expression import DualNode, ReducedCostNode
+from gems_craft.expression.expression import (
+    DualNode,
+    LowerBoundNode,
+    ReducedCostNode,
+    UpperBoundNode,
+)
 
 
 def test_comparison() -> None:
@@ -55,3 +60,8 @@ def test_abs_round_printer() -> None:
 def test_dual_reduced_cost_printer() -> None:
     assert visit(DualNode("balance"), PrinterVisitor()) == "dual(balance)"
     assert visit(ReducedCostNode("p"), PrinterVisitor()) == "reduced_cost(p)"
+
+
+def test_lower_upper_bound_printer() -> None:
+    assert visit(LowerBoundNode("x"), PrinterVisitor()) == "lower_bound(x)"
+    assert visit(UpperBoundNode("x"), PrinterVisitor()) == "upper_bound(x)"
diff --git a/tests/unittests/gems_craft/lib_parsing/test_lib_parsing.py b/tests/unittests/gems_craft/lib_parsing/test_lib_parsing.py
index e3930fb8..a2c9e08f 100644
--- a/tests/unittests/gems_craft/lib_parsing/test_lib_parsing.py
+++ b/tests/unittests/gems_craft/lib_parsing/test_lib_parsing.py
@@ -21,7 +21,9 @@
 )
 from gems_craft.expression.expression import (
     DualNode,
+    LowerBoundNode,
     ReducedCostNode,
+    UpperBoundNode,
     maximum,
     minimum,
     port_field,
@@ -307,6 +309,8 @@ def test_reduced_cost_in_objective_is_rejected() -> None:
 _UNRESTRICTED_EXPRS = [
     pytest.param("dual(balance)", DualNode("balance"), id="dual"),
     pytest.param("reduced_cost(x)", ReducedCostNode("x"), id="reduced_cost"),
+    pytest.param("lower_bound(x)", LowerBoundNode("x"), id="lower_bound"),
+    pytest.param("upper_bound(x)", UpperBoundNode("x"), id="upper_bound"),
     pytest.param("max(x, y)", maximum(var("x"), var("y")), id="max"),
     pytest.param("min(x, y)", minimum(var("x"), var("y")), id="min"),
     pytest.param("abs(x)", var("x").abs(), id="abs"),
@@ -654,3 +658,25 @@ def test_sum_connections_on_non_own_port_accepted() -> None:
         )
     )
     resolve_library([input_lib])  # must not raise
+
+
+# ---------------------------------------------------------------------------
+# Acceptance: lower_bound()/upper_bound() parse from YAML inside extra-outputs
+# ---------------------------------------------------------------------------
+
+
+def test_lower_upper_bound_in_extra_output_accepted() -> None:
+    """lower_bound(generation)/upper_bound(generation) are valid extra-output expressions."""
+    input_lib = parse_yaml_library(
+        io.StringIO(
+            _port_model_yaml(
+                extra_output_expr="lower_bound(generation) + upper_bound(generation)"
+            )
+        )
+    )
+    lib = resolve_library([input_lib])
+    eo = lib["test"].models["test.gen_model"].extra_outputs
+    assert eo is not None
+    assert expressions_equal(
+        eo["eo"], LowerBoundNode("generation") + UpperBoundNode("generation")
+    )
diff --git a/tests/unittests/gems_runner/expression/test_evaluation.py b/tests/unittests/gems_runner/expression/test_evaluation.py
index 72feee68..785ecb9b 100644
--- a/tests/unittests/gems_runner/expression/test_evaluation.py
+++ b/tests/unittests/gems_runner/expression/test_evaluation.py
@@ -27,7 +27,12 @@
     visit,
 )
 from gems_craft.expression.equality import expressions_equal
-from gems_craft.expression.expression import DualNode, ReducedCostNode
+from gems_craft.expression.expression import (
+    DualNode,
+    LowerBoundNode,
+    ReducedCostNode,
+    UpperBoundNode,
+)
 from gems_runner.expression import EvaluationContext, EvaluationVisitor, ValueProvider
 
 
@@ -110,3 +115,11 @@ def test_dual_reduced_cost_evaluation_raises() -> None:
         visit(DualNode("balance"), EvaluationVisitor(ctx))
     with pytest.raises(NotImplementedError, match="reduced_cost"):
         visit(ReducedCostNode("p"), EvaluationVisitor(ctx))
+
+
+def test_lower_upper_bound_evaluation_raises() -> None:
+    ctx = EvaluationContext()
+    with pytest.raises(NotImplementedError, match="lower_bound"):
+        visit(LowerBoundNode("x"), EvaluationVisitor(ctx))
+    with pytest.raises(NotImplementedError, match="upper_bound"):
+        visit(UpperBoundNode("x"), EvaluationVisitor(ctx))
diff --git a/tests/unittests/gems_runner/simulation/simulation_table_fakes.py b/tests/unittests/gems_runner/simulation/simulation_table_fakes.py
index 89027647..fea1a83b 100644
--- a/tests/unittests/gems_runner/simulation/simulation_table_fakes.py
+++ b/tests/unittests/gems_runner/simulation/simulation_table_fakes.py
@@ -78,6 +78,16 @@ def get_variable_solution(
             return None
         return self.linopy_model.solution.get(lv.name)
 
+    def get_variable_lower_bound(
+        self, model_id: object, var_name: str
+    ) -> Optional[xr.DataArray]:
+        return None
+
+    def get_variable_upper_bound(
+        self, model_id: object, var_name: str
+    ) -> Optional[xr.DataArray]:
+        return None
+
 
 def to_object_dtype(frame: pd.DataFrame) -> pd.DataFrame:
     """Cast every column to numpy object dtype, normalising all nulls to None."""
diff --git a/tests/unittests/gems_runner/simulation/test_integer_strategy.py b/tests/unittests/gems_runner/simulation/test_integer_strategy.py
index e9432ea1..9f4ada13 100644
--- a/tests/unittests/gems_runner/simulation/test_integer_strategy.py
+++ b/tests/unittests/gems_runner/simulation/test_integer_strategy.py
@@ -19,8 +19,9 @@
 """
 
 import pandas as pd
+import pytest
 
-from gems_craft.expression.expression import literal, param
+from gems_craft.expression.expression import LowerBoundNode, literal, param
 from gems_craft.expression.indexing_structure import IndexingStructure
 from gems_craft.model.model import model
 from gems_craft.model.parameter import float_parameter
@@ -30,6 +31,7 @@
 from gems_craft.study.parsing import HeuristicId, IntegerStrategy, IntegerStrategyId
 from gems_craft.study.system import Component
 from gems_runner.simulation import TimeBlock, build_problem
+from gems_runner.simulation.simulation_table import SimulationTableBuilder
 
 MIXED_MODEL = model(
     id="mixed_model",
@@ -179,3 +181,55 @@ def test_mixed_strategies_with_time_dependent_parameter_bound() -> None:
     assert (
         relaxed_var is not None and relaxed_var.upper.sel(component="c2").item() == 10.0
     )
+
+
+MIXED_MODEL_WITH_BOUND_OUTPUT = model(
+    id="mixed_model_with_bound_output",
+    variables=[
+        float_variable("generation", lower_bound=literal(0), upper_bound=literal(100)),
+    ],
+    extra_outputs={
+        "gen_lb": LowerBoundNode("generation"),
+    },
+)
+
+
+def test_lower_bound_extra_output_bypasses_merged_group_variable() -> None:
+    """lower_bound() extra-outputs must read the real per-component linopy
+    Variable, not the merged/detached _MergedGroupVariable copy built for
+    models split across relaxed/exact strategy groups — otherwise a
+    heuristic-style bound mutation on one component would not be visible (or
+    would leak across components).
+    """
+    system = System("test")
+    for comp_id, strategy in zip(
+        ["c1", "c2"], [IntegerStrategyId.EXACT, IntegerStrategyId.RELAXED]
+    ):
+        system.add_component(
+            Component(
+                model=MIXED_MODEL_WITH_BOUND_OUTPUT,
+                id=comp_id,
+                integer_strategy=IntegerStrategy(id=strategy),
+            )
+        )
+
+    problem = build_problem(
+        Study(system, DataBase()), TimeBlock(1, [0]), scenario_ids=[0]
+    )
+    problem.solve(solver_name="highs")
+
+    # Simulate what a heuristic does: mutate c2's bound directly via the real
+    # per-component Variable, bypassing the merged/detached copy.
+    c2_var = problem.get_component_variable(
+        "mixed_model_with_bound_output", "generation", "c2"
+    )
+    assert c2_var is not None
+    c2_var.lower.sel(component="c2")[:] = 42.0
+
+    st = SimulationTableBuilder().build(problem)
+
+    c2_lb = st.component("c2").output("gen_lb").value(time_index=0, scenario_index=0)
+    c1_lb = st.component("c1").output("gen_lb").value(time_index=0, scenario_index=0)
+
+    assert c2_lb == pytest.approx(42.0)
+    assert c1_lb == pytest.approx(0.0)
diff --git a/tests/unittests/gems_runner/simulation/test_simulation_table_extra_outputs.py b/tests/unittests/gems_runner/simulation/test_simulation_table_extra_outputs.py
index 723b2f59..2f3fabe2 100644
--- a/tests/unittests/gems_runner/simulation/test_simulation_table_extra_outputs.py
+++ b/tests/unittests/gems_runner/simulation/test_simulation_table_extra_outputs.py
@@ -225,6 +225,87 @@ def test_extra_output_min_on_variable() -> None:
     assert capped == pytest.approx(3.0)
 
 
+def test_extra_output_lower_bound_and_upper_bound() -> None:
+    """
+    lower_bound(x)/upper_bound(x) return the variable's current bounds post-solve.
+    """
+    from gems_craft.expression.expression import LowerBoundNode, UpperBoundNode, literal
+    from gems_craft.model.model import model
+    from gems_craft.model.variable import float_variable
+    from gems_craft.study import DataBase, Study, System, create_component
+    from gems_runner.simulation import TimeBlock, build_problem
+
+    SIMPLE_MODEL = model(
+        id="SIMPLE_BOUNDS",
+        variables=[float_variable("x", lower_bound=literal(2), upper_bound=literal(7))],
+        extra_outputs={
+            "lb": LowerBoundNode("x"),
+            "ub": UpperBoundNode("x"),
+        },
+    )
+
+    database = DataBase()
+    comp = create_component(model=SIMPLE_MODEL, id="comp_1")
+
+    system = System("test_bounds_extra")
+    system.add_component(comp)
+
+    problem = build_problem(
+        Study(system, database), TimeBlock(1, [0]), scenario_ids=list(range(1))
+    )
+    problem.solve(solver_name="highs")
+
+    df = SimulationTableBuilder().build(problem)
+    lb = df.component("comp_1").output("lb").value(time_index=0, scenario_index=0)
+    ub = df.component("comp_1").output("ub").value(time_index=0, scenario_index=0)
+    assert lb == pytest.approx(2.0)
+    assert ub == pytest.approx(7.0)
+
+
+def test_extra_output_bound_broadcasts_over_time() -> None:
+    """
+    upper_bound(x), when x's bound is a constant parameter, is broadcast over
+    every timestep of a time-varying variable in the resulting SimulationTable
+    (not just t=0).
+    """
+    from gems_craft.expression import param
+    from gems_craft.expression.expression import UpperBoundNode, literal
+    from gems_craft.expression.indexing_structure import IndexingStructure
+    from gems_craft.model.model import model
+    from gems_craft.model.parameter import float_parameter
+    from gems_craft.model.variable import float_variable
+    from gems_craft.study import ConstantData, DataBase, Study, System, create_component
+    from gems_runner.simulation import TimeBlock, build_problem
+
+    SIMPLE_MODEL = model(
+        id="SIMPLE_BOUND_BROADCAST",
+        parameters=[float_parameter("cap", structure=IndexingStructure(False, False))],
+        variables=[
+            float_variable("x", lower_bound=literal(0), upper_bound=param("cap"))
+        ],
+        extra_outputs={"x_ub": UpperBoundNode("x")},
+    )
+
+    database = DataBase()
+    comp = create_component(model=SIMPLE_MODEL, id="comp_1")
+    database.add_data("comp_1", "cap", ConstantData(7.0))
+
+    system = System("test_bound_broadcast")
+    system.add_component(comp)
+
+    problem = build_problem(
+        Study(system, database), TimeBlock(1, [0, 1, 2]), scenario_ids=list(range(1))
+    )
+    problem.solve(solver_name="highs")
+
+    df = SimulationTableBuilder().build(problem)
+    for t in range(3):
+        x_ub = (
+            df.component("comp_1").output("x_ub").value(time_index=t, scenario_index=0)
+        )
+        assert x_ub == pytest.approx(7.0), f"x_ub at t={t}: expected 7.0, got {x_ub}"
+
+
 def test_extra_output_comparison() -> None:
     """
     Comparison operators (>=, <=) are allowed in extra outputs and evaluated