Ctrl+K

Non-linear power flow after LOPF

Note

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

Non-linear power flow after LOPF#

In this example, the dispatch of generators is optimised using the linear OPF, then a non-linear power flow is run on the resulting dispatch.

Data sources#

Grid: based on SciGRID Version 0.2 which is based on OpenStreetMap.

Load size and location: based on Landkreise (NUTS 3) GDP and population.

Load time series: from ENTSO-E hourly data, scaled up uniformly by factor 1.12 (a simplification of the methodology in Schumacher, Hirth (2015)).

Conventional power plant capacities and locations: BNetzA list.

Wind and solar capacities and locations: EEG Stammdaten, based on http://www.energymap.info/download.html, which represents capacities at the end of 2014. Units without PLZ are removed.

Wind and solar time series: REatlas, Andresen et al, “Validation of Danish wind time series from a new global renewable energy atlas for energy system analysis,” Energy 93 (2015) 1074 - 1088.

Where SciGRID nodes have been split into 220kV and 380kV substations, all load and generation is attached to the 220kV substation.

Warnings#

The data behind the notebook is no longer supported. See PyPSA/pypsa-eur for a newer model that covers the whole of Europe.

This dataset is ONLY intended to demonstrate the capabilities of PyPSA and is NOT (yet) accurate enough to be used for research purposes.

