1-Node Capacity Expansion¶

In this example, we build a replica of model.energy. This tool calculates the cost of meeting a constant electricity demand from a combination of wind power, solar power and storage for different regions of the world. It includes capacity investments and dispatch optimisation. We deviate from model.energy by including an electricity demand profiles rather than a constant electricity demand.

In [2]:

import pandas as pd

import pypsa

import pandas as pd import pypsa

Techno-economic assumptions¶

We take techno-economic assumptions from the technology-data repository which collects assumptions on costs and efficiencies:

In [3]:

YEAR = 2030
url = f"https://raw.githubusercontent.com/PyPSA/technology-data/master/outputs/costs_{YEAR}.csv"
costs = pd.read_csv(url, index_col=[0, 1])
costs.loc[costs.unit.str.contains("/kW"), "value"] *= 1e3
costs = costs.value.unstack().fillna({"discount rate": 0.07, "lifetime": 20, "FOM": 0})

YEAR = 2030 url = f"https://raw.githubusercontent.com/PyPSA/technology-data/master/outputs/costs_{YEAR}.csv" costs = pd.read_csv(url, index_col=[0, 1]) costs.loc[costs.unit.str.contains("/kW"), "value"] *= 1e3 costs = costs.value.unstack().fillna({"discount rate": 0.07, "lifetime": 20, "FOM": 0})

Based on this, we calculate the marginal costs (€/MWh):

In [4]:

costs["marginal_cost"] = costs["VOM"] + costs["fuel"] / costs["efficiency"]

costs["marginal_cost"] = costs["VOM"] + costs["fuel"] / costs["efficiency"]

We also calculate the capital costs (i.e. annualised investment costs, €/MW/a or €/MWh/a for storage), using a small utility function to calculate the annuity factor to annualise investment costs based on is the discount rate $r$ and lifetime $n$.

In [5]:

from pypsa.costs import annuity

from pypsa.costs import annuity

In [6]:

a = costs.apply(lambda x: annuity(x["discount rate"], x["lifetime"]), axis=1)

a = costs.apply(lambda x: annuity(x["discount rate"], x["lifetime"]), axis=1)

In [7]:

costs["capital_cost"] = (a + costs["FOM"] / 100) * costs["investment"]

costs["capital_cost"] = (a + costs["FOM"] / 100) * costs["investment"]

Wind, solar and load time series¶

In [8]:

RESOLUTION = 3  # hours
url = "https://tubcloud.tu-berlin.de/s/9toBssWEdaLgHzq/download/time-series.csv"
ts = pd.read_csv(url, index_col=0, parse_dates=True)[::RESOLUTION]

RESOLUTION = 3 # hours url = "https://tubcloud.tu-berlin.de/s/9toBssWEdaLgHzq/download/time-series.csv" ts = pd.read_csv(url, index_col=0, parse_dates=True)[::RESOLUTION]

In [9]:

ts.head(3)

ts.head(3)

Out[9]:

	load_mw	pv_pu	wind_pu
timestamp
2019-01-01 00:00:00	5719.26	0.0	0.1846
2019-01-01 03:00:00	5474.74	0.0	0.3146
2019-01-01 06:00:00	5413.39	0.0	0.4957

Model initialisation¶

In [10]:

n = pypsa.Network()
n.add("Bus", "electricity", carrier="electricity")
n.set_snapshots(ts.index)

n = pypsa.Network() n.add("Bus", "electricity", carrier="electricity") n.set_snapshots(ts.index)

The weighting of the snapshots (e.g. how many hours they represent, see $w_t$ in problem formulation above) must be set in n.snapshot_weightings.

In [11]:

n.snapshot_weightings.loc[:, :] = RESOLUTION

n.snapshot_weightings.loc[:, :] = RESOLUTION

Adding carriers (this is only necessary for convenience in plotting later):

In [12]:

carriers = [
    "wind",
    "solar",
    "hydrogen storage",
    "battery storage",
    "load shedding",
    "electrolysis",
    "turbine",
    "electricity",
    "hydrogen",
]
colors = [
    "dodgerblue",
    "gold",
    "black",
    "yellowgreen",
    "darkorange",
    "magenta",
    "red",
    "grey",
    "grey",
]
n.add("Carrier", carriers, color=colors)

carriers = [ "wind", "solar", "hydrogen storage", "battery storage", "load shedding", "electrolysis", "turbine", "electricity", "hydrogen", ] colors = [ "dodgerblue", "gold", "black", "yellowgreen", "darkorange", "magenta", "red", "grey", "grey", ] n.add("Carrier", carriers, color=colors)

