Note

You can download this example as a Jupyter notebook or start it in interactive mode.

Islanded methanol 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.

[1]:

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:

[2]:

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\).

[3]:

def annuity(r: float, n: int) -> float:
    return r / (1.0 - 1.0 / (1.0 + r) ** n)

[4]:

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

[5]:

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

Wind and solar time series#

Go to model.energy to retrieve wind and solar capacity factor time series for other countries.

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

[7]:

ts.head(3)

[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.

[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

[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())

[9]:

Index(['wind', 'solar', 'hydrogen storage', 'battery storage 3h',
       'battery storage 6h', 'electrolysis', 'turbine', 'methanolisation',
       'direct air capture', 'co2 storage', 'methanol storage', 'electricity',
       'hydrogen', 'co2', 'methanol'],
      dtype='object')

Demand#

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

[10]:

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

[10]:

Index(['demand'], dtype='object')

Wind and solar#

Add the wind and solar generators:

[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,
)

[11]:

Index(['wind'], dtype='object')

[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,
)

[12]:

Index(['solar'], dtype='object')

Batteries#

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

[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,
    )

Hydrogen#

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

[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"],
)

[14]:

Index(['electrolysis'], dtype='object')

[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"],
)

[15]:

Index(['turbine'], dtype='object')

[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,
)

[16]:

Index(['hydrogen storage'], dtype='object')

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.

[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,
)

[17]:

Index(['direct air capture'], dtype='object')

[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,
)

[18]:

Index(['co2 storage'], dtype='object')

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.

[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,
)

[19]:

Index(['methanolisation'], dtype='object')

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)

[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,
)

[20]:

Index(['methanol storage'], dtype='object')

Optimisation#

[21]:

n.optimize(solver_name="highs")

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.io:Writing objective.
Writing constraints.: 100%|██████████| 23/23 [00:00<00:00, 59.38it/s]
Writing continuous variables.: 100%|██████████| 11/11 [00:00<00:00, 163.23it/s]
INFO:linopy.io: Writing time: 0.46s

Running HiGHS 1.10.0 (git hash: fd86653): Copyright (c) 2025 HiGHS under MIT licence terms
LP   linopy-problem-1q75vek9 has 85187 rows; 39323 cols; 175813 nonzeros
Coefficient ranges:
  Matrix [1e-03, 6e+00]
  Cost   [5e+00, 4e+05]
  Bound  [0e+00, 0e+00]
  RHS    [1e+02, 1e+02]
Presolving model
49130 rows, 36037 cols, 129918 nonzeros  0s
40394 rows, 29484 cols, 114630 nonzeros  0s
Dependent equations search running on 10920 equations with time limit of 1000.00s
Dependent equations search removed 0 rows and 0 nonzeros in 0.00s (limit = 1000.00s)
40394 rows, 29484 cols, 114630 nonzeros  0s
Presolve : Reductions: rows 40394(-44793); columns 29484(-9839); elements 114630(-61183)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     9.9929899048e+06 Pr: 6552(493028) 0s
      33564     1.6960257845e+08 Pr: 1536(797424); Du: 0(8.63386e-07) 5s

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

      39871     1.7400225927e+08 Pr: 0(0); Du: 0(2.02779e-12) 7s
Solving the original LP from the solution after postsolve
Model name          : linopy-problem-1q75vek9
Model status        : Optimal
Simplex   iterations: 39871
Objective value     :  1.7400225927e+08
Relative P-D gap    :  1.1304182405e-14
HiGHS run time      :          6.77
Writing the solution to /tmp/linopy-solve-0dkyz_h5.sol

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.

[21]:

('ok', 'optimal')

Model evaluation#

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

[22]:

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

[22]:

component    carrier
Generator    wind                  52.253560
             solar                 42.603818
Store        hydrogen storage      38.449031
Link         direct air capture    26.649610
             methanolisation        9.992990
StorageUnit  battery storage 6h     3.960301
Store        methanol storage       0.061370
Link         electrolysis           0.031580
dtype: float64

Costs per unit of fuel (€/MWh):

[23]:

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

[23]:

np.float64(174.00225926853003)

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

[24]:

n.statistics.optimal_capacity()

[24]:

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

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

[25]:

n.statistics.capacity_factor() * 100

[25]:

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

Curtailment for each technology (in TWh):

[26]:

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

[26]:

component    carrier
Generator    solar                 0.212877
             wind                  0.666425
StorageUnit  battery storage 6h    0.339661
dtype: float64

System operation on electricity, hydrogen, carbon dioxide and methanol sides (in MW):

[27]:

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

[27]:

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

../_images/examples_islanded-methanol-production_50_1.png

[28]:

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

[28]:

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

../_images/examples_islanded-methanol-production_51_1.png

[29]:

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

[29]:

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

../_images/examples_islanded-methanol-production_52_1.png

[30]:

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

[30]:

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

../_images/examples_islanded-methanol-production_53_1.png

Modifications: Biogenic Carbon Dioxide#

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

[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(solver_name="highs")

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.io:Writing objective.
Writing constraints.: 100%|██████████| 25/25 [00:00<00:00, 61.38it/s]
Writing continuous variables.: 100%|██████████| 11/11 [00:00<00:00, 152.93it/s]
INFO:linopy.io: Writing time: 0.49s

Running HiGHS 1.10.0 (git hash: fd86653): Copyright (c) 2025 HiGHS under MIT licence terms
LP   linopy-problem-8s2jqpzt has 89555 rows; 41507 cols; 182365 nonzeros
Coefficient ranges:
  Matrix [1e-03, 6e+00]
  Cost   [5e+00, 4e+05]
  Bound  [0e+00, 0e+00]
  RHS    [1e+02, 1e+03]
Presolving model
49130 rows, 38221 cols, 132102 nonzeros  0s
40394 rows, 31668 cols, 116814 nonzeros  0s
Dependent equations search running on 10920 equations with time limit of 1000.00s
Dependent equations search removed 0 rows and 0 nonzeros in 0.00s (limit = 1000.00s)
40394 rows, 31668 cols, 116814 nonzeros  0s
Presolve : Reductions: rows 40394(-49161); columns 31668(-9839); elements 116814(-65551)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     9.9929899048e+06 Pr: 6552(481227) 0s

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

      28538     1.5814488721e+08 Pr: 0(0); Du: 0(6.8408e-09) 3s
Solving the original LP from the solution after postsolve
Model name          : linopy-problem-8s2jqpzt
Model status        : Optimal
Simplex   iterations: 28538
Objective value     :  1.5814488721e+08
Relative P-D gap    :  1.2437665942e-14
HiGHS run time      :          3.28
Writing the solution to /tmp/linopy-solve-ukxonr5h.sol

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.

[31]:

('ok', 'optimal')

[32]:

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

[32]:

np.float64(158.14488720719996)

[33]:

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

[33]:

component  carrier          bus_carrier
Link       methanolisation  co2           -247320.54795
Generator  biogenic co2     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%)?

[34]:

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

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.io:Writing objective.
Writing constraints.: 100%|██████████| 25/25 [00:00<00:00, 60.98it/s]
Writing continuous variables.: 100%|██████████| 11/11 [00:00<00:00, 152.62it/s]
INFO:linopy.io: Writing time: 0.49s

Running HiGHS 1.10.0 (git hash: fd86653): Copyright (c) 2025 HiGHS under MIT licence terms
LP   linopy-problem-08z_b6_2 has 89555 rows; 41507 cols; 186733 nonzeros
Coefficient ranges:
  Matrix [1e-03, 6e+00]
  Cost   [5e+00, 4e+05]
  Bound  [0e+00, 0e+00]
  RHS    [1e+02, 1e+03]
Presolving model
53498 rows, 38221 cols, 138654 nonzeros  0s
40940 rows, 30030 cols, 127244 nonzeros  0s
Dependent equations search running on 9282 equations with time limit of 1000.00s
Dependent equations search removed 0 rows and 0 nonzeros in 0.00s (limit = 1000.00s)
40940 rows, 30030 cols, 127244 nonzeros  0s
Presolve : Reductions: rows 40940(-48615); columns 30030(-11477); elements 127244(-59489)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     9.9929899048e+06 Pr: 6552(792734); Du: 0(2.33813e-10) 0s

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

      19351     2.3951008302e+08 Pr: 0(0); Du: 0(3.31966e-11) 3s
Solving the original LP from the solution after postsolve
Model name          : linopy-problem-08z_b6_2
Model status        : Optimal
Simplex   iterations: 19351
Objective value     :  2.3951008302e+08
Relative P-D gap    :  6.0970869310e-15
HiGHS run time      :          3.55
Writing the solution to /tmp/linopy-solve-n4go8qv_.sol

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.

[34]:

('ok', 'optimal')

[35]:

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

[35]:

np.float64(239.51008301740998)

[36]:

n.statistics.optimal_capacity()

[36]:

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

[ ]:

[ ]:

Islanded methanol production

Contents

Islanded methanol production#

Techno-economic assumptions#

Wind and solar time series#

Build model#

Initialisation#

Demand#

Wind and solar#

Batteries#

Hydrogen#

Carbon Dioxide#

Methanol#

Optimisation#

Model evaluation#

Modifications: Biogenic Carbon Dioxide#

Modifications: Inflexible PtX#