diff --git a/docs/design/cluster_architecture.md b/docs/design/cluster_architecture.md new file mode 100644 index 000000000..e90aeccff --- /dev/null +++ b/docs/design/cluster_architecture.md @@ -0,0 +1,773 @@ +# Design Document: Cluster Architecture for flixopt + +## Executive Summary + +This document defines the architecture for cluster representation in flixopt using **true `(cluster, time)` dimensions**. + +### Key Decision: True Dimensions (Option B) + +**Chosen Approach:** +```python +# Clustered data structure: +data.dims = ('cluster', 'time', 'period', 'scenario') +data.shape = (9, 24, ...) # 9 clusters × 24 timesteps each +``` + +**Why True Dimensions?** +1. **Temporal constraints just work** - `x[:, 1:] - x[:, :-1]` naturally stays within clusters +2. **No boundary masking** - StorageModel, StatusModel constraints are clean and vectorized +3. **Plotting trivial** - existing `facet_col='cluster'` works automatically + +### Document Scope + +1. Current architecture analysis (Part 1) +2. Architectural options and recommendation (Part 2) +3. Impact on Features - StatusModel, StorageModel, etc. (Part 3) +4. Plotting improvements (Part 4) +5. Future considerations (Part 5) +6. Implementation roadmap (Part 6) + +--- + +## Part 1: Current Architecture Analysis + +### 1.1 Time Dimension Structure + +**Current Implementation:** +``` +time: (n_clusters × timesteps_per_cluster,) # Flat, e.g., (864,) for 9 clusters × 96 timesteps +``` + +**Key Properties:** +- `cluster_weight`: Shape `(time,)` with repeated values per cluster +- `timestep_duration`: Shape `(time,)` or scalar +- `aggregation_weight = timestep_duration × cluster_weight` + +**Cluster Tracking:** +- `cluster_start_positions`: Array of indices where each cluster begins +- `ClusterStructure`: Stores cluster_order, occurrences, n_clusters, timesteps_per_cluster + +### 1.2 Features Affected by Time Structure + +| Feature | Time Usage | Clustering Impact | +|---------|-----------|-------------------| +| **StatusModel** | `aggregation_weight` for active hours, `timestep_duration` for effects | Must sum correctly across clusters | +| **InvestmentModel** | Periodic (no time dim) | Unaffected by time structure | +| **PiecewiseModel** | Per-timestep lambda variables | Must preserve cluster structure | +| **ShareAllocationModel** | Uses `cluster_weight` explicitly | Directly depends on weight structure | +| **StorageModel** | Charge balance across time | Needs cluster boundary handling | +| **InterclusterStorageModel** | SOC_boundary linking | Uses cluster indices extensively | + +### 1.3 Current Plotting Structure + +**StatisticsPlotAccessor Methods:** +- `balance()`: Node flow visualization +- `storage()`: Dual-axis charge/discharge + SOC +- `heatmap()`: 2D time reshaping (days × hours) +- `duration_curve()`: Sorted load profiles +- `effects()`: Cost/emission breakdown + +**Clustering in Plots:** +- `ClusterStructure.plot()`: Shows cluster assignments +- Cluster weight applied when aggregating (`cluster_weight.sum('time')`) +- No visual separation between clusters in time series plots + +--- + +## Part 2: Architectural Options + +### 2.1 Option A: Enhanced Flat with Indexers + +Keep flat `time` dimension, add xarray indexer properties. + +**Pros:** Supports variable-length clusters +**Cons:** Every temporal constraint needs explicit boundary masking + +**NOT RECOMMENDED** - see Option B. + +### 2.2 Option B: True (cluster, time) Dimensions (RECOMMENDED) + +Reshape time to 2D when clustering is active: + +```python +# ═══════════════════════════════════════════════════════════════ +# DIMENSION STRUCTURE +# ═══════════════════════════════════════════════════════════════ + +# Non-clustered: +data.dims = ('time', 'period', 'scenario') +data.shape = (8760, ...) # Full year hourly + +# Clustered: +data.dims = ('cluster', 'time', 'period', 'scenario') +data.shape = (9, 24, ...) # 9 clusters × 24 timesteps each + +# Varying segment durations supported: +timestep_duration.dims = ('cluster', 'time') +timestep_duration.shape = (9, 24) # Different durations per segment per cluster +``` + +**Key Benefits - Temporal Constraints Just Work!** + +```python +# ═══════════════════════════════════════════════════════════════ +# STORAGE: Charge balance naturally within clusters +# ═══════════════════════════════════════════════════════════════ +# charge_state shape: (cluster, time+1, period, scenario) - extra timestep for boundaries +charge_state = ... # (9, 25, ...) + +# Balance constraint - NO MASKING NEEDED! +lhs = charge_state[:, 1:] - charge_state[:, :-1] * (1 - loss) - charge + discharge +# Shape: (cluster, time, period, scenario) = (9, 24, ...) + +# Delta per cluster (for inter-cluster linking): +delta_soc = charge_state[:, -1] - charge_state[:, 0] # Shape: (cluster, ...) + +# ═══════════════════════════════════════════════════════════════ +# STATUS: Uptime/downtime constraints stay within clusters +# ═══════════════════════════════════════════════════════════════ +status = ... # (cluster, time, ...) + +# State transitions - naturally per cluster! +activate = status[:, 1:] - status[:, :-1] # No boundary issues! + +# min_uptime constraint - works correctly, can't span clusters + +# ═══════════════════════════════════════════════════════════════ +# INTER-CLUSTER OPERATIONS +# ═══════════════════════════════════════════════════════════════ +# Select first/last timestep of each cluster: +at_start = data.isel(time=0) # Shape: (cluster, period, scenario) +at_end = data.isel(time=-1) # Shape: (cluster, period, scenario) + +# Compute per-cluster statistics: +mean_per_cluster = data.mean(dim='time') +max_per_cluster = data.max(dim='time') +``` + +**Varying Segment Durations (Future Segmentation):** + +```python +# Same NUMBER of segments per cluster, different DURATIONS: +timestep_duration = xr.DataArray( + [ + [2, 2, 1, 1, 2, 4], # Cluster 0: segments sum to 12h + [1, 3, 2, 2, 2, 2], # Cluster 1: segments sum to 12h + ... + ], + dims=['cluster', 'time'], + coords={'cluster': range(9), 'time': range(6)} +) + +# aggregation_weight still works: +aggregation_weight = timestep_duration * cluster_weight # (cluster, time) * (cluster,) +``` + +**Pros:** +- Temporal constraints naturally stay within clusters - NO MASKING! +- StatusModel uptime/downtime just works +- Storage balance is clean +- Much less code, fewer bugs +- Supports varying segment durations (same count, different lengths) + +**Cons:** +- More upfront refactoring +- All code paths need to handle `(cluster, time)` vs `(time,)` based on `is_clustered` + +### 2.3 Recommendation: Option B (True Dimensions) + +Given: +- Uniform timestep COUNT per cluster (tsam default) +- Variable segment DURATIONS supported via `timestep_duration[cluster, time]` +- Much cleaner constraint handling + +**Option B is the recommended choice.** + +--- + +## Part 3: Impact on Features + +### 3.1 StatusModel Impact - SOLVED BY TRUE DIMENSIONS + +**With `(cluster, time)` dimensions, temporal constraints naturally stay within clusters!** + +```python +status.dims = ('cluster', 'time', 'period', 'scenario') +status.shape = (9, 24, ...) + +# State transitions - per cluster, no boundary issues! +activate = status[:, 1:] - status[:, :-1] + +# min_uptime constraint operates within each cluster's time dimension +# Cannot accidentally span cluster boundaries +``` + +**What works automatically:** +- ✅ `min_uptime`, `min_downtime` - constraints stay within clusters +- ✅ `initial_status` - applies to each cluster's first timestep +- ✅ State transitions - naturally per cluster +- ✅ `active_hours` - uses `aggregation_weight` correctly +- ✅ `effects_per_startup` - counted per cluster, weighted by `cluster_weight` + +**Optional Enhancement - cluster_mode for special cases:** +```python +class StatusParameters: + # NEW: How to handle cluster boundaries (default: independent) + cluster_mode: Literal['independent', 'cyclic'] = 'independent' +``` + +| Mode | Behavior | +|------|----------| +| `independent` | Each cluster starts fresh (default, most common) | +| `cyclic` | `status[:, 0] == status[:, -1]` - status returns to start | + +### 3.2 StorageModel Impact - SIMPLIFIED + +**With `(cluster, time)` dimensions, storage constraints become trivial:** + +```python +charge_state.dims = ('cluster', 'time_extra', 'period', 'scenario') +charge_state.shape = (9, 25, ...) # 24 timesteps + 1 boundary per cluster + +# ═══════════════════════════════════════════════════════════════ +# Charge balance - NO MASKING! +# ═══════════════════════════════════════════════════════════════ +lhs = ( + charge_state[:, 1:] - + charge_state[:, :-1] * (1 - loss_rate) - + charge * eta_charge + + discharge / eta_discharge +) +self.add_constraints(lhs == 0, name='charge_balance') # Clean! + +# ═══════════════════════════════════════════════════════════════ +# Delta SOC per cluster (for inter-cluster linking) +# ═══════════════════════════════════════════════════════════════ +delta_soc = charge_state[:, -1] - charge_state[:, 0] # Shape: (cluster, ...) + +# ═══════════════════════════════════════════════════════════════ +# Cluster start constraint (relative SOC starts at 0) +# ═══════════════════════════════════════════════════════════════ +self.add_constraints(charge_state[:, 0] == 0, name='cluster_start') + +# ═══════════════════════════════════════════════════════════════ +# Cyclic constraint (optional) +# ═══════════════════════════════════════════════════════════════ +self.add_constraints(charge_state[:, 0] == charge_state[:, -1], name='cyclic') +``` + +**InterclusterStorageModel also simplified** - SOC_boundary linking uses clean slicing. + +### 3.3 ShareAllocationModel Impact + +**Current Code (features.py:624):** +```python +self._eq_total.lhs -= (self.total_per_timestep * self._model.cluster_weight).sum(dim='time') +``` + +**With Enhanced Helpers:** +No changes needed - `cluster_weight` structure preserved. + +### 3.4 PiecewiseModel Impact + +**Current Code:** Creates lambda variables per timestep. + +**With Enhanced Helpers:** +No changes needed - operates on flat time dimension. + +### 3.5 Summary: Models with True (cluster, time) Dimensions + +| Model | Cross-Timestep Constraints | With True Dims | Action Needed | +|-------|---------------------------|----------------|---------------| +| **StorageModel** | `cs[t] - cs[t-1]` | ✅ Just works | Simplify code | +| **StatusModel** | min_uptime, min_downtime | ✅ Just works | Optional cluster_mode | +| **consecutive_duration_tracking** | State machine | ✅ Just works | No changes | +| **state_transition_bounds** | `activate[t] - status[t-1]` | ✅ Just works | No changes | +| **PiecewiseModel** | Per-timestep only | ✅ Just works | No changes | +| **ShareAllocationModel** | Sum with cluster_weight | ✅ Just works | No changes | +| **InvestmentModel** | No time dimension | ✅ Just works | No changes | + +**Key Insight:** With true `(cluster, time)` dimensions, `x[:, 1:] - x[:, :-1]` naturally stays within clusters! + +--- + +## Part 4: Plotting Improvements + +### 4.1 Key Benefit of True Dimensions: Minimal Plotting Changes + +With true `(cluster, time)` dimensions, plotting becomes trivial because: +1. Data already has the right shape - no reshaping needed +2. Existing `facet_col='cluster'` parameter just works +3. Only minimal changes needed: auto-add cluster separators in combined views + +### 4.2 Proposed Approach: Leverage Existing Infrastructure + +#### 4.2.1 Use Existing facet_col Parameter + +**No new plot methods needed!** The existing infrastructure handles `cluster` dimension: + +```python +# ═══════════════════════════════════════════════════════════════ +# EXISTING API - works automatically with (cluster, time) dims! +# ═══════════════════════════════════════════════════════════════ +fs.statistics.plot.storage('Battery', facet_col='cluster') # One subplot per cluster +fs.statistics.plot.balance('Heat', facet_col='cluster') # One subplot per cluster +fs.statistics.plot.flows(..., facet_col='cluster') # Same pattern + +# Combine with other dimensions +fs.statistics.plot.balance('Heat', facet_col='cluster', facet_row='scenario') +``` + +#### 4.2.2 Auto-Add Cluster Separators (Small Change) + +For combined views (no faceting), add visual separators: + +```python +def _create_base_plot(self, data, **kwargs): + """Base plot creation - add cluster separators if combined view.""" + fig = ... # existing logic + + # Auto-add cluster separators if clustered and showing combined time + if self._fs.is_clustered and 'cluster' not in kwargs.get('facet_col', ''): + # Add subtle vertical lines between clusters + for cluster_idx in range(1, self._fs.clustering.n_clusters): + x_pos = cluster_idx * self._fs.clustering.timesteps_per_cluster + fig.add_vline(x=x_pos, line_dash='dot', opacity=0.3, line_color='gray') + + return fig +``` + +#### 4.2.3 Per-Cluster Statistics (Natural with True Dims) + +With `(cluster, time)` dimensions, aggregation is trivial: + +```python +# Mean per cluster - just use xarray +mean_per_cluster = data.mean(dim='time') # Shape: (cluster, ...) +max_per_cluster = data.max(dim='time') + +# Can plot directly +fs.statistics.plot.bar(data.mean('time'), x='cluster', title='Mean by Cluster') +``` + +#### 4.2.4 Heatmap (Already Correct Shape) + +With true dimensions, heatmaps work directly: + +```python +# Data already has (cluster, time) shape - heatmap just works! +def cluster_heatmap(self, variable): + data = self._get_variable(variable) + + # With (cluster, time) dims, no reshaping needed! + return self._plot_heatmap( + data, # Already (cluster, time, ...) + x='time', + y='cluster', + colorbar_title=variable + ) +``` + +### 4.3 Summary: Plotting Changes Required + +| Change | Scope | Complexity | +|--------|-------|------------| +| Auto cluster separators in base plot | ~10 lines in `_create_base_plot` | Low | +| Ensure facet_col='cluster' works | Should work already | None | +| Heatmap with cluster dim | Works automatically | None | +| No new plot methods needed | - | - | + +--- + +## Part 5: Future Considerations + +### 5.1 Variable Segment Durations (Out of Scope) + +tsam supports intra-period segmentation with variable segment durations. This could be supported in the future via: +- Integer-based `time` index (0, 1, 2, ...) instead of timestamps +- `timestep_duration[cluster, time]` array for variable durations per segment + +**Not implemented in initial version** - the architecture supports it, but it's not a priority. + +--- + +## Part 6: Implementation Roadmap + +### Phase 1: Core Dimension Refactoring (PRIORITY) + +**Goal:** Introduce true `(cluster, time)` dimensions throughout the codebase. + +**Tasks:** +1. Update `FlowSystem` to support `(cluster, time)` dimension structure when clustered +2. Add `is_clustered` property to `FlowSystem` +3. Update `Clustering` class with: + - `n_clusters: int` property + - `timesteps_per_cluster: int` property + - Coordinate accessors for cluster dimension +4. Update `cluster_weight` to have shape `(cluster,)` instead of `(time,)` +5. Update `timestep_duration` to have shape `(cluster, time)` when clustered +6. Update `aggregation_weight` computation to broadcast correctly + +**Files:** +- `flixopt/flow_system.py` - Core dimension handling +- `flixopt/clustering/base.py` - Updated Clustering class + +**Key Changes:** +```python +# FlowSystem property updates: +@property +def is_clustered(self) -> bool: + return self.clustering is not None + +@property +def cluster_weight(self) -> xr.DataArray: + if not self.is_clustered: + return xr.DataArray(1.0) + # Shape: (cluster,) - one weight per cluster + return xr.DataArray( + self.clustering.cluster_occurrences, + dims=['cluster'], + coords={'cluster': range(self.clustering.n_clusters)} + ) + +@property +def timestep_duration(self) -> xr.DataArray: + if not self.is_clustered: + return self._timestep_duration # Shape: (time,) or scalar + # Shape: (cluster, time) when clustered + return self._timestep_duration # Already 2D from clustering + +@property +def aggregation_weight(self) -> xr.DataArray: + return self.timestep_duration * self.cluster_weight # Broadcasting handles shapes +``` + +### Phase 2: Update Variable/Constraint Creation + +**Goal:** All variables and constraints use `(cluster, time)` dimensions when clustered. + +**Tasks:** +1. Update `create_variable` to use `(cluster, time, period, scenario)` dims when clustered +2. Update constraint generation in all models +3. Verify linopy handles multi-dimensional constraint arrays correctly +4. Add tests for both clustered and non-clustered paths + +**Files:** +- `flixopt/core.py` - Variable creation +- `flixopt/components.py` - StorageModel, other component models +- `flixopt/features.py` - StatusModel, other feature models + +**Key Pattern:** +```python +# Dimension-aware variable creation: +def _get_time_dims(self) -> list[str]: + if self.flow_system.is_clustered: + return ['cluster', 'time'] + return ['time'] + +def _get_time_coords(self) -> dict: + if self.flow_system.is_clustered: + return { + 'cluster': range(self.flow_system.clustering.n_clusters), + 'time': range(self.flow_system.clustering.timesteps_per_cluster) + } + return {'time': self.flow_system.time_coords} +``` + +### Phase 3: Simplify StorageModel and InterclusterStorageModel + +**Goal:** Leverage true dimensions for clean constraint generation. + +**Tasks:** +1. Simplify `StorageModel.charge_balance` - no boundary masking needed +2. Simplify delta SOC calculation: `charge_state[:, -1] - charge_state[:, 0]` +3. Simplify `InterclusterStorageModel` linking constraints +4. Update `intercluster_helpers.py` utilities + +**Files:** +- `flixopt/components.py` - StorageModel, InterclusterStorageModel +- `flixopt/clustering/intercluster_helpers.py` - Simplified helpers + +**Before/After:** +```python +# BEFORE (flat time with masking): +start_positions = clustering.cluster_start_positions +end_positions = start_positions[1:] - 1 +mask = _build_boundary_mask(...) +balance = charge_state.isel(time=slice(1, None)).where(~mask) - ... + +# AFTER (true dimensions): +# charge_state shape: (cluster, time+1, ...) +balance = ( + charge_state[:, 1:] - + charge_state[:, :-1] * (1 - loss_rate) - + charge * eta_charge + + discharge / eta_discharge +) +# No masking needed - constraints naturally stay within clusters! +``` + +### Phase 4: Update transform_accessor.cluster() + +**Goal:** Produce true `(cluster, time)` shaped data. + +**Tasks:** +1. Update `cluster()` to reshape time series to `(cluster, time)` +2. Generate proper coordinates for cluster dimension +3. Update `expand_solution()` to handle reverse transformation +4. Handle SOC_boundary expansion for inter-cluster storage + +**Files:** +- `flixopt/transform_accessor.py` - cluster() and expand_solution() + +**Key Implementation:** +```python +def cluster(self, n_clusters, cluster_duration, ...): + """Create clustered FlowSystem with (cluster, time) dimensions.""" + ... + # Reshape all time series: (flat_time,) → (cluster, time) + for key, ts in time_series.items(): + reshaped = ts.values.reshape(n_clusters, timesteps_per_cluster) + new_ts = xr.DataArray( + reshaped, + dims=['cluster', 'time'], + coords={'cluster': range(n_clusters), 'time': range(timesteps_per_cluster)} + ) + clustered_time_series[key] = new_ts + ... + +def expand_solution(self): + """Expand clustered solution back to original timeline.""" + expanded = {} + for var_name, var_data in self.solution.items(): + if 'cluster' in var_data.dims: + # Expand using cluster_order to map back to original periods + expanded[var_name] = self._expand_clustered_data(var_data) + else: + expanded[var_name] = var_data + return xr.Dataset(expanded) +``` + +### Phase 5: Plotting Integration + +**Goal:** Minimal changes - leverage existing infrastructure. + +**Tasks:** +1. Ensure `facet_col='cluster'` works with existing plot methods +2. Add auto cluster separators in combined time series views +3. Test heatmaps with `(cluster, time)` data + +**Files:** +- `flixopt/statistics_accessor.py` - Minor update to base plot method + +**Implementation:** +```python +# In _create_base_plot or similar: +def _add_cluster_separators(self, fig): + """Add subtle separators between clusters in combined view.""" + if self._fs.is_clustered: + for cluster_idx in range(1, self._fs.clustering.n_clusters): + x_pos = cluster_idx * self._fs.clustering.timesteps_per_cluster + fig.add_vline(x=x_pos, line_dash='dot', opacity=0.3) +``` + +### Phase Summary + +| Phase | Goal | Complexity | StatusModel Fix? | +|-------|------|------------|------------------| +| 1 | Core dimension refactoring | High | N/A (prep work) | +| 2 | Variable/constraint creation | Medium | ✅ Automatic | +| 3 | StorageModel simplification | Medium | N/A | +| 4 | transform_accessor updates | Medium | N/A | +| 5 | Plotting integration | Low | N/A | + +**Key Insight:** With true `(cluster, time)` dimensions, StatusModel and other temporal constraints **just work** without any special handling. The dimension structure naturally prevents constraints from spanning cluster boundaries. + +--- + +## Part 7: Testing Strategy + +### 7.1 Unit Tests + +```python +# Test cluster helpers +def test_cluster_labels_uniform(): + """Verify cluster_labels for uniform cluster lengths.""" + +def test_cluster_slices_variable(): + """Verify cluster_slices for variable cluster lengths.""" + +def test_boundaries_vary_by_period(): + """Verify boundary dispatch for different periods.""" +``` + +### 7.2 Integration Tests + +```python +# Test storage with different cluster modes +def test_storage_intercluster_with_helpers(): + """Verify intercluster storage using new helpers.""" + +def test_storage_variable_boundaries(): + """Verify storage with period-varying boundaries.""" +``` + +### 7.3 Plotting Tests + +```python +# Test new plot methods +def test_storage_by_cluster_facets(): + """Verify faceted cluster view.""" + +def test_cluster_heatmap(): + """Verify cluster heatmap rendering.""" +``` + +--- + +## Part 8: Decisions (Resolved) + +| Question | Decision | +|----------|----------| +| **Naming** | Use `cluster` as the dimension/coordinate name | +| **Indexer return type** | Always return proper multi-dimensional xarray DataArrays | +| **Segmentation** | tsam uniform segments only (sufficient for current needs) | +| **Backwards compatibility** | Not a concern - this is not released yet | + +--- + +## Appendix A: File Reference + +| File | Purpose | +|------|---------| +| `flixopt/clustering/base.py` | ClusterStructure, Clustering classes | +| `flixopt/clustering/intercluster_helpers.py` | SOC boundary utilities | +| `flixopt/flow_system.py` | FlowSystem with is_clustered property | +| `flixopt/transform_accessor.py` | cluster() method, solution expansion | +| `flixopt/components.py` | StorageModel, InterclusterStorageModel | +| `flixopt/features.py` | StatusModel, ShareAllocationModel | +| `flixopt/statistics_accessor.py` | Plotting methods | +| `flixopt/plotting.py` | Plot utilities | + +## Appendix B: Code Examples + +### B.1 Working with True (cluster, time) Dimensions + +```python +# ═══════════════════════════════════════════════════════════════ +# DIMENSION STRUCTURE +# ═══════════════════════════════════════════════════════════════ +# Non-clustered: +flow_rate.dims # ('time', 'period', 'scenario') +flow_rate.shape # (8760, ...) + +# Clustered: +flow_rate.dims # ('cluster', 'time', 'period', 'scenario') +flow_rate.shape # (9, 24, ...) # 9 clusters × 24 timesteps + +# ═══════════════════════════════════════════════════════════════ +# NATURAL CLUSTER BOUNDARY OPERATIONS +# ═══════════════════════════════════════════════════════════════ +# First/last timestep of each cluster - just use isel! +flow_at_start = flow_rate.isel(time=0) # Shape: (cluster, period, scenario) +flow_at_end = flow_rate.isel(time=-1) # Shape: (cluster, period, scenario) + +# Delta per cluster - trivial! +delta_per_cluster = flow_rate.isel(time=-1) - flow_rate.isel(time=0) + +# ═══════════════════════════════════════════════════════════════ +# TEMPORAL CONSTRAINTS - JUST WORK! +# ═══════════════════════════════════════════════════════════════ +# Storage balance - naturally stays within clusters +balance = charge_state[:, 1:] - charge_state[:, :-1] # No masking needed! + +# Status transitions - naturally per cluster +activate = status[:, 1:] - status[:, :-1] # No boundary issues! + +# ═══════════════════════════════════════════════════════════════ +# PER-CLUSTER AGGREGATION - use xarray directly +# ═══════════════════════════════════════════════════════════════ +mean_per_cluster = flow_rate.mean(dim='time') # Shape: (cluster, ...) +max_per_cluster = flow_rate.max(dim='time') +total_per_cluster = (flow_rate * timestep_duration).sum(dim='time') + +# ═══════════════════════════════════════════════════════════════ +# SELECT SPECIFIC WITHIN-CLUSTER TIMESTEP +# ═══════════════════════════════════════════════════════════════ +# Peak hour (hour 18) from each cluster +peak_values = flow_rate.isel(time=18) # Shape: (cluster, ...) + +# Multiple timesteps +morning_values = flow_rate.isel(time=slice(6, 12)) # Hours 6-11 from each cluster +``` + +### B.2 Storage Constraints with True Dimensions + +```python +# ═══════════════════════════════════════════════════════════════ +# charge_state has one extra timestep per cluster for boundaries +# ═══════════════════════════════════════════════════════════════ +# charge_state.dims = ('cluster', 'time_cs', 'period', 'scenario') +# charge_state.shape = (9, 25, ...) # 24 timesteps + 1 boundary + +# Charge balance - vectorized, no loops! +lhs = ( + charge_state[:, 1:] - # SOC at end of timestep + charge_state[:, :-1] * (1 - loss_rate) - # SOC at start, with loss + charge * eta_charge + # Charging adds energy + discharge / eta_discharge # Discharging removes energy +) +model.add_constraints(lhs == 0, name='charge_balance') + +# Delta SOC per cluster (for inter-cluster linking) +delta_soc = charge_state[:, -1] - charge_state[:, 0] # Shape: (cluster, ...) + +# Cluster start constraint (relative SOC starts at 0 within each cluster) +model.add_constraints(charge_state[:, 0] == 0, name='cluster_start') + +# Cyclic constraint (optional) +model.add_constraints( + charge_state[:, 0] == charge_state[:, -1], + name='cyclic' +) +``` + +### B.3 Cluster Plotting (Uses Existing API!) + +```python +# ═══════════════════════════════════════════════════════════════ +# FACET BY CLUSTER - uses existing facet_col parameter +# ═══════════════════════════════════════════════════════════════ +fs.statistics.plot.storage('Battery', facet_col='cluster') +fs.statistics.plot.balance('Heat', facet_col='cluster') +fs.statistics.plot.flows(..., facet_col='cluster') + +# ═══════════════════════════════════════════════════════════════ +# REGULAR PLOTS - auto-add cluster separators when clustered +# ═══════════════════════════════════════════════════════════════ +fs.statistics.plot.storage('Battery') # separators added automatically + +# ═══════════════════════════════════════════════════════════════ +# COMBINE WITH OTHER FACETS +# ═══════════════════════════════════════════════════════════════ +fs.statistics.plot.balance('Heat', facet_col='cluster', facet_row='scenario') +``` + +### B.4 Check Clustering Status and Access Properties + +```python +if flow_system.is_clustered: + clustering = flow_system.clustering + print(f"Clustered: {clustering.n_clusters} clusters × {clustering.timesteps_per_cluster} timesteps") + + # Dimension information + print(f"Data shape: (cluster={clustering.n_clusters}, time={clustering.timesteps_per_cluster})") + + # Cluster weights (how many original periods each cluster represents) + print(f"Cluster weights: {flow_system.cluster_weight.values}") + + # Aggregation weight (cluster_weight × timestep_duration) + print(f"Aggregation weight shape: {flow_system.aggregation_weight.shape}") +else: + print("Not clustered - full time resolution") +``` diff --git a/flixopt/features.py b/flixopt/features.py index bd43020a8..bb9864d64 100644 --- a/flixopt/features.py +++ b/flixopt/features.py @@ -214,13 +214,16 @@ def _do_modeling(self): self.add_variables(binary=True, short_name='startup', coords=self.get_coords()) self.add_variables(binary=True, short_name='shutdown', coords=self.get_coords()) + # Determine previous_state: None means relaxed (no constraint at t=0) + previous_state = self._previous_status.isel(time=-1) if self._previous_status is not None else None + BoundingPatterns.state_transition_bounds( self, state=self.status, activate=self.startup, deactivate=self.shutdown, name=f'{self.label_of_model}|switch', - previous_state=self._previous_status.isel(time=-1) if self._previous_status is not None else 0, + previous_state=previous_state, coord='time', ) @@ -231,12 +234,9 @@ def _do_modeling(self): coords=self._model.get_coords(('period', 'scenario')), short_name='startup_count', ) - # Apply cluster_weight to count startups correctly in clustered systems. - # A startup in a cluster with weight 10 represents 10 actual startups. # Sum over all temporal dimensions (time, and cluster if present) startup_temporal_dims = [d for d in self.startup.dims if d not in ('period', 'scenario')] - weighted_startup = self.startup * self._model.weights.get('cluster', 1.0) - self.add_constraints(count == weighted_startup.sum(startup_temporal_dims), short_name='startup_count') + self.add_constraints(count == self.startup.sum(startup_temporal_dims), short_name='startup_count') # 5. Consecutive active duration (uptime) using existing pattern if self.parameters.use_uptime_tracking: @@ -264,8 +264,19 @@ def _do_modeling(self): previous_duration=self._get_previous_downtime(), ) + # 7. Cyclic constraint for clustered systems + self._add_cluster_cyclic_constraint() + self._add_effects() + def _add_cluster_cyclic_constraint(self): + """For 'cyclic' cluster mode: each cluster's start status equals its end status.""" + if self._model.flow_system.clusters is not None and self.parameters.cluster_mode == 'cyclic': + self.add_constraints( + self.status.isel(time=0) == self.status.isel(time=-1), + short_name='cluster_cyclic', + ) + def _add_effects(self): """Add operational effects (use timestep_duration only, cluster_weight is applied when summing to total)""" if self.parameters.effects_per_active_hour: @@ -332,24 +343,22 @@ def downtime(self) -> linopy.Variable | None: def _get_previous_uptime(self): """Get previous uptime (consecutive active hours). - Returns 0 if no previous status is provided (assumes previously inactive). + Returns None if no previous status is provided (relaxed mode - no constraint at t=0). """ - hours_per_step = self._model.timestep_duration.isel(time=0).min().item() if self._previous_status is None: - return 0 - else: - return ModelingUtilities.compute_consecutive_hours_in_state(self._previous_status, hours_per_step) + return None # Relaxed mode + hours_per_step = self._model.timestep_duration.isel(time=0).min().item() + return ModelingUtilities.compute_consecutive_hours_in_state(self._previous_status, hours_per_step) def _get_previous_downtime(self): """Get previous downtime (consecutive inactive hours). - Returns one timestep duration if no previous status is provided (assumes previously inactive). + Returns None if no previous status is provided (relaxed mode - no constraint at t=0). """ - hours_per_step = self._model.timestep_duration.isel(time=0).min().item() if self._previous_status is None: - return hours_per_step - else: - return ModelingUtilities.compute_consecutive_hours_in_state(1 - self._previous_status, hours_per_step) + return None # Relaxed mode + hours_per_step = self._model.timestep_duration.isel(time=0).min().item() + return ModelingUtilities.compute_consecutive_hours_in_state(1 - self._previous_status, hours_per_step) class PieceModel(Submodel): diff --git a/flixopt/interface.py b/flixopt/interface.py index 13a9255da..227a63c7a 100644 --- a/flixopt/interface.py +++ b/flixopt/interface.py @@ -6,7 +6,7 @@ from __future__ import annotations import logging -from typing import TYPE_CHECKING, Any +from typing import TYPE_CHECKING, Any, Literal import numpy as np import pandas as pd @@ -1337,6 +1337,14 @@ class StatusParameters(Interface): force_startup_tracking: When True, creates startup variables even without explicit startup_limit constraint. Useful for tracking or reporting startup events without enforcing limits. + cluster_mode: How inter-timestep constraints are handled at cluster boundaries. + Only relevant when using ``transform.cluster()``. Options: + + - ``'relaxed'``: No constraint at cluster boundaries. Startups at the first + timestep of each cluster are not forced - the optimizer is free to choose. + This prevents clustering from inducing "phantom" startups. (default) + - ``'cyclic'``: Each cluster's final status equals its initial status. + Ensures consistent behavior within each representative period. Note: **Time Series Boundary Handling**: The final time period constraints for @@ -1472,6 +1480,7 @@ def __init__( max_downtime: Numeric_TPS | None = None, startup_limit: Numeric_PS | None = None, force_startup_tracking: bool = False, + cluster_mode: Literal['relaxed', 'cyclic'] = 'relaxed', ): self.effects_per_startup = effects_per_startup if effects_per_startup is not None else {} self.effects_per_active_hour = effects_per_active_hour if effects_per_active_hour is not None else {} @@ -1483,6 +1492,7 @@ def __init__( self.max_downtime = max_downtime self.startup_limit = startup_limit self.force_startup_tracking: bool = force_startup_tracking + self.cluster_mode = cluster_mode def transform_data(self) -> None: self.effects_per_startup = self._fit_effect_coords( diff --git a/flixopt/modeling.py b/flixopt/modeling.py index 6b81a0a4a..8f65349de 100644 --- a/flixopt/modeling.py +++ b/flixopt/modeling.py @@ -251,7 +251,7 @@ def consecutive_duration_tracking( maximum_duration: xr.DataArray | None = None, duration_dim: str = 'time', duration_per_step: int | float | xr.DataArray = None, - previous_duration: xr.DataArray = 0, + previous_duration: xr.DataArray | None = None, ) -> tuple[dict[str, linopy.Variable], dict[str, linopy.Constraint]]: """Creates consecutive duration tracking for a binary state variable. @@ -262,7 +262,7 @@ def consecutive_duration_tracking( duration[t] ≤ state[t] · M ∀t duration[t+1] ≤ duration[t] + duration_per_step[t] ∀t duration[t+1] ≥ duration[t] + duration_per_step[t] + (state[t+1] - 1) · M ∀t - duration[0] = (duration_per_step[0] + previous_duration) · state[0] + duration[0] = (duration_per_step[0] + previous_duration) · state[0] (if previous_duration is not None) If minimum_duration provided: duration[t] ≥ (state[t-1] - state[t]) · minimum_duration[t-1] ∀t > 0 @@ -278,16 +278,19 @@ def consecutive_duration_tracking( maximum_duration: Optional maximum consecutive duration (upper bound on duration variable) duration_dim: Dimension name to track duration along (default 'time') duration_per_step: Time increment per step in duration_dim - previous_duration: Initial duration value before first timestep (default 0) + previous_duration: Initial duration value before first timestep. If None (default), + the initial constraint is skipped ("relaxed" mode) - duration at t=0 is unconstrained. Returns: Tuple of (duration_variable, constraints_dict) - where constraints_dict contains: 'ub', 'forward', 'backward', 'initial', and optionally 'lb', 'initial_lb' + where constraints_dict contains: 'ub', 'forward', 'backward', and optionally 'initial', 'lb', 'initial_lb' """ if not isinstance(model, Submodel): raise ValueError('ModelingPrimitives.consecutive_duration_tracking() can only be used with a Submodel') - mega = duration_per_step.sum(duration_dim) + previous_duration # Big-M value + # Big-M value (use 0 for previous_duration if None/relaxed mode) + prev_for_mega = previous_duration if previous_duration is not None else 0 + mega = duration_per_step.sum(duration_dim) + prev_for_mega # Duration variable duration = model.add_variables( @@ -319,12 +322,14 @@ def consecutive_duration_tracking( name=f'{duration.name}|backward', ) - # Initial condition: duration[0] = (duration_per_step[0] + previous_duration) * state[0] - constraints['initial'] = model.add_constraints( - duration.isel({duration_dim: 0}) - == (duration_per_step.isel({duration_dim: 0}) + previous_duration) * state.isel({duration_dim: 0}), - name=f'{duration.name}|initial', - ) + # Initial condition (skip if previous_duration is None = relaxed mode) + if previous_duration is not None: + # duration[0] = (duration_per_step[0] + previous_duration) * state[0] + constraints['initial'] = model.add_constraints( + duration.isel({duration_dim: 0}) + == (duration_per_step.isel({duration_dim: 0}) + previous_duration) * state.isel({duration_dim: 0}), + name=f'{duration.name}|initial', + ) # Minimum duration constraint if provided if minimum_duration is not None: @@ -335,17 +340,18 @@ def consecutive_duration_tracking( name=f'{duration.name}|lb', ) - # Handle initial condition for minimum duration - prev = ( - float(previous_duration) - if not isinstance(previous_duration, xr.DataArray) - else float(previous_duration.max().item()) - ) - min0 = float(minimum_duration.isel({duration_dim: 0}).max().item()) - if prev > 0 and prev < min0: - constraints['initial_lb'] = model.add_constraints( - state.isel({duration_dim: 0}) == 1, name=f'{duration.name}|initial_lb' + # Handle initial condition for minimum duration (only if not relaxed mode) + if previous_duration is not None: + prev = ( + float(previous_duration) + if not isinstance(previous_duration, xr.DataArray) + else float(previous_duration.max().item()) ) + min0 = float(minimum_duration.isel({duration_dim: 0}).max().item()) + if prev > 0 and prev < min0: + constraints['initial_lb'] = model.add_constraints( + state.isel({duration_dim: 0}) == 1, name=f'{duration.name}|initial_lb' + ) variables = {'duration': duration} @@ -578,9 +584,9 @@ def state_transition_bounds( activate: linopy.Variable, deactivate: linopy.Variable, name: str, - previous_state: float | xr.DataArray = 0, + previous_state: float | xr.DataArray | None = None, coord: str = 'time', - ) -> tuple[linopy.Constraint, linopy.Constraint, linopy.Constraint]: + ) -> tuple[linopy.Constraint, linopy.Constraint | None, linopy.Constraint]: """Creates state transition constraints for binary state variables. Tracks transitions between active (1) and inactive (0) states using @@ -588,7 +594,7 @@ def state_transition_bounds( Mathematical formulation: activate[t] - deactivate[t] = state[t] - state[t-1] ∀t > 0 - activate[0] - deactivate[0] = state[0] - previous_state + activate[0] - deactivate[0] = state[0] - previous_state (if previous_state is not None) activate[t] + deactivate[t] ≤ 1 ∀t activate[t], deactivate[t] ∈ {0, 1} @@ -598,11 +604,13 @@ def state_transition_bounds( activate: Binary variable for transitions from inactive to active (0→1) deactivate: Binary variable for transitions from active to inactive (1→0) name: Base name for constraints - previous_state: State value before first timestep (default 0) + previous_state: State value before first timestep. If None (default), the initial + constraint is skipped ("relaxed" mode) - startup/shutdown at t=0 is not forced. coord: Time dimension name (default 'time') Returns: - Tuple of (transition_constraint, initial_constraint, mutex_constraint) + Tuple of (transition_constraint, initial_constraint, mutex_constraint). + initial_constraint is None if previous_state is None (relaxed mode). """ if not isinstance(model, Submodel): raise ValueError('BoundingPatterns.state_transition_bounds() can only be used with a Submodel') @@ -614,11 +622,14 @@ def state_transition_bounds( name=f'{name}|transition', ) - # Initial state transition for t = 0 - initial = model.add_constraints( - activate.isel({coord: 0}) - deactivate.isel({coord: 0}) == state.isel({coord: 0}) - previous_state, - name=f'{name}|initial', - ) + # Initial state transition for t = 0 (skip if previous_state is None = relaxed mode) + if previous_state is not None: + initial = model.add_constraints( + activate.isel({coord: 0}) - deactivate.isel({coord: 0}) == state.isel({coord: 0}) - previous_state, + name=f'{name}|initial', + ) + else: + initial = None # Relaxed mode: no constraint at t=0, startup/shutdown is "free" # At most one transition per timestep (mutual exclusivity) mutex = model.add_constraints(activate + deactivate <= 1, name=f'{name}|mutex') diff --git a/tests/deprecated/test_on_hours_computation.py b/tests/deprecated/test_on_hours_computation.py index 578fd7792..c74332565 100644 --- a/tests/deprecated/test_on_hours_computation.py +++ b/tests/deprecated/test_on_hours_computation.py @@ -9,7 +9,7 @@ class TestComputeConsecutiveDuration: """Tests for the compute_consecutive_hours_in_state static method.""" @pytest.mark.parametrize( - 'binary_values, hours_per_timestep, expected', + 'binary_values, timestep_duration, expected', [ # Case 1: Single timestep DataArrays (xr.DataArray([1], dims=['time']), 5, 5), @@ -26,22 +26,22 @@ class TestComputeConsecutiveDuration: (xr.DataArray([0, 1, 1, 1, 0, 0], dims=['time']), 1, 0), # ends with 0 ], ) - def test_compute_duration(self, binary_values, hours_per_timestep, expected): + def test_compute_duration(self, binary_values, timestep_duration, expected): """Test compute_consecutive_hours_in_state with various inputs.""" - result = ModelingUtilities.compute_consecutive_hours_in_state(binary_values, hours_per_timestep) + result = ModelingUtilities.compute_consecutive_hours_in_state(binary_values, timestep_duration) assert np.isclose(result, expected) @pytest.mark.parametrize( - 'binary_values, hours_per_timestep', + 'binary_values, timestep_duration', [ - # Case: hours_per_timestep must be scalar + # Case: timestep_duration must be scalar (xr.DataArray([1, 1, 1, 1, 1], dims=['time']), np.array([1, 2])), ], ) - def test_compute_duration_raises_error(self, binary_values, hours_per_timestep): + def test_compute_duration_raises_error(self, binary_values, timestep_duration): """Test error conditions.""" with pytest.raises(TypeError): - ModelingUtilities.compute_consecutive_hours_in_state(binary_values, hours_per_timestep) + ModelingUtilities.compute_consecutive_hours_in_state(binary_values, timestep_duration) class TestComputePreviousOnStates: diff --git a/tests/test_on_hours_computation.py b/tests/test_on_hours_computation.py index 578fd7792..c74332565 100644 --- a/tests/test_on_hours_computation.py +++ b/tests/test_on_hours_computation.py @@ -9,7 +9,7 @@ class TestComputeConsecutiveDuration: """Tests for the compute_consecutive_hours_in_state static method.""" @pytest.mark.parametrize( - 'binary_values, hours_per_timestep, expected', + 'binary_values, timestep_duration, expected', [ # Case 1: Single timestep DataArrays (xr.DataArray([1], dims=['time']), 5, 5), @@ -26,22 +26,22 @@ class TestComputeConsecutiveDuration: (xr.DataArray([0, 1, 1, 1, 0, 0], dims=['time']), 1, 0), # ends with 0 ], ) - def test_compute_duration(self, binary_values, hours_per_timestep, expected): + def test_compute_duration(self, binary_values, timestep_duration, expected): """Test compute_consecutive_hours_in_state with various inputs.""" - result = ModelingUtilities.compute_consecutive_hours_in_state(binary_values, hours_per_timestep) + result = ModelingUtilities.compute_consecutive_hours_in_state(binary_values, timestep_duration) assert np.isclose(result, expected) @pytest.mark.parametrize( - 'binary_values, hours_per_timestep', + 'binary_values, timestep_duration', [ - # Case: hours_per_timestep must be scalar + # Case: timestep_duration must be scalar (xr.DataArray([1, 1, 1, 1, 1], dims=['time']), np.array([1, 2])), ], ) - def test_compute_duration_raises_error(self, binary_values, hours_per_timestep): + def test_compute_duration_raises_error(self, binary_values, timestep_duration): """Test error conditions.""" with pytest.raises(TypeError): - ModelingUtilities.compute_consecutive_hours_in_state(binary_values, hours_per_timestep) + ModelingUtilities.compute_consecutive_hours_in_state(binary_values, timestep_duration) class TestComputePreviousOnStates: