Rolling Horizon¶

import timeit

import pandas as pd
import plotly.express as px
import plotly.graph_objects as go
import xarray as xr
from plotly.subplots import make_subplots

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-14 23:45:00']  # Two weeks
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'{len(timesteps)} timesteps ({len(timesteps) / 96:.0f} days), heat {heat_demand.min():.0f}-{heat_demand.max():.0f} MW'
)

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

    # Effects

    # Buses
    fs.add_elements(
        fx.Bus('Electricity'),
        fx.Bus('Heat'),
        fx.Bus('Gas'),
        fx.Bus('Coal'),
        fx.Effect('costs', '€', 'Total Costs', is_standard=True, is_objective=True),
        fx.Effect('CO2', 'kg', 'CO2 Emissions'),
        fx.linear_converters.CHP(
            'CHP',
            thermal_efficiency=0.58,
            electrical_efficiency=0.22,
            status_parameters=fx.StatusParameters(effects_per_startup=24000),
            electrical_flow=fx.Flow('P_el', bus='Electricity', size=200),
            thermal_flow=fx.Flow('Q_th', bus='Heat', size=200),
            fuel_flow=fx.Flow('Q_fu', bus='Coal', size=288, relative_minimum=87 / 288, previous_flow_rate=100),
        ),
        fx.linear_converters.Boiler(
            'Boiler',
            thermal_efficiency=0.85,
            thermal_flow=fx.Flow('Q_th', bus='Heat'),
            fuel_flow=fx.Flow(
                'Q_fu',
                bus='Gas',
                size=95,
                relative_minimum=12 / 95,
                previous_flow_rate=20,
                status_parameters=fx.StatusParameters(effects_per_startup=1000),
            ),
        ),
        fx.Storage(
            'Storage',
            capacity_in_flow_hours=684,
            initial_charge_state=137,
            minimal_final_charge_state=137,
            maximal_final_charge_state=158,
            eta_charge=1,
            eta_discharge=1,
            relative_loss_per_hour=0.001,
            prevent_simultaneous_charge_and_discharge=True,
            charging=fx.Flow('Charge', size=137, bus='Heat'),
            discharging=fx.Flow('Discharge', size=158, bus='Heat'),
        ),
        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})],
        ),
        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))],
        ),
        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()

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

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

start = timeit.default_timer()
fs_rolling = flow_system.copy()
segments = fs_rolling.optimize.rolling_horizon(
    solver,
    horizon=192,  # 2-day segments (192 timesteps at 15-min resolution)
    overlap=96,  # 1-day lookahead
)
time_rolling = timeit.default_timer() - start

print(f'Rolling ({len(segments)} segments): {time_rolling:.1f}s, {fs_rolling.solution["costs"].item():,.0f} €')

cost_full = fs_full.solution['costs'].item()
cost_rolling = fs_rolling.solution['costs'].item()
cost_gap = (cost_rolling - cost_full) / cost_full * 100

results = pd.DataFrame(
    {
        'Method': ['Full optimization', 'Rolling horizon'],
        'Time [s]': [time_full, time_rolling],
        'Cost [€]': [cost_full, cost_rolling],
        'Cost Gap [%]': [0.0, cost_gap],
    }
).set_index('Method')

results.style.format({'Time [s]': '{:.2f}', 'Cost [€]': '{:,.0f}', 'Cost Gap [%]': '{:.2f}'})

fs_full.statistics.plot.balance('Heat').figure.update_layout(title='Heat Balance (Full)')

fs_rolling.statistics.plot.balance('Heat').figure.update_layout(title='Heat Balance (Rolling)')

fig = make_subplots(
    rows=2, cols=1, shared_xaxes=True, vertical_spacing=0.1, subplot_titles=['Full Optimization', 'Rolling Horizon']
)

# Full optimization
charge_full = fs_full.solution['Storage|charge_state'].values[:-1]  # Drop final value
fig.add_trace(go.Scatter(x=timesteps, y=charge_full, name='Full', line=dict(color='blue')), row=1, col=1)

# Rolling horizon
charge_rolling = fs_rolling.solution['Storage|charge_state'].values[:-1]
fig.add_trace(go.Scatter(x=timesteps, y=charge_rolling, name='Rolling', line=dict(color='orange')), row=2, col=1)

fig.update_yaxes(title_text='Charge State [MWh]', row=1, col=1)
fig.update_yaxes(title_text='Charge State [MWh]', row=2, col=1)
fig.update_layout(height=400, showlegend=False)
fig.show()

print(f'{len(segments)} segments:')
for i, seg in enumerate(segments):
    print(
        f'  {i + 1}: {seg.timesteps[0]:%m-%d %H:%M} → {seg.timesteps[-1]:%m-%d %H:%M} | {seg.solution["costs"].item():,.0f} €'
    )

# Concatenate all segment solutions into one dataset (including overlaps)
ds = xr.concat([seg.solution for seg in segments], dim=pd.RangeIndex(len(segments), name='segment'), join='outer')

# Plot CHP thermal flow across all segments - each segment as a separate line
px.line(
    ds['Boiler(Q_th)|flow_rate'].to_pandas().T,
    labels={'value': 'Boiler Thermal Output [MW]', 'index': 'Timestep'},
)

px.line(
    ds['Storage|charge_state'].to_pandas().T,
    labels={'value': 'Storage Charge State [MW]', 'index': 'Timestep'},
)