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://www.pypsa.org/examples/scigrid-sclopf.ipynb).

  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

/home/docs/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.24.0/lib/python3.11/site-packages/pypsa/networkclustering.py:16: UserWarning: The namespace `pypsa.networkclustering` is deprecated and will be removed in PyPSA v0.24. Please use `pypsa.clustering.spatial instead`.
  warnings.warn(

[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.24.0. Read the release notes at https://pypsa.readthedocs.io/en/latest/release_notes.html to prepare your network for import.
INFO:pypsa.io:Imported network scigrid-de.nc has buses, generators, lines, loads, storage_units, transformers

Plot the distribution of the load and of generating tech

[3]:

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

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

/home/docs/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.24.0/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/v0.24.0/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 linear problem using Cbc solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 15/15 [00:00<00:00, 78.23it/s]
Writing continuous variables.: 100%|██████████| 7/7 [00:00<00:00, 253.01it/s]
INFO:linopy.io: Writing time: 0.23s
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



Welcome to the CBC MILP Solver
Version: 2.10.10
Build Date: Apr 19 2023

command line - cbc -printingOptions all -import /tmp/linopy-problem-xy4o4n1h.lp -solve -solu /tmp/linopy-solve-hbb7ed7x.sol (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 2374 (-21454) rows, 4321 (-5619) columns and 16165 (-28653) elements
Perturbing problem by 0.001% of 2449.2813 - largest nonzero change 0.00099323312 ( 0.01516664%) - largest zero change 0.00099115148
0  Obj 88145.87 Primal inf 5088423.3 (2169) Dual inf 643.19931 (4)
122  Obj 18794.556 Primal inf 3786129.8 (2140)
244  Obj 18794.556 Primal inf 3187310.3 (2069)
366  Obj 18796.14 Primal inf 2963428.8 (2038)
488  Obj 18796.715 Primal inf 2846998.6 (1990)
610  Obj 18797.754 Primal inf 7577693.4 (1993)
732  Obj 18799.448 Primal inf 2414939.2 (1876)
854  Obj 18800.396 Primal inf 2301237.2 (1802)
976  Obj 18802.076 Primal inf 2760961.2 (1744)
1098  Obj 18805.201 Primal inf 4747823.5 (1712)
1220  Obj 18806.81 Primal inf 5421341.6 (1636)
1342  Obj 18808.738 Primal inf 4412286.2 (1529)
1464  Obj 18811.035 Primal inf 2928402.3 (1413)
1586  Obj 22629.283 Primal inf 1758589.6 (1313)
1708  Obj 22632.404 Primal inf 2583963.1 (1218)
1830  Obj 22636.13 Primal inf 3940344.6 (1228)
1952  Obj 22640.767 Primal inf 3402892.1 (1090)
2074  Obj 22646.236 Primal inf 1648889 (893)
2196  Obj 121662.16 Primal inf 72850654 (740)
2318  Obj 361205.88 Primal inf 93351879 (927)
2440  Obj 670115.98 Primal inf 56792064 (597)
2562  Obj 843636.42 Primal inf 1207774.3 (439)
2684  Obj 1104652.9 Primal inf 389290.34 (370)
2806  Obj 1315871.3 Primal inf 14459.471 (178)
2928  Obj 1427248.6 Primal inf 24800.465 (127)
3050  Obj 1449427.4 Primal inf 40.13038 (18)
3085  Obj 1449694.6
3085  Obj 1449687.3 Dual inf 4.902528e-05 (2)
3088  Obj 1449687.3
Optimal - objective value 1449687.3
After Postsolve, objective 1449687.3, infeasibilities - dual 214.20597 (7), primal 1.7945026e-06 (7)
Presolved model was optimal, full model needs cleaning up
0  Obj 1449687.3 Dual inf 2.142059 (7)
End of values pass after 10 iterations
10  Obj 1449687.3
Optimal - objective value 1449687.3
Optimal objective 1449687.25 - 3098 iterations time 0.362, Presolve 0.02
Total time (CPU seconds):       0.55   (Wallclock seconds):       0.52


INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Kirchhoff-Voltage-Law 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 linear problem using Cbc solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 15/15 [00:00<00:00, 81.25it/s]
Writing continuous variables.: 100%|██████████| 7/7 [00:00<00:00, 244.63it/s]
INFO:linopy.io: Writing time: 0.23s
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 Kirchhoff-Voltage-Law 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-2shmfqqt.lp -solve -solu /tmp/linopy-solve-y49a7o_f.sol (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 2398 (-21430) rows, 4362 (-5578) columns and 16308 (-28510) elements
Perturbing problem by 0.001% of 2998.2956 - largest nonzero change 0.00099191028 ( 0.065266124%) - largest zero change 0.00098701061
0  Obj 81182.754 Primal inf 4978100.6 (2198) Dual inf 1089.1409 (4)
122  Obj 11831.44 Primal inf 3766653.8 (2155)
244  Obj 11831.44 Primal inf 3315166.8 (2103)
366  Obj 11831.75 Primal inf 2953190.7 (2055)
488  Obj 11834.627 Primal inf 3630596.4 (1984)
610  Obj 11835.54 Primal inf 2772497.4 (1963)
732  Obj 11837.297 Primal inf 3102886 (1943)
854  Obj 11837.881 Primal inf 2862875.9 (1869)
976  Obj 11839.34 Primal inf 3424482.3 (1824)
1098  Obj 11840.686 Primal inf 4303892.4 (1757)
1220  Obj 11843.134 Primal inf 3418470.4 (1656)
1342  Obj 11844.848 Primal inf 2359939 (1563)
1464  Obj 11847.643 Primal inf 2459170.5 (1448)
1586  Obj 11849.765 Primal inf 12799719 (1613)
1708  Obj 11853.249 Primal inf 10977995 (1526)
1830  Obj 11857.014 Primal inf 3682235.8 (1130)
1952  Obj 11861.46 Primal inf 2417292.3 (1201)
2074  Obj 11867.139 Primal inf 1282648.6 (887)
2196  Obj 20240.185 Primal inf 733162.46 (846)
2318  Obj 44373.571 Primal inf 1605040.2 (946)
2440  Obj 142737.02 Primal inf 2789905.1 (827)
2562  Obj 267523.29 Primal inf 21580745 (768)
2684  Obj 521893.06 Primal inf 9963160 (594)
2806  Obj 591547.19 Primal inf 10931492 (670)
2928  Obj 759444.72 Primal inf 10131743 (330)
3050  Obj 861917.81 Primal inf 12375.022 (131)
3172  Obj 872990.74 Primal inf 9277.2043 (81)
3266  Obj 873772.06
3266  Obj 873763.94 Dual inf 0.00078640246 (2)
3270  Obj 873763.9
Optimal - objective value 873763.9
After Postsolve, objective 873763.9, infeasibilities - dual 78.034695 (3), primal 9.1550862e-07 (3)
Presolved model was optimal, full model needs cleaning up
0  Obj 873763.9 Dual inf 0.78034665 (3)
End of values pass after 3 iterations
3  Obj 873763.9
Optimal - objective value 873763.9
Optimal objective 873763.9023 - 3273 iterations time 0.392, Presolve 0.03
Total time (CPU seconds):       0.60   (Wallclock seconds):       0.56


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 linear problem using Cbc solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 15/15 [00:00<00:00, 81.36it/s]
Writing continuous variables.: 100%|██████████| 7/7 [00:00<00:00, 249.29it/s]
INFO:linopy.io: Writing time: 0.23s
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 Kirchhoff-Voltage-Law 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-96k08g19.lp -solve -solu /tmp/linopy-solve-q6xkhn5d.sol (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 2458 (-21370) rows, 4430 (-5510) columns and 16535 (-28283) elements
Perturbing problem by 0.001% of 3268.792 - largest nonzero change 0.00099338525 ( 0.0062578095%) - largest zero change 0.00099328771
0  Obj -38.474412 Primal inf 5062095.5 (2261)
124  Obj -38.474412 Primal inf 3828339 (2220)
248  Obj -38.474412 Primal inf 3255443.4 (2172)
372  Obj -38.474412 Primal inf 2989752.7 (2133)
496  Obj -35.917673 Primal inf 2728944.5 (2069)
620  Obj -34.814278 Primal inf 11035160 (2056)
744  Obj -32.995215 Primal inf 7743498.5 (1971)
868  Obj -30.804241 Primal inf 8833459.4 (1907)
992  Obj -29.505447 Primal inf 4656231.2 (1775)
1116  Obj -28.18612 Primal inf 3589112.2 (1721)
1240  Obj -25.849398 Primal inf 31918935 (1698)
1364  Obj -24.940983 Primal inf 4631345.2 (1589)
1488  Obj -22.399434 Primal inf 6124081.4 (1497)
1612  Obj -20.074864 Primal inf 6354371.5 (1461)
1736  Obj 7643.0721 Primal inf 5274102.3 (1343)
1860  Obj 7646.345 Primal inf 12587994 (1217)
1984  Obj 7651.8656 Primal inf 4876545.5 (1078)
2108  Obj 7655.8583 Primal inf 4777023.8 (849)
2232  Obj 7719.2573 Primal inf 3442500.1 (885)
2356  Obj 13803.676 Primal inf 1578739.1 (869)
2480  Obj 95115.32 Primal inf 91602478 (1044)
2604  Obj 226156.69 Primal inf 1.4602469e+08 (896)
2728  Obj 369817.26 Primal inf 54731030 (612)
2852  Obj 494810.15 Primal inf 10646240 (553)
2976  Obj 550387.11 Primal inf 16648939 (485)
3100  Obj 707664.52 Primal inf 40817.296 (237)
3224  Obj 776753.89 Primal inf 12005.736 (138)
3348  Obj 790192.89 Primal inf 9241.6779 (47)
3399  Obj 790796.05
Optimal - objective value 790786.59
After Postsolve, objective 790786.59, infeasibilities - dual 18 (2), primal 9.8396307e-07 (2)
Presolved model was optimal, full model needs cleaning up
0  Obj 790786.59 Dual inf 0.1799998 (2)
End of values pass after 3 iterations
3  Obj 790786.59
Optimal - objective value 790786.59
Optimal objective 790786.5949 - 3402 iterations time 0.422, Presolve 0.02
Total time (CPU seconds):       0.61   (Wallclock seconds):       0.58


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 linear problem using Cbc solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 15/15 [00:00<00:00, 78.35it/s]
Writing continuous variables.: 100%|██████████| 7/7 [00:00<00:00, 247.36it/s]
INFO:linopy.io: Writing time: 0.23s
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 Kirchhoff-Voltage-Law 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-63oz_t6e.lp -solve -solu /tmp/linopy-solve-b040spvh.sol (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 2447 (-21381) rows, 4425 (-5515) columns and 16507 (-28311) elements
Perturbing problem by 0.001% of 2240.9252 - largest nonzero change 0.00099874891 ( 0.0087681826%) - largest zero change 0.0009971789
0  Obj 25305.903 Primal inf 5111220.9 (2264) Dual inf 265.06173 (1)
123  Obj 7968.0747 Primal inf 3951358.7 (2228)
246  Obj 7968.0747 Primal inf 3356526.8 (2165)
369  Obj 7968.4997 Primal inf 3789404 (2127)
492  Obj 7968.6176 Primal inf 3334611.5 (2064)
615  Obj 7973.7355 Primal inf 3116964.2 (2021)
738  Obj 7977.4642 Primal inf 2795152.9 (1918)
861  Obj 7979.4215 Primal inf 2035235.1 (1835)
984  Obj 7980.1229 Primal inf 1887713.3 (1746)
1107  Obj 7981.5095 Primal inf 3443003.4 (1733)
1230  Obj 7983.8865 Primal inf 5200853.2 (1708)
1353  Obj 7985.7809 Primal inf 2187227.3 (1573)
1476  Obj 7988.0269 Primal inf 5442040.4 (1581)
1599  Obj 23724.589 Primal inf 5203808.8 (1402)
1722  Obj 23727.555 Primal inf 7811350.2 (1260)
1845  Obj 23730.856 Primal inf 43660730 (1187)
1968  Obj 23733.434 Primal inf 16678941 (1209)
2091  Obj 23737.19 Primal inf 9128159.7 (997)
2214  Obj 37302.141 Primal inf 7097554.6 (992)
2337  Obj 107236.12 Primal inf 28621079 (857)
2460  Obj 388222.81 Primal inf 1.0426889e+08 (669)
2583  Obj 767754.56 Primal inf 363722.22 (478)
2706  Obj 1135755.4 Primal inf 29360.373 (313)
2829  Obj 1324229.3 Primal inf 14064.752 (197)
2952  Obj 1432668.8 Primal inf 4720.5706 (115)
3075  Obj 1455242 Primal inf 0.38744656 (6)
3081  Obj 1455243.4
Optimal - objective value 1455232.7
After Postsolve, objective 1455232.7, infeasibilities - dual 24.951498 (1), primal 3.1439549e-08 (1)
Presolved model was optimal, full model needs cleaning up
0  Obj 1455232.7 Dual inf 0.24951488 (1)
End of values pass after 2 iterations
2  Obj 1455232.7
Optimal - objective value 1455232.7
Optimal objective 1455232.678 - 3083 iterations time 0.362, Presolve 0.03
Total time (CPU seconds):       0.55   (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 linear problem using Cbc solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 15/15 [00:00<00:00, 81.54it/s]
Writing continuous variables.: 100%|██████████| 7/7 [00:00<00:00, 249.67it/s]
INFO:linopy.io: Writing time: 0.23s
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 Kirchhoff-Voltage-Law 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-35exxh_9.lp -solve -solu /tmp/linopy-solve-cig84kne.sol (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 2394 (-21434) rows, 4363 (-5577) columns and 16355 (-28463) elements
Perturbing problem by 0.001% of 2240.9252 - largest nonzero change 0.00099928784 ( 0.0096128225%) - largest zero change 0.00099738508
0  Obj 97175.877 Primal inf 5191777.3 (2229) Dual inf 1060.2469 (4)
122  Obj 27996.104 Primal inf 3837028.4 (2185)
244  Obj 27996.104 Primal inf 3703917.1 (2140)
366  Obj 27996.104 Primal inf 3794398.3 (2104)
488  Obj 27996.886 Primal inf 3642563.8 (2069)
610  Obj 27997.418 Primal inf 4074601.5 (2030)
732  Obj 28002.037 Primal inf 2826288.1 (1959)
854  Obj 28004.075 Primal inf 3405656.8 (1876)
976  Obj 28007.342 Primal inf 5020913.6 (1812)
1098  Obj 28008.405 Primal inf 3962477.4 (1747)
1220  Obj 28011.557 Primal inf 17729335 (1668)
1342  Obj 28013.737 Primal inf 28156712 (1615)
1464  Obj 59475.288 Primal inf 40830396 (1518)
1586  Obj 59479.732 Primal inf 46727574 (1484)
1708  Obj 59484.206 Primal inf 50196311 (1341)
1830  Obj 59487.743 Primal inf 56271009 (1244)
1952  Obj 69450.068 Primal inf 8.6717468e+08 (1351)
2074  Obj 98637.852 Primal inf 8.7428539e+08 (1369)
2196  Obj 162067.47 Primal inf 5.0468364e+08 (1156)
2318  Obj 393247.1 Primal inf 5.0926048e+08 (1068)
2440  Obj 525582.07 Primal inf 6.8827228e+08 (1165)
2562  Obj 1056621.3 Primal inf 3.0632943e+08 (670)
2684  Obj 1503763 Primal inf 763236.57 (417)
2806  Obj 2113545.3 Primal inf 75204.548 (400)
2928  Obj 2389767.8 Primal inf 30784.149 (259)
3050  Obj 2532521.2 Primal inf 8989.4248 (175)
3172  Obj 2626120.3 Primal inf 5996.7813 (137)
3294  Obj 2648715.8 Primal inf 33.612269 (29)
3320  Obj 2649040.1
Optimal - objective value 2649026.4
After Postsolve, objective 2649026.4, infeasibilities - dual 124.09584 (3), primal 7.3815688e-07 (3)
Presolved model was optimal, full model needs cleaning up
0  Obj 2649026.4 Dual inf 1.2409581 (3)
End of values pass after 4 iterations
4  Obj 2649026.4
Optimal - objective value 2649026.4
Optimal objective 2649026.448 - 3324 iterations time 0.412, Presolve 0.02
Total time (CPU seconds):       0.61   (Wallclock seconds):       0.57


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 linear problem using Cbc solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 15/15 [00:00<00:00, 80.90it/s]
Writing continuous variables.: 100%|██████████| 7/7 [00:00<00:00, 251.83it/s]
INFO:linopy.io: Writing time: 0.23s
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 Kirchhoff-Voltage-Law 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-tly91uhn.lp -solve -solu /tmp/linopy-solve-3id285mx.sol (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 2400 (-21428) rows, 4361 (-5579) columns and 16367 (-28451) elements
Perturbing problem by 0.001% of 6880.5027 - largest nonzero change 0.00098724516 ( 0.0042006337%) - largest zero change 0.00097642845
0  Obj 92706.97 Primal inf 5178889.1 (2227) Dual inf 2583.3603 (4)
123  Obj 23356.094 Primal inf 3897806.3 (2184)
246  Obj 23356.094 Primal inf 3592974.9 (2157)
369  Obj 23356.094 Primal inf 3689364.7 (2105)
492  Obj 23356.617 Primal inf 3579962.8 (2056)
615  Obj 23357.231 Primal inf 5364625.7 (2035)
738  Obj 23360.785 Primal inf 9121264.1 (1988)
861  Obj 23362.861 Primal inf 3526642.3 (1865)
984  Obj 23364.439 Primal inf 5539635.6 (1846)
1107  Obj 23365.598 Primal inf 4273554.4 (1749)
1230  Obj 23366.735 Primal inf 3991583.5 (1608)
1353  Obj 23368.187 Primal inf 3555483.7 (1512)
1476  Obj 40819.194 Primal inf 16632663 (1505)
1599  Obj 40820.659 Primal inf 22498717 (1579)
1722  Obj 40823.007 Primal inf 23240243 (1637)
1845  Obj 40825.072 Primal inf 10875948 (1280)
1968  Obj 41009.423 Primal inf 2411531.3 (1124)
2091  Obj 41234.773 Primal inf 3180208.6 (1017)
2214  Obj 41237.188 Primal inf 1264776.6 (851)
2337  Obj 66109.249 Primal inf 5097110.8 (999)
2460  Obj 421634.1 Primal inf 4254496.5 (816)
2583  Obj 1036119.8 Primal inf 144810.93 (431)
2706  Obj 1690865.1 Primal inf 26537.325 (342)
2829  Obj 1920365.2 Primal inf 10655.101 (214)
2952  Obj 2082573.4 Primal inf 5030.9035 (170)
3075  Obj 2130701.3 Primal inf 1569.5023 (96)
3198  Obj 2142584.3 Primal inf 32.942598 (28)
3241  Obj 2142962.3
Optimal - objective value 2142956.9
After Postsolve, objective 2142956.9, infeasibilities - dual 176.69321 (5), primal 1.4341023e-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 6 iterations
6  Obj 2142956.9
Optimal - objective value 2142956.9
Optimal objective 2142956.889 - 3247 iterations time 0.382, Presolve 0.02
Total time (CPU seconds):       0.58   (Wallclock seconds):       0.55


[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.002892
std        0.258575
min       -0.700000
25%       -0.129689
50%        0.003084
75%        0.122643
max        0.700000
dtype: float64

[16]:

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

/home/docs/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.24.0/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/v0.24.0/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.044446 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.044481 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.043828 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.043849 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.044255 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.044086 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.044562 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.044068 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.044169 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.043838 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.044449 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.044546 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.044386 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.044942 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.043778 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.043801 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.044272 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.045208 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.044738 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.044319 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.044060 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.044168 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.044939 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.043816 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.000434
std        0.260926
min       -0.813358
25%       -0.125478
50%        0.003032
75%        0.126695
max        0.827178
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.022840
std        2.373703
min      -12.158621
25%       -0.463026
50%        0.001586
75%        0.537443
max       17.959258
dtype: float64

Plot the reactive power

[29]:

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

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

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

/home/docs/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.24.0/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