Investment: Planning a Heating Plant Replacement¶

Investment models a singular build decision: the solver decides when to build and at what size. Once built, capacity is available for a limited number of periods (the lifetime), then retires.

This notebook demonstrates a realistic energy transition scenario: an aging CHP plant serves a district heating network, but rising gas prices and falling electricity prices shift the economics toward a heat pump. The solver decides when and what to build.

Key concepts: Investment with prior_size and lifetime · period-varying effects_per_flow_hour · technology transition dynamics

In [1]:

from datetime import datetime

import numpy as np
import pandas as pd
import plotly.graph_objects as go
import plotly.io as pio
import xarray as xr
from plotly.subplots import make_subplots

from fluxopt import Carrier, Converter, Effect, Flow, Investment, Port, optimize

pio.renderers.default = 'notebook_connected'

from datetime import datetime import numpy as np import pandas as pd import plotly.graph_objects as go import plotly.io as pio import xarray as xr from plotly.subplots import make_subplots from fluxopt import Carrier, Converter, Effect, Flow, Investment, Port, optimize pio.renderers.default = 'notebook_connected'

Setup¶

A district needs heat across four 5-year planning periods (2025--2040). Energy prices shift as the energy transition progresses:

Period	Gas (EUR/MWh)	Electricity (EUR/MWh)	Story
2025	30	120	Cheap gas era
2030	50	100	Carbon tax kicks in
2035	70	75	Approaching price parity
2040	90	60	Renewables dominate

We use 4 representative timesteps per period with dt=500 (each timestep represents 500 operating hours). NPV weights at 5% discount rate make the build-timing decision deterministic.

In [2]:

timesteps = [datetime(2025, 1, 15, h) for h in range(6, 10)]
periods = [2025, 2030, 2035, 2040]
dt = 500  # each representative timestep = 500 operating hours

# Heat demand profile [MW] -- 4 representative hours
demand = [8, 10, 9, 7]
max_demand = max(demand)

# NPV weights: 5% discount rate, 5-year blocks
r = 0.05
npv_weights = [round(sum(1 / (1 + r) ** (y0 + y) for y in range(5)), 2) for y0 in [0, 5, 10, 15]]

# Energy prices [EUR/MWh] -- vary by period
gas_price = xr.DataArray([30, 50, 70, 90], dims=['period'], coords={'period': periods})
elec_buy_price = xr.DataArray([120, 100, 75, 60], dims=['period'], coords={'period': periods})
elec_sell_price = elec_buy_price * 0.9  # feed-in at 90% of buy price

pd.DataFrame(
    {
        'gas': gas_price.values,
        'elec (buy)': elec_buy_price.values,
        'elec (sell)': elec_sell_price.values,
        'NPV weight': npv_weights,
    },
    index=pd.Index(periods, name='period'),
)

timesteps = [datetime(2025, 1, 15, h) for h in range(6, 10)] periods = [2025, 2030, 2035, 2040] dt = 500 # each representative timestep = 500 operating hours # Heat demand profile [MW] -- 4 representative hours demand = [8, 10, 9, 7] max_demand = max(demand) # NPV weights: 5% discount rate, 5-year blocks r = 0.05 npv_weights = [round(sum(1 / (1 + r) ** (y0 + y) for y in range(5)), 2) for y0 in [0, 5, 10, 15]] # Energy prices [EUR/MWh] -- vary by period gas_price = xr.DataArray([30, 50, 70, 90], dims=['period'], coords={'period': periods}) elec_buy_price = xr.DataArray([120, 100, 75, 60], dims=['period'], coords={'period': periods}) elec_sell_price = elec_buy_price * 0.9 # feed-in at 90% of buy price pd.DataFrame( { 'gas': gas_price.values, 'elec (buy)': elec_buy_price.values, 'elec (sell)': elec_sell_price.values, 'NPV weight': npv_weights, }, index=pd.Index(periods, name='period'), )

Out[2]:

	gas	elec (buy)	elec (sell)	NPV weight
period
2025	30	120	108.0	4.55
2030	50	100	90.0	3.56
2035	70	75	67.5	2.79
2040	90	60	54.0	2.19

In [3]:

x = [str(p) for p in periods]

fig = make_subplots(rows=1, cols=2, subplot_titles=('Energy Prices', 'Effective Heat Cost'))

# Left: raw carrier prices
fig.add_trace(
    go.Scatter(
        x=x,
        y=gas_price.values,
        name='Gas',
        mode='lines+markers',
        line={'color': '#ef553b', 'width': 3},
    ),
    row=1,
    col=1,
)
fig.add_trace(
    go.Scatter(
        x=x,
        y=elec_buy_price.values,
        name='Electricity (buy)',
        mode='lines+markers',
        line={'color': '#636efa', 'width': 3},
    ),
    row=1,
    col=1,
)
fig.add_trace(
    go.Scatter(
        x=x,
        y=elec_sell_price.values,
        name='Electricity (sell)',
        mode='lines+markers',
        line={'color': '#636efa', 'width': 2, 'dash': 'dash'},
    ),
    row=1,
    col=1,
)

# Right: effective cost per MWh of heat delivered
boiler_eff = gas_price.values / 0.90
hp_eff = elec_buy_price.values / 3.0
fig.add_trace(
    go.Scatter(
        x=x,
        y=boiler_eff,
        name='Boiler (gas/0.9)',
        mode='lines+markers',
        line={'color': '#ef553b', 'width': 2, 'dash': 'dot'},
    ),
    row=1,
    col=2,
)
fig.add_trace(
    go.Scatter(
        x=x,
        y=hp_eff,
        name='Heat pump (elec/3.0)',
        mode='lines+markers',
        line={'color': '#636efa', 'width': 2, 'dash': 'dot'},
    ),
    row=1,
    col=2,
)

fig.update_layout(
    height=350,
    template='plotly_white',
    margin={'l': 50, 'r': 20, 't': 40, 'b': 20},
    legend={'x': 0.01, 'y': 0.98},
)
fig.update_yaxes(title_text='EUR / MWh', row=1, col=1)
fig.update_yaxes(title_text='EUR / MWh_th', row=1, col=2)
fig

x = [str(p) for p in periods] fig = make_subplots(rows=1, cols=2, subplot_titles=('Energy Prices', 'Effective Heat Cost')) # Left: raw carrier prices fig.add_trace( go.Scatter( x=x, y=gas_price.values, name='Gas', mode='lines+markers', line={'color': '#ef553b', 'width': 3}, ), row=1, col=1, ) fig.add_trace( go.Scatter( x=x, y=elec_buy_price.values, name='Electricity (buy)', mode='lines+markers', line={'color': '#636efa', 'width': 3}, ), row=1, col=1, ) fig.add_trace( go.Scatter( x=x, y=elec_sell_price.values, name='Electricity (sell)', mode='lines+markers', line={'color': '#636efa', 'width': 2, 'dash': 'dash'}, ), row=1, col=1, ) # Right: effective cost per MWh of heat delivered boiler_eff = gas_price.values / 0.90 hp_eff = elec_buy_price.values / 3.0 fig.add_trace( go.Scatter( x=x, y=boiler_eff, name='Boiler (gas/0.9)', mode='lines+markers', line={'color': '#ef553b', 'width': 2, 'dash': 'dot'}, ), row=1, col=2, ) fig.add_trace( go.Scatter( x=x, y=hp_eff, name='Heat pump (elec/3.0)', mode='lines+markers', line={'color': '#636efa', 'width': 2, 'dash': 'dot'}, ), row=1, col=2, ) fig.update_layout( height=350, template='plotly_white', margin={'l': 50, 'r': 20, 't': 40, 'b': 20}, legend={'x': 0.01, 'y': 0.98}, ) fig.update_yaxes(title_text='EUR / MWh', row=1, col=1) fig.update_yaxes(title_text='EUR / MWh_th', row=1, col=2) fig

The lines cross: gas heating gets more expensive while electric heating gets cheaper. By 2030 the heat pump is already cheaper than the boiler per MWh of heat delivered.

