Clustering Storage Modes¶

import timeit

import pandas as pd
import plotly.graph_objects as go
from plotly.subplots import make_subplots

import flixopt as fx

fx.CONFIG.notebook()

from data.generate_example_systems import create_seasonal_storage_system

flow_system = create_seasonal_storage_system()
flow_system.connect_and_transform()  # Align all data as xarray

timesteps = flow_system.timesteps
print(f'FlowSystem: {len(timesteps)} timesteps ({len(timesteps) / 24:.0f} days)')
print(f'Components: {list(flow_system.components.keys())}')

# Visualize the seasonal patterns
solar_profile = flow_system.components['SolarThermal'].outputs[0].fixed_relative_profile
heat_demand = flow_system.components['HeatDemand'].inputs[0].fixed_relative_profile

# Compute daily averages using xarray resample
solar_daily = solar_profile.resample(time='1D').mean()
demand_daily = heat_demand.resample(time='1D').mean()

fig = make_subplots(rows=2, cols=1, shared_xaxes=True, vertical_spacing=0.1)
fig.add_trace(
    go.Scatter(x=solar_daily.time.values, y=solar_daily.values, name='Solar (daily avg)', fill='tozeroy'), row=1, col=1
)
fig.add_trace(
    go.Scatter(x=demand_daily.time.values, y=demand_daily.values, name='Heat Demand (daily avg)', fill='tozeroy'),
    row=2,
    col=1,
)
fig.update_layout(height=400, title='Seasonal Mismatch: Solar vs Heat Demand')
fig.update_xaxes(title_text='Day of Year', row=2, col=1)
fig.update_yaxes(title_text='Solar Profile', row=1, col=1)
fig.update_yaxes(title_text='Heat Demand [MW]', row=2, col=1)
fig.show()

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

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:.1f} seconds')
print(f'Total cost: {fs_full.solution["costs"].item():,.0f} EUR')
print('\nOptimized sizes:')
for name, size in fs_full.stats.sizes.items():
    print(f'  {name}: {float(size.item()):.2f}')

from tsam import ExtremeConfig

# Clustering parameters
N_CLUSTERS = 24  # 24 typical days for a full year
CLUSTER_DURATION = '1D'
PEAK_SERIES = ['HeatDemand(Q_th)|fixed_relative_profile']

# Storage modes to compare
storage_modes = ['independent', 'cyclic', 'intercluster', 'intercluster_cyclic']

results = {}
clustered_systems = {}

for mode in storage_modes:
    print(f'\n--- Mode: {mode} ---')

    # Create a copy and set the storage mode
    fs_copy = flow_system.copy()
    fs_copy.storages['SeasonalStorage'].cluster_mode = mode

    start = timeit.default_timer()
    fs_clustered = fs_copy.transform.cluster(
        n_clusters=N_CLUSTERS,
        cluster_duration=CLUSTER_DURATION,
        extremes=ExtremeConfig(method='new_cluster', max_value=PEAK_SERIES),
    )
    time_cluster = timeit.default_timer() - start

    start = timeit.default_timer()
    fs_clustered.optimize(solver)
    time_solve = timeit.default_timer() - start

    clustered_systems[mode] = fs_clustered

    results[mode] = {
        'Time [s]': time_cluster + time_solve,
        'Cost [EUR]': fs_clustered.solution['costs'].item(),
        'Solar [MW]': fs_clustered.stats.sizes.get('SolarThermal(Q_th)', 0),
        'Boiler [MW]': fs_clustered.stats.sizes.get('GasBoiler(Q_th)', 0),
        'Storage [MWh]': fs_clustered.stats.sizes.get('SeasonalStorage', 0),
    }

    # Handle xarray types
    for key in ['Solar [MW]', 'Boiler [MW]', 'Storage [MWh]']:
        val = results[mode][key]
        results[mode][key] = float(val.item()) if hasattr(val, 'item') else float(val)

    print(f'  Time: {results[mode]["Time [s]"]:.1f}s')
    print(f'  Cost: {results[mode]["Cost [EUR]"]:,.0f} EUR')
    print(f'  Storage: {results[mode]["Storage [MWh]"]:.0f} MWh')

# Add full optimization result for comparison
results['Full (baseline)'] = {
    'Time [s]': time_full,
    'Cost [EUR]': fs_full.solution['costs'].item(),
    'Solar [MW]': float(fs_full.stats.sizes.get('SolarThermal(Q_th)', 0).item()),
    'Boiler [MW]': float(fs_full.stats.sizes.get('GasBoiler(Q_th)', 0).item()),
    'Storage [MWh]': float(fs_full.stats.sizes.get('SeasonalStorage', 0).item()),
}

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

comparison.style.format(
    {
        'Time [s]': '{:.1f}',
        'Cost [EUR]': '{:,.0f}',
        'Solar [MW]': '{:.1f}',
        'Boiler [MW]': '{:.1f}',
        'Storage [MWh]': '{:.0f}',
        'Cost Gap [%]': '{:+.1f}',
        'Speedup': '{:.1f}x',
    }
)

# Plot storage charge state for each mode
fig = make_subplots(
    rows=len(storage_modes) + 1,
    cols=1,
    shared_xaxes=True,
    vertical_spacing=0.05,
    subplot_titles=['Full Optimization'] + [f'Mode: {m}' for m in storage_modes],
)

# Full optimization
soc_full = fs_full.solution['SeasonalStorage|charge_state']
fig.add_trace(go.Scatter(x=fs_full.timesteps, y=soc_full.values, name='Full', line=dict(width=0.8)), row=1, col=1)

# Expanded clustered solutions
for i, mode in enumerate(storage_modes, start=2):
    fs_exp = expanded_systems[mode]
    soc = fs_exp.solution['SeasonalStorage|charge_state']
    fig.add_trace(go.Scatter(x=fs_exp.timesteps, y=soc.values, name=mode, line=dict(width=0.8)), row=i, col=1)

fig.update_layout(height=800, title='Storage Charge State by Mode', showlegend=False)
for i in range(1, len(storage_modes) + 2):
    fig.update_yaxes(title_text='SOC [MWh]', row=i, col=1)
fig.show()

# Compare full optimization with the recommended intercluster_cyclic mode
comp = fx.Comparison([fs_full, expanded_systems['intercluster_cyclic']])
comp.stats.plot.balance('Heat')