Islanded Fuel Production¶

In this example, we build a simple model to assess the economics of islanded methanol production in Namibia (or Argentina). The model can optimise investment and operation of wind and solar, electrolysers, turbine, hydrogen and battery storage, direct air capture, CO$_2$ storage, methanolisation and methanol stores to supply a constant methanol demand of 1 TWh/a.

In [2]:

import pandas as pd

import pypsa
from pypsa.costs import annuity

import pandas as pd import pypsa from pypsa.costs import annuity

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})

# Let's also take a little more optimistic view on the costs of electrolysers
costs.loc["electrolysis", "investment"] = 500  # €/kW

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}) # Let's also take a little more optimistic view on the costs of electrolysers costs.loc["electrolysis", "investment"] = 500 # €/kW

We 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 [4]:

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 [5]:

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

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

Wind and solar time series¶

The wind and solar capacity factor time series have been retrieved from model.energy. Go there to find more time series for other countries and plug them in here.

In [6]:

RESOLUTION = 4  # hours
url = "https://model.energy/data/time-series-ca2bcb9e843aeb286cd6295854c885b6.csv"  # South of Argentina
# url = "https://model.energy/data/time-series-57f7bbcb5c4821506de052e52d022b48.csv" # Namibia
ts = pd.read_csv(url, index_col=0, parse_dates=True)[::RESOLUTION]

RESOLUTION = 4 # hours url = "https://model.energy/data/time-series-ca2bcb9e843aeb286cd6295854c885b6.csv" # South of Argentina # url = "https://model.energy/data/time-series-57f7bbcb5c4821506de052e52d022b48.csv" # Namibia ts = pd.read_csv(url, index_col=0, parse_dates=True)[::RESOLUTION]

In [7]:

ts.head(3)

ts.head(3)

Out[7]:

	onwind	solar
2011-01-01 00:00:00	0.910	0.032
2011-01-01 04:00:00	0.831	0.000
2011-01-01 08:00:00	0.758	0.000

Build model¶

Initialisation¶

Add buses, carriers and set snapshots.

In [8]:

n = pypsa.Network()
for carrier in ["electricity", "hydrogen", "co2", "methanol"]:
    n.add("Bus", carrier, carrier=carrier, unit="t/h" if carrier == "co2" else "MW")
n.set_snapshots(ts.index)
n.snapshot_weightings.loc[:, :] = RESOLUTION

n = pypsa.Network() for carrier in ["electricity", "hydrogen", "co2", "methanol"]: n.add("Bus", carrier, carrier=carrier, unit="t/h" if carrier == "co2" else "MW") n.set_snapshots(ts.index) n.snapshot_weightings.loc[:, :] = RESOLUTION

In [9]:

carriers = {
    "wind": "dodgerblue",
    "solar": "gold",
    "hydrogen storage": "blueviolet",
    "battery storage 3h": "yellowgreen",
    "battery storage 6h": "yellowgreen",
    "electrolysis": "magenta",
    "turbine": "darkorange",
    "methanolisation": "cyan",
    "direct air capture": "coral",
    "co2 storage": "black",
    "methanol storage": "cadetblue",
    "electricity": "grey",
    "hydrogen": "grey",
    "co2": "grey",
    "methanol": "grey",
}
n.add("Carrier", carriers.keys(), color=carriers.values())

carriers = { "wind": "dodgerblue", "solar": "gold", "hydrogen storage": "blueviolet", "battery storage 3h": "yellowgreen", "battery storage 6h": "yellowgreen", "electrolysis": "magenta", "turbine": "darkorange", "methanolisation": "cyan", "direct air capture": "coral", "co2 storage": "black", "methanol storage": "cadetblue", "electricity": "grey", "hydrogen": "grey", "co2": "grey", "methanol": "grey", } n.add("Carrier", carriers.keys(), color=carriers.values())

Demand¶

Add a constant methanol demand adding up to an annual target production of 1 TWh of methanol:

In [10]:

n.add(
    "Load",
    "demand",
    bus="methanol",
    p_set=1e6 / 8760,
)

n.add( "Load", "demand", bus="methanol", p_set=1e6 / 8760, )

Wind and solar¶

Add the wind and solar generators:

In [11]:

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

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

In [12]:

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

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

Batteries¶

Add a 3-hour and a 6-hour battery storage:

In [13]:

for max_hours in [3, 6]:
    n.add(
        "StorageUnit",
        f"battery storage {max_hours}h",
        bus="electricity",
        carrier=f"battery storage {max_hours}h",
        max_hours=max_hours,
        capital_cost=costs.at["battery inverter", "capital_cost"]
        + max_hours * 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,
    )

for max_hours in [3, 6]: n.add( "StorageUnit", f"battery storage {max_hours}h", bus="electricity", carrier=f"battery storage {max_hours}h", max_hours=max_hours, capital_cost=costs.at["battery inverter", "capital_cost"] + max_hours * 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, )

Hydrogen¶

Add electrolysers, hydrogen storage (steel tank), hydrogen turbine:

