Advanced models in pypsa: capacity expansion planning with sector coupling

Advanced models in `pypsa`: capacity expansion planning with sector coupling#

Note

If you have not yet set up Python on your computer, you can execute this tutorial in your browser via Google Colab. Click on the rocket in the top right corner and launch “Colab”. If that doesn’t work download the .ipynb file and import it in Google Colab.

Then install the following packages by executing the following command in a Jupyter cell at the top of the notebook.

!pip install pypsa pandas numpy matplotlib plotly highspy

Problem description#

To explore capacity expansion problem with sector-coupling options, let’s model an “greenfield” energy system model with:

electricity demand (here: Germany 2015)
several technologies for electricity generation (wind, solar, and gas for peak load)
hydrogen consumption (e.g., an offtaker in the industrial sector)
hydrogen production from electrolysis
hydrogen storage
hydrogen fuel cell

import pypsa
import pandas as pd

pd.options.plotting.backend = "plotly"

Prepare technology data#

At TU Berlin, we maintain a database (PyPSA/technology-data) which collects assumptions and projections for energy system technologies (such as costs, efficiencies, lifetimes, etc.) for given years, which we use for our research.

Reading this data into a useable pandas.DataFrame requires some pre-processing (e.g. converting units, setting defaults, re-arranging dimensions):

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.unit = costs.unit.str.replace("/kW", "/MW")

defaults = {
    "FOM": 0,
    "VOM": 0,
    "efficiency": 1,
    "fuel": 0,
    "investment": 0,
    "lifetime": 25,
    "CO2 intensity": 0,
    "discount rate": 0.07,
}
costs = costs.value.unstack().fillna(defaults)

costs.at["OCGT", "fuel"] = costs.at["gas", "fuel"]
# costs.at["CCGT", "fuel"] = costs.at["gas", "fuel"]
costs.at["OCGT", "CO2 intensity"] = costs.at["gas", "CO2 intensity"]
# costs.at["CCGT", "CO2 intensity"] = costs.at["gas", "CO2 intensity"]

# costs

Let’s also write a small utility function that calculates the annuity to annualise investment costs. The formula is

\[ a(r, n) = \frac{r}{1-(1+r)^{-n}} \]

where \(r\) is the discount rate and \(n\) is the lifetime.

def annuity_factor(r, n):
    return r / (1.0 - 1.0 / (1.0 + r) ** n)

The resulting annuity_factor (often called capital recovery factor) represents the factor that, when multiplied by the initial capital cost, yields the equivalent annual cost accounting for the time value of money:

annuity_factor(0.07, 20)

0.09439292574325567

Based on this, we can calculate the short-term marginal generation costs (€/MWh)

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

and the annualised investment costs (capital_cost in PyPSA terms, €/MW/a):

annuity = costs.apply(lambda x: annuity_factor(0.07, x["lifetime"]), axis=1)

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

Prepare load and renewable generation time-series#

We are also going to need some time series for wind, solar and load.

For now, we are going to recycle the time series we used in the pandas tutorial.

url = (
    "https://tubcloud.tu-berlin.de/s/pKttFadrbTKSJKF/download/time-series-lecture-2.csv"
)
ts = pd.read_csv(url, index_col=0, parse_dates=True)

# fallback if url is not available
# ts = pd.read_csv("../data/time-series-lecture-2.csv", index_col=0, parse_dates=True)

ts.head()

	load	onwind	offwind	prices
2015-01-01 00:00:00	41.151	0.1566	0.7030	NaN
2015-01-01 01:00:00	40.135	0.1659	0.6875	NaN
2015-01-01 02:00:00	39.106	0.1746	0.6535	NaN
2015-01-01 03:00:00	38.765	0.1745	0.6803	NaN
2015-01-01 04:00:00	38.941	0.1826	0.7272	NaN

Let’s also convert the load time series from GW to MW, the base unit of PyPSA:

ts.load *= 1e3

Optionally, we can downscale temporal resolution of the time series to save some computation time:

# here we sample only every other hour:
resolution = 4
ts = ts.resample(f"{resolution}h").first()

Build our energy system to simulate#

As always, let’s initialize an empty network:

n = pypsa.Network()

Then, we add a single electricity bus…

n.add("Bus", "electricity")

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

…and tell the pypsa.Network object n what the snapshots of the model will be using the utility function n.set_snapshots().

n.set_snapshots(ts.index)

n.snapshots