Adding load:

In [13]:

n.add(
    "Load",
    "demand",
    bus="electricity",
    p_set=ts.load_mw,
)

n.add( "Load", "demand", bus="electricity", p_set=ts.load_mw, )

Add a load shedding generator with high marginal cost of 2000 €/MWh:

In [14]:

n.add(
    "Generator",
    "load shedding",
    bus="electricity",
    carrier="load shedding",
    marginal_cost=2000,
    p_nom=ts.load_mw.max(),
)

n.add( "Generator", "load shedding", bus="electricity", carrier="load shedding", marginal_cost=2000, p_nom=ts.load_mw.max(), )

Adding the variable renewable generators works almost identically, but we also need to supply the capacity factors to the model via the attribute p_max_pu.

In [15]:

n.add(
    "Generator",
    "wind",
    bus="electricity",
    carrier="wind",
    p_max_pu=ts.wind_pu,
    capital_cost=costs.at["onwind", "capital_cost"],
    marginal_cost=costs.at["onwind", "marginal_cost"],
    p_nom_extendable=True,
)

n.add( "Generator", "wind", bus="electricity", carrier="wind", p_max_pu=ts.wind_pu, capital_cost=costs.at["onwind", "capital_cost"], marginal_cost=costs.at["onwind", "marginal_cost"], p_nom_extendable=True, )

In [16]:

n.add(
    "Generator",
    "solar",
    bus="electricity",
    carrier="solar",
    p_max_pu=ts.pv_pu,
    capital_cost=costs.at["solar", "capital_cost"],
    marginal_cost=costs.at["solar", "marginal_cost"],
    p_nom_extendable=True,
)

n.add( "Generator", "solar", bus="electricity", carrier="solar", p_max_pu=ts.pv_pu, capital_cost=costs.at["solar", "capital_cost"], marginal_cost=costs.at["solar", "marginal_cost"], p_nom_extendable=True, )

Adding a 3-hour battery:

In [17]:

n.add(
    "StorageUnit",
    "battery storage",
    bus="electricity",
    carrier="battery storage",
    max_hours=3,
    capital_cost=costs.at["battery inverter", "capital_cost"]
    + 3 * costs.at["battery storage", "capital_cost"],
    efficiency_store=costs.at["battery inverter", "efficiency"],
    efficiency_dispatch=costs.at["battery inverter", "efficiency"],
    p_nom_extendable=True,
    cyclic_state_of_charge=True,
)

n.add( "StorageUnit", "battery storage", bus="electricity", carrier="battery storage", max_hours=3, capital_cost=costs.at["battery inverter", "capital_cost"] + 3 * costs.at["battery storage", "capital_cost"], efficiency_store=costs.at["battery inverter", "efficiency"], efficiency_dispatch=costs.at["battery inverter", "efficiency"], p_nom_extendable=True, cyclic_state_of_charge=True, )

Adding a hydrogen storage system consisting of electrolyser, hydrogen turbine and underground storage, which can all be flexibly optimised:

In [18]:

n.add("Bus", "hydrogen", carrier="hydrogen")

n.add("Bus", "hydrogen", carrier="hydrogen")

In [19]:

n.add(
    "Link",
    "electrolysis",
    bus0="electricity",
    bus1="hydrogen",
    carrier="electrolysis",
    p_nom_extendable=True,
    efficiency=costs.at["electrolysis", "efficiency"],
    capital_cost=costs.at["electrolysis", "capital_cost"],
)

n.add( "Link", "electrolysis", bus0="electricity", bus1="hydrogen", carrier="electrolysis", p_nom_extendable=True, efficiency=costs.at["electrolysis", "efficiency"], capital_cost=costs.at["electrolysis", "capital_cost"], )

In [20]:

n.add(
    "Link",
    "turbine",
    bus0="hydrogen",
    bus1="electricity",
    carrier="turbine",
    p_nom_extendable=True,
    efficiency=costs.at["OCGT", "efficiency"],
    capital_cost=costs.at["OCGT", "capital_cost"] / costs.at["OCGT", "efficiency"],
)

n.add( "Link", "turbine", bus0="hydrogen", bus1="electricity", carrier="turbine", p_nom_extendable=True, efficiency=costs.at["OCGT", "efficiency"], capital_cost=costs.at["OCGT", "capital_cost"] / costs.at["OCGT", "efficiency"], )

In [21]:

n.add(
    "Store",
    "hydrogen storage",
    bus="hydrogen",
    carrier="hydrogen storage",
    capital_cost=costs.at["hydrogen storage underground", "capital_cost"],
    e_nom_extendable=True,
    e_cyclic=True,
)