System¶

Three technologies compete to supply district heat:

Technology	Efficiency	CAPEX	Lifetime	Notes
CHP	eta_el=0.35, eta_th=0.50	--	2 periods	Existing (prior_size=5 MW)
Gas boiler	eta=0.90	20k EUR/MW	2 periods	Investable
Heat pump	COP=3.0	200k EUR/MW	3 periods	Investable

The CHP co-generates electricity, which offsets grid imports while gas is cheap. After 2 periods it retires, and the solver must fill the gap.

Prices vary by period using xr.DataArray with a period dimension in effects_per_flow_hour:

In [4]:

result = optimize(
    timesteps=timesteps,
    carriers=[
        Carrier('gas'),
        Carrier('electricity'),
        Carrier('heat'),
        Carrier('ambient'),
    ],
    effects=[
        Effect('cost', unit='EUR', period_weights=npv_weights),
    ],
    ports=[
        Port(
            'gas_grid',
            imports=[
                Flow('gas', size=200, effects_per_flow_hour={'cost': gas_price}),
            ],
        ),
        Port(
            'elec_grid',
            imports=[Flow('electricity', short_id='buy', size=100, effects_per_flow_hour={'cost': elec_buy_price})],
            exports=[Flow('electricity', short_id='sell', size=100, effects_per_flow_hour={'cost': -elec_sell_price})],
        ),
        Port('ambient', imports=[Flow('ambient', size=100)]),
        Port(
            'district',
            exports=[
                Flow('heat', size=1, fixed_relative_profile=demand),
            ],
        ),
    ],
    converters=[
        Converter.chp(
            'chp',
            eta_el=0.35,
            eta_th=0.50,
            fuel_flow=Flow('gas', size=200),
            electrical_flow=Flow('electricity', size=100),
            thermal_flow=Flow(
                'heat',
                size=Investment(
                    0,
                    5,
                    mandatory=False,
                    prior_size=5,
                    lifetime=2,
                ),
            ),
        ),
        Converter.boiler(
            'boiler',
            thermal_efficiency=0.90,
            fuel_flow=Flow('gas', size=200),
            thermal_flow=Flow(
                'heat',
                size=Investment(
                    1,
                    15,
                    mandatory=False,
                    lifetime=2,
                    effects_per_size_at_build={'cost': 20_000},
                ),
            ),
        ),
        Converter.heat_pump(
            'hp',
            cop=3.0,
            electrical_flow=Flow('electricity', size=100),
            source_flow=Flow('ambient', size=100),
            thermal_flow=Flow(
                'heat',
                size=Investment(
                    1,
                    15,
                    mandatory=False,
                    lifetime=3,
                    effects_per_size_at_build={'cost': 200_000},
                ),
            ),
        ),
    ],
    periods=periods,
    period_weights=npv_weights,
    dt=dt,
    objective_effects='cost',
)

result = optimize( timesteps=timesteps, carriers=[ Carrier('gas'), Carrier('electricity'), Carrier('heat'), Carrier('ambient'), ], effects=[ Effect('cost', unit='EUR', period_weights=npv_weights), ], ports=[ Port( 'gas_grid', imports=[ Flow('gas', size=200, effects_per_flow_hour={'cost': gas_price}), ], ), Port( 'elec_grid', imports=[Flow('electricity', short_id='buy', size=100, effects_per_flow_hour={'cost': elec_buy_price})], exports=[Flow('electricity', short_id='sell', size=100, effects_per_flow_hour={'cost': -elec_sell_price})], ), Port('ambient', imports=[Flow('ambient', size=100)]), Port( 'district', exports=[ Flow('heat', size=1, fixed_relative_profile=demand), ], ), ], converters=[ Converter.chp( 'chp', eta_el=0.35, eta_th=0.50, fuel_flow=Flow('gas', size=200), electrical_flow=Flow('electricity', size=100), thermal_flow=Flow( 'heat', size=Investment( 0, 5, mandatory=False, prior_size=5, lifetime=2, ), ), ), Converter.boiler( 'boiler', thermal_efficiency=0.90, fuel_flow=Flow('gas', size=200), thermal_flow=Flow( 'heat', size=Investment( 1, 15, mandatory=False, lifetime=2, effects_per_size_at_build={'cost': 20_000}, ), ), ), Converter.heat_pump( 'hp', cop=3.0, electrical_flow=Flow('electricity', size=100), source_flow=Flow('ambient', size=100), thermal_flow=Flow( 'heat', size=Investment( 1, 15, mandatory=False, lifetime=3, effects_per_size_at_build={'cost': 200_000}, ), ), ), ], periods=periods, period_weights=npv_weights, dt=dt, objective_effects='cost', )