DatetimeIndex(['2015-01-01 00:00:00', '2015-01-01 04:00:00',
               '2015-01-01 08:00:00', '2015-01-01 12:00:00',
               '2015-01-01 16:00:00', '2015-01-01 20:00:00',
               '2015-01-02 00:00:00', '2015-01-02 04:00:00',
               '2015-01-02 08:00:00', '2015-01-02 12:00:00',
               ...
               '2015-12-30 08:00:00', '2015-12-30 12:00:00',
               '2015-12-30 16:00:00', '2015-12-30 20:00:00',
               '2015-12-31 00:00:00', '2015-12-31 04:00:00',
               '2015-12-31 08:00:00', '2015-12-31 12:00:00',
               '2015-12-31 16:00:00', '2015-12-31 20:00:00'],
              dtype='datetime64[ns]', name='snapshot', length=2190, freq='4h')

If we resampled the time series above, we need to adjust the weighting of the snapshots (i.e. how many hours they represent). We can do that with n.snapshot_weightings:

n.snapshot_weightings.head(3)

	objective	stores	generators
snapshot
2015-01-01 00:00:00	1.0	1.0	1.0
2015-01-01 04:00:00	1.0	1.0	1.0
2015-01-01 08:00:00	1.0	1.0	1.0

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

n.snapshot_weightings.head(3)

	objective	stores	generators
snapshot
2015-01-01 00:00:00	4.0	4.0	4.0
2015-01-01 04:00:00	4.0	4.0	4.0
2015-01-01 08:00:00	4.0	4.0	4.0

Adding components for electricity part#

Then, we add all the technologies we are going to include as carriers.

carriers = [
    "onwind",
    "offwind",
    "solar",
    "OCGT",
    "hydrogen storage underground",
]

n.add(
    "Carrier",
    carriers,
    color=["dodgerblue", "aquamarine", "gold", "indianred", "magenta"],
    co2_emissions=[costs.at[c, "CO2 intensity"] for c in carriers],
)

Index(['onwind', 'offwind', 'solar', 'OCGT', 'hydrogen storage underground'], dtype='object')

Next, we add the demand time series to the model.

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

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

Let’s have a check whether the data was read-in correctly.

n.loads_t.p_set.plot()

We are going to add one dispatchable generation technology to the model. This is an open-cycle gas turbine (OCGT) with CO\(_2\) emissions of 0.2 t/MWh\(_{th}\).

n.add(
    "Generator",
    "OCGT",
    bus="electricity",
    carrier="OCGT",
    capital_cost=costs.at["OCGT", "capital_cost"],
    marginal_cost=costs.at["OCGT", "marginal_cost"],
    efficiency=costs.at["OCGT", "efficiency"],
    p_nom_extendable=True,
)

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

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.

for tech in ["onwind", "offwind", "solar"]:
    n.add(
        "Generator",
        tech,
        bus="electricity",
        carrier=tech,
        p_max_pu=ts[tech],
        capital_cost=costs.at[tech, "capital_cost"],
        marginal_cost=costs.at[tech, "marginal_cost"],
        efficiency=costs.at[tech, "efficiency"],
        p_nom_extendable=True,
    )

# Making sure the capacity factors are read-in correctly
n.generators_t.p_max_pu.loc["2015-03"].plot()

Adding components for hydrogen part#

Add a dedicated Bus for the hydrogen energy carrier:

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

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

Add a Link for the hydrogen electrolysis:

n.add(
    "Link",
    "electrolysis",
    bus0="electricity",
    bus1="hydrogen",
    carrier="electrolysis",
    p_nom_extendable=True,
    efficiency=0.7,
    capital_cost=50e3,  # €/MW/a
)

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

Note

Some of the sector-coupling technologies we are going to add have multiple ouputs (e.g. CHP plants producing heat and power). PyPSA can automatically handle links have more than one input (bus0) and/or output (i.e. bus1, bus2, bus3) with a given efficieny (efficiency, efficiency2, efficiency3).

Add a Link for the fuel cell which reconverts hydrogen to electricity:

n.add(
    "Link",
    "fuel cell",
    bus0="hydrogen",
    bus1="electricity",
    carrier="fuel cell",
    p_nom_extendable=True,
    efficiency=0.5,
    capital_cost=120e3,  # €/MW/a
)

Index(['fuel cell'], dtype='object')

Add a Store for the hydrogen storage:

n.add(
    "Store",
    "hydrogen storage",
    bus="hydrogen",
    carrier="hydrogen storage",
    capital_cost=140,  # €/MWh/a
    e_nom_extendable=True,
    e_cyclic=True,  # cyclic state of charge
)

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

