From 09022298594fb6cf65fad9e874b86d80fec91706 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Thu, 15 Jan 2026 10:06:25 +0100
Subject: [PATCH 01/18]   Summary: Variable Registry Implementation

  Changes Made

  1. Added VariableCategory enum (structure.py:64-77)
    - STATE - For state variables like charge_state (interpolated within segments)
    - SEGMENT_TOTAL - For segment totals like effect contributions (divided by expansion divisor)
    - RATE - For rate variables like flow_rate (expanded as-is)
    - BINARY - For binary variables like status (expanded as-is)
    - OTHER - For uncategorized variables
  2. Added variable_categories registry to FlowSystemModel (structure.py:214)
    - Dictionary mapping variable names to their categories
  3. Modified add_variables() method (structure.py:388-396)
    - Added optional category parameter
    - Automatically registers variables with their category
  4. Updated variable creation calls:
    - components.py: Storage variables (charge_state as STATE, netto_discharge as RATE)
    - elements.py: Flow variables (flow_rate as RATE, status as BINARY)
    - features.py: Effect contributions (per_timestep as SEGMENT_TOTAL, temporal shares as SEGMENT_TOTAL, startup/shutdown as BINARY)
  5. Updated expand() method (transform_accessor.py:2074-2090)
    - Uses variable_categories registry to identify segment totals and state variables
    - Falls back to pattern matching for backwards compatibility with older FlowSystems

  Benefits

  - More robust categorization: Variables are categorized at creation time, not by pattern matching
  - Extensible: New variable types can easily be added with proper category
  - Backwards compatible: Old FlowSystems without categories still work via pattern matching fallback
---
 flixopt/components.py         |  9 +++++--
 flixopt/elements.py           | 16 ++++++++++--
 flixopt/features.py           | 20 ++++++++++++---
 flixopt/flow_system.py        | 16 +++++++++++-
 flixopt/structure.py          | 46 +++++++++++++++++++++++++++++++++--
 flixopt/transform_accessor.py | 28 ++++++++++++++++-----
 6 files changed, 119 insertions(+), 16 deletions(-)

diff --git a/flixopt/components.py b/flixopt/components.py
index 8c17bc6eb..f4cd775d3 100644
--- a/flixopt/components.py
+++ b/flixopt/components.py
@@ -17,7 +17,7 @@
 from .features import InvestmentModel, PiecewiseModel
 from .interface import InvestParameters, PiecewiseConversion, StatusParameters
 from .modeling import BoundingPatterns
-from .structure import FlowSystemModel, register_class_for_io
+from .structure import FlowSystemModel, VariableCategory, register_class_for_io
 
 if TYPE_CHECKING:
     import linopy
@@ -940,8 +940,13 @@ def _create_storage_variables(self):
             upper=ub,
             coords=self._model.get_coords(extra_timestep=True),
             short_name='charge_state',
+            category=VariableCategory.STATE,
+        )
+        self.add_variables(
+            coords=self._model.get_coords(),
+            short_name='netto_discharge',
+            category=VariableCategory.RATE,
         )
-        self.add_variables(coords=self._model.get_coords(), short_name='netto_discharge')
 
     def _add_netto_discharge_constraint(self):
         """Add constraint: netto_discharge = discharging - charging."""
diff --git a/flixopt/elements.py b/flixopt/elements.py
index 0cee53738..ae7d9055e 100644
--- a/flixopt/elements.py
+++ b/flixopt/elements.py
@@ -20,6 +20,7 @@
     Element,
     ElementModel,
     FlowSystemModel,
+    VariableCategory,
     register_class_for_io,
 )
 
@@ -672,6 +673,7 @@ def _do_modeling(self):
             upper=self.absolute_flow_rate_bounds[1],
             coords=self._model.get_coords(),
             short_name='flow_rate',
+            category=VariableCategory.RATE,
         )
 
         self._constraint_flow_rate()
@@ -726,7 +728,12 @@ def _do_modeling(self):
         self._create_shares()
 
     def _create_status_model(self):
