Screening Curves¶

Compute the long-term equilibrium power plant investment for a given load duration curve ($1000-1000z$ for $z \in [0,1]$) and a given set of generator investment options.

To read up on the theory behind screenin curves, see this lecture.

In [2]:

import numpy as np
import pandas as pd

import pypsa

import numpy as np import pandas as pd import pypsa

Generator marginal (m) and capital (c) costs are given in EUR/MWh - numbers chosen for simple answer.

In [3]:

generators = {
    "coal": {"m": 2, "c": 15},
    "gas": {"m": 12, "c": 10},
    "load-shedding": {"m": 1012, "c": 0},
}

generators = { "coal": {"m": 2, "c": 15}, "gas": {"m": 12, "c": 10}, "load-shedding": {"m": 1012, "c": 0}, }

The screening curve intersections are at 0.01 and 0.5.

In [4]:

x = np.linspace(0, 1, 101)
df = pd.DataFrame({k: v["c"] + x * v["m"] for k, v in generators.items()}, index=x)
df.plot(ylim=[0, 50], title="Screening Curve");

x = np.linspace(0, 1, 101) df = pd.DataFrame({k: v["c"] + x * v["m"] for k, v in generators.items()}, index=x) df.plot(ylim=[0, 50], title="Screening Curve");

In [5]:

n = pypsa.Network()

num_snapshots = 1001
n.snapshots = np.linspace(0, 1, num_snapshots)
n.snapshot_weightings = n.snapshot_weightings / num_snapshots

n.add("Bus", name="bus")

n.add("Load", name="load", bus="bus", p_set=1000 - 1000 * n.snapshots.values)

for gen in generators:
    n.add(
        "Generator",
        name=gen,
        bus="bus",
        p_nom_extendable=True,
        marginal_cost=generators[gen]["m"],
        capital_cost=generators[gen]["c"],
    )

n = pypsa.Network() num_snapshots = 1001 n.snapshots = np.linspace(0, 1, num_snapshots) n.snapshot_weightings = n.snapshot_weightings / num_snapshots n.add("Bus", name="bus") n.add("Load", name="load", bus="bus", p_set=1000 - 1000 * n.snapshots.values) for gen in generators: n.add( "Generator", name=gen, bus="bus", p_nom_extendable=True, marginal_cost=generators[gen]["m"], capital_cost=generators[gen]["c"], )

In [6]:

n.loads_t.p_set.plot.area(title="Load Duration Curve", ylabel="MW")

n.loads_t.p_set.plot.area(title="Load Duration Curve", ylabel="MW")

Out[6]:

<Axes: title={'center': 'Load Duration Curve'}, xlabel='snapshot', ylabel='MW'>

In [7]:

n.optimize()
n.objective

n.optimize() n.objective

/tmp/ipykernel_3927/542537002.py:1: FutureWarning: The default value of `include_objective_constant` will change from True to False in version 2.0. Set `include_objective_constant` explicitly to suppress this warning. Using False improves LP numerical conditioning by not including the objective constant as a variable.
  n.optimize()
WARNING:pypsa.consistency:The following buses have carriers which are not defined. Run n.sanitize() to add them. Components with undefined carriers:
Index(['bus'], dtype='object', name='name')

INFO:linopy.model: Solve problem using Highs solver

INFO:linopy.model:Solver options:
 - log_to_console: False

INFO:linopy.io: Writing time: 0.04s

INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 3006 primals, 7010 duals
Objective: 1.47e+04
Solver model: available
Solver message: Optimal

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-ext-p-lower, Generator-ext-p-upper were not assigned to the network.

Out[7]:

14706.193806194067

The capacity is set by total electricity required.

NB: No load shedding since all prices are below 10 000.

In [8]:

n.generators.p_nom_opt.round(2)

n.generators.p_nom_opt.round(2)

Out[8]:

name
coal             500.0
gas              490.0
load-shedding     10.0
Name: p_nom_opt, dtype: float64

In [9]:

n.buses_t.marginal_price.plot(title="Price Duration Curve")

n.buses_t.marginal_price.plot(title="Price Duration Curve")

Out[9]:

<Axes: title={'center': 'Price Duration Curve'}, xlabel='snapshot'>

The prices correspond either to VOLL (1012) for first 0.01 or the marginal costs (12 for 0.49 and 2 for 0.5)

Except for (infinitesimally small) points at the screening curve intersections, which correspond to changing the load duration near the intersection, so that capacity changes. This explains 7 = (12+10 - 15) (replacing coal with gas) and 22 = (12+10) (replacing load-shedding with gas).

