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

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

Plot the distribution of the load and of generating tech

[3]:

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

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

/home/docs/checkouts/readthedocs.org/user_builds/pypsa/conda/latest/lib/python3.11/site-packages/cartopy/mpl/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')
<__array_function__ internals>:200: RuntimeWarning: invalid value encountered in cast
<__array_function__ internals>:200: RuntimeWarning: invalid value encountered in cast
<__array_function__ internals>:200: RuntimeWarning: invalid value encountered in cast
<__array_function__ internals>:200: RuntimeWarning: invalid value encountered in cast
INFO:linopy.model: Solve linear problem using Cbc solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 15/15 [00:00<00:00, 64.21it/s]
Writing variables.: 100%|██████████| 7/7 [00:00<00:00, 231.66it/s]
INFO:linopy.io: Writing time: 0.3s
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-m02q7lbe.lp -solve -solu /tmp/linopy-solve-6zjwgoop.sol (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 2370 (-21458) rows, 4321 (-5619) columns and 15455 (-28591) elements
Perturbing problem by 0.001% of 2258.0586 - largest nonzero change 0.00099976141 ( 0.017785398%) - largest zero change 0.00099947237
0  Obj 88132.25 Primal inf 5731630.7 (2165) Dual inf 303.14607 (4)
122  Obj 18780.935 Primal inf 3741280.4 (2105)
244  Obj 18782.011 Primal inf 3156892.8 (2057)
366  Obj 18782.868 Primal inf 2974452.7 (2015)
488  Obj 18783.553 Primal inf 2896611.5 (1961)
610  Obj 18785.856 Primal inf 2633354.5 (1910)
732  Obj 18786.768 Primal inf 5489151.4 (1846)
854  Obj 18789.822 Primal inf 2603103.1 (1768)
976  Obj 18793.855 Primal inf 3580883.3 (1724)
1098  Obj 18799.182 Primal inf 2346268.8 (1609)
1220  Obj 18801.571 Primal inf 3652423.9 (1578)
1342  Obj 18804.272 Primal inf 7973281.5 (1484)
1464  Obj 18806.306 Primal inf 8234356.8 (1440)
1586  Obj 22625.93 Primal inf 13245722 (1401)
1708  Obj 22630.526 Primal inf 6509204.9 (1335)
1830  Obj 22635.378 Primal inf 4610184.4 (1164)
1952  Obj 22641.296 Primal inf 9984859.6 (1126)
2074  Obj 24886.471 Primal inf 8314771.3 (1002)
2196  Obj 86626.796 Primal inf 1464418.1 (833)
2318  Obj 170935.97 Primal inf 15776717 (854)
2440  Obj 317367.81 Primal inf 63700418 (723)
2562  Obj 474031.76 Primal inf 945173.56 (440)
2684  Obj 972991.86 Primal inf 992999 (378)
2806  Obj 1276987.6 Primal inf 22645.097 (190)
2928  Obj 1385632.5 Primal inf 13127.236 (119)
3050  Obj 1437639.7 Primal inf 949476.01 (628)
3172  Obj 1449323 Primal inf 116.40819 (21)
3213  Obj 1449696.3
3213  Obj 1449687.3 Dual inf 0.00094475941 (3)
3217  Obj 1449687.3
Optimal - objective value 1449687.3
After Postsolve, objective 1449687.3, infeasibilities - dual 214.20597 (7), primal 1.7945027e-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 - 3227 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')
<__array_function__ internals>:200: RuntimeWarning: invalid value encountered in cast
<__array_function__ internals>:200: RuntimeWarning: invalid value encountered in cast
<__array_function__ internals>:200: RuntimeWarning: invalid value encountered in cast
<__array_function__ internals>:200: RuntimeWarning: invalid value encountered in cast
INFO:linopy.model: Solve linear problem using Cbc solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 15/15 [00:00<00:00, 65.18it/s]
Writing variables.: 100%|██████████| 7/7 [00:00<00:00, 224.03it/s]
INFO:linopy.io: Writing time: 0.3s
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



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

command line - cbc -printingOptions all -import /tmp/linopy-problem-cymqu1b4.lp -solve -solu /tmp/linopy-solve-ajf21uvo.sol (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 2391 (-21437) rows, 4359 (-5581) columns and 15576 (-28470) elements
Perturbing problem by 0.001% of 2620.1345 - largest nonzero change 0.00099700892 ( 0.019493681%) - largest zero change 0.00099607519
0  Obj 81156.815 Primal inf 5995243.1 (2192) Dual inf 447.30519 (4)
122  Obj 11805.501 Primal inf 4003831 (2158)
244  Obj 11806.142 Primal inf 3894941.3 (2105)
366  Obj 11808.343 Primal inf 3489517.9 (2054)
488  Obj 11808.697 Primal inf 2879336.3 (1987)
610  Obj 11810.936 Primal inf 2805000.5 (1953)
732  Obj 11812.831 Primal inf 8250220.7 (1906)
854  Obj 11815.308 Primal inf 9179258.7 (1826)
976  Obj 11819.869 Primal inf 9834441 (1778)
1098  Obj 11823.145 Primal inf 10259354 (1723)
1220  Obj 11827.588 Primal inf 9027524.3 (1692)
1342  Obj 11831.177 Primal inf 5661830.6 (1640)
1464  Obj 11836.33 Primal inf 8649644.1 (1515)
1586  Obj 11840.748 Primal inf 25990687 (1446)
1708  Obj 11846.02 Primal inf 23035619 (1364)
1830  Obj 11852.455 Primal inf 7786896.5 (1193)
1952  Obj 11857.707 Primal inf 5105664.6 (1059)
2074  Obj 12212.649 Primal inf 5.0970035e+08 (1027)
2196  Obj 13492.323 Primal inf 1.1483789e+09 (1014)
2318  Obj 54457.313 Primal inf 48601189 (840)
2440  Obj 224493.74 Primal inf 77655088 (845)
2562  Obj 377861.99 Primal inf 13558376 (663)
2684  Obj 499790.01 Primal inf 222833.36 (421)
2806  Obj 663977.22 Primal inf 409323.95 (329)
2928  Obj 813446.17 Primal inf 36340.277 (165)
3050  Obj 857175.94 Primal inf 171789.94 (412)
3172  Obj 873625.51 Primal inf 94.924055 (21)
3209  Obj 873771.66
Optimal - objective value 873763.9
After Postsolve, objective 873763.9, infeasibilities - dual 78.034695 (3), primal 9.1550873e-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 - 3212 iterations time 0.372, Presolve 0.02
Total time (CPU seconds):       0.58   (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')
<__array_function__ internals>:200: RuntimeWarning: invalid value encountered in cast
<__array_function__ internals>:200: RuntimeWarning: invalid value encountered in cast
<__array_function__ internals>:200: RuntimeWarning: invalid value encountered in cast
<__array_function__ internals>:200: RuntimeWarning: invalid value encountered in cast
INFO:linopy.model: Solve linear problem using Cbc solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 15/15 [00:00<00:00, 65.30it/s]
Writing variables.: 100%|██████████| 7/7 [00:00<00:00, 224.02it/s]
INFO:linopy.io: Writing time: 0.3s
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



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

command line - cbc -printingOptions all -import /tmp/linopy-problem-flw5nwkg.lp -solve -solu /tmp/linopy-solve-dfqqgq80.sol (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 2459 (-21369) rows, 4434 (-5506) columns and 15837 (-28209) elements
Perturbing problem by 0.001% of 2745.1675 - largest nonzero change 0.00099629266 ( 0.010269742%) - largest zero change 0.00099225248
0  Obj -46.556636 Primal inf 5896437.2 (2262)
124  Obj -46.556636 Primal inf 3846928.8 (2202)
248  Obj -44.82498 Primal inf 3891477.9 (2161)
372  Obj -44.156006 Primal inf 3309202.3 (2099)
496  Obj -43.381401 Primal inf 3305884.9 (2059)
620  Obj -41.881269 Primal inf 3984973.3 (2013)
744  Obj -40.263946 Primal inf 3632772.1 (1932)
868  Obj -38.4766 Primal inf 2692217.7 (1845)
992  Obj -36.318847 Primal inf 5247973.4 (1836)
1116  Obj -34.308108 Primal inf 3554862.8 (1797)
1240  Obj -29.841983 Primal inf 6654365.9 (1644)
1364  Obj -27.452571 Primal inf 10423780 (1593)
1488  Obj -23.449055 Primal inf 6977159.9 (1500)
1612  Obj -20.90136 Primal inf 4160718.4 (1331)
1736  Obj 7645.0605 Primal inf 25154037 (1372)
1860  Obj 7648.9785 Primal inf 33121681 (1293)
1984  Obj 7653.6105 Primal inf 36524583 (1193)
2108  Obj 7664.7036 Primal inf 3595244.1 (1058)
2232  Obj 7669.8619 Primal inf 7969548.7 (1105)
2356  Obj 18319.399 Primal inf 12583904 (868)
2480  Obj 156737.75 Primal inf 5.9686087e+08 (621)
2604  Obj 301824.11 Primal inf 7.4529458e+08 (539)
2728  Obj 465315.97 Primal inf 53702.188 (283)
2852  Obj 598987.02 Primal inf 148166.16 (392)
2976  Obj 704133.04 Primal inf 49089.848 (226)
3100  Obj 768065.66 Primal inf 251485.75 (227)
3224  Obj 785981.22 Primal inf 10111.895 (80)
3308  Obj 790794.89
Optimal - objective value 790786.59
After Postsolve, objective 790786.59, infeasibilities - dual 7.9999999 (1), primal 4.8500652e-07 (1)
Presolved model was optimal, full model needs cleaning up
0  Obj 790786.59 Dual inf 0.0799999 (1)
End of values pass after 1 iterations
1  Obj 790786.59
Optimal - objective value 790786.59
Optimal objective 790786.5949 - 3309 iterations time 0.412, Presolve 0.03
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')
<__array_function__ internals>:200: RuntimeWarning: invalid value encountered in cast
<__array_function__ internals>:200: RuntimeWarning: invalid value encountered in cast
<__array_function__ internals>:200: RuntimeWarning: invalid value encountered in cast
<__array_function__ internals>:200: RuntimeWarning: invalid value encountered in cast
INFO:linopy.model: Solve linear problem using Cbc solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 15/15 [00:00<00:00, 65.47it/s]
Writing variables.: 100%|██████████| 7/7 [00:00<00:00, 228.83it/s]
INFO:linopy.io: Writing time: 0.3s
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



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

command line - cbc -printingOptions all -import /tmp/linopy-problem-g7bybe5q.lp -solve -solu /tmp/linopy-solve-stn_524j.sol (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 2446 (-21382) rows, 4424 (-5516) columns and 15796 (-28250) elements
Perturbing problem by 0.001% of 2193.4858 - largest nonzero change 0.00099891412 ( 0.012052428%) - largest zero change 0.00099818337
0  Obj 25296.904 Primal inf 5950415.8 (2263) Dual inf 118.67629 (1)
123  Obj 7959.0755 Primal inf 3978176.8 (2212)
246  Obj 7960.1705 Primal inf 3508620.9 (2157)
369  Obj 7961.9396 Primal inf 3123948.8 (2093)
492  Obj 7962.9138 Primal inf 2798936.5 (2048)
615  Obj 7964.9705 Primal inf 3357865.7 (1994)
738  Obj 7966.6686 Primal inf 3127552.2 (1948)
861  Obj 7971.3784 Primal inf 2301315.8 (1848)
984  Obj 7974.9175 Primal inf 2372425.7 (1783)
1107  Obj 7977.7474 Primal inf 1934878.7 (1700)
1230  Obj 7981.5501 Primal inf 3937741.9 (1699)
1353  Obj 7985.3938 Primal inf 18263103 (1584)
1476  Obj 7987.4254 Primal inf 3476458.5 (1569)
1599  Obj 12741.997 Primal inf 4981282.2 (1441)
1722  Obj 23727.135 Primal inf 2844103.7 (1342)
1845  Obj 23732.575 Primal inf 14172222 (1288)
1968  Obj 23737.962 Primal inf 5646282.6 (1395)
2091  Obj 23743.25 Primal inf 6416597.8 (1202)
2214  Obj 23869.961 Primal inf 5135444 (943)
2337  Obj 84853.49 Primal inf 724643 (705)
2460  Obj 329083.93 Primal inf 30185874 (846)
2583  Obj 434682.64 Primal inf 297958.23 (508)
2706  Obj 800514.73 Primal inf 405043.6 (494)
2829  Obj 1045136.2 Primal inf 423910.25 (319)
2952  Obj 1336163.6 Primal inf 20098.938 (181)
3075  Obj 1429673 Primal inf 4333.2787 (114)
3198  Obj 1455135.9 Primal inf 14.494187 (17)
3215  Obj 1455244.8
Optimal - objective value 1455232.7
After Postsolve, objective 1455232.7, infeasibilities - dual 25.49824 (1), primal 7.6227303e-08 (1)
Presolved model was optimal, full model needs cleaning up
0  Obj 1455232.7 Dual inf 0.30197075 (2)
End of values pass after 0 iterations
0  Obj 1455232.7 Dual inf 0.30197075 (2) w.o. free dual inf (1)
End of values pass after 12 iterations
12  Obj 1455232.7 Dual inf 0.046988453 (1) w.o. free dual inf (0)
13  Obj 1455232.7
Optimal - objective value 1455232.7
Optimal objective 1455232.678 - 3228 iterations time 0.402, Presolve 0.02
Total time (CPU seconds):       0.60   (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')
<__array_function__ internals>:200: RuntimeWarning: invalid value encountered in cast
<__array_function__ internals>:200: RuntimeWarning: invalid value encountered in cast
<__array_function__ internals>:200: RuntimeWarning: invalid value encountered in cast
<__array_function__ internals>:200: RuntimeWarning: invalid value encountered in cast
INFO:linopy.model: Solve linear problem using Cbc solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 15/15 [00:00<00:00, 63.89it/s]
Writing variables.: 100%|██████████| 7/7 [00:00<00:00, 225.45it/s]
INFO:linopy.io: Writing time: 0.3s
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



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

command line - cbc -printingOptions all -import /tmp/linopy-problem-cbwtdjz4.lp -solve -solu /tmp/linopy-solve-ody6fmwr.sol (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 2390 (-21438) rows, 4356 (-5584) columns and 15619 (-28427) elements
Perturbing problem by 0.001% of 2193.4858 - largest nonzero change 0.00099820297 ( 0.010152757%) - largest zero change 0.00099818337
0  Obj 97166.872 Primal inf 5972833.9 (2225) Dual inf 474.70514 (4)
122  Obj 27987.1 Primal inf 4466877.7 (2181)
244  Obj 27987.999 Primal inf 3523912.1 (2116)
366  Obj 27987.999 Primal inf 3307009.2 (2089)
488  Obj 27990.787 Primal inf 3429125.9 (2060)
610  Obj 27993.493 Primal inf 3698479 (1990)
732  Obj 27994.983 Primal inf 6970738.1 (1906)
854  Obj 28002.991 Primal inf 3701372.4 (1812)
976  Obj 28007.381 Primal inf 3743336.9 (1727)
1098  Obj 28009.916 Primal inf 11038730 (1714)
1220  Obj 28012.41 Primal inf 1919136.3 (1604)
1342  Obj 28015.155 Primal inf 2008027.9 (1498)
1464  Obj 44626.704 Primal inf 2004446.4 (1435)
1586  Obj 59481.392 Primal inf 1887845.6 (1288)
1708  Obj 59491.691 Primal inf 7485056.2 (1319)
1830  Obj 66967.925 Primal inf 7564881.8 (1248)
1952  Obj 87265.368 Primal inf 13515530 (1176)
2074  Obj 123053.26 Primal inf 28281720 (1264)
2196  Obj 194531.43 Primal inf 40238455 (1202)
2318  Obj 405537.74 Primal inf 7459420.9 (839)
2440  Obj 511907.66 Primal inf 9970654.1 (860)
2562  Obj 1194775.8 Primal inf 691561.84 (554)
2684  Obj 1884492.2 Primal inf 50732.517 (355)
2806  Obj 2231687.5 Primal inf 251590.65 (594)
2928  Obj 2462901.8 Primal inf 20933.126 (236)
3050  Obj 2592817.7 Primal inf 39739.072 (243)
3172  Obj 2623840.1 Primal inf 11215.152 (173)
3294  Obj 2645834.3 Primal inf 2552.0989 (97)
3382  Obj 2649041.9
Optimal - objective value 2649026.4
After Postsolve, objective 2649026.4, infeasibilities - dual 124.09584 (3), primal 7.3815671e-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 12 iterations
12  Obj 2649026.4
Optimal - objective value 2649026.4
Optimal objective 2649026.448 - 3394 iterations time 0.452, Presolve 0.02
Total time (CPU seconds):       0.66   (Wallclock seconds):       0.62


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')
<__array_function__ internals>:200: RuntimeWarning: invalid value encountered in cast
<__array_function__ internals>:200: RuntimeWarning: invalid value encountered in cast
<__array_function__ internals>:200: RuntimeWarning: invalid value encountered in cast
<__array_function__ internals>:200: RuntimeWarning: invalid value encountered in cast
INFO:linopy.model: Solve linear problem using Cbc solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 15/15 [00:00<00:00, 65.11it/s]
Writing variables.: 100%|██████████| 7/7 [00:00<00:00, 227.35it/s]
INFO:linopy.io: Writing time: 0.3s
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



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

command line - cbc -printingOptions all -import /tmp/linopy-problem-uw6x393b.lp -solve -solu /tmp/linopy-solve-5aa4ajzr.sol (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 2391 (-21437) rows, 4374 (-5566) columns and 15653 (-28393) elements
Perturbing problem by 0.001% of 8902.5915 - largest nonzero change 0.00097303945 ( 0.0039378659%) - largest zero change 0.00097207615
0  Obj 92710.129 Primal inf 6096896.4 (2218) Dual inf 1195.3238 (4)
122  Obj 23359.141 Primal inf 3993516.6 (2165)
244  Obj 23359.367 Primal inf 3279783 (2093)
366  Obj 23359.993 Primal inf 3144630.8 (2089)
488  Obj 23360.199 Primal inf 3082814.7 (2047)
610  Obj 23362.366 Primal inf 2788197.4 (1967)
732  Obj 23363.072 Primal inf 4921548.2 (1903)
854  Obj 23364.579 Primal inf 2683838.3 (1804)
976  Obj 23367.964 Primal inf 2691822.5 (1765)
1098  Obj 23368.91 Primal inf 4423665.5 (1640)
1220  Obj 23369.979 Primal inf 2721819.4 (1653)
1342  Obj 23371.108 Primal inf 2933038.8 (1571)
1464  Obj 23373.102 Primal inf 2246191.3 (1420)
1586  Obj 40824.39 Primal inf 1982682.2 (1277)
1708  Obj 40826.919 Primal inf 3724465.8 (1218)
1830  Obj 40829.913 Primal inf 2262983 (1192)
1952  Obj 40833.174 Primal inf 8634958.7 (1156)
2074  Obj 76839.224 Primal inf 5894395.7 (995)
2196  Obj 112207.39 Primal inf 1180179 (832)
2318  Obj 348102.06 Primal inf 841771.56 (686)
2440  Obj 678666.25 Primal inf 2318933.6 (749)
2562  Obj 949060.21 Primal inf 331056.5 (767)
2684  Obj 1436171 Primal inf 65058.383 (287)
2806  Obj 1845544.9 Primal inf 14568.194 (237)
2928  Obj 2045415.3 Primal inf 4098.1153 (157)
3050  Obj 2123896.8 Primal inf 38108.989 (361)
3172  Obj 2142760.5 Primal inf 7.4289097 (19)
3193  Obj 2142962.6
Optimal - objective value 2142956.9
After Postsolve, objective 2142956.9, infeasibilities - dual 138.1249 (4), primal 5.5034699e-07 (4)
Presolved model was optimal, full model needs cleaning up
0  Obj 2142956.9 Dual inf 1.3812486 (4)
End of values pass after 5 iterations
5  Obj 2142956.9
Optimal - objective value 2142956.9
Optimal objective 2142956.889 - 3198 iterations time 0.362, Presolve 0.02
Total time (CPU seconds):       0.54   (Wallclock seconds):       0.52


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

[16]:

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

/home/docs/checkouts/readthedocs.org/user_builds/pypsa/conda/latest/lib/python3.11/site-packages/cartopy/mpl/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

[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.045870 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.047222 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.046564 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.045386 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.045401 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.045926 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.045623 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.045584 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.045581 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.046002 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.044911 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.044681 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.045854 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.044866 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.045277 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.045268 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.044901 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.045786 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.045293 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.045301 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.044935 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.044666 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.044946 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.044325 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.000449
std        0.260873
min       -0.813322
25%       -0.125479
50%        0.003032
75%        0.126691
max        0.827140
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.022763
std        2.373578
min      -12.158659
25%       -0.463000
50%        0.001597
75%        0.537443
max       17.959258
dtype: float64

Plot the reactive power

[29]:

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

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

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

/home/docs/checkouts/readthedocs.org/user_builds/pypsa/conda/latest/lib/python3.11/site-packages/cartopy/mpl/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