Multi-Period Clustering with cluster()¶

import timeit

import numpy as np
import pandas as pd
import plotly.express as px

import flixopt as fx

fx.CONFIG.notebook()

from data.generate_example_systems import create_multiperiod_system

flow_system = create_multiperiod_system()

print(f'Timesteps: {len(flow_system.timesteps)} ({len(flow_system.timesteps) // 24} days)')
print(f'Periods: {list(flow_system.periods.values)}')
print(f'Scenarios: {list(flow_system.scenarios.values)}')
print(f'Scenario weights: {flow_system.scenario_weights.values}')
print(f'\nComponents: {list(flow_system.components.keys())}')

# Select first week only (168 hours)
flow_system = flow_system.transform.isel(time=slice(0, 168))

print(f'After isel: {len(flow_system.timesteps)} timesteps ({len(flow_system.timesteps) // 24} days)')

# Visualize demand scenarios (equal across periods)
heat_demand = flow_system.components['Building'].inputs[0].fixed_relative_profile

fig = px.line(heat_demand.to_dataframe('value').reset_index(), x='time', y='value', facet_row='scenario')

fig.update_layout(
    height=350,
    title='Heat Demand by Scenario (One Week)',
    xaxis_title='Time',
    yaxis_title='Heat Demand [kW]',
)
fig.show()

solver = fx.solvers.HighsSolver(mip_gap=0.01)

start = timeit.default_timer()
fs_full = flow_system.copy()
fs_full.name = 'Full Optimization'
fs_full.optimize(solver)
time_full = timeit.default_timer() - start

print(f'Full optimization: {time_full:.2f} seconds')
print(f'Total cost (objective): {fs_full.solution["objective"].item():,.0f} €')
print('\nOptimized sizes:')
for name, size in fs_full.stats.sizes.items():
    print(f'  {name}: {size.max().item():.1f}')

from tsam import ExtremeConfig

start = timeit.default_timer()

# Force inclusion of peak demand periods
peak_series = ['Building(Heat)|fixed_relative_profile']

# Cluster to 3 typical days (from 7 days)
fs_clustered = flow_system.transform.cluster(
    n_clusters=3,
    cluster_duration='1D',
    extremes=ExtremeConfig(method='replace', max_value=peak_series, preserve_n_clusters=True),
)

time_clustering = timeit.default_timer() - start

print(f'Clustering time: {time_clustering:.2f} seconds')
print(f'Reduced: {len(flow_system.timesteps)} → {len(fs_clustered.timesteps)} timesteps per period')
print('Total problem reduction: 7 days × 3 periods × 2 scenarios → 3 days × 3 × 2')

# Optimize the reduced system
start = timeit.default_timer()
fs_clustered.optimize(solver)
time_clustered = timeit.default_timer() - start

print(f'Clustered optimization: {time_clustered:.2f} seconds')
print(f'Total cost (objective): {fs_clustered.solution["objective"].item():,.0f} €')
print(f'\nSpeedup vs full: {time_full / (time_clustering + time_clustered):.1f}x')
print('\nOptimized sizes:')
for name, size in fs_clustered.stats.sizes.items():
    print(f'  {name}: {size.max().item():.1f}')

# Duration curves show how well the distribution is preserved per period/scenario
fs_clustered.clustering.plot.compare(
    kind='duration_curve',
)

# Heatmap shows cluster assignments - faceted by period and scenario
fs_clustered.clustering.plot.heatmap()

clustering = fs_clustered.clustering

print('Clustering Configuration:')
print(f'  Typical periods (clusters): {clustering.n_clusters}')
print(f'  Timesteps per cluster: {clustering.timesteps_per_cluster}')

# Access underlying results via xarray-like interface
print(f'\nClusteringResults dimensions: {clustering.results.dims}')
print(f'ClusteringResults coords: {clustering.results.coords}')

# The cluster_assignments shows which cluster each original day belongs to
# For multi-period systems, select a specific period/scenario combination
cluster_assignments = clustering.cluster_assignments.isel(period=0, scenario=0).values
day_names = ['Mon', 'Tue', 'Wed', 'Thu', 'Fri', 'Sat', 'Sun']

print('\nCluster assignments per day (period=2024, scenario=High):')
for i, cluster_id in enumerate(cluster_assignments):
    print(f'  {day_names[i]}: Cluster {cluster_id}')

