Aggregation¶

import timeit

import pandas as pd
import plotly.express as px
import xarray as xr

import flixopt as fx

fx.CONFIG.notebook()

# Load time series data (15-min resolution)
data = pd.read_csv('data/Zeitreihen2020.csv', index_col=0, parse_dates=True).sort_index()
data = data['2020-01-01':'2020-01-31 23:45:00']  # One month
data.index.name = 'time'  # Rename index for consistency

timesteps = data.index

# Extract profiles
electricity_demand = data['P_Netz/MW'].to_numpy()
heat_demand = data['Q_Netz/MW'].to_numpy()
electricity_price = data['Strompr.€/MWh'].to_numpy()
gas_price = data['Gaspr.€/MWh'].to_numpy()

print(f'Timesteps: {len(timesteps)} ({len(timesteps) / 96:.0f} days at 15-min resolution)')
print(f'Heat demand: {heat_demand.min():.1f} - {heat_demand.max():.1f} MW')
print(f'Electricity price: {electricity_price.min():.1f} - {electricity_price.max():.1f} €/MWh')

# Visualize first week
profiles = xr.Dataset(
    {
        'Heat Demand [MW]': xr.DataArray(heat_demand[:672], dims=['time'], coords={'time': timesteps[:672]}),
        'Electricity Price [€/MWh]': xr.DataArray(
            electricity_price[:672], dims=['time'], coords={'time': timesteps[:672]}
        ),
    }
)

df = profiles.to_dataframe().reset_index().melt(id_vars='time', var_name='variable', value_name='value')
fig = px.line(df, x='time', y='value', facet_col='variable', height=300)
fig.update_yaxes(matches=None, showticklabels=True)
fig.for_each_annotation(lambda a: a.update(text=a.text.split('=')[-1]))
fig

def build_system(timesteps, heat_demand, electricity_demand, electricity_price, gas_price):
    """Build a district heating system with CHP, boiler, and storage (with investment options)."""
    fs = fx.FlowSystem(timesteps)

    fs.add_elements(
        # Buses
        fx.Bus('Electricity'),
        fx.Bus('Heat'),
        fx.Bus('Gas'),
        fx.Bus('Coal'),
        # Effects
        fx.Effect('costs', '€', 'Total Costs', is_standard=True, is_objective=True),
        fx.Effect('CO2', 'kg', 'CO2 Emissions'),
        # CHP with investment optimization
        fx.linear_converters.CHP(
            'CHP',
            thermal_efficiency=0.58,
            electrical_efficiency=0.22,
            electrical_flow=fx.Flow('P_el', bus='Electricity', size=200),
            thermal_flow=fx.Flow(
                'Q_th',
                bus='Heat',
                size=fx.InvestParameters(
                    minimum_size=100,
                    maximum_size=300,
                    effects_of_investment_per_size={'costs': 10},
                ),
            ),
            fuel_flow=fx.Flow('Q_fu', bus='Coal'),
        ),
        # Gas Boiler with investment optimization
        fx.linear_converters.Boiler(
            'Boiler',
            thermal_efficiency=0.85,
            thermal_flow=fx.Flow(
                'Q_th',
                bus='Heat',
                size=fx.InvestParameters(
                    minimum_size=0,
                    maximum_size=150,
                    effects_of_investment_per_size={'costs': 5},
                ),
            ),
            fuel_flow=fx.Flow('Q_fu', bus='Gas'),
        ),
        # Thermal Storage with investment optimization
        fx.Storage(
            'Storage',
            capacity_in_flow_hours=fx.InvestParameters(
                minimum_size=0,
                maximum_size=1000,
                effects_of_investment_per_size={'costs': 0.5},
            ),
            initial_charge_state=0,
            eta_charge=1,
            eta_discharge=1,
            relative_loss_per_hour=0.001,
            charging=fx.Flow('Charge', size=137, bus='Heat'),
            discharging=fx.Flow('Discharge', size=158, bus='Heat'),
        ),
        # Fuel sources
        fx.Source(
            'GasGrid',
            outputs=[fx.Flow('Q_Gas', bus='Gas', size=1000, effects_per_flow_hour={'costs': gas_price, 'CO2': 0.3})],
        ),
        fx.Source(
            'CoalSupply',
            outputs=[fx.Flow('Q_Coal', bus='Coal', size=1000, effects_per_flow_hour={'costs': 4.6, 'CO2': 0.3})],
        ),
        # Electricity grid connection
        fx.Source(
            'GridBuy',
            outputs=[
                fx.Flow(
                    'P_el',
                    bus='Electricity',
                    size=1000,
                    effects_per_flow_hour={'costs': electricity_price + 0.5, 'CO2': 0.3},
                )
            ],
        ),
        fx.Sink(
            'GridSell',
            inputs=[fx.Flow('P_el', bus='Electricity', size=1000, effects_per_flow_hour=-(electricity_price - 0.5))],
        ),
        # Demands
        fx.Sink('HeatDemand', inputs=[fx.Flow('Q_th', bus='Heat', size=1, fixed_relative_profile=heat_demand)]),
        fx.Sink(
            'ElecDemand', inputs=[fx.Flow('P_el', bus='Electricity', size=1, fixed_relative_profile=electricity_demand)]
        ),
    )

    return fs