n.add( "Store", "hydrogen storage", bus="hydrogen", carrier="hydrogen storage", capital_cost=costs.at["hydrogen storage underground", "capital_cost"], e_nom_extendable=True, e_cyclic=True, )

Model run¶

In [22]:

n.optimize()

n.optimize()

/tmp/ipykernel_3500/1261279110.py:1: FutureWarning: The default value of `include_objective_constant` will change from True to False in version 2.0. Set `include_objective_constant` explicitly to suppress this warning. Using False improves LP numerical conditioning by not including the objective constant as a variable.
  n.optimize()

INFO:linopy.model: Solve problem using Highs solver

INFO:linopy.model:Solver options:
 - log_to_console: False

INFO:linopy.io:Writing objective.

Writing constraints.:   0%|          | 0/21 [00:00<?, ?it/s]

Writing constraints.:  90%|█████████ | 19/21 [00:00<00:00, 188.95it/s]

Writing constraints.: 100%|██████████| 21/21 [00:00<00:00, 182.03it/s]

Writing continuous variables.:   0%|          | 0/11 [00:00<?, ?it/s]

Writing continuous variables.: 100%|██████████| 11/11 [00:00<00:00, 422.57it/s]

INFO:linopy.io: Writing time: 0.16s

INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 29206 primals, 64246 duals
Objective: 1.01e+10
Solver model: available
Solver message: Optimal

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Generator-ext-p-lower, Generator-ext-p-upper, Link-ext-p-lower, Link-ext-p-upper, Store-ext-e-lower, Store-ext-e-upper, StorageUnit-ext-p_dispatch-lower, StorageUnit-ext-p_dispatch-upper, StorageUnit-ext-p_store-lower, StorageUnit-ext-p_store-upper, StorageUnit-ext-state_of_charge-lower, StorageUnit-ext-state_of_charge-upper, StorageUnit-energy_balance, Store-energy_balance were not assigned to the network.

Out[22]:

('ok', 'optimal')

Model evaluation¶

Total system cost by technology:

In [23]:

tsc = (
    pd.concat([n.statistics.capex(), n.statistics.opex()], axis=1).sum(axis=1).div(1e9)
)
tsc

tsc = ( pd.concat([n.statistics.capex(), n.statistics.opex()], axis=1).sum(axis=1).div(1e9) ) tsc

Out[23]:

component    carrier         
Generator    solar               1.684343
             wind                4.180646
Link         electrolysis        0.739954
             turbine             1.424007
StorageUnit  battery storage     1.124633
Store        hydrogen storage    0.688179
Generator    load shedding       0.290825
dtype: float64

In [24]:

tsc.sum()

tsc.sum()

Out[24]:

np.float64(10.132586933378578)

The optimised capacities in GW (GWh for Store component):

In [25]:

n.statistics.optimal_capacity().div(1e3)

n.statistics.optimal_capacity().div(1e3)

Out[25]:

component    carrier         
Generator    load shedding         10.901160
             solar                 26.089521
             wind                  32.583318
Link         electrolysis           3.118498
             turbine                9.692634
StorageUnit  battery storage       14.047947
Store        hydrogen storage    3672.190855
dtype: float64

Energy balances on electricity side (in TWh):

In [26]:

n.statistics.energy_balance(bus_carrier="electricity").sort_values().div(1e6)

n.statistics.energy_balance(bus_carrier="electricity").sort_values().div(1e6)

Out[26]:

component    carrier          bus_carrier
Load         -                electricity   -66.266089
Link         electrolysis     electricity   -15.656979
StorageUnit  battery storage  electricity    -0.731487
Generator    load shedding    electricity     0.145412
Link         turbine          electricity     3.990917
Generator    solar            electricity    25.926729
             wind             electricity    52.591497
dtype: float64

Energy balances plot as time series (in MW):

In [27]:

n.statistics.energy_balance.plot.area(linewidth=0, bus_carrier="electricity")

n.statistics.energy_balance.plot.area(linewidth=0, bus_carrier="electricity")

Out[27]:

(<Figure size 758.25x300 with 1 Axes>,
 <Axes: xlabel='snapshot', ylabel='Energy Balance []'>,
 <seaborn.axisgrid.FacetGrid at 0x7f34c273de80>)

Price time series for electricity and hydrogen:

In [28]:

n.buses_t.marginal_price.plot(figsize=(7, 2))

n.buses_t.marginal_price.plot(figsize=(7, 2))

Out[28]:

<Axes: xlabel='snapshot'>