Scenarios¶

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

import flixopt as fx

fx.CONFIG.notebook()

# Time horizon: one representative winter week
timesteps = pd.date_range('2024-01-15', periods=168, freq='h')  # 7 days

# Planning periods (years)
periods = pd.Index([2024, 2025, 2026], name='period')

# Scenarios with probabilities
scenarios = pd.Index(['Mild Winter', 'Harsh Winter'], name='scenario')
scenario_weights = np.array([0.6, 0.4])  # 60% mild, 40% harsh

print(f'Time dimension: {len(timesteps)} hours')
print(f'Periods: {list(periods)}')
print(f'Scenarios: {list(scenarios)}')
print(f'Scenario weights: {dict(zip(scenarios, scenario_weights, strict=False))}')

hours = np.arange(168)
hour_of_day = hours % 24

# Base daily pattern (kW): higher in morning/evening
daily_pattern = np.select(
    [
        (hour_of_day >= 6) & (hour_of_day < 9),  # Morning peak
        (hour_of_day >= 9) & (hour_of_day < 17),  # Daytime
        (hour_of_day >= 17) & (hour_of_day < 22),  # Evening peak
    ],
    [180, 120, 160],
    default=100,  # Night
).astype(float)

# Add random variation
np.random.seed(42)
noise = np.random.normal(0, 10, len(timesteps))

# Mild winter: lower demand
mild_demand = daily_pattern * 0.8 + noise
mild_demand = np.clip(mild_demand, 60, 200)

# Harsh winter: higher demand
harsh_demand = daily_pattern * 1.3 + noise * 1.5
harsh_demand = np.clip(harsh_demand, 100, 280)

# Create DataFrame with scenario columns (flixopt uses column names to match scenarios)
heat_demand = pd.DataFrame(
    {
        'Mild Winter': mild_demand,
        'Harsh Winter': harsh_demand,
    },
    index=timesteps,
)

print(f'Mild winter demand: {mild_demand.min():.0f} - {mild_demand.max():.0f} kW')
print(f'Harsh winter demand: {harsh_demand.min():.0f} - {harsh_demand.max():.0f} kW')

# Visualize demand scenarios with plotly
fig = px.line(
    heat_demand.iloc[:48],
    title='Heat Demand by Scenario (First 2 Days)',
    labels={'index': 'Time', 'value': 'kW', 'variable': 'Scenario'},
)
fig.update_traces(mode='lines')
fig

# Gas prices by period (€/kWh) - expected to rise
gas_prices = np.array([0.06, 0.08, 0.10])  # 2024, 2025, 2026

# Electricity sell prices by period (€/kWh) - CHP revenue
elec_prices = np.array([0.28, 0.34, 0.43])  # Rising with gas

print('Gas prices by period:')
for period, price in zip(periods, gas_prices, strict=False):
    print(f'  {period}: {price:.2f} €/kWh')

print('\nElectricity sell prices by period:')
for period, price in zip(periods, elec_prices, strict=False):
    print(f'  {period}: {price:.2f} €/kWh')

flow_system = fx.FlowSystem(
    timesteps=timesteps,
    periods=periods,
    scenarios=scenarios,
    scenario_weights=scenario_weights,
)
flow_system.add_carriers(
    fx.Carrier('gas', '#3498db', 'kW'),
    fx.Carrier('electricity', '#f1c40f', 'kW'),
    fx.Carrier('heat', '#e74c3c', 'kW'),
)

print(flow_system)

chp_size = flow_system.statistics.sizes['CHP(P_el)']
total_cost = flow_system.solution['costs']

print(f'Optimal CHP: {float(chp_size.max()):.0f} kW electrical ({float(chp_size.max()) * 0.50 / 0.35:.0f} kW thermal)')
print(f'Expected cost: {float(total_cost.sum()):.0f} €')

flow_system.statistics.plot.balance('Heat')

flow_system.statistics.plot.heatmap('CHP(Q_th)')

# View dimensions
flow_rates = flow_system.statistics.flow_rates
print('Flow rates dimensions:', dict(flow_rates.sizes))

# Plot flow rates
flow_system.statistics.plot.flows()

# CHP operation summary by scenario
chp_heat = flow_rates['CHP(Q_th)']

for scenario in scenarios:
    scenario_avg = float(chp_heat.sel(scenario=scenario).mean())
    scenario_max = float(chp_heat.sel(scenario=scenario).max())
    print(f'{scenario}: avg {scenario_avg:.0f} kW, max {scenario_max:.0f} kW')

# Select only the mild winter scenario
fs_mild = flow_system.transform.sel(scenario='Mild Winter')
fs_mild.optimize(fx.solvers.HighsSolver(mip_gap=0.01))

chp_size_mild = float(fs_mild.statistics.sizes['CHP(P_el)'].max())
chp_size_both = float(chp_size.max())

print(
    f'CHP sizing: {chp_size_mild:.0f} kW (mild only) vs {chp_size_both:.0f} kW (both scenarios) → +{chp_size_both - chp_size_mild:.0f} kW for uncertainty'
)

flow_system.statistics.plot.sankey.flows()