Multi-Carrier¶

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

import flixopt as fx

fx.CONFIG.notebook()

# One week, hourly
timesteps = pd.date_range('2024-02-05', periods=168, freq='h')
hours = np.arange(168)
hour_of_day = hours % 24

# Hospital electricity demand (kW)
# Base load + daily pattern (higher during day for equipment, lighting)
elec_base = 150  # 24/7 critical systems
elec_daily = 100 * np.sin((hour_of_day - 6) * np.pi / 12)  # Peak at noon
elec_daily = np.maximum(0, elec_daily)
electricity_demand = elec_base + elec_daily

# Hospital heat demand (kW)
# Higher in morning, drops during day, increases for hot water in evening
heat_pattern = np.select(
    [
        (hour_of_day >= 5) & (hour_of_day < 9),  # Morning warmup
        (hour_of_day >= 9) & (hour_of_day < 17),  # Daytime
        (hour_of_day >= 17) & (hour_of_day < 22),  # Evening
    ],
    [350, 250, 300],
    default=200,  # Night
)
heat_demand = heat_pattern.astype(float)

# Add random variation
np.random.seed(456)
electricity_demand += np.random.normal(0, 15, len(timesteps))
heat_demand += np.random.normal(0, 20, len(timesteps))
electricity_demand = np.clip(electricity_demand, 100, 300)
heat_demand = np.clip(heat_demand, 150, 400)

print(f'Electricity: {electricity_demand.min():.0f} - {electricity_demand.max():.0f} kW')
print(f'Heat: {heat_demand.min():.0f} - {heat_demand.max():.0f} kW')

# Visualize demands and prices with plotly - using xarray and faceting
profiles = xr.Dataset(
    {
        'Electricity Demand [kW]': xr.DataArray(electricity_demand, dims=['time'], coords={'time': timesteps}),
        'Heat Demand [kW]': xr.DataArray(heat_demand, dims=['time'], coords={'time': timesteps}),
        'Elec. Buy Price [€/kWh]': xr.DataArray(elec_buy_price, dims=['time'], coords={'time': timesteps}),
    }
)

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

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

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

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

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

total_costs = flow_system.solution['costs'].item()
total_co2 = flow_system.solution['CO2'].item()

# Energy flows
flow_rates = flow_system.statistics.flow_rates
grid_buy = flow_rates['GridBuy(Electricity)'].sum().item()
grid_sell = flow_rates['GridSell(Electricity)'].sum().item()
chp_elec = flow_rates['CHP(P_el)'].sum().item()
chp_heat = flow_rates['CHP(Q_th)'].sum().item()
boiler_heat = flow_rates['Boiler(Q_th)'].sum().item()

total_elec = electricity_demand.sum()
total_heat = heat_demand.sum()

# Display as compact summary
print(
    f'Electricity: {chp_elec:.0f} kWh CHP ({chp_elec / total_elec * 100:.0f}%) + {grid_buy:.0f} kWh grid, {grid_sell:.0f} kWh sold'
)
print(f'Heat: {chp_heat:.0f} kWh CHP ({chp_heat / total_heat * 100:.0f}%) + {boiler_heat:.0f} kWh boiler')
print(f'Costs: {total_costs:.2f} € | CO2: {total_co2:.0f} kg')

# Build system without CHP
fs_no_chp = fx.FlowSystem(timesteps)
fs_no_chp.add_carriers(
    fx.Carrier('gas', '#3498db', 'kW'),
    fx.Carrier('electricity', '#f1c40f', 'kW'),
    fx.Carrier('heat', '#e74c3c', 'kW'),
)
fs_no_chp.add_elements(
    fx.Bus('Electricity', carrier='electricity'),
    fx.Bus('Heat', carrier='heat'),
    fx.Bus('Gas', carrier='gas'),
    fx.Effect('costs', '€', 'Total Costs', is_standard=True, is_objective=True),
    fx.Effect('CO2', 'kg', 'CO2 Emissions'),
    fx.Source(
        'GasGrid',
        outputs=[fx.Flow('Gas', bus='Gas', size=1000, effects_per_flow_hour={'costs': gas_price, 'CO2': 0.2})],
    ),
    fx.Source(
        'GridBuy',
        outputs=[
            fx.Flow(
                'Electricity', bus='Electricity', size=500, effects_per_flow_hour={'costs': elec_buy_price, 'CO2': 0.4}
            )
        ],
    ),
    # Only boiler for heat
    fx.linear_converters.Boiler(
        'Boiler',
        thermal_efficiency=0.92,
        thermal_flow=fx.Flow('Q_th', bus='Heat', size=500),
        fuel_flow=fx.Flow('Q_fuel', bus='Gas'),
    ),
    fx.Sink(
        'HospitalElec', inputs=[fx.Flow('Load', bus='Electricity', size=1, fixed_relative_profile=electricity_demand)]
    ),
    fx.Sink('HospitalHeat', inputs=[fx.Flow('Load', bus='Heat', size=1, fixed_relative_profile=heat_demand)]),
)

fs_no_chp.optimize(fx.solvers.HighsSolver())

no_chp_costs = fs_no_chp.solution['costs'].item()
no_chp_co2 = fs_no_chp.solution['CO2'].item()

cost_saving = (no_chp_costs - total_costs) / no_chp_costs * 100
co2_saving = (no_chp_co2 - total_co2) / no_chp_co2 * 100
print(
    f'CHP saves {cost_saving:.1f}% costs ({no_chp_costs:.0f}→{total_costs:.0f} €) and {co2_saving:.1f}% CO2 ({no_chp_co2:.0f}→{total_co2:.0f} kg)'
)

flow_system.statistics.plot.sankey.flows()