Note: What remains unclear is what is causing l = 0... it should be 2.

In [10]:

n.buses_t.marginal_price.round(2).sum(axis=1).value_counts()

n.buses_t.marginal_price.round(2).sum(axis=1).value_counts()

Out[10]:

2.0       499
12.0      489
1012.0     10
22.0        1
7.0         1
0.0         1
Name: count, dtype: int64

In [11]:

n.generators_t.p.plot(ylim=[0, 600], title="Generation Dispatch")

n.generators_t.p.plot(ylim=[0, 600], title="Generation Dispatch")

Out[11]:

<Axes: title={'center': 'Generation Dispatch'}, xlabel='snapshot'>

Demonstrate zero-profit condition.

1. The total cost is given by

In [12]:

weights = n.snapshot_weightings.objective
(
    n.generators.p_nom_opt * n.generators.capital_cost
    + weights @ n.generators_t.p * n.generators.marginal_cost
)

weights = n.snapshot_weightings.objective ( n.generators.p_nom_opt * n.generators.capital_cost + weights @ n.generators_t.p * n.generators.marginal_cost )

Out[12]:

name
coal             8249.750250
gas              6400.839161
load-shedding      55.604396
dtype: float64

2. The total revenue by

In [13]:

weights @ n.generators_t.p.mul(n.buses_t.marginal_price["bus"], axis=0)

weights @ n.generators_t.p.mul(n.buses_t.marginal_price["bus"], axis=0)

Out[13]:

name
coal             8249.750250
gas              6400.839161
load-shedding      55.604396
Name: objective, dtype: float64

Now, take the capacities from the above long-term equilibrium, then disallow expansion.

Show that the resulting market prices are identical.

This holds in this example, but does not necessarily hold and breaks down in some circumstances (for example, when there is a lot of storage and inter-temporal shifting).

In [14]:

n.generators.p_nom_extendable = False
n.generators.p_nom = n.generators.p_nom_opt

n.generators.p_nom_extendable = False n.generators.p_nom = n.generators.p_nom_opt

In [15]:

n.optimize();

n.optimize();

/tmp/ipykernel_3927/1450685117.py:1: FutureWarning: The default value of `include_objective_constant` will change from True to False in version 2.0. Set `include_objective_constant` explicitly to suppress this warning. Using False improves LP numerical conditioning by not including the objective constant as a variable.
  n.optimize();
WARNING:pypsa.consistency:The following buses have carriers which are not defined. Run n.sanitize() to add them. Components with undefined carriers:
Index(['bus'], dtype='object', name='name')

WARNING:pypsa.consistency:The following sub_networks have carriers which are not defined. Run n.sanitize() to add them. Components with undefined carriers:
Index(['0'], dtype='object', name='name')

INFO:linopy.model: Solve problem using Highs solver

INFO:linopy.model:Solver options:
 - log_to_console: False

INFO:linopy.io: Writing time: 0.02s

INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 3003 primals, 7007 duals
Objective: 2.31e+03
Solver model: available
Solver message: Optimal

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper were not assigned to the network.

In [16]:

n.buses_t.marginal_price.plot(title="Price Duration Curve")

n.buses_t.marginal_price.plot(title="Price Duration Curve")

Out[16]:

<Axes: title={'center': 'Price Duration Curve'}, xlabel='snapshot'>

In [17]:

n.buses_t.marginal_price.sum(axis=1).value_counts()

n.buses_t.marginal_price.sum(axis=1).value_counts()

Out[17]:

2.0       500
12.0      490
1012.0     10
0.0         1
Name: count, dtype: int64

Demonstrate zero-profit condition. Differences are due to singular times, see above, not a problem

1. Total costs

In [18]:

(
    n.generators.p_nom * n.generators.capital_cost
    + weights @ n.generators_t.p * n.generators.marginal_cost
)

( n.generators.p_nom * n.generators.capital_cost + weights @ n.generators_t.p * n.generators.marginal_cost )

Out[18]:

name
coal             8249.750250
gas              6400.839161
load-shedding      55.604396
dtype: float64

2. Total revenue

In [19]:

weights @ n.generators_t.p.mul(n.buses_t.marginal_price["bus"], axis=0)

weights @ n.generators_t.p.mul(n.buses_t.marginal_price["bus"], axis=0)

Out[19]:

name
coal             8242.257742
gas              6395.944056
load-shedding      55.604396
Name: objective, dtype: float64