flow_system = build_system(timesteps, heat_demand, electricity_demand, electricity_price, gas_price)
print(f'System: {len(timesteps)} timesteps')

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

# Resample from 15min to 4h resolution
fs_resampled = flow_system.transform.resample('4h')

reduction = (1 - len(fs_resampled.timesteps) / len(flow_system.timesteps)) * 100
print(f'Resampled: {len(flow_system.timesteps)} → {len(fs_resampled.timesteps)} timesteps ({reduction:.0f}% reduction)')

# Optimize resampled system
start = timeit.default_timer()
fs_resampled.optimize(solver)
time_resampled = timeit.default_timer() - start

print(f'Resampled: {time_resampled:.1f}s, {fs_resampled.solution["costs"].item():,.0f} €')

# Stage 1: Sizing with resampled data
start = timeit.default_timer()
fs_sizing = flow_system.transform.resample('4h')
fs_sizing.optimize(solver)
time_stage1 = timeit.default_timer() - start

sizes = {k: float(v.item()) for k, v in fs_sizing.statistics.sizes.items()}
print(
    f'Stage 1 (sizing): {time_stage1:.1f}s → CHP {sizes["CHP(Q_th)"]:.0f}, Boiler {sizes["Boiler(Q_th)"]:.0f}, Storage {sizes["Storage"]:.0f}'
)

# Stage 2: Dispatch at full resolution with fixed sizes
start = timeit.default_timer()
fs_dispatch = flow_system.transform.fix_sizes(fs_sizing.statistics.sizes)
fs_dispatch.optimize(solver)
time_stage2 = timeit.default_timer() - start

print(
    f'Stage 2 (dispatch): {time_stage2:.1f}s, {fs_dispatch.solution["costs"].item():,.0f} € (total: {time_stage1 + time_stage2:.1f}s)'
)

start = timeit.default_timer()
fs_full = flow_system.copy()
fs_full.optimize(solver)
time_full = timeit.default_timer() - start

print(f'Full optimization: {time_full:.1f}s, {fs_full.solution["costs"].item():,.0f} €')

# Collect results
results = {
    'Full (baseline)': {
        'Time [s]': time_full,
        'Cost [€]': fs_full.solution['costs'].item(),
        'CHP Size [MW]': fs_full.statistics.sizes['CHP(Q_th)'].item(),
        'Boiler Size [MW]': fs_full.statistics.sizes['Boiler(Q_th)'].item(),
        'Storage Size [MWh]': fs_full.statistics.sizes['Storage'].item(),
    },
    'Resampled (4h)': {
        'Time [s]': time_resampled,
        'Cost [€]': fs_resampled.solution['costs'].item(),
        'CHP Size [MW]': fs_resampled.statistics.sizes['CHP(Q_th)'].item(),
        'Boiler Size [MW]': fs_resampled.statistics.sizes['Boiler(Q_th)'].item(),
        'Storage Size [MWh]': fs_resampled.statistics.sizes['Storage'].item(),
    },
    'Two-Stage': {
        'Time [s]': time_stage1 + time_stage2,
        'Cost [€]': fs_dispatch.solution['costs'].item(),
        'CHP Size [MW]': fs_dispatch.statistics.sizes['CHP(Q_th)'].item(),
        'Boiler Size [MW]': fs_dispatch.statistics.sizes['Boiler(Q_th)'].item(),
        'Storage Size [MWh]': fs_dispatch.statistics.sizes['Storage'].item(),
    },
}

comparison = pd.DataFrame(results).T

# Add relative metrics
baseline_cost = comparison.loc['Full (baseline)', 'Cost [€]']
baseline_time = comparison.loc['Full (baseline)', 'Time [s]']
comparison['Cost Gap [%]'] = ((comparison['Cost [€]'] - baseline_cost) / baseline_cost * 100).round(2)
comparison['Speedup'] = (baseline_time / comparison['Time [s]']).round(1)

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

# Full optimization heat balance
fs_full.statistics.plot.balance('Heat')

# Two-stage optimization heat balance
fs_dispatch.statistics.plot.balance('Heat')

fs_full.statistics.plot.sankey.flows()