Note

You can download this example as a Jupyter notebook or start it in interactive mode.

Redispatch Example with SciGRID network#

In this example, we compare a 2-stage market with an initial market clearing in two bidding zones with flow-based market coupling and a subsequent redispatch market (incl. curtailment) to an idealised nodal pricing scheme.

[1]:

import pypsa
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from pypsa.descriptors import get_switchable_as_dense as as_dense

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

[2]:

solver = "cbc"

Load example network#

[3]:

o = pypsa.examples.scigrid_de(from_master=True)
o.lines.s_max_pu = 0.7
o.lines.loc[["316", "527", "602"], "s_nom"] = 1715
o.set_snapshots([o.snapshots[12]])

WARNING:pypsa.io:Importing network from PyPSA version v0.17.1 while current version is v0.24.0. Read the release notes at https://pypsa.readthedocs.io/en/latest/release_notes.html to prepare your network for import.
INFO:pypsa.io:Imported network scigrid-de.nc has buses, generators, lines, loads, storage_units, transformers

[4]:

n = o.copy()  # for redispatch model
m = o.copy()  # for market model

[5]:

o.plot();

/home/docs/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.24.0/lib/python3.11/site-packages/cartopy/mpl/style.py:76: UserWarning: facecolor will have no effect as it has been defined as "never".
  warnings.warn('facecolor will have no effect as it has been '
../_images/examples_scigrid-redispatch_6_1.png

Solve original nodal market model o#

First, let us solve a nodal market using the original model o:

[6]:

o.optimize(solver_name=solver)

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 time: 0.12s
INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 2485 primals, 5957 duals
Objective: 3.01e+05
Solver model: not available
Solver message: Optimal - objective value 301209.38232509



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

command line - cbc -printingOptions all -import /tmp/linopy-problem-kk_hj2h8.lp -solve -solu /tmp/linopy-solve-qadvu5ri.sol (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 621 (-5336) rows, 1086 (-1399) columns and 3950 (-7069) elements
Perturbing problem by 0.001% of 2386.6826 - largest nonzero change 0.00098850976 ( 0.0079713351%) - largest zero change 0.0009873664
0  Obj -12.389216 Primal inf 1416482.1 (573)
87  Obj -11.955742 Primal inf 726036.62 (536)
174  Obj -10.558214 Primal inf 759152.92 (497)
261  Obj -9.4082981 Primal inf 663308.08 (469)
333  Obj -7.9420554 Primal inf 631600.28 (417)
412  Obj 3996.3541 Primal inf 555146.15 (322)
486  Obj 3999.023 Primal inf 1196052.8 (339)
553  Obj 4002.155 Primal inf 610168.03 (269)
635  Obj 4006.1589 Primal inf 118075.34 (160)
722  Obj 222031.83 Primal inf 129315.99 (79)
793  Obj 301212.25
Optimal - objective value 301209.38
After Postsolve, objective 301209.38, infeasibilities - dual 0 (0), primal 0 (0)
Optimal objective 301209.3823 - 793 iterations time 0.112, Presolve 0.02
Total time (CPU seconds):       0.15   (Wallclock seconds):       0.12


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

[6]:

('ok', 'optimal')

Costs are 301 k€.

Build market model m with two bidding zones#

For this example, we split the German transmission network into two market zones at latitude 51 degrees.

You can build any other market zones by providing an alternative mapping from bus to zone.

[7]:

zones = (n.buses.y > 51).map(lambda x: "North" if x else "South")

Next, we assign this mapping to the market model m.

We re-assign the buses of all generators and loads, and remove all transmission lines within each bidding zone.

Here, we assume that the bidding zones are coupled through the transmission lines that connect them.

[8]:

for c in m.iterate_components(m.one_port_components):
    c.df.bus = c.df.bus.map(zones)

for c in m.iterate_components(m.branch_components):
    c.df.bus0 = c.df.bus0.map(zones)
    c.df.bus1 = c.df.bus1.map(zones)
    internal = c.df.bus0 == c.df.bus1
    m.mremove(c.name, c.df.loc[internal].index)

m.mremove("Bus", m.buses.index)
m.madd("Bus", ["North", "South"]);

Now, we can solve the coupled market with two bidding zones.

[9]:

m.optimize(solver_name=solver)

INFO:linopy.model: Solve linear problem using Cbc solver
INFO:linopy.io: Writing time: 0.08s
INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 1561 primals, 3185 duals
Objective: 2.14e+05
Solver model: not available
Solver message: Optimal - objective value 213988.68595810


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-sx_wj2zy.lp -solve -solu /tmp/linopy-solve-q7he_gbd.sol (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 40 (-3145) rows, 410 (-1151) columns and 487 (-4342) elements
Perturbing problem by 0.001% of 212.59539 - largest nonzero change 0.00017578427 ( 0.0036987348%) - largest zero change 0.00015445146
0  Obj 0 Primal inf 11285.222 (1)
48  Obj 184184.9 Primal inf 1700.1029 (24)
86  Obj 213988.73
Optimal - objective value 213988.69
After Postsolve, objective 213988.69, infeasibilities - dual 0 (0), primal 0 (0)
Optimal objective 213988.686 - 86 iterations time 0.012, Presolve 0.01
Total time (CPU seconds):       0.04   (Wallclock seconds):       0.03


[9]:

('ok', 'optimal')

Costs are 214 k€, which is much lower than the 301 k€ of the nodal market.

This is because network restrictions apart from the North/South division are not taken into account yet.

We can look at the market clearing prices of each zone:

[10]:

m.buses_t.marginal_price

[10]:
Bus North South
snapshot
2011-01-01 12:00:00 8.0 25.0

Build redispatch model n#

Next, based on the market outcome with two bidding zones m, we build a secondary redispatch market n that rectifies transmission constraints through curtailment and ramping up/down thermal generators.

First, we fix the dispatch of generators to the results from the market simulation. (For simplicity, this example disregards storage units.)

[11]:

p = m.generators_t.p / m.generators.p_nom
n.generators_t.p_min_pu = p
n.generators_t.p_max_pu = p

Then, we add generators bidding into redispatch market using the following assumptions:

  • All generators can reduce their dispatch to zero. This includes also curtailment of renewables.

  • All generators can increase their dispatch to their available/nominal capacity.

  • No changes to the marginal costs, i.e. reducing dispatch lowers costs.

With these settings, the 2-stage market should result in the same cost as the nodal market.

[12]:

g_up = n.generators.copy()
g_down = n.generators.copy()

g_up.index = g_up.index.map(lambda x: x + " ramp up")
g_down.index = g_down.index.map(lambda x: x + " ramp down")

up = (
    as_dense(m, "Generator", "p_max_pu") * m.generators.p_nom - m.generators_t.p
).clip(0) / m.generators.p_nom
down = -m.generators_t.p / m.generators.p_nom

up.columns = up.columns.map(lambda x: x + " ramp up")
down.columns = down.columns.map(lambda x: x + " ramp down")

n.madd("Generator", g_up.index, p_max_pu=up, **g_up.drop("p_max_pu", axis=1))

n.madd(
    "Generator",
    g_down.index,
    p_min_pu=down,
    p_max_pu=0,
    **g_down.drop(["p_max_pu", "p_min_pu"], axis=1)
);

Now, let’s solve the redispatch market:

[13]:

n.optimize(solver_name=solver)

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 time: 0.15s
INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 5331 primals, 11649 duals
Objective: 3.01e+05
Solver model: not available
Solver message: Optimal - objective value 301209.38114435


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-1b9n7eip.lp -solve -solu /tmp/linopy-solve-e9pwm5f9.sol (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 624 (-11025) rows, 1326 (-4005) columns and 4213 (-15344) elements
Perturbing problem by 0.001% of 2420.52 - largest nonzero change 0.00098398274 ( 0.0085447571%) - largest zero change 0.00098270306
0  Obj 195154.11 Primal inf 1437635.3 (575) Dual inf 8051.4193 (158)
87  Obj -12.188894 Primal inf 733262.65 (537)
174  Obj -11.288618 Primal inf 1404232.7 (517)
261  Obj -9.6729363 Primal inf 1233892.3 (476)
348  Obj -7.4900378 Primal inf 520600.22 (395)
428  Obj 3996.6412 Primal inf 850063.86 (335)
501  Obj 3999.0691 Primal inf 340317.79 (290)
557  Obj 4034.6208 Primal inf 549962.33 (251)
634  Obj 4041.8437 Primal inf 822589.98 (278)
721  Obj 179168.61 Primal inf 17566.972 (89)
808  Obj 301212.15 Primal inf 0.14142613 (2)
809  Obj 301212.15
809  Obj 301209.38 Dual inf 6.4709158e-05 (4)
813  Obj 301209.38
Optimal - objective value 301209.38
After Postsolve, objective 301209.38, infeasibilities - dual 1532.6511 (103), primal 2.2654494e-05 (96)
Presolved model was optimal, full model needs cleaning up
0  Obj 301209.38 Primal inf 1.0448869e-06 (6) Dual inf 6.0000001e+08 (109)
End of values pass after 110 iterations
110  Obj 301209.38
Optimal - objective value 301209.38
Optimal objective 301209.3811 - 923 iterations time 0.092, Presolve 0.02
Total time (CPU seconds):       0.20   (Wallclock seconds):       0.17


[13]:

('ok', 'optimal')

And, as expected, the costs are the same as for the nodal market: 301 k€.

Now, we can plot both the market results of the 2 bidding zone market and the redispatch results:

[14]:

fig, axs = plt.subplots(
    1, 3, figsize=(20, 10), subplot_kw={"projection": ccrs.AlbersEqualArea()}
)

market = (
    n.generators_t.p[m.generators.index]
    .T.squeeze()
    .groupby(n.generators.bus)
    .sum()
    .div(2e4)
)
n.plot(ax=axs[0], bus_sizes=market, title="2 bidding zones market simulation")

redispatch_up = (
    n.generators_t.p.filter(like="ramp up")
    .T.squeeze()
    .groupby(n.generators.bus)
    .sum()
    .div(2e4)
)
n.plot(
    ax=axs[1], bus_sizes=redispatch_up, bus_colors="blue", title="Redispatch: ramp up"
)

redispatch_down = (
    n.generators_t.p.filter(like="ramp down")
    .T.squeeze()
    .groupby(n.generators.bus)
    .sum()
    .div(-2e4)
)
n.plot(
    ax=axs[2],
    bus_sizes=redispatch_down,
    bus_colors="red",
    title="Redispatch: ramp down / curtail",
);

/home/docs/checkouts/readthedocs.org/user_builds/pypsa/conda/v0.24.0/lib/python3.11/site-packages/cartopy/mpl/style.py:76: UserWarning: facecolor will have no effect as it has been defined as "never".
  warnings.warn('facecolor will have no effect as it has been '
../_images/examples_scigrid-redispatch_30_1.png

We can also read out the final dispatch of each generator:

[15]:

grouper = n.generators.index.str.split(" ramp", expand=True).get_level_values(0)

n.generators_t.p.groupby(grouper, axis=1).sum().squeeze()

[15]:

1 Gas                     0.000000
1 Hard Coal               0.000000
1 Solar                  11.326192
1 Wind Onshore            1.754375
100_220kV Solar          14.913326
                           ...
98 Wind Onshore          71.451646
99_220kV Gas              0.000000
99_220kV Hard Coal        0.000000
99_220kV Solar            8.246606
99_220kV Wind Onshore     3.432939
Name: 2011-01-01 12:00:00, Length: 1423, dtype: float64

Changing bidding strategies in redispatch market#

We can also formulate other bidding strategies or compensation mechanisms for the redispatch market.

For example, that ramping up a generator is twice as expensive.

[16]:

n.generators.loc[n.generators.index.str.contains("ramp up"), "marginal_cost"] *= 2

Or that generators need to be compensated for curtailing them or ramping them down at 50% of their marginal cost.

[17]:

n.generators.loc[n.generators.index.str.contains("ramp down"), "marginal_cost"] *= -0.5

In this way, the outcome should be more expensive than the ideal nodal market:

[18]:

n.optimize(solver_name=solver)

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 time: 0.15s
INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 5331 primals, 11649 duals
Objective: 4.79e+05
Solver model: not available
Solver message: Optimal - objective value 479003.12190570



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

command line - cbc -printingOptions all -import /tmp/linopy-problem-sa6vgbrq.lp -solve -solu /tmp/linopy-solve-86pyzflf.sol (default strategy 1)
Option for printingOptions changed from normal to all
Presolve 624 (-11025) rows, 1326 (-4005) columns and 4213 (-15344) elements
Perturbing problem by 0.001% of 4841.0401 - largest nonzero change 0.00053980221 ( 0.0024812559%) - largest zero change 0.0005285915
0  Obj 223384.43 Primal inf 1437635.3 (575)
87  Obj 223385.46 Primal inf 1968854.8 (535)
174  Obj 223387.25 Primal inf 758765.58 (522)
261  Obj 223389.5 Primal inf 873742.25 (466)
342  Obj 223421.63 Primal inf 2354311.3 (456)
429  Obj 230806.47 Primal inf 649143.58 (337)
493  Obj 231259.66 Primal inf 356589.93 (263)
563  Obj 231304.1 Primal inf 89728.959 (155)
632  Obj 231306.8 Primal inf 22873.786 (106)
719  Obj 320623.96 Primal inf 3340.972 (53)
806  Obj 479005.61
806  Obj 479003.13 Dual inf 0.00037331915 (7)
812  Obj 479003.12
Optimal - objective value 479003.12
After Postsolve, objective 479003.12, infeasibilities - dual 2851.6288 (92), primal 2.0051576e-05 (87)
Presolved model was optimal, full model needs cleaning up
0  Obj 479003.12 Primal inf 1.0448873e-06 (6) Dual inf 6.0000003e+08 (98)
End of values pass after 98 iterations
98  Obj 479003.12
Optimal - objective value 479003.12
Optimal objective 479003.1219 - 910 iterations time 0.112, Presolve 0.02
Total time (CPU seconds):       0.21   (Wallclock seconds):       0.17


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

[18]:

('ok', 'optimal')

Costs are now 502 k€ compared to 301 k€.