diff --git a/CHANGELOG.md b/CHANGELOG.md
index 20f2de7d7..f60497d80 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -53,7 +53,7 @@ Until here -->
 
 ## [5.1.0] - Upcoming
 
-**Summary**: Time-series clustering for faster optimization with configurable storage behavior across typical periods.
+**Summary**: Time-series clustering for faster optimization with configurable storage behavior across typical periods. Improved weights API with always-normalized scenario weights.
 
 ### ✨ Added
 
@@ -121,6 +121,23 @@ charge_state = fs_expanded.solution['SeasonalPit|charge_state']
     Use `'cyclic'` for short-term storage like batteries or hot water tanks where only daily patterns matter.
     Use `'independent'` for quick estimates when storage behavior isn't critical.
 
+### 💥 Breaking Changes
+
+- `FlowSystem.scenario_weights` are now always normalized to sum to 1 when set (including after `.sel()` subsetting)
+
+### ♻️ Changed
+
+- `FlowSystem.weights` returns `dict[str, xr.DataArray]` (unit weights instead of `1.0` float fallback)
+- `FlowSystemDimensions` type now includes `'cluster'`
+
+### 🗑️ Deprecated
+
+- `normalize_weights` parameter in `create_model()`, `build_model()`, `optimize()`
+
+### 🐛 Fixed
+
+- `temporal_weight` and `sum_temporal()` now use consistent implementation
+
 ### 👷 Development
 
 **New Test Suites for Clustering**:
diff --git a/flixopt/__init__.py b/flixopt/__init__.py
index 73784f2cd..474668b6f 100644
--- a/flixopt/__init__.py
+++ b/flixopt/__init__.py
@@ -31,7 +31,6 @@
 from .interface import InvestParameters, Piece, Piecewise, PiecewiseConversion, PiecewiseEffects, StatusParameters
 from .optimization import Optimization, SegmentedOptimization
 from .plot_result import PlotResult
-from .structure import TimeSeriesWeights
 
 __all__ = [
     'TimeSeriesData',
@@ -58,7 +57,6 @@
     'PiecewiseConversion',
     'PiecewiseEffects',
     'PlotResult',
-    'TimeSeriesWeights',
     'clustering',
     'plotting',
     'results',
diff --git a/flixopt/clustering/base.py b/flixopt/clustering/base.py
index 2c442e3d5..fadf28247 100644
--- a/flixopt/clustering/base.py
+++ b/flixopt/clustering/base.py
@@ -21,11 +21,11 @@
 from typing import TYPE_CHECKING, Any
 
 import numpy as np
+import pandas as pd
 import xarray as xr
 
 if TYPE_CHECKING:
     from ..color_processing import ColorType
-    from ..flow_system import FlowSystem
     from ..plot_result import PlotResult
     from ..statistics_accessor import SelectType
 
@@ -98,6 +98,31 @@ def __repr__(self) -> str:
             f')'
         )
 
+    def _create_reference_structure(self) -> tuple[dict, dict[str, xr.DataArray]]:
+        """Create reference structure for serialization."""
+        ref = {'__class__': self.__class__.__name__}
+        arrays = {}
+
+        # Store DataArrays with references
+        arrays[str(self.cluster_order.name)] = self.cluster_order
+        ref['cluster_order'] = f':::{self.cluster_order.name}'
+
+        arrays[str(self.cluster_occurrences.name)] = self.cluster_occurrences
+        ref['cluster_occurrences'] = f':::{self.cluster_occurrences.name}'
+
+        # Store scalar values
+        if isinstance(self.n_clusters, xr.DataArray):
+            n_clusters_name = self.n_clusters.name or 'n_clusters'
+            self.n_clusters = self.n_clusters.rename(n_clusters_name)
+            arrays[n_clusters_name] = self.n_clusters
+            ref['n_clusters'] = f':::{n_clusters_name}'
+        else:
+            ref['n_clusters'] = int(self.n_clusters)
+
+        ref['timesteps_per_cluster'] = self.timesteps_per_cluster
+
+        return ref, arrays
+
     @property
     def n_original_periods(self) -> int:
         """Number of original periods (before clustering)."""
@@ -281,10 +306,9 @@ def __post_init__(self):
             self.timestep_mapping = self.timestep_mapping.rename('timestep_mapping')
 
         # Ensure representative_weights is a DataArray
+        # Can be (cluster, time) for 2D structure or (time,) for flat structure
         if not isinstance(self.representative_weights, xr.DataArray):
-            self.representative_weights = xr.DataArray(
-                self.representative_weights, dims=['time'], name='representative_weights'
-            )
+            self.representative_weights = xr.DataArray(self.representative_weights, name='representative_weights')
         elif self.representative_weights.name is None:
             self.representative_weights = self.representative_weights.rename('representative_weights')
 
@@ -304,6 +328,37 @@ def __repr__(self) -> str:
             f')'
         )
 
+    def _create_reference_structure(self) -> tuple[dict, dict[str, xr.DataArray]]:
+        """Create reference structure for serialization."""
+        ref = {'__class__': self.__class__.__name__}
+        arrays = {}
+
+        # Store DataArrays with references
+        arrays[str(self.timestep_mapping.name)] = self.timestep_mapping
+        ref['timestep_mapping'] = f':::{self.timestep_mapping.name}'
+
+        arrays[str(self.representative_weights.name)] = self.representative_weights
+        ref['representative_weights'] = f':::{self.representative_weights.name}'
+
+        # Store scalar values
+        if isinstance(self.n_representatives, xr.DataArray):
+            n_rep_name = self.n_representatives.name or 'n_representatives'
+            self.n_representatives = self.n_representatives.rename(n_rep_name)
+            arrays[n_rep_name] = self.n_representatives
+            ref['n_representatives'] = f':::{n_rep_name}'
+        else:
+            ref['n_representatives'] = int(self.n_representatives)
+
+        # Store nested ClusterStructure if present
+        if self.cluster_structure is not None:
+            cs_ref, cs_arrays = self.cluster_structure._create_reference_structure()
+            ref['cluster_structure'] = cs_ref
+            arrays.update(cs_arrays)
+
+        # Skip aggregated_data and original_data - not needed for serialization
+
+        return ref, arrays
+
     @property
     def n_original_timesteps(self) -> int:
         """Number of original timesteps (before aggregation)."""
@@ -460,24 +515,50 @@ def validate(self) -> None:
         if max_idx >= n_rep:
             raise ValueError(f'timestep_mapping contains index {max_idx} but n_representatives is {n_rep}')
 
-        # Check weights length matches n_representatives
-        if len(self.representative_weights) != n_rep:
-            raise ValueError(
-                f'representative_weights has {len(self.representative_weights)} elements '
-                f'but n_representatives is {n_rep}'
-            )
+        # Check weights dimensions
+        # representative_weights should have (cluster,) dimension with n_clusters elements
+        # (plus optional period/scenario dimensions)
+        if self.cluster_structure is not None:
+            n_clusters = self.cluster_structure.n_clusters
+            if 'cluster' in self.representative_weights.dims:
+                weights_n_clusters = self.representative_weights.sizes['cluster']
+                if weights_n_clusters != n_clusters:
+                    raise ValueError(
+                        f'representative_weights has {weights_n_clusters} clusters '
+                        f'but cluster_structure has {n_clusters}'
+                    )
 
-        # Check weights sum roughly equals original timesteps
-        weight_sum = float(self.representative_weights.sum().values)
-        n_original = self.n_original_timesteps
-        if abs(weight_sum - n_original) > 1e-6:
-            # Warning only - some aggregation methods may not preserve this exactly
-            import warnings
+        # Check weights sum roughly equals number of original periods
+        # (each weight is how many original periods that cluster represents)
+        # Sum should be checked per period/scenario slice, not across all dimensions
+        if self.cluster_structure is not None:
+            n_original_periods = self.cluster_structure.n_original_periods
+            # Sum over cluster dimension only (keep period/scenario if present)
+            weight_sum_per_slice = self.representative_weights.sum(dim='cluster')
+            # Check each slice
+            if weight_sum_per_slice.size == 1:
+                # Simple case: no period/scenario
+                weight_sum = float(weight_sum_per_slice.values)
+                if abs(weight_sum - n_original_periods) > 1e-6:
+                    import warnings
+
+                    warnings.warn(
+                        f'representative_weights sum ({weight_sum}) does not match '
+                        f'n_original_periods ({n_original_periods})',
+                        stacklevel=2,
+                    )
+            else:
+                # Multi-dimensional: check each slice
+                for val in weight_sum_per_slice.values.flat:
+                    if abs(float(val) - n_original_periods) > 1e-6:
+                        import warnings
 
-            warnings.warn(
-                f'representative_weights sum ({weight_sum}) does not match n_original_timesteps ({n_original})',
-                stacklevel=2,
-            )
+                        warnings.warn(
+                            f'representative_weights sum per slice ({float(val)}) does not match '
+                            f'n_original_periods ({n_original_periods})',
+                            stacklevel=2,
+                        )
+                        break  # Only warn once
 
 
 class ClusteringPlotAccessor:
@@ -911,7 +992,6 @@ class Clustering:
 
     Attributes:
         result: The ClusterResult from the aggregation backend.
-        original_flow_system: Reference to the FlowSystem before aggregation.
         backend_name: Name of the aggregation backend used (e.g., 'tsam', 'manual').
 
     Example:
@@ -923,9 +1003,23 @@ class Clustering:
     """
 
     result: ClusterResult
-    original_flow_system: FlowSystem  # FlowSystem - avoid circular import
     backend_name: str = 'unknown'
 
+    def _create_reference_structure(self) -> tuple[dict, dict[str, xr.DataArray]]:
+        """Create reference structure for serialization."""
+        ref = {'__class__': self.__class__.__name__}
+        arrays = {}
+
+        # Store nested ClusterResult
+        result_ref, result_arrays = self.result._create_reference_structure()
+        ref['result'] = result_ref
+        arrays.update(result_arrays)
+
+        # Store scalar values
+        ref['backend_name'] = self.backend_name
+
+        return ref, arrays
+
     def __repr__(self) -> str:
         cs = self.result.cluster_structure
         if cs is not None:
@@ -1024,6 +1118,22 @@ def cluster_start_positions(self) -> np.ndarray:
         n_timesteps = self.n_clusters * self.timesteps_per_period
         return np.arange(0, n_timesteps, self.timesteps_per_period)
 
+    @property
+    def original_timesteps(self) -> pd.DatetimeIndex:
+        """Original timesteps before clustering.
+
+        Derived from the 'original_time' coordinate of timestep_mapping.
+
+        Raises:
+            KeyError: If 'original_time' coordinate is missing from timestep_mapping.
+        """
+        if 'original_time' not in self.result.timestep_mapping.coords:
+            raise KeyError(
+                "timestep_mapping is missing 'original_time' coordinate. "
+                'This may indicate corrupted or incompatible clustering results.'
+            )
+        return pd.DatetimeIndex(self.result.timestep_mapping.coords['original_time'].values)
+
 
 def create_cluster_structure_from_mapping(
     timestep_mapping: xr.DataArray,
@@ -1073,3 +1183,15 @@ def create_cluster_structure_from_mapping(
         n_clusters=n_clusters,
         timesteps_per_cluster=timesteps_per_cluster,
     )
+
+
+def _register_clustering_classes():
+    """Register clustering classes for IO.
+
+    Called from flow_system.py after all imports are complete to avoid circular imports.
+    """
+    from ..structure import CLASS_REGISTRY
+
+    CLASS_REGISTRY['ClusterStructure'] = ClusterStructure
+    CLASS_REGISTRY['ClusterResult'] = ClusterResult
+    CLASS_REGISTRY['Clustering'] = Clustering
diff --git a/flixopt/components.py b/flixopt/components.py
index 0f2a6077e..390fc6f02 100644
--- a/flixopt/components.py
+++ b/flixopt/components.py
@@ -1257,10 +1257,12 @@ def _absolute_charge_state_bounds(self) -> tuple[xr.DataArray, xr.DataArray]:
             return -np.inf, np.inf
         elif isinstance(self.element.capacity_in_flow_hours, InvestParameters):
             cap_max = self.element.capacity_in_flow_hours.maximum_or_fixed_size * relative_upper_bound
-            return -cap_max, cap_max
+            # Adding 0.0 converts -0.0 to 0.0 (linopy LP writer bug workaround)
+            return -cap_max + 0.0, cap_max + 0.0
         else:
             cap = self.element.capacity_in_flow_hours * relative_upper_bound
-            return -cap, cap
+            # Adding 0.0 converts -0.0 to 0.0 (linopy LP writer bug workaround)
+            return -cap + 0.0, cap + 0.0
 
     def _do_modeling(self):
         """Create storage model with inter-cluster linking constraints.