To model an industrial hydrogen offtaker, we add also a hydrogen demand to the hydrogen bus.

In the example below, we add a hydrogen demand such that it equals ~25% of the electricity demand (in MWh_H2)

n.add(
    "Load", "hydrogen demand", bus="hydrogen", carrier="hydrogen", p_set=19500
)  # MWh_H2/h

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

We are now ready to solve the model#

# n.optimize.create_model()

n.optimize(solver_name="highs")
# 66 seconds for 1H temporal resolution (8760 snapshots)
# 7 seconds for 4H temporal resolution (2190 snapshots)

WARNING:pypsa.consistency:The following loads have carriers which are not defined:
Index(['hydrogen demand'], dtype='object', name='Load')

WARNING:pypsa.consistency:The following links have carriers which are not defined:
Index(['electrolysis', 'fuel cell'], dtype='object', name='Link')

WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['electricity', 'hydrogen'], dtype='object', name='Bus')

WARNING:pypsa.consistency:The following stores have carriers which are not defined:
Index(['hydrogen storage'], dtype='object', name='Store')

WARNING:pypsa.consistency:The following loads have carriers which are not defined:
Index(['hydrogen demand'], dtype='object', name='Load')

WARNING:pypsa.consistency:The following links have carriers which are not defined:
Index(['electrolysis', 'fuel cell'], dtype='object', name='Link')

WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['electricity', 'hydrogen'], dtype='object', name='Bus')

WARNING:pypsa.consistency:The following stores have carriers which are not defined:
Index(['hydrogen storage'], dtype='object', name='Store')

INFO:linopy.model: Solve problem using Highs solver

INFO:linopy.io:Writing objective.

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

Writing constraints.:  71%|███████▏  | 10/14 [00:00<00:00, 82.73it/s]

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

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

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

INFO:linopy.io: Writing time: 0.26s

Running HiGHS 1.9.0 (git hash: fa40bdf): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e-04, 4e+00]
  Cost   [4e-02, 2e+05]
  Bound  [0e+00, 0e+00]
  RHS    [2e+04, 8e+04]
Presolving model
20850 rows, 16477 cols, 53790 nonzeros  0s
18660 rows, 14287 cols, 49410 nonzeros  0s
Presolve : Reductions: rows 18660(-18577); columns 14287(-3240); elements 49410(-21817)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     0.0000000000e+00 Pr: 4380(4.17793e+09) 0s

INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 17527 primals, 37237 duals
Objective: 4.75e+10
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, Store-energy_balance were not assigned to the network.

      18025     4.7521933007e+10 Pr: 0(0); Du: 0(3.95961e-13) 4s
Solving the original LP from the solution after postsolve
Model name          : linopy-problem-c24n0crh
Model status        : Optimal
Simplex   iterations: 18025
Objective value     :  4.7521933007e+10
Relative P-D gap    :  1.1639490831e-13
HiGHS run time      :          4.32
Writing the solution to /tmp/linopy-solve-m9hffbk1.sol

('ok', 'optimal')

Exploring model results#

The total system cost in billion Euros per year:

n.objective / 1e9

47.521933007318715

n.statistics() provides an informative overview of the model results.

See documentation: https://pypsa.readthedocs.io/en/stable/api/statistics.html

n.statistics()

		Optimal Capacity	Supply	Withdrawal	Energy Balance	Capacity Factor	Curtailment	Capital Expenditure	Operational Expenditure	Revenue	Market Value
Generator	OCGT	69495.97520	2.642126e+08	0.000000e+00	2.642126e+08	0.434000	3.445721e+08	3.316256e+09	1.709032e+10	2.040657e+10	77.235416
	offwind	97743.93570	2.985578e+08	0.000000e+00	2.985578e+08	0.348686	1.207735e+07	1.706181e+10	6.329425e+06	1.706814e+10	57.168634
	solar	147591.95957	1.601839e+08	0.000000e+00	1.601839e+08	0.123894	3.767638e+05	7.578379e+09	1.697949e+06	7.580077e+09	47.321103
