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.27.1. 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/style.py:76: 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/style.py:76: 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.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-247rcnlv.lp -solve -solu /tmp/linopy-solve-vdbffk51.sol (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 2410 (-21418) rows, 4349 (-5591) columns and 14945 (-28573) elements
Perturbing problem by 0.001% of 2769.2736 - largest nonzero change 0.00099944084 ( 0.015169862%) - largest zero change 0.00099884126
0  Obj 88139.312 Primal inf 5864837.9 (2205) Dual inf 371.77707 (4)
123  Obj 18787.997 Primal inf 3568080.9 (2162)
246  Obj 18787.997 Primal inf 3361257.3 (2111)
369  Obj 18788.186 Primal inf 3341077.2 (2089)
492  Obj 18789.315 Primal inf 2977841.4 (2024)
615  Obj 18789.773 Primal inf 5639654.1 (1985)
738  Obj 18790.821 Primal inf 2526002.5 (1912)
861  Obj 18792.007 Primal inf 2508693.3 (1856)
984  Obj 18795.794 Primal inf 2497340.3 (1825)
1107  Obj 18800.582 Primal inf 3076166.1 (1729)
1230  Obj 18802.281 Primal inf 3440737.3 (1642)
1353  Obj 18804.784 Primal inf 18097362 (1643)
1476  Obj 22623.938 Primal inf 20329479 (1532)
1599  Obj 22627.193 Primal inf 5547819.2 (1392)
1722  Obj 22631.563 Primal inf 4947738.4 (1389)
1845  Obj 22635.686 Primal inf 6232965.1 (1282)
1968  Obj 38971.764 Primal inf 8800510.8 (1281)
2091  Obj 39447.294 Primal inf 5584555.5 (1141)
2214  Obj 64731.652 Primal inf 93322443 (1273)
2337  Obj 106169.32 Primal inf 13596759 (899)
2460  Obj 399976.61 Primal inf 4000202.3 (650)
2583  Obj 633424.38 Primal inf 1785497.9 (483)
2706  Obj 851632.13 Primal inf 163208.91 (437)
2829  Obj 1238322.8 Primal inf 35698.282 (278)
2952  Obj 1372608.6 Primal inf 218380.93 (380)
3075  Obj 1437363.4 Primal inf 3831.7001 (87)
3183  Obj 1449694.1
3183  Obj 1449687.3 Dual inf 7.3445901e-06 (1)
3184  Obj 1449687.3
Optimal - objective value 1449687.3
After Postsolve, objective 1449687.3, infeasibilities - dual 214.20597 (7), primal 1.7945025e-06 (7)
Presolved model was optimal, full model needs cleaning up
0  Obj 1449687.3 Dual inf 11.934102 (14)
End of values pass after 0 iterations
0  Obj 1449687.3 Dual inf 11.934102 (14) w.o. free dual inf (13)
End of values pass after 9 iterations
9  Obj 1449687.3 Dual inf 9.7920427 (7) w.o. free dual inf (6)
10  Obj 1449687.3
Optimal - objective value 1449687.3
Optimal objective 1449687.25 - 3194 iterations time 0.402, Presolve 0.02
Total time (CPU seconds):       0.58   (Wallclock seconds):       0.55


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


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.19s
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-rpu6u4od.lp -solve -solu /tmp/linopy-solve-alkng133.sol (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 2431 (-21397) rows, 4383 (-5557) columns and 15035 (-28483) elements
Perturbing problem by 0.001% of 2769.2736 - largest nonzero change 0.00099975992 ( 0.02882223%) - largest zero change 0.00099884126
0  Obj 81167.008 Primal inf 5431406.8 (2232) Dual inf 472.76596 (4)
123  Obj 11815.693 Primal inf 3794399 (2194)
246  Obj 11817.42 Primal inf 3176889.7 (2126)
369  Obj 11817.42 Primal inf 4339087.9 (2106)
492  Obj 11817.782 Primal inf 3625771.3 (2038)
615  Obj 11818.158 Primal inf 2690857.3 (1995)
738  Obj 11818.397 Primal inf 4137095.3 (1975)
861  Obj 11821.208 Primal inf 2769123.7 (1896)
984  Obj 11825.196 Primal inf 2932228.3 (1866)
1107  Obj 11826.847 Primal inf 4564624.5 (1822)
1230  Obj 11830.849 Primal inf 3532837 (1731)
1353  Obj 11833.334 Primal inf 4914663.9 (1685)
1476  Obj 11838.113 Primal inf 8475518.5 (1643)
1599  Obj 11843.166 Primal inf 4409780.7 (1558)
1722  Obj 11845.68 Primal inf 2968620.2 (1341)
1845  Obj 11852.629 Primal inf 7518450.5 (1386)
1968  Obj 11859.488 Primal inf 6960064.5 (1266)
2091  Obj 11866.636 Primal inf 35155163 (1139)
2214  Obj 30386.531 Primal inf 3610136.1 (940)
2337  Obj 37606.279 Primal inf 1.0669928e+09 (1110)
2460  Obj 130967.29 Primal inf 1.4942478e+09 (1027)
2583  Obj 270868.45 Primal inf 1.9206231e+09 (1257)
2703  Obj 386891.79 Primal inf 19893994 (871)
2826  Obj 505436.52 Primal inf 385009.24 (484)
2949  Obj 693856.75 Primal inf 42411.174 (194)
3072  Obj 807999.23 Primal inf 66051.783 (258)
3195  Obj 858129.85 Primal inf 7833.3099 (106)
3318  Obj 873099.28 Primal inf 17004.704 (63)
3414  Obj 873770.27
3414  Obj 873763.91 Dual inf 0.00073281 (1)
3415  Obj 873763.9
Optimal - objective value 873763.9
After Postsolve, objective 873763.9, infeasibilities - dual 78.034695 (3), primal 9.1550868e-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 - 3418 iterations time 0.412, Presolve 0.02
Total time (CPU seconds):       0.59   (Wallclock seconds):       0.57


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.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-j7xnnbsk.lp -solve -solu /tmp/linopy-solve-rcr5bob4.sol (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 2499 (-21329) rows, 4452 (-5488) columns and 15275 (-28243) elements
Perturbing problem by 0.001% of 2795.0927 - largest nonzero change 0.00099540096 ( 0.010060591%) - largest zero change 0.00099424254
0  Obj -38.285862 Primal inf 5618748.7 (2302)
124  Obj -38.285862 Primal inf 3609101.7 (2258)
248  Obj -36.867459 Primal inf 3211565.3 (2201)
372  Obj -36.81767 Primal inf 3103386.6 (2154)
496  Obj -36.338414 Primal inf 2905460.5 (2099)
620  Obj -35.757067 Primal inf 2745952.4 (2056)
744  Obj -35.378888 Primal inf 2784044.8 (1998)
868  Obj -34.134538 Primal inf 2547903 (1932)
992  Obj -33.288863 Primal inf 2683463.6 (1890)
1116  Obj -31.990764 Primal inf 5805981.5 (1821)
1240  Obj -30.479519 Primal inf 2833035.8 (1741)
1364  Obj -27.491745 Primal inf 3240761.7 (1642)
1488  Obj -24.754853 Primal inf 2095504.9 (1603)
1612  Obj -22.521896 Primal inf 4732434.1 (1604)
1736  Obj -18.959295 Primal inf 3720821 (1364)
1860  Obj 7646.1989 Primal inf 13648030 (1314)
1984  Obj 7649.7809 Primal inf 8384731.6 (1153)
2108  Obj 7654.593 Primal inf 8271871.8 (1159)
2232  Obj 7664.5417 Primal inf 504909.26 (804)
2356  Obj 16020.088 Primal inf 5978928.4 (921)
2480  Obj 47302.137 Primal inf 8087937.3 (833)
2604  Obj 217774.25 Primal inf 238410.55 (449)
2728  Obj 333106.02 Primal inf 219710.77 (431)
2852  Obj 456473.59 Primal inf 22945.499 (230)
2976  Obj 634442.4 Primal inf 31148.281 (238)
3100  Obj 714245.75 Primal inf 13823.277 (183)
3224  Obj 784534.43 Primal inf 329590.31 (143)
3308  Obj 790793.87
Optimal - objective value 790786.59
After Postsolve, objective 790786.59, infeasibilities - dual 0 (0), primal 0 (0)
Optimal objective 790786.5949 - 3308 iterations time 0.352, Presolve 0.02
Total time (CPU seconds):       0.53   (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.19s
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-1kig23f9.lp -solve -solu /tmp/linopy-solve-tvoq5618.sol (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 2482 (-21346) rows, 4436 (-5504) columns and 15226 (-28292) elements
Perturbing problem by 0.001% of 2769.2736 - largest nonzero change 0.00099734427 ( 0.007753724%) - largest zero change 0.00099365468
0  Obj 25309.268 Primal inf 5735747 (2299) Dual inf 149.8287 (1)
124  Obj 7971.4393 Primal inf 3653069.7 (2244)
248  Obj 7972.316 Primal inf 3685697.3 (2198)
372  Obj 7972.7107 Primal inf 2911053.5 (2126)
496  Obj 7973.3922 Primal inf 2856878.8 (2081)
620  Obj 7973.8379 Primal inf 2862427.3 (2034)
744  Obj 7974.1401 Primal inf 3864344.1 (2023)
868  Obj 7975.9135 Primal inf 3232534.1 (1883)
992  Obj 7979.4422 Primal inf 3432307.9 (1852)
1116  Obj 7981.8436 Primal inf 4273206.2 (1772)
1240  Obj 7984.6598 Primal inf 4946246.1 (1648)
1364  Obj 7987.511 Primal inf 5412392.3 (1633)
1488  Obj 12741.204 Primal inf 4365504.8 (1540)
1612  Obj 23725.565 Primal inf 12264186 (1501)
1669  Obj 23726.823 Primal inf 17920289 (1454)
1793  Obj 23731.116 Primal inf 6945993.2 (1408)
1917  Obj 23735.636 Primal inf 6909640.2 (1379)
2041  Obj 23741.706 Primal inf 17957968 (1298)
2165  Obj 33717.512 Primal inf 40930250 (1110)
2289  Obj 46914.214 Primal inf 54065477 (1321)
2413  Obj 104078.13 Primal inf 12144516 (891)
2537  Obj 210368.51 Primal inf 5861073.4 (922)
2661  Obj 449360.14 Primal inf 7438010.6 (854)
2785  Obj 800790.64 Primal inf 64226.422 (396)
2909  Obj 1193673.9 Primal inf 3462122.4 (324)
3033  Obj 1355883.1 Primal inf 7626.8779 (152)
3157  Obj 1434975 Primal inf 3305300.6 (716)
3281  Obj 1454202 Primal inf 1375.0903 (44)
3322  Obj 1455240.5
Optimal - objective value 1455232.7
After Postsolve, objective 1455232.7, infeasibilities - dual 24.980602 (1), primal 1.1527299e-07 (1)
Presolved model was optimal, full model needs cleaning up
0  Obj 1455232.7 Dual inf 0.43422624 (3)
End of values pass after 0 iterations
0  Obj 1455232.7 Dual inf 0.43422624 (3) w.o. free dual inf (1)
End of values pass after 3 iterations
3  Obj 1455232.7 Dual inf 0.18442032 (2) w.o. free dual inf (0)
5  Obj 1455232.7
Optimal - objective value 1455232.7
Optimal objective 1455232.678 - 3327 iterations time 0.412, Presolve 0.02
Total time (CPU seconds):       0.61   (Wallclock seconds):       0.59


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.19s
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-5if5x3a0.lp -solve -solu /tmp/linopy-solve-ha4vj_su.sol (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 2422 (-21406) rows, 4366 (-5574) columns and 15046 (-28472) elements
Perturbing problem by 0.001% of 2769.2736 - largest nonzero change 0.00099182877 ( 0.0098766581%) - largest zero change 0.0009900962
0  Obj 97177.543 Primal inf 6318925.8 (2257) Dual inf 599.31475 (4)
123  Obj 27997.772 Primal inf 3754035.3 (2215)
246  Obj 27997.772 Primal inf 3583083.9 (2159)
369  Obj 27997.974 Primal inf 3256263.9 (2098)
492  Obj 27998.597 Primal inf 4305727.6 (2065)
615  Obj 27999.325 Primal inf 2980855 (2024)
738  Obj 28000.508 Primal inf 2815268.2 (1993)
861  Obj 28006.467 Primal inf 2454282.6 (1902)
984  Obj 28008.917 Primal inf 3603377.1 (1821)
1107  Obj 28010.134 Primal inf 2982064 (1727)
1230  Obj 28012.738 Primal inf 3334159.5 (1611)
1353  Obj 28015.383 Primal inf 3242843.5 (1496)
1476  Obj 36006.515 Primal inf 1794829.6 (1398)
1599  Obj 59483.208 Primal inf 3740181.7 (1416)
1722  Obj 59488.767 Primal inf 4390182.5 (1347)
1845  Obj 59492.107 Primal inf 5529566.4 (1278)
1968  Obj 59499.188 Primal inf 4782456.8 (1308)
2091  Obj 74508.256 Primal inf 14280303 (1169)
2214  Obj 268704.65 Primal inf 10753783 (1156)
2337  Obj 406253.8 Primal inf 2563580.2 (874)
2460  Obj 710161.11 Primal inf 1506214.8 (711)
2583  Obj 1260266.5 Primal inf 103672.19 (486)
2706  Obj 2000844.3 Primal inf 57329.445 (375)
2829  Obj 2303407.5 Primal inf 25961.045 (261)
2952  Obj 2536889.5 Primal inf 130893.02 (365)
3075  Obj 2611240.3 Primal inf 530994.26 (539)
3198  Obj 2623411.9 Primal inf 15619.127 (161)
3321  Obj 2642989.8 Primal inf 624.12941 (76)
3444  Obj 2648920.7 Primal inf 0.2075864 (8)
3452  Obj 2649040.5
Optimal - objective value 2649026.4
After Postsolve, objective 2649026.4, infeasibilities - dual 217.25558 (5), primal 2.0607825e-06 (5)
Presolved model was optimal, full model needs cleaning up
0  Obj 2649026.4 Primal inf 6.3133686e-07 (1) Dual inf 1e+08 (6)
End of values pass after 10 iterations
10  Obj 2649026.4
Optimal - objective value 2649026.4
Optimal objective 2649026.448 - 3462 iterations time 0.432, Presolve 0.02
Total time (CPU seconds):       0.61   (Wallclock seconds):       0.59


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 2649026.44753696


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-iu0wqrim.lp -solve -solu /tmp/linopy-solve-znkam15t.sol (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 2429 (-21399) rows, 4376 (-5564) columns and 15067 (-28451) elements
Perturbing problem by 0.001% of 8902.5915 - largest nonzero change 0.00099776449 ( 0.0031988742%) - largest zero change 0.00099484211
0  Obj 92713.603 Primal inf 6151288.6 (2256) Dual inf 1195.3238 (4)
123  Obj 23362.744 Primal inf 3550110.3 (2192)
246  Obj 23362.797 Primal inf 3218252 (2144)
369  Obj 23362.9 Primal inf 3283046.2 (2115)
492  Obj 23362.996 Primal inf 4296370.9 (2067)
615  Obj 23363.391 Primal inf 5178020.5 (2024)
738  Obj 23364.132 Primal inf 3245899.1 (1972)
861  Obj 23365.11 Primal inf 2714154.8 (1908)
984  Obj 23368.058 Primal inf 2412970.8 (1838)
1107  Obj 23369.595 Primal inf 2135699.7 (1726)
1230  Obj 23370.398 Primal inf 3787984.2 (1672)
1353  Obj 23371.887 Primal inf 2032569.4 (1594)
1476  Obj 23373.474 Primal inf 24737767 (1574)
1599  Obj 40824.572 Primal inf 5904967.5 (1514)
1722  Obj 40826.691 Primal inf 4501877 (1357)
1845  Obj 40829.378 Primal inf 5517405.1 (1311)
1968  Obj 52309.305 Primal inf 6452240.9 (1352)
2091  Obj 52894.244 Primal inf 1828920.5 (965)
2214  Obj 91595.631 Primal inf 1134984.7 (844)
2337  Obj 173876.96 Primal inf 24742815 (826)
2460  Obj 394331.46 Primal inf 416334.15 (624)
2583  Obj 955853.76 Primal inf 102055.66 (542)
2706  Obj 1453829.7 Primal inf 45134.535 (398)
2829  Obj 1809321 Primal inf 19955.012 (264)
2952  Obj 1988861.8 Primal inf 15606.542 (236)
3075  Obj 2101524.6 Primal inf 1799.5422 (115)
3198  Obj 2134897.5 Primal inf 1515.2003 (85)
3321  Obj 2142508.2 Primal inf 384.81958 (47)
3370  Obj 2142961.6
Optimal - objective value 2142956.9
After Postsolve, objective 2142956.9, infeasibilities - dual 137.15559 (4), primal 2.8487449e-07 (4)
Presolved model was optimal, full model needs cleaning up
0  Obj 2142956.9 Dual inf 1.3715555 (4)
End of values pass after 5 iterations
5  Obj 2142956.9
Optimal - objective value 2142956.9
Optimal objective 2142956.889 - 3375 iterations time 0.382, Presolve 0.02
Total time (CPU seconds):       0.56   (Wallclock seconds):       0.54


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_4796/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.002552
std        0.259205
min       -0.700000
25%       -0.129689
50%        0.003142
75%        0.121985
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/style.py:76: 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/style.py:76: 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_4796/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_4796/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.038000 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.037277 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.037432 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.036777 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.036448 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.037601 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.036881 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.036727 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.036493 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.037139 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.037670 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.037955 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.037040 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.037649 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.036689 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.036672 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.036554 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.036552 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.036487 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.036815 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.036780 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.036859 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.036656 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.036650 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.000812
std        0.261525
min       -0.813751
25%       -0.125460
50%        0.003045
75%        0.123852
max        0.827584
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.020187
std        2.385689
min      -12.158272
25%       -0.463104
50%        0.001597
75%        0.535470
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/style.py:76: 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