SciGRID Network¶
In this example, the dispatch of generators is optimised using the linear optimisation, then a non-linear power flow is run on the resulting dispatch. This example covers the German power system (roughly in the mid 2010s) for a single day at hourly resolution.
Data sources¶
Grid: Based on SciGRID which is based on OpenStreetMap.
Load size and location: Distributed proportional to NUTS3 GDP and population. In doubt, load and generation is attached to the 220kV substation.
Load time series: From ENTSO-E hourly data, scaled up uniformly by factor 1.12 (a simplification of the methodology in Schumacher and Hirth (2015)).
Conventional power plant capacities and locations: Based on list of BNetzA (German Federal Network Agency).
Wind and solar capacities and locations: Based on EEG Stammdaten, EnergyMap, which represents capacities at the end of 2014.
Wind and solar time series: Generated with REatlas tool based on Andresen et al. (2015).
Warnings¶
The data behind this example is no longer supported or updated. The dataset is only intended to demonstrate the capabilities of PyPSA and is not suitable for research purposes. Have a look at PyPSA-Eur for a newer grid model with active support.
Known problems include:
Rough approximations have been made for missing grid data, e.g. 220kV-380kV transformers and connections between close sub-stations missing from OSM.
There appears to be some unexpected congestion in parts of the network, which may mean for example that the load attachment method is inaccurate, particularly in regions with a high density of substations.
Attaching power plants to the nearest high voltage substation may not reflect reality.
The borders and neighbouring countries are not represented.
Hydroelectric power stations are not modelled accurately.
The marginal costs of generators are illustrative, not accurate.
Only the first day of 2011 is in the example, which is not representative.
The ENTSO-E total load for Germany may not be scaled correctly.
Biomass from the EEG reference data (Stammdaten) are not read in.
Power plant start up costs, ramping limits, and minimum part loads are not considered.
Input Data¶
In [2]:
Copied!
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pypsa
n = pypsa.examples.scigrid_de()
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pypsa
n = pypsa.examples.scigrid_de()
INFO:pypsa.network.io:Imported network 'SciGrid-DE' has buses, carriers, generators, lines, loads, storage_units, transformers
Plot the distribution of the load and of generating tech
In [3]:
Copied!
fig, ax = plt.subplots(
1,
1,
subplot_kw={"projection": ccrs.EqualEarth()},
)
load_distribution = n.loads_t.p_set.loc[n.snapshots[0]].groupby(n.loads.bus).sum()
n.plot(bus_size=load_distribution / 30000, ax=ax, title="Load distribution");
fig, ax = plt.subplots(
1,
1,
subplot_kw={"projection": ccrs.EqualEarth()},
)
load_distribution = n.loads_t.p_set.loc[n.snapshots[0]].groupby(n.loads.bus).sum()
n.plot(bus_size=load_distribution / 30000, ax=ax, title="Load distribution");
In [4]:
Copied!
n.generators.groupby("carrier")["p_nom"].sum().round(1)
n.generators.groupby("carrier")["p_nom"].sum().round(1)
Out[4]:
carrier Brown Coal 20879.5 Gas 23913.1 Geothermal 31.7 Hard Coal 25312.6 Multiple 152.7 Nuclear 12068.0 Oil 2710.2 Other 3027.8 Run of River 3999.1 Solar 37041.5 Storage Hydro 1445.0 Waste 1645.9 Wind Offshore 2973.5 Wind Onshore 37339.9 Name: p_nom, dtype: float64
In [5]:
Copied!
n.storage_units.groupby("carrier")["p_nom"].sum().round(1)
n.storage_units.groupby("carrier")["p_nom"].sum().round(1)
Out[5]:
carrier Pumped Hydro 9179.5 Name: p_nom, dtype: float64
In [6]:
Copied!
techs = ["Gas", "Brown Coal", "Hard Coal", "Wind Offshore", "Wind Onshore", "Solar"]
n_graphs = len(techs)
n_cols = 3
if n_graphs % n_cols == 0:
n_rows = n_graphs // n_cols
else:
n_rows = n_graphs // n_cols + 1
fig, axes = plt.subplots(
nrows=n_rows, ncols=n_cols, subplot_kw={"projection": ccrs.EqualEarth()}
)
size = 6
fig.set_size_inches(size * n_cols, size * n_rows)
for i, tech in enumerate(techs):
i_row = i // n_cols
i_col = i % n_cols
ax = axes[i_row, i_col]
gens = n.generators[n.generators.carrier == tech]
gen_distribution = (
gens.groupby("bus").sum()["p_nom"].reindex(n.buses.index, fill_value=0)
)
n.plot(ax=ax, bus_size=gen_distribution / 20000)
ax.set_title(tech)
fig.tight_layout()
techs = ["Gas", "Brown Coal", "Hard Coal", "Wind Offshore", "Wind Onshore", "Solar"]
n_graphs = len(techs)
n_cols = 3
if n_graphs % n_cols == 0:
n_rows = n_graphs // n_cols
else:
n_rows = n_graphs // n_cols + 1
fig, axes = plt.subplots(
nrows=n_rows, ncols=n_cols, subplot_kw={"projection": ccrs.EqualEarth()}
)
size = 6
fig.set_size_inches(size * n_cols, size * n_rows)
for i, tech in enumerate(techs):
i_row = i // n_cols
i_col = i % n_cols
ax = axes[i_row, i_col]
gens = n.generators[n.generators.carrier == tech]
gen_distribution = (
gens.groupby("bus").sum()["p_nom"].reindex(n.buses.index, fill_value=0)
)
n.plot(ax=ax, bus_size=gen_distribution / 20000)
ax.set_title(tech)
fig.tight_layout()
Rolling-Horizon Optimisation¶
Run optimisation (linear optimal power flow) on the first day of 2011.
To approximate $N-1$ security and allow room for reactive power flows, we do not allow any line to be loaded above 70% of their thermal rating:
In [7]:
Copied!
contingency_factor = 0.7
n.lines.s_max_pu = contingency_factor
contingency_factor = 0.7
n.lines.s_max_pu = contingency_factor
There are some infeasibilities without small extensions
In [8]:
Copied!
n.lines.loc[["316", "527", "602"], "s_nom"] = 1715
n.lines.loc[["316", "527", "602"], "s_nom"] = 1715
We are performing the optimisation for one day, 4 snapshots at a time.
In [9]:
Copied!
n.optimize.optimize_with_rolling_horizon(horizon=4, overlap=0);
n.optimize.optimize_with_rolling_horizon(horizon=4, overlap=0);
INFO:pypsa.optimization.abstract:Optimizing network for snapshot horizon [2011-01-01 00:00:00:2011-01-01 03:00:00] (1/6).
/home/docs/checkouts/readthedocs.org/user_builds/pypsa/checkouts/latest/pypsa/optimization/abstract.py:545: 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.
status, condition = n.optimize(sns, **kwargs)
WARNING:pypsa.consistency:The following transformers have zero r, which could break the linear load flow:
Index(['2', '5', '10', '12', '13', '15', '18', '20', '22', '24', '26', '30',
'32', '37', '42', '46', '52', '56', '61', '68', '69', '74', '78', '86',
'87', '94', '95', '96', '99', '100', '104', '105', '106', '107', '117',
'120', '123', '124', '125', '128', '129', '138', '143', '156', '157',
'159', '160', '165', '184', '191', '195', '201', '220', '231', '232',
'233', '236', '247', '248', '250', '251', '252', '261', '263', '264',
'267', '272', '279', '281', '282', '292', '303', '307', '308', '312',
'315', '317', '322', '332', '334', '336', '338', '351', '353', '360',
'362', '382', '384', '385', '391', '403', '404', '413', '421', '450',
'458'],
dtype='object', name='name')
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.model:Solver options: - log_to_console: False
INFO:linopy.io: Writing time: 0.09s
INFO:linopy.constants: Optimization successful: Status: ok Termination condition: optimal Solution: 9940 primals, 23828 duals Objective: 1.45e+06 Solver model: available Solver message: Optimal
INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Line-fix-s-lower, Line-fix-s-upper, Transformer-fix-s-lower, Transformer-fix-s-upper, StorageUnit-fix-p_dispatch-lower, StorageUnit-fix-p_dispatch-upper, StorageUnit-fix-p_store-lower, StorageUnit-fix-p_store-upper, StorageUnit-fix-state_of_charge-lower, StorageUnit-fix-state_of_charge-upper, Kirchhoff-Voltage-Law, StorageUnit-energy_balance were not assigned to the network.
INFO:pypsa.optimization.abstract:Optimizing network for snapshot horizon [2011-01-01 04:00:00:2011-01-01 07:00:00] (2/6).
/home/docs/checkouts/readthedocs.org/user_builds/pypsa/checkouts/latest/pypsa/optimization/abstract.py:545: 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. status, condition = n.optimize(sns, **kwargs)
WARNING:pypsa.consistency:The following transformers have zero r, which could break the linear load flow:
Index(['2', '5', '10', '12', '13', '15', '18', '20', '22', '24', '26', '30',
'32', '37', '42', '46', '52', '56', '61', '68', '69', '74', '78', '86',
'87', '94', '95', '96', '99', '100', '104', '105', '106', '107', '117',
'120', '123', '124', '125', '128', '129', '138', '143', '156', '157',
'159', '160', '165', '184', '191', '195', '201', '220', '231', '232',
'233', '236', '247', '248', '250', '251', '252', '261', '263', '264',
'267', '272', '279', '281', '282', '292', '303', '307', '308', '312',
'315', '317', '322', '332', '334', '336', '338', '351', '353', '360',
'362', '382', '384', '385', '391', '403', '404', '413', '421', '450',
'458'],
dtype='object', name='name')
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.model:Solver options: - log_to_console: False
INFO:linopy.io: Writing time: 0.08s
INFO:linopy.constants: Optimization successful: Status: ok Termination condition: optimal Solution: 9940 primals, 23828 duals Objective: 8.74e+05 Solver model: available Solver message: Optimal
INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Line-fix-s-lower, Line-fix-s-upper, Transformer-fix-s-lower, Transformer-fix-s-upper, StorageUnit-fix-p_dispatch-lower, StorageUnit-fix-p_dispatch-upper, StorageUnit-fix-p_store-lower, StorageUnit-fix-p_store-upper, StorageUnit-fix-state_of_charge-lower, StorageUnit-fix-state_of_charge-upper, Kirchhoff-Voltage-Law, StorageUnit-energy_balance were not assigned to the network.
INFO:pypsa.optimization.abstract:Optimizing network for snapshot horizon [2011-01-01 08:00:00:2011-01-01 11:00:00] (3/6).
/home/docs/checkouts/readthedocs.org/user_builds/pypsa/checkouts/latest/pypsa/optimization/abstract.py:545: 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.
status, condition = n.optimize(sns, **kwargs)
WARNING:pypsa.consistency:The following transformers have zero r, which could break the linear load flow:
Index(['2', '5', '10', '12', '13', '15', '18', '20', '22', '24', '26', '30',
'32', '37', '42', '46', '52', '56', '61', '68', '69', '74', '78', '86',
'87', '94', '95', '96', '99', '100', '104', '105', '106', '107', '117',
'120', '123', '124', '125', '128', '129', '138', '143', '156', '157',
'159', '160', '165', '184', '191', '195', '201', '220', '231', '232',
'233', '236', '247', '248', '250', '251', '252', '261', '263', '264',
'267', '272', '279', '281', '282', '292', '303', '307', '308', '312',
'315', '317', '322', '332', '334', '336', '338', '351', '353', '360',
'362', '382', '384', '385', '391', '403', '404', '413', '421', '450',
'458'],
dtype='object', name='name')
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.model:Solver options: - log_to_console: False
INFO:linopy.io: Writing time: 0.08s
INFO:linopy.constants: Optimization successful: Status: ok Termination condition: optimal Solution: 9940 primals, 23828 duals Objective: 7.91e+05 Solver model: available Solver message: Optimal
INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Line-fix-s-lower, Line-fix-s-upper, Transformer-fix-s-lower, Transformer-fix-s-upper, StorageUnit-fix-p_dispatch-lower, StorageUnit-fix-p_dispatch-upper, StorageUnit-fix-p_store-lower, StorageUnit-fix-p_store-upper, StorageUnit-fix-state_of_charge-lower, StorageUnit-fix-state_of_charge-upper, Kirchhoff-Voltage-Law, StorageUnit-energy_balance were not assigned to the network.
INFO:pypsa.optimization.abstract:Optimizing network for snapshot horizon [2011-01-01 12:00:00:2011-01-01 15:00:00] (4/6).
/home/docs/checkouts/readthedocs.org/user_builds/pypsa/checkouts/latest/pypsa/optimization/abstract.py:545: 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.
status, condition = n.optimize(sns, **kwargs)
WARNING:pypsa.consistency:The following transformers have zero r, which could break the linear load flow:
Index(['2', '5', '10', '12', '13', '15', '18', '20', '22', '24', '26', '30',
'32', '37', '42', '46', '52', '56', '61', '68', '69', '74', '78', '86',
'87', '94', '95', '96', '99', '100', '104', '105', '106', '107', '117',
'120', '123', '124', '125', '128', '129', '138', '143', '156', '157',
'159', '160', '165', '184', '191', '195', '201', '220', '231', '232',
'233', '236', '247', '248', '250', '251', '252', '261', '263', '264',
'267', '272', '279', '281', '282', '292', '303', '307', '308', '312',
'315', '317', '322', '332', '334', '336', '338', '351', '353', '360',
'362', '382', '384', '385', '391', '403', '404', '413', '421', '450',
'458'],
dtype='object', name='name')
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.model:Solver options: - log_to_console: False
INFO:linopy.io: Writing time: 0.08s
INFO:linopy.constants: Optimization successful: Status: ok Termination condition: optimal Solution: 9940 primals, 23828 duals Objective: 1.46e+06 Solver model: available Solver message: Optimal
INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Line-fix-s-lower, Line-fix-s-upper, Transformer-fix-s-lower, Transformer-fix-s-upper, StorageUnit-fix-p_dispatch-lower, StorageUnit-fix-p_dispatch-upper, StorageUnit-fix-p_store-lower, StorageUnit-fix-p_store-upper, StorageUnit-fix-state_of_charge-lower, StorageUnit-fix-state_of_charge-upper, Kirchhoff-Voltage-Law, StorageUnit-energy_balance were not assigned to the network.
INFO:pypsa.optimization.abstract:Optimizing network for snapshot horizon [2011-01-01 16:00:00:2011-01-01 19:00:00] (5/6).
/home/docs/checkouts/readthedocs.org/user_builds/pypsa/checkouts/latest/pypsa/optimization/abstract.py:545: 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.
status, condition = n.optimize(sns, **kwargs)
WARNING:pypsa.consistency:The following transformers have zero r, which could break the linear load flow:
Index(['2', '5', '10', '12', '13', '15', '18', '20', '22', '24', '26', '30',
'32', '37', '42', '46', '52', '56', '61', '68', '69', '74', '78', '86',
'87', '94', '95', '96', '99', '100', '104', '105', '106', '107', '117',
'120', '123', '124', '125', '128', '129', '138', '143', '156', '157',
'159', '160', '165', '184', '191', '195', '201', '220', '231', '232',
'233', '236', '247', '248', '250', '251', '252', '261', '263', '264',
'267', '272', '279', '281', '282', '292', '303', '307', '308', '312',
'315', '317', '322', '332', '334', '336', '338', '351', '353', '360',
'362', '382', '384', '385', '391', '403', '404', '413', '421', '450',
'458'],
dtype='object', name='name')
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.model:Solver options: - log_to_console: False
INFO:linopy.io: Writing time: 0.08s
INFO:linopy.constants: Optimization successful: Status: ok Termination condition: optimal Solution: 9940 primals, 23828 duals Objective: 2.65e+06 Solver model: available Solver message: Optimal
INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Line-fix-s-lower, Line-fix-s-upper, Transformer-fix-s-lower, Transformer-fix-s-upper, StorageUnit-fix-p_dispatch-lower, StorageUnit-fix-p_dispatch-upper, StorageUnit-fix-p_store-lower, StorageUnit-fix-p_store-upper, StorageUnit-fix-state_of_charge-lower, StorageUnit-fix-state_of_charge-upper, Kirchhoff-Voltage-Law, StorageUnit-energy_balance were not assigned to the network.
INFO:pypsa.optimization.abstract:Optimizing network for snapshot horizon [2011-01-01 20:00:00:2011-01-01 23:00:00] (6/6).
/home/docs/checkouts/readthedocs.org/user_builds/pypsa/checkouts/latest/pypsa/optimization/abstract.py:545: 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.
status, condition = n.optimize(sns, **kwargs)
WARNING:pypsa.consistency:The following transformers have zero r, which could break the linear load flow:
Index(['2', '5', '10', '12', '13', '15', '18', '20', '22', '24', '26', '30',
'32', '37', '42', '46', '52', '56', '61', '68', '69', '74', '78', '86',
'87', '94', '95', '96', '99', '100', '104', '105', '106', '107', '117',
'120', '123', '124', '125', '128', '129', '138', '143', '156', '157',
'159', '160', '165', '184', '191', '195', '201', '220', '231', '232',
'233', '236', '247', '248', '250', '251', '252', '261', '263', '264',
'267', '272', '279', '281', '282', '292', '303', '307', '308', '312',
'315', '317', '322', '332', '334', '336', '338', '351', '353', '360',
'362', '382', '384', '385', '391', '403', '404', '413', '421', '450',
'458'],
dtype='object', name='name')
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.model:Solver options: - log_to_console: False
INFO:linopy.io: Writing time: 0.08s
INFO:linopy.constants: Optimization successful: Status: ok Termination condition: optimal Solution: 9940 primals, 23828 duals Objective: 2.14e+06 Solver model: available Solver message: Optimal
INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Line-fix-s-lower, Line-fix-s-upper, Transformer-fix-s-lower, Transformer-fix-s-upper, StorageUnit-fix-p_dispatch-lower, StorageUnit-fix-p_dispatch-upper, StorageUnit-fix-p_store-lower, StorageUnit-fix-p_store-upper, StorageUnit-fix-state_of_charge-lower, StorageUnit-fix-state_of_charge-upper, Kirchhoff-Voltage-Law, StorageUnit-energy_balance were not assigned to the network.
Plot dispatch time series¶
In [10]:
Copied!
p_by_carrier = n.generators_t.p.T.groupby(n.generators.carrier).sum().T
to_drop = p_by_carrier.max()[p_by_carrier.max() < 1700].index
p_by_carrier.drop(to_drop, axis=1, inplace=True)
p_by_carrier.columns
p_by_carrier = n.generators_t.p.T.groupby(n.generators.carrier).sum().T
to_drop = p_by_carrier.max()[p_by_carrier.max() < 1700].index
p_by_carrier.drop(to_drop, axis=1, inplace=True)
p_by_carrier.columns
Out[10]:
Index(['Brown Coal', 'Gas', 'Hard Coal', 'Nuclear', 'Run of River', 'Solar',
'Wind Offshore', 'Wind Onshore'],
dtype='object', name='carrier')
In [11]:
Copied!
colors = {
"Brown Coal": "brown",
"Hard Coal": "k",
"Nuclear": "r",
"Run of River": "green",
"Wind Onshore": "blue",
"Solar": "yellow",
"Wind Offshore": "cyan",
"Waste": "orange",
"Gas": "orange",
}
# reorder
cols = [
"Nuclear",
"Run of River",
"Brown Coal",
"Hard Coal",
"Gas",
"Wind Offshore",
"Wind Onshore",
"Solar",
]
p_by_carrier = p_by_carrier[cols]
colors = {
"Brown Coal": "brown",
"Hard Coal": "k",
"Nuclear": "r",
"Run of River": "green",
"Wind Onshore": "blue",
"Solar": "yellow",
"Wind Offshore": "cyan",
"Waste": "orange",
"Gas": "orange",
}
# reorder
cols = [
"Nuclear",
"Run of River",
"Brown Coal",
"Hard Coal",
"Gas",
"Wind Offshore",
"Wind Onshore",
"Solar",
]
p_by_carrier = p_by_carrier[cols]
In [12]:
Copied!
c = [colors[col] for col in p_by_carrier.columns]
fig, ax = plt.subplots()
p_by_carrier.div(1e3).plot(kind="area", ax=ax, lw=0, color=c, alpha=0.7)
ax.legend(ncol=3, loc="upper left", bbox_to_anchor=(0, 1.02, 1, 0.2), frameon=False)
ax.set_ylabel("GW")
ax.set_xlabel("")
c = [colors[col] for col in p_by_carrier.columns]
fig, ax = plt.subplots()
p_by_carrier.div(1e3).plot(kind="area", ax=ax, lw=0, color=c, alpha=0.7)
ax.legend(ncol=3, loc="upper left", bbox_to_anchor=(0, 1.02, 1, 0.2), frameon=False)
ax.set_ylabel("GW")
ax.set_xlabel("")
Out[12]:
Text(0.5, 0, '')
Plot storage time series¶
In [13]:
Copied!
fig, ax = plt.subplots()
p_storage = n.storage_units_t.p.sum(axis=1)
state_of_charge = n.storage_units_t.state_of_charge.sum(axis=1)
p_storage.plot(label="Pumped hydro dispatch", ax=ax)
state_of_charge.plot(label="State of charge", ax=ax)
ax.axhline(0, color="k", lw=0.5, ls="--")
ax.legend()
ax.set_ylabel("MWh")
ax.set_xlabel("")
fig, ax = plt.subplots()
p_storage = n.storage_units_t.p.sum(axis=1)
state_of_charge = n.storage_units_t.state_of_charge.sum(axis=1)
p_storage.plot(label="Pumped hydro dispatch", ax=ax)
state_of_charge.plot(label="State of charge", ax=ax)
ax.axhline(0, color="k", lw=0.5, ls="--")
ax.legend()
ax.set_ylabel("MWh")
ax.set_xlabel("")
Out[13]:
Text(0.5, 0, '')
Line loading from optimisation¶
With the linear power flow, there is the following per unit loading at 3 AM:
In [14]:
Copied!
now = n.snapshots[4]
loading = n.lines_t.p0.loc[now] / n.lines.s_nom
loading.describe()
now = n.snapshots[4]
loading = n.lines_t.p0.loc[now] / n.lines.s_nom
loading.describe()
Out[14]:
count 852.000000 mean -0.003127 std 0.260218 min -0.700000 25% -0.127897 50% 0.003209 75% 0.121985 max 0.700000 dtype: float64
In [15]:
Copied!
fig, ax = plt.subplots(subplot_kw={"projection": ccrs.EqualEarth()})
n.plot(
ax=ax,
line_color=loading.abs(),
line_cmap="viridis",
title="Line loading",
bus_size=1e-3,
);
fig, ax = plt.subplots(subplot_kw={"projection": ccrs.EqualEarth()})
n.plot(
ax=ax,
line_color=loading.abs(),
line_cmap="viridis",
title="Line loading",
bus_size=1e-3,
);
Locational marginal prices¶
Let's have a look at the distribution of marginal prices:
In [16]:
Copied!
n.buses_t.marginal_price.loc[now].describe()
n.buses_t.marginal_price.loc[now].describe()
Out[16]:
count 585.000000 mean 15.737598 std 10.941995 min -10.397824 25% 6.992120 50% 15.841190 75% 25.048186 max 52.150120 Name: 2011-01-01 04:00:00, dtype: float64
In [17]:
Copied!
fig, ax = plt.subplots(subplot_kw={"projection": ccrs.PlateCarree()})
plt.hexbin(
n.buses.x,
n.buses.y,
gridsize=20,
C=n.buses_t.marginal_price.loc[now],
cmap="viridis",
zorder=-1,
)
n.plot(ax=ax, line_width=1, bus_size=0)
cb = plt.colorbar(location="right")
cb.set_label("Locational Marginal Price (€/MWh)")
fig, ax = plt.subplots(subplot_kw={"projection": ccrs.PlateCarree()})
plt.hexbin(
n.buses.x,
n.buses.y,
gridsize=20,
C=n.buses_t.marginal_price.loc[now],
cmap="viridis",
zorder=-1,
)
n.plot(ax=ax, line_width=1, bus_size=0)
cb = plt.colorbar(location="right")
cb.set_label("Locational Marginal Price (€/MWh)")
Curtailment¶
By considering how much power is available and how much is generated, you can see what share is curtailed:
In [18]:
Copied!
carrier = "Wind Onshore"
capacity = n.generators.groupby("carrier").sum().at[carrier, "p_nom"]
p_available = n.generators_t.p_max_pu.multiply(n.generators["p_nom"])
p_available_by_carrier = p_available.T.groupby(n.generators.carrier).sum().T
p_curtailed_by_carrier = p_available_by_carrier - p_by_carrier
carrier = "Wind Onshore"
capacity = n.generators.groupby("carrier").sum().at[carrier, "p_nom"]
p_available = n.generators_t.p_max_pu.multiply(n.generators["p_nom"])
p_available_by_carrier = p_available.T.groupby(n.generators.carrier).sum().T
p_curtailed_by_carrier = p_available_by_carrier - p_by_carrier
In [19]:
Copied!
p_df = pd.DataFrame(
{
carrier + " available": p_available_by_carrier[carrier],
carrier + " dispatched": p_by_carrier[carrier],
carrier + " curtailed": p_curtailed_by_carrier[carrier],
}
)
p_df[carrier + " capacity"] = capacity
p_df = pd.DataFrame(
{
carrier + " available": p_available_by_carrier[carrier],
carrier + " dispatched": p_by_carrier[carrier],
carrier + " curtailed": p_curtailed_by_carrier[carrier],
}
)
p_df[carrier + " capacity"] = capacity
In [20]:
Copied!
p_df.loc[p_df["Wind Onshore curtailed"] < 0, "Wind Onshore curtailed"] = 0
p_df.loc[p_df["Wind Onshore curtailed"] < 0, "Wind Onshore curtailed"] = 0
In [21]:
Copied!
fig, ax = plt.subplots()
p_df[[carrier + " dispatched", carrier + " curtailed"]].plot(kind="area", ax=ax, lw=0)
p_df[[carrier + " available", carrier + " capacity"]].plot(ax=ax)
ax.set_xlabel("")
ax.set_ylabel("Power [MW]")
ax.set_ylim([0, 40_000])
ax.legend()
fig, ax = plt.subplots()
p_df[[carrier + " dispatched", carrier + " curtailed"]].plot(kind="area", ax=ax, lw=0)
p_df[[carrier + " available", carrier + " capacity"]].plot(ax=ax)
ax.set_xlabel("")
ax.set_ylabel("Power [MW]")
ax.set_ylim([0, 40_000])
ax.legend()
Out[21]:
<matplotlib.legend.Legend at 0x7141147be490>
Non-Linear Power Flow¶
Now perform a full Newton-Raphson power flow on the first day. For the power flow, provide the setpoints from the optimisation.
In [22]:
Copied!
n.optimize.fix_optimal_dispatch()
n.optimize.fix_optimal_dispatch()
Set nearly all buses to PV, since we don't know what Q set points are. But we needs some PQ buses to ensure that the Jacobian is not singular.:
In [23]:
Copied!
n.generators.control = "PV"
f = n.generators[n.generators.bus == "492"]
n.generators.loc[f.index, "control"] = "PQ"
n.generators.control = "PV"
f = n.generators[n.generators.bus == "492"]
n.generators.loc[f.index, "control"] = "PQ"
Now, perform the non-linear PF.
In [24]:
Copied!
info = n.pf();
info = n.pf();
INFO:pypsa.network.power_flow:Performing non-linear load-flow on AC sub-network <pypsa.SubNetwork object at 0x71411f589f20> for snapshots DatetimeIndex(['2011-01-01 00:00:00', '2011-01-01 01:00:00',
'2011-01-01 02:00:00', '2011-01-01 03:00:00',
'2011-01-01 04:00:00', '2011-01-01 05:00:00',
'2011-01-01 06:00:00', '2011-01-01 07:00:00',
'2011-01-01 08:00:00', '2011-01-01 09:00:00',
'2011-01-01 10:00:00', '2011-01-01 11:00:00',
'2011-01-01 12:00:00', '2011-01-01 13:00:00',
'2011-01-01 14:00:00', '2011-01-01 15:00:00',
'2011-01-01 16:00:00', '2011-01-01 17:00:00',
'2011-01-01 18:00:00', '2011-01-01 19:00:00',
'2011-01-01 20:00:00', '2011-01-01 21:00:00',
'2011-01-01 22:00:00', '2011-01-01 23:00:00'],
dtype='datetime64[ns]', name='snapshot', freq=None)
Any failed to converge?
In [25]:
Copied!
(~info.converged).any().any()
(~info.converged).any().any()
Out[25]:
np.False_
With the non-linear load flow, there is the following per unit loading of the full thermal rating.
In [26]:
Copied!
(n.lines_t.p0.loc[now] / n.lines.s_nom).describe()
(n.lines_t.p0.loc[now] / n.lines.s_nom).describe()
Out[26]:
count 852.000000 mean 0.000242 std 0.262522 min -0.813731 25% -0.123203 50% 0.003104 75% 0.123850 max 0.827563 dtype: float64
Let's inspect the voltage angle differences across the lines have (in degrees)
In [27]:
Copied!
df = n.lines.copy()
for b in ["bus0", "bus1"]:
df = pd.merge(
df, n.buses_t.v_ang.loc[[now]].T, how="left", left_on=b, right_index=True
)
s = df[str(now) + "_x"] - df[str(now) + "_y"]
(s * 180 / np.pi).describe()
df = n.lines.copy()
for b in ["bus0", "bus1"]:
df = pd.merge(
df, n.buses_t.v_ang.loc[[now]].T, how="left", left_on=b, right_index=True
)
s = df[str(now) + "_x"] - df[str(now) + "_y"]
(s * 180 / np.pi).describe()
Out[27]:
count 852.000000 mean -0.022244 std 2.386535 min -12.158046 25% -0.463327 50% 0.001632 75% 0.534856 max 17.959258 dtype: float64
Plot the reactive power
In [28]:
Copied!
fig, ax = plt.subplots(subplot_kw={"projection": ccrs.EqualEarth()})
q = n.buses_t.q.loc[now]
bus_color = q.map(lambda x: "b" if x < 0 else "r")
n.plot(
bus_size=q.abs() / 1e4,
ax=ax,
bus_color=bus_color,
title="Reactive power feed-in\n(red = positive, blue = negative)",
);
fig, ax = plt.subplots(subplot_kw={"projection": ccrs.EqualEarth()})
q = n.buses_t.q.loc[now]
bus_color = q.map(lambda x: "b" if x < 0 else "r")
n.plot(
bus_size=q.abs() / 1e4,
ax=ax,
bus_color=bus_color,
title="Reactive power feed-in\n(red = positive, blue = negative)",
);
