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/v0.23.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.23.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')
<__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, 62.85it/s]
Writing variables.: 100%|██████████| 7/7 [00:00<00:00, 221.37it/s]
INFO:linopy.io: Writing time: 0.31s
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-ps78_0si.lp -solve -solu /tmp/linopy-solve-ock283r4.sol (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 2380 (-21448) rows, 4329 (-5611) columns and 15526 (-28700) elements
Perturbing problem by 0.001% of 2240.9252 - largest nonzero change 0.00099749507 ( 0.014846054%) - largest zero change 0.00099738508
0  Obj 88144.088 Primal inf 5298659.3 (2175) Dual inf 553.3846 (4)
122  Obj 18792.774 Primal inf 3913238.1 (2148)
244  Obj 18792.774 Primal inf 3676039.4 (2091)
366  Obj 18792.774 Primal inf 3021292.7 (2032)
488  Obj 18797.516 Primal inf 2662834.2 (1978)
610  Obj 18798.21 Primal inf 4473914.7 (1977)
732  Obj 18799.54 Primal inf 3264547.1 (1945)
854  Obj 18802.341 Primal inf 3442396 (1885)
976  Obj 18803.112 Primal inf 3907815.5 (1850)
1098  Obj 18806.652 Primal inf 68739944 (1740)
1220  Obj 18808.249 Primal inf 81732895 (1672)
1342  Obj 18810.458 Primal inf 35392259 (1545)
1464  Obj 22628.294 Primal inf 34747600 (1439)
1586  Obj 22631.335 Primal inf 39131896 (1369)
1708  Obj 22634.369 Primal inf 6592748.5 (1265)
1830  Obj 22638.666 Primal inf 3108955.9 (1163)
1952  Obj 22642.233 Primal inf 4769544.2 (1090)
2074  Obj 22681.019 Primal inf 11952750 (1127)
2196  Obj 173937.66 Primal inf 12025391 (1016)
2318  Obj 218125.34 Primal inf 4431494.9 (941)
2440  Obj 440392.46 Primal inf 1964679.1 (684)
2562  Obj 789101.36 Primal inf 295127.86 (595)
2684  Obj 1127777.1 Primal inf 37967.095 (281)
2806  Obj 1311398.2 Primal inf 11391.813 (175)
2928  Obj 1422205.5 Primal inf 3118.6274 (86)
3050  Obj 1449121.7 Primal inf 168.50409 (27)
3101  Obj 1449695.5
3101  Obj 1449687.3 Dual inf 0.00089529348 (2)
3104  Obj 1449687.3
Optimal - objective value 1449687.3
After Postsolve, objective 1449687.3, infeasibilities - dual 186.92063 (6), primal 1.3870423e-06 (6)
Presolved model was optimal, full model needs cleaning up
0  Obj 1449687.3 Dual inf 12.306893 (14)
End of values pass after 0 iterations
0  Obj 1449687.3 Dual inf 12.306893 (14) w.o. free dual inf (12)
End of values pass after 8 iterations
8  Obj 1449687.3 Dual inf 10.437687 (8) w.o. free dual inf (6)
10  Obj 1449687.3
Optimal - objective value 1449687.3
Optimal objective 1449687.25 - 3114 iterations time 0.382, Presolve 0.02
Total time (CPU seconds):       0.58   (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')
<__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.73it/s]
Writing variables.: 100%|██████████| 7/7 [00:00<00:00, 226.55it/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-h4deea38.lp -solve -solu /tmp/linopy-solve-83ijgo5z.sol (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 2403 (-21425) rows, 4362 (-5578) columns and 15744 (-28482) elements
Perturbing problem by 0.001% of 3070.478 - largest nonzero change 0.00099870555 ( 0.015179458%) - largest zero change 0.00099582028
0  Obj 81183.874 Primal inf 5108274.1 (2203) Dual inf 1043.2235 (4)
123  Obj 11832.559 Primal inf 3658988.9 (2161)
246  Obj 11832.559 Primal inf 3498053.5 (2103)
369  Obj 11833.961 Primal inf 3111221.3 (2076)
492  Obj 11834.125 Primal inf 3090398 (2006)
615  Obj 11835.057 Primal inf 3097681.1 (1940)
738  Obj 11837.729 Primal inf 2600358.1 (1878)
861  Obj 11840.049 Primal inf 2533201.7 (1821)
984  Obj 11842.114 Primal inf 22123968 (1784)
1107  Obj 11842.555 Primal inf 20665844 (1739)
1230  Obj 11845.314 Primal inf 13081951 (1656)
1353  Obj 11847.053 Primal inf 21997502 (1583)
1476  Obj 11849.762 Primal inf 24459627 (1476)
1599  Obj 11851.969 Primal inf 22285329 (1502)
1722  Obj 11855.302 Primal inf 28615222 (1437)
1845  Obj 11859.923 Primal inf 45814228 (1454)
1968  Obj 11865.676 Primal inf 15650037 (1015)
2091  Obj 26245.831 Primal inf 8519288.5 (1246)
2214  Obj 93471.772 Primal inf 7248911.4 (993)
2337  Obj 119576.58 Primal inf 1127730.5 (699)
2460  Obj 277547.31 Primal inf 2825341.6 (653)
2583  Obj 391310.18 Primal inf 13957155 (743)
2706  Obj 587466.74 Primal inf 69440.057 (310)
2829  Obj 746208.23 Primal inf 1250812.7 (445)
2952  Obj 834150.05 Primal inf 13948.832 (146)
3075  Obj 867929.52 Primal inf 1071.5131 (62)
3168  Obj 873772.41
3168  Obj 873763.9 Dual inf 0.0007078415 (1)
3169  Obj 873763.9
Optimal - objective value 873763.9
After Postsolve, objective 873763.9, infeasibilities - dual 78.034695 (3), primal 9.1550879e-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 - 3172 iterations time 0.392, Presolve 0.02
Total time (CPU seconds):       0.59   (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')
<__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.09it/s]
Writing variables.: 100%|██████████| 7/7 [00:00<00:00, 225.13it/s]
INFO:linopy.io: Writing time: 0.31s
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-6ah_fy3j.lp -solve -solu /tmp/linopy-solve-c9k12ini.sol (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 2466 (-21362) rows, 4425 (-5515) columns and 15973 (-28253) elements
Perturbing problem by 0.001% of 2240.9252 - largest nonzero change 0.00099117356 ( 0.0093201145%) - largest zero change 0.00098986715
0  Obj -44.584898 Primal inf 5220251.7 (2269)
124  Obj -44.584898 Primal inf 4103124.4 (2238)
248  Obj -44.584898 Primal inf 3420602.7 (2174)
372  Obj -40.256823 Primal inf 2983360.6 (2126)
496  Obj -37.582136 Primal inf 2700106.3 (2077)
620  Obj -36.185086 Primal inf 8567154 (2067)
744  Obj -34.92718 Primal inf 5142053.5 (2039)
868  Obj -32.380491 Primal inf 5594235.8 (1921)
992  Obj -30.416509 Primal inf 7927578.8 (1851)
1116  Obj -28.638856 Primal inf 12527217 (1763)
1240  Obj -26.146098 Primal inf 4812276.2 (1650)
1364  Obj -24.434687 Primal inf 9264882.4 (1604)
1488  Obj -22.005948 Primal inf 6585190.1 (1630)
1612  Obj -18.938605 Primal inf 12208133 (1496)
1736  Obj 7644.8319 Primal inf 7418446.4 (1378)
1860  Obj 7648.4993 Primal inf 2502777.6 (1162)
1984  Obj 7652.1304 Primal inf 1605687.6 (1032)
2108  Obj 7657.2909 Primal inf 2141736.5 (975)
2232  Obj 7727.0947 Primal inf 943120.44 (794)
2356  Obj 7786.6143 Primal inf 5789342.8 (969)
2480  Obj 82108.081 Primal inf 3747810.8 (740)
2604  Obj 251407.47 Primal inf 3237928.5 (660)
2728  Obj 393492.32 Primal inf 1255745.8 (538)
2852  Obj 605162.73 Primal inf 100552.45 (258)
2976  Obj 703797.58 Primal inf 10363.57 (169)
3100  Obj 779663.59 Primal inf 1661.7557 (82)
3224  Obj 790554.19 Primal inf 10188.046 (58)
3273  Obj 790797.17
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 - 3274 iterations time 0.402, Presolve 0.02
Total time (CPU seconds):       0.59   (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')
<__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, 62.64it/s]
Writing variables.: 100%|██████████| 7/7 [00:00<00:00, 223.73it/s]
INFO:linopy.io: Writing time: 0.31s
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-mdnoq1x6.lp -solve -solu /tmp/linopy-solve-l9t_va9l.sol (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 2450 (-21378) rows, 4425 (-5515) columns and 15937 (-28289) elements
Perturbing problem by 0.001% of 2240.9252 - largest nonzero change 0.00099880288 ( 0.010047233%) - largest zero change 0.0009976219
0  Obj 25308.261 Primal inf 5265409 (2267) Dual inf 247.5123 (1)
124  Obj 7970.4323 Primal inf 4163908.6 (2235)
248  Obj 7970.4323 Primal inf 3942473.4 (2187)
372  Obj 7970.4323 Primal inf 5913333.2 (2175)
496  Obj 7974.7219 Primal inf 3448768.4 (2098)
620  Obj 7978.9685 Primal inf 4425602.1 (2075)
744  Obj 7980.7471 Primal inf 3741758.5 (1998)
868  Obj 7983.9775 Primal inf 3025995.6 (1896)
992  Obj 7985.6157 Primal inf 3509923.7 (1836)
1116  Obj 7986.4856 Primal inf 2945161.6 (1764)
1240  Obj 7988.5425 Primal inf 5032595.6 (1744)
1364  Obj 7990.8858 Primal inf 4881470.3 (1654)
1488  Obj 12744.434 Primal inf 7775493.6 (1629)
1612  Obj 12746.795 Primal inf 4950056.1 (1493)
1736  Obj 23731.57 Primal inf 3495643.6 (1406)
1860  Obj 23735.001 Primal inf 36452432 (1403)
1984  Obj 23740.247 Primal inf 12407112 (1229)
2108  Obj 31295.75 Primal inf 1490729.7 (976)
2232  Obj 79549.243 Primal inf 1128898 (697)
2356  Obj 117110.24 Primal inf 5168696.9 (807)
2480  Obj 341872.75 Primal inf 1058342.1 (631)
2604  Obj 640245.65 Primal inf 1721167.3 (721)
2728  Obj 1010085.2 Primal inf 52753.257 (351)
2852  Obj 1262400.9 Primal inf 17774.604 (252)
2976  Obj 1402538.6 Primal inf 6379.0641 (139)
3100  Obj 1440540.1 Primal inf 2585.3527 (72)
3224  Obj 1455068.7 Primal inf 20.056235 (18)
3241  Obj 1455243.6
Optimal - objective value 1455232.7
After Postsolve, objective 1455232.7, infeasibilities - dual 24.259813 (1), primal 6.0436272e-07 (1)
Presolved model was optimal, full model needs cleaning up
0  Obj 1455232.7 Dual inf 0.29007703 (2)
End of values pass after 0 iterations
0  Obj 1455232.7 Dual inf 0.29007703 (2) w.o. free dual inf (1)
End of values pass after 1 iterations
1  Obj 1455232.7 Dual inf 0.047479 (1) w.o. free dual inf (0)
2  Obj 1455232.7
Optimal - objective value 1455232.7
Optimal objective 1455232.678 - 3243 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, 64.03it/s]
Writing variables.: 100%|██████████| 7/7 [00:00<00:00, 222.71it/s]
INFO:linopy.io: Writing time: 0.31s
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-dx218r9d.lp -solve -solu /tmp/linopy-solve-553vso5z.sol (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 2397 (-21431) rows, 4370 (-5570) columns and 15806 (-28420) elements
Perturbing problem by 0.001% of 2240.9252 - largest nonzero change 0.00099501233 ( 0.0099550754%) - largest zero change 0.00099429059
0  Obj 97176.272 Primal inf 5267954.4 (2232) Dual inf 990.04917 (4)
122  Obj 27996.499 Primal inf 3893927.7 (2188)
244  Obj 27996.499 Primal inf 3873699.5 (2149)
366  Obj 27996.499 Primal inf 4539384.7 (2120)
488  Obj 27997.447 Primal inf 5819207.4 (2077)
610  Obj 28006.083 Primal inf 4599924.5 (2004)
732  Obj 28008.02 Primal inf 3608164.8 (1921)
854  Obj 28009.228 Primal inf 4074478.4 (1867)
976  Obj 28009.714 Primal inf 3467745.5 (1787)
1098  Obj 28011.679 Primal inf 2952120.3 (1688)
1220  Obj 28014.085 Primal inf 4444161.8 (1658)
1342  Obj 28017.186 Primal inf 5525040.9 (1584)
1464  Obj 59480.236 Primal inf 2109331.4 (1413)
1586  Obj 59484.281 Primal inf 2296653 (1324)
1708  Obj 59488.006 Primal inf 2809372.2 (1280)
1830  Obj 59492.107 Primal inf 2808077.2 (1187)
1952  Obj 59497.442 Primal inf 3086809.6 (1044)
2074  Obj 73548.109 Primal inf 42404395 (1341)
2196  Obj 142289.61 Primal inf 83607795 (1386)
2318  Obj 467314.82 Primal inf 2.3388705e+08 (1133)
2440  Obj 1022511.1 Primal inf 14507336 (737)
2562  Obj 1523292.2 Primal inf 14426946 (567)
2684  Obj 1970336.7 Primal inf 323688.32 (324)
2806  Obj 2339169 Primal inf 37809.416 (292)
2928  Obj 2516521.7 Primal inf 17725.024 (193)
3050  Obj 2617878.7 Primal inf 46547.941 (224)
3172  Obj 2643414.1 Primal inf 1604.7586 (70)
3252  Obj 2649046.2
Optimal - objective value 2649026.4
After Postsolve, objective 2649026.4, infeasibilities - dual 165.36547 (4), primal 1.0356107e-06 (4)
Presolved model was optimal, full model needs cleaning up
0  Obj 2649026.4 Dual inf 1.6536543 (4)
End of values pass after 6 iterations
6  Obj 2649026.4
Optimal - objective value 2649026.4
Optimal objective 2649026.448 - 3258 iterations time 0.392, Presolve 0.02
Total time (CPU seconds):       0.55   (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')
<__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.16it/s]
Writing variables.: 100%|██████████| 7/7 [00:00<00:00, 220.72it/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-bt4es_28.lp -solve -solu /tmp/linopy-solve-kq81zmwk.sol (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 2403 (-21425) rows, 4376 (-5564) columns and 15832 (-28394) elements
Perturbing problem by 0.001% of 6880.5027 - largest nonzero change 0.00094531053 ( 0.0038408689%) - largest zero change 0.00093735584
0  Obj 92709.643 Primal inf 5294157.5 (2230) Dual inf 2421.7098 (4)
123  Obj 23358.678 Primal inf 3986691.9 (2187)
246  Obj 23358.678 Primal inf 3636007.5 (2158)
369  Obj 23358.678 Primal inf 3340400.7 (2095)
492  Obj 23359.019 Primal inf 3212825.3 (2043)
615  Obj 23359.78 Primal inf 3857273.7 (2034)
738  Obj 23361.59 Primal inf 3384191.3 (1959)
861  Obj 23363.047 Primal inf 3278099.8 (1924)
984  Obj 23364.324 Primal inf 3099244.2 (1829)
1107  Obj 23365.58 Primal inf 2832810.5 (1759)
1230  Obj 23366.8 Primal inf 5077389.2 (1718)
1353  Obj 23368.968 Primal inf 34944403 (1767)
1476  Obj 40819.856 Primal inf 39192825 (1679)
1599  Obj 40821.832 Primal inf 32531853 (1532)
1722  Obj 40824.498 Primal inf 6119840.8 (1482)
1845  Obj 40826.38 Primal inf 3320538.3 (1242)
1968  Obj 40828.583 Primal inf 6601247 (1287)
2091  Obj 41012.957 Primal inf 5889155.6 (1023)
2214  Obj 214570.46 Primal inf 5491671.7 (895)
2337  Obj 329648.53 Primal inf 39291819 (851)
2460  Obj 526049.21 Primal inf 1488711.7 (691)
2583  Obj 862846.79 Primal inf 248516.75 (399)
2706  Obj 1548968.3 Primal inf 167092.36 (302)
2829  Obj 1912777 Primal inf 11161.952 (206)
2952  Obj 2086971.7 Primal inf 2369.306 (116)
3075  Obj 2128633.5 Primal inf 656.14339 (75)
3176  Obj 2142962.1
Optimal - objective value 2142956.9
After Postsolve, objective 2142956.9, infeasibilities - dual 175.7239 (5), primal 1.1686308e-06 (5)
Presolved model was optimal, full model needs cleaning up
0  Obj 2142956.9 Dual inf 1.7572385 (5)
End of values pass after 8 iterations
8  Obj 2142956.9
Optimal - objective value 2142956.9
Optimal objective 2142956.889 - 3184 iterations time 0.372, Presolve 0.02
Total time (CPU seconds):       0.56   (Wallclock seconds):       0.53


[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.23.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.23.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.047050 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.046404 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.046934 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.046306 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.046268 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.046489 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.046864 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.047729 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.046389 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.046826 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.046838 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.046469 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.046797 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.046616 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.046969 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.048192 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.046718 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.047584 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.047398 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.047128 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.047003 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.047146 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.047751 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.048307 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.23.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