Running HiGHS 1.13.1 (git hash: 1d267d9): Copyright (c) 2026 under MIT licence terms
MIP linopy-problem-d1yb9h0y has 685 rows; 307 cols; 1216 nonzeros; 24 integer variables (24 binary)
Coefficient ranges:
  Matrix  [3e-01, 2e+05]
  Cost    [2e+00, 5e+00]
  Bound   [1e+00, 2e+01]
  RHS     [1e+00, 2e+02]
Presolving model
124 rows, 104 cols, 328 nonzeros  0s
110 rows, 82 cols, 300 nonzeros  0s
90 rows, 69 cols, 228 nonzeros  0s
87 rows, 62 cols, 215 nonzeros  0s
Presolve reductions: rows 87(-598); columns 62(-245); nonzeros 215(-1001) 

Solving MIP model with:
   87 rows
   62 cols (3 binary, 0 integer, 0 implied int., 59 continuous, 0 domain fixed)
   215 nonzeros

Src: B => Branching; C => Central rounding; F => Feasibility pump; H => Heuristic;
     I => Shifting; J => Feasibility jump; L => Sub-MIP; P => Empty MIP; R => Randomized rounding;
     S => Solve LP; T => Evaluate node; U => Unbounded; X => User solution; Y => HiGHS solution;
     Z => ZI Round; l => Trivial lower; p => Trivial point; u => Trivial upper; z => Trivial zero

        Nodes      |    B&B Tree     |            Objective Bounds              |  Dynamic Constraints |       Work      
Src  Proc. InQueue |  Leaves   Expl. | BestBound       BestSol              Gap |   Cuts   InLp Confl. | LpIters     Time

         0       0         0   0.00%   -160429650      inf                  inf        0      0      0         0     0.0s
 R       0       0         0   0.00%   4868612.231183  11809474.73118    58.77%        0      0      0        38     0.0s
 H       0       0         0   0.00%   8558816.273095  11018861.11111    22.33%       96      8      0       129     0.0s
         1       0         1 100.00%   11018861.11111  11018861.11111     0.00%       96      8      2       168     0.0s

Solving report
  Model             linopy-problem-d1yb9h0y
  Status            Optimal
  Primal bound      11018861.1111
  Dual bound        11018861.1111
  Gap               0% (tolerance: 0.01%)
  P-D integral      0.00524592766546
  Solution status   feasible
                    11018861.1111 (objective)
                    0 (bound viol.)
                    2.22044604925e-16 (int. viol.)
                    0 (row viol.)
  Timing            0.03
  Max sub-MIP depth 1
  Nodes             1
  Repair LPs        0
  LP iterations     168
                    0 (strong br.)
                    25 (separation)
                    105 (heuristics)

