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

Plot the distribution of the load and of generating tech

[3]:

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

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

/home/docs/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.25.1/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.25.1/lib/python3.11/site-packages/cartopy/mpl/style.py:76: UserWarning: facecolor will have no effect as it has been defined as "never".
  warnings.warn('facecolor will have no effect as it has been '
../_images/examples_scigrid-lopf-then-pf_7_1.png

Run Linear Optimal Power Flow on the first day of 2011.

To approximate n-1 security and allow room for reactive power flows, don’t allow any line to be loaded above 70% of their thermal rating

[7]:

contingency_factor = 0.7
network.lines.s_max_pu = contingency_factor

There are some infeasibilities without small extensions

[8]:

network.lines.loc[["316", "527", "602"], "s_nom"] = 1715

We performing a linear OPF for one day, 4 snapshots at a time.

[9]:

group_size = 4
network.storage_units.state_of_charge_initial = 0.0

for i in range(int(24 / group_size)):
    # set the initial state of charge based on previous round
    if i:
        network.storage_units.state_of_charge_initial = (
            network.storage_units_t.state_of_charge.loc[
                network.snapshots[group_size * i - 1]
            ]
        )
    network.optimize(
        network.snapshots[group_size * i : group_size * i + group_size],
        solver_name="cbc",
    )

WARNING:pypsa.components:The following transformers have zero r, which could break the linear load flow:
Index(['2', '5', '10', '12', '13', '15', '18', '20', '22', '24', '26', '30',
       '32', '37', '42', '46', '52', '56', '61', '68', '69', '74', '78', '86',
       '87', '94', '95', '96', '99', '100', '104', '105', '106', '107', '117',
       '120', '123', '124', '125', '128', '129', '138', '143', '156', '157',
       '159', '160', '165', '184', '191', '195', '201', '220', '231', '232',
       '233', '236', '247', '248', '250', '251', '252', '261', '263', '264',
       '267', '272', '279', '281', '282', '292', '303', '307', '308', '312',
       '315', '317', '322', '332', '334', '336', '338', '351', '353', '360',
       '362', '382', '384', '385', '391', '403', '404', '413', '421', '450',
       '458'],
      dtype='object', name='Transformer')
WARNING:pypsa.components:The following transformers have zero r, which could break the linear load flow:
Index(['2', '5', '10', '12', '13', '15', '18', '20', '22', '24', '26', '30',
       '32', '37', '42', '46', '52', '56', '61', '68', '69', '74', '78', '86',
       '87', '94', '95', '96', '99', '100', '104', '105', '106', '107', '117',
       '120', '123', '124', '125', '128', '129', '138', '143', '156', '157',
       '159', '160', '165', '184', '191', '195', '201', '220', '231', '232',
       '233', '236', '247', '248', '250', '251', '252', '261', '263', '264',
       '267', '272', '279', '281', '282', '292', '303', '307', '308', '312',
       '315', '317', '322', '332', '334', '336', '338', '351', '353', '360',
       '362', '382', '384', '385', '391', '403', '404', '413', '421', '450',
       '458'],
      dtype='object', name='Transformer')
INFO:linopy.model: Solve linear problem using Cbc solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 15/15 [00:00<00:00, 78.48it/s]
Writing continuous variables.: 100%|██████████| 7/7 [00:00<00:00, 250.50it/s]
INFO:linopy.io: Writing time: 0.23s
INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 9940 primals, 23828 duals
Objective: 1.45e+06
Solver model: not available
Solver message: Optimal - objective value 1449687.25014697


INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Line-fix-s-lower, Line-fix-s-upper, Transformer-fix-s-lower, Transformer-fix-s-upper, StorageUnit-fix-p_dispatch-lower, StorageUnit-fix-p_dispatch-upper, StorageUnit-fix-p_store-lower, StorageUnit-fix-p_store-upper, StorageUnit-fix-state_of_charge-lower, StorageUnit-fix-state_of_charge-upper, Kirchhoff-Voltage-Law, StorageUnit-energy_balance were not assigned to the network.

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

command line - cbc -printingOptions all -import /tmp/linopy-problem-cqv665cx.lp -solve -solu /tmp/linopy-solve-7brkwo9w.sol (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 2396 (-21432) rows, 4338 (-5602) columns and 17225 (-28825) elements
Perturbing problem by 0.001% of 2327.2654 - largest nonzero change 0.00099751112 ( 0.017779791%) - largest zero change 0.00099725216
0  Obj 88150.3 Primal inf 6561637.9 (2191) Dual inf 574.70586 (4)
122  Obj 18798.986 Primal inf 3533540.2 (2138)
244  Obj 18798.986 Primal inf 3347990.9 (2084)
366  Obj 18799.233 Primal inf 3288071.4 (2052)
488  Obj 18799.55 Primal inf 2965523.6 (1982)
610  Obj 18801.755 Primal inf 3041883.8 (1949)
732  Obj 18808.073 Primal inf 2442282.4 (1855)
854  Obj 18809.966 Primal inf 4461839.1 (1826)
976  Obj 18810.776 Primal inf 3966631.1 (1759)
1098  Obj 18811.843 Primal inf 4848582.3 (1752)
1220  Obj 18813.338 Primal inf 5956633 (1689)
1342  Obj 18815.102 Primal inf 7122586.2 (1665)
1464  Obj 22632.85 Primal inf 12814617 (1530)
1586  Obj 22634.612 Primal inf 44752656 (1528)
1708  Obj 22637.095 Primal inf 12110012 (1487)
1830  Obj 22638.993 Primal inf 5316640.1 (1231)
1952  Obj 22642.546 Primal inf 6856722.3 (1185)
2074  Obj 56217.442 Primal inf 3820580.2 (1184)
2196  Obj 73809.885 Primal inf 2163606.4 (885)
2318  Obj 174061.57 Primal inf 2044603.5 (795)
2440  Obj 447675.03 Primal inf 1906557.2 (656)
2562  Obj 479474.4 Primal inf 2549636.3 (542)
2684  Obj 719428.6 Primal inf 102884.02 (390)
2806  Obj 976321.97 Primal inf 233452.77 (482)
2928  Obj 1332105.1 Primal inf 67307.555 (200)
3050  Obj 1427979.3 Primal inf 1543.1247 (78)
3172  Obj 1449452.1 Primal inf 54.277877 (15)
3191  Obj 1449691.9
3191  Obj 1449687.3 Dual inf 7.3445901e-06 (1)
3192  Obj 1449687.3
Optimal - objective value 1449687.3
After Postsolve, objective 1449687.3, infeasibilities - dual 215.9672 (7), primal 1.5133967e-06 (7)
Presolved model was optimal, full model needs cleaning up
0  Obj 1449687.3 Primal inf 1.2635423e-07 (1) Dual inf 1e+08 (8)
End of values pass after 12 iterations
12  Obj 1449687.3
Optimal - objective value 1449687.3
Optimal objective 1449687.25 - 3204 iterations time 0.392, Presolve 0.02
Total time (CPU seconds):       0.59   (Wallclock seconds):       0.57


WARNING:pypsa.components:The following transformers have zero r, which could break the linear load flow:
Index(['2', '5', '10', '12', '13', '15', '18', '20', '22', '24', '26', '30',
       '32', '37', '42', '46', '52', '56', '61', '68', '69', '74', '78', '86',
       '87', '94', '95', '96', '99', '100', '104', '105', '106', '107', '117',
       '120', '123', '124', '125', '128', '129', '138', '143', '156', '157',
       '159', '160', '165', '184', '191', '195', '201', '220', '231', '232',
       '233', '236', '247', '248', '250', '251', '252', '261', '263', '264',
       '267', '272', '279', '281', '282', '292', '303', '307', '308', '312',
       '315', '317', '322', '332', '334', '336', '338', '351', '353', '360',
       '362', '382', '384', '385', '391', '403', '404', '413', '421', '450',
       '458'],
      dtype='object', name='Transformer')
WARNING:pypsa.components:The following transformers have zero r, which could break the linear load flow:
Index(['2', '5', '10', '12', '13', '15', '18', '20', '22', '24', '26', '30',
       '32', '37', '42', '46', '52', '56', '61', '68', '69', '74', '78', '86',
       '87', '94', '95', '96', '99', '100', '104', '105', '106', '107', '117',
       '120', '123', '124', '125', '128', '129', '138', '143', '156', '157',
       '159', '160', '165', '184', '191', '195', '201', '220', '231', '232',
       '233', '236', '247', '248', '250', '251', '252', '261', '263', '264',
       '267', '272', '279', '281', '282', '292', '303', '307', '308', '312',
       '315', '317', '322', '332', '334', '336', '338', '351', '353', '360',
       '362', '382', '384', '385', '391', '403', '404', '413', '421', '450',
       '458'],
      dtype='object', name='Transformer')
INFO:linopy.model: Solve linear problem using Cbc solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 15/15 [00:00<00:00, 80.65it/s]
Writing continuous variables.: 100%|██████████| 7/7 [00:00<00:00, 255.80it/s]
INFO:linopy.io: Writing time: 0.23s
INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 9940 primals, 23828 duals
Objective: 8.74e+05
Solver model: not available
Solver message: Optimal - objective value 873763.90232907


INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Line-fix-s-lower, Line-fix-s-upper, Transformer-fix-s-lower, Transformer-fix-s-upper, StorageUnit-fix-p_dispatch-lower, StorageUnit-fix-p_dispatch-upper, StorageUnit-fix-p_store-lower, StorageUnit-fix-p_store-upper, StorageUnit-fix-state_of_charge-lower, StorageUnit-fix-state_of_charge-upper, Kirchhoff-Voltage-Law, StorageUnit-energy_balance were not assigned to the network.

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

command line - cbc -printingOptions all -import /tmp/linopy-problem-yqy0nijg.lp -solve -solu /tmp/linopy-solve-g4bh3ufi.sol (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 2419 (-21409) rows, 4368 (-5572) columns and 17502 (-28548) elements
Perturbing problem by 0.001% of 3070.478 - largest nonzero change 0.00099526042 ( 0.019389923%) - largest zero change 0.00099451166
0  Obj 81193.865 Primal inf 5494255.5 (2219) Dual inf 1043.2235 (4)
123  Obj 11842.551 Primal inf 3830534.6 (2167)
246  Obj 11842.88 Primal inf 3043037.5 (2111)
369  Obj 11842.88 Primal inf 3000798.7 (2086)
492  Obj 11845.363 Primal inf 3243530.5 (2066)
615  Obj 11845.594 Primal inf 3116325.2 (1987)
738  Obj 11846.197 Primal inf 3106767.1 (1925)
861  Obj 11848.884 Primal inf 3550933.5 (1877)
984  Obj 11849.433 Primal inf 3956501 (1840)
1107  Obj 11850.964 Primal inf 2608199.3 (1722)
1230  Obj 11851.894 Primal inf 4413841.3 (1691)
1353  Obj 11853.518 Primal inf 6986413.9 (1658)
1476  Obj 11855.124 Primal inf 4771432 (1543)
1599  Obj 11857.272 Primal inf 17230040 (1528)
1722  Obj 11859.199 Primal inf 14911702 (1439)
1845  Obj 11861.172 Primal inf 7868857.1 (1346)
1968  Obj 11864.542 Primal inf 7633831.8 (1314)
2091  Obj 15366.27 Primal inf 7203194.1 (1206)
2214  Obj 21034.295 Primal inf 4788815.4 (967)
2337  Obj 29296.872 Primal inf 1090346.8 (717)
2460  Obj 115228.94 Primal inf 7471645.7 (818)
2583  Obj 294554.12 Primal inf 385292.79 (477)
2706  Obj 544453.14 Primal inf 170573.67 (432)
2829  Obj 722449.78 Primal inf 78653.084 (393)
2952  Obj 828889.96 Primal inf 10859.173 (140)
3075  Obj 857225.17 Primal inf 9904.9891 (111)
3198  Obj 873222.99 Primal inf 973.51594 (29)
3263  Obj 873767.79
Optimal - objective value 873763.9
After Postsolve, objective 873763.9, infeasibilities - dual 78.034695 (3), primal 9.155089e-07 (3)
Presolved model was optimal, full model needs cleaning up
0  Obj 873763.9 Dual inf 0.78034665 (3)
End of values pass after 3 iterations
3  Obj 873763.9
Optimal - objective value 873763.9
Optimal objective 873763.9023 - 3266 iterations time 0.422, Presolve 0.02
Total time (CPU seconds):       0.62   (Wallclock seconds):       0.59


WARNING:pypsa.components:The following transformers have zero r, which could break the linear load flow:
Index(['2', '5', '10', '12', '13', '15', '18', '20', '22', '24', '26', '30',
       '32', '37', '42', '46', '52', '56', '61', '68', '69', '74', '78', '86',
       '87', '94', '95', '96', '99', '100', '104', '105', '106', '107', '117',
       '120', '123', '124', '125', '128', '129', '138', '143', '156', '157',
       '159', '160', '165', '184', '191', '195', '201', '220', '231', '232',
       '233', '236', '247', '248', '250', '251', '252', '261', '263', '264',
       '267', '272', '279', '281', '282', '292', '303', '307', '308', '312',
       '315', '317', '322', '332', '334', '336', '338', '351', '353', '360',
       '362', '382', '384', '385', '391', '403', '404', '413', '421', '450',
       '458'],
      dtype='object', name='Transformer')
WARNING:pypsa.components:The following transformers have zero r, which could break the linear load flow:
Index(['2', '5', '10', '12', '13', '15', '18', '20', '22', '24', '26', '30',
       '32', '37', '42', '46', '52', '56', '61', '68', '69', '74', '78', '86',
       '87', '94', '95', '96', '99', '100', '104', '105', '106', '107', '117',
       '120', '123', '124', '125', '128', '129', '138', '143', '156', '157',
       '159', '160', '165', '184', '191', '195', '201', '220', '231', '232',
       '233', '236', '247', '248', '250', '251', '252', '261', '263', '264',
       '267', '272', '279', '281', '282', '292', '303', '307', '308', '312',
       '315', '317', '322', '332', '334', '336', '338', '351', '353', '360',
       '362', '382', '384', '385', '391', '403', '404', '413', '421', '450',
       '458'],
      dtype='object', name='Transformer')
INFO:linopy.model: Solve linear problem using Cbc solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 15/15 [00:00<00:00, 80.14it/s]
Writing continuous variables.: 100%|██████████| 7/7 [00:00<00:00, 251.53it/s]
INFO:linopy.io: Writing time: 0.23s
INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 9940 primals, 23828 duals
Objective: 7.91e+05
Solver model: not available
Solver message: Optimal - objective value 790786.59488191


INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Line-fix-s-lower, Line-fix-s-upper, Transformer-fix-s-lower, Transformer-fix-s-upper, StorageUnit-fix-p_dispatch-lower, StorageUnit-fix-p_dispatch-upper, StorageUnit-fix-p_store-lower, StorageUnit-fix-p_store-upper, StorageUnit-fix-state_of_charge-lower, StorageUnit-fix-state_of_charge-upper, Kirchhoff-Voltage-Law, StorageUnit-energy_balance were not assigned to the network.

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

command line - cbc -printingOptions all -import /tmp/linopy-problem-f7_bm3zp.lp -solve -solu /tmp/linopy-solve-etrw6amp.sol (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 2482 (-21346) rows, 4437 (-5503) columns and 17760 (-28290) elements
Perturbing problem by 0.001% of 2327.2654 - largest nonzero change 0.00099940449 ( 0.0089548564%) - largest zero change 0.00099290601
0  Obj -35.947231 Primal inf 5817841 (2285)
124  Obj -35.947231 Primal inf 3617850.8 (2236)
248  Obj -35.947231 Primal inf 3215341.9 (2164)
372  Obj -35.947231 Primal inf 3006826.5 (2117)
496  Obj -29.01886 Primal inf 2708219.4 (2076)
620  Obj -28.673879 Primal inf 3891582.8 (2073)
744  Obj -26.704935 Primal inf 2745783.9 (1997)
868  Obj -24.891039 Primal inf 2846365.7 (1933)
992  Obj -23.690271 Primal inf 2804675.9 (1844)
1116  Obj -22.242021 Primal inf 5328384.7 (1843)
1240  Obj -20.756623 Primal inf 2388080.3 (1730)
1364  Obj -19.156596 Primal inf 10867715 (1734)
1488  Obj -17.052352 Primal inf 6881646.6 (1643)
1612  Obj -15.253495 Primal inf 5079510 (1467)
1736  Obj 7648.308 Primal inf 3919663.6 (1295)
1860  Obj 7651.097 Primal inf 1587548.8 (1214)
1984  Obj 7653.5283 Primal inf 1742262.8 (1093)
2108  Obj 7655.9015 Primal inf 17312403 (1290)
2232  Obj 7724.7418 Primal inf 1183320.8 (833)
2356  Obj 7728.3851 Primal inf 753564.38 (571)
2480  Obj 88058.258 Primal inf 500047.18 (438)
2604  Obj 298462.98 Primal inf 115223.05 (422)
2728  Obj 506800.66 Primal inf 143882.05 (475)
2852  Obj 619572.93 Primal inf 30741.362 (276)
2976  Obj 692808.13 Primal inf 26716.3 (206)
3100  Obj 769053.42 Primal inf 6375.2674 (125)
3224  Obj 790780.39 Primal inf 90.297361 (11)
3237  Obj 790794.5
Optimal - objective value 790786.59
After Postsolve, objective 790786.59, infeasibilities - dual 0 (0), primal 0 (0)
Optimal objective 790786.5949 - 3237 iterations time 0.362, Presolve 0.02
Total time (CPU seconds):       0.56   (Wallclock seconds):       0.53


WARNING:pypsa.components:The following transformers have zero r, which could break the linear load flow:
Index(['2', '5', '10', '12', '13', '15', '18', '20', '22', '24', '26', '30',
       '32', '37', '42', '46', '52', '56', '61', '68', '69', '74', '78', '86',
       '87', '94', '95', '96', '99', '100', '104', '105', '106', '107', '117',
       '120', '123', '124', '125', '128', '129', '138', '143', '156', '157',
       '159', '160', '165', '184', '191', '195', '201', '220', '231', '232',
       '233', '236', '247', '248', '250', '251', '252', '261', '263', '264',
       '267', '272', '279', '281', '282', '292', '303', '307', '308', '312',
       '315', '317', '322', '332', '334', '336', '338', '351', '353', '360',
       '362', '382', '384', '385', '391', '403', '404', '413', '421', '450',
       '458'],
      dtype='object', name='Transformer')
WARNING:pypsa.components:The following transformers have zero r, which could break the linear load flow:
Index(['2', '5', '10', '12', '13', '15', '18', '20', '22', '24', '26', '30',
       '32', '37', '42', '46', '52', '56', '61', '68', '69', '74', '78', '86',
       '87', '94', '95', '96', '99', '100', '104', '105', '106', '107', '117',
       '120', '123', '124', '125', '128', '129', '138', '143', '156', '157',
       '159', '160', '165', '184', '191', '195', '201', '220', '231', '232',
       '233', '236', '247', '248', '250', '251', '252', '261', '263', '264',
       '267', '272', '279', '281', '282', '292', '303', '307', '308', '312',
       '315', '317', '322', '332', '334', '336', '338', '351', '353', '360',
       '362', '382', '384', '385', '391', '403', '404', '413', '421', '450',
       '458'],
      dtype='object', name='Transformer')
INFO:linopy.model: Solve linear problem using Cbc solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 15/15 [00:00<00:00, 78.58it/s]
Writing continuous variables.: 100%|██████████| 7/7 [00:00<00:00, 256.72it/s]
INFO:linopy.io: Writing time: 0.23s
INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 9940 primals, 23828 duals
Objective: 1.46e+06
Solver model: not available
Solver message: Optimal - objective value 1455232.67827293


INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Line-fix-s-lower, Line-fix-s-upper, Transformer-fix-s-lower, Transformer-fix-s-upper, StorageUnit-fix-p_dispatch-lower, StorageUnit-fix-p_dispatch-upper, StorageUnit-fix-p_store-lower, StorageUnit-fix-p_store-upper, StorageUnit-fix-state_of_charge-lower, StorageUnit-fix-state_of_charge-upper, Kirchhoff-Voltage-Law, StorageUnit-energy_balance were not assigned to the network.

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

command line - cbc -printingOptions all -import /tmp/linopy-problem-r55yr3by.lp -solve -solu /tmp/linopy-solve-3siuuln9.sol (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 2466 (-21362) rows, 4427 (-5513) columns and 17708 (-28342) elements
Perturbing problem by 0.001% of 2327.2654 - largest nonzero change 0.00099659414 ( 0.0084437175%) - largest zero change 0.00099290601
0  Obj 25318 Primal inf 5953235.3 (2283) Dual inf 257.04866 (1)
124  Obj 7980.1713 Primal inf 3678218.1 (2234)
248  Obj 7980.1713 Primal inf 3161230.9 (2154)
372  Obj 7980.1713 Primal inf 3159607.2 (2122)
496  Obj 7982.4593 Primal inf 3085333.4 (2087)
620  Obj 7984.9022 Primal inf 3765702.3 (2070)
744  Obj 7985.9546 Primal inf 11500073 (2026)
868  Obj 7989.4851 Primal inf 5994065.4 (1891)
992  Obj 7990.7625 Primal inf 2861796.1 (1806)
1116  Obj 7992.2846 Primal inf 9681223 (1816)
1240  Obj 7993.803 Primal inf 4222481.3 (1674)
1364  Obj 7995.6091 Primal inf 6182873.1 (1667)
1488  Obj 7997.6587 Primal inf 9164733.4 (1674)
1612  Obj 23733.32 Primal inf 8787498.2 (1554)
1736  Obj 23735.4 Primal inf 4137439.9 (1384)
1860  Obj 23737.272 Primal inf 1671432 (1195)
1984  Obj 32180.004 Primal inf 2642054.8 (1192)
2108  Obj 34252.04 Primal inf 534045.84 (736)
2232  Obj 61378.031 Primal inf 630584.17 (807)
2356  Obj 175374.62 Primal inf 700687.43 (751)
2480  Obj 394537.26 Primal inf 354882.56 (590)
2604  Obj 698709.87 Primal inf 110466.2 (415)
2728  Obj 1128187.8 Primal inf 42190.259 (327)
2852  Obj 1294953 Primal inf 24474.535 (226)
2976  Obj 1410172.6 Primal inf 12908.517 (128)
3100  Obj 1448466.5 Primal inf 3506.4167 (78)
3193  Obj 1455240
Optimal - objective value 1455232.7
After Postsolve, objective 1455232.7, infeasibilities - dual 24.259813 (1), primal 6.0436269e-07 (1)
Presolved model was optimal, full model needs cleaning up
0  Obj 1455232.7 Dual inf 0.28898438 (2)
End of values pass after 0 iterations
0  Obj 1455232.7 Dual inf 0.28898438 (2) w.o. free dual inf (1)
End of values pass after 1 iterations
1  Obj 1455232.7 Dual inf 0.046386348 (1) w.o. free dual inf (0)
2  Obj 1455232.7
Optimal - objective value 1455232.7
Optimal objective 1455232.678 - 3195 iterations time 0.402, 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')
INFO:linopy.model: Solve linear problem using Cbc solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 15/15 [00:00<00:00, 80.76it/s]
Writing continuous variables.: 100%|██████████| 7/7 [00:00<00:00, 246.60it/s]
INFO:linopy.io: Writing time: 0.23s
INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 9940 primals, 23828 duals
Objective: 2.65e+06
Solver model: not available
Solver message: Optimal - objective value 2649026.44753696


INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Line-fix-s-lower, Line-fix-s-upper, Transformer-fix-s-lower, Transformer-fix-s-upper, StorageUnit-fix-p_dispatch-lower, StorageUnit-fix-p_dispatch-upper, StorageUnit-fix-p_store-lower, StorageUnit-fix-p_store-upper, StorageUnit-fix-state_of_charge-lower, StorageUnit-fix-state_of_charge-upper, Kirchhoff-Voltage-Law, StorageUnit-energy_balance were not assigned to the network.

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

command line - cbc -printingOptions all -import /tmp/linopy-problem-k49vvuzj.lp -solve -solu /tmp/linopy-solve-zwb70h7x.sol (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 2413 (-21415) rows, 4373 (-5567) columns and 17577 (-28473) elements
Perturbing problem by 0.001% of 2327.2654 - largest nonzero change 0.00099909145 ( 0.0082541046%) - largest zero change 0.00099876841
0  Obj 97183.555 Primal inf 6522465.4 (2248) Dual inf 1028.1946 (4)
123  Obj 28003.782 Primal inf 3713931.9 (2187)
246  Obj 28003.782 Primal inf 3443959.2 (2127)
369  Obj 28003.782 Primal inf 4098306.5 (2087)
492  Obj 28004.679 Primal inf 3403566.6 (2040)
615  Obj 28005.071 Primal inf 4309864.7 (2008)
738  Obj 28010.794 Primal inf 3240961.7 (1942)
861  Obj 28013.925 Primal inf 11770243 (1871)
984  Obj 28015.306 Primal inf 6287239.9 (1819)
1107  Obj 28016.582 Primal inf 3207929.1 (1725)
1230  Obj 28018.722 Primal inf 2963397.2 (1618)
1353  Obj 28021.1 Primal inf 19211549 (1568)
1476  Obj 59483.261 Primal inf 4440951.4 (1402)
1599  Obj 59485.129 Primal inf 2.7045896e+08 (1412)
1722  Obj 59489.804 Primal inf 90875714 (1361)
1845  Obj 59494.725 Primal inf 38184749 (1135)
1968  Obj 59504.152 Primal inf 916463.94 (915)
2091  Obj 130785.08 Primal inf 1206319.2 (903)
2214  Obj 344478.12 Primal inf 1.4278501e+08 (798)
2337  Obj 665395.07 Primal inf 2.0601155e+09 (792)
2460  Obj 1069728.2 Primal inf 1.4649021e+09 (798)
2583  Obj 1701403.9 Primal inf 59045004 (427)
2706  Obj 2141705.6 Primal inf 1.4640981e+08 (325)
2829  Obj 2401937.7 Primal inf 19895.418 (248)
2952  Obj 2591218.9 Primal inf 7624.9249 (149)
3075  Obj 2633322.5 Primal inf 16130.894 (156)
3198  Obj 2648504.4 Primal inf 38.482883 (41)
3233  Obj 2649041.8
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 14 iterations
14  Obj 2649026.4
Optimal - objective value 2649026.4
Optimal objective 2649026.448 - 3247 iterations time 0.382, Presolve 0.02
Total time (CPU seconds):       0.58   (Wallclock seconds):       0.56


WARNING:pypsa.components:The following transformers have zero r, which could break the linear load flow:
Index(['2', '5', '10', '12', '13', '15', '18', '20', '22', '24', '26', '30',
       '32', '37', '42', '46', '52', '56', '61', '68', '69', '74', '78', '86',
       '87', '94', '95', '96', '99', '100', '104', '105', '106', '107', '117',
       '120', '123', '124', '125', '128', '129', '138', '143', '156', '157',
       '159', '160', '165', '184', '191', '195', '201', '220', '231', '232',
       '233', '236', '247', '248', '250', '251', '252', '261', '263', '264',
       '267', '272', '279', '281', '282', '292', '303', '307', '308', '312',
       '315', '317', '322', '332', '334', '336', '338', '351', '353', '360',
       '362', '382', '384', '385', '391', '403', '404', '413', '421', '450',
       '458'],
      dtype='object', name='Transformer')
WARNING:pypsa.components:The following transformers have zero r, which could break the linear load flow:
Index(['2', '5', '10', '12', '13', '15', '18', '20', '22', '24', '26', '30',
       '32', '37', '42', '46', '52', '56', '61', '68', '69', '74', '78', '86',
       '87', '94', '95', '96', '99', '100', '104', '105', '106', '107', '117',
       '120', '123', '124', '125', '128', '129', '138', '143', '156', '157',
       '159', '160', '165', '184', '191', '195', '201', '220', '231', '232',
       '233', '236', '247', '248', '250', '251', '252', '261', '263', '264',
       '267', '272', '279', '281', '282', '292', '303', '307', '308', '312',
       '315', '317', '322', '332', '334', '336', '338', '351', '353', '360',
       '362', '382', '384', '385', '391', '403', '404', '413', '421', '450',
       '458'],
      dtype='object', name='Transformer')
INFO:linopy.model: Solve linear problem using Cbc solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 15/15 [00:00<00:00, 79.27it/s]
Writing continuous variables.: 100%|██████████| 7/7 [00:00<00:00, 253.75it/s]
INFO:linopy.io: Writing time: 0.23s
INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 9940 primals, 23828 duals
Objective: 2.14e+06
Solver model: not available
Solver message: Optimal - objective value 2142956.88872811


INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Line-fix-s-lower, Line-fix-s-upper, Transformer-fix-s-lower, Transformer-fix-s-upper, StorageUnit-fix-p_dispatch-lower, StorageUnit-fix-p_dispatch-upper, StorageUnit-fix-p_store-lower, StorageUnit-fix-p_store-upper, StorageUnit-fix-state_of_charge-lower, StorageUnit-fix-state_of_charge-upper, Kirchhoff-Voltage-Law, StorageUnit-energy_balance were not assigned to the network.

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

command line - cbc -printingOptions all -import /tmp/linopy-problem-k6kfcn14.lp -solve -solu /tmp/linopy-solve-obvr6v58.sol (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 2419 (-21409) rows, 4370 (-5570) columns and 17602 (-28448) elements
Perturbing problem by 0.001% of 6880.5027 - largest nonzero change 0.00099919989 ( 0.0042006337%) - largest zero change 0.00099447956
0  Obj 92712.384 Primal inf 6319244.8 (2246) Dual inf 2421.7098 (4)
123  Obj 23361.441 Primal inf 4040192.4 (2206)
246  Obj 23361.441 Primal inf 3262836.9 (2154)
369  Obj 23361.441 Primal inf 2950810.8 (2079)
492  Obj 23362.526 Primal inf 3054064.2 (2052)
615  Obj 23363.285 Primal inf 2726340.5 (2002)
738  Obj 23364.506 Primal inf 2450128.5 (1918)
861  Obj 23365.55 Primal inf 3743974 (1873)
984  Obj 23366.807 Primal inf 3471253.9 (1820)
1107  Obj 23367.491 Primal inf 3430092.7 (1776)
1230  Obj 23368.564 Primal inf 3619425.4 (1668)
1353  Obj 23370.054 Primal inf 4148296 (1634)
1476  Obj 40820.437 Primal inf 3656378.9 (1460)
1599  Obj 40821.89 Primal inf 6246520.5 (1392)
1722  Obj 40823.744 Primal inf 13163866 (1275)
1845  Obj 40824.91 Primal inf 4691517.3 (1215)
1968  Obj 40827.697 Primal inf 25939416 (1116)
2091  Obj 60976.185 Primal inf 10484092 (990)
2214  Obj 185987.22 Primal inf 22073928 (890)
2337  Obj 343497.23 Primal inf 9168656.6 (789)
2460  Obj 806246.97 Primal inf 3719368.6 (599)
2583  Obj 1337036.3 Primal inf 141960.96 (425)
2706  Obj 1726216.5 Primal inf 28189.683 (295)
2829  Obj 1877231.4 Primal inf 14460.612 (250)
2952  Obj 2029306 Primal inf 26484.117 (311)
3075  Obj 2110041.7 Primal inf 2402.7327 (117)
3198  Obj 2135814.4 Primal inf 534.25763 (72)
3296  Obj 2142961.2
Optimal - objective value 2142956.9
After Postsolve, objective 2142956.9, infeasibilities - dual 137.65508 (4), primal 1.0570064e-06 (4)
Presolved model was optimal, full model needs cleaning up
0  Obj 2142956.9 Dual inf 1.3765504 (4)
End of values pass after 5 iterations
5  Obj 2142956.9
Optimal - objective value 2142956.9
Optimal objective 2142956.889 - 3301 iterations time 0.402, Presolve 0.02
Total time (CPU seconds):       0.59   (Wallclock seconds):       0.56


[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.258607
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.25.1/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.25.1/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.043731 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.044072 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.043288 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.044316 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.042980 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.044468 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.043098 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.043186 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.043176 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.042891 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.044288 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.044505 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.042956 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.042895 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.042972 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.043036 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.043306 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.043199 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.043259 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.042990 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.043490 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.043404 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.042988 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.043140 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.000435
std        0.260958
min       -0.813383
25%       -0.125478
50%        0.003032
75%        0.126690
max        0.827203
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.022968
std        2.373833
min      -12.158423
25%       -0.462981
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.25.1/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