diff --git a/flixopt/core.py b/flixopt/core.py
index fdcab029b..3d456fff1 100644
--- a/flixopt/core.py
+++ b/flixopt/core.py
@@ -15,7 +15,7 @@
 
 logger = logging.getLogger('flixopt')
 
-FlowSystemDimensions = Literal['time', 'period', 'scenario']
+FlowSystemDimensions = Literal['time', 'cluster', 'period', 'scenario']
 """Possible dimensions of a FlowSystem."""
 
 
diff --git a/flixopt/elements.py b/flixopt/elements.py
index ba2b72f80..0cee53738 100644
--- a/flixopt/elements.py
+++ b/flixopt/elements.py
@@ -677,14 +677,10 @@ def _do_modeling(self):
         self._constraint_flow_rate()
 
         # Total flow hours tracking (per period)
-        # Sum over all temporal dimensions (time, and cluster if present)
-        weighted_flow = self.flow_rate * self._model.aggregation_weight
-        # Get temporal_dims from aggregation_weight (not weighted_flow which has linopy's _term dim)
-        temporal_dims = [d for d in self._model.aggregation_weight.dims if d not in ('period', 'scenario')]
         ModelingPrimitives.expression_tracking_variable(
             model=self,
             name=f'{self.label_full}|total_flow_hours',
-            tracked_expression=weighted_flow.sum(temporal_dims),
+            tracked_expression=self._model.sum_temporal(self.flow_rate),
             bounds=(
                 self.element.flow_hours_min if self.element.flow_hours_min is not None else 0,
                 self.element.flow_hours_max if self.element.flow_hours_max is not None else None,
@@ -841,9 +837,8 @@ def _create_bounds_for_load_factor(self):
         # Get the size (either from element or investment)
         size = self.investment.size if self.with_investment else self.element.size
 
-        # Sum over all temporal dimensions (time, and cluster if present)
-        temporal_dims = [d for d in self._model.aggregation_weight.dims if d not in ('period', 'scenario')]
-        total_hours = self._model.aggregation_weight.sum(temporal_dims)
+        # Total hours in the period (sum of temporal weights)
+        total_hours = self._model.temporal_weight.sum(self._model.temporal_dims)
 
         # Maximum load factor constraint
         if self.element.load_factor_max is not None:
@@ -959,9 +954,7 @@ def _do_modeling(self):
 
         # Add virtual supply/demand to balance and penalty if needed
         if self.element.allows_imbalance:
-            imbalance_penalty = np.multiply(
-                self._model.aggregation_weight, self.element.imbalance_penalty_per_flow_hour
-            )
+            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'
diff --git a/flixopt/features.py b/flixopt/features.py
index e0a018a7f..289640ddd 100644
--- a/flixopt/features.py
+++ b/flixopt/features.py
@@ -196,22 +196,14 @@ def _do_modeling(self):
             inactive = self.add_variables(binary=True, short_name='inactive', coords=self._model.get_coords())
             self.add_constraints(self.status + inactive == 1, short_name='complementary')
 
-        # 3. Total duration tracking using existing pattern
-        # Sum over all temporal dimensions (time, and cluster if present)
-        weighted_status = self.status * self._model.aggregation_weight
-        # Get temporal_dims from aggregation_weight (not weighted_status which has linopy's _term dim)
-        temporal_dims = [d for d in self._model.aggregation_weight.dims if d not in ('period', 'scenario')]
-        agg_weight_sum = self._model.aggregation_weight.sum(temporal_dims)
+        # 3. Total duration tracking
+        total_hours = self._model.temporal_weight.sum(self._model.temporal_dims)
         ModelingPrimitives.expression_tracking_variable(
             self,
-            tracked_expression=weighted_status.sum(temporal_dims),
+            tracked_expression=self._model.sum_temporal(self.status),
             bounds=(
                 self.parameters.active_hours_min if self.parameters.active_hours_min is not None else 0,
-                self.parameters.active_hours_max
-                if self.parameters.active_hours_max is not None
-                else agg_weight_sum.max().item()
-                if hasattr(agg_weight_sum, 'max')
-                else agg_weight_sum,
+                self.parameters.active_hours_max if self.parameters.active_hours_max is not None else total_hours,
             ),
             short_name='active_hours',
             coords=['period', 'scenario'],
@@ -631,10 +623,8 @@ def _do_modeling(self):
 
             # Add it to the total (cluster_weight handles cluster representation, defaults to 1.0)
             # Sum over all temporal dimensions (time, and cluster if present)
-            weighted_per_timestep = self.total_per_timestep * self._model.cluster_weight
-            # Get temporal_dims from total_per_timestep (linopy Variable) - its coords are the actual dims
-            temporal_dims = [d for d in self.total_per_timestep.dims if d not in ('period', 'scenario')]
-            self._eq_total.lhs -= weighted_per_timestep.sum(dim=temporal_dims)
+            weighted_per_timestep = self.total_per_timestep * self._model.weights.get('cluster', 1.0)
+            self._eq_total.lhs -= weighted_per_timestep.sum(dim=self._model.temporal_dims)
 
     def add_share(
         self,
diff --git a/flixopt/flow_system.py b/flixopt/flow_system.py
index c10a1defb..7e12029ed 100644
--- a/flixopt/flow_system.py
+++ b/flixopt/flow_system.py
@@ -40,11 +40,15 @@
 
     from .clustering import Clustering
     from .solvers import _Solver
-    from .structure import TimeSeriesWeights
     from .types import Effect_TPS, Numeric_S, Numeric_TPS, NumericOrBool
 
 from .carrier import Carrier, CarrierContainer
 
+# Register clustering classes for IO (deferred to avoid circular imports)
+from .clustering.base import _register_clustering_classes
+
+_register_clustering_classes()
+
 logger = logging.getLogger('flixopt')
 
 
@@ -69,10 +73,10 @@ class FlowSystem(Interface, CompositeContainerMixin[Element]):
         scenario_weights: The weights of each scenario. If None, all scenarios have the same weight (normalized to 1).
             Period weights are always computed internally from the period index (like timestep_duration for time).
             The final `weights` array (accessible via `flow_system.model.objective_weights`) is computed as period_weights × normalized_scenario_weights, with normalization applied to the scenario weights by default.
-        cluster_weight: Weight for each timestep representing cluster representation count.
-            If None (default), all timesteps have weight 1.0. Used by cluster() to specify
-            how many original timesteps each cluster represents. Combined with timestep_duration
-            via aggregation_weight for proper time aggregation in clustered models.
+        cluster_weight: Weight for each cluster.
+            If None (default), all clusters have weight 1.0. Used by cluster() to specify
+            how many original timesteps each cluster represents. Multiply with timestep_duration
+            for proper time aggregation in clustered models.
         scenario_independent_sizes: Controls whether investment sizes are equalized across scenarios.
             - True: All sizes are shared/equalized across scenarios
             - False: All sizes are optimized separately per scenario
@@ -201,10 +205,13 @@ def __init__(
         # Cluster weight for cluster() optimization (default 1.0)
         # Represents how many original timesteps each cluster represents
         # May have period/scenario dimensions if cluster() was used with those
-        self.cluster_weight = self.fit_to_model_coords(
-            'cluster_weight',
-            np.ones(len(self.timesteps)) if cluster_weight is None else cluster_weight,
-            dims=['time', 'period', 'scenario'],  # Gracefully ignores dims not present
+        self.cluster_weight: xr.DataArray | None = (
+            self.fit_to_model_coords(
+                'cluster_weight',
+                cluster_weight,
+            )
+            if cluster_weight is not None
+            else None
         )
 
         self.scenario_weights = scenario_weights  # Use setter
@@ -502,6 +509,9 @@ def _update_period_metadata(
         period index. This ensures period metadata stays synchronized with the actual
         periods after operations like selection.
 
+        When the period dimension is dropped (single value selected), this method
+        removes the scalar coordinate, period_weights DataArray, and cleans up attributes.
+
         This is analogous to _update_time_metadata() for time-related metadata.
 
         Args:
@@ -513,7 +523,16 @@ def _update_period_metadata(
             The same dataset with updated period-related attributes and data variables
         """
         new_period_index = dataset.indexes.get('period')
-        if new_period_index is not None and len(new_period_index) >= 1:
+
+        if new_period_index is None:
+            # Period dimension was dropped (single value selected)
+            if 'period' in dataset.coords:
+                dataset = dataset.drop_vars('period')
+            dataset = dataset.drop_vars(['period_weights'], errors='ignore')
+            dataset.attrs.pop('weight_of_last_period', None)
+            return dataset
+
+        if len(new_period_index) >= 1:
             # Reuse stored weight_of_last_period when not explicitly overridden.
             # This is essential for single-period subsets where it cannot be inferred from intervals.
             if weight_of_last_period is None:
@@ -542,6 +561,9 @@ def _update_scenario_metadata(cls, dataset: xr.Dataset) -> xr.Dataset:
         Recomputes or removes scenario weights. This ensures scenario metadata stays synchronized with the actual
         scenarios after operations like selection.
 
+        When the scenario dimension is dropped (single value selected), this method
+        removes the scalar coordinate, scenario_weights DataArray, and cleans up attributes.
+
         This is analogous to _update_period_metadata() for time-related metadata.
 
         Args:
@@ -551,7 +573,16 @@ def _update_scenario_metadata(cls, dataset: xr.Dataset) -> xr.Dataset:
             The same dataset with updated scenario-related attributes and data variables
         """
         new_scenario_index = dataset.indexes.get('scenario')
-        if new_scenario_index is None or len(new_scenario_index) <= 1:
+
+        if new_scenario_index is None:
+            # Scenario dimension was dropped (single value selected)
+            if 'scenario' in dataset.coords:
+                dataset = dataset.drop_vars('scenario')
+            dataset = dataset.drop_vars(['scenario_weights'], errors='ignore')
+            dataset.attrs.pop('scenario_weights', None)
+            return dataset
+
+        if len(new_scenario_index) <= 1:
             dataset.attrs.pop('scenario_weights', None)
 
         return dataset
@@ -645,13 +676,21 @@ def to_dataset(self, include_solution: bool = True) -> xr.Dataset:
                 carriers_structure[name] = carrier_ref
             ds.attrs['carriers'] = json.dumps(carriers_structure)
 
-        # Include cluster info for clustered FlowSystems
+        # Include cluster info for clustered FlowSystems (old structure)
         if self.clusters is not None:
             ds.attrs['is_clustered'] = True
             ds.attrs['n_clusters'] = len(self.clusters)
             ds.attrs['timesteps_per_cluster'] = len(self.timesteps)
             ds.attrs['timestep_duration'] = float(self.timestep_duration.mean())
 
+        # Serialize Clustering object if present (new structure)
+        if self.clustering is not None:
+            clustering_ref, clustering_arrays = self.clustering._create_reference_structure()
+            # Add clustering arrays with prefix
+            for name, arr in clustering_arrays.items():
+                ds[f'clustering|{name}'] = arr
+            ds.attrs['clustering'] = json.dumps(clustering_ref)
+
         # Add version info
         ds.attrs['flixopt_version'] = __version__
 
@@ -708,6 +747,11 @@ def from_dataset(cls, ds: xr.Dataset) -> FlowSystem:
                 else None
             )
 
+        # Resolve scenario_weights only if scenario dimension exists
+        scenario_weights = None
+        if ds.indexes.get('scenario') is not None and 'scenario_weights' in reference_structure:
+            scenario_weights = cls._resolve_dataarray_reference(reference_structure['scenario_weights'], arrays_dict)
+
         # Create FlowSystem instance with constructor parameters
         flow_system = cls(
             timesteps=ds.indexes['time'],
@@ -717,9 +761,7 @@ def from_dataset(cls, ds: xr.Dataset) -> FlowSystem:
             hours_of_last_timestep=reference_structure.get('hours_of_last_timestep'),
             hours_of_previous_timesteps=reference_structure.get('hours_of_previous_timesteps'),
             weight_of_last_period=reference_structure.get('weight_of_last_period'),
-            scenario_weights=cls._resolve_dataarray_reference(reference_structure['scenario_weights'], arrays_dict)
-            if 'scenario_weights' in reference_structure
-            else None,
+            scenario_weights=scenario_weights,
             cluster_weight=cluster_weight_for_constructor,
             scenario_independent_sizes=reference_structure.get('scenario_independent_sizes', True),
             scenario_independent_flow_rates=reference_structure.get('scenario_independent_flow_rates', False),
@@ -765,6 +807,19 @@ def from_dataset(cls, ds: xr.Dataset) -> FlowSystem:
                 carrier = cls._resolve_reference_structure(carrier_data, {})
                 flow_system._carriers.add(carrier)
 
+        # Restore Clustering object if present
+        if 'clustering' in reference_structure:
+            clustering_structure = json.loads(reference_structure['clustering'])
+            # Collect clustering arrays (prefixed with 'clustering|')
+            clustering_arrays = {}
+            for name, arr in ds.data_vars.items():
+                if name.startswith('clustering|'):
+                    # Remove 'clustering|' prefix (11 chars)
+                    arr_name = name[11:]
+                    clustering_arrays[arr_name] = arr
+            clustering = cls._resolve_reference_structure(clustering_structure, clustering_arrays)
+            flow_system.clustering = clustering
+
         # Reconnect network to populate bus inputs/outputs (not stored in NetCDF).
         flow_system.connect_and_transform()
 
@@ -1061,6 +1116,7 @@ def connect_and_transform(self):
         self._connect_network()
         self._register_missing_carriers()
         self._assign_element_colors()
+
         for element in chain(self.components.values(), self.effects.values(), self.buses.values()):
             element.transform_data()
 
@@ -1274,22 +1330,29 @@ def flow_carriers(self) -> dict[str, str]:
 
         return self._flow_carriers
 
-    def create_model(self, normalize_weights: bool = True) -> FlowSystemModel:
+    def create_model(self, normalize_weights: bool | None = None) -> FlowSystemModel:
         """
         Create a linopy model from the FlowSystem.
 
         Args:
-            normalize_weights: Whether to automatically normalize the weights (periods and scenarios) to sum up to 1 when solving.
+            normalize_weights: Deprecated. Scenario weights are now always normalized in FlowSystem.
         """
+        if normalize_weights is not None:
+            warnings.warn(
+                f'\n\nnormalize_weights parameter is deprecated and will be removed in {DEPRECATION_REMOVAL_VERSION}. '
+                'Scenario weights are now always normalized when set on FlowSystem.\n',
+                DeprecationWarning,
+                stacklevel=2,
+            )
         if not self.connected_and_transformed:
             raise RuntimeError(
                 'FlowSystem is not connected_and_transformed. Call FlowSystem.connect_and_transform() first.'
             )
         # System integrity was already validated in connect_and_transform()
-        self.model = FlowSystemModel(self, normalize_weights)
+        self.model = FlowSystemModel(self)
         return self.model
 
-    def build_model(self, normalize_weights: bool = True) -> FlowSystem:
+    def build_model(self, normalize_weights: bool | None = None) -> FlowSystem:
         """
         Build the optimization model for this FlowSystem.
 
@@ -1303,7 +1366,7 @@ def build_model(self, normalize_weights: bool = True) -> FlowSystem:
         before solving.
 
         Args:
-            normalize_weights: Whether to normalize scenario/period weights to sum to 1.
+            normalize_weights: Deprecated. Scenario weights are now always normalized in FlowSystem.
 
         Returns:
             Self, for method chaining.
@@ -1313,8 +1376,15 @@ def build_model(self, normalize_weights: bool = True) -> FlowSystem:
             >>> print(flow_system.model.variables)  # Inspect variables before solving
             >>> flow_system.solve(solver)
         """
+        if normalize_weights is not None:
+            warnings.warn(
+                f'\n\nnormalize_weights parameter is deprecated and will be removed in {DEPRECATION_REMOVAL_VERSION}. '
+                'Scenario weights are now always normalized when set on FlowSystem.\n',
+                DeprecationWarning,
+                stacklevel=2,
+            )
         self.connect_and_transform()
-        self.create_model(normalize_weights)
+        self.create_model()
 
         self.model.do_modeling()
 
@@ -1862,27 +1932,85 @@ def storages(self) -> ElementContainer[Storage]:
             self._storages_cache = ElementContainer(storages, element_type_name='storages', truncate_repr=10)
         return self._storages_cache
 
+    @property
+    def dims(self) -> list[str]:
+        """Active dimension names.
+
+        Returns:
+            List of active dimension names in order.
+
+        Example:
+            >>> fs.dims
+            ['time']  # simple case
+            >>> fs_clustered.dims
+            ['cluster', 'time', 'period', 'scenario']  # full case
+        """
+        result = []
+        if self.clusters is not None:
+            result.append('cluster')
+        result.append('time')
+        if self.periods is not None:
+            result.append('period')
+        if self.scenarios is not None:
+            result.append('scenario')
+        return result
+
+    @property
+    def indexes(self) -> dict[str, pd.Index]:
+        """Indexes for active dimensions.
+
+        Returns:
+            Dict mapping dimension names to pandas Index objects.
+
+        Example:
+            >>> fs.indexes['time']
+            DatetimeIndex(['2024-01-01', ...], dtype='datetime64[ns]', name='time')
+        """
+        result: dict[str, pd.Index] = {}
+        if self.clusters is not None:
+            result['cluster'] = self.clusters
+        result['time'] = self.timesteps
+        if self.periods is not None:
+            result['period'] = self.periods
+        if self.scenarios is not None:
+            result['scenario'] = self.scenarios
+        return result
+
+    @property
+    def temporal_dims(self) -> list[str]:
+        """Temporal dimensions for summing over time.
+
+        Returns ['time', 'cluster'] for clustered systems, ['time'] otherwise.
+        """
+        if self.clusters is not None:
+            return ['time', 'cluster']
+        return ['time']
+
+    @property
+    def temporal_weight(self) -> xr.DataArray:
+        """Combined temporal weight (timestep_duration × cluster_weight).
+
+        Use for converting rates to totals before summing.
+        Note: cluster_weight is used even without a clusters dimension.
+        """
+        # Use cluster_weight directly if set, otherwise check weights dict, fallback to 1.0
+        cluster_weight = self.weights.get('cluster', self.cluster_weight if self.cluster_weight is not None else 1.0)
+        return self.weights['time'] * cluster_weight
+
     @property
     def coords(self) -> dict[FlowSystemDimensions, pd.Index]:
         """Active coordinates for variable creation.
 
+        .. deprecated::
+            Use :attr:`indexes` instead.
+
         Returns a dict of dimension names to coordinate arrays. When clustered,
         includes 'cluster' dimension before 'time'.
 
         Returns:
             Dict mapping dimension names to coordinate arrays.
         """
-        active_coords: dict[str, pd.Index] = {}
-
-        if self.clusters is not None:
-            active_coords['cluster'] = self.clusters
-        active_coords['time'] = self.timesteps
-
-        if self.periods is not None:
-            active_coords['period'] = self.periods
-        if self.scenarios is not None:
-            active_coords['scenario'] = self.scenarios
-        return active_coords
+        return self.indexes
 
     @property
     def _use_true_cluster_dims(self) -> bool:
@@ -1928,14 +2056,15 @@ def scenario_weights(self) -> xr.DataArray | None:
     @scenario_weights.setter
     def scenario_weights(self, value: Numeric_S | None) -> None:
         """
-        Set scenario weights.
+        Set scenario weights (always normalized to sum to 1).
 
         Args:
-            value: Scenario weights to set (will be converted to DataArray with 'scenario' dimension)
-                or None to clear weights.
+            value: Scenario weights to set (will be converted to DataArray with 'scenario' dimension
+                and normalized to sum to 1), or None to clear weights.
 
         Raises:
             ValueError: If value is not None and no scenarios are defined in the FlowSystem.
+            ValueError: If weights sum to zero (cannot normalize).
         """
         if value is None:
             self._scenario_weights = None
@@ -1947,48 +2076,65 @@ def scenario_weights(self, value: Numeric_S | None) -> None:
                 'Either define scenarios in FlowSystem(scenarios=...) or set scenario_weights to None.'
             )
 
-        self._scenario_weights = self.fit_to_model_coords('scenario_weights', value, dims=['scenario'])
+        weights = self.fit_to_model_coords('scenario_weights', value, dims=['scenario'])
 
-    @property
-    def weights(self) -> TimeSeriesWeights:
-        """Unified weighting system for time series aggregation.
+        # Normalize to sum to 1
+        norm = weights.sum('scenario')
+        if np.isclose(norm, 0.0).any():
+            raise ValueError('scenario_weights sum to 0; cannot normalize.')
+        self._scenario_weights = weights / norm
 
-        Returns a TimeSeriesWeights object providing a clean, unified interface
-        for all weight types used in flixopt. This is the recommended way to
-        access weights for new code (PyPSA-inspired design).
+    def _unit_weight(self, dim: str) -> xr.DataArray:
+        """Create a unit weight DataArray (all 1.0) for a dimension."""
+        index = self.indexes[dim]
+        return xr.DataArray(
+            np.ones(len(index), dtype=float),
+            coords={dim: index},
+            dims=[dim],
+            name=f'{dim}_weight',
+        )
 
-        The temporal weight combines timestep_duration and cluster_weight,
-        which is the proper weight for summing over time.
+    @property
+    def weights(self) -> dict[str, xr.DataArray]:
+        """Weights for active dimensions (unit weights if not explicitly set).
 
         Returns:
-            TimeSeriesWeights with temporal, period, and scenario weights.
+            Dict mapping dimension names to weight DataArrays.
+            Keys match :attr:`dims` and :attr:`indexes`.
 
         Example:
-            >>> weights = flow_system.weights
-            >>> weighted_total = (flow_rate * weights.temporal).sum('time')
-            >>> # Or use the convenience method:
-            >>> weighted_total = weights.sum_over_time(flow_rate)
+            >>> fs.weights['time']  # timestep durations
+            >>> fs.weights['cluster']  # cluster weights (unit if not set)
         """
-        from .structure import TimeSeriesWeights
+        result: dict[str, xr.DataArray] = {'time': self.timestep_duration}
+        if self.clusters is not None:
+            result['cluster'] = self.cluster_weight if self.cluster_weight is not None else self._unit_weight('cluster')
+        if self.periods is not None:
+            result['period'] = self.period_weights if self.period_weights is not None else self._unit_weight('period')
+        if self.scenarios is not None:
+            result['scenario'] = (
+                self.scenario_weights if self.scenario_weights is not None else self._unit_weight('scenario')
+            )
+        return result
 
-        return TimeSeriesWeights(
-            temporal=self.timestep_duration * self.cluster_weight,
-            period=self.period_weights,
-            scenario=self._scenario_weights,
-        )
+    def sum_temporal(self, data: xr.DataArray) -> xr.DataArray:
+        """Sum data over temporal dimensions with full temporal weighting.
 
-    @property
-    def aggregation_weight(self) -> xr.DataArray:
-        """Combined weight for time aggregation.
+        Applies both timestep_duration and cluster_weight, then sums over temporal dimensions.
+        Use this to convert rates to totals (e.g., flow_rate → total_energy).
+
+        Args:
+            data: Data with time dimension (and optionally cluster).
+                  Typically a rate (e.g., flow_rate in MW, status as 0/1).
 
-        Combines timestep_duration (physical duration) and cluster_weight (cluster representation).
-        Use this for proper time aggregation in clustered models.
+        Returns:
+            Data summed over temporal dims with full temporal weighting applied.
 
-        Note:
-            This is equivalent to `weights.temporal`. The unified TimeSeriesWeights
-            interface (via `flow_system.weights`) is recommended for new code.
+        Example:
+            >>> total_energy = fs.sum_temporal(flow_rate)  # MW → MWh total
+            >>> active_hours = fs.sum_temporal(status)  # count → hours
         """
-        return self.timestep_duration * self.cluster_weight
+        return (data * self.temporal_weight).sum(self.temporal_dims)
 
     @property
     def is_clustered(self) -> bool:
diff --git a/flixopt/optimization.py b/flixopt/optimization.py
index 6a1a87ce1..0b567387f 100644
--- a/flixopt/optimization.py
+++ b/flixopt/optimization.py
@@ -82,7 +82,7 @@ def _initialize_optimization_common(
     name: str,
     flow_system: FlowSystem,
     folder: pathlib.Path | None = None,
-    normalize_weights: bool = True,
+    normalize_weights: bool | None = None,
 ) -> None:
     """
     Shared initialization logic for all optimization types.
@@ -95,7 +95,7 @@ def _initialize_optimization_common(
         name: Name of the optimization
         flow_system: FlowSystem to optimize
         folder: Directory for saving results
-        normalize_weights: Whether to normalize scenario weights
+        normalize_weights: Deprecated. Scenario weights are now always normalized in FlowSystem.
     """
     obj.name = name
 
@@ -106,7 +106,8 @@ def _initialize_optimization_common(
         )
         flow_system = flow_system.copy()
 
-    obj.normalize_weights = normalize_weights
+    # normalize_weights is deprecated but kept for backwards compatibility
+    obj.normalize_weights = True  # Always True now
 
     flow_system._used_in_optimization = True
 
@@ -186,7 +187,7 @@ def do_modeling(self) -> Optimization:
         t_start = timeit.default_timer()
         self.flow_system.connect_and_transform()
 
-        self.model = self.flow_system.create_model(self.normalize_weights)
+        self.model = self.flow_system.create_model()
         self.model.do_modeling()
 
         self.durations['modeling'] = round(timeit.default_timer() - t_start, 2)
diff --git a/flixopt/optimize_accessor.py b/flixopt/optimize_accessor.py
index f88cdf982..7aee930a4 100644
--- a/flixopt/optimize_accessor.py
+++ b/flixopt/optimize_accessor.py
@@ -53,7 +53,7 @@ def __init__(self, flow_system: FlowSystem) -> None:
         """
         self._fs = flow_system
 
-    def __call__(self, solver: _Solver, normalize_weights: bool = True) -> FlowSystem:
+    def __call__(self, solver: _Solver, normalize_weights: bool | None = None) -> FlowSystem:
         """
         Build and solve the optimization model in one step.
 
@@ -64,7 +64,7 @@ def __call__(self, solver: _Solver, normalize_weights: bool = True) -> FlowSyste
 
         Args:
             solver: The solver to use (e.g., HighsSolver, GurobiSolver).
-            normalize_weights: Whether to normalize scenario/period weights to sum to 1.
+            normalize_weights: Deprecated. Scenario weights are now always normalized in FlowSystem.
 
         Returns:
             The FlowSystem, for method chaining.
@@ -85,7 +85,18 @@ def __call__(self, solver: _Solver, normalize_weights: bool = True) -> FlowSyste
 
             >>> solution = flow_system.optimize(solver).solution
         """
-        self._fs.build_model(normalize_weights)
+        if normalize_weights is not None:
+            import warnings
+
+            from .config import DEPRECATION_REMOVAL_VERSION
+
+            warnings.warn(
+                f'\n\nnormalize_weights parameter is deprecated and will be removed in {DEPRECATION_REMOVAL_VERSION}. '
+                'Scenario weights are now always normalized when set on FlowSystem.\n',
+                DeprecationWarning,
+                stacklevel=2,
+            )
+        self._fs.build_model()
         self._fs.solve(solver)
         return self._fs
 
diff --git a/flixopt/statistics_accessor.py b/flixopt/statistics_accessor.py
index bee26a0e2..73d115df0 100644
--- a/flixopt/statistics_accessor.py
+++ b/flixopt/statistics_accessor.py
@@ -847,7 +847,7 @@ def get_contributor_type(contributor: str) -> str:
                             # For total mode, sum temporal over time (apply cluster_weight for proper weighting)
                             # Sum over all temporal dimensions (time, and cluster if present)
                             if mode == 'total' and current_mode == 'temporal' and 'time' in da.dims:
-                                weighted = da * self._fs.cluster_weight
+                                weighted = da * self._fs.weights.get('cluster', 1.0)
                                 temporal_dims = [d for d in weighted.dims if d not in ('period', 'scenario')]
                                 da = weighted.sum(temporal_dims)
                             if share_total is None:
diff --git a/flixopt/structure.py b/flixopt/structure.py
index 7996565e8..4b7734199 100644
--- a/flixopt/structure.py
+++ b/flixopt/structure.py
@@ -43,92 +43,6 @@
 CLASS_REGISTRY = {}
 
 
-@dataclass
-class TimeSeriesWeights:
-    """Unified weighting system for time series aggregation (PyPSA-inspired).
-
-    This class provides a clean, unified interface for time series weights,
-    combining the various weight types used in flixopt into a single object.
-
-    Attributes:
-        temporal: Combined weight for temporal operations (timestep_duration × cluster_weight).
-            Applied to all time-summing operations. dims: [time] or [time, period, scenario]
-        period: Weight for each period in multi-period optimization.
-            dims: [period] or None
-        scenario: Weight for each scenario in stochastic optimization.
-            dims: [scenario] or None
-        objective: Optional override weight for objective function calculations.
-            If None, uses temporal weight. dims: [time] or [time, period, scenario]
-        storage: Optional override weight for storage balance equations.
-            If None, uses temporal weight. dims: [time] or [time, period, scenario]
-
-    Example:
-        >>> # Access via FlowSystem
-        >>> weights = flow_system.weights
-        >>> weighted_sum = (flow_rate * weights.temporal).sum('time')
-        >>>
-        >>> # With period/scenario weighting
-        >>> total = weighted_sum * weights.period * weights.scenario
-
-    Note:
-        For backwards compatibility, the existing properties (cluster_weight,
-        timestep_duration, aggregation_weight) are still available on FlowSystem
-        and FlowSystemModel.
-    """
-
-    temporal: xr.DataArray
-    period: xr.DataArray | None = None
-    scenario: xr.DataArray | None = None
-    objective: xr.DataArray | None = None
-    storage: xr.DataArray | None = None
-
-    def __post_init__(self):
-        """Validate weights."""
-        if not isinstance(self.temporal, xr.DataArray):
-            raise TypeError('temporal must be an xarray DataArray')
-        if 'time' not in self.temporal.dims:
-            raise ValueError("temporal must have 'time' dimension")
-
-    @property
-    def effective_objective(self) -> xr.DataArray:
-        """Get effective objective weight (override or temporal)."""
-        return self.objective if self.objective is not None else self.temporal
-
-    @property
-    def effective_storage(self) -> xr.DataArray:
-        """Get effective storage weight (override or temporal)."""
-        return self.storage if self.storage is not None else self.temporal
-
-    def sum_over_time(self, data: xr.DataArray) -> xr.DataArray:
-        """Sum data over time dimension with proper weighting.
-
-        Args:
-            data: DataArray with 'time' dimension.
-
-        Returns:
-            Data summed over time with temporal weighting applied.
-        """
-        if 'time' not in data.dims:
-            return data
-        return (data * self.temporal).sum('time')
-
-    def apply_period_scenario_weights(self, data: xr.DataArray) -> xr.DataArray:
-        """Apply period and scenario weights to data.
-
-        Args:
-            data: DataArray, optionally with 'period' and/or 'scenario' dims.
-
-        Returns:
-            Data with period and scenario weights applied.
-        """
-        result = data
-        if self.period is not None and 'period' in data.dims:
-            result = result * self.period
-        if self.scenario is not None and 'scenario' in data.dims:
-            result = result * self.scenario
-        return result
-
-
 def register_class_for_io(cls):
     """Register a class for serialization/deserialization."""
     name = cls.__name__
@@ -176,13 +90,11 @@ class FlowSystemModel(linopy.Model, SubmodelsMixin):
 
     Args:
         flow_system: The flow_system that is used to create the model.
-        normalize_weights: Whether to automatically normalize the weights to sum up to 1 when solving.
     """
 
-    def __init__(self, flow_system: FlowSystem, normalize_weights: bool):
+    def __init__(self, flow_system: FlowSystem):
         super().__init__(force_dim_names=True)
         self.flow_system = flow_system
-        self.normalize_weights = normalize_weights
         self.effects: EffectCollectionModel | None = None
         self.submodels: Submodels = Submodels({})
 
@@ -314,53 +226,63 @@ def hours_of_previous_timesteps(self):
         return self.flow_system.hours_of_previous_timesteps
 
     @property
-    def cluster_weight(self) -> xr.DataArray:
-        """Cluster weight for cluster() optimization.
+    def dims(self) -> list[str]:
+        """Active dimension names."""
+        return self.flow_system.dims
 
-        Represents how many original timesteps each cluster represents.
-        Default is 1.0 for all timesteps.
+    @property
+    def indexes(self) -> dict[str, pd.Index]:
+        """Indexes for active dimensions."""
+        return self.flow_system.indexes
+
+    @property
+    def weights(self) -> dict[str, xr.DataArray]:
+        """Weights for active dimensions (unit weights if not set).
+
+        Scenario weights are always normalized (handled by FlowSystem).
         """
-        return self.flow_system.cluster_weight
+        return self.flow_system.weights
 
     @property
-    def aggregation_weight(self) -> xr.DataArray:
-        """Combined weight for time aggregation.
+    def temporal_dims(self) -> list[str]:
+        """Temporal dimensions for summing over time.
 
-        Combines timestep_duration (physical duration) and cluster_weight (cluster representation).
-        Use this for proper time aggregation in clustered models.
+        Returns ['time', 'cluster'] for clustered systems, ['time'] otherwise.
         """
-        return self.timestep_duration * self.cluster_weight
+        return self.flow_system.temporal_dims
+
+    @property
+    def temporal_weight(self) -> xr.DataArray:
+        """Combined temporal weight (timestep_duration × cluster_weight)."""
+        return self.flow_system.temporal_weight
+
+    def sum_temporal(self, data: xr.DataArray) -> xr.DataArray:
+        """Sum data over temporal dimensions with full temporal weighting.
+
+        Example:
+            >>> total_energy = model.sum_temporal(flow_rate)
+        """
+        return self.flow_system.sum_temporal(data)
 
     @property
     def scenario_weights(self) -> xr.DataArray:
         """
-        Scenario weights of model. With optional normalization.
+        Scenario weights of model (always normalized, via FlowSystem).
+
+        Returns unit weights if no scenarios defined or no explicit weights set.
         """
         if self.flow_system.scenarios is None:
             return xr.DataArray(1)
 
         if self.flow_system.scenario_weights is None:
-            scenario_weights = xr.DataArray(
-                np.ones(self.flow_system.scenarios.size, dtype=float),
-                coords={'scenario': self.flow_system.scenarios},
-                dims=['scenario'],
-                name='scenario_weights',
-            )
-        else:
-            scenario_weights = self.flow_system.scenario_weights
-
-        if not self.normalize_weights:
-            return scenario_weights
+            return self.flow_system._unit_weight('scenario')
 
-        norm = scenario_weights.sum('scenario')
-        if np.isclose(norm, 0.0).any():
-            raise ValueError('FlowSystemModel.scenario_weights: weights sum to 0; cannot normalize.')
-        return scenario_weights / norm
+        return self.flow_system.scenario_weights
 
     @property
     def objective_weights(self) -> xr.DataArray:
         """
-        Objective weights of model. With optional normalization of scenario weights.
+        Objective weights of model (period_weights × scenario_weights).
         """
         period_weights = self.flow_system.effects.objective_effect.submodel.period_weights
         scenario_weights = self.scenario_weights
diff --git a/flixopt/transform_accessor.py b/flixopt/transform_accessor.py
index 93a1e2247..dc63dcacb 100644
--- a/flixopt/transform_accessor.py
+++ b/flixopt/transform_accessor.py
@@ -801,7 +801,11 @@ def _build_cluster_weight_for_key(key: tuple) -> xr.DataArray:
                     da = TimeSeriesData.from_dataarray(da.assign_attrs(original_da.attrs))
                 ds_new_vars[name] = da
 
-        ds_new = xr.Dataset(ds_new_vars, attrs=ds.attrs)
+        # Copy attrs but remove cluster_weight - the clustered FlowSystem gets its own
+        # cluster_weight set after from_dataset (original reference has wrong shape)
+        new_attrs = dict(ds.attrs)
+        new_attrs.pop('cluster_weight', None)
+        ds_new = xr.Dataset(ds_new_vars, attrs=new_attrs)
         ds_new.attrs['timesteps_per_cluster'] = timesteps_per_cluster
         ds_new.attrs['timestep_duration'] = dt
         ds_new.attrs['n_clusters'] = actual_n_clusters
@@ -848,6 +852,9 @@ def _build_cluster_occurrences_for_key(key: tuple) -> np.ndarray:
             timestep_mapping_slices = {}
             cluster_occurrences_slices = {}
 
+            # Use renamed timesteps as coordinates for multi-dimensional case
+            original_timesteps_coord = self._fs.timesteps.rename('original_time')
+
             for p in periods:
                 for s in scenarios:
                     key = (p, s)
@@ -855,7 +862,10 @@ def _build_cluster_occurrences_for_key(key: tuple) -> np.ndarray:
                         cluster_orders[key], dims=['original_period'], name='cluster_order'
                     )
                     timestep_mapping_slices[key] = xr.DataArray(
-                        _build_timestep_mapping_for_key(key), dims=['original_time'], name='timestep_mapping'
+                        _build_timestep_mapping_for_key(key),
+                        dims=['original_time'],
+                        coords={'original_time': original_timesteps_coord},
+                        name='timestep_mapping',
                     )
                     cluster_occurrences_slices[key] = xr.DataArray(
                         _build_cluster_occurrences_for_key(key), dims=['cluster'], name='cluster_occurrences'
@@ -874,8 +884,13 @@ def _build_cluster_occurrences_for_key(key: tuple) -> np.ndarray:
         else:
             # Simple case: single (None, None) slice
             cluster_order_da = xr.DataArray(cluster_orders[first_key], dims=['original_period'], name='cluster_order')
+            # Use renamed timesteps as coordinates
+            original_timesteps_coord = self._fs.timesteps.rename('original_time')
             timestep_mapping_da = xr.DataArray(
-                _build_timestep_mapping_for_key(first_key), dims=['original_time'], name='timestep_mapping'
+                _build_timestep_mapping_for_key(first_key),
+                dims=['original_time'],
+                coords={'original_time': original_timesteps_coord},
+                name='timestep_mapping',
             )
             cluster_occurrences_da = xr.DataArray(
                 _build_cluster_occurrences_for_key(first_key), dims=['cluster'], name='cluster_occurrences'
@@ -888,16 +903,17 @@ def _build_cluster_occurrences_for_key(key: tuple) -> np.ndarray:
             timesteps_per_cluster=timesteps_per_cluster,
         )
 
-        # Create representative_weights in flat format for ClusterResult compatibility
-        # This repeats each cluster's weight for all timesteps within that cluster
-        def _build_flat_weights_for_key(key: tuple) -> xr.DataArray:
+        # Create representative_weights with (cluster,) dimension only
+        # Each cluster has one weight (same for all timesteps within it)
+        def _build_cluster_weights_for_key(key: tuple) -> xr.DataArray:
             occurrences = cluster_occurrences_all[key]
-            weights = np.repeat([occurrences.get(c, 1) for c in range(actual_n_clusters)], timesteps_per_cluster)
-            return xr.DataArray(weights, dims=['time'], name='representative_weights')
+            # Shape: (n_clusters,) - one weight per cluster
+            weights = np.array([occurrences.get(c, 1) for c in range(actual_n_clusters)])
+            return xr.DataArray(weights, dims=['cluster'], name='representative_weights')
 
-        flat_weights_slices = {key: _build_flat_weights_for_key(key) for key in cluster_occurrences_all}
+        weights_slices = {key: _build_cluster_weights_for_key(key) for key in cluster_occurrences_all}
         representative_weights = self._combine_slices_to_dataarray_generic(
-            flat_weights_slices, ['time'], periods, scenarios, 'representative_weights'
+            weights_slices, ['cluster'], periods, scenarios, 'representative_weights'
         )
 
         aggregation_result = ClusterResult(
@@ -911,7 +927,6 @@ def _build_flat_weights_for_key(key: tuple) -> xr.DataArray:
 
         reduced_fs.clustering = Clustering(
             result=aggregation_result,
-            original_flow_system=self._fs,
             backend_name='tsam',
         )
 
@@ -1127,19 +1142,20 @@ def expand_solution(self) -> FlowSystem:
             raise ValueError('No cluster structure available for expansion.')
 
         timesteps_per_cluster = cluster_structure.timesteps_per_cluster
-        original_fs: FlowSystem = info.original_flow_system
         n_clusters = (
             int(cluster_structure.n_clusters)
             if isinstance(cluster_structure.n_clusters, (int, np.integer))
             else int(cluster_structure.n_clusters.values)
         )
-        has_periods = original_fs.periods is not None
-        has_scenarios = original_fs.scenarios is not None
 
-        periods = list(original_fs.periods) if has_periods else [None]
-        scenarios = list(original_fs.scenarios) if has_scenarios else [None]
+        # Get original timesteps from clustering, but periods/scenarios from the FlowSystem
+        # (the clustered FlowSystem preserves the same periods/scenarios)
+        original_timesteps = info.original_timesteps
+        has_periods = self._fs.periods is not None
+        has_scenarios = self._fs.scenarios is not None
 
-        original_timesteps = original_fs.timesteps
+        periods = list(self._fs.periods) if has_periods else [None]
+        scenarios = list(self._fs.scenarios) if has_scenarios else [None]
         n_original_timesteps = len(original_timesteps)
         n_reduced_timesteps = n_clusters * timesteps_per_cluster
 
@@ -1151,11 +1167,23 @@ def expand_da(da: xr.DataArray) -> xr.DataArray:
 
         # 1. Expand FlowSystem data (with cluster_weight set to 1.0 for all timesteps)
         reduced_ds = self._fs.to_dataset(include_solution=False)
-        expanded_ds = xr.Dataset(
-            {name: expand_da(da) for name, da in reduced_ds.data_vars.items() if name != 'cluster_weight'},
-            attrs=reduced_ds.attrs,
-        )
-        expanded_ds.attrs['timestep_duration'] = original_fs.timestep_duration.values.tolist()
+        # Filter out cluster-related variables and copy attrs without clustering info
+        data_vars = {
+            name: expand_da(da)
+            for name, da in reduced_ds.data_vars.items()
+            if name != 'cluster_weight' and not name.startswith('clustering|')
+        }
+        attrs = {
+            k: v
+            for k, v in reduced_ds.attrs.items()
+            if k not in ('is_clustered', 'n_clusters', 'timesteps_per_cluster', 'clustering')
+        }
+        expanded_ds = xr.Dataset(data_vars, attrs=attrs)
+        # Compute timestep_duration from original timesteps
+        # Add extra timestep for duration calculation (assume same interval as last)
+        original_timesteps_extra = FlowSystem._create_timesteps_with_extra(original_timesteps, None)
+        timestep_duration = FlowSystem.calculate_timestep_duration(original_timesteps_extra)
+        expanded_ds.attrs['timestep_duration'] = timestep_duration.values.tolist()
 
         # Create cluster_weight with value 1.0 for all timesteps (no weighting needed for expanded)
         # Use _combine_slices_to_dataarray for consistent multi-dim handling
@@ -1201,8 +1229,8 @@ def expand_da(da: xr.DataArray) -> xr.DataArray:
             soc_boundary_per_timestep = soc_boundary_per_timestep.assign_coords(time=original_timesteps)
 
             # Apply self-discharge decay to SOC_boundary based on time within period
-            # Get the storage's relative_loss_per_hour from original flow system
-            storage = original_fs.storages[storage_name]
+            # Get the storage's relative_loss_per_hour from the clustered flow system
+            storage = self._fs.storages.get(storage_name)
             if storage is not None:
                 # Time within period for each timestep (0, 1, 2, ..., timesteps_per_cluster-1, 0, 1, ...)
                 time_within_period = np.arange(n_original_timesteps) % timesteps_per_cluster
diff --git a/tests/deprecated/test_scenarios.py b/tests/deprecated/test_scenarios.py
index 65ea62d81..2699647ad 100644
--- a/tests/deprecated/test_scenarios.py
+++ b/tests/deprecated/test_scenarios.py
@@ -341,12 +341,14 @@ def test_scenarios_selection(flow_system_piecewise_conversion_scenarios):
 
     assert flow_system.scenarios.equals(flow_system_full.scenarios[0:2])
 
-    np.testing.assert_allclose(flow_system.scenario_weights.values, flow_system_full.scenario_weights[0:2])
+    # Scenario weights are always normalized - subset is re-normalized to sum to 1
+    subset_weights = flow_system_full.scenario_weights[0:2]
+    expected_normalized = subset_weights / subset_weights.sum()
+    np.testing.assert_allclose(flow_system.scenario_weights.values, expected_normalized.values)
 
-    # Optimize using new API with normalize_weights=False
+    # Optimize using new API
     flow_system.optimize(
         fx.solvers.GurobiSolver(mip_gap=0.01, time_limit_seconds=60),
-        normalize_weights=False,
     )
 
     # Penalty has same structure as other effects: 'Penalty' is the total, 'Penalty(temporal)' and 'Penalty(periodic)' are components
@@ -769,7 +771,10 @@ def test_weights_selection():
 
     # Verify weights are correctly sliced
     assert fs_subset.scenarios.equals(pd.Index(['base', 'high'], name='scenario'))
-    np.testing.assert_allclose(fs_subset.scenario_weights.values, custom_scenario_weights[[0, 2]])
+    # Scenario weights are always normalized - subset is re-normalized to sum to 1
+    subset_weights = np.array([0.3, 0.2])  # Original weights for selected scenarios
+    expected_normalized = subset_weights / subset_weights.sum()
+    np.testing.assert_allclose(fs_subset.scenario_weights.values, expected_normalized)
 
     # Verify weights are 1D with just scenario dimension (no period dimension)
     assert fs_subset.scenario_weights.dims == ('scenario',)
diff --git a/tests/test_cluster_reduce_expand.py b/tests/test_cluster_reduce_expand.py
index af2864563..7072fe22e 100644
--- a/tests/test_cluster_reduce_expand.py
+++ b/tests/test_cluster_reduce_expand.py
@@ -292,8 +292,9 @@ def test_cluster_with_scenarios(timesteps_8_days, scenarios_2):
     assert info is not None
     assert info.result.cluster_structure is not None
     assert info.result.cluster_structure.n_clusters == 2
-    # Original FlowSystem had scenarios
-    assert info.original_flow_system.scenarios is not None
+    # Clustered FlowSystem preserves scenarios
+    assert fs_reduced.scenarios is not None
+    assert len(fs_reduced.scenarios) == 2
 
 
 def test_cluster_and_expand_with_scenarios(solver_fixture, timesteps_8_days, scenarios_2):
@@ -465,8 +466,8 @@ def test_storage_cluster_mode_intercluster_cyclic(self, solver_fixture, timestep
         assert 'cluster_boundary' in soc_boundary.dims
 
         # First and last SOC_boundary values should be equal (cyclic constraint)
-        first_soc = float(soc_boundary.isel(cluster_boundary=0).values)
-        last_soc = float(soc_boundary.isel(cluster_boundary=-1).values)
+        first_soc = soc_boundary.isel(cluster_boundary=0).item()
+        last_soc = soc_boundary.isel(cluster_boundary=-1).item()
         assert_allclose(first_soc, last_soc, rtol=1e-6)
 
 
@@ -543,15 +544,15 @@ def test_expanded_charge_state_matches_manual_calculation(self, solver_fixture,
         # Manual verification for first few timesteps of first period
         p = 0  # First period
         cluster = int(cluster_order[p])
-        soc_b = float(soc_boundary.isel(cluster_boundary=p).values)
+        soc_b = soc_boundary.isel(cluster_boundary=p).item()
 
         for t in [0, 5, 12, 23]:
             global_t = p * timesteps_per_cluster + t
-            delta_e = float(cs_clustered.isel(cluster=cluster, time=t).values)
+            delta_e = cs_clustered.isel(cluster=cluster, time=t).item()
             decay = (1 - loss_rate) ** t
             expected = soc_b * decay + delta_e
             expected_clipped = max(0.0, expected)
-            actual = float(cs_expanded.isel(time=global_t).values)
+            actual = cs_expanded.isel(time=global_t).item()
 
             assert_allclose(
                 actual,
diff --git a/tests/test_clustering/test_base.py b/tests/test_clustering/test_base.py
index a6c4d8cc7..9c63f25f6 100644
--- a/tests/test_clustering/test_base.py
+++ b/tests/test_clustering/test_base.py
@@ -152,7 +152,6 @@ def test_creation(self):
 
         info = Clustering(
             result=result,
-            original_flow_system=None,  # Would be FlowSystem in practice
             backend_name='tsam',
         )
 
diff --git a/tests/test_clustering/test_integration.py b/tests/test_clustering/test_integration.py
index 587e39160..2d04a51c1 100644
--- a/tests/test_clustering/test_integration.py
+++ b/tests/test_clustering/test_integration.py
@@ -5,85 +5,121 @@
 import pytest
 import xarray as xr
 
-from flixopt import FlowSystem, TimeSeriesWeights
+from flixopt import FlowSystem
 
 
-class TestTimeSeriesWeights:
-    """Tests for TimeSeriesWeights class."""
+class TestWeights:
+    """Tests for FlowSystem.weights dict property."""
 
-    def test_creation(self):
-        """Test TimeSeriesWeights creation."""
-        temporal = xr.DataArray([1.0, 1.0, 1.0], dims=['time'])
-        weights = TimeSeriesWeights(temporal=temporal)
+    def test_weights_is_dict(self):
+        """Test weights returns a dict."""
+        fs = FlowSystem(timesteps=pd.date_range('2024-01-01', periods=24, freq='h'))
+        weights = fs.weights
+
+        assert isinstance(weights, dict)
+        assert 'time' in weights
+
+    def test_time_weight(self):
+        """Test weights['time'] returns timestep_duration."""
+        fs = FlowSystem(timesteps=pd.date_range('2024-01-01', periods=24, freq='h'))
+        weights = fs.weights
+
+        # For hourly data, timestep_duration is 1.0
+        assert float(weights['time'].mean()) == 1.0
+
+    def test_cluster_not_in_weights_when_non_clustered(self):
+        """Test weights doesn't have 'cluster' key for non-clustered systems."""
+        fs = FlowSystem(timesteps=pd.date_range('2024-01-01', periods=24, freq='h'))
+        weights = fs.weights
+
+        # Non-clustered: 'cluster' not in weights
+        assert 'cluster' not in weights
 
-        assert 'time' in weights.temporal.dims
-        assert float(weights.temporal.sum().values) == 3.0
+    def test_temporal_dims_non_clustered(self):
+        """Test temporal_dims is ['time'] for non-clustered systems."""
+        fs = FlowSystem(timesteps=pd.date_range('2024-01-01', periods=24, freq='h'))
+
+        assert fs.temporal_dims == ['time']
 
-    def test_invalid_no_time_dim(self):
-        """Test error when temporal has no time dimension."""
-        temporal = xr.DataArray([1.0, 1.0], dims=['other'])
+    def test_temporal_weight(self):
+        """Test temporal_weight returns time * cluster."""
+        fs = FlowSystem(timesteps=pd.date_range('2024-01-01', periods=24, freq='h'))
 
-        with pytest.raises(ValueError, match='time'):
-            TimeSeriesWeights(temporal=temporal)
+        expected = fs.weights['time'] * fs.weights.get('cluster', 1.0)
+        xr.testing.assert_equal(fs.temporal_weight, expected)
 
-    def test_sum_over_time(self):
-        """Test sum_over_time convenience method."""
-        temporal = xr.DataArray([2.0, 3.0, 1.0], dims=['time'], coords={'time': [0, 1, 2]})
-        weights = TimeSeriesWeights(temporal=temporal)
+    def test_sum_temporal(self):
+        """Test sum_temporal applies full temporal weighting (time * cluster) and sums."""
+        fs = FlowSystem(timesteps=pd.date_range('2024-01-01', periods=3, freq='h'))
 
-        data = xr.DataArray([10.0, 20.0, 30.0], dims=['time'], coords={'time': [0, 1, 2]})
-        result = weights.sum_over_time(data)
+        # Input is a rate (e.g., flow_rate in MW)
+        data = xr.DataArray([10.0, 20.0, 30.0], dims=['time'], coords={'time': fs.timesteps})
 
-        # 10*2 + 20*3 + 30*1 = 20 + 60 + 30 = 110
-        assert float(result.values) == 110.0
+        result = fs.sum_temporal(data)
 
-    def test_effective_objective(self):
-        """Test effective_objective property."""
-        temporal = xr.DataArray([1.0, 1.0], dims=['time'])
-        objective = xr.DataArray([2.0, 2.0], dims=['time'])
+        # For hourly non-clustered: temporal = time * cluster = 1.0 * 1.0 = 1.0
+        # result = sum(data * temporal) = sum(data) = 60
+        assert float(result.values) == 60.0
 
-        # Without override
-        weights1 = TimeSeriesWeights(temporal=temporal)
-        assert np.array_equal(weights1.effective_objective.values, temporal.values)
 
-        # With override
-        weights2 = TimeSeriesWeights(temporal=temporal, objective=objective)
-        assert np.array_equal(weights2.effective_objective.values, objective.values)
+class TestFlowSystemDimsIndexesWeights:
+    """Tests for FlowSystem.dims, .indexes, .weights properties."""
 
+    def test_dims_property(self):
+        """Test that FlowSystem.dims returns active dimension names."""
+        fs = FlowSystem(timesteps=pd.date_range('2024-01-01', periods=24, freq='h'))
 
-class TestFlowSystemWeightsProperty:
-    """Tests for FlowSystem.weights property."""
+        dims = fs.dims
+        assert dims == ['time']
 
-    def test_weights_property_exists(self):
-        """Test that FlowSystem has weights property."""
+    def test_indexes_property(self):
+        """Test that FlowSystem.indexes returns active indexes."""
         fs = FlowSystem(timesteps=pd.date_range('2024-01-01', periods=24, freq='h'))
 
-        weights = fs.weights
-        assert isinstance(weights, TimeSeriesWeights)
+        indexes = fs.indexes
+        assert isinstance(indexes, dict)
+        assert 'time' in indexes
+        assert len(indexes['time']) == 24
 
-    def test_weights_temporal_equals_aggregation_weight(self):
-        """Test that weights.temporal equals aggregation_weight."""
+    def test_weights_keys_match_dims(self):
+        """Test that weights.keys() is subset of dims (only 'time' for simple case)."""
         fs = FlowSystem(timesteps=pd.date_range('2024-01-01', periods=24, freq='h'))
 
-        weights = fs.weights
-        aggregation_weight = fs.aggregation_weight
+        # For non-clustered, weights only has 'time'
+        assert set(fs.weights.keys()) == {'time'}
 
-        np.testing.assert_array_almost_equal(weights.temporal.values, aggregation_weight.values)
+    def test_temporal_weight_calculation(self):
+        """Test that temporal_weight = timestep_duration * cluster_weight."""
+        fs = FlowSystem(timesteps=pd.date_range('2024-01-01', periods=24, freq='h'))
+
+        expected = fs.timestep_duration * 1.0  # cluster is 1.0 for non-clustered
+
+        np.testing.assert_array_almost_equal(fs.temporal_weight.values, expected.values)
 
     def test_weights_with_cluster_weight(self):
-        """Test weights property includes cluster_weight."""
+        """Test weights property includes cluster_weight when provided."""
         # Create FlowSystem with custom cluster_weight
         timesteps = pd.date_range('2024-01-01', periods=24, freq='h')
-        cluster_weight = np.array([2.0] * 12 + [1.0] * 12)  # First 12h weighted 2x
+        cluster_weight = xr.DataArray(
+            np.array([2.0] * 12 + [1.0] * 12),
+            dims=['time'],
+            coords={'time': timesteps},
+        )
 
         fs = FlowSystem(timesteps=timesteps, cluster_weight=cluster_weight)
 
         weights = fs.weights
 
-        # temporal = timestep_duration * cluster_weight
-        # timestep_duration is 1h for all, so temporal = cluster_weight
+        # cluster weight should be in weights (FlowSystem has cluster_weight set)
+        # But note: 'cluster' only appears in weights if clusters dimension exists
+        # Since we didn't set clusters, 'cluster' won't be in weights
+        # The cluster_weight is applied via temporal_weight
+        assert 'cluster' not in weights  # No cluster dimension
+
+        # temporal_weight = timestep_duration * cluster_weight
+        # timestep_duration is 1h for all
         expected = 1.0 * cluster_weight
-        np.testing.assert_array_almost_equal(weights.temporal.values, expected)
+        np.testing.assert_array_almost_equal(fs.temporal_weight.values, expected.values)
 
 
 class TestClusterMethod:
diff --git a/tests/test_clustering_io.py b/tests/test_clustering_io.py
new file mode 100644
index 000000000..483cdc447
--- /dev/null
+++ b/tests/test_clustering_io.py
@@ -0,0 +1,241 @@
+"""Tests for clustering serialization and deserialization."""
+
+import numpy as np
+import pandas as pd
+import pytest
+
+import flixopt as fx
+
+
+@pytest.fixture
+def simple_system_24h():
+    """Create a simple flow system with 24 hourly timesteps."""
+    timesteps = pd.date_range('2023-01-01', periods=24, freq='h')
+
+    fs = fx.FlowSystem(timesteps)
+    fs.add_elements(
+        fx.Bus('heat'),
+        fx.Effect('costs', unit='EUR', description='costs', is_objective=True, is_standard=True),
+    )
+    fs.add_elements(
+        fx.Sink('demand', inputs=[fx.Flow('in', bus='heat', fixed_relative_profile=np.ones(24), size=10)]),
+        fx.Source('source', outputs=[fx.Flow('out', bus='heat', size=50, effects_per_flow_hour={'costs': 0.05})]),
+    )
+    return fs
+
+
+@pytest.fixture
+def simple_system_8_days():
+    """Create a simple flow system with 8 days of hourly timesteps."""
+    timesteps = pd.date_range('2023-01-01', periods=8 * 24, freq='h')
+
+    # Create varying demand profile with different patterns for different days
+    # 4 "weekdays" with high demand, 4 "weekend" days with low demand
+    hourly_pattern = np.sin(np.linspace(0, 2 * np.pi, 24)) * 0.5 + 0.5
+    weekday_profile = hourly_pattern * 1.5  # Higher demand
+    weekend_profile = hourly_pattern * 0.5  # Lower demand
+    demand_profile = np.concatenate(
+        [
+            weekday_profile,
+            weekday_profile,
+            weekday_profile,
+            weekday_profile,
+            weekend_profile,
+            weekend_profile,
+            weekend_profile,
+            weekend_profile,
+        ]
+    )
+
+    fs = fx.FlowSystem(timesteps)
+    fs.add_elements(
+        fx.Bus('heat'),
+        fx.Effect('costs', unit='EUR', description='costs', is_objective=True, is_standard=True),
+    )
+    fs.add_elements(
+        fx.Sink('demand', inputs=[fx.Flow('in', bus='heat', fixed_relative_profile=demand_profile, size=10)]),
+        fx.Source('source', outputs=[fx.Flow('out', bus='heat', size=50, effects_per_flow_hour={'costs': 0.05})]),
+    )
+    return fs
+
+
+class TestClusteringRoundtrip:
+    """Test that clustering survives dataset roundtrip."""
+
+    def test_clustering_to_dataset_has_clustering_attrs(self, simple_system_8_days):
+        """Clustered FlowSystem dataset should have clustering info."""
+        fs = simple_system_8_days
+        fs_clustered = fs.transform.cluster(n_clusters=2, cluster_duration='1D')
+
+        ds = fs_clustered.to_dataset(include_solution=False)
+
+        # Check that clustering attrs are present
+        assert 'clustering' in ds.attrs
+
+        # Check that clustering arrays are present with prefix
+        clustering_vars = [name for name in ds.data_vars if name.startswith('clustering|')]
+        assert len(clustering_vars) > 0
+
+    def test_clustering_roundtrip_preserves_clustering_object(self, simple_system_8_days):
+        """Clustering object should be restored after roundtrip."""
+        fs = simple_system_8_days
+        fs_clustered = fs.transform.cluster(n_clusters=2, cluster_duration='1D')
+
+        # Roundtrip
+        ds = fs_clustered.to_dataset(include_solution=False)
+        fs_restored = fx.FlowSystem.from_dataset(ds)
+
+        # Clustering should be restored
+        assert fs_restored.clustering is not None
+        assert fs_restored.clustering.backend_name == 'tsam'
+
+    def test_clustering_roundtrip_preserves_n_clusters(self, simple_system_8_days):
+        """Number of clusters should be preserved after roundtrip."""
+        fs = simple_system_8_days
+        fs_clustered = fs.transform.cluster(n_clusters=2, cluster_duration='1D')
+
+        ds = fs_clustered.to_dataset(include_solution=False)
+        fs_restored = fx.FlowSystem.from_dataset(ds)
+
+        assert fs_restored.clustering.n_clusters == 2
+
+    def test_clustering_roundtrip_preserves_timesteps_per_cluster(self, simple_system_8_days):
+        """Timesteps per cluster should be preserved after roundtrip."""
+        fs = simple_system_8_days
+        fs_clustered = fs.transform.cluster(n_clusters=2, cluster_duration='1D')
+
+        ds = fs_clustered.to_dataset(include_solution=False)
+        fs_restored = fx.FlowSystem.from_dataset(ds)
+
+        assert fs_restored.clustering.timesteps_per_cluster == 24
+
+    def test_clustering_roundtrip_preserves_original_timesteps(self, simple_system_8_days):
+        """Original timesteps should be preserved after roundtrip."""
+        fs = simple_system_8_days
+        fs_clustered = fs.transform.cluster(n_clusters=2, cluster_duration='1D')
+        original_timesteps = fs_clustered.clustering.original_timesteps
+
+        ds = fs_clustered.to_dataset(include_solution=False)
+        fs_restored = fx.FlowSystem.from_dataset(ds)
+
+        pd.testing.assert_index_equal(fs_restored.clustering.original_timesteps, original_timesteps)
+
+    def test_clustering_roundtrip_preserves_timestep_mapping(self, simple_system_8_days):
+        """Timestep mapping should be preserved after roundtrip."""
+        fs = simple_system_8_days
+        fs_clustered = fs.transform.cluster(n_clusters=2, cluster_duration='1D')
+        original_mapping = fs_clustered.clustering.timestep_mapping.values.copy()
+
+        ds = fs_clustered.to_dataset(include_solution=False)
+        fs_restored = fx.FlowSystem.from_dataset(ds)
+
+        np.testing.assert_array_equal(fs_restored.clustering.timestep_mapping.values, original_mapping)
+
+
+class TestClusteringWithSolutionRoundtrip:
+    """Test that clustering with solution survives roundtrip."""
+
+    def test_expand_solution_after_roundtrip(self, simple_system_8_days, solver_fixture):
+        """expand_solution should work after loading from dataset."""
+        fs = simple_system_8_days
+        fs_clustered = fs.transform.cluster(n_clusters=2, cluster_duration='1D')
+
+        # Solve
+        fs_clustered.optimize(solver_fixture)
+
+        # Roundtrip
+        ds = fs_clustered.to_dataset(include_solution=True)
+        fs_restored = fx.FlowSystem.from_dataset(ds)
+
+        # expand_solution should work
+        fs_expanded = fs_restored.transform.expand_solution()
+
+        # Check expanded FlowSystem has correct number of timesteps
+        assert len(fs_expanded.timesteps) == 8 * 24
+
+    def test_expand_solution_after_netcdf_roundtrip(self, simple_system_8_days, tmp_path, solver_fixture):
+        """expand_solution should work after loading from NetCDF file."""
+        fs = simple_system_8_days
+        fs_clustered = fs.transform.cluster(n_clusters=2, cluster_duration='1D')
+
+        # Solve
+        fs_clustered.optimize(solver_fixture)
+
+        # Save to NetCDF
+        nc_path = tmp_path / 'clustered.nc'
+        fs_clustered.to_netcdf(nc_path)
+
+        # Load from NetCDF
+        fs_restored = fx.FlowSystem.from_netcdf(nc_path)
+
+        # expand_solution should work
+        fs_expanded = fs_restored.transform.expand_solution()
+
+        # Check expanded FlowSystem has correct number of timesteps
+        assert len(fs_expanded.timesteps) == 8 * 24
+
+
+class TestClusteringDerivedProperties:
+    """Test derived properties on Clustering object."""
+
+    def test_original_timesteps_property(self, simple_system_8_days):
+        """original_timesteps property should return correct DatetimeIndex."""
+        fs = simple_system_8_days
+        original_timesteps = fs.timesteps
+
+        fs_clustered = fs.transform.cluster(n_clusters=2, cluster_duration='1D')
+
+        # Check values are equal (name attribute may differ)
+        pd.testing.assert_index_equal(
+            fs_clustered.clustering.original_timesteps,
+            original_timesteps,
+            check_names=False,
+        )
+
+    def test_simple_system_has_no_periods_or_scenarios(self, simple_system_8_days):
+        """Clustered simple system should preserve that it has no periods/scenarios."""
+        fs = simple_system_8_days
+        fs_clustered = fs.transform.cluster(n_clusters=2, cluster_duration='1D')
+
+        # FlowSystem without periods/scenarios should remain so after clustering
+        assert fs_clustered.periods is None
+        assert fs_clustered.scenarios is None
+
+
+class TestClusteringWithScenarios:
+    """Test clustering IO with scenarios."""
+
+    @pytest.fixture
+    def system_with_scenarios(self):
+        """Create a flow system with scenarios."""
+        timesteps = pd.date_range('2023-01-01', periods=4 * 24, freq='h')
+        scenarios = pd.Index(['Low', 'High'], name='scenario')
+
+        # Create varying demand profile for clustering
+        demand_profile = np.tile(np.sin(np.linspace(0, 2 * np.pi, 24)) * 0.5 + 0.5, 4)
+
+        fs = fx.FlowSystem(timesteps, scenarios=scenarios)
+        fs.add_elements(
+            fx.Bus('heat'),
+            fx.Effect('costs', unit='EUR', description='costs', is_objective=True, is_standard=True),
+        )
+        fs.add_elements(
+            fx.Sink('demand', inputs=[fx.Flow('in', bus='heat', fixed_relative_profile=demand_profile, size=10)]),
+            fx.Source('source', outputs=[fx.Flow('out', bus='heat', size=50, effects_per_flow_hour={'costs': 0.05})]),
+        )
+        return fs
+
+    def test_clustering_roundtrip_preserves_scenarios(self, system_with_scenarios):
+        """Scenarios should be preserved after clustering and roundtrip."""
+        fs = system_with_scenarios
+        fs_clustered = fs.transform.cluster(n_clusters=2, cluster_duration='1D')
+
+        ds = fs_clustered.to_dataset(include_solution=False)
+        fs_restored = fx.FlowSystem.from_dataset(ds)
+
+        # Scenarios should be preserved in the FlowSystem itself
+        pd.testing.assert_index_equal(
+            fs_restored.scenarios,
+            pd.Index(['Low', 'High'], name='scenario'),
+            check_names=False,
+        )
diff --git a/tests/test_io_conversion.py b/tests/test_io_conversion.py
index 33bda8c91..7775f987a 100644
--- a/tests/test_io_conversion.py
+++ b/tests/test_io_conversion.py
@@ -762,6 +762,11 @@ def test_v4_reoptimized_objective_matches_original(self, result_name):
         new_objective = float(fs.solution['objective'].item())
         new_effect_total = float(fs.solution[objective_effect_label].sum().item())
 
+        # Skip comparison for scenarios test case - scenario weights are now always normalized,
+        # which changes the objective value when loading old results with non-normalized weights
+        if result_name == '04_scenarios':
+            pytest.skip('Scenario weights are now always normalized - old results have different weights')
+
         # Verify objective matches (within tolerance)
         assert new_objective == pytest.approx(old_objective, rel=1e-5, abs=1), (
             f'Objective mismatch for {result_name}: new={new_objective}, old={old_objective}'
diff --git a/tests/test_scenarios.py b/tests/test_scenarios.py
index 65ea62d81..2699647ad 100644
--- a/tests/test_scenarios.py
+++ b/tests/test_scenarios.py
@@ -341,12 +341,14 @@ def test_scenarios_selection(flow_system_piecewise_conversion_scenarios):
 
     assert flow_system.scenarios.equals(flow_system_full.scenarios[0:2])
 
-    np.testing.assert_allclose(flow_system.scenario_weights.values, flow_system_full.scenario_weights[0:2])
+    # Scenario weights are always normalized - subset is re-normalized to sum to 1
+    subset_weights = flow_system_full.scenario_weights[0:2]
+    expected_normalized = subset_weights / subset_weights.sum()
+    np.testing.assert_allclose(flow_system.scenario_weights.values, expected_normalized.values)
 
-    # Optimize using new API with normalize_weights=False
+    # Optimize using new API
     flow_system.optimize(
         fx.solvers.GurobiSolver(mip_gap=0.01, time_limit_seconds=60),
-        normalize_weights=False,
     )
 
     # Penalty has same structure as other effects: 'Penalty' is the total, 'Penalty(temporal)' and 'Penalty(periodic)' are components
@@ -769,7 +771,10 @@ def test_weights_selection():
 
     # Verify weights are correctly sliced
     assert fs_subset.scenarios.equals(pd.Index(['base', 'high'], name='scenario'))
-    np.testing.assert_allclose(fs_subset.scenario_weights.values, custom_scenario_weights[[0, 2]])
+    # Scenario weights are always normalized - subset is re-normalized to sum to 1
+    subset_weights = np.array([0.3, 0.2])  # Original weights for selected scenarios
+    expected_normalized = subset_weights / subset_weights.sum()
+    np.testing.assert_allclose(fs_subset.scenario_weights.values, expected_normalized)
 
     # Verify weights are 1D with just scenario dimension (no period dimension)
     assert fs_subset.scenario_weights.dims == ('scenario',)
diff --git a/tests/test_sel_isel_single_selection.py b/tests/test_sel_isel_single_selection.py
new file mode 100644
index 000000000..4d84ced51
--- /dev/null
+++ b/tests/test_sel_isel_single_selection.py
@@ -0,0 +1,193 @@
+"""Tests for sel/isel with single period/scenario selection."""
+
+import numpy as np
+import pandas as pd
+import pytest
+
+import flixopt as fx
+
+
+@pytest.fixture
+def fs_with_scenarios():
+    """FlowSystem with scenarios for testing single selection."""
+    timesteps = pd.date_range('2023-01-01', periods=24, freq='h')
+    scenarios = pd.Index(['A', 'B', 'C'], name='scenario')
+    scenario_weights = np.array([0.5, 0.3, 0.2])
+
+    fs = fx.FlowSystem(timesteps, scenarios=scenarios, scenario_weights=scenario_weights)
+    fs.add_elements(
+        fx.Bus('heat'),
+        fx.Effect('costs', unit='EUR', description='costs', is_objective=True, is_standard=True),
+    )
+    fs.add_elements(
+        fx.Sink('demand', inputs=[fx.Flow('in', bus='heat', fixed_relative_profile=np.ones(24), size=10)]),
+        fx.Source('source', outputs=[fx.Flow('out', bus='heat', size=50, effects_per_flow_hour={'costs': 0.05})]),
+    )
+    return fs
+
+
+@pytest.fixture
+def fs_with_periods():
+    """FlowSystem with periods for testing single selection."""
+    timesteps = pd.date_range('2023-01-01', periods=24, freq='h')
+    periods = pd.Index([2020, 2030, 2040], name='period')
+
+    fs = fx.FlowSystem(timesteps, periods=periods, weight_of_last_period=10)
+    fs.add_elements(
+        fx.Bus('heat'),
+        fx.Effect('costs', unit='EUR', description='costs', is_objective=True, is_standard=True),
+    )
+    fs.add_elements(
+        fx.Sink('demand', inputs=[fx.Flow('in', bus='heat', fixed_relative_profile=np.ones(24), size=10)]),
+        fx.Source('source', outputs=[fx.Flow('out', bus='heat', size=50, effects_per_flow_hour={'costs': 0.05})]),
+    )
+    return fs
+
+
+@pytest.fixture
+def fs_with_periods_and_scenarios():
+    """FlowSystem with both periods and scenarios."""
+    timesteps = pd.date_range('2023-01-01', periods=24, freq='h')
+    periods = pd.Index([2020, 2030], name='period')
+    scenarios = pd.Index(['Low', 'High'], name='scenario')
+
+    fs = fx.FlowSystem(timesteps, periods=periods, scenarios=scenarios, weight_of_last_period=10)
+    fs.add_elements(
+        fx.Bus('heat'),
+        fx.Effect('costs', unit='EUR', description='costs', is_objective=True, is_standard=True),
+    )
+    fs.add_elements(
+        fx.Sink('demand', inputs=[fx.Flow('in', bus='heat', fixed_relative_profile=np.ones(24), size=10)]),
+        fx.Source('source', outputs=[fx.Flow('out', bus='heat', size=50, effects_per_flow_hour={'costs': 0.05})]),
+    )
+    return fs
+
+
+class TestIselSingleScenario:
+    """Test isel with single scenario selection."""
+
+    def test_isel_single_scenario_drops_dimension(self, fs_with_scenarios):
+        """Selecting a single scenario with isel should drop the scenario dimension."""
+        fs_selected = fs_with_scenarios.transform.isel(scenario=0)
+
+        assert fs_selected.scenarios is None
+        assert 'scenario' not in fs_selected.to_dataset().dims
+
+    def test_isel_single_scenario_removes_scenario_weights(self, fs_with_scenarios):
+        """scenario_weights should be removed when scenario dimension is dropped."""
+        fs_selected = fs_with_scenarios.transform.isel(scenario=0)
+
+        ds = fs_selected.to_dataset()
+        assert 'scenario_weights' not in ds.data_vars
+        assert 'scenario_weights' not in ds.attrs
+
+    def test_isel_single_scenario_preserves_time(self, fs_with_scenarios):
+        """Time dimension should be preserved."""
+        fs_selected = fs_with_scenarios.transform.isel(scenario=0)
+
+        assert len(fs_selected.timesteps) == 24
+
+    def test_isel_single_scenario_roundtrip(self, fs_with_scenarios):
+        """FlowSystem should survive to_dataset/from_dataset roundtrip after single selection."""
+        fs_selected = fs_with_scenarios.transform.isel(scenario=0)
+
+        ds = fs_selected.to_dataset()
+        fs_restored = fx.FlowSystem.from_dataset(ds)
+
+        assert fs_restored.scenarios is None
+        assert len(fs_restored.timesteps) == 24
+
+
+class TestSelSingleScenario:
+    """Test sel with single scenario selection."""
+
+    def test_sel_single_scenario_drops_dimension(self, fs_with_scenarios):
+        """Selecting a single scenario with sel should drop the scenario dimension."""
+        fs_selected = fs_with_scenarios.transform.sel(scenario='B')
+
+        assert fs_selected.scenarios is None
+
+
+class TestIselSinglePeriod:
+    """Test isel with single period selection."""
+
+    def test_isel_single_period_drops_dimension(self, fs_with_periods):
+        """Selecting a single period with isel should drop the period dimension."""
+        fs_selected = fs_with_periods.transform.isel(period=0)
+
+        assert fs_selected.periods is None
+        assert 'period' not in fs_selected.to_dataset().dims
+
+    def test_isel_single_period_removes_period_weights(self, fs_with_periods):
+        """period_weights should be removed when period dimension is dropped."""
+        fs_selected = fs_with_periods.transform.isel(period=0)
+
+        ds = fs_selected.to_dataset()
+        assert 'period_weights' not in ds.data_vars
+        assert 'weight_of_last_period' not in ds.attrs
+
+    def test_isel_single_period_roundtrip(self, fs_with_periods):
+        """FlowSystem should survive roundtrip after single period selection."""
+        fs_selected = fs_with_periods.transform.isel(period=0)
+
+        ds = fs_selected.to_dataset()
+        fs_restored = fx.FlowSystem.from_dataset(ds)
+
+        assert fs_restored.periods is None
+
+
+class TestSelSinglePeriod:
+    """Test sel with single period selection."""
+
+    def test_sel_single_period_drops_dimension(self, fs_with_periods):
+        """Selecting a single period with sel should drop the period dimension."""
+        fs_selected = fs_with_periods.transform.sel(period=2030)
+
+        assert fs_selected.periods is None
+
+
+class TestMixedSelection:
+    """Test mixed selections (single + multiple)."""
+
+    def test_single_period_multiple_scenarios(self, fs_with_periods_and_scenarios):
+        """Single period but multiple scenarios should only drop period."""
+        fs_selected = fs_with_periods_and_scenarios.transform.isel(period=0)
+
+        assert fs_selected.periods is None
+        assert fs_selected.scenarios is not None
+        assert len(fs_selected.scenarios) == 2
+
+    def test_multiple_periods_single_scenario(self, fs_with_periods_and_scenarios):
+        """Multiple periods but single scenario should only drop scenario."""
+        fs_selected = fs_with_periods_and_scenarios.transform.isel(scenario=0)
+
+        assert fs_selected.periods is not None
+        assert len(fs_selected.periods) == 2
+        assert fs_selected.scenarios is None
+
+    def test_single_period_single_scenario(self, fs_with_periods_and_scenarios):
+        """Single period and single scenario should drop both."""
+        fs_selected = fs_with_periods_and_scenarios.transform.isel(period=0, scenario=0)
+
+        assert fs_selected.periods is None
+        assert fs_selected.scenarios is None
+
+
+class TestSliceSelection:
+    """Test that slice selection preserves dimensions."""
+
+    def test_slice_scenarios_preserves_dimension(self, fs_with_scenarios):
+        """Slice selection should preserve dimension even with 1 element."""
+        # Select a slice that results in 2 elements
+        fs_selected = fs_with_scenarios.transform.isel(scenario=slice(0, 2))
+
+        assert fs_selected.scenarios is not None
+        assert len(fs_selected.scenarios) == 2
+
+    def test_list_selection_preserves_dimension(self, fs_with_scenarios):
+        """List selection should preserve dimension even with 1 element."""
+        fs_selected = fs_with_scenarios.transform.isel(scenario=[0])
+
+        # List selection should preserve dimension
+        assert fs_selected.scenarios is not None
+        assert len(fs_selected.scenarios) == 1