/home/docs/checkouts/readthedocs.org/user_builds/fluxopt/envs/v0.0.8/lib/python3.13/site-packages/linopy/common.py:492: UserWarning: Coordinates across variables not equal. Perform outer join.
  warn(

Results¶

Investment Timeline¶

When does each technology get built, and how long is it active?

In [5]:

sol = result.solution

techs = {
    'CHP': 'chp(heat)',
    'Boiler': 'boiler(heat)',
    'Heat pump': 'hp(heat)',
}
colors = {'CHP': '#2ca02c', 'Boiler': '#ef553b', 'Heat pump': '#636efa'}

# Investment status table
rows = []
for label, fid in techs.items():
    active = sol['invest--active'].sel(flow=fid).values
    build = sol['invest--build'].sel(flow=fid).values
    size = float(sol['invest--size'].sel(flow=fid).values)
    for i, p in enumerate(periods):
        rows.append(
            {
                'Technology': label,
                'Period': p,
                'Active': bool(active[i] > 0.5),
                'Build': bool(build[i] > 0.5),
                'Size (MW)': round(size, 1),
            }
        )

pd.DataFrame(rows).pivot(index='Period', columns='Technology', values='Active')

sol = result.solution techs = { 'CHP': 'chp(heat)', 'Boiler': 'boiler(heat)', 'Heat pump': 'hp(heat)', } colors = {'CHP': '#2ca02c', 'Boiler': '#ef553b', 'Heat pump': '#636efa'} # Investment status table rows = [] for label, fid in techs.items(): active = sol['invest--active'].sel(flow=fid).values build = sol['invest--build'].sel(flow=fid).values size = float(sol['invest--size'].sel(flow=fid).values) for i, p in enumerate(periods): rows.append( { 'Technology': label, 'Period': p, 'Active': bool(active[i] > 0.5), 'Build': bool(build[i] > 0.5), 'Size (MW)': round(size, 1), } ) pd.DataFrame(rows).pivot(index='Period', columns='Technology', values='Active')

Out[5]:

Technology	Boiler	CHP	Heat pump
Period
2025	True	True	False
2030	True	True	False
2035	False	False	True
2040	False	False	True

In [6]:

fig = go.Figure()
x = [str(p) for p in periods]

for label, fid in techs.items():
    active = sol['invest--active'].sel(flow=fid).values
    build = sol['invest--build'].sel(flow=fid).values
    size = float(sol['invest--size'].sel(flow=fid).values)

    fig.add_trace(
        go.Bar(
            x=x,
            y=active * size,
            name=f'{label} ({size:.0f} MW)',
            marker_color=colors[label],
            opacity=0.7,
        )
    )

    # Mark build events with a star
    build_x = [str(periods[i]) for i in range(len(periods)) if build[i] > 0.5]
    if build_x:
        fig.add_trace(
            go.Scatter(
                x=build_x,
                y=[size] * len(build_x),
                mode='markers',
                marker={'size': 16, 'symbol': 'star', 'color': colors[label], 'line': {'width': 1.5, 'color': 'black'}},
                showlegend=False,
            )
        )

fig.update_layout(
    barmode='stack',
    height=350,
    template='plotly_white',
    title='Active Capacity by Period',
    margin={'l': 50, 'r': 20, 't': 50, 'b': 20},
    yaxis={'title': 'Capacity (MW)'},
    xaxis={'title': 'Period'},
)
fig

fig = go.Figure() x = [str(p) for p in periods] for label, fid in techs.items(): active = sol['invest--active'].sel(flow=fid).values build = sol['invest--build'].sel(flow=fid).values size = float(sol['invest--size'].sel(flow=fid).values) fig.add_trace( go.Bar( x=x, y=active * size, name=f'{label} ({size:.0f} MW)', marker_color=colors[label], opacity=0.7, ) ) # Mark build events with a star build_x = [str(periods[i]) for i in range(len(periods)) if build[i] > 0.5] if build_x: fig.add_trace( go.Scatter( x=build_x, y=[size] * len(build_x), mode='markers', marker={'size': 16, 'symbol': 'star', 'color': colors[label], 'line': {'width': 1.5, 'color': 'black'}}, showlegend=False, ) ) fig.update_layout( barmode='stack', height=350, template='plotly_white', title='Active Capacity by Period', margin={'l': 50, 'r': 20, 't': 50, 'b': 20}, yaxis={'title': 'Capacity (MW)'}, xaxis={'title': 'Period'}, ) fig

Heat Supply Mix¶

Who supplies the heat in each period? The stacked bars show average heat production by technology.

In [7]:

rates = sol['flow--rate']
x = [str(p) for p in periods]

fig = go.Figure()
for label, fid in techs.items():
    avg_rate = rates.sel(flow=fid).mean('time').values
    fig.add_trace(go.Bar(x=x, y=avg_rate, name=label, marker_color=colors[label]))

fig.update_layout(
    barmode='stack',
    height=350,
    template='plotly_white',
    title='Heat Supply Mix by Period',
    margin={'l': 50, 'r': 20, 't': 50, 'b': 20},
    yaxis={'title': 'Average Heat Output (MW)'},
    xaxis={'title': 'Period'},
)
fig

rates = sol['flow--rate'] x = [str(p) for p in periods] fig = go.Figure() for label, fid in techs.items(): avg_rate = rates.sel(flow=fid).mean('time').values fig.add_trace(go.Bar(x=x, y=avg_rate, name=label, marker_color=colors[label])) fig.update_layout( barmode='stack', height=350, template='plotly_white', title='Heat Supply Mix by Period', margin={'l': 50, 'r': 20, 't': 50, 'b': 20}, yaxis={'title': 'Average Heat Output (MW)'}, xaxis={'title': 'Period'}, ) fig

Cost Breakdown¶

In [8]:

totals = result.effect_totals.sel(effect='cost')

pd.DataFrame(
    {
        'cost/period (EUR)': totals.values.round(0),
        'NPV weight': npv_weights,
        'weighted (EUR)': (totals.values * np.array(npv_weights)).round(0),
    },
    index=pd.Index(periods, name='period'),
).assign(objective=result.objective)

totals = result.effect_totals.sel(effect='cost') pd.DataFrame( { 'cost/period (EUR)': totals.values.round(0), 'NPV weight': npv_weights, 'weighted (EUR)': (totals.values * np.array(npv_weights)).round(0), }, index=pd.Index(periods, name='period'), ).assign(objective=result.objective)

Out[8]:

	cost/period (EUR)	NPV weight	weighted (EUR)	objective
period
2025	177333.0	4.55	806867.0	1.101886e+07
2030	758889.0	3.56	2701644.0	1.101886e+07
2035	2425000.0	2.79	6765750.0	1.101886e+07
2040	340000.0	2.19	744600.0	1.101886e+07

Insights¶

1. CHP + boiler early: The existing CHP covers most demand in 2025 while gas is cheap and electricity sales generate revenue. A cheap gas boiler fills the remaining gap.

2. Heat pump takes over: From 2030, the heat pump dominates as electricity prices fall. The CHP and boiler retire after 2 periods, leaving the heat pump as the sole supplier.

3. CAPEX vs OPEX trade-off: The boiler has low CAPEX (20k/MW) but relies on increasingly expensive gas. The heat pump has high CAPEX (200k/MW) but benefits from falling electricity prices and a longer lifetime (3 periods) that amortizes the investment.

4. Period-varying prices: Passing xr.DataArray with a period dimension to effects_per_flow_hour lets the solver see the full price trajectory when making build decisions.

How it works¶

Investment adds three decision variables per flow:

Variable	Dims	Meaning
invest_size	(flow,)	Capacity chosen once (period-independent)
invest_build	(flow, period)	Binary: build in this period?
invest_active	(flow, period)	Binary: operational in this period?

The existing flow_size variable is linked to invest_size * active via big-M constraints, so it equals the invested capacity when active and zero when retired. Flow rate bounds (P <= flow_size * rel_ub) apply automatically.

Key constraints:

• Build-once: sum_p build[p] == 1 (mandatory) or <= 1 (optional)
• Active logic: active[p] = sum_{tau: p in [tau, tau+lifetime)} build[tau]
• Prior capacity: active[p] = 1 for p < lifetime when prior_size > 0

Key API patterns:

• Investment(prior_size=...) -- pre-existing capacity with limited lifetime
• Investment(lifetime=...) -- solver decides build period; capacity expires
• effects_per_size_at_build -- one-time CAPEX charged in the build period
• effects_per_flow_hour with xr.DataArray(period=...) -- period-varying prices
• NPV period_weights on Effect -- discount future costs