diff --git a/flixopt/components.py b/flixopt/components.py
index 481135d1c..f9c59485b 100644
--- a/flixopt/components.py
+++ b/flixopt/components.py
@@ -4,6 +4,7 @@
 
 from __future__ import annotations
 
+import functools
 import logging
 import warnings
 from typing import TYPE_CHECKING, Literal
@@ -1102,7 +1103,7 @@ def _absolute_charge_state_bounds(self) -> tuple[xr.DataArray, xr.DataArray]:
                 relative_upper_bound * cap,
             )
 
-    @property
+    @functools.cached_property
     def _relative_charge_state_bounds(self) -> tuple[xr.DataArray, xr.DataArray]:
         """
         Get relative charge state bounds with final timestep values.
@@ -1152,7 +1153,9 @@ def _relative_charge_state_bounds(self) -> tuple[xr.DataArray, xr.DataArray]:
             # Original is scalar - broadcast to full time range (constant value)
             max_bounds = rel_max.expand_dims(time=timesteps_extra)
 
-        return min_bounds, max_bounds
+        # Ensure both bounds have matching dimensions (broadcast once here,
+        # so downstream code doesn't need to handle dimension mismatches)
+        return xr.broadcast(min_bounds, max_bounds)
 
     @property
     def _investment(self) -> InvestmentModel | None:
diff --git a/flixopt/elements.py b/flixopt/elements.py
index e2def702d..791596b28 100644
--- a/flixopt/elements.py
+++ b/flixopt/elements.py
@@ -4,6 +4,7 @@
 
 from __future__ import annotations
 
+import functools
 import logging
 from typing import TYPE_CHECKING
 
@@ -866,11 +867,13 @@ def _create_bounds_for_load_factor(self):
                 short_name='load_factor_min',
             )
 
-    @property
+    @functools.cached_property
     def relative_flow_rate_bounds(self) -> tuple[xr.DataArray, xr.DataArray]:
         if self.element.fixed_relative_profile is not None:
             return self.element.fixed_relative_profile, self.element.fixed_relative_profile
-        return self.element.relative_minimum, self.element.relative_maximum
+        # Ensure both bounds have matching dimensions (broadcast once here,
+        # so downstream code doesn't need to handle dimension mismatches)
+        return xr.broadcast(self.element.relative_minimum, self.element.relative_maximum)
 
     @property
     def absolute_flow_rate_bounds(self) -> tuple[xr.DataArray, xr.DataArray]:
diff --git a/flixopt/modeling.py b/flixopt/modeling.py
index 3adce5338..ff84c808f 100644
--- a/flixopt/modeling.py
+++ b/flixopt/modeling.py
@@ -76,6 +76,27 @@ def _scalar_safe_reduce(data: xr.DataArray | Any, dim: str, method: str = 'mean'
     return data
 
 
+def _xr_allclose(a: xr.DataArray, b: xr.DataArray, rtol: float = 1e-5, atol: float = 1e-8) -> bool:
+    """Check if two DataArrays are element-wise equal within tolerance.
+
+    Args:
+        a: First DataArray
+        b: Second DataArray
+        rtol: Relative tolerance (default matches np.allclose)
+        atol: Absolute tolerance (default matches np.allclose)
+
+    Returns:
+        True if all elements are close (including matching NaN positions)
+    """
+    # Fast path: same dims and shape - use numpy directly
+    if a.dims == b.dims and a.shape == b.shape:
+        return np.allclose(a.values, b.values, rtol=rtol, atol=atol, equal_nan=True)
+
+    # Slow path: broadcast to common shape, then use numpy
+    a_bc, b_bc = xr.broadcast(a, b)
+    return np.allclose(a_bc.values, b_bc.values, rtol=rtol, atol=atol, equal_nan=True)
+
+
 class ModelingUtilitiesAbstract:
     """Utility functions for modeling - leveraging xarray for temporal data"""
 
@@ -546,7 +567,7 @@ def bounds_with_state(
         lower_bound, upper_bound = bounds
         name = name or f'{variable.name}'
 
-        if np.allclose(lower_bound, upper_bound, atol=1e-10, equal_nan=True):
+        if _xr_allclose(lower_bound, upper_bound):
             fix_constraint = model.add_constraints(variable == state * upper_bound, name=f'{name}|fix')
             return [fix_constraint]
 
@@ -588,7 +609,7 @@ def scaled_bounds(
         rel_lower, rel_upper = relative_bounds
         name = name or f'{variable.name}'
 
-        if np.allclose(rel_lower, rel_upper, atol=1e-10, equal_nan=True):
+        if _xr_allclose(rel_lower, rel_upper):
             return [model.add_constraints(variable == scaling_variable * rel_lower, name=f'{name}|fixed')]
 
         upper_constraint = model.add_constraints(variable <= scaling_variable * rel_upper, name=f'{name}|ub')