In [14]:

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 [15]:

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 [16]:

tech = "hydrogen storage tank type 1 including compressor"

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

tech = "hydrogen storage tank type 1 including compressor" n.add( "Store", "hydrogen storage", bus="hydrogen", carrier="hydrogen storage", capital_cost=costs.at[tech, "capital_cost"], e_nom_extendable=True, e_cyclic=True, )

Carbon Dioxide¶

Add liquid carbon dioxide storage and direct air capture, assuming for simplicity an electricity demand 2 MWh/tCO2 for heat and electricity needs of the process. More detailed modelling would also model the heat supply for the direct air capture.

In [17]:

electricity_input = 2  # MWh/tCO2

n.add(
    "Link",
    "direct air capture",
    bus0="electricity",
    bus1="co2",
    carrier="direct air capture",
    p_nom_extendable=True,
    efficiency=1 / electricity_input,
    capital_cost=costs.at["direct air capture", "capital_cost"] / electricity_input,
)

electricity_input = 2 # MWh/tCO2 n.add( "Link", "direct air capture", bus0="electricity", bus1="co2", carrier="direct air capture", p_nom_extendable=True, efficiency=1 / electricity_input, capital_cost=costs.at["direct air capture", "capital_cost"] / electricity_input, )

In [18]:

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

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

Methanol¶

Add methanolisation unit, which takes hydrogen, electricity and carbon dioxide as input, and methanol storage. Efficiencies and capital costs need to be expressed in units of bus0.

In [19]:

eff_h2 = 1 / costs.at["methanolisation", "hydrogen-input"]

n.add(
    "Link",
    "methanolisation",
    bus0="hydrogen",
    bus1="methanol",
    bus2="electricity",
    bus3="co2",
    carrier="methanolisation",
    p_nom_extendable=True,
    capital_cost=costs.at["methanolisation", "capital_cost"] * eff_h2,
    efficiency=eff_h2,
    efficiency2=-costs.at["methanolisation", "electricity-input"] * eff_h2,
    efficiency3=-costs.at["methanolisation", "carbondioxide-input"] * eff_h2,
)

eff_h2 = 1 / costs.at["methanolisation", "hydrogen-input"] n.add( "Link", "methanolisation", bus0="hydrogen", bus1="methanol", bus2="electricity", bus3="co2", carrier="methanolisation", p_nom_extendable=True, capital_cost=costs.at["methanolisation", "capital_cost"] * eff_h2, efficiency=eff_h2, efficiency2=-costs.at["methanolisation", "electricity-input"] * eff_h2, efficiency3=-costs.at["methanolisation", "carbondioxide-input"] * eff_h2, )

Costs for g is given in €/m³. We need to convert it to €/MWh, using the volumetric energy density of methanol (15.6 MJ/L).

In [20]:

capital_cost = costs.at[
    "General liquid hydrocarbon storage (crude)", "capital_cost"
] / (15.6 * 1000 / 3600)

n.add(
    "Store",
    "methanol storage",
    bus="co2",
    carrier="methanol storage",
    capital_cost=capital_cost,
    e_nom_extendable=True,
    e_cyclic=True,
)

capital_cost = costs.at[ "General liquid hydrocarbon storage (crude)", "capital_cost" ] / (15.6 * 1000 / 3600) n.add( "Store", "methanol storage", bus="co2", carrier="methanol storage", capital_cost=capital_cost, e_nom_extendable=True, e_cyclic=True, )

Optimisation¶

In [21]:

n.optimize()

n.optimize()

/tmp/ipykernel_4999/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()
WARNING:pypsa.consistency:Encountered nan's in static data for columns ['efficiency3', 'efficiency2'] of component 'Link'.

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/19 [00:00<?, ?it/s]

Writing constraints.:  89%|████████▉ | 17/19 [00:00<00:00, 159.84it/s]

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

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

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

INFO:linopy.io: Writing time: 0.17s

INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 39323 primals, 85187 duals
Objective: 2.24e+08
Solver model: available
Solver message: Optimal

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints 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[21]:

('ok', 'optimal')

Model evaluation¶

Total system cost by technology (we only added CAPEX components):

In [22]:

n.statistics.capex().div(1e6).sort_values(ascending=False)  # mn€/a

n.statistics.capex().div(1e6).sort_values(ascending=False) # mn€/a

Out[22]:

component    carrier           
Generator    wind                  65.960159
             solar                 53.567253
Store        hydrogen storage      48.579899
Link         direct air capture    33.507477
             methanolisation       16.916982
StorageUnit  battery storage 6h     5.003787
Store        methanol storage       0.077775
Link         electrolysis           0.031580
dtype: float64

Costs per unit of fuel (€/MWh):

In [23]:

n.statistics.capex().sum() / (8760 * n.loads.p_set.sum())

n.statistics.capex().sum() / (8760 * n.loads.p_set.sum())

Out[23]:

np.float64(223.64491163334)

The optimised capacities in MW (MWh for Store component):

In [24]:

n.statistics.optimal_capacity()

n.statistics.optimal_capacity()

Out[24]:

component    carrier           
Generator    solar                   829.72636
             wind                    514.08343
Link         direct air capture       61.73482
             electrolysis            502.02806
             methanolisation         129.90868
StorageUnit  battery storage 6h       38.80922
Store        hydrogen storage       8820.73540
             methanol storage      13477.63523
dtype: float64

Utilisation rates and capacity factors for each technology (in percent):

In [25]:

n.statistics.capacity_factor() * 100

n.statistics.capacity_factor() * 100

Out[25]:

component    carrier           
Generator    solar                   9.451986
             wind                   42.436697
Link         direct air capture     91.716474
             electrolysis           41.622610
             methanolisation       100.000000
Load         -                            NaN
StorageUnit  battery storage 6h      4.509109
Store        hydrogen storage       19.240184
             methanol storage       52.682847
dtype: float64

Curtailment for each technology (in TWh):

In [26]:

n.statistics.curtailment().div(1e6)

n.statistics.curtailment().div(1e6)

Out[26]:

component    carrier           
Generator    solar                 0.213441
             wind                  0.665860
Link         direct air capture    0.044674
             electrolysis          2.560267
StorageUnit  battery storage 6h    0.339661
dtype: float64

System operation on electricity, hydrogen, carbon dioxide and methanol sides (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 780.25x300 with 1 Axes>,
 <Axes: xlabel='snapshot', ylabel='Energy Balance [MW]'>,
 <seaborn.axisgrid.FacetGrid at 0x730b18652a50>)

In [28]:

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

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

Out[28]:

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

In [29]:

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

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

Out[29]:

(<Figure size 773.25x300 with 1 Axes>,
 <Axes: xlabel='snapshot', ylabel='Energy Balance [t/h]'>,
 <seaborn.axisgrid.FacetGrid at 0x730af4dfbc50>)

In [30]:

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

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

Out[30]:

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

Modifications: Biogenic Carbon Dioxide¶

How are costs affected by the availability of biogenic carbon dioxide costed at 50 €/t?

In [31]:

n.add(
    "Generator",
    "biogenic co2",
    bus="co2",
    carrier="biogenic co2",
    p_nom=1000,  # non-binding
    marginal_cost=50,
)
n.add(
    "Carrier",
    "biogenic co2",
    color="forestgreen",
)
n.optimize()

n.add( "Generator", "biogenic co2", bus="co2", carrier="biogenic co2", p_nom=1000, # non-binding marginal_cost=50, ) n.add( "Carrier", "biogenic co2", color="forestgreen", ) n.optimize()

/tmp/ipykernel_4999/4026235608.py:14: 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()
WARNING:pypsa.consistency:Encountered nan's in static data for columns ['efficiency3', 'efficiency2'] of component 'Link'.

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, 174.11it/s]

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

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

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

INFO:linopy.io: Writing time: 0.17s

INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 41507 primals, 89555 duals
Objective: 2.01e+08
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[31]:

('ok', 'optimal')

In [32]:

(n.statistics.capex().sum() + n.statistics.opex().sum()) / 1e6  # €/MWh

(n.statistics.capex().sum() + n.statistics.opex().sum()) / 1e6 # €/MWh

Out[32]:

np.float64(188.15502293155)

In [33]:

n.statistics.energy_balance(bus_carrier="co2")

n.statistics.energy_balance(bus_carrier="co2")

Out[33]:

component  carrier          bus_carrier
Link       methanolisation  co2           -247320.54795
dtype: float64

Modifications: Inflexible PtX¶

How are costs affected by limited flexibility of electrolyer and methanolisation units (e.g. with a minimum part load of 80%)?

In [34]:

n.links.loc["electrolysis", "p_min_pu"] = 0.8
n.links.loc["methanolisation", "p_min_pu"] = 0.8
n.optimize()

n.links.loc["electrolysis", "p_min_pu"] = 0.8 n.links.loc["methanolisation", "p_min_pu"] = 0.8 n.optimize()

/tmp/ipykernel_4999/1351952013.py:3: 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()

WARNING:pypsa.consistency:Encountered nan's in static data for columns ['efficiency3', 'efficiency2'] of component 'Link'.

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, 177.60it/s]

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

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

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

INFO:linopy.io: Writing time: 0.16s

INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 41507 primals, 89555 duals
Objective: 3.03e+08
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[34]:

('ok', 'optimal')

In [35]:

(n.statistics.capex().sum() + n.statistics.opex().sum()) / 1e6  # €/MWh

(n.statistics.capex().sum() + n.statistics.opex().sum()) / 1e6 # €/MWh

Out[35]:

np.float64(303.27748910232)

In [36]:

n.statistics.optimal_capacity()

n.statistics.optimal_capacity()

Out[36]:

component    carrier           
Generator    biogenic co2          1000.00000
             solar                  780.11055
             wind                  1010.92272
Link         electrolysis           211.83481
             methanolisation        129.90868
             turbine                 15.38233
StorageUnit  battery storage 6h     607.51087
Store        hydrogen storage      2418.75950
dtype: float64