Link	electrolysis	46714.11681	1.708200e+08	2.440286e+08	-7.320857e+07	0.596332	0.000000e+00	2.335706e+09	0.000000e+00	2.335706e+09	13.673492
Load	-	0.00000	0.000000e+00	4.789257e+08	-4.789257e+08	NaN	0.000000e+00	0.000000e+00	0.000000e+00	-3.142989e+10	NaN
Load	hydrogen	0.00000	0.000000e+00	1.708200e+08	-1.708200e+08	NaN	0.000000e+00	0.000000e+00	0.000000e+00	-1.609205e+10	NaN
Store	hydrogen storage	938822.02685	3.675079e+07	3.675079e+07	0.000000e+00	0.088991	0.000000e+00	1.314351e+08	0.000000e+00	1.314351e+08	3.576388

We can use n.statisitcs() to get a quick overview of optimised capacities across all components:

n.statistics.expanded_capacity().div(1e3).round(1)  # GW

component  carrier         
Link       electrolysis         46.7
Generator  OCGT                 69.5
           offwind              97.7
           solar               147.6
Store      hydrogen storage    938.8
dtype: float64

You can also use n.statistics() to promptly get an energy balance for the complete system or even any specific bus:

n.statistics.energy_balance(aggregate_time=False)

		snapshot	2015-01-01 00:00:00	2015-01-01 04:00:00	2015-01-01 08:00:00	2015-01-01 12:00:00	2015-01-01 16:00:00	2015-01-01 20:00:00	2015-01-02 00:00:00	2015-01-02 04:00:00	2015-01-02 08:00:00	2015-01-02 12:00:00	...	2015-12-30 08:00:00	2015-12-30 12:00:00	2015-12-30 16:00:00	2015-12-30 20:00:00	2015-12-31 00:00:00	2015-12-31 04:00:00	2015-12-31 08:00:00	2015-12-31 12:00:00	2015-12-31 16:00:00	2015-12-31 20:00:00
component	carrier	bus_carrier
Link	electrolysis	AC	-27562.98680	-32138.39004	-40074.51748	-46714.11681	-34774.57879	-46345.52625	-46714.11681	-46714.11681	-42621.69190	-46714.11681	...	-26843.95389	-46714.11681	-25168.36394	-33945.88902	-42405.86367	-39956.80787	-30172.89347	-46714.11681	-3053.24524	NaN
	fuel cell	hydrogen	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
	electrolysis	hydrogen	19294.09076	22496.87303	28052.16223	32699.88177	24342.20515	32441.86838	32699.88177	32699.88177	29835.18433	32699.88177	...	18790.76772	32699.88177	17617.85476	23762.12231	29684.10457	27969.76551	21121.02543	32699.88177	2137.27167	NaN
	fuel cell	AC	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
Generator	OCGT	AC	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	6588.61467
	offwind	AC	68713.98680	71079.39004	76513.95286	63267.54439	88184.57879	95075.52625	87291.11681	89096.11681	96219.13030	88983.96447	...	73268.85420	71457.24329	83991.36394	86004.88902	83189.86367	80462.80787	74256.06795	75183.49011	57698.24524	39322.38533
	onwind	AC	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	...	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
	solar	AC	NaN	NaN	6523.56461	30610.57241	NaN	NaN	NaN	NaN	2907.56160	16294.15234	...	6981.09969	29488.87352	NaN	NaN	NaN	NaN	6139.82552	22271.62670	NaN	NaN
Store	hydrogen storage	hydrogen	205.90924	-2996.87303	-8552.16223	-13199.88177	-4842.20515	-12941.86838	-13199.88177	-13199.88177	-10335.18433	-13199.88177	...	709.23228	-13199.88177	1882.14524	-4262.12231	-10184.10457	-8469.76551	-1621.02543	-13199.88177	17362.72833	19500.00000
Load	-	AC	-41151.00000	-38941.00000	-42963.00000	-47164.00000	-53410.00000	-48730.00000	-40577.00000	-42382.00000	-56505.00000	-58564.00000	...	-53406.00000	-54232.00000	-58823.00000	-52059.00000	-40784.00000	-40506.00000	-50223.00000	-50741.00000	-54645.00000	-45911.00000
Load	hydrogen	hydrogen	-19500.00000	-19500.00000	-19500.00000	-19500.00000	-19500.00000	-19500.00000	-19500.00000	-19500.00000	-19500.00000	-19500.00000	...	-19500.00000	-19500.00000	-19500.00000	-19500.00000	-19500.00000	-19500.00000	-19500.00000	-19500.00000	-19500.00000	-19500.00000

11 rows × 2190 columns

n.statistics.energy_balance(aggregate_time=False, bus_carrier="AC").div(1e3).groupby(
    "carrier"
).sum().T.plot()

n.statistics.energy_balance(aggregate_time=False, bus_carrier="hydrogen").div(
    1e3
).groupby("carrier").sum().T.plot()

