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.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.25.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.25.0/lib/python3.11/site-packages/cartopy/mpl/style.py:76: UserWarning: facecolor will have no effect as it has been defined as "never".
  warnings.warn('facecolor will have no effect as it has been '
../_images/examples_scigrid-lopf-then-pf_7_1.png

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

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

[7]:

contingency_factor = 0.7
network.lines.s_max_pu = contingency_factor

There are some infeasibilities without small extensions

[8]:

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

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

[9]:

group_size = 4
network.storage_units.state_of_charge_initial = 0.0

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

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


INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Kirchhoff-Voltage-Law were not assigned to the network.

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

command line - cbc -printingOptions all -import /tmp/linopy-problem-9qjz8stv.lp -solve -solu /tmp/linopy-solve-8bfa3ml9.sol (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 2410 (-21418) rows, 4353 (-5587) columns and 17359 (-28575) elements
Perturbing problem by 0.001% of 2454.8041 - largest nonzero change 0.00099832053 ( 0.0092323128%) - largest zero change 0.00099672906
0  Obj 88140.86 Primal inf 6084716.3 (2205) Dual inf 329.5593 (4)
123  Obj 18789.546 Primal inf 4306398.5 (2161)
246  Obj 18789.546 Primal inf 3583236.4 (2121)
369  Obj 18789.546 Primal inf 3184697.5 (2066)
492  Obj 18790.354 Primal inf 3131813.4 (2013)
615  Obj 18791.526 Primal inf 3082292.7 (1967)
738  Obj 18792.758 Primal inf 2907960.9 (1920)
861  Obj 18795.184 Primal inf 3086091.6 (1874)
984  Obj 18798.6 Primal inf 2495927.6 (1740)
1107  Obj 18800.036 Primal inf 3733422.1 (1762)
1230  Obj 18802.139 Primal inf 4756395.8 (1716)
1353  Obj 18805.77 Primal inf 44396497 (1691)
1476  Obj 18808.922 Primal inf 17745727 (1575)
1599  Obj 22628.151 Primal inf 14163577 (1668)
1722  Obj 22631.967 Primal inf 24635561 (1568)
1845  Obj 22636.99 Primal inf 37338688 (1531)
1968  Obj 22641.946 Primal inf 32318699 (1327)
2091  Obj 44542.392 Primal inf 95108608 (1095)
2214  Obj 51234.636 Primal inf 2.5121574e+08 (938)
2337  Obj 99672.587 Primal inf 4.0164952e+08 (1113)
2460  Obj 315705.62 Primal inf 1153201.7 (730)
2583  Obj 535270.17 Primal inf 3690717.4 (578)
2706  Obj 798556.22 Primal inf 836799.29 (411)
2829  Obj 1185911.4 Primal inf 318644.43 (343)
2952  Obj 1381557 Primal inf 64700.221 (261)
3075  Obj 1440029.2 Primal inf 30578.853 (130)
3198  Obj 1449431.8 Primal inf 361.22998 (20)
3233  Obj 1449695.1
3233  Obj 1449687.3 Dual inf 7.3445901e-06 (1)
3234  Obj 1449687.3
Optimal - objective value 1449687.3
After Postsolve, objective 1449687.3, infeasibilities - dual 243.25253 (8), primal 1.9208568e-06 (8)
Presolved model was optimal, full model needs cleaning up
0  Obj 1449687.3 Primal inf 1.2635402e-07 (1) Dual inf 1.0000001e+08 (16)
End of values pass after 0 iterations
0  Obj 1449687.3 Primal inf 1.2635402e-07 (1) Dual inf 1.0000001e+08 (16) w.o. free dual inf (15)
End of values pass after 12 iterations
12  Obj 1449687.3 Dual inf 9.80985 (7) w.o. free dual inf (6)
13  Obj 1449687.3
Optimal - objective value 1449687.3
Optimal objective 1449687.25 - 3247 iterations time 0.412, 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')
INFO:linopy.model: Solve linear problem using Cbc solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 15/15 [00:00<00:00, 73.77it/s]
Writing continuous variables.: 100%|██████████| 7/7 [00:00<00:00, 231.70it/s]
INFO:linopy.io: Writing time: 0.25s
INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 9940 primals, 23828 duals
Objective: 8.74e+05
Solver model: not available
Solver message: Optimal - objective value 873763.90232907


INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Kirchhoff-Voltage-Law were not assigned to the network.

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

command line - cbc -printingOptions all -import /tmp/linopy-problem-w4ym3tan.lp -solve -solu /tmp/linopy-solve-srkclulj.sol (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 2431 (-21397) rows, 4386 (-5554) columns and 17475 (-28459) elements
Perturbing problem by 0.001% of 2620.1345 - largest nonzero change 0.00099700892 ( 0.031892336%) - largest zero change 0.00099607519
0  Obj 81173.04 Primal inf 5698083.6 (2232) Dual inf 447.30519 (4)
123  Obj 11821.726 Primal inf 4260799.9 (2198)
246  Obj 11823.35 Primal inf 3793102.7 (2123)
369  Obj 11823.353 Primal inf 3410096 (2088)
492  Obj 11823.351 Primal inf 3013960 (2044)
615  Obj 11823.654 Primal inf 2745108.5 (1961)
738  Obj 11824.746 Primal inf 2883575.2 (1954)
861  Obj 11826.736 Primal inf 3440808 (1909)
984  Obj 11831.309 Primal inf 5080654.5 (1869)
1107  Obj 11836.078 Primal inf 7969727 (1768)
1230  Obj 11837.736 Primal inf 7329787.3 (1787)
1353  Obj 11839.213 Primal inf 10436725 (1720)
1476  Obj 11842.938 Primal inf 55518068 (1675)
1599  Obj 11848.627 Primal inf 36436888 (1465)
1722  Obj 11853.046 Primal inf 9397408.5 (1368)
1845  Obj 11858.402 Primal inf 24102585 (1322)
1968  Obj 11861.553 Primal inf 5505447.6 (1261)
2091  Obj 11867.184 Primal inf 3773037.5 (1043)
2214  Obj 11873.132 Primal inf 2145380.3 (842)
2337  Obj 49922.262 Primal inf 22587915 (918)
2460  Obj 138341.25 Primal inf 68780172 (883)
2583  Obj 254219.54 Primal inf 278285.49 (424)
2706  Obj 486281.1 Primal inf 156438.95 (352)
2829  Obj 675879.54 Primal inf 28143.183 (225)
2952  Obj 815631.56 Primal inf 58970.849 (173)
3075  Obj 868070.92 Primal inf 30727.814 (146)
3198  Obj 873765.08 Primal inf 49.566208 (11)
3210  Obj 873769.92
3210  Obj 873763.94 Dual inf 9.5224768e-05 (1)
3213  Obj 873763.9
Optimal - objective value 873763.9
After Postsolve, objective 873763.9, infeasibilities - dual 78.034695 (3), primal 9.1550868e-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 - 3216 iterations time 0.392, Presolve 0.02
Total time (CPU seconds):       0.57   (Wallclock seconds):       0.55


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


INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Kirchhoff-Voltage-Law were not assigned to the network.

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

command line - cbc -printingOptions all -import /tmp/linopy-problem-4585664p.lp -solve -solu /tmp/linopy-solve-_tw1gff3.sol (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 2502 (-21326) rows, 4459 (-5481) columns and 17714 (-28220) elements
Perturbing problem by 0.001% of 2745.1675 - largest nonzero change 0.00099862309 ( 0.0089469271%) - largest zero change 0.00099692494
0  Obj -41.834805 Primal inf 5679983.1 (2305)
125  Obj -41.834805 Primal inf 4337926.1 (2257)
250  Obj -41.428796 Primal inf 3913959.1 (2195)
375  Obj -40.250334 Primal inf 3325637.2 (2133)
500  Obj -40.250334 Primal inf 3263029.3 (2097)
625  Obj -40.065742 Primal inf 3158994.9 (2060)
750  Obj -39.893192 Primal inf 3131828.6 (2005)
875  Obj -37.965938 Primal inf 2946158.4 (1934)
1000  Obj -35.188331 Primal inf 2251450.6 (1834)
1125  Obj -33.795251 Primal inf 3064673.8 (1834)
1250  Obj -30.638947 Primal inf 2026927.6 (1647)
1375  Obj -28.147552 Primal inf 2109590.2 (1596)
1500  Obj -24.617826 Primal inf 4653146.9 (1617)
1625  Obj -19.979951 Primal inf 8066609.3 (1552)
1750  Obj 7644.3653 Primal inf 3458681.6 (1388)
1875  Obj 7648.4023 Primal inf 4452502.1 (1376)
2000  Obj 7652.5276 Primal inf 20669445 (1357)
2125  Obj 7656.8458 Primal inf 2709617.2 (1065)
2250  Obj 7719.0474 Primal inf 3106411.9 (1017)
2375  Obj 40037.343 Primal inf 689966.92 (666)
2500  Obj 139484.31 Primal inf 454034.64 (477)
2625  Obj 330454.53 Primal inf 146013.35 (415)
2750  Obj 496259.6 Primal inf 54926.512 (330)
2875  Obj 636585.27 Primal inf 79484.24 (329)
3000  Obj 731192.39 Primal inf 10241.738 (170)
3125  Obj 772479.73 Primal inf 2118.7682 (96)
3250  Obj 790754.18 Primal inf 16.386815 (10)
3261  Obj 790794.91
Optimal - objective value 790786.59
After Postsolve, objective 790786.59, infeasibilities - dual 0 (0), primal 0 (0)
Optimal objective 790786.5949 - 3261 iterations time 0.382, Presolve 0.02
Total time (CPU seconds):       0.56   (Wallclock seconds):       0.54


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


INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Kirchhoff-Voltage-Law were not assigned to the network.

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

command line - cbc -printingOptions all -import /tmp/linopy-problem-3_lwpxsa.lp -solve -solu /tmp/linopy-solve-4t9acxd9.sol (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 2485 (-21343) rows, 4451 (-5489) columns and 17674 (-28260) elements
Perturbing problem by 0.001% of 2454.8041 - largest nonzero change 0.00099435818 ( 0.0075378327%) - largest zero change 0.00098881386
0  Obj 25303.664 Primal inf 7227382.1 (2302) Dual inf 132.81465 (1)
124  Obj 7965.835 Primal inf 4379234.7 (2270)
248  Obj 7965.835 Primal inf 3659509.9 (2207)
372  Obj 7967.0413 Primal inf 3308482 (2138)
496  Obj 7967.3651 Primal inf 2972298 (2087)
620  Obj 7968.0961 Primal inf 3335007.6 (2067)
744  Obj 7969.2247 Primal inf 3257972.2 (2014)
868  Obj 7971.0334 Primal inf 2963495.4 (1966)
992  Obj 7976.3347 Primal inf 5301274.4 (1888)
1116  Obj 7978.2399 Primal inf 16306251 (1870)
1240  Obj 7981.6505 Primal inf 7071746.7 (1837)
1364  Obj 7983.2165 Primal inf 6524511 (1769)
1488  Obj 7987.6663 Primal inf 15186037 (1674)
1612  Obj 12742.507 Primal inf 4399942.2 (1548)
1736  Obj 23727.809 Primal inf 3408276.5 (1421)
1860  Obj 23731.922 Primal inf 3709798.1 (1285)
1984  Obj 23738.957 Primal inf 2613298.5 (1174)
2108  Obj 23744.269 Primal inf 29528825 (1193)
2232  Obj 23752.416 Primal inf 34120116 (1105)
2356  Obj 142051.69 Primal inf 32369688 (929)
2480  Obj 313481.99 Primal inf 1266069.8 (777)
2604  Obj 407138.38 Primal inf 3144248.9 (579)
2728  Obj 591707.49 Primal inf 1240886.6 (470)
2852  Obj 956738.49 Primal inf 182762.16 (364)
2976  Obj 1176145.3 Primal inf 24861.808 (268)
3100  Obj 1369622.3 Primal inf 10091.07 (163)
3224  Obj 1442408.2 Primal inf 3093.7711 (83)
3348  Obj 1455239.1 Primal inf 0.10959895 (4)
3352  Obj 1455243.1
Optimal - objective value 1455232.7
After Postsolve, objective 1455232.7, infeasibilities - dual 49.951498 (2), primal 2.7246272e-07 (2)
Presolved model was optimal, full model needs cleaning up
0  Obj 1455232.7 Dual inf 0.70043366 (4)
End of values pass after 0 iterations
0  Obj 1455232.7 Dual inf 0.70043366 (4) w.o. free dual inf (2)
End of values pass after 2 iterations
2  Obj 1455232.7 Dual inf 0.20091888 (2) w.o. free dual inf (0)
4  Obj 1455232.7
Optimal - objective value 1455232.7
Optimal objective 1455232.678 - 3356 iterations time 0.422, Presolve 0.02
Total time (CPU seconds):       0.65   (Wallclock seconds):       0.60


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, 72.71it/s]
Writing continuous variables.: 100%|██████████| 7/7 [00:00<00:00, 227.94it/s]
INFO:linopy.io: Writing time: 0.25s
INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 9940 primals, 23828 duals
Objective: 2.65e+06
Solver model: not available
Solver message: Optimal - objective value 2649026.44753696


INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Kirchhoff-Voltage-Law were not assigned to the network.

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

command line - cbc -printingOptions all -import /tmp/linopy-problem-efiacw75.lp -solve -solu /tmp/linopy-solve-ufovuw9s.sol (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 2426 (-21402) rows, 4392 (-5548) columns and 17527 (-28407) elements
Perturbing problem by 0.001% of 2454.8041 - largest nonzero change 0.00099831964 ( 0.010521421%) - largest zero change 0.00099672906
0  Obj 97176.228 Primal inf 7542643.9 (2261) Dual inf 531.25856 (4)
123  Obj 27996.457 Primal inf 4844128.8 (2231)
246  Obj 27996.457 Primal inf 3994473.7 (2171)
369  Obj 27996.457 Primal inf 3237270.4 (2130)
492  Obj 27997.595 Primal inf 3016801.2 (2080)
615  Obj 27998.188 Primal inf 3903735.8 (2022)
738  Obj 27999.594 Primal inf 2850935.8 (1953)
861  Obj 28002.511 Primal inf 3003972.6 (1890)
984  Obj 28008.363 Primal inf 5444607.3 (1863)
1107  Obj 28011.395 Primal inf 6511075.8 (1843)
1230  Obj 28014.715 Primal inf 5022781.5 (1736)
1353  Obj 28018.975 Primal inf 3968364.4 (1558)
1476  Obj 59481.633 Primal inf 3183728.1 (1460)
1599  Obj 59483.902 Primal inf 7314146.3 (1485)
1722  Obj 59487.872 Primal inf 9048018.3 (1343)
1845  Obj 59495.417 Primal inf 1.2070481e+08 (1477)
1968  Obj 59502.242 Primal inf 3355483.3 (1310)
2091  Obj 60374.728 Primal inf 6992183.2 (1239)
2214  Obj 191146.06 Primal inf 7036687.8 (1121)
2337  Obj 366005.97 Primal inf 4085402.9 (938)
2460  Obj 742084.09 Primal inf 3244596.7 (757)
2583  Obj 1286456.9 Primal inf 6.4834068e+09 (689)
2706  Obj 1901304.7 Primal inf 3638375.1 (458)
2829  Obj 2231098.6 Primal inf 27168.641 (252)
2952  Obj 2477209.8 Primal inf 42619.637 (227)
3075  Obj 2599924.4 Primal inf 8646.9262 (165)
3198  Obj 2639675.3 Primal inf 3490.0862 (116)
3321  Obj 2648898.2 Primal inf 13.850264 (23)
3345  Obj 2649044.1
Optimal - objective value 2649026.4
After Postsolve, objective 2649026.4, infeasibilities - dual 124.09584 (3), primal 7.3815671e-07 (3)
Presolved model was optimal, full model needs cleaning up
0  Obj 2649026.4 Dual inf 1.2409581 (3)
End of values pass after 4 iterations
4  Obj 2649026.4
Optimal - objective value 2649026.4
Optimal objective 2649026.448 - 3349 iterations time 0.412, Presolve 0.02
Total time (CPU seconds):       0.60   (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, 75.59it/s]
Writing continuous variables.: 100%|██████████| 7/7 [00:00<00:00, 237.69it/s]
INFO:linopy.io: Writing time: 0.24s
INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 9940 primals, 23828 duals
Objective: 2.14e+06
Solver model: not available
Solver message: Optimal - objective value 2142956.88872811


INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Kirchhoff-Voltage-Law were not assigned to the network.

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

command line - cbc -printingOptions all -import /tmp/linopy-problem-y1pqfhy2.lp -solve -solu /tmp/linopy-solve-_de7zo5q.sol (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 2431 (-21397) rows, 4386 (-5554) columns and 17537 (-28397) elements
Perturbing problem by 0.001% of 8902.5915 - largest nonzero change 0.00097370854 ( 0.0040748628%) - largest zero change 0.00096865049
0  Obj 92710.211 Primal inf 7647212.1 (2258) Dual inf 1195.3238 (4)
123  Obj 23359.333 Primal inf 3973273.9 (2221)
246  Obj 23359.333 Primal inf 3468408.5 (2158)
369  Obj 23359.514 Primal inf 3345991 (2122)
492  Obj 23359.542 Primal inf 3815704.4 (2083)
615  Obj 23359.577 Primal inf 5540384.5 (2051)
738  Obj 23360.633 Primal inf 3197917.2 (1949)
861  Obj 23361.651 Primal inf 4195108.4 (1924)
984  Obj 23363.08 Primal inf 7528182.4 (1876)
1107  Obj 23365.819 Primal inf 2522607.9 (1762)
1230  Obj 23367.181 Primal inf 1884859.8 (1650)
1353  Obj 23369.035 Primal inf 2081224.7 (1639)
1476  Obj 23371.141 Primal inf 3218384.1 (1572)
1599  Obj 40822.702 Primal inf 9731638.5 (1540)
1722  Obj 40825.123 Primal inf 8562290.5 (1350)
1845  Obj 40828.263 Primal inf 4514591.2 (1271)
1968  Obj 40830.105 Primal inf 1203051 (1129)
2091  Obj 40832.225 Primal inf 2607965.5 (1022)
2214  Obj 59388.409 Primal inf 42960916 (1393)
2337  Obj 220024.87 Primal inf 3879726.3 (977)
2460  Obj 345323.49 Primal inf 705229.29 (858)
2583  Obj 822658.91 Primal inf 1377450.9 (624)
2706  Obj 1408231 Primal inf 98704.041 (347)
2829  Obj 1822469.2 Primal inf 52241.387 (434)
2952  Obj 1991560.2 Primal inf 10665.086 (243)
3075  Obj 2104981.3 Primal inf 1572.1526 (113)
3198  Obj 2139923.1 Primal inf 220.96157 (49)
3317  Obj 2142962.9
Optimal - objective value 2142956.9
After Postsolve, objective 2142956.9, infeasibilities - dual 214.76203 (6), primal 1.5457266e-06 (6)
Presolved model was optimal, full model needs cleaning up
0  Obj 2142956.9 Dual inf 2.1476197 (6)
End of values pass after 7 iterations
7  Obj 2142956.9
Optimal - objective value 2142956.9
Optimal objective 2142956.889 - 3324 iterations time 0.432, Presolve 0.02
Total time (CPU seconds):       0.63   (Wallclock seconds):       0.59


[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.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.25.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.046904 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.045363 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.049335 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.045262 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.045631 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.045453 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.045372 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.045028 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.045739 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.044974 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.045147 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.046248 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.044705 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.044878 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.044741 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.044976 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.045192 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.045137 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.045367 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.046620 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.046184 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.046924 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.044996 seconds
INFO:pypsa.pf:Newton-Raphson solved in 4 iterations with error of 0.000000 in 0.044976 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.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