Known problems include:

  1. Rough approximations have been made for missing grid data, e.g. 220kV-380kV transformers and connections between close sub-stations missing from OSM.

  2. There appears to be some unexpected congestion in parts of the network, which may mean for example that the load attachment method (by Voronoi cell overlap with Landkreise) isn’t working, particularly in regions with a high density of substations.

  3. Attaching power plants to the nearest high voltage substation may not reflect reality.

  4. There is no proper n-1 security in the calculations - this can either be simulated with a blanket e.g. 70% reduction in thermal limits (as done here) or a proper security constrained OPF (see e.g. http://pypsa.readthedocs.io/en/latest/examples/scigrid-sclopf.html).

  5. The borders and neighbouring countries are not represented.

  6. Hydroelectric power stations are not modelled accurately.

  7. The marginal costs are illustrative, not accurate.

  8. Only the first day of 2011 is in the github dataset, which is not representative. The full year of 2011 can be downloaded at http://www.pypsa.org/examples/scigrid-with-load-gen-trafos-2011.zip.

  9. The ENTSO-E total load for Germany may not be scaled correctly; it is scaled up uniformly by factor 1.12 (a simplification of the methodology in Schumacher, Hirth (2015), which suggests monthly factors).

  10. Biomass from the EEG Stammdaten are not read in at the moment.

  11. Power plant start up costs, ramping limits/costs, minimum loading rates are not considered.

[1]:

import pypsa
import numpy as np
import pandas as pd
import os
import matplotlib.pyplot as plt
import cartopy.crs as ccrs

ERROR 1: PROJ: proj_create_from_database: Open of /home/docs/checkouts/readthedocs.org/user_builds/pypsa/conda/latest/share/proj failed

[2]:

network = pypsa.examples.scigrid_de(from_master=True)

WARNING:pypsa.io:Importing network from PyPSA version v0.17.1 while current version is v0.28.0. Read the release notes at https://pypsa.readthedocs.io/en/latest/release_notes.html to prepare your network for import.
INFO:pypsa.io:Imported network scigrid-de.nc has buses, generators, lines, loads, storage_units, transformers

Plot the distribution of the load and of generating tech

[3]:

fig, ax = plt.subplots(
    1, 1, subplot_kw={"projection": ccrs.EqualEarth()}, figsize=(8, 8)
)

load_distribution = (
    network.loads_t.p_set.loc[network.snapshots[0]].groupby(network.loads.bus).sum()
)
network.plot(bus_sizes=1e-5 * load_distribution, ax=ax, title="Load distribution");

/home/docs/checkouts/readthedocs.org/user_builds/pypsa/conda/latest/lib/python3.11/site-packages/cartopy/mpl/feature_artist.py:144: UserWarning: facecolor will have no effect as it has been defined as "never".
  warnings.warn('facecolor will have no effect as it has been '
../_images/examples_scigrid-lopf-then-pf_4_1.png 
[4]:

network.generators.groupby("carrier")["p_nom"].sum()

[4]:

carrier
Brown Coal       20879.500000
Gas              23913.130000
Geothermal          31.700000
Hard Coal        25312.600000
Multiple           152.700000
Nuclear          12068.000000
Oil               2710.200000
Other             3027.800000
Run of River      3999.100000
Solar            37041.524779
Storage Hydro     1445.000000
Waste             1645.900000
Wind Offshore     2973.500000
Wind Onshore     37339.895329
Name: p_nom, dtype: float64

[5]:

network.storage_units.groupby("carrier")["p_nom"].sum()

[5]:

carrier
Pumped Hydro    9179.5
Name: p_nom, dtype: float64

[6]:

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 = network.generators[network.generators.carrier == tech]
    gen_distribution = (
        gens.groupby("bus").sum()["p_nom"].reindex(network.buses.index, fill_value=0.0)
    )
    network.plot(ax=ax, bus_sizes=2e-5 * gen_distribution)
    ax.set_title(tech)
fig.tight_layout()

/home/docs/checkouts/readthedocs.org/user_builds/pypsa/conda/latest/lib/python3.11/site-packages/cartopy/mpl/feature_artist.py:144: UserWarning: facecolor will have no effect as it has been defined as "never".
  warnings.warn('facecolor will have no effect as it has been '
/home/docs/checkouts/readthedocs.org/user_builds/pypsa/conda/latest/lib/python3.11/site-packages/cartopy/mpl/feature_artist.py:144: UserWarning: facecolor will have no effect as it has been defined as "never".
  warnings.warn('facecolor will have no effect as it has been '
/home/docs/checkouts/readthedocs.org/user_builds/pypsa/conda/latest/lib/python3.11/site-packages/cartopy/mpl/feature_artist.py:144: UserWarning: facecolor will have no effect as it has been defined as "never".
  warnings.warn('facecolor will have no effect as it has been '
/home/docs/checkouts/readthedocs.org/user_builds/pypsa/conda/latest/lib/python3.11/site-packages/cartopy/mpl/feature_artist.py:144: UserWarning: facecolor will have no effect as it has been defined as "never".
  warnings.warn('facecolor will have no effect as it has been '
/home/docs/checkouts/readthedocs.org/user_builds/pypsa/conda/latest/lib/python3.11/site-packages/cartopy/mpl/feature_artist.py:144: UserWarning: facecolor will have no effect as it has been defined as "never".
  warnings.warn('facecolor will have no effect as it has been '
/home/docs/checkouts/readthedocs.org/user_builds/pypsa/conda/latest/lib/python3.11/site-packages/cartopy/mpl/feature_artist.py:144: UserWarning: facecolor will have no effect as it has been defined as "never".
  warnings.warn('facecolor will have no effect as it has been '
../_images/examples_scigrid-lopf-then-pf_7_1.png

Run Linear Optimal Power Flow on the first day of 2011.

To approximate n-1 security and allow room for reactive power flows, don’t allow any line to be loaded above 70% of their thermal rating

[7]:

contingency_factor = 0.7
network.lines.s_max_pu = contingency_factor

There are some infeasibilities without small extensions

[8]:

network.lines.loc[["316", "527", "602"], "s_nom"] = 1715

We performing a linear OPF for one day, 4 snapshots at a time.

[9]:

group_size = 4
network.storage_units.state_of_charge_initial = 0.0

for i in range(int(24 / group_size)):
    # set the initial state of charge based on previous round
    if i:
        network.storage_units.state_of_charge_initial = (
            network.storage_units_t.state_of_charge.loc[
                network.snapshots[group_size * i - 1]
            ]
        )
    network.optimize(
        network.snapshots[group_size * i : group_size * i + group_size],
        solver_name="cbc",
    )

WARNING:pypsa.components: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='Transformer')
WARNING:pypsa.components: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='Transformer')
INFO:linopy.model: Solve problem using Cbc solver
INFO:linopy.io: Writing time: 0.17s
INFO:linopy.solvers:Welcome to the CBC MILP Solver
Version: 2.10.10
Build Date: Apr 19 2023

command line - cbc -printingOptions all -import /tmp/linopy-problem-8gbb84oa.lp -solve -solu /tmp/linopy-solve-_a4k_o2l.sol (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 2429 (-21399) rows, 4357 (-5583) columns and 15341 (-28565) elements
Perturbing problem by 0.001% of 2769.2736 - largest nonzero change 0.00099246172 ( 0.015099382%) - largest zero change 0.00099142717
0  Obj 88131.77 Primal inf 4920261.9 (2224) Dual inf 374.16895 (4)
123  Obj 18780.455 Primal inf 3391702.9 (2173)
246  Obj 18780.455 Primal inf 3046627.4 (2112)
369  Obj 18780.455 Primal inf 2895580.8 (2095)
492  Obj 18780.772 Primal inf 3094322.3 (2056)
615  Obj 18782.245 Primal inf 2875702.9 (1988)
738  Obj 18783.485 Primal inf 3065605.6 (1958)
861  Obj 18783.961 Primal inf 3176424.6 (1912)
984  Obj 18789.18 Primal inf 2455791.8 (1793)
1107  Obj 18791.617 Primal inf 3951671.5 (1752)
1230  Obj 18795.35 Primal inf 6158655.2 (1748)
1353  Obj 18797.413 Primal inf 21602140 (1620)
1476  Obj 22617.311 Primal inf 21429797 (1500)
1599  Obj 22620.718 Primal inf 30347474 (1366)
1722  Obj 22626.175 Primal inf 51192620 (1242)
1845  Obj 22636.272 Primal inf 53452705 (1197)
1968  Obj 24965.395 Primal inf 1.2429637e+08 (1303)
2091  Obj 67923.979 Primal inf 23572094 (1055)
2214  Obj 200877.69 Primal inf 1.5632781e+08 (968)
2337  Obj 381603.01 Primal inf 1.2251578e+08 (649)
2460  Obj 497008.72 Primal inf 269221.95 (543)
2583  Obj 907814.35 Primal inf 94615.812 (358)
2706  Obj 1134071.5 Primal inf 40615.798 (310)
2829  Obj 1339938.1 Primal inf 45147.776 (245)
2952  Obj 1419000 Primal inf 3171.3744 (93)
3075  Obj 1448893.7 Primal inf 261.52701 (35)
3142  Obj 1449690.8
3142  Obj 1449687.3 Dual inf 6.3634756e-06 (1)
3144  Obj 1449687.3
Optimal - objective value 1449687.3
After Postsolve, objective 1449687.3, infeasibilities - dual 186.92063 (6), primal 1.3870423e-06 (6)
Presolved model was optimal, full model needs cleaning up
0  Obj 1449687.3 Dual inf 12.773483 (15)
End of values pass after 0 iterations
0  Obj 1449687.3 Dual inf 12.773483 (15) w.o. free dual inf (13)
End of values pass after 0 iterations
0  Obj 1449687.3 Dual inf 12.773483 (15) w.o. free dual inf (12)
End of values pass after 10 iterations
10  Obj 1449687.3 Dual inf 10.904277 (9) w.o. free dual inf (6)
13  Obj 1449687.3
Optimal - objective value 1449687.3
Optimal objective 1449687.25 - 3157 iterations time 0.382, Presolve 0.03
Total time (CPU seconds):       0.59   (Wallclock seconds):       0.52


INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 9940 primals, 23828 duals
Objective: 1.45e+06
Solver model: not available
Solver message: Optimal - objective value 1449687.25014697


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.
WARNING:pypsa.components: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='Transformer')
WARNING:pypsa.components: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='Transformer')
INFO:linopy.model: Solve problem using Cbc solver
INFO:linopy.io: Writing time: 0.17s
INFO:linopy.solvers:Welcome to the CBC MILP Solver
Version: 2.10.10
Build Date: Apr 19 2023

command line - cbc -printingOptions all -import /tmp/linopy-problem-cusci0ko.lp -solve -solu /tmp/linopy-solve-2igqt7vq.sol (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 2451 (-21377) rows, 4387 (-5553) columns and 15471 (-28435) elements
Perturbing problem by 0.001% of 2769.2736 - largest nonzero change 0.00099246172 ( 0.014954013%) - largest zero change 0.00099142717
0  Obj 81174.124 Primal inf 4944643.9 (2252) Dual inf 476.35378 (4)
124  Obj 11822.81 Primal inf 3577799.5 (2197)
248  Obj 11823.847 Primal inf 3089038.7 (2148)
372  Obj 11824.185 Primal inf 2687198 (2094)
496  Obj 11825.23 Primal inf 3245529.8 (2076)
620  Obj 11826.05 Primal inf 3887022.3 (2023)
744  Obj 11826.51 Primal inf 4047441.3 (2012)
868  Obj 11828.548 Primal inf 4221995.3 (1940)
992  Obj 11830.658 Primal inf 18626562 (1880)
1116  Obj 11832.699 Primal inf 16559669 (1813)
1240  Obj 11835.71 Primal inf 55927543 (1736)
1364  Obj 11838.701 Primal inf 4040983.3 (1611)
1488  Obj 11843.389 Primal inf 1998223.5 (1433)
1612  Obj 11849.733 Primal inf 5826938.4 (1576)
1736  Obj 11854.117 Primal inf 3195252.1 (1404)
1860  Obj 11858.459 Primal inf 1625623.4 (1170)
1984  Obj 11862.808 Primal inf 2424545.4 (1122)
2108  Obj 11868.555 Primal inf 2055960.9 (981)
2232  Obj 30821.999 Primal inf 1081704.4 (891)
2356  Obj 57488.864 Primal inf 4935718.2 (944)
2480  Obj 157151.82 Primal inf 4498932.6 (726)
2604  Obj 340564.67 Primal inf 107587.92 (397)
2728  Obj 513792.38 Primal inf 48916.892 (285)
2852  Obj 666145.84 Primal inf 753271.2 (359)
2976  Obj 780216.62 Primal inf 16038.716 (162)
3100  Obj 862337.25 Primal inf 2187.7859 (82)
3224  Obj 873762.56 Primal inf 15.459476 (7)
3241  Obj 873769.94
Optimal - objective value 873763.9
After Postsolve, objective 873763.9, infeasibilities - dual 78.034695 (3), primal 9.1550862e-07 (3)
Presolved model was optimal, full model needs cleaning up
0  Obj 873763.9 Dual inf 0.78034665 (3)
End of values pass after 3 iterations
3  Obj 873763.9
Optimal - objective value 873763.9
Optimal objective 873763.9023 - 3244 iterations time 0.402, Presolve 0.04
Total time (CPU seconds):       0.62   (Wallclock seconds):       0.53


INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 9940 primals, 23828 duals
Objective: 8.74e+05
Solver model: not available
Solver message: Optimal - objective value 873763.90232907


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.
WARNING:pypsa.components: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='Transformer')
WARNING:pypsa.components: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='Transformer')
INFO:linopy.model: Solve problem using Cbc solver
INFO:linopy.io: Writing time: 0.17s
INFO:linopy.solvers:Welcome to the CBC MILP Solver
Version: 2.10.10
Build Date: Apr 19 2023

command line - cbc -printingOptions all -import /tmp/linopy-problem-g3h6e7wu.lp -solve -solu /tmp/linopy-solve-sdksgxqz.sol (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 2520 (-21308) rows, 4463 (-5477) columns and 15683 (-28223) elements
Perturbing problem by 0.001% of 3148.2846 - largest nonzero change 0.00099819563 ( 0.0069530294%) - largest zero change 0.00099758325
0  Obj -37.859231 Primal inf 5101322.9 (2320)
125  Obj -37.859231 Primal inf 3414946.6 (2269)
250  Obj -37.236804 Primal inf 3530110.8 (2210)
375  Obj -36.71281 Primal inf 3442912.3 (2156)
500  Obj -36.544476 Primal inf 3115512.6 (2112)
625  Obj -36.17333 Primal inf 3705271.2 (2066)
750  Obj -35.614977 Primal inf 3004088.7 (2018)
875  Obj -35.350996 Primal inf 3887020.7 (1996)
1000  Obj -33.981271 Primal inf 3343218.2 (1941)
1125  Obj -31.717471 Primal inf 2483275.3 (1808)
1250  Obj -30.102391 Primal inf 1912510.8 (1722)
1375  Obj -28.011261 Primal inf 1979451.8 (1628)
1500  Obj -26.044779 Primal inf 2329222.5 (1550)
1625  Obj 7639.0388 Primal inf 6759859.7 (1555)
1750  Obj 7642.2875 Primal inf 22736410 (1567)
1875  Obj 7645.6722 Primal inf 39703221 (1483)
2000  Obj 7648.6847 Primal inf 10450888 (1308)
2125  Obj 7658.9477 Primal inf 29501422 (1175)
2250  Obj 7662.5535 Primal inf 8822329.7 (1124)
2375  Obj 7668.0361 Primal inf 5668588.7 (698)
2500  Obj 59225.32 Primal inf 325223.44 (586)
2625  Obj 131756.97 Primal inf 138426.08 (463)
2750  Obj 388348.75 Primal inf 68788.219 (393)
2875  Obj 499863.1 Primal inf 576740 (358)
3000  Obj 632724.89 Primal inf 118701.55 (312)
3125  Obj 746299.75 Primal inf 14720.878 (174)
3250  Obj 783521.6 Primal inf 2252.4206 (71)
3372  Obj 790793.8
Optimal - objective value 790786.59
After Postsolve, objective 790786.59, infeasibilities - dual 0 (0), primal 0 (0)
Optimal objective 790786.5949 - 3372 iterations time 0.372, Presolve 0.03
Total time (CPU seconds):       0.60   (Wallclock seconds):       0.51


INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 9940 primals, 23828 duals
Objective: 7.91e+05
Solver model: not available
Solver message: Optimal - objective value 790786.59488191


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.
WARNING:pypsa.components: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='Transformer')
WARNING:pypsa.components: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='Transformer')
INFO:linopy.model: Solve problem using Cbc solver
INFO:linopy.io: Writing time: 0.17s
INFO:linopy.solvers:Welcome to the CBC MILP Solver
Version: 2.10.10
Build Date: Apr 19 2023

command line - cbc -printingOptions all -import /tmp/linopy-problem-eekrsw1u.lp -solve -solu /tmp/linopy-solve-k0s0732s.sol (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 2497 (-21331) rows, 4449 (-5491) columns and 15642 (-28264) elements
Perturbing problem by 0.001% of 2769.2736 - largest nonzero change 0.00099944084 ( 0.008627573%) - largest zero change 0.00099884126
0  Obj 25309.974 Primal inf 5242665.2 (2315) Dual inf 151.02464 (1)
124  Obj 7972.1452 Primal inf 3496955 (2255)
248  Obj 7972.1452 Primal inf 3234550.5 (2208)
372  Obj 7972.8452 Primal inf 3003954.6 (2146)
496  Obj 7973.1893 Primal inf 2583174.5 (2109)
620  Obj 7973.6596 Primal inf 2884825.9 (2069)
744  Obj 7973.9487 Primal inf 3553314.1 (2050)
868  Obj 7975.2235 Primal inf 7348484.2 (2031)
992  Obj 7981.065 Primal inf 7701551.1 (1914)
1116  Obj 7982.4484 Primal inf 10525904 (1861)
1240  Obj 7984.5622 Primal inf 4686591.9 (1735)
1364  Obj 12738.193 Primal inf 5593207.6 (1659)
1488  Obj 12740.427 Primal inf 7540643.4 (1640)
1612  Obj 23724.89 Primal inf 6989181.2 (1539)
1736  Obj 23728.525 Primal inf 4269760.6 (1312)
1860  Obj 23732.443 Primal inf 4003323.4 (1177)
1984  Obj 23737.232 Primal inf 6557694.9 (1298)
2108  Obj 23740.723 Primal inf 11695870 (1338)
2232  Obj 23870.813 Primal inf 14979844 (1203)
2356  Obj 99205.314 Primal inf 5743530.5 (963)
2480  Obj 133195.14 Primal inf 9262255.4 (985)
2604  Obj 297144.26 Primal inf 23955882 (810)
2728  Obj 484259.31 Primal inf 2333847.8 (514)
2852  Obj 925240.06 Primal inf 68855.099 (329)
2976  Obj 1276166.9 Primal inf 14169.551 (206)
3100  Obj 1406228.3 Primal inf 35327.859 (173)
3224  Obj 1454698.7 Primal inf 824.92415 (36)
3261  Obj 1455241.5
Optimal - objective value 1455232.7
After Postsolve, objective 1455232.7, infeasibilities - dual 0 (0), primal 0 (0)
Optimal objective 1455232.678 - 3261 iterations time 0.362, Presolve 0.03
Total time (CPU seconds):       0.58   (Wallclock seconds):       0.49


INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 9940 primals, 23828 duals
Objective: 1.46e+06
Solver model: not available
Solver message: Optimal - objective value 1455232.67827293


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.
WARNING:pypsa.components: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='Transformer')
WARNING:pypsa.components: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='Transformer')
INFO:linopy.model: Solve problem using Cbc solver
INFO:linopy.io: Writing time: 0.18s
INFO:linopy.solvers:Welcome to the CBC MILP Solver
Version: 2.10.10
Build Date: Apr 19 2023

command line - cbc -printingOptions all -import /tmp/linopy-problem-9a0gt5je.lp -solve -solu /tmp/linopy-solve-gwyx_54i.sol (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 2434 (-21394) rows, 4388 (-5552) columns and 15512 (-28394) elements
Perturbing problem by 0.001% of 2769.2736 - largest nonzero change 0.00099975902 ( 0.009641762%) - largest zero change 0.00099920441
0  Obj 97177.304 Primal inf 5426833 (2270) Dual inf 604.09852 (4)
123  Obj 27997.531 Primal inf 3474041 (2224)
246  Obj 27997.531 Primal inf 3271740.7 (2189)
369  Obj 27997.531 Primal inf 3855753.1 (2140)
492  Obj 27997.897 Primal inf 3735841.6 (2086)
615  Obj 27998.371 Primal inf 3124968.2 (2020)
738  Obj 27998.812 Primal inf 3638331.3 (1986)
861  Obj 28000.397 Primal inf 2314356.9 (1870)
984  Obj 28008.413 Primal inf 3597507.7 (1821)
1107  Obj 28012.427 Primal inf 2380464.6 (1714)
1230  Obj 28014.319 Primal inf 8511238.5 (1700)
1353  Obj 59476.18 Primal inf 6096965.6 (1570)
1476  Obj 59479.151 Primal inf 9035886.9 (1457)
1599  Obj 59482.544 Primal inf 8341277 (1378)
1722  Obj 59486.062 Primal inf 7041870.2 (1279)
1845  Obj 59491.501 Primal inf 18741033 (1538)
1968  Obj 73342.067 Primal inf 6221362.2 (1165)
2091  Obj 127515.18 Primal inf 75398656 (1261)
2214  Obj 260744.69 Primal inf 5.3327485e+08 (1038)
2337  Obj 539356.8 Primal inf 6.1773447e+08 (765)
2460  Obj 880279.04 Primal inf 1.206367e+08 (711)
2583  Obj 1175776.6 Primal inf 18096431 (600)
2706  Obj 1601376.9 Primal inf 33945667 (442)
2829  Obj 2084791.1 Primal inf 58113.371 (350)
2952  Obj 2398552.5 Primal inf 20031.596 (242)
3075  Obj 2591167.8 Primal inf 5425.5718 (144)
3198  Obj 2635716.6 Primal inf 3109.0375 (107)
3321  Obj 2647305.6 Primal inf 256.39086 (44)
3368  Obj 2647916.6
Optimal - objective value 2647902.5
After Postsolve, objective 2647902.5, infeasibilities - dual 217.1748 (5), primal 1.6669476e-06 (5)
Presolved model was optimal, full model needs cleaning up
0  Obj 2647902.5 Primal inf 6.3133663e-07 (1) Dual inf 1e+08 (6)
End of values pass after 9 iterations
9  Obj 2647902.5
Optimal - objective value 2647902.5
Optimal objective 2647902.52 - 3377 iterations time 0.462, Presolve 0.04
Total time (CPU seconds):       0.68   (Wallclock seconds):       0.58


INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 9940 primals, 23828 duals
Objective: 2.65e+06
Solver model: not available
Solver message: Optimal - objective value 2647902.51963847


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.
WARNING:pypsa.components: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='Transformer')
WARNING:pypsa.components: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='Transformer')
INFO:linopy.model: Solve problem using Cbc solver
INFO:linopy.io: Writing time: 0.17s
INFO:linopy.solvers:Welcome to the CBC MILP Solver
Version: 2.10.10
Build Date: Apr 19 2023

command line - cbc -printingOptions all -import /tmp/linopy-problem-pbxmj4k9.lp -solve -solu /tmp/linopy-solve-ytk2kdvm.sol (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 2442 (-21386) rows, 4392 (-5548) columns and 15524 (-28382) elements
Perturbing problem by 0.001% of 9068.3747 - largest nonzero change 0.0009413617 ( 0.0033681695%) - largest zero change 0.00094111171
0  Obj 92714.595 Primal inf 5406628.1 (2269) Dual inf 1225.4155 (4)
123  Obj 23363.589 Primal inf 3455547 (2205)
246  Obj 23363.589 Primal inf 3130258.2 (2161)
369  Obj 23363.715 Primal inf 3858547.6 (2152)
492  Obj 23363.925 Primal inf 4360359.5 (2095)
615  Obj 23364.181 Primal inf 3534015.5 (2048)
738  Obj 23364.753 Primal inf 2613352.5 (1980)
861  Obj 23367.792 Primal inf 3552238.4 (1919)
984  Obj 23368.355 Primal inf 3534460.1 (1843)
1107  Obj 23369.923 Primal inf 2132424.4 (1717)
1230  Obj 23370.784 Primal inf 2059777.6 (1647)
1353  Obj 23371.676 Primal inf 1705599.3 (1529)
1476  Obj 40822.922 Primal inf 2675337 (1537)
1599  Obj 40824.168 Primal inf 8371287.9 (1463)
1722  Obj 40826.952 Primal inf 4096460.2 (1395)
1845  Obj 40829.66 Primal inf 4355786.4 (1338)
1968  Obj 40832.089 Primal inf 1774061 (1234)
2091  Obj 42856.307 Primal inf 1704135.1 (1136)
2195  Obj 138986.1 Primal inf 9.2979695e+08 (1186)
2318  Obj 421904.27 Primal inf 77148776 (1150)
2441  Obj 630756.14 Primal inf 1.0802175e+09 (858)
2564  Obj 1069512.5 Primal inf 2.4410129e+08 (615)
2687  Obj 1386353.6 Primal inf 224602.48 (410)
2810  Obj 1614619.7 Primal inf 88773.98 (344)
2933  Obj 1838260.4 Primal inf 21680.684 (210)
3056  Obj 1935716.9 Primal inf 20004.402 (181)
3179  Obj 2007790.5 Primal inf 5887.2821 (97)
3302  Obj 2099341.4 Primal inf 1242.7627 (54)
3425  Obj 2142605.2 Primal inf 92.334755 (34)
3459  Obj 2142961.1
Optimal - objective value 2142956.9
After Postsolve, objective 2142956.9, infeasibilities - dual 99.086763 (3), primal 1.7325076e-07 (3)
Presolved model was optimal, full model needs cleaning up
0  Obj 2142956.9 Dual inf 0.99086733 (3)
End of values pass after 7 iterations
7  Obj 2142956.9
Optimal - objective value 2142956.9
Optimal objective 2142956.889 - 3466 iterations time 0.412, Presolve 0.04
Total time (CPU seconds):       0.63   (Wallclock seconds):       0.55


INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 9940 primals, 23828 duals
Objective: 2.14e+06
Solver model: not available
Solver message: Optimal - objective value 2142956.88872811


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.

[10]:

p_by_carrier = network.generators_t.p.groupby(network.generators.carrier, axis=1).sum()
p_by_carrier.drop(
    (p_by_carrier.max()[p_by_carrier.max() < 1700.0]).index, axis=1, inplace=True
)
p_by_carrier.columns

/tmp/ipykernel_4806/901685434.py:1: FutureWarning: DataFrame.groupby with axis=1 is deprecated. Do `frame.T.groupby(...)` without axis instead.
  p_by_carrier = network.generators_t.p.groupby(network.generators.carrier, axis=1).sum()

[10]:

Index(['Brown Coal', 'Gas', 'Hard Coal', 'Nuclear', 'Run of River', 'Solar',
       'Wind Offshore', 'Wind Onshore'],
      dtype='object', name='carrier')

[11]:

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]

[12]:

c = [colors[col] for col in p_by_carrier.columns]

fig, ax = plt.subplots(figsize=(12, 6))
(p_by_carrier / 1e3).plot(kind="area", ax=ax, linewidth=4, color=c, alpha=0.7)
ax.legend(ncol=4, loc="upper left")
ax.set_ylabel("GW")
ax.set_xlabel("")
fig.tight_layout()
../_images/examples_scigrid-lopf-then-pf_16_0.png 
[13]:

fig, ax = plt.subplots(figsize=(12, 6))

p_storage = network.storage_units_t.p.sum(axis=1)
state_of_charge = network.storage_units_t.state_of_charge.sum(axis=1)
p_storage.plot(label="Pumped hydro dispatch", ax=ax, linewidth=3)
state_of_charge.plot(label="State of charge", ax=ax, linewidth=3)

ax.legend()
ax.grid()
ax.set_ylabel("MWh")
ax.set_xlabel("")
fig.tight_layout()
../_images/examples_scigrid-lopf-then-pf_17_0.png 
[14]:

now = network.snapshots[4]

With the linear load flow, there is the following per unit loading:

[15]:

loading = network.lines_t.p0.loc[now] / network.lines.s_nom
loading.describe()

[15]:

count    852.000000
mean      -0.002892
std        0.258607
min       -0.700000
25%       -0.129689
50%        0.003084
75%        0.122643
max        0.700000
dtype: float64

[16]:

fig, ax = plt.subplots(subplot_kw={"projection": ccrs.EqualEarth()}, figsize=(9, 9))
network.plot(
    ax=ax,
    line_colors=abs(loading),
    line_cmap=plt.cm.jet,
    title="Line loading",
    bus_sizes=1e-3,
    bus_alpha=0.7,
)
fig.tight_layout();

/home/docs/checkouts/readthedocs.org/user_builds/pypsa/conda/latest/lib/python3.11/site-packages/cartopy/mpl/feature_artist.py:144: UserWarning: facecolor will have no effect as it has been defined as "never".
  warnings.warn('facecolor will have no effect as it has been '
../_images/examples_scigrid-lopf-then-pf_21_1.png

Let’s have a look at the marginal prices

[17]:

network.buses_t.marginal_price.loc[now].describe()

[17]:

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

[18]:

fig, ax = plt.subplots(subplot_kw={"projection": ccrs.PlateCarree()}, figsize=(8, 8))

plt.hexbin(
    network.buses.x,
    network.buses.y,
    gridsize=20,
    C=network.buses_t.marginal_price.loc[now],
    cmap=plt.cm.jet,
    zorder=3,
)
network.plot(ax=ax, line_widths=pd.Series(0.5, network.lines.index), bus_sizes=0)

cb = plt.colorbar(location="bottom")
cb.set_label("Locational Marginal Price (EUR/MWh)")
fig.tight_layout()

/home/docs/checkouts/readthedocs.org/user_builds/pypsa/conda/latest/lib/python3.11/site-packages/cartopy/mpl/feature_artist.py:144: UserWarning: facecolor will have no effect as it has been defined as "never".
  warnings.warn('facecolor will have no effect as it has been '
../_images/examples_scigrid-lopf-then-pf_24_1.png

Curtailment variable#

By considering how much power is available and how much is generated, you can see what share is curtailed:

[19]:

carrier = "Wind Onshore"

capacity = network.generators.groupby("carrier").sum().at[carrier, "p_nom"]
p_available = network.generators_t.p_max_pu.multiply(network.generators["p_nom"])
p_available_by_carrier = p_available.groupby(network.generators.carrier, axis=1).sum()
p_curtailed_by_carrier = p_available_by_carrier - p_by_carrier

/tmp/ipykernel_4806/2097778438.py:5: FutureWarning: DataFrame.groupby with axis=1 is deprecated. Do `frame.T.groupby(...)` without axis instead.
  p_available_by_carrier = p_available.groupby(network.generators.carrier, axis=1).sum()

[20]:

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

[21]:

p_df["Wind Onshore curtailed"][p_df["Wind Onshore curtailed"] < 0.0] = 0.0

/tmp/ipykernel_4806/1546145974.py:1: FutureWarning: ChainedAssignmentError: behaviour will change in pandas 3.0!
You are setting values through chained assignment. Currently this works in certain cases, but when using Copy-on-Write (which will become the default behaviour in pandas 3.0) this will never work to update the original DataFrame or Series, because the intermediate object on which we are setting values will behave as a copy.
A typical example is when you are setting values in a column of a DataFrame, like:

df["col"][row_indexer] = value

Use `df.loc[row_indexer, "col"] = values` instead, to perform the assignment in a single step and ensure this keeps updating the original `df`.

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy

  p_df["Wind Onshore curtailed"][p_df["Wind Onshore curtailed"] < 0.0] = 0.0

[22]:

fig, ax = plt.subplots(figsize=(10, 4))
p_df[[carrier + " dispatched", carrier + " curtailed"]].plot(
    kind="area", ax=ax, linewidth=3
)
p_df[[carrier + " available", carrier + " capacity"]].plot(ax=ax, linewidth=3)

ax.set_xlabel("")
ax.set_ylabel("Power [MW]")
ax.set_ylim([0, 40000])
ax.legend()
fig.tight_layout()
../_images/examples_scigrid-lopf-then-pf_29_0.png

Non-Linear Power Flow#

Now perform a full Newton-Raphson power flow on the first hour. For the PF, set the P to the optimised P.

[23]:

network.generators_t.p_set = network.generators_t.p
network.storage_units_t.p_set = network.storage_units_t.p

Set all buses to PV, since we don’t know what Q set points are

[24]:

network.generators.control = "PV"

# Need some PQ buses so that Jacobian doesn't break
f = network.generators[network.generators.bus == "492"]
network.generators.loc[f.index, "control"] = "PQ"

Now, perform the non-linear PF.

[25]:

info = network.pf();

INFO:pypsa.pf:Performing non-linear load-flow on AC sub-network SubNetwork 0 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)
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.033733 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.033595 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.033771 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.033657 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.033154 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.033438 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.032925 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.032873 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.033531 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.033713 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.033082 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.033127 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.033061 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.033537 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.033695 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.032922 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.033554 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.033694 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.033304 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.033150 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.033128 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.033046 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.033302 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.033164 seconds

Any failed to converge?

[26]:

(~info.converged).any().any()

[26]:

False

With the non-linear load flow, there is the following per unit loading of the full thermal rating.

[27]:

(network.lines_t.p0.loc[now] / network.lines.s_nom).describe()

[27]:

count    852.000000
mean       0.000435
std        0.260958
min       -0.813383
25%       -0.125478
50%        0.003032
75%        0.126690
max        0.827203
dtype: float64

Let’s inspect the voltage angle differences across the lines have (in degrees)

[28]:

df = network.lines.copy()

for b in ["bus0", "bus1"]:
    df = pd.merge(
        df, network.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()

[28]:

count    852.000000
mean      -0.022968
std        2.373833
min      -12.158423
25%       -0.462981
50%        0.001586
75%        0.537443
max       17.959258
dtype: float64

Plot the reactive power

[29]:

fig, ax = plt.subplots(subplot_kw={"projection": ccrs.EqualEarth()}, figsize=(9, 9))

q = network.buses_t.q.loc[now]
bus_colors = pd.Series("r", network.buses.index)
bus_colors[q < 0.0] = "b"

network.plot(
    bus_sizes=1e-4 * abs(q),
    ax=ax,
    bus_colors=bus_colors,
    title="Reactive power feed-in (red=+ve, blue=-ve)",
);

/home/docs/checkouts/readthedocs.org/user_builds/pypsa/conda/latest/lib/python3.11/site-packages/cartopy/mpl/feature_artist.py:144: UserWarning: facecolor will have no effect as it has been defined as "never".
  warnings.warn('facecolor will have no effect as it has been '
../_images/examples_scigrid-lopf-then-pf_43_1.png
Contents