-        status = self.add_variables(binary=True, short_name='status', coords=self._model.get_coords())
+        status = self.add_variables(
+            binary=True,
+            short_name='status',
+            coords=self._model.get_coords(),
+            category=VariableCategory.BINARY,
+        )
         self.add_submodels(
             StatusModel(
                 model=self._model,
@@ -1028,7 +1035,12 @@ def _do_modeling(self):
 
         # Create component status variable and StatusModel if needed
         if self.element.status_parameters:
-            status = self.add_variables(binary=True, short_name='status', coords=self._model.get_coords())
+            status = self.add_variables(
+                binary=True,
+                short_name='status',
+                coords=self._model.get_coords(),
+                category=VariableCategory.BINARY,
+            )
             if len(all_flows) == 1:
                 self.add_constraints(status == all_flows[0].submodel.status.status, short_name='status')
             else:
diff --git a/flixopt/features.py b/flixopt/features.py
index bb9864d64..64137d1a8 100644
--- a/flixopt/features.py
+++ b/flixopt/features.py
@@ -11,7 +11,7 @@
 import numpy as np
 
 from .modeling import BoundingPatterns, ModelingPrimitives, ModelingUtilities
-from .structure import FlowSystemModel, Submodel
+from .structure import FlowSystemModel, Submodel, VariableCategory
 
 if TYPE_CHECKING:
     from collections.abc import Collection
@@ -211,8 +211,18 @@ def _do_modeling(self):
 
         # 4. Switch tracking using existing pattern
         if self.parameters.use_startup_tracking:
-            self.add_variables(binary=True, short_name='startup', coords=self.get_coords())
-            self.add_variables(binary=True, short_name='shutdown', coords=self.get_coords())
+            self.add_variables(
+                binary=True,
+                short_name='startup',
+                coords=self.get_coords(),
+                category=VariableCategory.BINARY,
+            )
+            self.add_variables(
+                binary=True,
+                short_name='shutdown',
+                coords=self.get_coords(),
+                category=VariableCategory.BINARY,
+            )
 
             # Determine previous_state: None means relaxed (no constraint at t=0)
             previous_state = self._previous_status.isel(time=-1) if self._previous_status is not None else None
@@ -629,6 +639,7 @@ def _do_modeling(self):
                 upper=np.inf if (self._max_per_hour is None) else self._max_per_hour * self._model.timestep_duration,
                 coords=self._model.get_coords(self._dims),
                 short_name='per_timestep',
+                category=VariableCategory.SEGMENT_TOTAL,
             )
 
             self._eq_total_per_timestep = self.add_constraints(self.total_per_timestep == 0, short_name='per_timestep')
@@ -668,10 +679,13 @@ def add_share(
         if name in self.shares:
             self.share_constraints[name].lhs -= expression
         else:
+            # Temporal shares (with 'time' dim) are segment totals that need division
+            category = VariableCategory.SEGMENT_TOTAL if 'time' in dims else None
             self.shares[name] = self.add_variables(
                 coords=self._model.get_coords(dims),
                 name=f'{name}->{self.label_full}',
                 short_name=name,
+                category=category,
             )
 
             self.share_constraints[name] = self.add_constraints(
diff --git a/flixopt/flow_system.py b/flixopt/flow_system.py
index 3d7fe3acd..9478e3034 100644
--- a/flixopt/flow_system.py
+++ b/flixopt/flow_system.py
@@ -29,7 +29,14 @@
 from .elements import Bus, Component, Flow
 from .optimize_accessor import OptimizeAccessor
 from .statistics_accessor import StatisticsAccessor
-from .structure import CompositeContainerMixin, Element, ElementContainer, FlowSystemModel, Interface
+from .structure import (
+    CompositeContainerMixin,
+    Element,
+    ElementContainer,
+    FlowSystemModel,
+    Interface,
+    VariableCategory,
+)
 from .topology_accessor import TopologyAccessor
 from .transform_accessor import TransformAccessor
 
@@ -249,6 +256,10 @@ def __init__(
         # Solution dataset - populated after optimization or loaded from file
         self._solution: xr.Dataset | None = None
 
+        # Variable categories for segment expansion handling
+        # Populated when model is built, used by transform.expand()
+        self._variable_categories: dict[str, VariableCategory] = {}
+
         # Aggregation info - populated by transform.cluster()
         self.clustering: Clustering | None = None
 
@@ -1605,6 +1616,9 @@ def solve(self, solver: _Solver) -> FlowSystem:
         # Store solution on FlowSystem for direct Element access
         self.solution = self.model.solution
 
+        # Copy variable categories for segment expansion handling
+        self._variable_categories = self.model.variable_categories.copy()
+
         logger.info(f'Optimization solved successfully. Objective: {self.model.objective.value:.4f}')
 
         return self
diff --git a/flixopt/structure.py b/flixopt/structure.py
index 67a02750f..3511763b8 100644
--- a/flixopt/structure.py
+++ b/flixopt/structure.py
@@ -13,6 +13,7 @@
 import warnings
 from dataclasses import dataclass
 from difflib import get_close_matches
+from enum import Enum
 from typing import (
     TYPE_CHECKING,
     Any,
@@ -40,6 +41,27 @@
 
 logger = logging.getLogger('flixopt')
 
+
+class VariableCategory(Enum):
+    """Category of a solution variable for segment expansion handling.
+
+    When expanding segmented clustered systems back to hourly resolution,
+    different variable types require different handling:
+
+    - STATE: Interpolated between segment boundaries (e.g., charge_state)
+    - SEGMENT_TOTAL: Divided by segment duration (e.g., effect contributions)
+    - RATE: Expanded as-is, already averaged by tsam (e.g., flow_rate)
+    - BINARY: Constant within segment, cannot interpolate (e.g., status)
+    - OTHER: Default, no special handling
+    """
+
+    STATE = 'state'
+    SEGMENT_TOTAL = 'segment_total'
+    RATE = 'rate'
+    BINARY = 'binary'
+    OTHER = 'other'
+
+
 CLASS_REGISTRY = {}
 
 
@@ -97,6 +119,7 @@ def __init__(self, flow_system: FlowSystem):
         self.flow_system = flow_system
         self.effects: EffectCollectionModel | None = None
         self.submodels: Submodels = Submodels({})
+        self.variable_categories: dict[str, VariableCategory] = {}
 
     def do_modeling(self):
         # Create all element models
@@ -1603,8 +1626,22 @@ def __init__(self, model: FlowSystemModel, label_of_element: str, label_of_model
         logger.debug(f'Creating {self.__class__.__name__}  "{self.label_full}"')
         self._do_modeling()
 
-    def add_variables(self, short_name: str = None, **kwargs) -> linopy.Variable:
-        """Create and register a variable in one step"""
+    def add_variables(
+        self,
+        short_name: str = None,
+        category: VariableCategory = None,
+        **kwargs,
+    ) -> linopy.Variable:
+        """Create and register a variable in one step.
+
+        Args:
+            short_name: Short name for the variable (used as suffix in full name).
+            category: Category for segment expansion handling. See VariableCategory.
+            **kwargs: Additional arguments passed to linopy.Model.add_variables().
+
+        Returns:
+            The created linopy Variable.
+        """
         if kwargs.get('name') is None:
             if short_name is None:
                 raise ValueError('Short name must be provided when no name is given')
@@ -1612,6 +1649,11 @@ def add_variables(self, short_name: str = None, **kwargs) -> linopy.Variable:
 
         variable = self._model.add_variables(**kwargs)
         self.register_variable(variable, short_name)
+
+        # Register category in FlowSystemModel for segment expansion handling
+        if category is not None:
+            self._model.variable_categories[variable.name] = category
+
         return variable
 
     def add_constraints(self, expression, short_name: str = None, **kwargs) -> linopy.Constraint:
diff --git a/flixopt/transform_accessor.py b/flixopt/transform_accessor.py
index 57fae3fec..148e3111e 100644
--- a/flixopt/transform_accessor.py
+++ b/flixopt/transform_accessor.py
@@ -16,6 +16,8 @@
 import pandas as pd
 import xarray as xr
 
+from .structure import VariableCategory
+
 if TYPE_CHECKING:
     from tsam.config import ClusterConfig, ExtremeConfig, SegmentConfig
 
@@ -2069,9 +2071,23 @@ def expand(self) -> FlowSystem:
         # For segmented systems: build expansion divisor and identify segment total variables
         expansion_divisor = None
         segment_total_vars: set[str] = set()
+        variable_categories = getattr(self._fs, '_variable_categories', {})
         if clustering.is_segmented:
             expansion_divisor = clustering.build_expansion_divisor(original_time=original_timesteps)
-            segment_total_vars = self._build_segment_total_varnames()
+            # Build segment total vars using registry first, fall back to pattern matching
+            segment_total_vars = {
+                name for name, cat in variable_categories.items() if cat == VariableCategory.SEGMENT_TOTAL
+            }
+            # Fall back to pattern matching for backwards compatibility (old FlowSystems without categories)
+            if not segment_total_vars:
+                segment_total_vars = self._build_segment_total_varnames()
+
+        def _is_state_variable(var_name: str) -> bool:
+            """Check if a variable is a state variable (should be interpolated)."""
+            if var_name in variable_categories:
+                return variable_categories[var_name] == VariableCategory.STATE
+            # Fall back to pattern matching for backwards compatibility
+            return var_name.endswith('|charge_state')
 
         def expand_da(da: xr.DataArray, var_name: str = '', is_solution: bool = False) -> xr.DataArray:
             """Expand a DataArray from clustered to original timesteps.
@@ -2085,9 +2101,9 @@ def expand_da(da: xr.DataArray, var_name: str = '', is_solution: bool = False) -
             if 'time' not in da.dims:
                 return da.copy()
 
-            # For charge_state in segmented systems: interpolate within segments
-            # to show the actual charge trajectory as storage charges/discharges
-            if var_name.endswith('|charge_state') and 'cluster' in da.dims and clustering.is_segmented:
+            # For state variables (like charge_state) in segmented systems: interpolate within segments
+            # to show the actual state trajectory as the storage charges/discharges
+            if _is_state_variable(var_name) and 'cluster' in da.dims and clustering.is_segmented:
                 expanded = self._interpolate_charge_state_segmented(da, clustering, original_timesteps)
                 # Append the extra timestep value (final charge state)
                 cluster_assignments = clustering.cluster_assignments
@@ -2109,8 +2125,8 @@ def expand_da(da: xr.DataArray, var_name: str = '', is_solution: bool = False) -
             if is_solution and expansion_divisor is not None and var_name in segment_total_vars:
                 expanded = expanded / expansion_divisor
 
-            # For charge_state with cluster dim (non-segmented), append the extra timestep value
-            if var_name.endswith('|charge_state') and 'cluster' in da.dims:
+            # For state variables (like charge_state) with cluster dim (non-segmented), append the extra timestep value
+            if _is_state_variable(var_name) and 'cluster' in da.dims:
                 cluster_assignments = clustering.cluster_assignments
                 if cluster_assignments.ndim == 1:
                     last_cluster = int(cluster_assignments[last_original_cluster_idx])

From 6fa5a24aabcebf398bb76dcb8cb37edbd5a5aaa0 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Thu, 15 Jan 2026 10:34:22 +0100
Subject: [PATCH 02/18]   Summary: Fine-Grained Variable Categories
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

  New Categories (structure.py:45-103)

  class VariableCategory(Enum):
      # State variables
      CHARGE_STATE, SOC_BOUNDARY

      # Rate/Power variables
      FLOW_RATE, NETTO_DISCHARGE, VIRTUAL_FLOW

      # Binary state
      STATUS, INACTIVE

      # Binary events
      STARTUP, SHUTDOWN

      # Effect variables
      PER_TIMESTEP, SHARE, TOTAL, TOTAL_OVER_PERIODS

      # Investment
      SIZE, INVESTED

      # Counting/Duration
      STARTUP_COUNT, DURATION

      # Piecewise linearization
      INSIDE_PIECE, LAMBDA0, LAMBDA1, ZERO_POINT

      # Other
      OTHER

  Logical Groupings for Expansion

  EXPAND_INTERPOLATE = {CHARGE_STATE}      # Interpolate between boundaries
  EXPAND_DIVIDE = {PER_TIMESTEP, SHARE}    # Divide by expansion factor
  # Default: repeat within segment

  Files Modified
  ┌───────────────────────┬─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┐
  │         File          │                                                          Variables Updated                                                          │
  ├───────────────────────┼─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┤
  │ components.py         │ charge_state, netto_discharge, SOC_boundary                                                                                         │
  ├───────────────────────┼─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┤
  │ elements.py           │ flow_rate, status, virtual_supply, virtual_demand                                                                                   │
  ├───────────────────────┼─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┤
  │ features.py           │ size, invested, inactive, startup, shutdown, startup_count, inside_piece, lambda0, lambda1, zero_point, total, per_timestep, shares │
  ├───────────────────────┼─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┤
  │ effects.py            │ total, total_over_periods                                                                                                           │
  ├───────────────────────┼─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┤
  │ modeling.py           │ duration                                                                                                                            │
  ├───────────────────────┼─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┤
  │ transform_accessor.py │ Updated to use EXPAND_INTERPOLATE and EXPAND_DIVIDE groupings                                                                       │
  └───────────────────────┴─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┘
  Test Results

  - All 83 cluster/expand tests pass
  - Variable categories correctly populated and grouped
---
 flixopt/components.py         |  5 ++-
 flixopt/effects.py            | 12 +++++-
 flixopt/elements.py           | 16 +++++---
 flixopt/features.py           | 23 +++++++++---
 flixopt/modeling.py           |  3 +-
 flixopt/structure.py          | 69 ++++++++++++++++++++++++++++-------
 flixopt/transform_accessor.py |  8 ++--
 7 files changed, 103 insertions(+), 33 deletions(-)

diff --git a/flixopt/components.py b/flixopt/components.py
index f4cd775d3..56a0f5632 100644
--- a/flixopt/components.py
+++ b/flixopt/components.py
@@ -940,12 +940,12 @@ def _create_storage_variables(self):
             upper=ub,
             coords=self._model.get_coords(extra_timestep=True),
             short_name='charge_state',
-            category=VariableCategory.STATE,
+            category=VariableCategory.CHARGE_STATE,
         )
         self.add_variables(
             coords=self._model.get_coords(),
             short_name='netto_discharge',
-            category=VariableCategory.RATE,
+            category=VariableCategory.NETTO_DISCHARGE,
         )
 
     def _add_netto_discharge_constraint(self):
@@ -1347,6 +1347,7 @@ def _add_intercluster_linking(self) -> None:
             coords=boundary_coords,
             dims=boundary_dims,
             short_name='SOC_boundary',
+            category=VariableCategory.SOC_BOUNDARY,
         )
 
         # 3. Link SOC_boundary to investment size
diff --git a/flixopt/effects.py b/flixopt/effects.py
index 3a2322988..b32a4edd8 100644
--- a/flixopt/effects.py
+++ b/flixopt/effects.py
@@ -17,7 +17,15 @@
 
 from .core import PlausibilityError
 from .features import ShareAllocationModel
-from .structure import Element, ElementContainer, ElementModel, FlowSystemModel, Submodel, register_class_for_io
+from .structure import (
+    Element,
+    ElementContainer,
+    ElementModel,
+    FlowSystemModel,
+    Submodel,
+    VariableCategory,
+    register_class_for_io,
+)
 
 if TYPE_CHECKING:
     from collections.abc import Iterator
@@ -377,6 +385,7 @@ def _do_modeling(self):
             upper=self.element.maximum_total if self.element.maximum_total is not None else np.inf,
             coords=self._model.get_coords(['period', 'scenario']),
             name=self.label_full,
+            category=VariableCategory.TOTAL,
         )
 
         self.add_constraints(
@@ -394,6 +403,7 @@ def _do_modeling(self):
                 upper=self.element.maximum_over_periods if self.element.maximum_over_periods is not None else np.inf,
                 coords=self._model.get_coords(['scenario']),
                 short_name='total_over_periods',
+                category=VariableCategory.TOTAL_OVER_PERIODS,
             )
 
             self.add_constraints(self.total_over_periods == weighted_total, short_name='total_over_periods')
diff --git a/flixopt/elements.py b/flixopt/elements.py
index ae7d9055e..0fa620bed 100644
--- a/flixopt/elements.py
+++ b/flixopt/elements.py
@@ -673,7 +673,7 @@ def _do_modeling(self):
             upper=self.absolute_flow_rate_bounds[1],
             coords=self._model.get_coords(),
             short_name='flow_rate',
-            category=VariableCategory.RATE,
+            category=VariableCategory.FLOW_RATE,
         )
 
         self._constraint_flow_rate()
@@ -732,7 +732,7 @@ def _create_status_model(self):
             binary=True,
             short_name='status',
             coords=self._model.get_coords(),
-            category=VariableCategory.BINARY,
+            category=VariableCategory.STATUS,
         )
         self.add_submodels(
             StatusModel(
@@ -964,11 +964,17 @@ def _do_modeling(self):
             imbalance_penalty = self.element.imbalance_penalty_per_flow_hour * self._model.timestep_duration
 
             self.virtual_supply = self.add_variables(
-                lower=0, coords=self._model.get_coords(), short_name='virtual_supply'
+                lower=0,
+                coords=self._model.get_coords(),
+                short_name='virtual_supply',
+                category=VariableCategory.VIRTUAL_FLOW,
             )
 
             self.virtual_demand = self.add_variables(
-                lower=0, coords=self._model.get_coords(), short_name='virtual_demand'
+                lower=0,
+                coords=self._model.get_coords(),
+                short_name='virtual_demand',
+                category=VariableCategory.VIRTUAL_FLOW,
             )
 
             # Σ(inflows) + virtual_supply = Σ(outflows) + virtual_demand
@@ -1039,7 +1045,7 @@ def _do_modeling(self):
                 binary=True,
                 short_name='status',
                 coords=self._model.get_coords(),
-                category=VariableCategory.BINARY,
+                category=VariableCategory.STATUS,
             )
             if len(all_flows) == 1:
                 self.add_constraints(status == all_flows[0].submodel.status.status, short_name='status')
diff --git a/flixopt/features.py b/flixopt/features.py
index 64137d1a8..af1ff1471 100644
--- a/flixopt/features.py
+++ b/flixopt/features.py
@@ -69,6 +69,7 @@ def _create_variables_and_constraints(self):
             lower=size_min if self.parameters.mandatory else 0,
             upper=size_max,
             coords=self._model.get_coords(['period', 'scenario']),
+            category=VariableCategory.SIZE,
         )
 
         if not self.parameters.mandatory:
@@ -76,6 +77,7 @@ def _create_variables_and_constraints(self):
                 binary=True,
                 coords=self._model.get_coords(['period', 'scenario']),
                 short_name='invested',
+                category=VariableCategory.INVESTED,
             )
             BoundingPatterns.bounds_with_state(
                 self,
@@ -193,7 +195,12 @@ def _do_modeling(self):
         # Create a separate binary 'inactive' variable when needed for downtime tracking or explicit use
         # When not needed, the expression (1 - self.status) can be used instead
         if self.parameters.use_downtime_tracking:
-            inactive = self.add_variables(binary=True, short_name='inactive', coords=self._model.get_coords())
+            inactive = self.add_variables(
+                binary=True,
+                short_name='inactive',
+                coords=self._model.get_coords(),
+                category=VariableCategory.INACTIVE,
+            )
             self.add_constraints(self.status + inactive == 1, short_name='complementary')
 
         # 3. Total duration tracking
@@ -215,13 +222,13 @@ def _do_modeling(self):
                 binary=True,
                 short_name='startup',
                 coords=self.get_coords(),
-                category=VariableCategory.BINARY,
+                category=VariableCategory.STARTUP,
             )
             self.add_variables(
                 binary=True,
                 short_name='shutdown',
                 coords=self.get_coords(),
-                category=VariableCategory.BINARY,
+                category=VariableCategory.SHUTDOWN,
             )
 
             # Determine previous_state: None means relaxed (no constraint at t=0)
@@ -243,6 +250,7 @@ def _do_modeling(self):
                     upper=self.parameters.startup_limit,
                     coords=self._model.get_coords(('period', 'scenario')),
                     short_name='startup_count',
+                    category=VariableCategory.STARTUP_COUNT,
                 )
                 # Sum over all temporal dimensions (time, and cluster if present)
                 startup_temporal_dims = [d for d in self.startup.dims if d not in ('period', 'scenario')]
@@ -397,12 +405,14 @@ def _do_modeling(self):
             binary=True,
             short_name='inside_piece',
             coords=self._model.get_coords(dims=self.dims),
+            category=VariableCategory.INSIDE_PIECE,
         )
         self.lambda0 = self.add_variables(
             lower=0,
             upper=1,
             short_name='lambda0',
             coords=self._model.get_coords(dims=self.dims),
+            category=VariableCategory.LAMBDA0,
         )
 
         self.lambda1 = self.add_variables(
@@ -410,6 +420,7 @@ def _do_modeling(self):
             upper=1,
             short_name='lambda1',
             coords=self._model.get_coords(dims=self.dims),
+            category=VariableCategory.LAMBDA1,
         )
 
         # Create constraints
@@ -505,6 +516,7 @@ def _do_modeling(self):
                     coords=self._model.get_coords(self.dims),
                     binary=True,
                     short_name='zero_point',
+                    category=VariableCategory.ZERO_POINT,
                 )
                 rhs = self.zero_point
             else:
@@ -629,6 +641,7 @@ def _do_modeling(self):
             coords=self._model.get_coords([dim for dim in self._dims if dim != 'time']),
             name=self.label_full,
             short_name='total',
+            category=VariableCategory.TOTAL,
         )
         # eq: sum = sum(share_i) # skalar
         self._eq_total = self.add_constraints(self.total == 0, name=self.label_full)
@@ -639,7 +652,7 @@ def _do_modeling(self):
                 upper=np.inf if (self._max_per_hour is None) else self._max_per_hour * self._model.timestep_duration,
                 coords=self._model.get_coords(self._dims),
                 short_name='per_timestep',
-                category=VariableCategory.SEGMENT_TOTAL,
+                category=VariableCategory.PER_TIMESTEP,
             )
 
             self._eq_total_per_timestep = self.add_constraints(self.total_per_timestep == 0, short_name='per_timestep')
@@ -680,7 +693,7 @@ def add_share(
             self.share_constraints[name].lhs -= expression
         else:
             # Temporal shares (with 'time' dim) are segment totals that need division
-            category = VariableCategory.SEGMENT_TOTAL if 'time' in dims else None
+            category = VariableCategory.SHARE if 'time' in dims else None
             self.shares[name] = self.add_variables(
                 coords=self._model.get_coords(dims),
                 name=f'{name}->{self.label_full}',
diff --git a/flixopt/modeling.py b/flixopt/modeling.py
index f22ffbc04..31511c90e 100644
--- a/flixopt/modeling.py
+++ b/flixopt/modeling.py
@@ -5,7 +5,7 @@
 import xarray as xr
 
 from .config import CONFIG
-from .structure import Submodel
+from .structure import Submodel, VariableCategory
 
 logger = logging.getLogger('flixopt')
 
@@ -303,6 +303,7 @@ def consecutive_duration_tracking(
             coords=state.coords,
             name=name,
             short_name=short_name,
+            category=VariableCategory.DURATION,
         )
 
         constraints = {}
diff --git a/flixopt/structure.py b/flixopt/structure.py
index 3511763b8..8ac18e76a 100644
--- a/flixopt/structure.py
+++ b/flixopt/structure.py
@@ -43,23 +43,64 @@
 
 
 class VariableCategory(Enum):
-    """Category of a solution variable for segment expansion handling.
+    """Fine-grained variable categories - names mirror variable names.
 
-    When expanding segmented clustered systems back to hourly resolution,
-    different variable types require different handling:
-
-    - STATE: Interpolated between segment boundaries (e.g., charge_state)
-    - SEGMENT_TOTAL: Divided by segment duration (e.g., effect contributions)
-    - RATE: Expanded as-is, already averaged by tsam (e.g., flow_rate)
-    - BINARY: Constant within segment, cannot interpolate (e.g., status)
-    - OTHER: Default, no special handling
+    Each variable type has its own category for precise handling during
+    segment expansion and statistics calculation.
     """
 
-    STATE = 'state'
-    SEGMENT_TOTAL = 'segment_total'
-    RATE = 'rate'
-    BINARY = 'binary'
-    OTHER = 'other'
+    # === State variables ===
+    CHARGE_STATE = 'charge_state'  # Storage SOC (interpolate between boundaries)
+    SOC_BOUNDARY = 'soc_boundary'  # Intercluster SOC boundaries
+
+    # === Rate/Power variables ===
+    FLOW_RATE = 'flow_rate'  # Flow rate (kW)
+    NETTO_DISCHARGE = 'netto_discharge'  # Storage net discharge
+    VIRTUAL_FLOW = 'virtual_flow'  # Bus penalty slack variables
+
+    # === Binary state ===
+    STATUS = 'status'  # On/off status (persists through segment)
+    INACTIVE = 'inactive'  # Complementary inactive status
+
+    # === Binary events ===
+    STARTUP = 'startup'  # Startup event
+    SHUTDOWN = 'shutdown'  # Shutdown event
+
+    # === Effect variables ===
+    PER_TIMESTEP = 'per_timestep'  # Effect per timestep
+    SHARE = 'share'  # All temporal contributions (flow, active, startup)
+    TOTAL = 'total'  # Effect total (per period/scenario)
+    TOTAL_OVER_PERIODS = 'total_over_periods'  # Effect total over all periods
+
+    # === Investment ===
+    SIZE = 'size'  # Investment size
+    INVESTED = 'invested'  # Invested yes/no binary
+
+    # === Counting/Duration ===
+    STARTUP_COUNT = 'startup_count'  # Count of startups
+    DURATION = 'duration'  # Duration tracking (uptime/downtime)
+
+    # === Piecewise linearization ===
+    INSIDE_PIECE = 'inside_piece'  # Binary segment selection
+    LAMBDA0 = 'lambda0'  # Interpolation weight
+    LAMBDA1 = 'lambda1'  # Interpolation weight
+    ZERO_POINT = 'zero_point'  # Zero point handling
+
+    # === Other ===
+    OTHER = 'other'  # Uncategorized
+
+
+# === Logical Groupings for Segment Expansion ===
+# Default behavior (not listed): repeat value within segment
+
+EXPAND_INTERPOLATE: set[VariableCategory] = {VariableCategory.CHARGE_STATE}
+"""State variables that should be interpolated between segment boundaries."""
+
+EXPAND_DIVIDE: set[VariableCategory] = {VariableCategory.PER_TIMESTEP, VariableCategory.SHARE}
+"""Segment totals that should be divided by expansion factor to preserve sums."""
+
+EXPAND_FIRST_TIMESTEP: set[VariableCategory] = {VariableCategory.STARTUP, VariableCategory.SHUTDOWN}
+"""Binary events that should appear only at the first timestep of the segment."""
 
 
 CLASS_REGISTRY = {}
diff --git a/flixopt/transform_accessor.py b/flixopt/transform_accessor.py
index 148e3111e..909004517 100644
--- a/flixopt/transform_accessor.py
+++ b/flixopt/transform_accessor.py
@@ -16,7 +16,7 @@
 import pandas as pd
 import xarray as xr
 
-from .structure import VariableCategory
+from .structure import EXPAND_DIVIDE, EXPAND_INTERPOLATE
 
 if TYPE_CHECKING:
     from tsam.config import ClusterConfig, ExtremeConfig, SegmentConfig
@@ -2075,9 +2075,7 @@ def expand(self) -> FlowSystem:
         if clustering.is_segmented:
             expansion_divisor = clustering.build_expansion_divisor(original_time=original_timesteps)
             # Build segment total vars using registry first, fall back to pattern matching
-            segment_total_vars = {
-                name for name, cat in variable_categories.items() if cat == VariableCategory.SEGMENT_TOTAL
-            }
+            segment_total_vars = {name for name, cat in variable_categories.items() if cat in EXPAND_DIVIDE}
             # Fall back to pattern matching for backwards compatibility (old FlowSystems without categories)
             if not segment_total_vars:
                 segment_total_vars = self._build_segment_total_varnames()
@@ -2085,7 +2083,7 @@ def expand(self) -> FlowSystem:
         def _is_state_variable(var_name: str) -> bool:
             """Check if a variable is a state variable (should be interpolated)."""
             if var_name in variable_categories:
-                return variable_categories[var_name] == VariableCategory.STATE
+                return variable_categories[var_name] in EXPAND_INTERPOLATE
             # Fall back to pattern matching for backwards compatibility
             return var_name.endswith('|charge_state')
 

From 47969a70c23e7620b0d016d7221a86781e0f3805 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Thu, 15 Jan 2026 10:36:48 +0100
Subject: [PATCH 03/18] Add IO for variable categories

---
 flixopt/flow_system.py | 14 ++++++++++++++
 1 file changed, 14 insertions(+)

diff --git a/flixopt/flow_system.py b/flixopt/flow_system.py
index 9478e3034..155d8379b 100644
--- a/flixopt/flow_system.py
+++ b/flixopt/flow_system.py
@@ -747,6 +747,12 @@ def to_dataset(self, include_solution: bool = True, include_original_data: bool
                 ds[f'clustering|{name}'] = arr
             ds.attrs['clustering'] = json.dumps(clustering_ref)
 
+        # Serialize variable categories for segment expansion handling
+        if self._variable_categories:
+            # Convert enum values to strings for JSON serialization
+            categories_dict = {name: cat.value for name, cat in self._variable_categories.items()}
+            ds.attrs['variable_categories'] = json.dumps(categories_dict)
+
         # Add version info
         ds.attrs['flixopt_version'] = __version__
 
@@ -909,6 +915,14 @@ def from_dataset(cls, ds: xr.Dataset) -> FlowSystem:
             if hasattr(clustering, 'representative_weights'):
                 flow_system.cluster_weight = clustering.representative_weights
 
+        # Restore variable categories if present
+        if 'variable_categories' in reference_structure:
+            categories_dict = json.loads(reference_structure['variable_categories'])
+            # Convert string values back to VariableCategory enum
+            flow_system._variable_categories = {
+                name: VariableCategory(value) for name, value in categories_dict.items()
+            }
+
         # Reconnect network to populate bus inputs/outputs (not stored in NetCDF).
         flow_system.connect_and_transform()
 

From f7e066987cd36b31cd66a1e8eed8bd3368aab641 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Thu, 15 Jan 2026 11:07:13 +0100
Subject: [PATCH 04/18]   The refactoring is complete. Here's what was
 accomplished:
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

  Changes Made

  1. Added combine_slices() utility to flixopt/clustering/base.py (lines 52-107)
    - Simple function that stacks dict of {(dim_values): np.ndarray} into a DataArray
    - Much cleaner than the previous reverse-concat pattern
  2. Refactored 3 methods to use the new utility:
    - Clustering.expand_data() - reduced from ~25 to ~12 lines
    - Clustering.build_expansion_divisor() - reduced from ~35 to ~20 lines
    - TransformAccessor._interpolate_charge_state_segmented() - reduced from ~43 to ~27 lines
  3. Added 4 unit tests for combine_slices() in tests/test_cluster_reduce_expand.py

  Results
  ┌───────────────────────────────────┬──────────┬────────────────────────┐
  │              Metric               │  Before  │         After          │
  ├───────────────────────────────────┼──────────┼────────────────────────┤
  │ Complex reverse-concat blocks     │ 3        │ 0                      │
  ├───────────────────────────────────┼──────────┼────────────────────────┤
  │ Lines of dimension iteration code │ ~100     │ ~60                    │
  ├───────────────────────────────────┼──────────┼────────────────────────┤
  │ Test coverage                     │ 83 tests │ 87 tests (all passing) │
  └───────────────────────────────────┴──────────┴────────────────────────┘
  The Pattern Change

  Before (complex reverse-concat):
  result_arrays = slices
  for dim in reversed(extra_dims):
      grouped = {}
      for key, arr in result_arrays.items():
          rest_key = key[:-1] if len(key) > 1 else ()
          grouped.setdefault(rest_key, []).append(arr)
      result_arrays = {k: xr.concat(v, dim=...) for k, v in grouped.items()}
  result = list(result_arrays.values())[0].transpose('time', ...)

  After (simple combine):
  return combine_slices(slices, extra_dims, dim_coords, 'time', output_coord, attrs)
---
 flixopt/clustering/base.py          | 107 ++++++++++++++++++----------
 flixopt/transform_accessor.py       |  30 ++------
 tests/test_cluster_reduce_expand.py |  83 +++++++++++++++++++++
 3 files changed, 158 insertions(+), 62 deletions(-)

diff --git a/flixopt/clustering/base.py b/flixopt/clustering/base.py
index 322bc71e4..fec150bcb 100644
--- a/flixopt/clustering/base.py
+++ b/flixopt/clustering/base.py
@@ -49,6 +49,64 @@ def _select_dims(da: xr.DataArray, period: Any = None, scenario: Any = None) ->
     return da
 
 
+def combine_slices(
+    slices: dict[tuple, np.ndarray],
+    extra_dims: list[str],
+    dim_coords: dict[str, list],
+    output_dim: str,
+    output_coord: Any,
+    attrs: dict | None = None,
+) -> xr.DataArray:
+    """Combine {(dim_values): 1D_array} dict into a DataArray.
+
+    This utility simplifies the common pattern of iterating over extra dimensions
+    (like period, scenario), processing each slice, and combining results.
+
+    Args:
+        slices: Dict mapping dimension value tuples to 1D numpy arrays.
+            Keys are tuples like ('period1', 'scenario1') matching extra_dims order.
+        extra_dims: Dimension names in order (e.g., ['period', 'scenario']).
+        dim_coords: Dict mapping dimension names to coordinate values.
+        output_dim: Name of the output dimension (typically 'time').
+        output_coord: Coordinate values for output dimension.
+        attrs: Optional DataArray attributes.
+
+    Returns:
+        DataArray with dims [output_dim, *extra_dims].
+
+    Example:
+        >>> slices = {
+        ...     ('P1', 'base'): np.array([1, 2, 3]),
+        ...     ('P1', 'high'): np.array([4, 5, 6]),
+        ...     ('P2', 'base'): np.array([7, 8, 9]),
+        ...     ('P2', 'high'): np.array([10, 11, 12]),
+        ... }
+        >>> result = combine_slices(
+        ...     slices,
+        ...     extra_dims=['period', 'scenario'],
+        ...     dim_coords={'period': ['P1', 'P2'], 'scenario': ['base', 'high']},
+        ...     output_dim='time',
+        ...     output_coord=[0, 1, 2],
+        ... )
+        >>> result.dims
+        ('time', 'period', 'scenario')
+    """
+    n_output = len(next(iter(slices.values())))
+    shape = [n_output] + [len(dim_coords[d]) for d in extra_dims]
+    data = np.empty(shape)
+
+    for combo in np.ndindex(*shape[1:]):
+        key = tuple(dim_coords[d][i] for d, i in zip(extra_dims, combo, strict=True))
+        data[(slice(None),) + combo] = slices[key]
+
+    return xr.DataArray(
+        data,
+        dims=[output_dim] + extra_dims,
+        coords={output_dim: output_coord, **dim_coords},
+        attrs=attrs or {},
+    )
+
+
 def _cluster_occurrences(cr: TsamClusteringResult) -> np.ndarray:
     """Compute cluster occurrences from ClusteringResult."""
     counts = Counter(cr.cluster_assignments)
@@ -890,32 +948,19 @@ def _expand_slice(mapping: np.ndarray, data: xr.DataArray) -> np.ndarray:
                 attrs=aggregated.attrs,
             )
 
-        # Multi-dimensional: expand each slice and recombine
+        # Multi-dimensional: expand each slice and combine
         dim_coords = {d: list(timestep_mapping.coords[d].values) for d in extra_dims}
-        expanded_slices = {}
+        slices = {}
         for combo in np.ndindex(*[len(v) for v in dim_coords.values()]):
             selector = {d: dim_coords[d][i] for d, i in zip(extra_dims, combo, strict=True)}
+            key = tuple(selector.values())
             mapping = _select_dims(timestep_mapping, **selector).values
             data_slice = (
                 _select_dims(aggregated, **selector) if any(d in aggregated.dims for d in selector) else aggregated
             )
-            expanded_slices[tuple(selector.values())] = xr.DataArray(
-                _expand_slice(mapping, data_slice),
-                coords={'time': original_time},
-                dims=['time'],
-            )
+            slices[key] = _expand_slice(mapping, data_slice)
 
-        # Concatenate along extra dimensions
-        result_arrays = expanded_slices
-        for dim in reversed(extra_dims):
-            dim_vals = dim_coords[dim]
-            grouped = {}
-            for key, arr in result_arrays.items():
-                rest_key = key[:-1] if len(key) > 1 else ()
-                grouped.setdefault(rest_key, []).append(arr)
-            result_arrays = {k: xr.concat(v, dim=pd.Index(dim_vals, name=dim)) for k, v in grouped.items()}
-        result = list(result_arrays.values())[0]
-        return result.transpose('time', ...).assign_attrs(aggregated.attrs)
+        return combine_slices(slices, extra_dims, dim_coords, 'time', original_time, aggregated.attrs)
 
     def build_expansion_divisor(
         self,
@@ -990,12 +1035,13 @@ def _build_divisor_for_result(
                 name='expansion_divisor',
             )
 
-        # Multi-dimensional: build divisor for each period/scenario
+        # Multi-dimensional: build divisor for each period/scenario and combine
         dim_coords = self.results.coords
-        divisor_slices = {}
+        slices = {}
 
         for combo in np.ndindex(*[len(v) for v in dim_coords.values()]):
             selector = {d: dim_coords[d][i] for d, i in zip(extra_dims, combo, strict=True)}
+            key = tuple(selector.values())
 
             # Get cluster assignments for this period/scenario
             if any(d in cluster_assignments.dims for d in selector):
@@ -1005,26 +1051,9 @@ def _build_divisor_for_result(
 
             # Get the clustering result for this period/scenario
             result = self.results.sel(**selector)
-            divisor_values = _build_divisor_for_result(assignments, result)
-
-            divisor_slices[tuple(selector.values())] = xr.DataArray(
-                divisor_values,
-                dims=['time'],
-                coords={'time': original_time},
-            )
+            slices[key] = _build_divisor_for_result(assignments, result)
 
-        # Concatenate along extra dimensions
-        result_arrays = divisor_slices
-        for dim in reversed(extra_dims):
-            dim_vals = dim_coords[dim]
-            grouped = {}
-            for key, arr in result_arrays.items():
-                rest_key = key[:-1] if len(key) > 1 else ()
-                grouped.setdefault(rest_key, []).append(arr)
-            result_arrays = {k: xr.concat(v, dim=pd.Index(dim_vals, name=dim)) for k, v in grouped.items()}
-
-        result = list(result_arrays.values())[0]
-        return result.transpose('time', ...).assign_attrs({'name': 'expansion_divisor'})
+        return combine_slices(slices, extra_dims, dim_coords, 'time', original_time, {'name': 'expansion_divisor'})
 
     def get_result(
         self,
diff --git a/flixopt/transform_accessor.py b/flixopt/transform_accessor.py
index 909004517..c13a00402 100644
--- a/flixopt/transform_accessor.py
+++ b/flixopt/transform_accessor.py
@@ -1904,17 +1904,18 @@ def _interpolate_slice(
                 attrs=da.attrs,
             )
 
-        # Multi-dimensional case
+        # Multi-dimensional case: interpolate each slice and combine
+        from .clustering.base import _select_dims, combine_slices
+
         dim_coords = clustering.results.coords
-        interpolated_slices = {}
+        slices = {}
 
         for combo in np.ndindex(*[len(v) for v in dim_coords.values()]):
             selector = {d: dim_coords[d][i] for d, i in zip(extra_dims, combo, strict=True)}
+            key = tuple(selector.values())
 
             # Get cluster assignments for this period/scenario
             if any(d in cluster_assignments.dims for d in selector):
-                from .clustering.base import _select_dims
-
                 assignments = _select_dims(cluster_assignments, **selector).values
             else:
                 assignments = cluster_assignments.values
@@ -1927,26 +1928,9 @@ def _interpolate_slice(
 
             # Get clustering result for this period/scenario
             result = clustering.results.sel(**selector)
-            interpolated = _interpolate_slice(da_slice.values, assignments, result)
-
-            interpolated_slices[tuple(selector.values())] = xr.DataArray(
-                interpolated,
-                dims=['time'],
-                coords={'time': original_timesteps},
-            )
+            slices[key] = _interpolate_slice(da_slice.values, assignments, result)
 
-        # Concatenate along extra dimensions
-        result_arrays = interpolated_slices
-        for dim in reversed(extra_dims):
-            dim_vals = dim_coords[dim]
-            grouped = {}
-            for key, arr in result_arrays.items():
-                rest_key = key[:-1] if len(key) > 1 else ()
-                grouped.setdefault(rest_key, []).append(arr)
-            result_arrays = {k: xr.concat(v, dim=pd.Index(dim_vals, name=dim)) for k, v in grouped.items()}
-
-        result = list(result_arrays.values())[0]
-        return result.transpose('time', ...).assign_attrs(da.attrs)
+        return combine_slices(slices, extra_dims, dim_coords, 'time', original_timesteps, da.attrs)
 
     def expand(self) -> FlowSystem:
         """Expand a clustered FlowSystem back to full original timesteps.
diff --git a/tests/test_cluster_reduce_expand.py b/tests/test_cluster_reduce_expand.py
index 9c119ee2d..c187358d7 100644
--- a/tests/test_cluster_reduce_expand.py
+++ b/tests/test_cluster_reduce_expand.py
@@ -1393,3 +1393,86 @@ def test_segmented_expand_after_load(self, solver_fixture, timesteps_8_days, tmp
             fs_segmented.solution['objective'].item(),
             rtol=1e-6,
         )
+
+
+class TestCombineSlices:
+    """Tests for the combine_slices utility function."""
+
+    def test_single_dim(self):
+        """Test combining slices with a single extra dimension."""
+        from flixopt.clustering.base import combine_slices
+
+        slices = {
+            ('A',): np.array([1.0, 2.0, 3.0]),
+            ('B',): np.array([4.0, 5.0, 6.0]),
+        }
+        result = combine_slices(
+            slices,
+            extra_dims=['x'],
+            dim_coords={'x': ['A', 'B']},
+            output_dim='time',
+            output_coord=[0, 1, 2],
+        )
+
+        assert result.dims == ('time', 'x')
+        assert result.shape == (3, 2)
+        assert result.sel(x='A').values.tolist() == [1.0, 2.0, 3.0]
+        assert result.sel(x='B').values.tolist() == [4.0, 5.0, 6.0]
+
+    def test_two_dims(self):
+        """Test combining slices with two extra dimensions."""
+        from flixopt.clustering.base import combine_slices
+
+        slices = {
+            ('P1', 'base'): np.array([1.0, 2.0]),
+            ('P1', 'high'): np.array([3.0, 4.0]),
+            ('P2', 'base'): np.array([5.0, 6.0]),
+            ('P2', 'high'): np.array([7.0, 8.0]),
+        }
+        result = combine_slices(
+            slices,
+            extra_dims=['period', 'scenario'],
+            dim_coords={'period': ['P1', 'P2'], 'scenario': ['base', 'high']},
+            output_dim='time',
+            output_coord=[0, 1],
+        )
+
+        assert result.dims == ('time', 'period', 'scenario')
+        assert result.shape == (2, 2, 2)
+        assert result.sel(period='P1', scenario='base').values.tolist() == [1.0, 2.0]
+        assert result.sel(period='P2', scenario='high').values.tolist() == [7.0, 8.0]
+
+    def test_attrs_propagation(self):
+        """Test that attrs are propagated to the result."""
+        from flixopt.clustering.base import combine_slices
+
+        slices = {('A',): np.array([1.0, 2.0])}
+        result = combine_slices(
+            slices,
+            extra_dims=['x'],
+            dim_coords={'x': ['A']},
+            output_dim='time',
+            output_coord=[0, 1],
+            attrs={'units': 'kW', 'description': 'power'},
+        )
+
+        assert result.attrs['units'] == 'kW'
+        assert result.attrs['description'] == 'power'
+
+    def test_datetime_coords(self):
+        """Test with pandas DatetimeIndex as output coordinates."""
+        from flixopt.clustering.base import combine_slices
+
+        time_index = pd.date_range('2020-01-01', periods=3, freq='h')
+        slices = {('A',): np.array([1.0, 2.0, 3.0])}
+        result = combine_slices(
+            slices,
+            extra_dims=['x'],
+            dim_coords={'x': ['A']},
+            output_dim='time',
+            output_coord=time_index,
+        )
+
+        assert result.dims == ('time', 'x')
+        assert len(result.coords['time']) == 3
+        assert result.coords['time'][0].values == time_index[0]

From 77ea9dbc370686f240ad6e81d2231c13b27cb996 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Thu, 15 Jan 2026 11:20:33 +0100
Subject: [PATCH 05/18]  Here's what we accomplished:
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

  1. Fully Vectorized expand_data()

  Before (~65 lines with loops):
  for combo in np.ndindex(*[len(v) for v in dim_coords.values()]):
      selector = {...}
      mapping = _select_dims(timestep_mapping, **selector).values
      data_slice = _select_dims(aggregated, **selector)
      slices[key] = _expand_slice(mapping, data_slice)
  return combine_slices(slices, ...)

  After (~25 lines, fully vectorized):
  timestep_mapping = self.timestep_mapping  # Already multi-dimensional!
  cluster_indices = timestep_mapping // time_dim_size
  time_indices = timestep_mapping % time_dim_size
  expanded = aggregated.isel(cluster=cluster_indices, time=time_indices)
  # xarray handles broadcasting across period/scenario automatically

  2. build_expansion_divisor() and _interpolate_charge_state_segmented()

  These still use combine_slices() because they need per-result segment data (segment_assignments, segment_durations) which isn't available as concatenated Clustering properties yet.

  Current State
  ┌───────────────────────────────────────┬─────────────────┬─────────────────────────────────┐
  │                Method                 │   Vectorized?   │   Uses Clustering Properties    │
  ├───────────────────────────────────────┼─────────────────┼─────────────────────────────────┤
  │ expand_data()                         │ Yes             │ timestep_mapping (fully)        │
  ├───────────────────────────────────────┼─────────────────┼─────────────────────────────────┤
  │ build_expansion_divisor()             │ No (small loop) │ cluster_assignments (partially) │
  ├───────────────────────────────────────┼─────────────────┼─────────────────────────────────┤
  │ _interpolate_charge_state_segmented() │ No (small loop) │ cluster_assignments (partially) │
  └───────────────────────────────────────┴─────────────────┴─────────────────────────────────┘
---
 flixopt/clustering/base.py | 64 ++++++++++++++------------------------
 1 file changed, 23 insertions(+), 41 deletions(-)

diff --git a/flixopt/clustering/base.py b/flixopt/clustering/base.py
index fec150bcb..13b24d05d 100644
--- a/flixopt/clustering/base.py
+++ b/flixopt/clustering/base.py
@@ -905,7 +905,8 @@ def expand_data(
         """Expand aggregated data back to original timesteps.
 
         Uses the timestep_mapping to map each original timestep to its
-        representative value from the aggregated data.
+        representative value from the aggregated data. Fully vectorized using
+        xarray's advanced indexing - no loops over period/scenario dimensions.
 
         Args:
             aggregated: DataArray with aggregated (cluster, time) or (time,) dimension.
@@ -917,50 +918,31 @@ def expand_data(
         if original_time is None:
             original_time = self.original_timesteps
 
-        timestep_mapping = self.timestep_mapping
-        has_cluster_dim = 'cluster' in aggregated.dims
+        timestep_mapping = self.timestep_mapping  # Already multi-dimensional DataArray
 
-        # For segmented systems, the time dimension size is n_segments, not timesteps_per_cluster.
-        # The timestep_mapping uses timesteps_per_cluster for creating indices, but when
-        # indexing into aggregated data with (cluster, time) shape, we need the actual
-        # time dimension size.
-        if has_cluster_dim and self.is_segmented and self.n_segments is not None:
-            time_dim_size = self.n_segments
+        if 'cluster' not in aggregated.dims:
+            # No cluster dimension: use mapping directly as time index
+            expanded = aggregated.isel(time=timestep_mapping)
         else:
-            time_dim_size = self.timesteps_per_cluster
-
-        def _expand_slice(mapping: np.ndarray, data: xr.DataArray) -> np.ndarray:
-            """Expand a single slice using the mapping."""
-            if has_cluster_dim:
-                cluster_ids = mapping // time_dim_size
-                time_within = mapping % time_dim_size
-                return data.values[cluster_ids, time_within]
-            return data.values[mapping]
-
-        # Simple case: no period/scenario dimensions
-        extra_dims = [d for d in timestep_mapping.dims if d != 'original_time']
-        if not extra_dims:
-            expanded_values = _expand_slice(timestep_mapping.values, aggregated)
-            return xr.DataArray(
-                expanded_values,
-                coords={'time': original_time},
-                dims=['time'],
-                attrs=aggregated.attrs,
-            )
+            # Has cluster dimension: compute cluster and time indices from mapping
+            # For segmented systems, time dimension is n_segments, not timesteps_per_cluster
+            if self.is_segmented and self.n_segments is not None:
+                time_dim_size = self.n_segments
+            else:
+                time_dim_size = self.timesteps_per_cluster
 
-        # Multi-dimensional: expand each slice and combine
-        dim_coords = {d: list(timestep_mapping.coords[d].values) for d in extra_dims}
-        slices = {}
-        for combo in np.ndindex(*[len(v) for v in dim_coords.values()]):
-            selector = {d: dim_coords[d][i] for d, i in zip(extra_dims, combo, strict=True)}
-            key = tuple(selector.values())
-            mapping = _select_dims(timestep_mapping, **selector).values
-            data_slice = (
-                _select_dims(aggregated, **selector) if any(d in aggregated.dims for d in selector) else aggregated
-            )
-            slices[key] = _expand_slice(mapping, data_slice)
+            cluster_indices = timestep_mapping // time_dim_size
+            time_indices = timestep_mapping % time_dim_size
+
+            # xarray's advanced indexing handles broadcasting across period/scenario dims
+            expanded = aggregated.isel(cluster=cluster_indices, time=time_indices)
+
+        # Clean up: drop coordinate artifacts from isel, then rename original_time -> time
+        # The isel operation may leave 'cluster' and 'time' as non-dimension coordinates
+        expanded = expanded.drop_vars(['cluster', 'time'], errors='ignore')
+        expanded = expanded.rename({'original_time': 'time'}).assign_coords(time=original_time)
 
-        return combine_slices(slices, extra_dims, dim_coords, 'time', original_time, aggregated.attrs)
+        return expanded.transpose('time', ...).assign_attrs(aggregated.attrs)
 
     def build_expansion_divisor(
         self,

From b146d3e071f5f4ee920ae891fda9b963747502e1 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Thu, 15 Jan 2026 11:30:09 +0100
Subject: [PATCH 06/18]   Completed:

  1. _interpolate_charge_state_segmented() - Fully vectorized from ~110 lines to ~55 lines
    - Uses clustering.timestep_mapping for indexing
    - Uses clustering.results.segment_assignments, segment_durations, and position_within_segment
    - Single xarray expression instead of triple-nested loops

  Previously completed (from before context limit):
  - Added segment_assignments multi-dimensional property to ClusteringResults
  - Added segment_durations multi-dimensional property to ClusteringResults
  - Added position_within_segment property to ClusteringResults
  - Vectorized expand_data()
  - Vectorized build_expansion_divisor()

  Test results: All 130 tests pass (87 cluster/expand + 43 IO tests)

  The combine_slices utility function is still available in clustering/base.py if needed in the future, but all the main dimension-handling methods now use xarray's vectorized advanced indexing instead of the loop-based slice-and-combine pattern.
---
 flixopt/clustering/base.py    | 186 ++++++++++++++++++++--------------
 flixopt/transform_accessor.py | 128 +++++++----------------
 2 files changed, 148 insertions(+), 166 deletions(-)

diff --git a/flixopt/clustering/base.py b/flixopt/clustering/base.py
index 13b24d05d..67c2ca94f 100644
--- a/flixopt/clustering/base.py
+++ b/flixopt/clustering/base.py
@@ -413,25 +413,36 @@ def segment_assignments(self) -> xr.DataArray | None:
         Only available if segmentation was configured during clustering.
 
         Returns:
-            DataArray with dims [cluster, time] or None if no segmentation.
+            DataArray with dims [cluster, time] or [cluster, time, period?, scenario?].
+            Returns None if no segmentation.
         """
         first = self._first_result
         if first.segment_assignments is None:
             return None
 
+        periods = self._get_dim_values('period')
+        scenarios = self._get_dim_values('scenario')
+        timesteps_per_cluster = first.n_timesteps_per_period
+
         if not self.dim_names:
-            # segment_assignments is tuple of tuples: (cluster0_assignments, cluster1_assignments, ...)
+            # Simple case: segment_assignments is tuple of tuples
             data = np.array(first.segment_assignments)
             return xr.DataArray(
                 data,
                 dims=['cluster', 'time'],
-                coords={'cluster': range(self.n_clusters)},
+                coords={'cluster': range(self.n_clusters), 'time': range(timesteps_per_cluster)},
                 name='segment_assignments',
             )
 
-        # Multi-dim case would need more complex handling
-        # For now, return None for multi-dim
-        return None
+        # Multi-dimensional case: build array for each period/scenario
+        return self._build_multi_dim_array(
+            lambda cr: np.array(cr.segment_assignments),
+            base_dims=['cluster', 'time'],
+            base_coords={'cluster': range(self.n_clusters), 'time': range(timesteps_per_cluster)},
+            periods=periods,
+            scenarios=scenarios,
+            name='segment_assignments',
+        )
 
     @property
     def segment_durations(self) -> xr.DataArray | None:
@@ -440,18 +451,25 @@ def segment_durations(self) -> xr.DataArray | None:
         Only available if segmentation was configured during clustering.
 
         Returns:
-            DataArray with dims [cluster, segment] or None if no segmentation.
+            DataArray with dims [cluster, segment] or [cluster, segment, period?, scenario?].
+            Returns None if no segmentation.
         """
         first = self._first_result
         if first.segment_durations is None:
             return None
 
+        periods = self._get_dim_values('period')
+        scenarios = self._get_dim_values('scenario')
+        n_segments = first.n_segments
+
+        def _get_padded_durations(cr: TsamClusteringResult) -> np.ndarray:
+            """Get segment durations, padded with NaN for ragged arrays."""
+            durations = cr.segment_durations
+            return np.array([list(d) + [np.nan] * (n_segments - len(d)) for d in durations])
+
         if not self.dim_names:
-            # segment_durations is tuple of tuples: (cluster0_durations, cluster1_durations, ...)
-            # Each cluster may have different segment counts, so we need to handle ragged arrays
-            durations = first.segment_durations
-            n_segments = first.n_segments
-            data = np.array([list(d) + [np.nan] * (n_segments - len(d)) for d in durations])
+            # Simple case: segment_durations is tuple of tuples
+            data = _get_padded_durations(first)
             return xr.DataArray(
                 data,
                 dims=['cluster', 'segment'],
@@ -460,7 +478,15 @@ def segment_durations(self) -> xr.DataArray | None:
                 attrs={'units': 'hours'},
             )
 
-        return None
+        # Multi-dimensional case: build array for each period/scenario
+        return self._build_multi_dim_array(
+            _get_padded_durations,
+            base_dims=['cluster', 'segment'],
+            base_coords={'cluster': range(self.n_clusters), 'segment': range(n_segments)},
+            periods=periods,
+            scenarios=scenarios,
+            name='segment_durations',
+        )
 
     @property
     def segment_centers(self) -> xr.DataArray | None:
@@ -478,6 +504,59 @@ def segment_centers(self) -> xr.DataArray | None:
         # tsam's segment_centers may be None even with segments configured
         return None
 
+    @property
+    def position_within_segment(self) -> xr.DataArray | None:
+        """Position of each timestep within its segment (0-indexed).
+
+        For each (cluster, time) position, returns how many timesteps into the
+        segment that position is. Used for interpolation within segments.
+
+        Returns:
+            DataArray with dims [cluster, time] or [cluster, time, period?, scenario?].
+            Returns None if no segmentation.
+        """
+        segment_assignments = self.segment_assignments
+        if segment_assignments is None:
+            return None
+
+        def _compute_positions(seg_assigns: np.ndarray) -> np.ndarray:
+            """Compute position within segment for each (cluster, time)."""
+            n_clusters, n_times = seg_assigns.shape
+            positions = np.zeros_like(seg_assigns)
+            for c in range(n_clusters):
+                pos = 0
+                prev_seg = -1
+                for t in range(n_times):
+                    seg = seg_assigns[c, t]
+                    if seg != prev_seg:
+                        pos = 0
+                        prev_seg = seg
+                    positions[c, t] = pos
+                    pos += 1
+            return positions
+
+        # Handle extra dimensions by applying _compute_positions to each slice
+        extra_dims = [d for d in segment_assignments.dims if d not in ('cluster', 'time')]
+
+        if not extra_dims:
+            positions = _compute_positions(segment_assignments.values)
+            return xr.DataArray(
+                positions,
+                dims=['cluster', 'time'],
+                coords=segment_assignments.coords,
+                name='position_within_segment',
+            )
+
+        # Multi-dimensional case: compute for each period/scenario slice
+        result = xr.apply_ufunc(
+            _compute_positions,
+            segment_assignments,
+            input_core_dims=[['cluster', 'time']],
+            output_core_dims=[['cluster', 'time']],
+            vectorize=True,
+        )
+        return result.rename('position_within_segment')
+
     # === Serialization ===
 
     def to_dict(self) -> dict:
@@ -957,6 +1036,8 @@ def build_expansion_divisor(
         For each original timestep, returns the number of original timesteps that map
         to the same (cluster, segment) - i.e., the segment duration in timesteps.
 
+        Fully vectorized using xarray's advanced indexing - no loops over period/scenario.
+
         Args:
             original_time: Original time coordinates. Defaults to self.original_timesteps.
 
@@ -970,72 +1051,29 @@ def build_expansion_divisor(
         if original_time is None:
             original_time = self.original_timesteps
 
-        n_original = self.n_original_clusters
-        timesteps_per_cluster = self.timesteps_per_cluster
-        cluster_assignments = self.cluster_assignments  # Maps original clusters to typical clusters
-
-        def _build_divisor_for_result(
-            assignments: np.ndarray,
-            clustering_result: TsamClusteringResult,
-        ) -> np.ndarray:
-            """Build divisor for a single period/scenario using its clustering result."""
-            n_timesteps = n_original * timesteps_per_cluster
-            divisor = np.empty(n_timesteps)
+        timestep_mapping = self.timestep_mapping  # Already multi-dimensional
+        segment_assignments = self.results.segment_assignments  # [cluster, time, period?, scenario?]
+        segment_durations = self.results.segment_durations  # [cluster, segment, period?, scenario?]
 
-            # Get segment info from the specific clustering result
-            # segment_durations: tuple of tuples, one per cluster
-            # segment_assignments: tuple of tuples, one per cluster
-            seg_durations_per_cluster = clustering_result.segment_durations
-            seg_assignments_per_cluster = clustering_result.segment_assignments
+        # Decode cluster and time indices from timestep_mapping
+        # For segmented systems, time dimension is n_segments
+        time_dim_size = self.n_segments
+        cluster_indices = timestep_mapping // time_dim_size
+        time_indices = timestep_mapping % time_dim_size
 
-            for orig_cluster_idx in range(n_original):
-                typical_cluster_idx = int(assignments[orig_cluster_idx])
+        # Step 1: Get segment index for each original timestep
+        # segment_assignments[cluster, time] -> segment index
+        seg_indices = segment_assignments.isel(cluster=cluster_indices, time=time_indices)
 
-                # Get segment durations and assignments for this typical cluster
-                seg_durs = seg_durations_per_cluster[typical_cluster_idx]
-                seg_assigns = seg_assignments_per_cluster[typical_cluster_idx]
-
-                # For each timestep within this cluster
-                cluster_start = orig_cluster_idx * timesteps_per_cluster
-                for t in range(timesteps_per_cluster):
-                    seg_idx = int(seg_assigns[t])
-                    divisor[cluster_start + t] = seg_durs[seg_idx]
-
-            return divisor
-
-        # Handle extra dimensions (period, scenario)
-        extra_dims = self.results.dim_names
-
-        if not extra_dims:
-            # Simple case: no period/scenario dimensions
-            result = self.results.sel()  # Get the single result
-            divisor_values = _build_divisor_for_result(cluster_assignments.values, result)
-            return xr.DataArray(
-                divisor_values,
-                dims=['time'],
-                coords={'time': original_time},
-                name='expansion_divisor',
-            )
-
-        # Multi-dimensional: build divisor for each period/scenario and combine
-        dim_coords = self.results.coords
-        slices = {}
-
-        for combo in np.ndindex(*[len(v) for v in dim_coords.values()]):
-            selector = {d: dim_coords[d][i] for d, i in zip(extra_dims, combo, strict=True)}
-            key = tuple(selector.values())
-
-            # Get cluster assignments for this period/scenario
-            if any(d in cluster_assignments.dims for d in selector):
-                assignments = _select_dims(cluster_assignments, **selector).values
-            else:
-                assignments = cluster_assignments.values
+        # Step 2: Get duration for each segment
+        # segment_durations[cluster, segment] -> duration
+        divisor = segment_durations.isel(cluster=cluster_indices, segment=seg_indices)
 
-            # Get the clustering result for this period/scenario
-            result = self.results.sel(**selector)
-            slices[key] = _build_divisor_for_result(assignments, result)
+        # Clean up coordinates and rename
+        divisor = divisor.drop_vars(['cluster', 'time', 'segment'], errors='ignore')
+        divisor = divisor.rename({'original_time': 'time'}).assign_coords(time=original_time)
 
-        return combine_slices(slices, extra_dims, dim_coords, 'time', original_time, {'name': 'expansion_divisor'})
+        return divisor.transpose('time', ...).rename('expansion_divisor')
 
     def get_result(
         self,
diff --git a/flixopt/transform_accessor.py b/flixopt/transform_accessor.py
index c13a00402..9c0e7e89e 100644
--- a/flixopt/transform_accessor.py
+++ b/flixopt/transform_accessor.py
@@ -1831,6 +1831,8 @@ def _interpolate_charge_state_segmented(
         this method interpolates between start and end boundary values to show the
         actual charge trajectory as the storage charges/discharges.
 
+        Uses vectorized xarray operations via Clustering class properties.
+
         Args:
             da: charge_state DataArray with dims (cluster, time) where time has n_segments+1 entries.
             clustering: Clustering object with segment info.
@@ -1839,98 +1841,40 @@ def _interpolate_charge_state_segmented(
         Returns:
             Interpolated charge_state with dims (time, ...) for original timesteps.
         """
-        timesteps_per_cluster = clustering.timesteps_per_cluster
-        n_original_clusters = clustering.n_original_clusters
-        cluster_assignments = clustering.cluster_assignments
-
-        # Get segment assignments and durations from clustering results
-        extra_dims = clustering.results.dim_names
-
-        def _interpolate_slice(
-            charge_state_data: np.ndarray,
-            assignments: np.ndarray,
-            clustering_result,
-        ) -> np.ndarray:
-            """Interpolate charge_state for a single period/scenario slice."""
-            n_timesteps = n_original_clusters * timesteps_per_cluster
-            result = np.zeros(n_timesteps)
-
-            seg_assignments = clustering_result.segment_assignments
-            seg_durations = clustering_result.segment_durations
-
-            for orig_cluster_idx in range(n_original_clusters):
-                typical_cluster_idx = int(assignments[orig_cluster_idx])
-                cluster_seg_assigns = seg_assignments[typical_cluster_idx]
-                cluster_seg_durs = seg_durations[typical_cluster_idx]
-
-                # Build cumulative positions within cluster for interpolation
-                for t in range(timesteps_per_cluster):
-                    seg_idx = int(cluster_seg_assigns[t])
-                    # Count how many timesteps before this one are in the same segment
-                    seg_start_t = 0
-                    for prev_t in range(t):
-                        if cluster_seg_assigns[prev_t] == seg_idx:
-                            break
-                        seg_start_t = prev_t + 1
-                    t_within_seg = t - seg_start_t
-                    seg_duration = cluster_seg_durs[seg_idx]
-
-                    # Interpolation factor: position within segment (0 to 1)
-                    # At t_within_seg=0, factor=0 (start of segment)
-                    # At t_within_seg=seg_duration-1, factor approaches 1 (end of segment)
-                    if seg_duration > 1:
-                        factor = (t_within_seg + 0.5) / seg_duration
-                    else:
-                        factor = 0.5
-
-                    # Get start and end boundary values
-                    start_val = charge_state_data[typical_cluster_idx, seg_idx]
-                    end_val = charge_state_data[typical_cluster_idx, seg_idx + 1]
-
-                    # Linear interpolation
-                    result[orig_cluster_idx * timesteps_per_cluster + t] = start_val + (end_val - start_val) * factor
-
-            return result
-
-        # Handle extra dimensions (period, scenario)
-        if not extra_dims:
-            # Simple case: no period/scenario dimensions
-            result = clustering.results.sel()
-            interpolated = _interpolate_slice(da.values, cluster_assignments.values, result)
-            return xr.DataArray(
-                interpolated,
-                dims=['time'],
-                coords={'time': original_timesteps},
-                attrs=da.attrs,
-            )
-
-        # Multi-dimensional case: interpolate each slice and combine
-        from .clustering.base import _select_dims, combine_slices
-
-        dim_coords = clustering.results.coords
-        slices = {}
-
-        for combo in np.ndindex(*[len(v) for v in dim_coords.values()]):
-            selector = {d: dim_coords[d][i] for d, i in zip(extra_dims, combo, strict=True)}
-            key = tuple(selector.values())
-
-            # Get cluster assignments for this period/scenario
-            if any(d in cluster_assignments.dims for d in selector):
-                assignments = _select_dims(cluster_assignments, **selector).values
-            else:
-                assignments = cluster_assignments.values
-
-            # Get charge_state data for this period/scenario
-            da_slice = da
-            for dim, val in selector.items():
-                if dim in da.dims:
-                    da_slice = da_slice.sel({dim: val})
-
-            # Get clustering result for this period/scenario
-            result = clustering.results.sel(**selector)
-            slices[key] = _interpolate_slice(da_slice.values, assignments, result)
-
-        return combine_slices(slices, extra_dims, dim_coords, 'time', original_timesteps, da.attrs)
+        # Get multi-dimensional properties from Clustering
+        timestep_mapping = clustering.timestep_mapping
+        segment_assignments = clustering.results.segment_assignments
+        segment_durations = clustering.results.segment_durations
+        position_within_segment = clustering.results.position_within_segment
+
+        # Decode timestep_mapping into cluster and time indices
+        time_dim_size = clustering.timesteps_per_cluster
+        cluster_indices = timestep_mapping // time_dim_size
+        time_indices = timestep_mapping % time_dim_size
+
+        # Get segment index and position for each original timestep
+        seg_indices = segment_assignments.isel(cluster=cluster_indices, time=time_indices)
+        positions = position_within_segment.isel(cluster=cluster_indices, time=time_indices)
+        durations = segment_durations.isel(cluster=cluster_indices, segment=seg_indices)
+
+        # Calculate interpolation factor: position within segment (0 to 1)
+        # At position=0, factor=0.5/duration (start of segment)
+        # At position=duration-1, factor approaches 1 (end of segment)
+        factor = xr.where(durations > 1, (positions + 0.5) / durations, 0.5)
+
+        # Get start and end boundary values from charge_state
+        # charge_state has dims (cluster, time) where time = segment boundaries (n_segments+1)
+        start_vals = da.isel(cluster=cluster_indices, time=seg_indices)
+        end_vals = da.isel(cluster=cluster_indices, time=seg_indices + 1)
+
+        # Linear interpolation
+        interpolated = start_vals + (end_vals - start_vals) * factor
+
+        # Clean up coordinate artifacts and rename
+        interpolated = interpolated.drop_vars(['cluster', 'time', 'segment'], errors='ignore')
+        interpolated = interpolated.rename({'original_time': 'time'}).assign_coords(time=original_timesteps)
+
+        return interpolated.transpose('time', ...).assign_attrs(da.attrs)
 
     def expand(self) -> FlowSystem:
         """Expand a clustered FlowSystem back to full original timesteps.

From ed221cb9b8d93554e95c28b3a9a08a71a729b61e Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Thu, 15 Jan 2026 13:56:47 +0100
Subject: [PATCH 07/18] All simplifications complete! Here's a summary of what
 we cleaned up:
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

  Summary of Simplifications

  1. expand_da() in transform_accessor.py

  - Extracted duplicate "append extra timestep" logic into _append_final_state() helper
  - Reduced from ~50 lines to ~25 lines
  - Eliminated code duplication

  2. _build_multi_dim_array() → _build_property_array() in clustering/base.py

  - Replaced 6 conditional branches with unified np.ndindex() pattern
  - Now handles both simple and multi-dimensional cases in one method
  - Reduced from ~50 lines to ~25 lines
  - Preserves dtype (fixed integer indexing bug)

  3. Property boilerplate in ClusteringResults

  - 5 properties (cluster_assignments, cluster_occurrences, cluster_centers, segment_assignments, segment_durations) now use the unified _build_property_array()
  - Each property reduced from ~25 lines to ~8 lines
  - Total: ~165 lines → ~85 lines

  4. _build_timestep_mapping() in Clustering

  - Simplified to single call using _build_property_array()
  - Reduced from ~16 lines to ~9 lines

  Total lines removed: ~150+ lines of duplicated/complex code
---
 flixopt/clustering/base.py    | 224 ++++++++--------------------------
 flixopt/transform_accessor.py |  60 ++++-----
 2 files changed, 77 insertions(+), 207 deletions(-)

diff --git a/flixopt/clustering/base.py b/flixopt/clustering/base.py
index 67c2ca94f..6c7af3c11 100644
--- a/flixopt/clustering/base.py
+++ b/flixopt/clustering/base.py
@@ -324,167 +324,82 @@ def n_segments(self) -> int | None:
 
     @property
     def cluster_assignments(self) -> xr.DataArray:
-        """Build multi-dimensional cluster_assignments DataArray.
+        """Maps each original cluster to its typical cluster index.
 
         Returns:
-            DataArray with dims [original_cluster] or [original_cluster, period?, scenario?].
+            DataArray with dims [original_cluster, period?, scenario?].
         """
-        if not self.dim_names:
-            # Simple case: no extra dimensions
-            # Note: Don't include coords - they cause issues when used as isel() indexer
-            return xr.DataArray(
-                np.array(self._results[()].cluster_assignments),
-                dims=['original_cluster'],
-                name='cluster_assignments',
-            )
-
-        # Multi-dimensional case
-        # Note: Don't include coords - they cause issues when used as isel() indexer
-        periods = self._get_dim_values('period')
-        scenarios = self._get_dim_values('scenario')
-
-        return self._build_multi_dim_array(
+        # Note: No coords on original_cluster - they cause issues when used as isel() indexer
+        return self._build_property_array(
             lambda cr: np.array(cr.cluster_assignments),
             base_dims=['original_cluster'],
-            base_coords={},  # No coords on original_cluster
-            periods=periods,
-            scenarios=scenarios,
             name='cluster_assignments',
         )
 
     @property
     def cluster_occurrences(self) -> xr.DataArray:
-        """Build multi-dimensional cluster_occurrences DataArray.
+        """How many original clusters map to each typical cluster.
 
         Returns:
-            DataArray with dims [cluster] or [cluster, period?, scenario?].
+            DataArray with dims [cluster, period?, scenario?].
         """
-        if not self.dim_names:
-            return xr.DataArray(
-                _cluster_occurrences(self._results[()]),
-                dims=['cluster'],
-                coords={'cluster': range(self.n_clusters)},
-                name='cluster_occurrences',
-            )
-
-        periods = self._get_dim_values('period')
-        scenarios = self._get_dim_values('scenario')
-
-        return self._build_multi_dim_array(
+        return self._build_property_array(
             _cluster_occurrences,
             base_dims=['cluster'],
             base_coords={'cluster': range(self.n_clusters)},
-            periods=periods,
-            scenarios=scenarios,
             name='cluster_occurrences',
         )
 
     @property
     def cluster_centers(self) -> xr.DataArray:
-        """Which original period is the representative (center) for each cluster.
+        """Which original cluster is the representative (center) for each typical cluster.
 
         Returns:
-            DataArray with dims [cluster] containing original period indices.
+            DataArray with dims [cluster, period?, scenario?].
         """
-        if not self.dim_names:
-            return xr.DataArray(
-                np.array(self._results[()].cluster_centers),
-                dims=['cluster'],
-                coords={'cluster': range(self.n_clusters)},
-                name='cluster_centers',
-            )
-
-        periods = self._get_dim_values('period')
-        scenarios = self._get_dim_values('scenario')
-
-        return self._build_multi_dim_array(
+        return self._build_property_array(
             lambda cr: np.array(cr.cluster_centers),
             base_dims=['cluster'],
             base_coords={'cluster': range(self.n_clusters)},
-            periods=periods,
-            scenarios=scenarios,
             name='cluster_centers',
         )
 
     @property
     def segment_assignments(self) -> xr.DataArray | None:
-        """For each timestep within a cluster, which intra-period segment it belongs to.
-
-        Only available if segmentation was configured during clustering.
+        """For each timestep within a cluster, which segment it belongs to.
 
         Returns:
-            DataArray with dims [cluster, time] or [cluster, time, period?, scenario?].
-            Returns None if no segmentation.
+            DataArray with dims [cluster, time, period?, scenario?], or None if not segmented.
         """
-        first = self._first_result
-        if first.segment_assignments is None:
+        if self._first_result.segment_assignments is None:
             return None
-
-        periods = self._get_dim_values('period')
-        scenarios = self._get_dim_values('scenario')
-        timesteps_per_cluster = first.n_timesteps_per_period
-
-        if not self.dim_names:
-            # Simple case: segment_assignments is tuple of tuples
-            data = np.array(first.segment_assignments)
-            return xr.DataArray(
-                data,
-                dims=['cluster', 'time'],
-                coords={'cluster': range(self.n_clusters), 'time': range(timesteps_per_cluster)},
-                name='segment_assignments',
-            )
-
-        # Multi-dimensional case: build array for each period/scenario
-        return self._build_multi_dim_array(
+        timesteps = self._first_result.n_timesteps_per_period
+        return self._build_property_array(
             lambda cr: np.array(cr.segment_assignments),
             base_dims=['cluster', 'time'],
-            base_coords={'cluster': range(self.n_clusters), 'time': range(timesteps_per_cluster)},
-            periods=periods,
-            scenarios=scenarios,
+            base_coords={'cluster': range(self.n_clusters), 'time': range(timesteps)},
             name='segment_assignments',
         )
 
     @property
     def segment_durations(self) -> xr.DataArray | None:
-        """Duration of each intra-period segment in hours.
-
-        Only available if segmentation was configured during clustering.
+        """Duration of each segment in timesteps.
 
         Returns:
-            DataArray with dims [cluster, segment] or [cluster, segment, period?, scenario?].
-            Returns None if no segmentation.
+            DataArray with dims [cluster, segment, period?, scenario?], or None if not segmented.
         """
-        first = self._first_result
-        if first.segment_durations is None:
+        if self._first_result.segment_durations is None:
             return None
-
-        periods = self._get_dim_values('period')
-        scenarios = self._get_dim_values('scenario')
-        n_segments = first.n_segments
+        n_segments = self._first_result.n_segments
 
         def _get_padded_durations(cr: TsamClusteringResult) -> np.ndarray:
-            """Get segment durations, padded with NaN for ragged arrays."""
-            durations = cr.segment_durations
-            return np.array([list(d) + [np.nan] * (n_segments - len(d)) for d in durations])
-
-        if not self.dim_names:
-            # Simple case: segment_durations is tuple of tuples
-            data = _get_padded_durations(first)
-            return xr.DataArray(
-                data,
-                dims=['cluster', 'segment'],
-                coords={'cluster': range(self.n_clusters), 'segment': range(n_segments)},
-                name='segment_durations',
-                attrs={'units': 'hours'},
-            )
+            """Pad ragged segment durations to uniform shape."""
+            return np.array([list(d) + [np.nan] * (n_segments - len(d)) for d in cr.segment_durations])
 
-        # Multi-dimensional case: build array for each period/scenario
-        return self._build_multi_dim_array(
+        return self._build_property_array(
             _get_padded_durations,
             base_dims=['cluster', 'segment'],
             base_coords={'cluster': range(self.n_clusters), 'segment': range(n_segments)},
-            periods=periods,
-            scenarios=scenarios,
             name='segment_durations',
         )
 
@@ -605,58 +520,41 @@ def _get_dim_values(self, dim: str) -> list | None:
         idx = self._dim_names.index(dim)
         return sorted(set(k[idx] for k in self._results.keys()))
 
-    def _build_multi_dim_array(
+    def _build_property_array(
         self,
         get_data: callable,
         base_dims: list[str],
-        base_coords: dict,
-        periods: list | None,
-        scenarios: list | None,
-        name: str,
+        base_coords: dict | None = None,
+        name: str | None = None,
     ) -> xr.DataArray:
-        """Build a multi-dimensional DataArray from per-result data."""
-        has_periods = periods is not None
-        has_scenarios = scenarios is not None
-
-        slices = {}
-        if has_periods and has_scenarios:
-            for p in periods:
-                for s in scenarios:
-                    slices[(p, s)] = xr.DataArray(
-                        get_data(self._results[(p, s)]),
-                        dims=base_dims,
-                        coords=base_coords,
-                    )
-        elif has_periods:
-            for p in periods:
-                slices[(p,)] = xr.DataArray(
-                    get_data(self._results[(p,)]),
-                    dims=base_dims,
-                    coords=base_coords,
-                )
-        elif has_scenarios:
-            for s in scenarios:
-                slices[(s,)] = xr.DataArray(
-                    get_data(self._results[(s,)]),
-                    dims=base_dims,
-                    coords=base_coords,
-                )
-
-        # Combine slices into multi-dimensional array
-        if has_periods and has_scenarios:
-            period_arrays = []
-            for p in periods:
-                scenario_arrays = [slices[(p, s)] for s in scenarios]
-                period_arrays.append(xr.concat(scenario_arrays, dim=pd.Index(scenarios, name='scenario')))
-            result = xr.concat(period_arrays, dim=pd.Index(periods, name='period'))
-        elif has_periods:
-            result = xr.concat([slices[(p,)] for p in periods], dim=pd.Index(periods, name='period'))
-        else:
-            result = xr.concat([slices[(s,)] for s in scenarios], dim=pd.Index(scenarios, name='scenario'))
+        """Build a DataArray property, handling both single and multi-dimensional cases."""
+        base_coords = base_coords or {}
+        periods = self._get_dim_values('period')
+        scenarios = self._get_dim_values('scenario')
 
-        # Ensure base dims come first
-        dim_order = base_dims + [d for d in result.dims if d not in base_dims]
-        return result.transpose(*dim_order).rename(name)
+        # Build list of (dim_name, values) for dimensions that exist
+        extra_dims = []
+        if periods is not None:
+            extra_dims.append(('period', periods))
+        if scenarios is not None:
+            extra_dims.append(('scenario', scenarios))
+
+        # Simple case: no extra dimensions
+        if not extra_dims:
+            return xr.DataArray(get_data(self._results[()]), dims=base_dims, coords=base_coords, name=name)
+
+        # Multi-dimensional: stack data for each combination
+        first_data = get_data(next(iter(self._results.values())))
+        shape = list(first_data.shape) + [len(vals) for _, vals in extra_dims]
+        data = np.empty(shape, dtype=first_data.dtype)  # Preserve dtype
+
+        for combo in np.ndindex(*[len(vals) for _, vals in extra_dims]):
+            key = tuple(extra_dims[i][1][idx] for i, idx in enumerate(combo))
+            data[(...,) + combo] = get_data(self._results[key])
+
+        dims = base_dims + [dim_name for dim_name, _ in extra_dims]
+        coords = {**base_coords, **{dim_name: vals for dim_name, vals in extra_dims}}
+        return xr.DataArray(data, dims=dims, coords=coords, name=name)
 
     @staticmethod
     def _key_to_str(key: tuple) -> str:
@@ -1180,24 +1078,10 @@ def _build_timestep_mapping(self) -> xr.DataArray:
         """Build timestep_mapping DataArray."""
         n_original = len(self.original_timesteps)
         original_time_coord = self.original_timesteps.rename('original_time')
-
-        if not self.dim_names:
-            # Simple case: no extra dimensions
-            mapping = _build_timestep_mapping(self.results[()], n_original)
-            return xr.DataArray(
-                mapping,
-                dims=['original_time'],
-                coords={'original_time': original_time_coord},
-                name='timestep_mapping',
-            )
-
-        # Multi-dimensional case: combine slices into multi-dim array
-        return self.results._build_multi_dim_array(
+        return self.results._build_property_array(
             lambda cr: _build_timestep_mapping(cr, n_original),
             base_dims=['original_time'],
             base_coords={'original_time': original_time_coord},
-            periods=self.results._get_dim_values('period'),
-            scenarios=self.results._get_dim_values('scenario'),
             name='timestep_mapping',
         )
 
diff --git a/flixopt/transform_accessor.py b/flixopt/transform_accessor.py
index 9c0e7e89e..8b78540e1 100644
--- a/flixopt/transform_accessor.py
+++ b/flixopt/transform_accessor.py
@@ -2015,54 +2015,40 @@ def _is_state_variable(var_name: str) -> bool:
             # Fall back to pattern matching for backwards compatibility
             return var_name.endswith('|charge_state')
 
+        def _append_final_state(expanded: xr.DataArray, da: xr.DataArray) -> xr.DataArray:
+            """Append final state value from original data to expanded data."""
+            cluster_assignments = clustering.cluster_assignments
+            if cluster_assignments.ndim == 1:
+                last_cluster = int(cluster_assignments.values[last_original_cluster_idx])
+                extra_val = da.isel(cluster=last_cluster, time=-1)
+            else:
+                last_clusters = cluster_assignments.isel(original_cluster=last_original_cluster_idx)
+                extra_val = da.isel(cluster=last_clusters, time=-1)
+            extra_val = extra_val.drop_vars(['cluster', 'time'], errors='ignore')
+            extra_val = extra_val.expand_dims(time=[original_timesteps_extra[-1]])
+            return xr.concat([expanded, extra_val], dim='time')
+
         def expand_da(da: xr.DataArray, var_name: str = '', is_solution: bool = False) -> xr.DataArray:
-            """Expand a DataArray from clustered to original timesteps.
-
-            Args:
-                da: DataArray to expand.
-                var_name: Variable name for segment total lookup.
-                is_solution: True if this is a solution variable (may need segment correction).
-                    FlowSystem data (is_solution=False) is never corrected for segments.
-            """
+            """Expand a DataArray from clustered to original timesteps."""
             if 'time' not in da.dims:
                 return da.copy()
 
-            # For state variables (like charge_state) in segmented systems: interpolate within segments
-            # to show the actual state trajectory as the storage charges/discharges
-            if _is_state_variable(var_name) and 'cluster' in da.dims and clustering.is_segmented:
+            is_state = _is_state_variable(var_name) and 'cluster' in da.dims
+
+            # State variables in segmented systems: interpolate within segments
+            if is_state and clustering.is_segmented:
                 expanded = self._interpolate_charge_state_segmented(da, clustering, original_timesteps)
-                # Append the extra timestep value (final charge state)
-                cluster_assignments = clustering.cluster_assignments
-                if cluster_assignments.ndim == 1:
-                    last_cluster = int(cluster_assignments[last_original_cluster_idx])
-                    extra_val = da.isel(cluster=last_cluster, time=-1)
-                else:
-                    last_clusters = cluster_assignments.isel(original_cluster=last_original_cluster_idx)
-                    extra_val = da.isel(cluster=last_clusters, time=-1)
-                extra_val = extra_val.drop_vars(['cluster', 'time'], errors='ignore')
-                extra_val = extra_val.expand_dims(time=[original_timesteps_extra[-1]])
-                expanded = xr.concat([expanded, extra_val], dim='time')
-                return expanded
+                return _append_final_state(expanded, da)
 
             expanded = clustering.expand_data(da, original_time=original_timesteps)
 
-            # For segmented systems: divide segment totals by expansion divisor
-            # ONLY for solution variables explicitly identified as segment totals
+            # Segment totals: divide by expansion divisor
             if is_solution and expansion_divisor is not None and var_name in segment_total_vars:
                 expanded = expanded / expansion_divisor
 
-            # For state variables (like charge_state) with cluster dim (non-segmented), append the extra timestep value
-            if _is_state_variable(var_name) and 'cluster' in da.dims:
-                cluster_assignments = clustering.cluster_assignments
-                if cluster_assignments.ndim == 1:
-                    last_cluster = int(cluster_assignments[last_original_cluster_idx])
-                    extra_val = da.isel(cluster=last_cluster, time=-1)
-                else:
-                    last_clusters = cluster_assignments.isel(original_cluster=last_original_cluster_idx)
-                    extra_val = da.isel(cluster=last_clusters, time=-1)
-                extra_val = extra_val.drop_vars(['cluster', 'time'], errors='ignore')
-                extra_val = extra_val.expand_dims(time=[original_timesteps_extra[-1]])
-                expanded = xr.concat([expanded, extra_val], dim='time')
+            # State variables: append final state
+            if is_state:
+                expanded = _append_final_state(expanded, da)
 
             return expanded
 

From 64ce6ce2fa8bfc6f4d2c46b9c6eb276a6f291633 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Thu, 15 Jan 2026 14:20:50 +0100
Subject: [PATCH 08/18] Removed the unnecessary lookup and use segment_indices
 directl

---
 flixopt/clustering/base.py | 15 +++++----------
 1 file changed, 5 insertions(+), 10 deletions(-)

diff --git a/flixopt/clustering/base.py b/flixopt/clustering/base.py
index 6c7af3c11..51d6a7fa4 100644
--- a/flixopt/clustering/base.py
+++ b/flixopt/clustering/base.py
@@ -950,22 +950,17 @@ def build_expansion_divisor(
             original_time = self.original_timesteps
 
         timestep_mapping = self.timestep_mapping  # Already multi-dimensional
-        segment_assignments = self.results.segment_assignments  # [cluster, time, period?, scenario?]
         segment_durations = self.results.segment_durations  # [cluster, segment, period?, scenario?]
 
-        # Decode cluster and time indices from timestep_mapping
-        # For segmented systems, time dimension is n_segments
+        # Decode cluster and segment indices from timestep_mapping
+        # For segmented systems, encoding is: cluster_id * n_segments + segment_idx
         time_dim_size = self.n_segments
         cluster_indices = timestep_mapping // time_dim_size
-        time_indices = timestep_mapping % time_dim_size
+        segment_indices = timestep_mapping % time_dim_size  # This IS the segment index
 
-        # Step 1: Get segment index for each original timestep
-        # segment_assignments[cluster, time] -> segment index
-        seg_indices = segment_assignments.isel(cluster=cluster_indices, time=time_indices)
-
-        # Step 2: Get duration for each segment
+        # Get duration for each segment directly
         # segment_durations[cluster, segment] -> duration
-        divisor = segment_durations.isel(cluster=cluster_indices, segment=seg_indices)
+        divisor = segment_durations.isel(cluster=cluster_indices, segment=segment_indices)
 
         # Clean up coordinates and rename
         divisor = divisor.drop_vars(['cluster', 'time', 'segment'], errors='ignore')

From ff43e32f3f5ba055157684cc737ea54fd1c85c76 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Thu, 15 Jan 2026 14:41:29 +0100
Subject: [PATCH 09/18] The IO roundtrip fix is working correctly. Here's a
 summary of what was fixed:
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

  Summary

  The IO roundtrip bug was caused by representative_weights (a variable with only ('cluster',) dimension) being copied as-is during expansion, which caused the cluster dimension to incorrectly persist in the expanded dataset.

  Fix applied in transform_accessor.py:2063-2065:
  # Skip cluster-only vars (no time dim) - they don't make sense after expansion
  if da.dims == ('cluster',):
      continue

  This skips variables that have only a cluster dimension (and no time dimension) during expansion, as these variables don't make sense after the clustering structure is removed.

  Test results:
  - All 87 tests in test_cluster_reduce_expand.py pass ✓
  - All 43 tests in test_clustering_io.py pass ✓
  - Manual IO roundtrip test passes ✓
  - Tests with different segment counts (3, 6) pass ✓
  - Tests with 2-hour timesteps pass ✓
---
 flixopt/transform_accessor.py | 14 +++++++++-----
 1 file changed, 9 insertions(+), 5 deletions(-)

diff --git a/flixopt/transform_accessor.py b/flixopt/transform_accessor.py
index 8b78540e1..0eeb0a4ae 100644
--- a/flixopt/transform_accessor.py
+++ b/flixopt/transform_accessor.py
@@ -2055,11 +2055,15 @@ def expand_da(da: xr.DataArray, var_name: str = '', is_solution: bool = False) -
         # 1. Expand FlowSystem data
         reduced_ds = self._fs.to_dataset(include_solution=False)
         clustering_attrs = {'is_clustered', 'n_clusters', 'timesteps_per_cluster', 'clustering', 'cluster_weight'}
-        data_vars = {
-            name: expand_da(da, name)
-            for name, da in reduced_ds.data_vars.items()
-            if name != 'cluster_weight' and not name.startswith('clustering|')
-        }
+        skip_vars = {'cluster_weight', 'timestep_duration'}  # These have special handling
+        data_vars = {}
+        for name, da in reduced_ds.data_vars.items():
+            if name in skip_vars or name.startswith('clustering|'):
+                continue
+            # Skip cluster-only vars (no time dim) - they don't make sense after expansion
+            if da.dims == ('cluster',):
+                continue
+            data_vars[name] = expand_da(da, name)
         attrs = {k: v for k, v in reduced_ds.attrs.items() if k not in clustering_attrs}
         expanded_ds = xr.Dataset(data_vars, attrs=attrs)
 

From 3bdcd1d9f5e4105eb0c4333cdde72fba8546ebf2 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Thu, 15 Jan 2026 14:44:24 +0100
Subject: [PATCH 10/18] Updated condition in transform_accessor.py:2063-2066:

  # Skip vars with cluster dim but no time dim - they don't make sense after expansion
  # (e.g., representative_weights with dims ('cluster',) or ('cluster', 'period'))
  if 'cluster' in da.dims and 'time' not in da.dims:
      continue

  This correctly handles:
  - ('cluster',) - simple cluster-only variables like cluster_weight
  - ('cluster', 'period') - cluster variables with period dimension
  - ('cluster', 'scenario') - cluster variables with scenario dimension
  - ('cluster', 'period', 'scenario') - cluster variables with both

  Variables with both cluster and time dimensions (like timestep_duration with dims ('cluster', 'time')) are correctly expanded since they contain time-series data that needs to be mapped back to original timesteps.
---
 flixopt/transform_accessor.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/flixopt/transform_accessor.py b/flixopt/transform_accessor.py
index 0eeb0a4ae..1a3f16d8a 100644
--- a/flixopt/transform_accessor.py
+++ b/flixopt/transform_accessor.py
@@ -2060,8 +2060,9 @@ def expand_da(da: xr.DataArray, var_name: str = '', is_solution: bool = False) -
         for name, da in reduced_ds.data_vars.items():
             if name in skip_vars or name.startswith('clustering|'):
                 continue
-            # Skip cluster-only vars (no time dim) - they don't make sense after expansion
-            if da.dims == ('cluster',):
+            # Skip vars with cluster dim but no time dim - they don't make sense after expansion
+            # (e.g., representative_weights with dims ('cluster',) or ('cluster', 'period'))
+            if 'cluster' in da.dims and 'time' not in da.dims:
                 continue
             data_vars[name] = expand_da(da, name)
         attrs = {k: v for k, v in reduced_ds.attrs.items() if k not in clustering_attrs}

From 94e0a1d8108dc9e1100b513edece2049ef21b14b Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Thu, 15 Jan 2026 16:01:56 +0100
Subject: [PATCH 11/18]   Summary of Fixes
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

  1. clustering/base.py - combine_slices() hardening (lines 52-118)

  - Added validation for empty input: if not slices: raise ValueError("slices cannot be empty")
  - Capture first array and preserve dtype: first = next(iter(slices.values())) → np.empty(shape, dtype=first.dtype)
  - Clearer error on missing keys with try/except: raise KeyError(f"Missing slice for key {key} (extra_dims={extra_dims})")

  2. flow_system.py - Variable categories cleanup and safe enum restoration

  - Added self._variable_categories.clear() in _invalidate_model() (line 1692) to prevent stale categories from being reused
  - Hardened VariableCategory restoration (lines 922-930) with try/except to handle unknown/renamed enum values gracefully with a warning instead of crashing

  3. transform_accessor.py - Correct timestep_mapping decode for segmented systems (lines 1850-1857)

  - For segmented systems, now uses clustering.n_segments instead of clustering.timesteps_per_cluster as the divisor
  - This matches the encoding logic in expand_data() and build_expansion_divisor()
---
 flixopt/clustering/base.py    | 17 ++++++++++++++---
 flixopt/flow_system.py        | 15 +++++++++++----
 flixopt/transform_accessor.py |  6 +++++-
 3 files changed, 30 insertions(+), 8 deletions(-)

diff --git a/flixopt/clustering/base.py b/flixopt/clustering/base.py
index 51d6a7fa4..2c6c70b6b 100644
--- a/flixopt/clustering/base.py
+++ b/flixopt/clustering/base.py
@@ -74,6 +74,10 @@ def combine_slices(
     Returns:
         DataArray with dims [output_dim, *extra_dims].
 
+    Raises:
+        ValueError: If slices is empty.
+        KeyError: If a required key is missing from slices.
+
     Example:
         >>> slices = {
         ...     ('P1', 'base'): np.array([1, 2, 3]),
@@ -91,13 +95,20 @@ def combine_slices(
         >>> result.dims
         ('time', 'period', 'scenario')
     """
-    n_output = len(next(iter(slices.values())))
+    if not slices:
+        raise ValueError('slices cannot be empty')
+
+    first = next(iter(slices.values()))
+    n_output = len(first)
     shape = [n_output] + [len(dim_coords[d]) for d in extra_dims]
-    data = np.empty(shape)
+    data = np.empty(shape, dtype=first.dtype)
 
     for combo in np.ndindex(*shape[1:]):
         key = tuple(dim_coords[d][i] for d, i in zip(extra_dims, combo, strict=True))
-        data[(slice(None),) + combo] = slices[key]
+        try:
+            data[(slice(None),) + combo] = slices[key]
+        except KeyError:
+            raise KeyError(f'Missing slice for key {key} (extra_dims={extra_dims})') from None
 
     return xr.DataArray(
         data,
diff --git a/flixopt/flow_system.py b/flixopt/flow_system.py
index 155d8379b..89466e7b6 100644
--- a/flixopt/flow_system.py
+++ b/flixopt/flow_system.py
@@ -918,10 +918,16 @@ def from_dataset(cls, ds: xr.Dataset) -> FlowSystem:
         # Restore variable categories if present
         if 'variable_categories' in reference_structure:
             categories_dict = json.loads(reference_structure['variable_categories'])
-            # Convert string values back to VariableCategory enum
-            flow_system._variable_categories = {
-                name: VariableCategory(value) for name, value in categories_dict.items()
-            }
+            # Convert string values back to VariableCategory enum with safe fallback
+            restored_categories = {}
+            for name, value in categories_dict.items():
+                try:
+                    restored_categories[name] = VariableCategory(value)
+                except ValueError:
+                    # Unknown category value (e.g., renamed/removed enum) - skip it
+                    # The variable will be treated as uncategorized during expansion
+                    logger.warning(f'Unknown VariableCategory value "{value}" for "{name}", skipping')
+            flow_system._variable_categories = restored_categories
 
         # Reconnect network to populate bus inputs/outputs (not stored in NetCDF).
         flow_system.connect_and_transform()
@@ -1689,6 +1695,7 @@ def _invalidate_model(self) -> None:
         self._connected_and_transformed = False
         self._topology = None  # Invalidate topology accessor (and its cached colors)
         self._flow_carriers = None  # Invalidate flow-to-carrier mapping
+        self._variable_categories.clear()  # Clear stale categories for segment expansion
         for element in self.values():
             element.submodel = None
             element._variable_names = []
diff --git a/flixopt/transform_accessor.py b/flixopt/transform_accessor.py
index 1a3f16d8a..83078c3b7 100644
--- a/flixopt/transform_accessor.py
+++ b/flixopt/transform_accessor.py
@@ -1848,7 +1848,11 @@ def _interpolate_charge_state_segmented(
         position_within_segment = clustering.results.position_within_segment
 
         # Decode timestep_mapping into cluster and time indices
-        time_dim_size = clustering.timesteps_per_cluster
+        # For segmented systems, use n_segments as the divisor (matches expand_data/build_expansion_divisor)
+        if clustering.is_segmented and clustering.n_segments is not None:
+            time_dim_size = clustering.n_segments
+        else:
+            time_dim_size = clustering.timesteps_per_cluster
         cluster_indices = timestep_mapping // time_dim_size
         time_indices = timestep_mapping % time_dim_size
 

From 3cddeb8f674e90fdf1330f7e45557ead0b6b9af1 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Thu, 15 Jan 2026 16:03:19 +0100
Subject: [PATCH 12/18]  Added test_segmented_total_effects_match_solution to
 TestSegmentation class

---
 tests/test_cluster_reduce_expand.py | 45 +++++++++++++++++++++++++++++
 1 file changed, 45 insertions(+)

diff --git a/tests/test_cluster_reduce_expand.py b/tests/test_cluster_reduce_expand.py
index c187358d7..d6c991783 100644
--- a/tests/test_cluster_reduce_expand.py
+++ b/tests/test_cluster_reduce_expand.py
@@ -1180,6 +1180,51 @@ def test_segmented_timestep_mapping_uses_segment_assignments(self, timesteps_8_d
         assert mapping.min() >= 0
         assert mapping.max() <= max_valid_idx
 
+    @pytest.mark.parametrize('freq', ['1h', '2h'])
+    def test_segmented_total_effects_match_solution(self, solver_fixture, freq):
+        """Test that total_effects matches solution Cost after expand with segmentation.
+
+        This is a regression test for the bug where expansion_divisor was computed
+        incorrectly for segmented systems, causing total_effects to not match the
+        solution's objective value.
+        """
+        from tsam.config import SegmentConfig
+
+        # Create system with specified timestep frequency
+        n_timesteps = 72 if freq == '1h' else 36  # 3 days worth
+        timesteps = pd.date_range('2024-01-01', periods=n_timesteps, freq=freq)
+        fs = fx.FlowSystem(timesteps=timesteps)
+
+        # Minimal components: effect + source + sink with varying demand
+        fs.add_elements(fx.Effect('Cost', unit='EUR', is_objective=True))
+        fs.add_elements(fx.Bus('Heat'))
+        fs.add_elements(
+            fx.Source(
+                'Boiler',
+                outputs=[fx.Flow('Q', bus='Heat', size=100, effects_per_flow_hour={'Cost': 50})],
+            )
+        )
+        demand_profile = np.tile([0.5, 1], n_timesteps // 2)
+        fs.add_elements(
+            fx.Sink('Demand', inputs=[fx.Flow('Q', bus='Heat', size=50, fixed_relative_profile=demand_profile)])
+        )
+
+        # Cluster with segments -> solve -> expand
+        fs_clustered = fs.transform.cluster(
+            n_clusters=2,
+            cluster_duration='1D',
+            segments=SegmentConfig(n_segments=4),
+        )
+        fs_clustered.optimize(solver_fixture)
+        fs_expanded = fs_clustered.transform.expand()
+
+        # Validate: total_effects must match solution objective
+        computed = fs_expanded.statistics.total_effects['Cost'].sum('contributor')
+        expected = fs_expanded.solution['Cost']
+        assert np.allclose(computed.values, expected.values, rtol=1e-5), (
+            f'total_effects mismatch: computed={float(computed):.2f}, expected={float(expected):.2f}'
+        )
+
 
 class TestSegmentationWithStorage:
     """Tests for segmentation combined with storage components."""

From dd671f2f47e0d99f63dea4a0976ba1e695a592f5 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Thu, 15 Jan 2026 19:42:02 +0100
Subject: [PATCH 13/18] Added all remaining tsam.aggregate() paramaters and
 missing type hint

---
 flixopt/structure.py          |  2 +-
 flixopt/transform_accessor.py | 38 ++++++++++++++++++++++++++++-------
 2 files changed, 32 insertions(+), 8 deletions(-)

diff --git a/flixopt/structure.py b/flixopt/structure.py
index 8ac18e76a..7d76c760c 100644
--- a/flixopt/structure.py
+++ b/flixopt/structure.py
@@ -1671,7 +1671,7 @@ def add_variables(
         self,
         short_name: str = None,
         category: VariableCategory = None,
-        **kwargs,
+        **kwargs: Any,
     ) -> linopy.Variable:
         """Create and register a variable in one step.
 
diff --git a/flixopt/transform_accessor.py b/flixopt/transform_accessor.py
index 83078c3b7..3f7841085 100644
--- a/flixopt/transform_accessor.py
+++ b/flixopt/transform_accessor.py
@@ -1172,6 +1172,10 @@ def cluster(
         cluster: ClusterConfig | None = None,
         extremes: ExtremeConfig | None = None,
         segments: SegmentConfig | None = None,
+        preserve_column_means: bool = True,
+        rescale_exclude_columns: list[str] | None = None,
+        round_decimals: int | None = None,
+        numerical_tolerance: float = 1e-13,
         **tsam_kwargs: Any,
     ) -> FlowSystem:
         """
@@ -1212,8 +1216,19 @@ def cluster(
             segments: Optional tsam ``SegmentConfig`` object specifying intra-period
                 segmentation. Segments divide each cluster period into variable-duration
                 sub-segments. Example: ``SegmentConfig(n_segments=4)``.
-            **tsam_kwargs: Additional keyword arguments passed to ``tsam.aggregate()``.
-                See tsam documentation for all options (e.g., ``preserve_column_means``).
+            preserve_column_means: Rescale typical periods so each column's weighted mean
+                matches the original data's mean. Ensures total energy/load is preserved
+                when weights represent occurrence counts. Default is True.
+            rescale_exclude_columns: Column names to exclude from rescaling when
+                ``preserve_column_means=True``. Useful for binary/indicator columns (0/1 values)
+                that should not be rescaled.
+            round_decimals: Round output values to this many decimal places.
+                If None (default), no rounding is applied.
+            numerical_tolerance: Tolerance for numerical precision issues. Controls when
+                warnings are raised for aggregated values exceeding original time series bounds.
+                Default is 1e-13.
+            **tsam_kwargs: Additional keyword arguments passed to ``tsam.aggregate()``
+                for forward compatibility. See tsam documentation for all options.
 
         Returns:
             A new FlowSystem with reduced timesteps (only typical clusters).
@@ -1302,11 +1317,16 @@ def cluster(
 
         # Validate tsam_kwargs doesn't override explicit parameters
         reserved_tsam_keys = {
-            'n_periods',
-            'period_hours',
-            'resolution',
-            'cluster',  # ClusterConfig object (weights are passed through this)
-            'extremes',  # ExtremeConfig object
+            'n_clusters',
+            'period_duration',  # exposed as cluster_duration
+            'timestep_duration',  # computed automatically
+            'cluster',
+            'segments',
+            'extremes',
+            'preserve_column_means',
+            'rescale_exclude_columns',
+            'round_decimals',
+            'numerical_tolerance',
         }
         conflicts = reserved_tsam_keys & set(tsam_kwargs.keys())
         if conflicts:
@@ -1359,6 +1379,10 @@ def cluster(
                         cluster=cluster_config,
                         extremes=extremes,
                         segments=segments,
+                        preserve_column_means=preserve_column_means,
+                        rescale_exclude_columns=rescale_exclude_columns,
+                        round_decimals=round_decimals,
+                        numerical_tolerance=numerical_tolerance,
                         **tsam_kwargs,
                     )
 

From fe3ae385e5316dbae656f8b0b81dd2cad83d2b39 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Thu, 15 Jan 2026 19:57:34 +0100
Subject: [PATCH 14/18] Updated expression_tracking_variable
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

  modeling.py:200-242 - Added category: VariableCategory = None parameter and passed it to both add_variables calls.

  Updated Callers
  ┌─────────────┬──────┬─────────────────────────┬────────────────────┐
  │    File     │ Line │        Variable         │      Category      │
  ├─────────────┼──────┼─────────────────────────┼────────────────────┤
  │ features.py │ 208  │ active_hours            │ TOTAL              │
  ├─────────────┼──────┼─────────────────────────┼────────────────────┤
  │ elements.py │ 682  │ total_flow_hours        │ TOTAL              │
  ├─────────────┼──────┼─────────────────────────┼────────────────────┤
  │ elements.py │ 709  │ flow_hours_over_periods │ TOTAL_OVER_PERIODS │
  └─────────────┴──────┴─────────────────────────┴────────────────────┘
  All expression tracking variables now properly register their categories for segment expansion handling. The pattern is consistent: callers specify the appropriate category based on what the tracked expression represents.
---
 flixopt/elements.py | 2 ++
 flixopt/features.py | 1 +
 flixopt/modeling.py | 7 ++++++-
 3 files changed, 9 insertions(+), 1 deletion(-)

diff --git a/flixopt/elements.py b/flixopt/elements.py
index 0fa620bed..c80c263d1 100644
--- a/flixopt/elements.py
+++ b/flixopt/elements.py
@@ -689,6 +689,7 @@ def _do_modeling(self):
             ),
             coords=['period', 'scenario'],
             short_name='total_flow_hours',
+            category=VariableCategory.TOTAL,
         )
 
         # Weighted sum over all periods constraint
@@ -719,6 +720,7 @@ def _do_modeling(self):
                 ),
                 coords=['scenario'],
                 short_name='flow_hours_over_periods',
+                category=VariableCategory.TOTAL_OVER_PERIODS,
             )
 
         # Load factor constraints
diff --git a/flixopt/features.py b/flixopt/features.py
index af1ff1471..9ef400392 100644
--- a/flixopt/features.py
+++ b/flixopt/features.py
@@ -214,6 +214,7 @@ def _do_modeling(self):
             ),
             short_name='active_hours',
             coords=['period', 'scenario'],
+            category=VariableCategory.TOTAL,
         )
 
         # 4. Switch tracking using existing pattern
diff --git a/flixopt/modeling.py b/flixopt/modeling.py
index 31511c90e..591702cd0 100644
--- a/flixopt/modeling.py
+++ b/flixopt/modeling.py
@@ -204,6 +204,7 @@ def expression_tracking_variable(
         short_name: str = None,
         bounds: tuple[xr.DataArray, xr.DataArray] = None,
         coords: str | list[str] | None = None,
+        category: VariableCategory = None,
     ) -> tuple[linopy.Variable, linopy.Constraint]:
         """Creates a variable constrained to equal a given expression.
 
@@ -218,6 +219,7 @@ def expression_tracking_variable(
             short_name: Short name for display purposes
             bounds: Optional (lower_bound, upper_bound) tuple for the tracker variable
             coords: Coordinate dimensions for the variable (None uses all model coords)
+            category: Category for segment expansion handling. See VariableCategory.
 
         Returns:
             Tuple of (tracker_variable, tracking_constraint)
@@ -226,7 +228,9 @@ def expression_tracking_variable(
             raise ValueError('ModelingPrimitives.expression_tracking_variable() can only be used with a Submodel')
 
         if not bounds:
-            tracker = model.add_variables(name=name, coords=model.get_coords(coords), short_name=short_name)
+            tracker = model.add_variables(
+                name=name, coords=model.get_coords(coords), short_name=short_name, category=category
+            )
         else:
             tracker = model.add_variables(
                 lower=bounds[0] if bounds[0] is not None else -np.inf,
@@ -234,6 +238,7 @@ def expression_tracking_variable(
                 name=name,
                 coords=model.get_coords(coords),
                 short_name=short_name,
+                category=category,
             )
 
         # Constraint: tracker = expression

From 9df0a8535772fcc4bc170d7619de16971faae031 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Thu, 15 Jan 2026 20:19:48 +0100
Subject: [PATCH 15/18]   Added to flow_system.py
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

  variable_categories property (line 1672):
  @property
  def variable_categories(self) -> dict[str, VariableCategory]:
      """Variable categories for filtering and segment expansion."""
      return self._variable_categories

  get_variables_by_category() method (line 1681):
  def get_variables_by_category(
      self, *categories: VariableCategory, from_solution: bool = True
  ) -> list[str]:
      """Get variable names matching any of the specified categories."""

  Updated in statistics_accessor.py
  ┌───────────────┬──────────────────────────────────────────┬──────────────────────────────────────────────────┐
  │   Property    │                  Before                  │                      After                       │
  ├───────────────┼──────────────────────────────────────────┼──────────────────────────────────────────────────┤
  │ flow_rates    │ endswith('|flow_rate')                   │ get_variables_by_category(FLOW_RATE)             │
  ├───────────────┼──────────────────────────────────────────┼──────────────────────────────────────────────────┤
  │ flow_sizes    │ endswith('|size') + flow_labels check    │ get_variables_by_category(SIZE) + flow_labels    │
  ├───────────────┼──────────────────────────────────────────┼──────────────────────────────────────────────────┤
  │ storage_sizes │ endswith('|size') + storage_labels check │ get_variables_by_category(SIZE) + storage_labels │
  ├───────────────┼──────────────────────────────────────────┼──────────────────────────────────────────────────┤
  │ charge_states │ endswith('|charge_state')                │ get_variables_by_category(CHARGE_STATE)          │
  └───────────────┴──────────────────────────────────────────┴──────────────────────────────────────────────────┘
  Benefits

  1. Single source of truth - Categories defined once in VariableCategory enum
  2. Easier maintenance - Adding new variable types only requires updating one place
  3. Type safety - Using enum values instead of magic strings
  4. Flexible filtering - Can filter by multiple categories: get_variables_by_category(SIZE, INVESTED)
  5. Consistent naming - Uses rsplit('|', 1)[0] instead of replace('|suffix', '') for label extraction
---
 flixopt/flow_system.py         | 33 +++++++++++++++++++++++++++++++++
 flixopt/statistics_accessor.py | 29 +++++++++++++----------------
 2 files changed, 46 insertions(+), 16 deletions(-)

diff --git a/flixopt/flow_system.py b/flixopt/flow_system.py
index 89466e7b6..9bf79292a 100644
--- a/flixopt/flow_system.py
+++ b/flixopt/flow_system.py
@@ -1669,6 +1669,39 @@ def solution(self, value: xr.Dataset | None) -> None:
         self._solution = value
         self._statistics = None  # Invalidate cached statistics
 
+    @property
+    def variable_categories(self) -> dict[str, VariableCategory]:
+        """Variable categories for filtering and segment expansion.
+
+        Returns:
+            Dict mapping variable names to their VariableCategory.
+        """
+        return self._variable_categories
+
+    def get_variables_by_category(self, *categories: VariableCategory, from_solution: bool = True) -> list[str]:
+        """Get variable names matching any of the specified categories.
+
+        Args:
+            *categories: One or more VariableCategory values to filter by.
+            from_solution: If True, only return variables present in solution.
+                If False, return all registered variables matching categories.
+
+        Returns:
+            List of variable names matching any of the specified categories.
+
+        Example:
+            >>> fs.get_variables_by_category(VariableCategory.FLOW_RATE)
+            ['Boiler(Q_th)|flow_rate', 'CHP(Q_th)|flow_rate', ...]
+            >>> fs.get_variables_by_category(VariableCategory.SIZE, VariableCategory.INVESTED)
+            ['Boiler(Q_th)|size', 'Boiler(Q_th)|invested', ...]
+        """
+        category_set = set(categories)
+        matching = [name for name, cat in self._variable_categories.items() if cat in category_set]
+        if from_solution and self._solution is not None:
+            solution_vars = set(self._solution.data_vars)
+            matching = [v for v in matching if v in solution_vars]
+        return matching
+
     @property
     def is_locked(self) -> bool:
         """Check if the FlowSystem is locked (has a solution).
diff --git a/flixopt/statistics_accessor.py b/flixopt/statistics_accessor.py
index 90ad875b7..632112d15 100644
--- a/flixopt/statistics_accessor.py
+++ b/flixopt/statistics_accessor.py
@@ -31,6 +31,7 @@
 from .color_processing import ColorType, hex_to_rgba, process_colors
 from .config import CONFIG
 from .plot_result import PlotResult
+from .structure import VariableCategory
 
 if TYPE_CHECKING:
     from .flow_system import FlowSystem
@@ -523,12 +524,12 @@ def flow_rates(self) -> xr.Dataset:
         """
         self._require_solution()
         if self._flow_rates is None:
-            flow_rate_vars = [v for v in self._fs.solution.data_vars if v.endswith('|flow_rate')]
+            flow_rate_vars = self._fs.get_variables_by_category(VariableCategory.FLOW_RATE)
             flow_carriers = self._fs.flow_carriers  # Cached lookup
             carrier_units = self.carrier_units  # Cached lookup
             data_vars = {}
             for v in flow_rate_vars:
-                flow_label = v.replace('|flow_rate', '')
+                flow_label = v.rsplit('|', 1)[0]  # Extract label from 'label|flow_rate'
                 da = self._fs.solution[v].copy()
                 # Add carrier and unit as attributes
                 carrier = flow_carriers.get(flow_label)
@@ -568,10 +569,10 @@ def flow_sizes(self) -> xr.Dataset:
         self._require_solution()
         if self._flow_sizes is None:
             flow_labels = set(self._fs.flows.keys())
-            size_vars = [
-                v for v in self._fs.solution.data_vars if v.endswith('|size') and v.replace('|size', '') in flow_labels
-            ]
-            self._flow_sizes = xr.Dataset({v.replace('|size', ''): self._fs.solution[v] for v in size_vars})
+            size_vars = self._fs.get_variables_by_category(VariableCategory.SIZE)
+            # Filter to only flow-related sizes
+            flow_size_vars = [v for v in size_vars if v.rsplit('|', 1)[0] in flow_labels]
+            self._flow_sizes = xr.Dataset({v.rsplit('|', 1)[0]: self._fs.solution[v] for v in flow_size_vars})
         return self._flow_sizes
 
     @property
@@ -580,12 +581,10 @@ def storage_sizes(self) -> xr.Dataset:
         self._require_solution()
         if self._storage_sizes is None:
             storage_labels = set(self._fs.storages.keys())
-            size_vars = [
-                v
-                for v in self._fs.solution.data_vars
-                if v.endswith('|size') and v.replace('|size', '') in storage_labels
-            ]
-            self._storage_sizes = xr.Dataset({v.replace('|size', ''): self._fs.solution[v] for v in size_vars})
+            size_vars = self._fs.get_variables_by_category(VariableCategory.SIZE)
+            # Filter to only storage-related sizes
+            storage_size_vars = [v for v in size_vars if v.rsplit('|', 1)[0] in storage_labels]
+            self._storage_sizes = xr.Dataset({v.rsplit('|', 1)[0]: self._fs.solution[v] for v in storage_size_vars})
         return self._storage_sizes
 
     @property
@@ -600,10 +599,8 @@ def charge_states(self) -> xr.Dataset:
         """All storage charge states as a Dataset with storage labels as variable names."""
         self._require_solution()
         if self._charge_states is None:
-            charge_vars = [v for v in self._fs.solution.data_vars if v.endswith('|charge_state')]
-            self._charge_states = xr.Dataset(
-                {v.replace('|charge_state', ''): self._fs.solution[v] for v in charge_vars}
-            )
+            charge_vars = self._fs.get_variables_by_category(VariableCategory.CHARGE_STATE)
+            self._charge_states = xr.Dataset({v.rsplit('|', 1)[0]: self._fs.solution[v] for v in charge_vars})
         return self._charge_states
 
     @property

From cf59335f1b31c6ae0651b9b0a7a14948f429f7aa Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Thu, 15 Jan 2026 20:21:09 +0100
Subject: [PATCH 16/18] Ensure backwards compatability

---
 flixopt/flow_system.py | 13 ++++++++++++-
 1 file changed, 12 insertions(+), 1 deletion(-)

diff --git a/flixopt/flow_system.py b/flixopt/flow_system.py
index 9bf79292a..98b4357e0 100644
--- a/flixopt/flow_system.py
+++ b/flixopt/flow_system.py
@@ -1696,7 +1696,18 @@ def get_variables_by_category(self, *categories: VariableCategory, from_solution
             ['Boiler(Q_th)|size', 'Boiler(Q_th)|invested', ...]
         """
         category_set = set(categories)
-        matching = [name for name, cat in self._variable_categories.items() if cat in category_set]
+
+        if self._variable_categories:
+            # Use registered categories
+            matching = [name for name, cat in self._variable_categories.items() if cat in category_set]
+        elif self._solution is not None:
+            # Fallback for old files without categories: match by suffix pattern
+            # Category values match the variable suffix (e.g., FLOW_RATE.value = 'flow_rate')
+            suffixes = tuple(f'|{cat.value}' for cat in category_set)
+            matching = [v for v in self._solution.data_vars if v.endswith(suffixes)]
+        else:
+            matching = []
+
         if from_solution and self._solution is not None:
             solution_vars = set(self._solution.data_vars)
             matching = [v for v in matching if v in solution_vars]

From 7acec0c2a6281ccf225fd04a044b730794037004 Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Thu, 15 Jan 2026 20:43:34 +0100
Subject: [PATCH 17/18]   Summary of Changes
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

  1. New SIZE Sub-Categories (structure.py)

  - Added FLOW_SIZE and STORAGE_SIZE to differentiate flow vs storage investments
  - Kept SIZE for backward compatibility

  2. InvestmentModel Updated (features.py)

  - Added size_category parameter to InvestmentModel.__init__()
  - Callers now specify the appropriate category

  3. Variable Registrations Updated

  - elements.py: FlowModel uses FLOW_SIZE
  - components.py: StorageModel uses STORAGE_SIZE (2 locations)

  4. Statistics Accessor Simplified (statistics_accessor.py)

  - flow_sizes: Now uses get_variables_by_category(FLOW_SIZE) directly
  - storage_sizes: Now uses get_variables_by_category(STORAGE_SIZE) directly
  - No more filtering by element labels after getting SIZE variables

  5. Backward-Compatible Fallback (flow_system.py)

  - get_variables_by_category() handles old files:
    - FLOW_SIZE → matches |size suffix + flow labels
    - STORAGE_SIZE → matches |size suffix + storage labels

  6. SOC Boundary Pattern Matching Replaced (transform_accessor.py)

  - Changed from endswith('|SOC_boundary') to get_variables_by_category(SOC_BOUNDARY)

  7. Effect Variables Verified

  - PER_TIMESTEP ✓ (features.py:659)
  - SHARE ✓ (features.py:700 for temporal shares)
  - TOTAL / TOTAL_OVER_PERIODS ✓ (multiple locations)

  8. Documentation Updated

  - _build_segment_total_varnames() marked as backwards-compatibility fallback

  Benefits

  - Cleaner code: No more string manipulation to filter by element type
  - Type safety: Using enum values instead of magic strings
  - Single source of truth: Categories defined once, used everywhere
  - Backward compatible: Old files still work via fallback logic
---
 flixopt/components.py          |  2 ++
 flixopt/elements.py            |  1 +
 flixopt/features.py            |  5 ++++-
 flixopt/flow_system.py         | 23 +++++++++++++++++++++--
 flixopt/statistics_accessor.py | 10 ++--------
 flixopt/structure.py           |  4 +++-
 flixopt/transform_accessor.py  | 13 +++++++------
 7 files changed, 40 insertions(+), 18 deletions(-)

diff --git a/flixopt/components.py b/flixopt/components.py
index 56a0f5632..a3bd20014 100644
--- a/flixopt/components.py
+++ b/flixopt/components.py
@@ -977,6 +977,7 @@ def _add_investment_model(self):
                     label_of_element=self.label_of_element,
                     label_of_model=self.label_of_element,
                     parameters=self.element.capacity_in_flow_hours,
+                    size_category=VariableCategory.STORAGE_SIZE,
                 ),
                 short_name='investment',
             )
@@ -1291,6 +1292,7 @@ def _add_investment_model(self):
                     label_of_element=self.label_of_element,
                     label_of_model=self.label_of_element,
                     parameters=self.element.capacity_in_flow_hours,
+                    size_category=VariableCategory.STORAGE_SIZE,
                 ),
                 short_name='investment',
             )
diff --git a/flixopt/elements.py b/flixopt/elements.py
index c80c263d1..e2def702d 100644
--- a/flixopt/elements.py
+++ b/flixopt/elements.py
@@ -755,6 +755,7 @@ def _create_investment_model(self):
                 label_of_element=self.label_of_element,
                 parameters=self.element.size,
                 label_of_model=self.label_of_element,
+                size_category=VariableCategory.FLOW_SIZE,
             ),
             'investment',
         )
diff --git a/flixopt/features.py b/flixopt/features.py
index 9ef400392..e85636435 100644
--- a/flixopt/features.py
+++ b/flixopt/features.py
@@ -37,6 +37,7 @@ class InvestmentModel(Submodel):
         label_of_element: The label of the parent (Element). Used to construct the full label of the model.
         parameters: The parameters of the feature model.
         label_of_model: The label of the model. This is needed to construct the full label of the model.
+        size_category: Category for the size variable (FLOW_SIZE, STORAGE_SIZE, or SIZE for generic).
     """
 
     parameters: InvestParameters
@@ -47,9 +48,11 @@ def __init__(
         label_of_element: str,
         parameters: InvestParameters,
         label_of_model: str | None = None,
+        size_category: VariableCategory = VariableCategory.SIZE,
     ):
         self.piecewise_effects: PiecewiseEffectsModel | None = None
         self.parameters = parameters
+        self._size_category = size_category
         super().__init__(model, label_of_element=label_of_element, label_of_model=label_of_model)
 
     def _do_modeling(self):
@@ -69,7 +72,7 @@ def _create_variables_and_constraints(self):
             lower=size_min if self.parameters.mandatory else 0,
             upper=size_max,
             coords=self._model.get_coords(['period', 'scenario']),
-            category=VariableCategory.SIZE,
+            category=self._size_category,
         )
 
         if not self.parameters.mandatory:
diff --git a/flixopt/flow_system.py b/flixopt/flow_system.py
index 98b4357e0..e32d1b2ad 100644
--- a/flixopt/flow_system.py
+++ b/flixopt/flow_system.py
@@ -1703,8 +1703,27 @@ def get_variables_by_category(self, *categories: VariableCategory, from_solution
         elif self._solution is not None:
             # Fallback for old files without categories: match by suffix pattern
             # Category values match the variable suffix (e.g., FLOW_RATE.value = 'flow_rate')
-            suffixes = tuple(f'|{cat.value}' for cat in category_set)
-            matching = [v for v in self._solution.data_vars if v.endswith(suffixes)]
+            matching = []
+            for cat in category_set:
+                # Handle new sub-categories that map to old |size suffix
+                if cat == VariableCategory.FLOW_SIZE:
+                    flow_labels = set(self.flows.keys())
+                    matching.extend(
+                        v
+                        for v in self._solution.data_vars
+                        if v.endswith('|size') and v.rsplit('|', 1)[0] in flow_labels
+                    )
+                elif cat == VariableCategory.STORAGE_SIZE:
+                    storage_labels = set(self.storages.keys())
+                    matching.extend(
+                        v
+                        for v in self._solution.data_vars
+                        if v.endswith('|size') and v.rsplit('|', 1)[0] in storage_labels
+                    )
+                else:
+                    # Standard suffix matching
+                    suffix = f'|{cat.value}'
+                    matching.extend(v for v in self._solution.data_vars if v.endswith(suffix))
         else:
             matching = []
 
diff --git a/flixopt/statistics_accessor.py b/flixopt/statistics_accessor.py
index 632112d15..0092d4989 100644
--- a/flixopt/statistics_accessor.py
+++ b/flixopt/statistics_accessor.py
@@ -568,10 +568,7 @@ def flow_sizes(self) -> xr.Dataset:
         """Flow sizes as a Dataset with flow labels as variable names."""
         self._require_solution()
         if self._flow_sizes is None:
-            flow_labels = set(self._fs.flows.keys())
-            size_vars = self._fs.get_variables_by_category(VariableCategory.SIZE)
-            # Filter to only flow-related sizes
-            flow_size_vars = [v for v in size_vars if v.rsplit('|', 1)[0] in flow_labels]
+            flow_size_vars = self._fs.get_variables_by_category(VariableCategory.FLOW_SIZE)
             self._flow_sizes = xr.Dataset({v.rsplit('|', 1)[0]: self._fs.solution[v] for v in flow_size_vars})
         return self._flow_sizes
 
@@ -580,10 +577,7 @@ def storage_sizes(self) -> xr.Dataset:
         """Storage capacity sizes as a Dataset with storage labels as variable names."""
         self._require_solution()
         if self._storage_sizes is None:
-            storage_labels = set(self._fs.storages.keys())
-            size_vars = self._fs.get_variables_by_category(VariableCategory.SIZE)
-            # Filter to only storage-related sizes
-            storage_size_vars = [v for v in size_vars if v.rsplit('|', 1)[0] in storage_labels]
+            storage_size_vars = self._fs.get_variables_by_category(VariableCategory.STORAGE_SIZE)
             self._storage_sizes = xr.Dataset({v.rsplit('|', 1)[0]: self._fs.solution[v] for v in storage_size_vars})
         return self._storage_sizes
 
diff --git a/flixopt/structure.py b/flixopt/structure.py
index 7d76c760c..6dea3aa19 100644
--- a/flixopt/structure.py
+++ b/flixopt/structure.py
@@ -73,7 +73,9 @@ class VariableCategory(Enum):
     TOTAL_OVER_PERIODS = 'total_over_periods'  # Effect total over all periods
 
     # === Investment ===
-    SIZE = 'size'  # Investment size
+    SIZE = 'size'  # Generic investment size (for backwards compatibility)
+    FLOW_SIZE = 'flow_size'  # Flow investment size
+    STORAGE_SIZE = 'storage_size'  # Storage capacity size
     INVESTED = 'invested'  # Invested yes/no binary
 
     # === Counting/Duration ===
diff --git a/flixopt/transform_accessor.py b/flixopt/transform_accessor.py
index 3f7841085..a81de0aa9 100644
--- a/flixopt/transform_accessor.py
+++ b/flixopt/transform_accessor.py
@@ -16,7 +16,7 @@
 import pandas as pd
 import xarray as xr
 
-from .structure import EXPAND_DIVIDE, EXPAND_INTERPOLATE
+from .structure import EXPAND_DIVIDE, EXPAND_INTERPOLATE, VariableCategory
 
 if TYPE_CHECKING:
     from tsam.config import ClusterConfig, ExtremeConfig, SegmentConfig
@@ -1697,7 +1697,7 @@ def _combine_intercluster_charge_states(
             n_original_clusters: Number of original clusters before aggregation.
         """
         n_original_timesteps_extra = len(original_timesteps_extra)
-        soc_boundary_vars = [name for name in reduced_solution.data_vars if name.endswith('|SOC_boundary')]
+        soc_boundary_vars = self._fs.get_variables_by_category(VariableCategory.SOC_BOUNDARY)
 
         for soc_boundary_name in soc_boundary_vars:
             storage_name = soc_boundary_name.rsplit('|', 1)[0]
@@ -1800,15 +1800,16 @@ def _apply_soc_decay(
         return soc_boundary_per_timestep * decay_da
 
     def _build_segment_total_varnames(self) -> set[str]:
-        """Build the set of solution variable names that represent segment totals.
+        """Build segment total variable names - BACKWARDS COMPATIBILITY FALLBACK.
+
+        This method is only used when variable_categories is empty (old FlowSystems
+        saved before category registration was implemented). New FlowSystems use
+        the VariableCategory registry with EXPAND_DIVIDE categories (PER_TIMESTEP, SHARE).
 
         For segmented systems, these variables contain values that are summed over
         segments. When expanded to hourly resolution, they need to be divided by
         segment duration to get correct hourly rates.
 
-        Derives variable names directly from FlowSystem structure (effects, flows,
-        components) rather than pattern matching, ensuring robustness.
-
         Returns:
             Set of variable names that should be divided by expansion divisor.
         """

From d3ce0647b7f3fdc7705fb422c7afda6c1d8391fa Mon Sep 17 00:00:00 2001
From: FBumann <117816358+FBumann@users.noreply.github.com>
Date: Fri, 16 Jan 2026 12:28:46 +0100
Subject: [PATCH 18/18] perf: Keep data in minimal form (no pre-broadcasting)
 (#575)
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

* Add collapse to IO

* Temp

* 1. Added dimension reduction functions in io.py:
    - _reduce_constant_dims(): Reduces dimensions where all values are constant
    - _expand_reduced_dims(): Restores original dimensions with correct order
  2. Added reduce_constants parameter to FlowSystem.to_dataset():
    - Default: False (opt-in)
    - When True: Reduces constant dimensions for memory efficiency
  3. Updated FlowSystem.from_dataset():
    - Expands both collapsed arrays (from NetCDF) and reduced dimensions
  4. Key fixes:
    - Store original dimension order in attrs to preserve (cluster, time) vs (time, cluster) ordering
    - Skip solution variables in reduction (they're prefixed with solution|)

  The optimization is opt-in via to_dataset(reduce_constants=True). For file storage, save_dataset_to_netcdf still collapses constants to scalars by default.

* Set reduce_constants=True by default in to_dataset()

Benchmarks show 99%+ reduction in memory and file size for multi-period
models, with faster IO speeds. The optimization is now enabled by default.

Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>

* Feature/speedup io+internals (#576)

* Add safe isel wrappers for dimension-independent operations

- Add _scalar_safe_isel_drop() to modeling.py for selecting from
  potentially reduced arrays (handles both scalar and array cases)
- Update Storage validation to use _scalar_safe_isel for
  relative_minimum/maximum_charge_state access at time=0
- Update StorageModel._relative_charge_state_bounds to handle
  reduced arrays properly with dimension checks

These wrappers prepare the codebase for future optimization where
_expand_reduced_dims() could be removed from from_dataset(),
keeping data compact throughout the system.

Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>

* Remove pre-broadcasting from DataConverter - data stays minimal

Major change: DataConverter.to_dataarray() now validates dimensions
instead of broadcasting to full target dimensions. This keeps data
compact (scalars stay scalar, 1D arrays stay 1D) and lets linopy
handle broadcasting at variable creation time.

Changes:
- core.py: Replace _broadcast_dataarray_to_target_specification with
  _validate_dataarray_dims in to_dataarray() method
- components.py: Fix _relative_charge_state_bounds to handle scalar
  inputs that have no time dimension (expand to timesteps_extra)
- conftest.py: Add assert_dims_compatible helper for tests
- test_*.py: Update dimension assertions to use subset checking

Memory impact (8760 timesteps × 20 periods):
- Before: 24.1 MB dataset size
- After: 137.1 KB dataset size

Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>

* Remove redundant dimension reduction - data is minimal from start

Since DataConverter no longer broadcasts, these optimizations are redundant:
- Remove _reduce_constant_dims() and _expand_reduced_dims() from io.py
- Remove reduce_constants parameter from to_dataset()
- Remove _expand_reduced_dims() call from from_dataset()

Also ensure Dataset always has model coordinates (time, period, scenario,
cluster) even if no data variable uses them. This is important for Dataset
validity and downstream operations.

Update benchmark_io.py to reflect the simplified API.

Memory footprint now:
- 24 timesteps: 568 bytes (23/24 vars are scalar)
- 168 timesteps: 2.8 KB
- 8760 timesteps × 20 periods: ~137 KB (vs 24 MB with old broadcasting)

Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>

* ⏺ Done. I've removed the _collapse_constant_arrays / _expand_collapsed_arrays optimization from the codebase. Here's a summary of the changes:

  Files modified:

  1. flixopt/io.py:
    - Removed COLLAPSED_VAR_PREFIX constant
    - Removed _collapse_constant_arrays() function (~45 lines)
    - Removed _expand_collapsed_arrays() function (~40 lines)
    - Removed collapse_constants parameter from save_dataset_to_netcdf()
    - Removed expansion logic from load_dataset_from_netcdf()
  2. flixopt/flow_system.py:
    - Removed import and call to _expand_collapsed_arrays in from_dataset()
  3. benchmark_io.py:
    - Simplified benchmark_netcdf() to only benchmark single save/load (no more comparison)
    - Removed collapse_constants parameter from roundtrip function

  Rationale: Since data is now kept in minimal form throughout the system (scalars stay scalars, 1D arrays stay 1D), there's no need for the extra collapsing/expanding layer when saving to NetCDF. The file sizes are naturally small because the data isn't expanded to full dimensions.

  Test results:
  - All 10 IO tests pass
  - All 4 integration tests pass
  - Benchmark runs successfully

* Summary of Fixes

  1. _dataset_resample handling all-scalar data (transform_accessor.py)

  When all data variables are scalars (no time dimension), the resample function now:
  - Creates a dummy variable to resample the time coordinate
  - Preserves all original scalar data variables
  - Preserves all non-time coordinates (period, scenario, cluster)

  2. _scalar_safe_reduce helper (modeling.py)

  Added a new helper function that safely performs aggregation operations (mean/sum/etc) over a dimension:
  - Returns reduced data if the dimension exists
  - Returns data unchanged if the dimension doesn't exist (scalar data)

  3. Updated Storage intercluster linking (components.py)

  Used _scalar_safe_reduce for:
  - relative_loss_per_hour.mean('time')
  - timestep_duration.mean('time')

  4. Updated clustering expansion (transform_accessor.py)

  Used _scalar_safe_reduce for relative_loss_per_hour.mean('time')

  5. Fixed dimension order in expand_data (clustering/base.py)

  When expanding clustered data back to original timesteps:
  - Added transpose to ensure (cluster, time) dimension order before numpy indexing
  - Fixes IndexError when dimensions were in different order

* Added reduction on loading old datasets

* Summary of changes made:

  1. flixopt/structure.py: Added dimension transpose in _ensure_coords() to ensure correct dimension order when data already has all dims but in wrong order
  2. tests/test_io_conversion.py: Updated test_returns_same_object → test_returns_equivalent_dataset since convert_old_dataset now creates a new dataset when reducing constants
  3. tests/deprecated/test_flow.py: Updated 6 assertions to expect minimal dimensions ('time',) instead of broadcasting to all model coords
  4. tests/deprecated/test_component.py: Updated 2 assertions to expect minimal dimensions ('time',)
  5. tests/deprecated/test_linear_converter.py: Updated 1 assertion to expect minimal dimensions ('time',)

  The key change is that data now stays in minimal form - a 1D time-varying array stays as ('time',) dimensions rather than being pre-broadcast to ('time', 'scenario') or other full model dimensions. Broadcasting happens at the linopy interface layer in FlowSystemModel.add_variables() when needed.

---------

Co-authored-by: Claude Opus 4.5 <noreply@anthropic.com>

---------

Co-authored-by: Claude Opus 4.5 <noreply@anthropic.com>
---
 flixopt/clustering/base.py                |    3 +
 flixopt/components.py                     |   61 +-
 flixopt/core.py                           |   54 +-
 flixopt/flow_system.py                    |   15 +
 flixopt/io.py                             |   67 +-
 flixopt/modeling.py                       |   70 +-
 flixopt/structure.py                      |   56 +
 flixopt/transform_accessor.py             |   30 +-
 tests/conftest.py                         |   19 +
 tests/deprecated/test_component.py        |    6 +-
 tests/deprecated/test_dataconverter.py    | 1262 ---------------------
 tests/deprecated/test_flow.py             |   30 +-
 tests/deprecated/test_linear_converter.py |    3 +-
 tests/test_component.py                   |    5 +-
 tests/test_dataconverter.py               |  366 +++---
 tests/test_flow.py                        |   26 +-
 tests/test_io_conversion.py               |    8 +-
 tests/test_linear_converter.py            |    4 +-
 18 files changed, 554 insertions(+), 1531 deletions(-)
 delete mode 100644 tests/deprecated/test_dataconverter.py

diff --git a/flixopt/clustering/base.py b/flixopt/clustering/base.py
index 3e27697e6..10224a1a5 100644
--- a/flixopt/clustering/base.py
+++ b/flixopt/clustering/base.py
@@ -439,6 +439,9 @@ def _expand_slice(mapping: np.ndarray, data: xr.DataArray) -> np.ndarray:
             if has_cluster_dim:
                 cluster_ids = mapping // timesteps_per_cluster
                 time_within = mapping % timesteps_per_cluster
+                # Ensure dimension order is (cluster, time) for correct indexing
+                if data.dims != ('cluster', 'time'):
+                    data = data.transpose('cluster', 'time')
                 return data.values[cluster_ids, time_within]
             return data.values[mapping]
 
diff --git a/flixopt/components.py b/flixopt/components.py
index b720dd0ba..7cb5b9fc4 100644
--- a/flixopt/components.py
+++ b/flixopt/components.py
@@ -16,7 +16,7 @@
 from .elements import Component, ComponentModel, Flow
 from .features import InvestmentModel, PiecewiseModel
 from .interface import InvestParameters, PiecewiseConversion, StatusParameters
-from .modeling import BoundingPatterns
+from .modeling import BoundingPatterns, _scalar_safe_isel, _scalar_safe_isel_drop, _scalar_safe_reduce
 from .structure import FlowSystemModel, register_class_for_io
 
 if TYPE_CHECKING:
@@ -570,8 +570,12 @@ def _plausibility_checks(self) -> None:
             # Initial charge state should not constrain investment decision
             # If initial > (min_cap * rel_max), investment is forced to increase capacity
             # If initial < (max_cap * rel_min), investment is forced to decrease capacity
-            min_initial_at_max_capacity = maximum_capacity * self.relative_minimum_charge_state.isel(time=0)
-            max_initial_at_min_capacity = minimum_capacity * self.relative_maximum_charge_state.isel(time=0)
+            min_initial_at_max_capacity = maximum_capacity * _scalar_safe_isel(
+                self.relative_minimum_charge_state, {'time': 0}
+            )
+            max_initial_at_min_capacity = minimum_capacity * _scalar_safe_isel(
+                self.relative_maximum_charge_state, {'time': 0}
+            )
 
             # Only perform numeric comparisons if using a numeric initial_charge_state
             if not initial_equals_final and self.initial_charge_state is not None:
@@ -1100,24 +1104,47 @@ def _relative_charge_state_bounds(self) -> tuple[xr.DataArray, xr.DataArray]:
         Returns:
             Tuple of (minimum_bounds, maximum_bounds) DataArrays extending to final timestep
         """
-        final_coords = {'time': [self._model.flow_system.timesteps_extra[-1]]}
+        timesteps_extra = self._model.flow_system.timesteps_extra
+
+        # Get the original bounds (may be scalar or have time dim)
+        rel_min = self.element.relative_minimum_charge_state
+        rel_max = self.element.relative_maximum_charge_state
 
         # Get final minimum charge state
         if self.element.relative_minimum_final_charge_state is None:
-            min_final = self.element.relative_minimum_charge_state.isel(time=-1, drop=True)
+            min_final_value = _scalar_safe_isel_drop(rel_min, 'time', -1)
         else:
-            min_final = self.element.relative_minimum_final_charge_state
-        min_final = min_final.expand_dims('time').assign_coords(time=final_coords['time'])
+            min_final_value = self.element.relative_minimum_final_charge_state
 
         # Get final maximum charge state
         if self.element.relative_maximum_final_charge_state is None:
-            max_final = self.element.relative_maximum_charge_state.isel(time=-1, drop=True)
+            max_final_value = _scalar_safe_isel_drop(rel_max, 'time', -1)
+        else:
+            max_final_value = self.element.relative_maximum_final_charge_state
+
+        # Build bounds arrays for timesteps_extra (includes final timestep)
+        # Handle case where original data may be scalar (no time dim)
+        if 'time' in rel_min.dims:
+            # Original has time dim - concat with final value
+            min_final_da = (
+                min_final_value.expand_dims('time') if 'time' not in min_final_value.dims else min_final_value
+            )
+            min_final_da = min_final_da.assign_coords(time=[timesteps_extra[-1]])
+            min_bounds = xr.concat([rel_min, min_final_da], dim='time')
+        else:
+            # Original is scalar - broadcast to full time range (constant value)
+            min_bounds = rel_min.expand_dims(time=timesteps_extra)
+
+        if 'time' in rel_max.dims:
+            # Original has time dim - concat with final value
+            max_final_da = (
+                max_final_value.expand_dims('time') if 'time' not in max_final_value.dims else max_final_value
+            )
+            max_final_da = max_final_da.assign_coords(time=[timesteps_extra[-1]])
+            max_bounds = xr.concat([rel_max, max_final_da], dim='time')
         else:
-            max_final = self.element.relative_maximum_final_charge_state
-        max_final = max_final.expand_dims('time').assign_coords(time=final_coords['time'])
-        # Concatenate with original bounds
-        min_bounds = xr.concat([self.element.relative_minimum_charge_state, min_final], dim='time')
-        max_bounds = xr.concat([self.element.relative_maximum_charge_state, max_final], dim='time')
+            # Original is scalar - broadcast to full time range (constant value)
+            max_bounds = rel_max.expand_dims(time=timesteps_extra)
 
         return min_bounds, max_bounds
 
@@ -1477,8 +1504,8 @@ def _add_linking_constraints(
         # relative_loss_per_hour is per-hour, so we need hours = timesteps * duration
         # Use mean over time (linking operates at period level, not timestep)
         # Keep as DataArray to respect per-period/scenario values
-        rel_loss = self.element.relative_loss_per_hour.mean('time')
-        hours_per_cluster = timesteps_per_cluster * self._model.timestep_duration.mean('time')
+        rel_loss = _scalar_safe_reduce(self.element.relative_loss_per_hour, 'time', 'mean')
+        hours_per_cluster = timesteps_per_cluster * _scalar_safe_reduce(self._model.timestep_duration, 'time', 'mean')
         decay_n = (1 - rel_loss) ** hours_per_cluster
 
         lhs = soc_after - soc_before * decay_n - delta_soc_ordered
@@ -1522,8 +1549,8 @@ def _add_combined_bound_constraints(
         # Get self-discharge rate for decay calculation
         # relative_loss_per_hour is per-hour, so we need to convert offsets to hours
         # Keep as DataArray to respect per-period/scenario values
-        rel_loss = self.element.relative_loss_per_hour.mean('time')
-        mean_timestep_duration = self._model.timestep_duration.mean('time')
+        rel_loss = _scalar_safe_reduce(self.element.relative_loss_per_hour, 'time', 'mean')
+        mean_timestep_duration = _scalar_safe_reduce(self._model.timestep_duration, 'time', 'mean')
 
         sample_offsets = [0, timesteps_per_cluster // 2, timesteps_per_cluster - 1]
 
diff --git a/flixopt/core.py b/flixopt/core.py
index 3d456fff1..02cfc2940 100644
--- a/flixopt/core.py
+++ b/flixopt/core.py
@@ -366,6 +366,55 @@ def _broadcast_dataarray_to_target_specification(
         broadcasted = source_data.broadcast_like(target_template)
         return broadcasted.transpose(*target_dims)
 
+    @staticmethod
+    def _validate_dataarray_dims(
+        data: xr.DataArray, target_coords: dict[str, pd.Index], target_dims: tuple[str, ...]
+    ) -> xr.DataArray:
+        """
+        Validate that DataArray dims are a subset of target dims (without broadcasting).
+
+        This method validates compatibility without expanding to full dimensions,
+        allowing data to remain in compact form. Broadcasting happens later at
+        the linopy interface (FlowSystemModel.add_variables).
+
+        Also transposes data to canonical dimension order (matching target_dims order).
+
+        Args:
+            data: DataArray to validate
+            target_coords: Target coordinates {dim_name: coordinate_index}
+            target_dims: Target dimension names in canonical order
+
+        Returns:
+            DataArray with validated dims, transposed to canonical order
+
+        Raises:
+            ConversionError: If data has dimensions not in target_dims,
+                           or coordinate values don't match
+        """
+        # Validate: all data dimensions must exist in target
+        extra_dims = set(data.dims) - set(target_dims)
+        if extra_dims:
+            raise ConversionError(f'Data has dimensions {extra_dims} not in target dimensions {target_dims}')
+
+        # Validate: coordinate compatibility for overlapping dimensions
+        for dim in data.dims:
+            if dim in data.coords and dim in target_coords:
+                data_coords = data.coords[dim]
+                target_coords_for_dim = target_coords[dim]
+
+                if not np.array_equal(data_coords.values, target_coords_for_dim.values):
+                    raise ConversionError(
+                        f'Coordinate mismatch for dimension "{dim}". Data and target coordinates have different values.'
+                    )
+
+        # Transpose to canonical dimension order (subset of target_dims that data has)
+        if data.dims:
+            canonical_order = tuple(d for d in target_dims if d in data.dims)
+            if data.dims != canonical_order:
+                data = data.transpose(*canonical_order)
+
+        return data
+
     @classmethod
     def to_dataarray(
         cls,
@@ -480,8 +529,9 @@ def to_dataarray(
                 f'Unsupported data type: {type(data).__name__}. Supported types: {", ".join(supported_types)}'
             )
 
-        # Broadcast intermediate result to target specification
-        return cls._broadcast_dataarray_to_target_specification(intermediate, validated_coords, target_dims)
+        # Validate dims are compatible (no broadcasting - data stays compact)
+        # Broadcasting happens at FlowSystemModel.add_variables() via _ensure_coords
+        return cls._validate_dataarray_dims(intermediate, validated_coords, target_dims)
 
     @staticmethod
     def _validate_and_prepare_target_coordinates(
diff --git a/flixopt/flow_system.py b/flixopt/flow_system.py
index 775789963..2a9d2c383 100644
--- a/flixopt/flow_system.py
+++ b/flixopt/flow_system.py
@@ -632,6 +632,10 @@ def to_dataset(self, include_solution: bool = True) -> xr.Dataset:
         Convert the FlowSystem to an xarray Dataset.
         Ensures FlowSystem is connected before serialization.
 
+        Data is stored in minimal form (scalars stay scalar, 1D arrays stay 1D) without
+        broadcasting to full model dimensions. This provides significant memory savings
+        for multi-period and multi-scenario models.
+
         If a solution is present and `include_solution=True`, it will be included
         in the dataset with variable names prefixed by 'solution|' to avoid conflicts
         with FlowSystem configuration variables. Solution time coordinates are renamed
@@ -687,6 +691,17 @@ def to_dataset(self, include_solution: bool = True) -> xr.Dataset:
         # Add version info
         ds.attrs['flixopt_version'] = __version__
 
+        # Ensure model coordinates are always present in the Dataset
+        # (even if no data variable uses them, they define the model structure)
+        model_coords = {'time': self.timesteps}
+        if self.periods is not None:
+            model_coords['period'] = self.periods
+        if self.scenarios is not None:
+            model_coords['scenario'] = self.scenarios
+        if self.clusters is not None:
+            model_coords['cluster'] = self.clusters
+        ds = ds.assign_coords(model_coords)
+
         return ds
 
     @classmethod
diff --git a/flixopt/io.py b/flixopt/io.py
index 164a9fb2e..22347ec08 100644
--- a/flixopt/io.py
+++ b/flixopt/io.py
@@ -548,6 +548,7 @@ def save_dataset_to_netcdf(
         raise ValueError(f'Invalid file extension for path {path}. Only .nc and .nc4 are supported')
 
     ds = ds.copy(deep=True)
+
     ds.attrs = {'attrs': json.dumps(ds.attrs)}
 
     # Convert all DataArray attrs to JSON strings
@@ -572,6 +573,49 @@ def save_dataset_to_netcdf(
         )
 
 
+def _reduce_constant_arrays(ds: xr.Dataset) -> xr.Dataset:
+    """
+    Reduce constant dimensions in arrays for more efficient storage.
+
+    For each array, checks each dimension and removes it if values are constant
+    along that dimension. This handles cases like:
+    - Shape (8760,) all identical → scalar
+    - Shape (8760, 2) constant along time → shape (2,)
+    - Shape (8760, 2, 3) constant along time → shape (2, 3)
+
+    This is useful for datasets saved with older versions where data was
+    broadcast to full dimensions.
+
+    Args:
+        ds: Dataset with potentially constant arrays.
+
+    Returns:
+        Dataset with constant dimensions reduced.
+    """
+    new_data_vars = {}
+
+    for name, da in ds.data_vars.items():
+        if not da.dims or da.size == 0:
+            new_data_vars[name] = da
+            continue
+
+        # Try to reduce each dimension
+        reduced = da
+        for dim in list(da.dims):
+            if dim not in reduced.dims:
+                continue  # Already removed
+            # Check if constant along this dimension
+            first_slice = reduced.isel({dim: 0})
+            is_constant = (reduced == first_slice).all()
+            if is_constant:
+                # Remove this dimension by taking first slice
+                reduced = first_slice
+
+        new_data_vars[name] = reduced
+
+    return xr.Dataset(new_data_vars, coords=ds.coords, attrs=ds.attrs)
+
+
 def load_dataset_from_netcdf(path: str | pathlib.Path) -> xr.Dataset:
     """
     Load a dataset from a netcdf file. Load all attrs from 'attrs' attributes.
@@ -741,20 +785,25 @@ def convert_old_dataset(
     ds: xr.Dataset,
     key_renames: dict[str, str] | None = None,
     value_renames: dict[str, dict] | None = None,
+    reduce_constants: bool = True,
 ) -> xr.Dataset:
-    """Convert an old FlowSystem dataset to use new parameter names.
+    """Convert an old FlowSystem dataset to the current format.
+
+    This function performs two conversions:
+    1. Renames parameters in the reference structure to current naming conventions
+    2. Reduces constant arrays to minimal dimensions (e.g., broadcasted scalars back to scalars)
 
-    This function updates the reference structure in a dataset's attrs to use
-    the current parameter naming conventions. This is useful for loading
-    FlowSystem files saved with older versions of flixopt.
+    This is useful for loading FlowSystem files saved with older versions of flixopt.
 
     Args:
-        ds: The dataset to convert (will be modified in place)
+        ds: The dataset to convert
         key_renames: Custom key renames to apply. If None, uses PARAMETER_RENAMES.
         value_renames: Custom value renames to apply. If None, uses VALUE_RENAMES.
+        reduce_constants: If True (default), reduce constant arrays to minimal dimensions.
+            Old files may have scalars broadcasted to full (time, period, scenario) shape.
 
     Returns:
-        The converted dataset (same object, modified in place)
+        The converted dataset
 
     Examples:
         Convert an old netCDF file to new format:
@@ -765,7 +814,7 @@ def convert_old_dataset(
         # Load old file
         ds = io.load_dataset_from_netcdf('old_flow_system.nc4')
 
-        # Convert parameter names
+        # Convert to current format
         ds = io.convert_old_dataset(ds)
 
         # Now load as FlowSystem
@@ -782,6 +831,10 @@ def convert_old_dataset(
     # Convert the attrs (reference_structure)
     ds.attrs = _rename_keys_recursive(ds.attrs, key_renames, value_renames)
 
+    # Reduce constant arrays to minimal dimensions
+    if reduce_constants:
+        ds = _reduce_constant_arrays(ds)
+
     return ds
 
 
diff --git a/flixopt/modeling.py b/flixopt/modeling.py
index f22ffbc04..a0abeec77 100644
--- a/flixopt/modeling.py
+++ b/flixopt/modeling.py
@@ -1,4 +1,5 @@
 import logging
+from typing import Any
 
 import linopy
 import numpy as np
@@ -10,6 +11,71 @@
 logger = logging.getLogger('flixopt')
 
 
+def _scalar_safe_isel(data: xr.DataArray | Any, indexers: dict) -> xr.DataArray | Any:
+    """Apply isel if data has the required dimensions, otherwise return data as-is.
+
+    This allows parameters to remain compact (scalar or lower-dimensional) while still
+    being usable in constraint expressions that use .isel() for slicing.
+
+    Args:
+        data: DataArray or scalar value
+        indexers: Dictionary of {dim: indexer} for isel
+
+    Returns:
+        Sliced DataArray if dims exist, otherwise original data
+    """
+    if not isinstance(data, xr.DataArray):
+        return data
+    # Only apply isel if data has all the required dimensions
+    if all(dim in data.dims for dim in indexers):
+        return data.isel(indexers)
+    return data
+
+
+def _scalar_safe_isel_drop(data: xr.DataArray | Any, dim: str, index: int) -> xr.DataArray | Any:
+    """Apply isel with drop=True if data has the dimension, otherwise return data as-is.
+
+    Useful for cases like selecting the last value of a potentially reduced array:
+    - If data has time dimension: returns data.isel(time=-1, drop=True)
+    - If data is reduced (no time dimension): returns data unchanged (already represents constant)
+
+    Args:
+        data: DataArray or scalar value
+        dim: Dimension name to select from
+        index: Index to select (e.g., -1 for last, 0 for first)
+
+    Returns:
+        Selected value with dimension dropped if dim exists, otherwise original data
+    """
+    if not isinstance(data, xr.DataArray):
+        return data
+    if dim in data.dims:
+        return data.isel({dim: index}, drop=True)
+    return data
+
+
+def _scalar_safe_reduce(data: xr.DataArray | Any, dim: str, method: str = 'mean') -> xr.DataArray | Any:
+    """Apply reduction (mean/sum/etc) over dimension if it exists, otherwise return data as-is.
+
+    Useful for aggregating over time dimension when data may be scalar (constant):
+    - If data has time dimension: returns getattr(data, method)(dim)
+    - If data is reduced (no time dimension): returns data unchanged (already represents constant)
+
+    Args:
+        data: DataArray or scalar value
+        dim: Dimension name to reduce over
+        method: Reduction method ('mean', 'sum', 'min', 'max', etc.)
+
+    Returns:
+        Reduced value if dim exists, otherwise original data
+    """
+    if not isinstance(data, xr.DataArray):
+        return data
+    if dim in data.dims:
+        return getattr(data, method)(dim)
+    return data
+
+
 class ModelingUtilitiesAbstract:
     """Utility functions for modeling - leveraging xarray for temporal data"""
 
@@ -340,7 +406,7 @@ def consecutive_duration_tracking(
             constraints['lb'] = model.add_constraints(
                 duration
                 >= (state.isel({duration_dim: slice(None, -1)}) - state.isel({duration_dim: slice(1, None)}))
-                * minimum_duration.isel({duration_dim: slice(None, -1)}),
+                * _scalar_safe_isel(minimum_duration, {duration_dim: slice(None, -1)}),
                 name=f'{duration.name}|lb',
             )
 
@@ -351,7 +417,7 @@ def consecutive_duration_tracking(
                     if not isinstance(previous_duration, xr.DataArray)
                     else float(previous_duration.max().item())
                 )
-                min0 = float(minimum_duration.isel({duration_dim: 0}).max().item())
+                min0 = float(_scalar_safe_isel(minimum_duration, {duration_dim: 0}).max().item())
                 if prev > 0 and prev < min0:
                     constraints['initial_lb'] = model.add_constraints(
                         state.isel({duration_dim: 0}) == 1, name=f'{duration.name}|initial_lb'
diff --git a/flixopt/structure.py b/flixopt/structure.py
index 67a02750f..5333d37ae 100644
--- a/flixopt/structure.py
+++ b/flixopt/structure.py
@@ -40,6 +40,44 @@
 
 logger = logging.getLogger('flixopt')
 
+
+def _ensure_coords(
+    data: xr.DataArray | float | int,
+    coords: xr.Coordinates | dict,
+) -> xr.DataArray | float:
+    """Broadcast data to coords if needed.
+
+    This is used at the linopy interface to ensure bounds are properly broadcasted
+    to the target variable shape. Linopy needs at least one bound to have all
+    dimensions to determine the variable shape.
+
+    Note: Infinity values (-inf, inf) are kept as scalars because linopy uses
+    special checks like `if (lower != -inf)` that fail with DataArrays.
+    """
+    # Handle both dict and xr.Coordinates
+    if isinstance(coords, dict):
+        coord_dims = list(coords.keys())
+    else:
+        coord_dims = list(coords.dims)
+
+    # Keep infinity values as scalars (linopy uses them for special checks)
+    if not isinstance(data, xr.DataArray):
+        if np.isinf(data):
+            return data
+        # Finite scalar - create full DataArray
+        return xr.DataArray(data, coords=coords, dims=coord_dims)
+
+    if set(data.dims) == set(coord_dims):
+        # Has all dims - ensure correct order
+        if data.dims != tuple(coord_dims):
+            return data.transpose(*coord_dims)
+        return data
+
+    # Broadcast to full coords (broadcast_like ensures correct dim order)
+    template = xr.DataArray(coords=coords, dims=coord_dims)
+    return data.broadcast_like(template)
+
+
 CLASS_REGISTRY = {}
 
 
@@ -98,6 +136,24 @@ def __init__(self, flow_system: FlowSystem):
         self.effects: EffectCollectionModel | None = None
         self.submodels: Submodels = Submodels({})
 
+    def add_variables(
+        self,
+        lower: xr.DataArray | float = -np.inf,
+        upper: xr.DataArray | float = np.inf,
+        coords: xr.Coordinates | None = None,
+        **kwargs,
+    ) -> linopy.Variable:
+        """Override to ensure bounds are broadcasted to coords shape.
+
+        Linopy uses the union of all DataArray dimensions to determine variable shape.
+        This override ensures at least one bound has all target dimensions when coords
+        is provided, allowing internal data to remain compact (scalars, 1D arrays).
+        """
+        if coords is not None:
+            lower = _ensure_coords(lower, coords)
+            upper = _ensure_coords(upper, coords)
+        return super().add_variables(lower=lower, upper=upper, coords=coords, **kwargs)
+
     def do_modeling(self):
         # Create all element models
         self.effects = self.flow_system.effects.create_model(self)
diff --git a/flixopt/transform_accessor.py b/flixopt/transform_accessor.py
index 07381cf5f..854b23525 100644
--- a/flixopt/transform_accessor.py
+++ b/flixopt/transform_accessor.py
@@ -16,6 +16,8 @@
 import pandas as pd
 import xarray as xr
 
+from .modeling import _scalar_safe_reduce
+
 if TYPE_CHECKING:
     from .flow_system import FlowSystem
 
@@ -366,6 +368,26 @@ def _dataset_resample(
         time_var_names = [v for v in dataset.data_vars if 'time' in dataset[v].dims]
         non_time_var_names = [v for v in dataset.data_vars if v not in time_var_names]
 
+        # Handle case where no data variables have time dimension (all scalars)
+        # We still need to resample the time coordinate itself
+        if not time_var_names:
+            if 'time' not in dataset.coords:
+                raise ValueError('Dataset has no time dimension to resample')
+            # Create a dummy variable to resample the time coordinate
+            dummy = xr.DataArray(
+                np.zeros(len(dataset.coords['time'])), dims=['time'], coords={'time': dataset.coords['time']}
+            )
+            dummy_ds = xr.Dataset({'__dummy__': dummy})
+            resampled_dummy = getattr(dummy_ds.resample(time=freq, **kwargs), method)()
+            # Get the resampled time coordinate
+            resampled_time = resampled_dummy.coords['time']
+            # Create result with all original scalar data and resampled time coordinate
+            # Keep all existing coordinates (period, scenario, etc.) except time which gets resampled
+            result = dataset.copy()
+            result = result.assign_coords(time=resampled_time)
+            result.attrs.update(original_attrs)
+            return FlowSystem._update_time_metadata(result, hours_of_last_timestep, hours_of_previous_timesteps)
+
         time_dataset = dataset[time_var_names]
         resampled_time_dataset = cls._resample_by_dimension_groups(time_dataset, freq, method, **kwargs)
 
@@ -399,6 +421,12 @@ def _dataset_resample(
         else:
             result = resampled_time_dataset
 
+        # Preserve all original coordinates that aren't 'time' (e.g., period, scenario, cluster)
+        # These may be lost during merge if no data variable uses them
+        for coord_name, coord_val in dataset.coords.items():
+            if coord_name != 'time' and coord_name not in result.coords:
+                result = result.assign_coords({coord_name: coord_val})
+
         result.attrs.update(original_attrs)
         return FlowSystem._update_time_metadata(result, hours_of_last_timestep, hours_of_previous_timesteps)
 
@@ -1347,7 +1375,7 @@ def _apply_soc_decay(
         )
 
         # Decay factor: (1 - loss)^t
-        loss_value = storage.relative_loss_per_hour.mean('time')
+        loss_value = _scalar_safe_reduce(storage.relative_loss_per_hour, 'time', 'mean')
         if not np.any(loss_value.values > 0):
             return soc_boundary_per_timestep
 
diff --git a/tests/conftest.py b/tests/conftest.py
index ee2c0f4e2..9923af896 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -829,6 +829,25 @@ def assert_sets_equal(set1: Iterable, set2: Iterable, msg=''):
         raise AssertionError(error_msg)
 
 
+def assert_dims_compatible(data: xr.DataArray, model_coords: tuple[str, ...], msg: str = ''):
+    """Assert that data dimensions are a subset of model coordinates (compatible with broadcasting).
+
+    Parameters in flixopt now stay in minimal form (scalar, 1D, etc.) and are broadcast
+    at the linopy interface. This helper verifies that data dims are valid for the model.
+
+    Args:
+        data: DataArray to check
+        model_coords: Tuple of model coordinate names (from model.get_coords())
+        msg: Optional message for assertion error
+    """
+    extra_dims = set(data.dims) - set(model_coords)
+    if extra_dims:
+        error = f'Data has dimensions {extra_dims} not in model coordinates {model_coords}'
+        if msg:
+            error = f'{msg}: {error}'
+        raise AssertionError(error)
+
+
 # ============================================================================
 # PLOTTING CLEANUP FIXTURES
 # ============================================================================
diff --git a/tests/deprecated/test_component.py b/tests/deprecated/test_component.py
index 765a08c5a..f81ca270e 100644
--- a/tests/deprecated/test_component.py
+++ b/tests/deprecated/test_component.py
@@ -131,7 +131,8 @@ def test_on_with_multiple_flows(self, basic_flow_system_linopy_coords, coords_co
 
         upper_bound_flow_rate = outputs[1].relative_maximum
 
-        assert upper_bound_flow_rate.dims == tuple(model.get_coords())
+        # Data stays in minimal form (1D array stays 1D)
+        assert upper_bound_flow_rate.dims == ('time',)
 
         assert_var_equal(
             model['TestComponent(Out2)|flow_rate'],
@@ -311,7 +312,8 @@ def test_previous_states_with_multiple_flows(self, basic_flow_system_linopy_coor
 
         upper_bound_flow_rate = outputs[1].relative_maximum
 
-        assert upper_bound_flow_rate.dims == tuple(model.get_coords())
+        # Data stays in minimal form (1D array stays 1D)
+        assert upper_bound_flow_rate.dims == ('time',)
 
         assert_var_equal(
             model['TestComponent(Out2)|flow_rate'],
diff --git a/tests/deprecated/test_dataconverter.py b/tests/deprecated/test_dataconverter.py
deleted file mode 100644
index a5774fd6b..000000000
--- a/tests/deprecated/test_dataconverter.py
+++ /dev/null
@@ -1,1262 +0,0 @@
-import numpy as np
-import pandas as pd
-import pytest
-import xarray as xr
-
-from flixopt.core import (  # Adjust this import to match your project structure
-    ConversionError,
-    DataConverter,
-    TimeSeriesData,
-)
-
-
-@pytest.fixture
-def time_coords():
-    return pd.date_range('2024-01-01', periods=5, freq='D', name='time')
-
-
-@pytest.fixture
-def scenario_coords():
-    return pd.Index(['baseline', 'high', 'low'], name='scenario')
-
-
-@pytest.fixture
-def region_coords():
-    return pd.Index(['north', 'south', 'east'], name='region')
-
-
-@pytest.fixture
-def standard_coords():
-    """Standard coordinates with unique lengths for easy testing."""
-    return {
-        'time': pd.date_range('2024-01-01', periods=5, freq='D', name='time'),  # length 5
-        'scenario': pd.Index(['A', 'B', 'C'], name='scenario'),  # length 3
-        'region': pd.Index(['north', 'south'], name='region'),  # length 2
-    }
-
-
-class TestScalarConversion:
-    """Test scalar data conversions with different coordinate configurations."""
-
-    def test_scalar_no_coords(self):
-        """Scalar without coordinates should create 0D DataArray."""
-        result = DataConverter.to_dataarray(42)
-        assert result.shape == ()
-        assert result.dims == ()
-        assert result.item() == 42
-
-    def test_scalar_single_coord(self, time_coords):
-        """Scalar with single coordinate should broadcast."""
-        result = DataConverter.to_dataarray(42, coords={'time': time_coords})
-        assert result.shape == (5,)
-        assert result.dims == ('time',)
-        assert np.all(result.values == 42)
-
-    def test_scalar_multiple_coords(self, time_coords, scenario_coords):
-        """Scalar with multiple coordinates should broadcast to all."""
-        result = DataConverter.to_dataarray(42, coords={'time': time_coords, 'scenario': scenario_coords})
-        assert result.shape == (5, 3)
-        assert result.dims == ('time', 'scenario')
-        assert np.all(result.values == 42)
-
-    def test_numpy_scalars(self, time_coords):
-        """Test numpy scalar types."""
-        for scalar in [np.int32(42), np.int64(42), np.float32(42.5), np.float64(42.5)]:
-            result = DataConverter.to_dataarray(scalar, coords={'time': time_coords})
-            assert result.shape == (5,)
-            assert np.all(result.values == scalar.item())
-
-    def test_scalar_many_dimensions(self, standard_coords):
-        """Scalar should broadcast to any number of dimensions."""
-        coords = {**standard_coords, 'technology': pd.Index(['solar', 'wind'], name='technology')}
-
-        result = DataConverter.to_dataarray(42, coords=coords)
-        assert result.shape == (5, 3, 2, 2)
-        assert result.dims == ('time', 'scenario', 'region', 'technology')
-        assert np.all(result.values == 42)
-
-
-class TestOneDimensionalArrayConversion:
-    """Test 1D numpy array and pandas Series conversions."""
-
-    def test_1d_array_no_coords(self):
-        """1D array without coords should fail unless single element."""
-        # Multi-element fails
-        with pytest.raises(ConversionError):
-            DataConverter.to_dataarray(np.array([1, 2, 3]))
-
-        # Single element succeeds
-        result = DataConverter.to_dataarray(np.array([42]))
-        assert result.shape == ()
-        assert result.item() == 42
-
-    def test_1d_array_matching_coord(self, time_coords):
-        """1D array matching coordinate length should work."""
-        arr = np.array([10, 20, 30, 40, 50])
-        result = DataConverter.to_dataarray(arr, coords={'time': time_coords})
-        assert result.shape == (5,)
-        assert result.dims == ('time',)
-        assert np.array_equal(result.values, arr)
-
-    def test_1d_array_mismatched_coord(self, time_coords):
-        """1D array not matching coordinate length should fail."""
-        arr = np.array([10, 20, 30])  # Length 3, time_coords has length 5
-        with pytest.raises(ConversionError):
-            DataConverter.to_dataarray(arr, coords={'time': time_coords})
-
-    def test_1d_array_broadcast_to_multiple_coords(self, time_coords, scenario_coords):
-        """1D array should broadcast to matching dimension."""
-        # Array matching time dimension
-        time_arr = np.array([10, 20, 30, 40, 50])
-        result = DataConverter.to_dataarray(time_arr, coords={'time': time_coords, 'scenario': scenario_coords})
-        assert result.shape == (5, 3)
-        assert result.dims == ('time', 'scenario')
-
-        # Each scenario should have the same time values
-        for scenario in scenario_coords:
-            assert np.array_equal(result.sel(scenario=scenario).values, time_arr)
-
-        # Array matching scenario dimension
-        scenario_arr = np.array([100, 200, 300])
-        result = DataConverter.to_dataarray(scenario_arr, coords={'time': time_coords, 'scenario': scenario_coords})
-        assert result.shape == (5, 3)
-        assert result.dims == ('time', 'scenario')
-
-        # Each time should have the same scenario values
-        for time in time_coords:
-            assert np.array_equal(result.sel(time=time).values, scenario_arr)
-
-    def test_1d_array_ambiguous_length(self):
-        """Array length matching multiple dimensions should fail."""
-        # Both dimensions have length 3
-        coords_3x3 = {
-            'time': pd.date_range('2024-01-01', periods=3, freq='D', name='time'),
-            'scenario': pd.Index(['A', 'B', 'C'], name='scenario'),
-        }
-        arr = np.array([1, 2, 3])
-
-        with pytest.raises(ConversionError, match='matches multiple dimension'):
-            DataConverter.to_dataarray(arr, coords=coords_3x3)
-
-    def test_1d_array_broadcast_to_many_dimensions(self, standard_coords):
-        """1D array should broadcast to many dimensions."""
-        # Array matching time dimension
-        time_arr = np.array([10, 20, 30, 40, 50])
-        result = DataConverter.to_dataarray(time_arr, coords=standard_coords)
-
-        assert result.shape == (5, 3, 2)
-        assert result.dims == ('time', 'scenario', 'region')
-
-        # Check broadcasting - all scenarios and regions should have same time values
-        for scenario in standard_coords['scenario']:
-            for region in standard_coords['region']:
-                assert np.array_equal(result.sel(scenario=scenario, region=region).values, time_arr)
-
-
-class TestSeriesConversion:
-    """Test pandas Series conversions."""
-
-    def test_series_no_coords(self):
-        """Series without coords should fail unless single element."""
-        # Multi-element fails
-        series = pd.Series([1, 2, 3])
-        with pytest.raises(ConversionError):
-            DataConverter.to_dataarray(series)
-
-        # Single element succeeds
-        single_series = pd.Series([42])
-        result = DataConverter.to_dataarray(single_series)
-        assert result.shape == ()
-        assert result.item() == 42
-
-    def test_series_matching_index(self, time_coords, scenario_coords):
-        """Series with matching index should work."""
-        # Time-indexed series
-        time_series = pd.Series([10, 20, 30, 40, 50], index=time_coords)
-        result = DataConverter.to_dataarray(time_series, coords={'time': time_coords})
-        assert result.shape == (5,)
-        assert result.dims == ('time',)
-        assert np.array_equal(result.values, time_series.values)
-
-        # Scenario-indexed series
-        scenario_series = pd.Series([100, 200, 300], index=scenario_coords)
-        result = DataConverter.to_dataarray(scenario_series, coords={'scenario': scenario_coords})
-        assert result.shape == (3,)
-        assert result.dims == ('scenario',)
-        assert np.array_equal(result.values, scenario_series.values)
-
-    def test_series_mismatched_index(self, time_coords):
-        """Series with non-matching index should fail."""
-        wrong_times = pd.date_range('2025-01-01', periods=5, freq='D', name='time')
-        series = pd.Series([10, 20, 30, 40, 50], index=wrong_times)
-
-        with pytest.raises(ConversionError):
-            DataConverter.to_dataarray(series, coords={'time': time_coords})
-
-    def test_series_broadcast_to_multiple_coords(self, time_coords, scenario_coords):
-        """Series should broadcast to non-matching dimensions."""
-        # Time series broadcast to scenarios
-        time_series = pd.Series([10, 20, 30, 40, 50], index=time_coords)
-        result = DataConverter.to_dataarray(time_series, coords={'time': time_coords, 'scenario': scenario_coords})
-        assert result.shape == (5, 3)
-
-        for scenario in scenario_coords:
-            assert np.array_equal(result.sel(scenario=scenario).values, time_series.values)
-
-    def test_series_wrong_dimension(self, time_coords, region_coords):
-        """Series indexed by dimension not in coords should fail."""
-        wrong_series = pd.Series([1, 2, 3], index=region_coords)
-
-        with pytest.raises(ConversionError):
-            DataConverter.to_dataarray(wrong_series, coords={'time': time_coords})
-
-    def test_series_broadcast_to_many_dimensions(self, standard_coords):
-        """Series should broadcast to many dimensions."""
-        time_series = pd.Series([100, 200, 300, 400, 500], index=standard_coords['time'])
-        result = DataConverter.to_dataarray(time_series, coords=standard_coords)
-
-        assert result.shape == (5, 3, 2)
-        assert result.dims == ('time', 'scenario', 'region')
-
-        # Check that all non-time dimensions have the same time series values
-        for scenario in standard_coords['scenario']:
-            for region in standard_coords['region']:
-                assert np.array_equal(result.sel(scenario=scenario, region=region).values, time_series.values)
-
-
-class TestDataFrameConversion:
-    """Test pandas DataFrame conversions."""
-
-    def test_single_column_dataframe(self, time_coords):
-        """Single-column DataFrame should work like Series."""
-        df = pd.DataFrame({'value': [10, 20, 30, 40, 50]}, index=time_coords)
-        result = DataConverter.to_dataarray(df, coords={'time': time_coords})
-
-        assert result.shape == (5,)
-        assert result.dims == ('time',)
-        assert np.array_equal(result.values, df['value'].values)
-
-    def test_multi_column_dataframe_accepted(self, time_coords, scenario_coords):
-        """Multi-column DataFrame should now be accepted and converted via numpy array path."""
-        df = pd.DataFrame(
-            {'value1': [10, 20, 30, 40, 50], 'value2': [15, 25, 35, 45, 55], 'value3': [12, 22, 32, 42, 52]},
-            index=time_coords,
-        )
-
-        # Should work by converting to numpy array (5x3) and matching to time x scenario
-        result = DataConverter.to_dataarray(df, coords={'time': time_coords, 'scenario': scenario_coords})
-
-        assert result.shape == (5, 3)
-        assert result.dims == ('time', 'scenario')
-        assert np.array_equal(result.values, df.to_numpy())
-
-    def test_empty_dataframe_rejected(self, time_coords):
-        """Empty DataFrame should be rejected."""
-        df = pd.DataFrame(index=time_coords)  # No columns
-
-        with pytest.raises(ConversionError, match='DataFrame must have at least one column'):
-            DataConverter.to_dataarray(df, coords={'time': time_coords})
-
-    def test_dataframe_broadcast(self, time_coords, scenario_coords):
-        """Single-column DataFrame should broadcast like Series."""
-        df = pd.DataFrame({'power': [10, 20, 30, 40, 50]}, index=time_coords)
-        result = DataConverter.to_dataarray(df, coords={'time': time_coords, 'scenario': scenario_coords})
-
-        assert result.shape == (5, 3)
-        for scenario in scenario_coords:
-            assert np.array_equal(result.sel(scenario=scenario).values, df['power'].values)
-
-
-class TestMultiDimensionalArrayConversion:
-    """Test multi-dimensional numpy array conversions."""
-
-    def test_2d_array_unique_dimensions(self, standard_coords):
-        """2D array with unique dimension lengths should work."""
-        # 5x3 array should map to time x scenario
-        data_2d = np.random.rand(5, 3)
-        result = DataConverter.to_dataarray(
-            data_2d, coords={'time': standard_coords['time'], 'scenario': standard_coords['scenario']}
-        )
-
-        assert result.shape == (5, 3)
-        assert result.dims == ('time', 'scenario')
-        assert np.array_equal(result.values, data_2d)
-
-        # 3x5 array should map to scenario x time
-        data_2d_flipped = np.random.rand(3, 5)
-        result_flipped = DataConverter.to_dataarray(
-            data_2d_flipped, coords={'time': standard_coords['time'], 'scenario': standard_coords['scenario']}
-        )
-
-        assert result_flipped.shape == (5, 3)
-        assert result_flipped.dims == ('time', 'scenario')
-        assert np.array_equal(result_flipped.values.transpose(), data_2d_flipped)
-
-    def test_2d_array_broadcast_to_3d(self, standard_coords):
-        """2D array should broadcast to additional dimensions when using partial matching."""
-        # With improved integration, 2D array (5x3) should match time×scenario and broadcast to region
-        data_2d = np.random.rand(5, 3)
-        result = DataConverter.to_dataarray(data_2d, coords=standard_coords)
-
-        assert result.shape == (5, 3, 2)
-        assert result.dims == ('time', 'scenario', 'region')
-
-        # Check that all regions have the same time x scenario data
-        for region in standard_coords['region']:
-            assert np.array_equal(result.sel(region=region).values, data_2d)
-
-    def test_3d_array_unique_dimensions(self, standard_coords):
-        """3D array with unique dimension lengths should work."""
-        # 5x3x2 array should map to time x scenario x region
-        data_3d = np.random.rand(5, 3, 2)
-        result = DataConverter.to_dataarray(data_3d, coords=standard_coords)
-
-        assert result.shape == (5, 3, 2)
-        assert result.dims == ('time', 'scenario', 'region')
-        assert np.array_equal(result.values, data_3d)
-
-    def test_3d_array_different_permutation(self, standard_coords):
-        """3D array with different dimension order should work."""
-        # 2x5x3 array should map to region x time x scenario
-        data_3d = np.random.rand(2, 5, 3)
-        result = DataConverter.to_dataarray(data_3d, coords=standard_coords)
-
-        assert result.shape == (5, 3, 2)
-        assert result.dims == ('time', 'scenario', 'region')
-        assert np.array_equal(result.transpose('region', 'time', 'scenario').values, data_3d)
-
-    def test_4d_array_unique_dimensions(self):
-        """4D array with unique dimension lengths should work."""
-        coords = {
-            'time': pd.date_range('2024-01-01', periods=2, freq='D', name='time'),  # length 2
-            'scenario': pd.Index(['A', 'B', 'C'], name='scenario'),  # length 3
-            'region': pd.Index(['north', 'south', 'east', 'west'], name='region'),  # length 4
-            'technology': pd.Index(['solar', 'wind', 'gas', 'coal', 'hydro'], name='technology'),  # length 5
-        }
-
-        # 3x5x2x4 array should map to scenario x technology x time x region
-        data_4d = np.random.rand(3, 5, 2, 4)
-        result = DataConverter.to_dataarray(data_4d, coords=coords)
-
-        assert result.shape == (2, 3, 4, 5)
-        assert result.dims == ('time', 'scenario', 'region', 'technology')
-        assert np.array_equal(result.transpose('scenario', 'technology', 'time', 'region').values, data_4d)
-
-    def test_2d_array_ambiguous_dimensions_error(self):
-        """2D array with ambiguous dimension lengths should fail."""
-        # Both dimensions have length 3
-        coords_ambiguous = {
-            'scenario': pd.Index(['A', 'B', 'C'], name='scenario'),  # length 3
-            'region': pd.Index(['north', 'south', 'east'], name='region'),  # length 3
-        }
-
-        data_2d = np.random.rand(3, 3)
-        with pytest.raises(ConversionError, match='matches multiple dimension combinations'):
-            DataConverter.to_dataarray(data_2d, coords=coords_ambiguous)
-
-    def test_multid_array_no_coords(self):
-        """Multi-D arrays without coords should fail unless scalar."""
-        # Multi-element fails
-        data_2d = np.random.rand(2, 3)
-        with pytest.raises(ConversionError, match='Cannot convert multi-element array without target dimensions'):
-            DataConverter.to_dataarray(data_2d)
-
-        # Single element succeeds
-        single_element = np.array([[42]])
-        result = DataConverter.to_dataarray(single_element)
-        assert result.shape == ()
-        assert result.item() == 42
-
-    def test_array_no_matching_dimensions_error(self, standard_coords):
-        """Array with no matching dimension lengths should fail."""
-        # 7x8 array - no dimension has length 7 or 8
-        data_2d = np.random.rand(7, 8)
-        coords_2d = {
-            'time': standard_coords['time'],  # length 5
-            'scenario': standard_coords['scenario'],  # length 3
-        }
-
-        with pytest.raises(ConversionError, match='cannot be mapped to any combination'):
-            DataConverter.to_dataarray(data_2d, coords=coords_2d)
-
-    def test_multid_array_special_values(self, standard_coords):
-        """Multi-D arrays should preserve special values."""
-        # Create 2D array with special values
-        data_2d = np.array(
-            [[1.0, np.nan, 3.0], [np.inf, 5.0, -np.inf], [7.0, 8.0, 9.0], [10.0, np.nan, 12.0], [13.0, 14.0, np.inf]]
-        )
-
-        result = DataConverter.to_dataarray(
-            data_2d, coords={'time': standard_coords['time'], 'scenario': standard_coords['scenario']}
-        )
-
-        assert result.shape == (5, 3)
-        assert np.array_equal(np.isnan(result.values), np.isnan(data_2d))
-        assert np.array_equal(np.isinf(result.values), np.isinf(data_2d))
-
-    def test_multid_array_dtype_preservation(self, standard_coords):
-        """Multi-D arrays should preserve data types."""
-        # Integer array
-        int_data = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9], [10, 11, 12], [13, 14, 15]], dtype=np.int32)
-
-        result_int = DataConverter.to_dataarray(
-            int_data, coords={'time': standard_coords['time'], 'scenario': standard_coords['scenario']}
-        )
-
-        assert result_int.dtype == np.int32
-        assert np.array_equal(result_int.values, int_data)
-
-        # Boolean array
-        bool_data = np.array(
-            [[True, False, True], [False, True, False], [True, True, False], [False, False, True], [True, False, True]]
-        )
-
-        result_bool = DataConverter.to_dataarray(
-            bool_data, coords={'time': standard_coords['time'], 'scenario': standard_coords['scenario']}
-        )
-
-        assert result_bool.dtype == bool
-        assert np.array_equal(result_bool.values, bool_data)
-
-
-class TestDataArrayConversion:
-    """Test xarray DataArray conversions."""
-
-    def test_compatible_dataarray(self, time_coords):
-        """Compatible DataArray should pass through."""
-        original = xr.DataArray([10, 20, 30, 40, 50], coords={'time': time_coords}, dims='time')
-        result = DataConverter.to_dataarray(original, coords={'time': time_coords})
-
-        assert result.shape == (5,)
-        assert result.dims == ('time',)
-        assert np.array_equal(result.values, original.values)
-
-        # Should be a copy
-        result[0] = 999
-        assert original[0].item() == 10
-
-    def test_incompatible_dataarray_coords(self, time_coords):
-        """DataArray with wrong coordinates should fail."""
-        wrong_times = pd.date_range('2025-01-01', periods=5, freq='D', name='time')
-        original = xr.DataArray([10, 20, 30, 40, 50], coords={'time': wrong_times}, dims='time')
-
-        with pytest.raises(ConversionError):
-            DataConverter.to_dataarray(original, coords={'time': time_coords})
-
-    def test_incompatible_dataarray_dims(self, time_coords):
-        """DataArray with wrong dimensions should fail."""
-        original = xr.DataArray([10, 20, 30, 40, 50], coords={'wrong_dim': range(5)}, dims='wrong_dim')
-
-        with pytest.raises(ConversionError):
-            DataConverter.to_dataarray(original, coords={'time': time_coords})
-
-    def test_dataarray_broadcast(self, time_coords, scenario_coords):
-        """DataArray should broadcast to additional dimensions."""
-        # 1D time DataArray to 2D time+scenario
-        original = xr.DataArray([10, 20, 30, 40, 50], coords={'time': time_coords}, dims='time')
-        result = DataConverter.to_dataarray(original, coords={'time': time_coords, 'scenario': scenario_coords})
-
-        assert result.shape == (5, 3)
-        assert result.dims == ('time', 'scenario')
-
-        for scenario in scenario_coords:
-            assert np.array_equal(result.sel(scenario=scenario).values, original.values)
-
-    def test_scalar_dataarray_broadcast(self, time_coords, scenario_coords):
-        """Scalar DataArray should broadcast to all dimensions."""
-        scalar_da = xr.DataArray(42)
-        result = DataConverter.to_dataarray(scalar_da, coords={'time': time_coords, 'scenario': scenario_coords})
-
-        assert result.shape == (5, 3)
-        assert np.all(result.values == 42)
-
-    def test_2d_dataarray_broadcast_to_more_dimensions(self, standard_coords):
-        """DataArray should broadcast to additional dimensions."""
-        # Start with 2D DataArray
-        original = xr.DataArray(
-            [[10, 20, 30], [40, 50, 60], [70, 80, 90], [100, 110, 120], [130, 140, 150]],
-            coords={'time': standard_coords['time'], 'scenario': standard_coords['scenario']},
-            dims=('time', 'scenario'),
-        )
-
-        # Broadcast to 3D
-        result = DataConverter.to_dataarray(original, coords=standard_coords)
-
-        assert result.shape == (5, 3, 2)
-        assert result.dims == ('time', 'scenario', 'region')
-
-        # Check that all regions have the same time+scenario values
-        for region in standard_coords['region']:
-            assert np.array_equal(result.sel(region=region).values, original.values)
-
-
-class TestTimeSeriesDataConversion:
-    """Test TimeSeriesData conversions."""
-
-    def test_timeseries_data_basic(self, time_coords):
-        """TimeSeriesData should work like DataArray."""
-        data_array = xr.DataArray([10, 20, 30, 40, 50], coords={'time': time_coords}, dims='time')
-        ts_data = TimeSeriesData(data_array, clustering_group='test')
-
-        result = DataConverter.to_dataarray(ts_data, coords={'time': time_coords})
-
-        assert result.shape == (5,)
-        assert result.dims == ('time',)
-        assert np.array_equal(result.values, [10, 20, 30, 40, 50])
-
-    def test_timeseries_data_broadcast(self, time_coords, scenario_coords):
-        """TimeSeriesData should broadcast to additional dimensions."""
-        data_array = xr.DataArray([10, 20, 30, 40, 50], coords={'time': time_coords}, dims='time')
-        ts_data = TimeSeriesData(data_array)
-
-        result = DataConverter.to_dataarray(ts_data, coords={'time': time_coords, 'scenario': scenario_coords})
-
-        assert result.shape == (5, 3)
-        for scenario in scenario_coords:
-            assert np.array_equal(result.sel(scenario=scenario).values, [10, 20, 30, 40, 50])
-
-
-class TestAsDataArrayAlias:
-    """Test that to_dataarray works as an alias for to_dataarray."""
-
-    def test_to_dataarray_is_alias(self, time_coords, scenario_coords):
-        """to_dataarray should work identically to to_dataarray."""
-        # Test with scalar
-        result_to = DataConverter.to_dataarray(42, coords={'time': time_coords})
-        result_as = DataConverter.to_dataarray(42, coords={'time': time_coords})
-        assert np.array_equal(result_to.values, result_as.values)
-        assert result_to.dims == result_as.dims
-        assert result_to.shape == result_as.shape
-
-        # Test with array
-        arr = np.array([10, 20, 30, 40, 50])
-        result_to_arr = DataConverter.to_dataarray(arr, coords={'time': time_coords})
-        result_as_arr = DataConverter.to_dataarray(arr, coords={'time': time_coords})
-        assert np.array_equal(result_to_arr.values, result_as_arr.values)
-        assert result_to_arr.dims == result_as_arr.dims
-
-        # Test with Series
-        series = pd.Series([100, 200, 300, 400, 500], index=time_coords)
-        result_to_series = DataConverter.to_dataarray(series, coords={'time': time_coords, 'scenario': scenario_coords})
-        result_as_series = DataConverter.to_dataarray(series, coords={'time': time_coords, 'scenario': scenario_coords})
-        assert np.array_equal(result_to_series.values, result_as_series.values)
-        assert result_to_series.dims == result_as_series.dims
-
-
-class TestCustomDimensions:
-    """Test with custom dimension names beyond time/scenario."""
-
-    def test_custom_single_dimension(self, region_coords):
-        """Test with custom dimension name."""
-        result = DataConverter.to_dataarray(42, coords={'region': region_coords})
-        assert result.shape == (3,)
-        assert result.dims == ('region',)
-        assert np.all(result.values == 42)
-
-    def test_custom_multiple_dimensions(self):
-        """Test with multiple custom dimensions."""
-        products = pd.Index(['A', 'B'], name='product')
-        technologies = pd.Index(['solar', 'wind', 'gas'], name='technology')
-
-        # Array matching technology dimension
-        arr = np.array([100, 150, 80])
-        result = DataConverter.to_dataarray(arr, coords={'product': products, 'technology': technologies})
-
-        assert result.shape == (2, 3)
-        assert result.dims == ('product', 'technology')
-
-        # Should broadcast across products
-        for product in products:
-            assert np.array_equal(result.sel(product=product).values, arr)
-
-    def test_mixed_dimension_types(self):
-        """Test mixing time dimension with custom dimensions."""
-        time_coords = pd.date_range('2024-01-01', periods=3, freq='D', name='time')
-        regions = pd.Index(['north', 'south'], name='region')
-
-        # Time series should broadcast to regions
-        time_series = pd.Series([10, 20, 30], index=time_coords)
-        result = DataConverter.to_dataarray(time_series, coords={'time': time_coords, 'region': regions})
-
-        assert result.shape == (3, 2)
-        assert result.dims == ('time', 'region')
-
-    def test_custom_dimensions_complex(self):
-        """Test complex scenario with custom dimensions."""
-        coords = {
-            'product': pd.Index(['A', 'B'], name='product'),
-            'factory': pd.Index(['F1', 'F2', 'F3'], name='factory'),
-            'quarter': pd.Index(['Q1', 'Q2', 'Q3', 'Q4'], name='quarter'),
-        }
-
-        # Array matching factory dimension
-        factory_arr = np.array([100, 200, 300])
-        result = DataConverter.to_dataarray(factory_arr, coords=coords)
-
-        assert result.shape == (2, 3, 4)
-        assert result.dims == ('product', 'factory', 'quarter')
-
-        # Check broadcasting
-        for product in coords['product']:
-            for quarter in coords['quarter']:
-                slice_data = result.sel(product=product, quarter=quarter)
-                assert np.array_equal(slice_data.values, factory_arr)
-
-
-class TestValidation:
-    """Test coordinate validation."""
-
-    def test_empty_coords(self):
-        """Empty coordinates should work for scalars."""
-        result = DataConverter.to_dataarray(42, coords={})
-        assert result.shape == ()
-        assert result.item() == 42
-
-    def test_invalid_coord_type(self):
-        """Non-pandas Index coordinates should fail."""
-        with pytest.raises(ConversionError):
-            DataConverter.to_dataarray(42, coords={'time': [1, 2, 3]})
-
-    def test_empty_coord_index(self):
-        """Empty coordinate index should fail."""
-        empty_index = pd.Index([], name='time')
-        with pytest.raises(ConversionError):
-            DataConverter.to_dataarray(42, coords={'time': empty_index})
-
-    def test_time_coord_validation(self):
-        """Time coordinates must be DatetimeIndex."""
-        # Non-datetime index with name 'time' should fail
-        wrong_time = pd.Index([1, 2, 3], name='time')
-        with pytest.raises(ConversionError, match='DatetimeIndex'):
-            DataConverter.to_dataarray(42, coords={'time': wrong_time})
-
-    def test_coord_naming(self, time_coords):
-        """Coordinates should be auto-renamed to match dimension."""
-        # Unnamed time index should be renamed
-        unnamed_time = time_coords.rename(None)
-        result = DataConverter.to_dataarray(42, coords={'time': unnamed_time})
-        assert result.coords['time'].name == 'time'
-
-
-class TestErrorHandling:
-    """Test error handling and edge cases."""
-
-    def test_unsupported_data_types(self, time_coords):
-        """Unsupported data types should fail with clear messages."""
-        unsupported = ['string', object(), None, {'dict': 'value'}, [1, 2, 3]]
-
-        for data in unsupported:
-            with pytest.raises(ConversionError):
-                DataConverter.to_dataarray(data, coords={'time': time_coords})
-
-    def test_dimension_mismatch_messages(self, time_coords, scenario_coords):
-        """Error messages should be informative."""
-        # Array with wrong length
-        wrong_arr = np.array([1, 2])  # Length 2, but no dimension has length 2
-        with pytest.raises(ConversionError, match='does not match any target dimension lengths'):
-            DataConverter.to_dataarray(wrong_arr, coords={'time': time_coords, 'scenario': scenario_coords})
-
-    def test_multidimensional_array_dimension_count_mismatch(self, standard_coords):
-        """Array with wrong number of dimensions should fail with clear error."""
-        # 4D array with 3D coordinates
-        data_4d = np.random.rand(5, 3, 2, 4)
-        with pytest.raises(ConversionError, match='cannot be mapped to any combination'):
-            DataConverter.to_dataarray(data_4d, coords=standard_coords)
-
-    def test_error_message_quality(self, standard_coords):
-        """Error messages should include helpful information."""
-        # Wrong shape array
-        data_2d = np.random.rand(7, 8)
-        coords_2d = {
-            'time': standard_coords['time'],  # length 5
-            'scenario': standard_coords['scenario'],  # length 3
-        }
-
-        try:
-            DataConverter.to_dataarray(data_2d, coords=coords_2d)
-            raise AssertionError('Should have raised ConversionError')
-        except ConversionError as e:
-            error_msg = str(e)
-            assert 'Array shape (7, 8)' in error_msg
-            assert 'target coordinate lengths:' in error_msg
-
-
-class TestDataIntegrity:
-    """Test data copying and integrity."""
-
-    def test_array_copy_independence(self, time_coords):
-        """Converted arrays should be independent copies."""
-        original_arr = np.array([10, 20, 30, 40, 50])
-        result = DataConverter.to_dataarray(original_arr, coords={'time': time_coords})
-
-        # Modify result
-        result[0] = 999
-
-        # Original should be unchanged
-        assert original_arr[0] == 10
-
-    def test_series_copy_independence(self, time_coords):
-        """Converted Series should be independent copies."""
-        original_series = pd.Series([10, 20, 30, 40, 50], index=time_coords)
-        result = DataConverter.to_dataarray(original_series, coords={'time': time_coords})
-
-        # Modify result
-        result[0] = 999
-
-        # Original should be unchanged
-        assert original_series.iloc[0] == 10
-
-    def test_dataframe_copy_independence(self, time_coords):
-        """Converted DataFrames should be independent copies."""
-        original_df = pd.DataFrame({'value': [10, 20, 30, 40, 50]}, index=time_coords)
-        result = DataConverter.to_dataarray(original_df, coords={'time': time_coords})
-
-        # Modify result
-        result[0] = 999
-
-        # Original should be unchanged
-        assert original_df.loc[time_coords[0], 'value'] == 10
-
-    def test_multid_array_copy_independence(self, standard_coords):
-        """Multi-D arrays should be independent copies."""
-        original_data = np.random.rand(5, 3)
-        result = DataConverter.to_dataarray(
-            original_data, coords={'time': standard_coords['time'], 'scenario': standard_coords['scenario']}
-        )
-
-        # Modify result
-        result[0, 0] = 999
-
-        # Original should be unchanged
-        assert original_data[0, 0] != 999
-
-
-class TestBooleanValues:
-    """Test handling of boolean values and arrays."""
-
-    def test_scalar_boolean_to_dataarray(self, time_coords):
-        """Scalar boolean values should work with to_dataarray."""
-        result_true = DataConverter.to_dataarray(True, coords={'time': time_coords})
-        assert result_true.shape == (5,)
-        assert result_true.dtype == bool
-        assert np.all(result_true.values)
-
-        result_false = DataConverter.to_dataarray(False, coords={'time': time_coords})
-        assert result_false.shape == (5,)
-        assert result_false.dtype == bool
-        assert not np.any(result_false.values)
-
-    def test_numpy_boolean_scalar(self, time_coords):
-        """Numpy boolean scalars should work."""
-        result_np_true = DataConverter.to_dataarray(np.bool_(True), coords={'time': time_coords})
-        assert result_np_true.shape == (5,)
-        assert result_np_true.dtype == bool
-        assert np.all(result_np_true.values)
-
-        result_np_false = DataConverter.to_dataarray(np.bool_(False), coords={'time': time_coords})
-        assert result_np_false.shape == (5,)
-        assert result_np_false.dtype == bool
-        assert not np.any(result_np_false.values)
-
-    def test_boolean_array_to_dataarray(self, time_coords):
-        """Boolean arrays should work with to_dataarray."""
-        bool_arr = np.array([True, False, True, False, True])
-        result = DataConverter.to_dataarray(bool_arr, coords={'time': time_coords})
-        assert result.shape == (5,)
-        assert result.dims == ('time',)
-        assert result.dtype == bool
-        assert np.array_equal(result.values, bool_arr)
-
-    def test_boolean_no_coords(self):
-        """Boolean scalar without coordinates should create 0D DataArray."""
-        result = DataConverter.to_dataarray(True)
-        assert result.shape == ()
-        assert result.dims == ()
-        assert result.item()
-
-        result_as = DataConverter.to_dataarray(False)
-        assert result_as.shape == ()
-        assert result_as.dims == ()
-        assert not result_as.item()
-
-    def test_boolean_multidimensional_broadcast(self, standard_coords):
-        """Boolean values should broadcast to multiple dimensions."""
-        result = DataConverter.to_dataarray(True, coords=standard_coords)
-        assert result.shape == (5, 3, 2)
-        assert result.dims == ('time', 'scenario', 'region')
-        assert result.dtype == bool
-        assert np.all(result.values)
-
-        result_as = DataConverter.to_dataarray(False, coords=standard_coords)
-        assert result_as.shape == (5, 3, 2)
-        assert result_as.dims == ('time', 'scenario', 'region')
-        assert result_as.dtype == bool
-        assert not np.any(result_as.values)
-
-    def test_boolean_series(self, time_coords):
-        """Boolean Series should work."""
-        bool_series = pd.Series([True, False, True, False, True], index=time_coords)
-        result = DataConverter.to_dataarray(bool_series, coords={'time': time_coords})
-        assert result.shape == (5,)
-        assert result.dtype == bool
-        assert np.array_equal(result.values, bool_series.values)
-
-        result_as = DataConverter.to_dataarray(bool_series, coords={'time': time_coords})
-        assert result_as.shape == (5,)
-        assert result_as.dtype == bool
-        assert np.array_equal(result_as.values, bool_series.values)
-
-    def test_boolean_dataframe(self, time_coords):
-        """Boolean DataFrame should work."""
-        bool_df = pd.DataFrame({'values': [True, False, True, False, True]}, index=time_coords)
-        result = DataConverter.to_dataarray(bool_df, coords={'time': time_coords})
-        assert result.shape == (5,)
-        assert result.dtype == bool
-        assert np.array_equal(result.values, bool_df['values'].values)
-
-        result_as = DataConverter.to_dataarray(bool_df, coords={'time': time_coords})
-        assert result_as.shape == (5,)
-        assert result_as.dtype == bool
-        assert np.array_equal(result_as.values, bool_df['values'].values)
-
-    def test_multidimensional_boolean_array(self, standard_coords):
-        """Multi-dimensional boolean arrays should work."""
-        bool_data = np.array(
-            [[True, False, True], [False, True, False], [True, True, False], [False, False, True], [True, False, True]]
-        )
-        result = DataConverter.to_dataarray(
-            bool_data, coords={'time': standard_coords['time'], 'scenario': standard_coords['scenario']}
-        )
-        assert result.shape == (5, 3)
-        assert result.dtype == bool
-        assert np.array_equal(result.values, bool_data)
-
-        result_as = DataConverter.to_dataarray(
-            bool_data, coords={'time': standard_coords['time'], 'scenario': standard_coords['scenario']}
-        )
-        assert result_as.shape == (5, 3)
-        assert result_as.dtype == bool
-        assert np.array_equal(result_as.values, bool_data)
-
-
-class TestSpecialValues:
-    """Test handling of special numeric values."""
-
-    def test_nan_values(self, time_coords):
-        """NaN values should be preserved."""
-        arr_with_nan = np.array([1, np.nan, 3, np.nan, 5])
-        result = DataConverter.to_dataarray(arr_with_nan, coords={'time': time_coords})
-
-        assert np.array_equal(np.isnan(result.values), np.isnan(arr_with_nan))
-        assert np.array_equal(result.values[~np.isnan(result.values)], arr_with_nan[~np.isnan(arr_with_nan)])
-
-    def test_infinite_values(self, time_coords):
-        """Infinite values should be preserved."""
-        arr_with_inf = np.array([1, np.inf, 3, -np.inf, 5])
-        result = DataConverter.to_dataarray(arr_with_inf, coords={'time': time_coords})
-
-        assert np.array_equal(result.values, arr_with_inf)
-
-    def test_boolean_values(self, time_coords):
-        """Boolean values should be preserved."""
-        bool_arr = np.array([True, False, True, False, True])
-        result = DataConverter.to_dataarray(bool_arr, coords={'time': time_coords})
-
-        assert result.dtype == bool
-        assert np.array_equal(result.values, bool_arr)
-
-    def test_mixed_numeric_types(self, time_coords):
-        """Mixed integer/float should become float."""
-        mixed_arr = np.array([1, 2.5, 3, 4.5, 5])
-        result = DataConverter.to_dataarray(mixed_arr, coords={'time': time_coords})
-
-        assert np.issubdtype(result.dtype, np.floating)
-        assert np.array_equal(result.values, mixed_arr)
-
-    def test_special_values_in_multid_arrays(self, standard_coords):
-        """Special values should be preserved in multi-D arrays and broadcasting."""
-        # Array with NaN and inf
-        special_arr = np.array([1, np.nan, np.inf, -np.inf, 5])
-        result = DataConverter.to_dataarray(special_arr, coords=standard_coords)
-
-        assert result.shape == (5, 3, 2)
-
-        # Check that special values are preserved in all broadcasts
-        for scenario in standard_coords['scenario']:
-            for region in standard_coords['region']:
-                slice_data = result.sel(scenario=scenario, region=region)
-                assert np.array_equal(np.isnan(slice_data.values), np.isnan(special_arr))
-                assert np.array_equal(np.isinf(slice_data.values), np.isinf(special_arr))
-
-
-class TestAdvancedBroadcasting:
-    """Test advanced broadcasting scenarios and edge cases."""
-
-    def test_partial_dimension_matching_with_broadcasting(self, standard_coords):
-        """Test that partial dimension matching works with the improved integration."""
-        # 1D array matching one dimension should broadcast to all target dimensions
-        time_arr = np.array([10, 20, 30, 40, 50])  # matches time (length 5)
-        result = DataConverter.to_dataarray(time_arr, coords=standard_coords)
-
-        assert result.shape == (5, 3, 2)
-        assert result.dims == ('time', 'scenario', 'region')
-
-        # Verify broadcasting
-        for scenario in standard_coords['scenario']:
-            for region in standard_coords['region']:
-                assert np.array_equal(result.sel(scenario=scenario, region=region).values, time_arr)
-
-    def test_complex_multid_scenario(self):
-        """Complex real-world scenario with multi-D array and broadcasting."""
-        # Energy system data: time x technology, broadcast to regions
-        coords = {
-            'time': pd.date_range('2024-01-01', periods=24, freq='h', name='time'),  # 24 hours
-            'technology': pd.Index(['solar', 'wind', 'gas', 'coal'], name='technology'),  # 4 technologies
-            'region': pd.Index(['north', 'south', 'east'], name='region'),  # 3 regions
-        }
-
-        # Capacity factors: 24 x 4 (will broadcast to 24 x 4 x 3)
-        capacity_factors = np.random.rand(24, 4)
-
-        result = DataConverter.to_dataarray(capacity_factors, coords=coords)
-
-        assert result.shape == (24, 4, 3)
-        assert result.dims == ('time', 'technology', 'region')
-        assert isinstance(result.indexes['time'], pd.DatetimeIndex)
-
-        # Verify broadcasting: all regions should have same time×technology data
-        for region in coords['region']:
-            assert np.array_equal(result.sel(region=region).values, capacity_factors)
-
-    def test_ambiguous_length_handling(self):
-        """Test handling of ambiguous length scenarios across different data types."""
-        # All dimensions have length 3
-        coords_3x3x3 = {
-            'time': pd.date_range('2024-01-01', periods=3, freq='D', name='time'),
-            'scenario': pd.Index(['A', 'B', 'C'], name='scenario'),
-            'region': pd.Index(['X', 'Y', 'Z'], name='region'),
-        }
-
-        # 1D array - should fail
-        arr_1d = np.array([1, 2, 3])
-        with pytest.raises(ConversionError, match='matches multiple dimension'):
-            DataConverter.to_dataarray(arr_1d, coords=coords_3x3x3)
-
-        # 2D array - should fail
-        arr_2d = np.random.rand(3, 3)
-        with pytest.raises(ConversionError, match='matches multiple dimension'):
-            DataConverter.to_dataarray(arr_2d, coords=coords_3x3x3)
-
-        # 3D array - should fail
-        arr_3d = np.random.rand(3, 3, 3)
-        with pytest.raises(ConversionError, match='matches multiple dimension'):
-            DataConverter.to_dataarray(arr_3d, coords=coords_3x3x3)
-
-    def test_mixed_broadcasting_scenarios(self):
-        """Test various broadcasting scenarios with different input types."""
-        coords = {
-            'time': pd.date_range('2024-01-01', periods=4, freq='D', name='time'),  # length 4
-            'scenario': pd.Index(['A', 'B'], name='scenario'),  # length 2
-            'region': pd.Index(['north', 'south', 'east'], name='region'),  # length 3
-            'product': pd.Index(['X', 'Y', 'Z', 'W', 'V'], name='product'),  # length 5
-        }
-
-        # Scalar to 4D
-        scalar_result = DataConverter.to_dataarray(42, coords=coords)
-        assert scalar_result.shape == (4, 2, 3, 5)
-        assert np.all(scalar_result.values == 42)
-
-        # 1D array (length 4, matches time) to 4D
-        arr_1d = np.array([10, 20, 30, 40])
-        arr_result = DataConverter.to_dataarray(arr_1d, coords=coords)
-        assert arr_result.shape == (4, 2, 3, 5)
-        # Verify broadcasting
-        for scenario in coords['scenario']:
-            for region in coords['region']:
-                for product in coords['product']:
-                    assert np.array_equal(
-                        arr_result.sel(scenario=scenario, region=region, product=product).values, arr_1d
-                    )
-
-        # 2D array (4x2, matches time×scenario) to 4D
-        arr_2d = np.random.rand(4, 2)
-        arr_2d_result = DataConverter.to_dataarray(arr_2d, coords=coords)
-        assert arr_2d_result.shape == (4, 2, 3, 5)
-        # Verify broadcasting
-        for region in coords['region']:
-            for product in coords['product']:
-                assert np.array_equal(arr_2d_result.sel(region=region, product=product).values, arr_2d)
-
-
-class TestAmbiguousDimensionLengthHandling:
-    """Test that DataConverter correctly raises errors when multiple dimensions have the same length."""
-
-    def test_1d_array_ambiguous_dimensions_simple(self):
-        """Test 1D array with two dimensions of same length should fail."""
-        # Both dimensions have length 3
-        coords_ambiguous = {
-            'scenario': pd.Index(['A', 'B', 'C'], name='scenario'),  # length 3
-            'region': pd.Index(['north', 'south', 'east'], name='region'),  # length 3
-        }
-
-        arr_1d = np.array([1, 2, 3])  # length 3 - matches both dimensions
-
-        with pytest.raises(ConversionError, match='matches multiple dimension'):
-            DataConverter.to_dataarray(arr_1d, coords=coords_ambiguous)
-
-    def test_1d_array_ambiguous_dimensions_complex(self):
-        """Test 1D array with multiple dimensions of same length."""
-        # Three dimensions have length 4
-        coords_4x4x4 = {
-            'time': pd.date_range('2024-01-01', periods=4, freq='D', name='time'),  # length 4
-            'scenario': pd.Index(['A', 'B', 'C', 'D'], name='scenario'),  # length 4
-            'region': pd.Index(['north', 'south', 'east', 'west'], name='region'),  # length 4
-            'product': pd.Index(['X', 'Y'], name='product'),  # length 2 - unique
-        }
-
-        # Array matching the ambiguous length
-        arr_1d = np.array([10, 20, 30, 40])  # length 4 - matches time, scenario, region
-
-        with pytest.raises(ConversionError, match='matches multiple dimension'):
-            DataConverter.to_dataarray(arr_1d, coords=coords_4x4x4)
-
-        # Array matching the unique length should work
-        arr_1d_unique = np.array([100, 200])  # length 2 - matches only product
-        result = DataConverter.to_dataarray(arr_1d_unique, coords=coords_4x4x4)
-        assert result.shape == (4, 4, 4, 2)  # broadcast to all dimensions
-        assert result.dims == ('time', 'scenario', 'region', 'product')
-
-    def test_2d_array_ambiguous_dimensions_both_same(self):
-        """Test 2D array where both dimensions have the same ambiguous length."""
-        # All dimensions have length 3
-        coords_3x3x3 = {
-            'time': pd.date_range('2024-01-01', periods=3, freq='D', name='time'),  # length 3
-            'scenario': pd.Index(['A', 'B', 'C'], name='scenario'),  # length 3
-            'region': pd.Index(['X', 'Y', 'Z'], name='region'),  # length 3
-        }
-
-        # 3x3 array - could be any combination of the three dimensions
-        arr_2d = np.random.rand(3, 3)
-
-        with pytest.raises(ConversionError, match='matches multiple dimension'):
-            DataConverter.to_dataarray(arr_2d, coords=coords_3x3x3)
-
-    def test_2d_array_one_dimension_ambiguous(self):
-        """Test 2D array where only one dimension length is ambiguous."""
-        coords_mixed = {
-            'time': pd.date_range('2024-01-01', periods=5, freq='D', name='time'),  # length 5 - unique
-            'scenario': pd.Index(['A', 'B', 'C'], name='scenario'),  # length 3
-            'region': pd.Index(['X', 'Y', 'Z'], name='region'),  # length 3 - same as scenario
-            'product': pd.Index(['P1', 'P2'], name='product'),  # length 2 - unique
-        }
-
-        # 5x3 array - first dimension clearly maps to time (unique length 5)
-        # but second dimension could be scenario or region (both length 3)
-        arr_5x3 = np.random.rand(5, 3)
-
-        with pytest.raises(ConversionError, match='matches multiple dimension'):
-            DataConverter.to_dataarray(arr_5x3, coords=coords_mixed)
-
-        # 5x2 array should work - dimensions are unambiguous
-        arr_5x2 = np.random.rand(5, 2)
-        result = DataConverter.to_dataarray(
-            arr_5x2, coords={'time': coords_mixed['time'], 'product': coords_mixed['product']}
-        )
-        assert result.shape == (5, 2)
-        assert result.dims == ('time', 'product')
-
-    def test_3d_array_all_dimensions_ambiguous(self):
-        """Test 3D array where all dimension lengths are ambiguous."""
-        # All dimensions have length 2
-        coords_2x2x2x2 = {
-            'scenario': pd.Index(['A', 'B'], name='scenario'),  # length 2
-            'region': pd.Index(['north', 'south'], name='region'),  # length 2
-            'technology': pd.Index(['solar', 'wind'], name='technology'),  # length 2
-            'product': pd.Index(['X', 'Y'], name='product'),  # length 2
-        }
-
-        # 2x2x2 array - could be any combination of 3 dimensions from the 4 available
-        arr_3d = np.random.rand(2, 2, 2)
-
-        with pytest.raises(ConversionError, match='matches multiple dimension'):
-            DataConverter.to_dataarray(arr_3d, coords=coords_2x2x2x2)
-
-    def test_3d_array_partial_ambiguity(self):
-        """Test 3D array with partial dimension ambiguity."""
-        coords_partial = {
-            'time': pd.date_range('2024-01-01', periods=4, freq='D', name='time'),  # length 4 - unique
-            'scenario': pd.Index(['A', 'B', 'C'], name='scenario'),  # length 3
-            'region': pd.Index(['X', 'Y', 'Z'], name='region'),  # length 3 - same as scenario
-            'technology': pd.Index(['solar', 'wind'], name='technology'),  # length 2 - unique
-        }
-
-        # 4x3x2 array - first and third dimensions are unique, middle is ambiguous
-        # This should still fail because middle dimension (length 3) could be scenario or region
-        arr_4x3x2 = np.random.rand(4, 3, 2)
-
-        with pytest.raises(ConversionError, match='matches multiple dimension'):
-            DataConverter.to_dataarray(arr_4x3x2, coords=coords_partial)
-
-    def test_pandas_series_ambiguous_dimensions(self):
-        """Test pandas Series with ambiguous dimension lengths."""
-        coords_ambiguous = {
-            'scenario': pd.Index(['A', 'B', 'C'], name='scenario'),  # length 3
-            'region': pd.Index(['north', 'south', 'east'], name='region'),  # length 3
-        }
-
-        # Series with length 3 but index that doesn't match either coordinate exactly
-        generic_series = pd.Series([10, 20, 30], index=[0, 1, 2])
-
-        # Should fail because length matches multiple dimensions and index doesn't match any
-        with pytest.raises(ConversionError, match='Series index does not match any target dimension coordinates'):
-            DataConverter.to_dataarray(generic_series, coords=coords_ambiguous)
-
-        # Series with index that matches one of the ambiguous coordinates should work
-        scenario_series = pd.Series([10, 20, 30], index=coords_ambiguous['scenario'])
-        result = DataConverter.to_dataarray(scenario_series, coords=coords_ambiguous)
-        assert result.shape == (3, 3)  # should broadcast to both dimensions
-        assert result.dims == ('scenario', 'region')
-
-    def test_edge_case_many_same_lengths(self):
-        """Test edge case with many dimensions having the same length."""
-        # Five dimensions all have length 2
-        coords_many = {
-            'dim1': pd.Index(['A', 'B'], name='dim1'),
-            'dim2': pd.Index(['X', 'Y'], name='dim2'),
-            'dim3': pd.Index(['P', 'Q'], name='dim3'),
-            'dim4': pd.Index(['M', 'N'], name='dim4'),
-            'dim5': pd.Index(['U', 'V'], name='dim5'),
-        }
-
-        # 1D array
-        arr_1d = np.array([1, 2])
-        with pytest.raises(ConversionError, match='matches multiple dimension'):
-            DataConverter.to_dataarray(arr_1d, coords=coords_many)
-
-        # 2D array
-        arr_2d = np.random.rand(2, 2)
-        with pytest.raises(ConversionError, match='matches multiple dimension'):
-            DataConverter.to_dataarray(arr_2d, coords=coords_many)
-
-        # 3D array
-        arr_3d = np.random.rand(2, 2, 2)
-        with pytest.raises(ConversionError, match='matches multiple dimension'):
-            DataConverter.to_dataarray(arr_3d, coords=coords_many)
-
-    def test_mixed_lengths_with_duplicates(self):
-        """Test mixed scenario with some duplicate and some unique lengths."""
-        coords_mixed = {
-            'time': pd.date_range('2024-01-01', periods=8, freq='D', name='time'),  # length 8 - unique
-            'scenario': pd.Index(['A', 'B', 'C'], name='scenario'),  # length 3
-            'region': pd.Index(['X', 'Y', 'Z'], name='region'),  # length 3 - same as scenario
-            'technology': pd.Index(['solar'], name='technology'),  # length 1 - unique
-            'product': pd.Index(['P1', 'P2', 'P3', 'P4', 'P5'], name='product'),  # length 5 - unique
-        }
-
-        # Arrays with unique lengths should work
-        arr_8 = np.arange(8)
-        result_8 = DataConverter.to_dataarray(arr_8, coords=coords_mixed)
-        assert result_8.dims == ('time', 'scenario', 'region', 'technology', 'product')
-
-        arr_1 = np.array([42])
-        result_1 = DataConverter.to_dataarray(arr_1, coords={'technology': coords_mixed['technology']})
-        assert result_1.shape == (1,)
-
-        arr_5 = np.arange(5)
-        result_5 = DataConverter.to_dataarray(arr_5, coords={'product': coords_mixed['product']})
-        assert result_5.shape == (5,)
-
-        # Arrays with ambiguous length should fail
-        arr_3 = np.array([1, 2, 3])  # matches both scenario and region
-        with pytest.raises(ConversionError, match='matches multiple dimension'):
-            DataConverter.to_dataarray(arr_3, coords=coords_mixed)
-
-    def test_dataframe_with_ambiguous_dimensions(self):
-        """Test DataFrame handling with ambiguous dimensions."""
-        coords_ambiguous = {
-            'scenario': pd.Index(['A', 'B', 'C'], name='scenario'),  # length 3
-            'region': pd.Index(['X', 'Y', 'Z'], name='region'),  # length 3
-        }
-
-        # Multi-column DataFrame with ambiguous dimensions
-        df = pd.DataFrame({'col1': [1, 2, 3], 'col2': [4, 5, 6], 'col3': [7, 8, 9]})  # 3x3 DataFrame
-
-        # Should fail due to ambiguous dimensions
-        with pytest.raises(ConversionError, match='matches multiple dimension'):
-            DataConverter.to_dataarray(df, coords=coords_ambiguous)
-
-    def test_error_message_quality_for_ambiguous_dimensions(self):
-        """Test that error messages for ambiguous dimensions are helpful."""
-        coords_ambiguous = {
-            'scenario': pd.Index(['A', 'B', 'C'], name='scenario'),
-            'region': pd.Index(['north', 'south', 'east'], name='region'),
-            'technology': pd.Index(['solar', 'wind', 'gas'], name='technology'),
-        }
-
-        # 1D array case
-        arr_1d = np.array([1, 2, 3])
-        try:
-            DataConverter.to_dataarray(arr_1d, coords=coords_ambiguous)
-            raise AssertionError('Should have raised ConversionError')
-        except ConversionError as e:
-            error_msg = str(e)
-            assert 'matches multiple dimension' in error_msg
-            assert 'scenario' in error_msg
-            assert 'region' in error_msg
-            assert 'technology' in error_msg
-
-        # 2D array case
-        arr_2d = np.random.rand(3, 3)
-        try:
-            DataConverter.to_dataarray(arr_2d, coords=coords_ambiguous)
-            raise AssertionError('Should have raised ConversionError')
-        except ConversionError as e:
-            error_msg = str(e)
-            assert 'matches multiple dimension combinations' in error_msg
-            assert '(3, 3)' in error_msg
-
-    def test_ambiguous_with_broadcasting_target(self):
-        """Test ambiguous dimensions when target includes broadcasting."""
-        coords_ambiguous_plus = {
-            'time': pd.date_range('2024-01-01', periods=5, freq='D', name='time'),  # length 5
-            'scenario': pd.Index(['A', 'B', 'C'], name='scenario'),  # length 3
-            'region': pd.Index(['X', 'Y', 'Z'], name='region'),  # length 3 - same as scenario
-            'technology': pd.Index(['solar', 'wind'], name='technology'),  # length 2
-        }
-
-        # 1D array with ambiguous length, but targeting broadcast scenario
-        arr_3 = np.array([10, 20, 30])  # length 3, matches scenario and region
-
-        # Should fail even though it would broadcast to other dimensions
-        with pytest.raises(ConversionError, match='matches multiple dimension'):
-            DataConverter.to_dataarray(arr_3, coords=coords_ambiguous_plus)
-
-        # 2D array with one ambiguous dimension
-        arr_5x3 = np.random.rand(5, 3)  # 5 is unique (time), 3 is ambiguous (scenario/region)
-
-        with pytest.raises(ConversionError, match='matches multiple dimension'):
-            DataConverter.to_dataarray(arr_5x3, coords=coords_ambiguous_plus)
-
-    def test_time_dimension_ambiguity(self):
-        """Test ambiguity specifically involving time dimension."""
-        # Create scenario where time has same length as another dimension
-        coords_time_ambiguous = {
-            'time': pd.date_range('2024-01-01', periods=3, freq='D', name='time'),  # length 3
-            'scenario': pd.Index(['base', 'high', 'low'], name='scenario'),  # length 3 - same as time
-            'region': pd.Index(['north', 'south'], name='region'),  # length 2 - unique
-        }
-
-        # Time-indexed series should work even with ambiguous lengths (index matching takes precedence)
-        time_series = pd.Series([100, 200, 300], index=coords_time_ambiguous['time'])
-        result = DataConverter.to_dataarray(time_series, coords=coords_time_ambiguous)
-        assert result.shape == (3, 3, 2)
-        assert result.dims == ('time', 'scenario', 'region')
-
-        # But generic array with length 3 should still fail
-        generic_array = np.array([100, 200, 300])
-        with pytest.raises(ConversionError, match='matches multiple dimension'):
-            DataConverter.to_dataarray(generic_array, coords=coords_time_ambiguous)
-
-
-if __name__ == '__main__':
-    pytest.main()
diff --git a/tests/deprecated/test_flow.py b/tests/deprecated/test_flow.py
index 487a96c25..69922482a 100644
--- a/tests/deprecated/test_flow.py
+++ b/tests/deprecated/test_flow.py
@@ -69,8 +69,9 @@ def test_flow(self, basic_flow_system_linopy_coords, coords_config):
             model.add_variables(lower=10, upper=1000, coords=model.get_coords(['period', 'scenario'])),
         )
 
-        assert flow.relative_minimum.dims == tuple(model.get_coords())
-        assert flow.relative_maximum.dims == tuple(model.get_coords())
+        # Data stays in minimal form (not broadcast to all model dimensions)
+        assert flow.relative_minimum.dims == ('time',)  # Only time dimension
+        assert flow.relative_maximum.dims == ('time',)  # Only time dimension
 
         assert_var_equal(
             flow.submodel.flow_rate,
@@ -182,8 +183,9 @@ def test_flow_invest(self, basic_flow_system_linopy_coords, coords_config):
             model.add_variables(lower=20, upper=100, coords=model.get_coords(['period', 'scenario'])),
         )
 
-        assert flow.relative_minimum.dims == tuple(model.get_coords())
-        assert flow.relative_maximum.dims == tuple(model.get_coords())
+        # Data stays in minimal form (not broadcast to all model dimensions)
+        assert flow.relative_minimum.dims == ('time',)  # Only time dimension
+        assert flow.relative_maximum.dims == ('time',)  # Only time dimension
 
         # flow_rate
         assert_var_equal(
@@ -247,8 +249,9 @@ def test_flow_invest_optional(self, basic_flow_system_linopy_coords, coords_conf
             model.add_variables(binary=True, coords=model.get_coords(['period', 'scenario'])),
         )
 
-        assert flow.relative_minimum.dims == tuple(model.get_coords())
-        assert flow.relative_maximum.dims == tuple(model.get_coords())
+        # Data stays in minimal form (not broadcast to all model dimensions)
+        assert flow.relative_minimum.dims == ('time',)  # Only time dimension
+        assert flow.relative_maximum.dims == ('time',)  # Only time dimension
 
         # flow_rate
         assert_var_equal(
@@ -322,8 +325,9 @@ def test_flow_invest_optional_wo_min_size(self, basic_flow_system_linopy_coords,
             model.add_variables(binary=True, coords=model.get_coords(['period', 'scenario'])),
         )
 
-        assert flow.relative_minimum.dims == tuple(model.get_coords())
-        assert flow.relative_maximum.dims == tuple(model.get_coords())
+        # Data stays in minimal form (not broadcast to all model dimensions)
+        assert flow.relative_minimum.dims == ('time',)  # Only time dimension
+        assert flow.relative_maximum.dims == ('time',)  # Only time dimension
 
         # flow_rate
         assert_var_equal(
@@ -390,8 +394,9 @@ def test_flow_invest_wo_min_size_non_optional(self, basic_flow_system_linopy_coo
             model.add_variables(lower=1e-5, upper=100, coords=model.get_coords(['period', 'scenario'])),
         )
 
-        assert flow.relative_minimum.dims == tuple(model.get_coords())
-        assert flow.relative_maximum.dims == tuple(model.get_coords())
+        # Data stays in minimal form (not broadcast to all model dimensions)
+        assert flow.relative_minimum.dims == ('time',)  # Only time dimension
+        assert flow.relative_maximum.dims == ('time',)  # Only time dimension
 
         # flow_rate
         assert_var_equal(
@@ -629,8 +634,9 @@ def test_effects_per_active_hour(self, basic_flow_system_linopy_coords, coords_c
         costs_per_running_hour = flow.status_parameters.effects_per_active_hour['costs']
         co2_per_running_hour = flow.status_parameters.effects_per_active_hour['CO2']
 
-        assert costs_per_running_hour.dims == tuple(model.get_coords())
-        assert co2_per_running_hour.dims == tuple(model.get_coords())
+        # Data stays in minimal form (1D array stays 1D)
+        assert costs_per_running_hour.dims == ('time',)
+        assert co2_per_running_hour.dims == ('time',)
 
         assert_conequal(
             model.constraints['Sink(Wärme)->costs(temporal)'],
diff --git a/tests/deprecated/test_linear_converter.py b/tests/deprecated/test_linear_converter.py
index d20d104d0..76a45553e 100644
--- a/tests/deprecated/test_linear_converter.py
+++ b/tests/deprecated/test_linear_converter.py
@@ -282,7 +282,8 @@ def test_edge_case_time_varying_conversion(self, basic_flow_system_linopy_coords
 
         factor = converter.conversion_factors[0]['electricity']
 
-        assert factor.dims == tuple(model.get_coords())
+        # Data stays in minimal form (1D array stays 1D)
+        assert factor.dims == ('time',)
 
         # Verify the constraint has the time-varying coefficient
         assert_conequal(
diff --git a/tests/test_component.py b/tests/test_component.py
index 626a64272..c5ebd34a3 100644
--- a/tests/test_component.py
+++ b/tests/test_component.py
@@ -7,6 +7,7 @@
 from .conftest import (
     assert_almost_equal_numeric,
     assert_conequal,
+    assert_dims_compatible,
     assert_sets_equal,
     assert_var_equal,
     create_linopy_model,
@@ -131,7 +132,7 @@ def test_on_with_multiple_flows(self, basic_flow_system_linopy_coords, coords_co
 
         upper_bound_flow_rate = outputs[1].relative_maximum
 
-        assert upper_bound_flow_rate.dims == tuple(model.get_coords())
+        assert_dims_compatible(upper_bound_flow_rate, tuple(model.get_coords()))
 
         assert_var_equal(
             model['TestComponent(Out2)|flow_rate'],
@@ -311,7 +312,7 @@ def test_previous_states_with_multiple_flows(self, basic_flow_system_linopy_coor
 
         upper_bound_flow_rate = outputs[1].relative_maximum
 
-        assert upper_bound_flow_rate.dims == tuple(model.get_coords())
+        assert_dims_compatible(upper_bound_flow_rate, tuple(model.get_coords()))
 
         assert_var_equal(
             model['TestComponent(Out2)|flow_rate'],
diff --git a/tests/test_dataconverter.py b/tests/test_dataconverter.py
index a5774fd6b..f9f2df889 100644
--- a/tests/test_dataconverter.py
+++ b/tests/test_dataconverter.py
@@ -46,34 +46,38 @@ def test_scalar_no_coords(self):
         assert result.item() == 42
 
     def test_scalar_single_coord(self, time_coords):
-        """Scalar with single coordinate should broadcast."""
+        """Scalar with single coordinate stays as scalar (no broadcasting at conversion time)."""
         result = DataConverter.to_dataarray(42, coords={'time': time_coords})
-        assert result.shape == (5,)
-        assert result.dims == ('time',)
-        assert np.all(result.values == 42)
+        # Scalars stay as scalars - broadcasting happens at linopy interface
+        assert result.shape == ()
+        assert result.dims == ()
+        assert result.item() == 42
 
     def test_scalar_multiple_coords(self, time_coords, scenario_coords):
-        """Scalar with multiple coordinates should broadcast to all."""
+        """Scalar with multiple coordinates stays as scalar (no broadcasting at conversion time)."""
         result = DataConverter.to_dataarray(42, coords={'time': time_coords, 'scenario': scenario_coords})
-        assert result.shape == (5, 3)
-        assert result.dims == ('time', 'scenario')
-        assert np.all(result.values == 42)
+        # Scalars stay as scalars - broadcasting happens at linopy interface
+        assert result.shape == ()
+        assert result.dims == ()
+        assert result.item() == 42
 
     def test_numpy_scalars(self, time_coords):
-        """Test numpy scalar types."""
+        """Test numpy scalar types stay as scalars."""
         for scalar in [np.int32(42), np.int64(42), np.float32(42.5), np.float64(42.5)]:
             result = DataConverter.to_dataarray(scalar, coords={'time': time_coords})
-            assert result.shape == (5,)
-            assert np.all(result.values == scalar.item())
+            # Scalars stay as scalars - broadcasting happens at linopy interface
+            assert result.shape == ()
+            assert result.item() == scalar.item()
 
     def test_scalar_many_dimensions(self, standard_coords):
-        """Scalar should broadcast to any number of dimensions."""
+        """Scalar stays as scalar regardless of target dimensions."""
         coords = {**standard_coords, 'technology': pd.Index(['solar', 'wind'], name='technology')}
 
         result = DataConverter.to_dataarray(42, coords=coords)
-        assert result.shape == (5, 3, 2, 2)
-        assert result.dims == ('time', 'scenario', 'region', 'technology')
-        assert np.all(result.values == 42)
+        # Scalars stay as scalars - broadcasting happens at linopy interface
+        assert result.shape == ()
+        assert result.dims == ()
+        assert result.item() == 42
 
 
 class TestOneDimensionalArrayConversion:
@@ -105,26 +109,20 @@ def test_1d_array_mismatched_coord(self, time_coords):
             DataConverter.to_dataarray(arr, coords={'time': time_coords})
 
     def test_1d_array_broadcast_to_multiple_coords(self, time_coords, scenario_coords):
-        """1D array should broadcast to matching dimension."""
-        # Array matching time dimension
+        """1D array stays in minimal form, matching only one dimension (no broadcast at conversion)."""
+        # Array matching time dimension - stays as 1D with time dim only
         time_arr = np.array([10, 20, 30, 40, 50])
         result = DataConverter.to_dataarray(time_arr, coords={'time': time_coords, 'scenario': scenario_coords})
-        assert result.shape == (5, 3)
-        assert result.dims == ('time', 'scenario')
-
-        # Each scenario should have the same time values
-        for scenario in scenario_coords:
-            assert np.array_equal(result.sel(scenario=scenario).values, time_arr)
+        assert result.shape == (5,)
+        assert result.dims == ('time',)
+        assert np.array_equal(result.values, time_arr)
 
-        # Array matching scenario dimension
+        # Array matching scenario dimension - stays as 1D with scenario dim only
         scenario_arr = np.array([100, 200, 300])
         result = DataConverter.to_dataarray(scenario_arr, coords={'time': time_coords, 'scenario': scenario_coords})
-        assert result.shape == (5, 3)
-        assert result.dims == ('time', 'scenario')
-
-        # Each time should have the same scenario values
-        for time in time_coords:
-            assert np.array_equal(result.sel(time=time).values, scenario_arr)
+        assert result.shape == (3,)
+        assert result.dims == ('scenario',)
+        assert np.array_equal(result.values, scenario_arr)
 
     def test_1d_array_ambiguous_length(self):
         """Array length matching multiple dimensions should fail."""
@@ -139,18 +137,14 @@ def test_1d_array_ambiguous_length(self):
             DataConverter.to_dataarray(arr, coords=coords_3x3)
 
     def test_1d_array_broadcast_to_many_dimensions(self, standard_coords):
-        """1D array should broadcast to many dimensions."""
-        # Array matching time dimension
+        """1D array stays in minimal form with matched dimension only (no broadcast at conversion)."""
+        # Array matching time dimension - stays as 1D with time dim only
         time_arr = np.array([10, 20, 30, 40, 50])
         result = DataConverter.to_dataarray(time_arr, coords=standard_coords)
 
-        assert result.shape == (5, 3, 2)
-        assert result.dims == ('time', 'scenario', 'region')
-
-        # Check broadcasting - all scenarios and regions should have same time values
-        for scenario in standard_coords['scenario']:
-            for region in standard_coords['region']:
-                assert np.array_equal(result.sel(scenario=scenario, region=region).values, time_arr)
+        assert result.shape == (5,)
+        assert result.dims == ('time',)
+        assert np.array_equal(result.values, time_arr)
 
 
 class TestSeriesConversion:
@@ -194,14 +188,13 @@ def test_series_mismatched_index(self, time_coords):
             DataConverter.to_dataarray(series, coords={'time': time_coords})
 
     def test_series_broadcast_to_multiple_coords(self, time_coords, scenario_coords):
-        """Series should broadcast to non-matching dimensions."""
-        # Time series broadcast to scenarios
+        """Series stays in minimal form with matched dimension only (no broadcast at conversion)."""
+        # Time series stays as 1D with time dim only
         time_series = pd.Series([10, 20, 30, 40, 50], index=time_coords)
         result = DataConverter.to_dataarray(time_series, coords={'time': time_coords, 'scenario': scenario_coords})
-        assert result.shape == (5, 3)
-
-        for scenario in scenario_coords:
-            assert np.array_equal(result.sel(scenario=scenario).values, time_series.values)
+        assert result.shape == (5,)
+        assert result.dims == ('time',)
+        assert np.array_equal(result.values, time_series.values)
 
     def test_series_wrong_dimension(self, time_coords, region_coords):
         """Series indexed by dimension not in coords should fail."""
@@ -211,17 +204,13 @@ def test_series_wrong_dimension(self, time_coords, region_coords):
             DataConverter.to_dataarray(wrong_series, coords={'time': time_coords})
 
     def test_series_broadcast_to_many_dimensions(self, standard_coords):
-        """Series should broadcast to many dimensions."""
+        """Series stays in minimal form with matched dimension only (no broadcast at conversion)."""
         time_series = pd.Series([100, 200, 300, 400, 500], index=standard_coords['time'])
         result = DataConverter.to_dataarray(time_series, coords=standard_coords)
 
-        assert result.shape == (5, 3, 2)
-        assert result.dims == ('time', 'scenario', 'region')
-
-        # Check that all non-time dimensions have the same time series values
-        for scenario in standard_coords['scenario']:
-            for region in standard_coords['region']:
-                assert np.array_equal(result.sel(scenario=scenario, region=region).values, time_series.values)
+        assert result.shape == (5,)
+        assert result.dims == ('time',)
+        assert np.array_equal(result.values, time_series.values)
 
 
 class TestDataFrameConversion:
@@ -258,13 +247,14 @@ def test_empty_dataframe_rejected(self, time_coords):
             DataConverter.to_dataarray(df, coords={'time': time_coords})
 
     def test_dataframe_broadcast(self, time_coords, scenario_coords):
-        """Single-column DataFrame should broadcast like Series."""
+        """Single-column DataFrame stays 1D (no broadcasting at conversion time)."""
         df = pd.DataFrame({'power': [10, 20, 30, 40, 50]}, index=time_coords)
         result = DataConverter.to_dataarray(df, coords={'time': time_coords, 'scenario': scenario_coords})
 
-        assert result.shape == (5, 3)
-        for scenario in scenario_coords:
-            assert np.array_equal(result.sel(scenario=scenario).values, df['power'].values)
+        # Data stays in minimal form - 1D along time
+        assert result.shape == (5,)
+        assert result.dims == ('time',)
+        assert np.array_equal(result.values, df['power'].values)
 
 
 class TestMultiDimensionalArrayConversion:
@@ -282,7 +272,7 @@ def test_2d_array_unique_dimensions(self, standard_coords):
         assert result.dims == ('time', 'scenario')
         assert np.array_equal(result.values, data_2d)
 
-        # 3x5 array should map to scenario x time
+        # 3x5 array should map to scenario x time, then transpose to canonical (time, scenario) order
         data_2d_flipped = np.random.rand(3, 5)
         result_flipped = DataConverter.to_dataarray(
             data_2d_flipped, coords={'time': standard_coords['time'], 'scenario': standard_coords['scenario']}
@@ -292,18 +282,15 @@ def test_2d_array_unique_dimensions(self, standard_coords):
         assert result_flipped.dims == ('time', 'scenario')
         assert np.array_equal(result_flipped.values.transpose(), data_2d_flipped)
 
-    def test_2d_array_broadcast_to_3d(self, standard_coords):
-        """2D array should broadcast to additional dimensions when using partial matching."""
-        # With improved integration, 2D array (5x3) should match time×scenario and broadcast to region
+    def test_2d_array_stays_2d(self, standard_coords):
+        """2D array stays 2D even when more coords are provided (no broadcasting)."""
+        # 2D array (5x3) matches time×scenario, stays 2D (no broadcast to region)
         data_2d = np.random.rand(5, 3)
         result = DataConverter.to_dataarray(data_2d, coords=standard_coords)
 
-        assert result.shape == (5, 3, 2)
-        assert result.dims == ('time', 'scenario', 'region')
-
-        # Check that all regions have the same time x scenario data
-        for region in standard_coords['region']:
-            assert np.array_equal(result.sel(region=region).values, data_2d)
+        assert result.shape == (5, 3)
+        assert result.dims == ('time', 'scenario')
+        assert np.array_equal(result.values, data_2d)
 
     def test_3d_array_unique_dimensions(self, standard_coords):
         """3D array with unique dimension lengths should work."""
@@ -316,8 +303,8 @@ def test_3d_array_unique_dimensions(self, standard_coords):
         assert np.array_equal(result.values, data_3d)
 
     def test_3d_array_different_permutation(self, standard_coords):
-        """3D array with different dimension order should work."""
-        # 2x5x3 array should map to region x time x scenario
+        """3D array with different dimension order should be transposed to canonical order."""
+        # 2x5x3 array should map to region x time x scenario, then transpose to canonical order
         data_3d = np.random.rand(2, 5, 3)
         result = DataConverter.to_dataarray(data_3d, coords=standard_coords)
 
@@ -450,28 +437,26 @@ def test_incompatible_dataarray_dims(self, time_coords):
         with pytest.raises(ConversionError):
             DataConverter.to_dataarray(original, coords={'time': time_coords})
 
-    def test_dataarray_broadcast(self, time_coords, scenario_coords):
-        """DataArray should broadcast to additional dimensions."""
-        # 1D time DataArray to 2D time+scenario
+    def test_dataarray_no_broadcast(self, time_coords, scenario_coords):
+        """DataArray stays in minimal form (no broadcasting at conversion time)."""
+        # 1D time DataArray stays 1D even with additional coords
         original = xr.DataArray([10, 20, 30, 40, 50], coords={'time': time_coords}, dims='time')
         result = DataConverter.to_dataarray(original, coords={'time': time_coords, 'scenario': scenario_coords})
 
-        assert result.shape == (5, 3)
-        assert result.dims == ('time', 'scenario')
-
-        for scenario in scenario_coords:
-            assert np.array_equal(result.sel(scenario=scenario).values, original.values)
+        assert result.shape == (5,)
+        assert result.dims == ('time',)
+        assert np.array_equal(result.values, original.values)
 
-    def test_scalar_dataarray_broadcast(self, time_coords, scenario_coords):
-        """Scalar DataArray should broadcast to all dimensions."""
+    def test_scalar_dataarray_stays_scalar(self, time_coords, scenario_coords):
+        """Scalar DataArray stays scalar (no broadcasting at conversion time)."""
         scalar_da = xr.DataArray(42)
         result = DataConverter.to_dataarray(scalar_da, coords={'time': time_coords, 'scenario': scenario_coords})
 
-        assert result.shape == (5, 3)
-        assert np.all(result.values == 42)
+        assert result.shape == ()
+        assert result.item() == 42
 
-    def test_2d_dataarray_broadcast_to_more_dimensions(self, standard_coords):
-        """DataArray should broadcast to additional dimensions."""
+    def test_2d_dataarray_stays_2d(self, standard_coords):
+        """DataArray stays in minimal form (no broadcasting at conversion time)."""
         # Start with 2D DataArray
         original = xr.DataArray(
             [[10, 20, 30], [40, 50, 60], [70, 80, 90], [100, 110, 120], [130, 140, 150]],
@@ -479,15 +464,12 @@ def test_2d_dataarray_broadcast_to_more_dimensions(self, standard_coords):
             dims=('time', 'scenario'),
         )
 
-        # Broadcast to 3D
+        # Stays 2D (no broadcast to 3D)
         result = DataConverter.to_dataarray(original, coords=standard_coords)
 
-        assert result.shape == (5, 3, 2)
-        assert result.dims == ('time', 'scenario', 'region')
-
-        # Check that all regions have the same time+scenario values
-        for region in standard_coords['region']:
-            assert np.array_equal(result.sel(region=region).values, original.values)
+        assert result.shape == (5, 3)
+        assert result.dims == ('time', 'scenario')
+        assert np.array_equal(result.values, original.values)
 
 
 class TestTimeSeriesDataConversion:
@@ -504,16 +486,17 @@ def test_timeseries_data_basic(self, time_coords):
         assert result.dims == ('time',)
         assert np.array_equal(result.values, [10, 20, 30, 40, 50])
 
-    def test_timeseries_data_broadcast(self, time_coords, scenario_coords):
-        """TimeSeriesData should broadcast to additional dimensions."""
+    def test_timeseries_data_stays_minimal(self, time_coords, scenario_coords):
+        """TimeSeriesData stays in minimal form (no broadcasting at conversion time)."""
         data_array = xr.DataArray([10, 20, 30, 40, 50], coords={'time': time_coords}, dims='time')
         ts_data = TimeSeriesData(data_array)
 
         result = DataConverter.to_dataarray(ts_data, coords={'time': time_coords, 'scenario': scenario_coords})
 
-        assert result.shape == (5, 3)
-        for scenario in scenario_coords:
-            assert np.array_equal(result.sel(scenario=scenario).values, [10, 20, 30, 40, 50])
+        # Stays 1D (time only)
+        assert result.shape == (5,)
+        assert result.dims == ('time',)
+        assert np.array_equal(result.values, [10, 20, 30, 40, 50])
 
 
 class TestAsDataArrayAlias:
@@ -547,60 +530,52 @@ class TestCustomDimensions:
     """Test with custom dimension names beyond time/scenario."""
 
     def test_custom_single_dimension(self, region_coords):
-        """Test with custom dimension name."""
+        """Scalar stays scalar even with custom dimension."""
         result = DataConverter.to_dataarray(42, coords={'region': region_coords})
-        assert result.shape == (3,)
-        assert result.dims == ('region',)
-        assert np.all(result.values == 42)
+        assert result.shape == ()
+        assert result.dims == ()
+        assert result.item() == 42
 
     def test_custom_multiple_dimensions(self):
-        """Test with multiple custom dimensions."""
+        """Array with matching dimension stays 1D (no broadcasting)."""
         products = pd.Index(['A', 'B'], name='product')
         technologies = pd.Index(['solar', 'wind', 'gas'], name='technology')
 
-        # Array matching technology dimension
+        # Array matching technology dimension stays 1D
         arr = np.array([100, 150, 80])
         result = DataConverter.to_dataarray(arr, coords={'product': products, 'technology': technologies})
 
-        assert result.shape == (2, 3)
-        assert result.dims == ('product', 'technology')
-
-        # Should broadcast across products
-        for product in products:
-            assert np.array_equal(result.sel(product=product).values, arr)
+        assert result.shape == (3,)
+        assert result.dims == ('technology',)
+        assert np.array_equal(result.values, arr)
 
     def test_mixed_dimension_types(self):
-        """Test mixing time dimension with custom dimensions."""
+        """Time series stays 1D (no broadcasting to regions)."""
         time_coords = pd.date_range('2024-01-01', periods=3, freq='D', name='time')
         regions = pd.Index(['north', 'south'], name='region')
 
-        # Time series should broadcast to regions
         time_series = pd.Series([10, 20, 30], index=time_coords)
         result = DataConverter.to_dataarray(time_series, coords={'time': time_coords, 'region': regions})
 
-        assert result.shape == (3, 2)
-        assert result.dims == ('time', 'region')
+        # Stays 1D along time
+        assert result.shape == (3,)
+        assert result.dims == ('time',)
 
     def test_custom_dimensions_complex(self):
-        """Test complex scenario with custom dimensions."""
+        """Array with matching dimension stays 1D (no broadcasting)."""
         coords = {
             'product': pd.Index(['A', 'B'], name='product'),
             'factory': pd.Index(['F1', 'F2', 'F3'], name='factory'),
             'quarter': pd.Index(['Q1', 'Q2', 'Q3', 'Q4'], name='quarter'),
         }
 
-        # Array matching factory dimension
+        # Array matching factory dimension stays 1D
         factory_arr = np.array([100, 200, 300])
         result = DataConverter.to_dataarray(factory_arr, coords=coords)
 
-        assert result.shape == (2, 3, 4)
-        assert result.dims == ('product', 'factory', 'quarter')
-
-        # Check broadcasting
-        for product in coords['product']:
-            for quarter in coords['quarter']:
-                slice_data = result.sel(product=product, quarter=quarter)
-                assert np.array_equal(slice_data.values, factory_arr)
+        assert result.shape == (3,)
+        assert result.dims == ('factory',)
+        assert np.array_equal(result.values, factory_arr)
 
 
 class TestValidation:
@@ -631,11 +606,11 @@ def test_time_coord_validation(self):
             DataConverter.to_dataarray(42, coords={'time': wrong_time})
 
     def test_coord_naming(self, time_coords):
-        """Coordinates should be auto-renamed to match dimension."""
-        # Unnamed time index should be renamed
-        unnamed_time = time_coords.rename(None)
-        result = DataConverter.to_dataarray(42, coords={'time': unnamed_time})
-        assert result.coords['time'].name == 'time'
+        """Scalar with coords stays scalar but validates coords."""
+        # Scalar stays scalar regardless of coords
+        result = DataConverter.to_dataarray(42, coords={'time': time_coords})
+        assert result.shape == ()
+        assert result.item() == 42
 
 
 class TestErrorHandling:
@@ -735,28 +710,28 @@ class TestBooleanValues:
     """Test handling of boolean values and arrays."""
 
     def test_scalar_boolean_to_dataarray(self, time_coords):
-        """Scalar boolean values should work with to_dataarray."""
+        """Scalar boolean values stay scalar (no broadcasting)."""
         result_true = DataConverter.to_dataarray(True, coords={'time': time_coords})
-        assert result_true.shape == (5,)
+        assert result_true.shape == ()
         assert result_true.dtype == bool
-        assert np.all(result_true.values)
+        assert result_true.item() is True
 
         result_false = DataConverter.to_dataarray(False, coords={'time': time_coords})
-        assert result_false.shape == (5,)
+        assert result_false.shape == ()
         assert result_false.dtype == bool
-        assert not np.any(result_false.values)
+        assert result_false.item() is False
 
     def test_numpy_boolean_scalar(self, time_coords):
-        """Numpy boolean scalars should work."""
+        """Numpy boolean scalars stay scalar (no broadcasting)."""
         result_np_true = DataConverter.to_dataarray(np.bool_(True), coords={'time': time_coords})
-        assert result_np_true.shape == (5,)
+        assert result_np_true.shape == ()
         assert result_np_true.dtype == bool
-        assert np.all(result_np_true.values)
+        assert result_np_true.item() is True
 
         result_np_false = DataConverter.to_dataarray(np.bool_(False), coords={'time': time_coords})
-        assert result_np_false.shape == (5,)
+        assert result_np_false.shape == ()
         assert result_np_false.dtype == bool
-        assert not np.any(result_np_false.values)
+        assert result_np_false.item() is False
 
     def test_boolean_array_to_dataarray(self, time_coords):
         """Boolean arrays should work with to_dataarray."""
@@ -779,19 +754,17 @@ def test_boolean_no_coords(self):
         assert result_as.dims == ()
         assert not result_as.item()
 
-    def test_boolean_multidimensional_broadcast(self, standard_coords):
-        """Boolean values should broadcast to multiple dimensions."""
+    def test_boolean_scalar_stays_scalar(self, standard_coords):
+        """Boolean scalars stay scalar (no broadcasting)."""
         result = DataConverter.to_dataarray(True, coords=standard_coords)
-        assert result.shape == (5, 3, 2)
-        assert result.dims == ('time', 'scenario', 'region')
+        assert result.shape == ()
         assert result.dtype == bool
-        assert np.all(result.values)
+        assert result.item() is True
 
         result_as = DataConverter.to_dataarray(False, coords=standard_coords)
-        assert result_as.shape == (5, 3, 2)
-        assert result_as.dims == ('time', 'scenario', 'region')
+        assert result_as.shape == ()
         assert result_as.dtype == bool
-        assert not np.any(result_as.values)
+        assert result_as.item() is False
 
     def test_boolean_series(self, time_coords):
         """Boolean Series should work."""
@@ -873,60 +846,51 @@ def test_mixed_numeric_types(self, time_coords):
         assert np.issubdtype(result.dtype, np.floating)
         assert np.array_equal(result.values, mixed_arr)
 
-    def test_special_values_in_multid_arrays(self, standard_coords):
-        """Special values should be preserved in multi-D arrays and broadcasting."""
+    def test_special_values_in_1d_arrays(self, standard_coords):
+        """Special values should be preserved in arrays (no broadcasting)."""
         # Array with NaN and inf
         special_arr = np.array([1, np.nan, np.inf, -np.inf, 5])
         result = DataConverter.to_dataarray(special_arr, coords=standard_coords)
 
-        assert result.shape == (5, 3, 2)
+        # Stays 1D (no broadcasting)
+        assert result.shape == (5,)
+        assert result.dims == ('time',)
 
-        # Check that special values are preserved in all broadcasts
-        for scenario in standard_coords['scenario']:
-            for region in standard_coords['region']:
-                slice_data = result.sel(scenario=scenario, region=region)
-                assert np.array_equal(np.isnan(slice_data.values), np.isnan(special_arr))
-                assert np.array_equal(np.isinf(slice_data.values), np.isinf(special_arr))
+        # Check that special values are preserved
+        assert np.array_equal(np.isnan(result.values), np.isnan(special_arr))
+        assert np.array_equal(np.isinf(result.values), np.isinf(special_arr))
 
 
 class TestAdvancedBroadcasting:
-    """Test advanced broadcasting scenarios and edge cases."""
+    """Test advanced scenarios (no broadcasting - data stays minimal)."""
 
-    def test_partial_dimension_matching_with_broadcasting(self, standard_coords):
-        """Test that partial dimension matching works with the improved integration."""
-        # 1D array matching one dimension should broadcast to all target dimensions
+    def test_partial_dimension_matching_stays_1d(self, standard_coords):
+        """1D array stays 1D (no broadcasting to additional dimensions)."""
         time_arr = np.array([10, 20, 30, 40, 50])  # matches time (length 5)
         result = DataConverter.to_dataarray(time_arr, coords=standard_coords)
 
-        assert result.shape == (5, 3, 2)
-        assert result.dims == ('time', 'scenario', 'region')
-
-        # Verify broadcasting
-        for scenario in standard_coords['scenario']:
-            for region in standard_coords['region']:
-                assert np.array_equal(result.sel(scenario=scenario, region=region).values, time_arr)
+        # Stays 1D
+        assert result.shape == (5,)
+        assert result.dims == ('time',)
+        assert np.array_equal(result.values, time_arr)
 
     def test_complex_multid_scenario(self):
-        """Complex real-world scenario with multi-D array and broadcasting."""
-        # Energy system data: time x technology, broadcast to regions
+        """2D array stays 2D (no broadcasting to additional dimensions)."""
         coords = {
             'time': pd.date_range('2024-01-01', periods=24, freq='h', name='time'),  # 24 hours
             'technology': pd.Index(['solar', 'wind', 'gas', 'coal'], name='technology'),  # 4 technologies
             'region': pd.Index(['north', 'south', 'east'], name='region'),  # 3 regions
         }
 
-        # Capacity factors: 24 x 4 (will broadcast to 24 x 4 x 3)
+        # Capacity factors: 24 x 4 stays 2D (no broadcast to 24 x 4 x 3)
         capacity_factors = np.random.rand(24, 4)
 
         result = DataConverter.to_dataarray(capacity_factors, coords=coords)
 
-        assert result.shape == (24, 4, 3)
-        assert result.dims == ('time', 'technology', 'region')
+        assert result.shape == (24, 4)
+        assert result.dims == ('time', 'technology')
         assert isinstance(result.indexes['time'], pd.DatetimeIndex)
-
-        # Verify broadcasting: all regions should have same time×technology data
-        for region in coords['region']:
-            assert np.array_equal(result.sel(region=region).values, capacity_factors)
+        assert np.array_equal(result.values, capacity_factors)
 
     def test_ambiguous_length_handling(self):
         """Test handling of ambiguous length scenarios across different data types."""
@@ -952,8 +916,8 @@ def test_ambiguous_length_handling(self):
         with pytest.raises(ConversionError, match='matches multiple dimension'):
             DataConverter.to_dataarray(arr_3d, coords=coords_3x3x3)
 
-    def test_mixed_broadcasting_scenarios(self):
-        """Test various broadcasting scenarios with different input types."""
+    def test_no_broadcasting_scenarios(self):
+        """Data stays in minimal form (no broadcasting to additional dimensions)."""
         coords = {
             'time': pd.date_range('2024-01-01', periods=4, freq='D', name='time'),  # length 4
             'scenario': pd.Index(['A', 'B'], name='scenario'),  # length 2
@@ -961,31 +925,24 @@ def test_mixed_broadcasting_scenarios(self):
             'product': pd.Index(['X', 'Y', 'Z', 'W', 'V'], name='product'),  # length 5
         }
 
-        # Scalar to 4D
+        # Scalar stays scalar
         scalar_result = DataConverter.to_dataarray(42, coords=coords)
-        assert scalar_result.shape == (4, 2, 3, 5)
-        assert np.all(scalar_result.values == 42)
+        assert scalar_result.shape == ()
+        assert scalar_result.item() == 42
 
-        # 1D array (length 4, matches time) to 4D
+        # 1D array (length 4, matches time) stays 1D
         arr_1d = np.array([10, 20, 30, 40])
         arr_result = DataConverter.to_dataarray(arr_1d, coords=coords)
-        assert arr_result.shape == (4, 2, 3, 5)
-        # Verify broadcasting
-        for scenario in coords['scenario']:
-            for region in coords['region']:
-                for product in coords['product']:
-                    assert np.array_equal(
-                        arr_result.sel(scenario=scenario, region=region, product=product).values, arr_1d
-                    )
-
-        # 2D array (4x2, matches time×scenario) to 4D
+        assert arr_result.shape == (4,)
+        assert arr_result.dims == ('time',)
+        assert np.array_equal(arr_result.values, arr_1d)
+
+        # 2D array (4x2, matches time×scenario) stays 2D
         arr_2d = np.random.rand(4, 2)
         arr_2d_result = DataConverter.to_dataarray(arr_2d, coords=coords)
-        assert arr_2d_result.shape == (4, 2, 3, 5)
-        # Verify broadcasting
-        for region in coords['region']:
-            for product in coords['product']:
-                assert np.array_equal(arr_2d_result.sel(region=region, product=product).values, arr_2d)
+        assert arr_2d_result.shape == (4, 2)
+        assert arr_2d_result.dims == ('time', 'scenario')
+        assert np.array_equal(arr_2d_result.values, arr_2d)
 
 
 class TestAmbiguousDimensionLengthHandling:
@@ -1020,11 +977,11 @@ def test_1d_array_ambiguous_dimensions_complex(self):
         with pytest.raises(ConversionError, match='matches multiple dimension'):
             DataConverter.to_dataarray(arr_1d, coords=coords_4x4x4)
 
-        # Array matching the unique length should work
+        # Array matching the unique length should work (stays 1D, no broadcasting)
         arr_1d_unique = np.array([100, 200])  # length 2 - matches only product
         result = DataConverter.to_dataarray(arr_1d_unique, coords=coords_4x4x4)
-        assert result.shape == (4, 4, 4, 2)  # broadcast to all dimensions
-        assert result.dims == ('time', 'scenario', 'region', 'product')
+        assert result.shape == (2,)  # stays 1D
+        assert result.dims == ('product',)
 
     def test_2d_array_ambiguous_dimensions_both_same(self):
         """Test 2D array where both dimensions have the same ambiguous length."""
@@ -1111,11 +1068,11 @@ def test_pandas_series_ambiguous_dimensions(self):
         with pytest.raises(ConversionError, match='Series index does not match any target dimension coordinates'):
             DataConverter.to_dataarray(generic_series, coords=coords_ambiguous)
 
-        # Series with index that matches one of the ambiguous coordinates should work
+        # Series with index that matches one of the ambiguous coordinates should work (stays 1D)
         scenario_series = pd.Series([10, 20, 30], index=coords_ambiguous['scenario'])
         result = DataConverter.to_dataarray(scenario_series, coords=coords_ambiguous)
-        assert result.shape == (3, 3)  # should broadcast to both dimensions
-        assert result.dims == ('scenario', 'region')
+        assert result.shape == (3,)  # stays 1D
+        assert result.dims == ('scenario',)
 
     def test_edge_case_many_same_lengths(self):
         """Test edge case with many dimensions having the same length."""
@@ -1153,10 +1110,10 @@ def test_mixed_lengths_with_duplicates(self):
             'product': pd.Index(['P1', 'P2', 'P3', 'P4', 'P5'], name='product'),  # length 5 - unique
         }
 
-        # Arrays with unique lengths should work
+        # Arrays with unique lengths should work (stays minimal, no broadcasting)
         arr_8 = np.arange(8)
         result_8 = DataConverter.to_dataarray(arr_8, coords=coords_mixed)
-        assert result_8.dims == ('time', 'scenario', 'region', 'technology', 'product')
+        assert result_8.dims == ('time',)  # Stays 1D, matches time dimension
 
         arr_1 = np.array([42])
         result_1 = DataConverter.to_dataarray(arr_1, coords={'technology': coords_mixed['technology']})
@@ -1247,10 +1204,11 @@ def test_time_dimension_ambiguity(self):
         }
 
         # Time-indexed series should work even with ambiguous lengths (index matching takes precedence)
+        # Stays minimal - no broadcasting to other dimensions
         time_series = pd.Series([100, 200, 300], index=coords_time_ambiguous['time'])
         result = DataConverter.to_dataarray(time_series, coords=coords_time_ambiguous)
-        assert result.shape == (3, 3, 2)
-        assert result.dims == ('time', 'scenario', 'region')
+        assert result.shape == (3,)  # Stays 1D
+        assert result.dims == ('time',)  # Matches time via index
 
         # But generic array with length 3 should still fail
         generic_array = np.array([100, 200, 300])
diff --git a/tests/test_flow.py b/tests/test_flow.py
index 275963348..aa75b3c66 100644
--- a/tests/test_flow.py
+++ b/tests/test_flow.py
@@ -4,7 +4,7 @@
 
 import flixopt as fx
 
-from .conftest import assert_conequal, assert_sets_equal, assert_var_equal, create_linopy_model
+from .conftest import assert_conequal, assert_dims_compatible, assert_sets_equal, assert_var_equal, create_linopy_model
 
 
 class TestFlowModel:
@@ -69,8 +69,8 @@ def test_flow(self, basic_flow_system_linopy_coords, coords_config):
             model.add_variables(lower=10, upper=1000, coords=model.get_coords(['period', 'scenario'])),
         )
 
-        assert flow.relative_minimum.dims == tuple(model.get_coords())
-        assert flow.relative_maximum.dims == tuple(model.get_coords())
+        assert_dims_compatible(flow.relative_minimum, tuple(model.get_coords()))
+        assert_dims_compatible(flow.relative_maximum, tuple(model.get_coords()))
 
         assert_var_equal(
             flow.submodel.flow_rate,
@@ -182,8 +182,8 @@ def test_flow_invest(self, basic_flow_system_linopy_coords, coords_config):
             model.add_variables(lower=20, upper=100, coords=model.get_coords(['period', 'scenario'])),
         )
 
-        assert flow.relative_minimum.dims == tuple(model.get_coords())
-        assert flow.relative_maximum.dims == tuple(model.get_coords())
+        assert_dims_compatible(flow.relative_minimum, tuple(model.get_coords()))
+        assert_dims_compatible(flow.relative_maximum, tuple(model.get_coords()))
 
         # flow_rate
         assert_var_equal(
@@ -247,8 +247,8 @@ def test_flow_invest_optional(self, basic_flow_system_linopy_coords, coords_conf
             model.add_variables(binary=True, coords=model.get_coords(['period', 'scenario'])),
         )
 
-        assert flow.relative_minimum.dims == tuple(model.get_coords())
-        assert flow.relative_maximum.dims == tuple(model.get_coords())
+        assert_dims_compatible(flow.relative_minimum, tuple(model.get_coords()))
+        assert_dims_compatible(flow.relative_maximum, tuple(model.get_coords()))
 
         # flow_rate
         assert_var_equal(
@@ -322,8 +322,8 @@ def test_flow_invest_optional_wo_min_size(self, basic_flow_system_linopy_coords,
             model.add_variables(binary=True, coords=model.get_coords(['period', 'scenario'])),
         )
 
-        assert flow.relative_minimum.dims == tuple(model.get_coords())
-        assert flow.relative_maximum.dims == tuple(model.get_coords())
+        assert_dims_compatible(flow.relative_minimum, tuple(model.get_coords()))
+        assert_dims_compatible(flow.relative_maximum, tuple(model.get_coords()))
 
         # flow_rate
         assert_var_equal(
@@ -390,8 +390,8 @@ def test_flow_invest_wo_min_size_non_optional(self, basic_flow_system_linopy_coo
             model.add_variables(lower=1e-5, upper=100, coords=model.get_coords(['period', 'scenario'])),
         )
 
-        assert flow.relative_minimum.dims == tuple(model.get_coords())
-        assert flow.relative_maximum.dims == tuple(model.get_coords())
+        assert_dims_compatible(flow.relative_minimum, tuple(model.get_coords()))
+        assert_dims_compatible(flow.relative_maximum, tuple(model.get_coords()))
 
         # flow_rate
         assert_var_equal(
@@ -629,8 +629,8 @@ def test_effects_per_active_hour(self, basic_flow_system_linopy_coords, coords_c
         costs_per_running_hour = flow.status_parameters.effects_per_active_hour['costs']
         co2_per_running_hour = flow.status_parameters.effects_per_active_hour['CO2']
 
-        assert costs_per_running_hour.dims == tuple(model.get_coords())
-        assert co2_per_running_hour.dims == tuple(model.get_coords())
+        assert_dims_compatible(costs_per_running_hour, tuple(model.get_coords()))
+        assert_dims_compatible(co2_per_running_hour, tuple(model.get_coords()))
 
         assert_conequal(
             model.constraints['Sink(Wärme)->costs(temporal)'],
diff --git a/tests/test_io_conversion.py b/tests/test_io_conversion.py
index 7775f987a..dffba1dfc 100644
--- a/tests/test_io_conversion.py
+++ b/tests/test_io_conversion.py
@@ -204,12 +204,12 @@ def test_custom_renames(self):
         assert 'custom_new' in result.attrs
         assert 'custom_old' not in result.attrs
 
-    def test_returns_same_object(self):
-        """Test that the function modifies and returns the same dataset object."""
+    def test_returns_equivalent_dataset(self):
+        """Test that the function converts and returns equivalent dataset."""
         ds = xr.Dataset(attrs={'minimum_operation': 100})
         result = convert_old_dataset(ds)
-        # Note: attrs are modified in place, so the object should be the same
-        assert result is ds
+        # Check that attrs are converted
+        assert result.attrs == {'minimum_temporal': 100}
 
 
 class TestConvertOldNetcdf:
diff --git a/tests/test_linear_converter.py b/tests/test_linear_converter.py
index d20d104d0..c8fc3fb52 100644
--- a/tests/test_linear_converter.py
+++ b/tests/test_linear_converter.py
@@ -4,7 +4,7 @@
 
 import flixopt as fx
 
-from .conftest import assert_conequal, assert_var_equal, create_linopy_model
+from .conftest import assert_conequal, assert_dims_compatible, assert_var_equal, create_linopy_model
 
 
 class TestLinearConverterModel:
@@ -282,7 +282,7 @@ def test_edge_case_time_varying_conversion(self, basic_flow_system_linopy_coords
 
         factor = converter.conversion_factors[0]['electricity']
 
-        assert factor.dims == tuple(model.get_coords())
+        assert_dims_compatible(factor, tuple(model.get_coords()))
 
         # Verify the constraint has the time-varying coefficient
         assert_conequal(