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
[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.25.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 objective.
Writing constraints.: 100%|██████████| 15/15 [00:00<00:00, 74.15it/s]
Writing continuous variables.: 100%|██████████| 7/7 [00:00<00:00, 200.28it/s]
INFO:linopy.io: Writing time: 0.25s
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.
Welcome to the CBC MILP Solver
Version: 2.10.10
Build Date: Apr 19 2023

command line - cbc -printingOptions all -import /tmp/linopy-problem-zjteiovs.lp -solve -solu /tmp/linopy-solve-bl4c0s4f.sol (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 2382 (-21446) rows, 4329 (-5611) columns and 16900 (-28614) elements
Perturbing problem by 0.001% of 2559.136 - largest nonzero change 0.00099541521 ( 0.015141654%) - largest zero change 0.00099500976
0  Obj 88148.835 Primal inf 5208378.9 (2177) Dual inf 479.04969 (4)
122  Obj 18797.52 Primal inf 4223375.2 (2150)
244  Obj 18797.52 Primal inf 3599932.2 (2102)
366  Obj 18797.52 Primal inf 3102746.1 (2056)
488  Obj 18797.52 Primal inf 2965133.6 (2003)
610  Obj 18799.174 Primal inf 3388932 (1973)
732  Obj 18801.261 Primal inf 3191502 (1880)
854  Obj 18802.116 Primal inf 3203737.1 (1833)
976  Obj 18803.079 Primal inf 4673642 (1761)
1098  Obj 18804.402 Primal inf 4043080.5 (1698)
1220  Obj 18807.651 Primal inf 15795526 (1698)
1342  Obj 18809.442 Primal inf 35418111 (1643)
1464  Obj 18812.921 Primal inf 4.306715e+08 (1679)
1586  Obj 22632.671 Primal inf 40483375 (1446)
1708  Obj 23932.133 Primal inf 16642510 (1382)
1830  Obj 24641.898 Primal inf 2.0220361e+08 (1561)
1952  Obj 24645.061 Primal inf 2.0148198e+08 (1286)
2074  Obj 24648.422 Primal inf 1.7698456e+08 (1116)
2136  Obj 42558.579 Primal inf 38573457 (1025)
2258  Obj 74927.9 Primal inf 39448559 (1282)
2380  Obj 189048.91 Primal inf 1.2709321e+08 (1138)
2502  Obj 364493.95 Primal inf 10368509 (684)
2624  Obj 621165.04 Primal inf 8382325.4 (558)
2746  Obj 977591.95 Primal inf 1.6783572e+08 (406)
2868  Obj 1259217.5 Primal inf 171016.93 (222)
2990  Obj 1391267.3 Primal inf 20681.372 (143)
3112  Obj 1439830.7 Primal inf 1584.6131 (58)
3197  Obj 1449695.5
3197  Obj 1449687.3 Dual inf 4.902528e-05 (2)
3200  Obj 1449687.3
Optimal - objective value 1449687.3
After Postsolve, objective 1449687.3, infeasibilities - dual 243.25253 (8), primal 1.9208569e-06 (8)
Presolved model was optimal, full model needs cleaning up
0  Obj 1449687.3 Primal inf 1.2635434e-07 (1) Dual inf 1e+08 (9)
End of values pass after 10 iterations
10  Obj 1449687.3
Optimal - objective value 1449687.3
Optimal objective 1449687.25 - 3210 iterations time 0.372, Presolve 0.02
Total time (CPU seconds):       0.57   (Wallclock seconds):       0.53

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 objective.
Writing constraints.: 100%|██████████| 15/15 [00:00<00:00, 75.81it/s]
Writing continuous variables.: 100%|██████████| 7/7 [00:00<00:00, 199.94it/s]
INFO:linopy.io: Writing time: 0.25s
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.
Welcome to the CBC MILP Solver
Version: 2.10.10
Build Date: Apr 19 2023

command line - cbc -printingOptions all -import /tmp/linopy-problem-66sa5hkd.lp -solve -solu /tmp/linopy-solve-kyql6x6b.sol (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 2402 (-21426) rows, 4364 (-5576) columns and 17020 (-28494) elements
Perturbing problem by 0.001% of 3241.309 - largest nonzero change 0.0009995126 ( 0.019358106%) - largest zero change 0.00099917067
0  Obj 81187.89 Primal inf 5093389.6 (2202) Dual inf 810.74934 (4)
123  Obj 11836.575 Primal inf 3989004.5 (2166)
246  Obj 11836.575 Primal inf 3564147 (2128)
369  Obj 11836.944 Primal inf 3217172.6 (2084)
492  Obj 11836.944 Primal inf 2991892.3 (2001)
615  Obj 11838.489 Primal inf 3247742.5 (1969)
738  Obj 11840.725 Primal inf 4001291.5 (1925)
861  Obj 11843.438 Primal inf 3013909.1 (1853)
984  Obj 11844.296 Primal inf 3029120.5 (1847)
1107  Obj 11845.694 Primal inf 3607572.1 (1710)
1230  Obj 11848.374 Primal inf 14429361 (1701)
1353  Obj 11851.101 Primal inf 2938562.5 (1531)
1476  Obj 11853.469 Primal inf 10789109 (1634)
1599  Obj 11855.547 Primal inf 10944208 (1584)
1722  Obj 11859.127 Primal inf 5639308.4 (1384)
1845  Obj 11862.844 Primal inf 4633652.9 (1400)
1968  Obj 11866.649 Primal inf 15384887 (1446)
2091  Obj 11871.071 Primal inf 45702441 (1473)
2214  Obj 11874.675 Primal inf 38412472 (1182)
2337  Obj 25783.045 Primal inf 7275034.4 (850)
2460  Obj 167549.02 Primal inf 4197040.4 (675)
2583  Obj 386921.9 Primal inf 13519416 (577)
2706  Obj 472955.25 Primal inf 3889122.4 (492)
2829  Obj 745213.27 Primal inf 363698.78 (236)
2952  Obj 835312.17 Primal inf 66600.393 (200)
3075  Obj 870693.36 Primal inf 5293.1068 (67)
3198  Obj 873766.43 Primal inf 11.783275 (7)
3212  Obj 873770.12
Optimal - objective value 873763.9
After Postsolve, objective 873763.9, infeasibilities - dual 78.034695 (3), primal 9.155089e-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 - 3215 iterations time 0.402, Presolve 0.02
Total time (CPU seconds):       0.57   (Wallclock seconds):       0.55

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 objective.
Writing constraints.: 100%|██████████| 15/15 [00:00<00:00, 75.45it/s]
Writing continuous variables.: 100%|██████████| 7/7 [00:00<00:00, 202.20it/s]
INFO:linopy.io: Writing time: 0.25s
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.
Welcome to the CBC MILP Solver
Version: 2.10.10
Build Date: Apr 19 2023

command line - cbc -printingOptions all -import /tmp/linopy-problem-wv1y8_1x.lp -solve -solu /tmp/linopy-solve-c9a7y_i5.sol (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 2462 (-21366) rows, 4429 (-5511) columns and 17238 (-28276) elements
Perturbing problem by 0.001% of 2559.136 - largest nonzero change 0.00099999381 ( 0.0083925711%) - largest zero change 0.00099994984
0  Obj -37.237207 Primal inf 5103732.3 (2265)
124  Obj -37.237207 Primal inf 3974916.8 (2221)
248  Obj -37.237207 Primal inf 3510635.1 (2195)
372  Obj -37.237207 Primal inf 3178253.9 (2138)
496  Obj -37.004607 Primal inf 2968157.5 (2089)
620  Obj -36.150684 Primal inf 3071626.1 (2045)
744  Obj -34.163641 Primal inf 2862331.6 (1980)
868  Obj -33.005464 Primal inf 5532250.5 (1945)
992  Obj -30.89407 Primal inf 8185840.3 (1931)
1116  Obj -28.84108 Primal inf 4026430.8 (1846)
1240  Obj -26.419961 Primal inf 5797766.5 (1765)
1364  Obj -23.616529 Primal inf 6801973.5 (1707)
1488  Obj -20.742335 Primal inf 10714438 (1627)
1612  Obj -17.78221 Primal inf 9666620.9 (1524)
1736  Obj 7645.7357 Primal inf 5111770.6 (1316)
1860  Obj 7649.1136 Primal inf 4147646.8 (1212)
1984  Obj 7652.7673 Primal inf 6877204.1 (1186)
2108  Obj 7664.5278 Primal inf 16606892 (981)
2232  Obj 7667.9284 Primal inf 7337855.8 (1063)
2356  Obj 61885.799 Primal inf 6920551 (938)
2480  Obj 127660.29 Primal inf 1814655.6 (743)
2604  Obj 333897.76 Primal inf 250486.28 (386)
2728  Obj 439888.09 Primal inf 367232.62 (439)
2852  Obj 651189.63 Primal inf 150584.94 (305)
2976  Obj 752893.13 Primal inf 67866.475 (173)
3100  Obj 786248.7 Primal inf 3730.9278 (71)
3206  Obj 790797.27
Optimal - objective value 790786.59
After Postsolve, objective 790786.59, infeasibilities - dual 0 (0), primal 0 (0)
Optimal objective 790786.5949 - 3206 iterations time 0.342, Presolve 0.02
Total time (CPU seconds):       0.55   (Wallclock seconds):       0.50

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 objective.
Writing constraints.: 100%|██████████| 15/15 [00:00<00:00, 74.46it/s]
Writing continuous variables.: 100%|██████████| 7/7 [00:00<00:00, 200.54it/s]
INFO:linopy.io: Writing time: 0.25s
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.
Welcome to the CBC MILP Solver
Version: 2.10.10
Build Date: Apr 19 2023

command line - cbc -printingOptions all -import /tmp/linopy-problem-y8fe28m5.lp -solve -solu /tmp/linopy-solve-kb6atxcr.sol (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 2451 (-21377) rows, 4440 (-5500) columns and 17226 (-28288) elements
Perturbing problem by 0.001% of 2559.136 - largest nonzero change 0.00099999381 ( 0.0080608671%) - largest zero change 0.00099994984
0  Obj 25313.96 Primal inf 5194791.2 (2268) Dual inf 206.2013 (1)
124  Obj 7976.1317 Primal inf 4101795.9 (2225)
248  Obj 7976.1317 Primal inf 3367222.8 (2180)
372  Obj 7976.1317 Primal inf 3782832.4 (2134)
496  Obj 7976.6954 Primal inf 5076042.7 (2101)
620  Obj 7979.0837 Primal inf 3450981.8 (2070)
744  Obj 7981.1975 Primal inf 3101908.5 (1982)
868  Obj 7982.825 Primal inf 2798363.8 (1892)
992  Obj 7983.9478 Primal inf 3221266.4 (1820)
1116  Obj 7986.2941 Primal inf 3145952.3 (1808)
1240  Obj 7989.2541 Primal inf 2522541.3 (1688)
1364  Obj 7992.0921 Primal inf 4986757.3 (1655)
1488  Obj 7994.6565 Primal inf 16769866 (1754)
1612  Obj 12749.202 Primal inf 32448990 (1460)
1736  Obj 23734.329 Primal inf 4244445.1 (1378)
1860  Obj 23738.033 Primal inf 10174688 (1272)
1984  Obj 23741.412 Primal inf 3913900.2 (1055)
2108  Obj 23777.471 Primal inf 14186809 (1226)
2232  Obj 101826.77 Primal inf 3527449.3 (845)
2356  Obj 162533.65 Primal inf 2945595.9 (640)
2480  Obj 285275.01 Primal inf 396002.45 (563)
2604  Obj 748236.26 Primal inf 249039.81 (508)
2728  Obj 1092656.3 Primal inf 87554.249 (341)
2852  Obj 1290628.8 Primal inf 225146.24 (341)
2976  Obj 1406589.6 Primal inf 157224.3 (230)
3100  Obj 1442647.7 Primal inf 4675.3356 (109)
3223  Obj 1455244.6
Optimal - objective value 1455232.7
After Postsolve, objective 1455232.7, infeasibilities - dual 0 (0), primal 0 (0)
Optimal objective 1455232.678 - 3223 iterations time 0.362, Presolve 0.02
Total time (CPU seconds):       0.56   (Wallclock seconds):       0.52

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 objective.
Writing constraints.: 100%|██████████| 15/15 [00:00<00:00, 74.07it/s]
Writing continuous variables.: 100%|██████████| 7/7 [00:00<00:00, 204.26it/s]
INFO:linopy.io: Writing time: 0.25s
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.
Welcome to the CBC MILP Solver
Version: 2.10.10
Build Date: Apr 19 2023

command line - cbc -printingOptions all -import /tmp/linopy-problem-5yoaxfsd.lp -solve -solu /tmp/linopy-solve-a67z6e81.sol (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 2398 (-21430) rows, 4374 (-5566) columns and 17074 (-28440) elements
Perturbing problem by 0.001% of 2559.136 - largest nonzero change 0.00099747629 ( 0.0087001389%) - largest zero change 0.0009969372
0  Obj 97182.151 Primal inf 5533521.7 (2233) Dual inf 824.80515 (4)
122  Obj 28002.378 Primal inf 4044220.3 (2183)
244  Obj 28002.378 Primal inf 3635437.6 (2153)
366  Obj 28002.378 Primal inf 3752879.9 (2091)
488  Obj 28003.75 Primal inf 3823462.7 (2064)
610  Obj 28004.776 Primal inf 5559411.6 (2024)
732  Obj 28007.205 Primal inf 3474532.5 (1958)
854  Obj 28009.005 Primal inf 5010859.5 (1929)
976  Obj 28011.279 Primal inf 4358441.6 (1826)
1098  Obj 28014.56 Primal inf 4835711 (1806)
1220  Obj 28017.865 Primal inf 17237744 (1691)
1342  Obj 28019.43 Primal inf 45234292 (1676)
1464  Obj 59482.687 Primal inf 11263859 (1557)
1586  Obj 59485.29 Primal inf 9143920.3 (1350)
1708  Obj 59489.767 Primal inf 3565540.2 (1231)
1830  Obj 59493.366 Primal inf 6990649.8 (1207)
1952  Obj 87201.54 Primal inf 13765050 (1185)
2074  Obj 128086.18 Primal inf 5380793.3 (1079)
2196  Obj 277662.31 Primal inf 1475091.7 (895)
2318  Obj 410445.7 Primal inf 845555.69 (671)
2440  Obj 1144207.4 Primal inf 541350.92 (572)
2562  Obj 1538189 Primal inf 6750439.9 (556)
2684  Obj 1772545.6 Primal inf 537571.34 (637)
2806  Obj 2453761.4 Primal inf 141258.6 (390)
2928  Obj 2583267.7 Primal inf 4174.8612 (142)
3050  Obj 2636197.5 Primal inf 9505.3827 (109)
3172  Obj 2648912.7 Primal inf 0.21975397 (8)
3180  Obj 2649040.4
Optimal - objective value 2649026.4
After Postsolve, objective 2649026.4, infeasibilities - dual 258.5252 (6), primal 2.3582362e-06 (6)
Presolved model was optimal, full model needs cleaning up
0  Obj 2649026.4 Primal inf 6.3133674e-07 (1) Dual inf 1e+08 (7)
End of values pass after 7 iterations
7  Obj 2649026.4
Optimal - objective value 2649026.4
Optimal objective 2649026.448 - 3187 iterations time 0.382, Presolve 0.02
Total time (CPU seconds):       0.56   (Wallclock seconds):       0.54

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 objective.
Writing constraints.: 100%|██████████| 15/15 [00:00<00:00, 74.39it/s]
Writing continuous variables.: 100%|██████████| 7/7 [00:00<00:00, 200.96it/s]
INFO:linopy.io: Writing time: 0.25s
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.
Welcome to the CBC MILP Solver
Version: 2.10.10
Build Date: Apr 19 2023

command line - cbc -printingOptions all -import /tmp/linopy-problem-5mt6sddy.lp -solve -solu /tmp/linopy-solve-t004bhu5.sol (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 2405 (-21423) rows, 4377 (-5563) columns and 17094 (-28420) elements
Perturbing problem by 0.001% of 10344.616 - largest nonzero change 0.00099119089 ( 0.0040686248%) - largest zero change 0.0009891535
0  Obj 92714.292 Primal inf 5404576 (2232) Dual inf 1936.5966 (4)
123  Obj 23363.44 Primal inf 4037887.6 (2184)
246  Obj 23363.44 Primal inf 3854125.5 (2158)
369  Obj 23363.44 Primal inf 5330363.4 (2118)
492  Obj 23363.615 Primal inf 3607375.8 (2065)
615  Obj 23363.94 Primal inf 3767050.8 (2009)
738  Obj 23365.621 Primal inf 3374128.2 (1938)
861  Obj 23366.56 Primal inf 4694872 (1906)
984  Obj 23367.668 Primal inf 13193522 (1854)
1107  Obj 23368.394 Primal inf 23152034 (1792)
1230  Obj 23369.556 Primal inf 7427100.1 (1650)
1353  Obj 23370.584 Primal inf 7279777.8 (1605)
1476  Obj 23371.877 Primal inf 2226961.6 (1447)
1599  Obj 40822.46 Primal inf 2361538.8 (1379)
1722  Obj 40823.944 Primal inf 1881779.5 (1259)
1845  Obj 40826.325 Primal inf 1592837.2 (1219)
1968  Obj 40828.268 Primal inf 1224388.3 (1068)
2091  Obj 40830.688 Primal inf 1417660.7 (958)
2214  Obj 41015.479 Primal inf 1044677.4 (911)
2337  Obj 229375.22 Primal inf 1137201.9 (860)
2460  Obj 534844.03 Primal inf 822495.53 (711)
2583  Obj 844177.62 Primal inf 135806.31 (469)
2706  Obj 1485379.9 Primal inf 37075.653 (287)
2829  Obj 1924833 Primal inf 10962.854 (207)
2952  Obj 2075480.1 Primal inf 3313.2508 (126)
3075  Obj 2133993.1 Primal inf 280.0322 (62)
3177  Obj 2142962.4
Optimal - objective value 2142956.9
After Postsolve, objective 2142956.9, infeasibilities - dual 176.69321 (5), primal 1.4341022e-06 (5)
Presolved model was optimal, full model needs cleaning up
0  Obj 2142956.9 Dual inf 1.7669316 (5)
End of values pass after 9 iterations
9  Obj 2142956.9
Optimal - objective value 2142956.9
Optimal objective 2142956.889 - 3186 iterations time 0.352, Presolve 0.02
Total time (CPU seconds):       0.53   (Wallclock seconds):       0.51

[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
[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.002503
std        0.259138
min       -0.700000
25%       -0.127897
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/numpy/core/fromnumeric.py:3504: RuntimeWarning: Mean of empty slice.
  return _methods._mean(a, axis=axis, dtype=dtype,
/home/docs/checkouts/readthedocs.org/user_builds/pypsa/conda/latest/lib/python3.11/site-packages/numpy/core/_methods.py:129: RuntimeWarning: invalid value encountered in scalar divide
  ret = ret.dtype.type(ret / rcount)
/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
[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
[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.044525 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.141684 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.043783 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.043761 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.045038 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.044794 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.043262 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.043236 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.044595 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.043636 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.044842 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.044736 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.044647 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.044950 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.043842 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.044403 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.045235 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.044824 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.045365 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.045097 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.045618 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.044458 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.043844 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.043637 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.000861
std        0.261455
min       -0.813714
25%       -0.123196
50%        0.003046
75%        0.123850
max        0.827546
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.020094
std        2.385565
min      -12.158311
25%       -0.463078
50%        0.001626
75%        0.535464
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