Possibly, we are also interested in the total emissions:

emissions = (
    n.generators_t.p
    / n.generators.efficiency
    * n.generators.carrier.map(n.carriers.co2_emissions)
)  # t/h

n.snapshot_weightings.generators @ emissions.sum(axis=1).div(1e6)  # Mt

127.59537376038045

Adding emission limits#

The gas power plant offers sufficient and cheap enough backup capacity to run in periods of low wind and solar generation. But what happens if this source of flexibility disappears. Let’s model a 100% renewable electricity system by adding a CO\(_2\) emission limit as global constraint:

n.add(
    "GlobalConstraint",
    "CO2Limit",
    carrier_attribute="co2_emissions",
    sense="<=",
    constant=0,
)

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

When we run the model now…

n.optimize(solver_name="highs")

Show code cell output Hide code cell output

WARNING:pypsa.consistency:The following loads have carriers which are not defined:
Index(['hydrogen demand'], dtype='object', name='Load')

WARNING:pypsa.consistency:The following links have carriers which are not defined:
Index(['electrolysis', 'fuel cell'], dtype='object', name='Link')

WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['electricity', 'hydrogen'], dtype='object', name='Bus')

WARNING:pypsa.consistency:The following stores have carriers which are not defined:
Index(['hydrogen storage'], dtype='object', name='Store')

WARNING:pypsa.consistency:The following loads have carriers which are not defined:
Index(['hydrogen demand'], dtype='object', name='Load')

WARNING:pypsa.consistency:The following links have carriers which are not defined:
Index(['electrolysis', 'fuel cell'], dtype='object', name='Link')

WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['electricity', 'hydrogen'], dtype='object', name='Bus')

WARNING:pypsa.consistency:The following stores have carriers which are not defined:
Index(['hydrogen storage'], dtype='object', name='Store')

INFO:linopy.model: Solve problem using Highs solver

INFO:linopy.io:Writing objective.

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

Writing constraints.:  67%|██████▋   | 10/15 [00:00<00:00, 86.54it/s]

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

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

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

INFO:linopy.io: Writing time: 0.25s

Running HiGHS 1.9.0 (git hash: fa40bdf): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e-04, 4e+00]
  Cost   [4e-02, 2e+05]
  Bound  [0e+00, 0e+00]
  RHS    [2e+04, 8e+04]
Presolving model
18660 rows, 14286 cols, 47220 nonzeros  0s
16470 rows, 12096 cols, 42840 nonzeros  0s
Presolve : Reductions: rows 16470(-20768); columns 12096(-5431); elements 42840(-30577)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     0.0000000000e+00 Pr: 4380(3.74436e+09) 0s

INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 17527 primals, 37238 duals
Objective: 7.53e+10
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, Store-energy_balance were not assigned to the network.

      11188     7.5338336224e+10 Pr: 0(0); Du: 0(1.20154e-13) 4s
Solving the original LP from the solution after postsolve
Model name          : linopy-problem-_ja4mbm3
Model status        : Optimal
Simplex   iterations: 11188
Objective value     :  7.5338336224e+10
Relative P-D gap    :  1.4380116123e-14
HiGHS run time      :          3.77
Writing the solution to /tmp/linopy-solve-afczry3b.sol

('ok', 'optimal')

The total system cost (billion Euros per year) is now higher than before:

n.objective / 1e9

75.33833622431514

The optimal capacity mix now does not include any gas power plants and includes onshore wind and fuel cell technologies:

n.statistics.expanded_capacity().div(1e3).round(1)  # GW

component  carrier         
Link       electrolysis          145.0
           fuel cell             135.2
Generator  offwind               101.8
           onwind                111.9
           solar                 332.5
Store      hydrogen storage    38417.0
dtype: float64

Fuel cell technology steps in hours with low wind and solar generation:

n.statistics.energy_balance(aggregate_time=False, bus_carrier="AC").div(1e3).groupby(
    "carrier"
).sum().T.plot()

n.statistics.energy_balance(aggregate_time=False, bus_carrier="hydrogen").div(
    1e3
).groupby("carrier").sum().T.plot()

n.stores_t.e.plot()

Total emissions are now zero:

emissions = (
    n.generators_t.p
    / n.generators.efficiency
    * n.generators.carrier.map(n.carriers.co2_emissions)
)  # t/h

n.snapshot_weightings.generators @ emissions.sum(axis=1).div(1e6)  # Mt

0.0