# Cluster occurrences (how many original days each cluster represents)
unique, counts = np.unique(cluster_assignments, return_counts=True)
print('\nCluster weights (days represented):')
for cluster_id, count in zip(unique, counts, strict=True):
    print(f'  Cluster {cluster_id}: {count} days')

# Stage 1 already done - apply safety margin
SAFETY_MARGIN = 1.10  # 10% buffer for multi-period uncertainty

sizes_with_margin = {name: size.max().item() * SAFETY_MARGIN for name, size in fs_clustered.stats.sizes.items()}

print('Stage 1: Sizing with clustering')
print(f'  Time: {time_clustering + time_clustered:.2f} seconds')
print(f'  Cost estimate: {fs_clustered.solution["objective"].item():,.0f} €')
print(f'\nSizes with {(SAFETY_MARGIN - 1) * 100:.0f}% safety margin:')
for name, size in sizes_with_margin.items():
    original = fs_clustered.stats.sizes[name].max().item()
    print(f'  {name}: {original:.1f} → {size:.1f}')

# Stage 2: Full resolution dispatch with fixed sizes
print('Stage 2: Full resolution dispatch')
start = timeit.default_timer()

fs_dispatch = flow_system.transform.fix_sizes(sizes_with_margin)
fs_dispatch.name = 'Two-Stage'
fs_dispatch.optimize(solver)

time_dispatch = timeit.default_timer() - start

print(f'  Time: {time_dispatch:.2f} seconds')
print(f'  Actual cost: {fs_dispatch.solution["objective"].item():,.0f} €')

# Total comparison
total_two_stage = time_clustering + time_clustered + time_dispatch
print(f'\nTotal two-stage time: {total_two_stage:.2f} seconds')
print(f'Speedup vs full: {time_full / total_two_stage:.1f}x')

results = {
    'Full (baseline)': {
        'Time [s]': time_full,
        'Cost [€]': fs_full.solution['objective'].item(),
        'Boiler': fs_full.stats.sizes['Boiler(Heat)'].max().item(),
        'Storage': fs_full.stats.sizes['ThermalStorage'].max().item(),
    },
    'Clustered (3 days)': {
        'Time [s]': time_clustering + time_clustered,
        'Cost [€]': fs_clustered.solution['objective'].item(),
        'Boiler': fs_clustered.stats.sizes['Boiler(Heat)'].max().item(),
        'Storage': fs_clustered.stats.sizes['ThermalStorage'].max().item(),
    },
    'Two-Stage': {
        'Time [s]': total_two_stage,
        'Cost [€]': fs_dispatch.solution['objective'].item(),
        'Boiler': sizes_with_margin['Boiler(Heat)'],
        'Storage': sizes_with_margin['ThermalStorage'],
    },
}

comparison = pd.DataFrame(results).T
baseline_cost = comparison.loc['Full (baseline)', 'Cost [€]']
baseline_time = comparison.loc['Full (baseline)', 'Time [s]']
comparison['Cost Gap [%]'] = ((comparison['Cost [€]'] - baseline_cost) / abs(baseline_cost) * 100).round(2)
comparison['Speedup'] = (baseline_time / comparison['Time [s]']).round(1)

comparison.style.format(
    {
        'Time [s]': '{:.2f}',
        'Cost [€]': '{:,.0f}',
        'Boiler': '{:.1f}',
        'Storage': '{:.0f}',
        'Cost Gap [%]': '{:.2f}',
        'Speedup': '{:.1f}x',
    }
)

# Plot flow rates with automatic faceting by period and scenario
fs_full.stats.plot.flows(component='Boiler')

# Side-by-side comparison using the Comparison class
comp = fx.Comparison([fs_full, fs_dispatch])
comp.stats.plot.balance('Heat')

# Expand the clustered solution
fs_expanded = fs_clustered.transform.expand()

print(f'Expanded: {len(fs_clustered.timesteps)} → {len(fs_expanded.timesteps)} timesteps')
print(f'Cost (objective): {fs_expanded.solution["objective"].item():,.0f} €')

# Compare expanded solution - shows the repeated cluster patterns
fs_expanded.stats.plot.flows(component='Boiler')