Note
You can download this example as a Jupyter notebook or start it in interactive mode.
Stochastic optimisation example#
This example optimises a single node with only a wind, solar, gas and lignite generator under uncertainty about the gas price.
In Stage 1 decisions are made about capacities of the generators while the gas price is unknown.
First we solve assuming knowledge about the gas price, then stochastically according to the probability distribution of gas prices.
We then show that the average total cost of the system of the stochastically optimised capacities is lower than the means of the solutions from the deterministically determined capacities.
Required data#
For this example, we need solar and wind generation time-series. For convenience, we will be fetching the time-series data directly from the renewables.ninja server. An arbitrary example of Germany’s data is retrieved.
The fetched files: - PV (1985-2016, SARAH) (6.37 MB) - Wind (Current fleet, onshore/offshore separate, MERRA-2) (13.93 MB)
Dependencies#
[1]:
from io import StringIO
import matplotlib.pyplot as plt
import pandas as pd
import requests
from linopy.expressions import merge
from xarray import DataArray
import pypsa
from pypsa.descriptors import (
get_bounds_pu,
nominal_attrs,
)
from pypsa.descriptors import get_switchable_as_dense as get_as_dense
from pypsa.optimization.common import reindex
[2]:
%matplotlib inline
Retrieve PV & Wind data#
[3]:
urls = {
"solar_pu": "https://www.renewables.ninja/country_downloads/DE/ninja_pv_country_DE_sarah_corrected.csv",
"wind_pu": "https://www.renewables.ninja/country_downloads/DE/ninja_wind_country_DE_current-merra-2_corrected.csv",
}
[4]:
def fetch_timeseries_data(url):
"""Fetch the timeseries data from the renewable.ninja server"""
response = requests.get(url)
response.raise_for_status() # Raise an error for bad responses
return pd.read_csv(
StringIO(response.text), skiprows=2, parse_dates=["time"], index_col="time"
)["national"]
[5]:
solar_pu = fetch_timeseries_data(urls["solar_pu"])
wind_pu = fetch_timeseries_data(urls["wind_pu"])
Major settings#
[6]:
scenarios = ["low", "med", "high"]
# this just determines the default scenario when building stochastic model
base_scenario = "low"
# in EUR/MWh_th
gas_prices = {"low": 40, "med": 70, "high": 100}
probability = {"low": 0.4, "med": 0.3, "high": 0.3}
[7]:
# years for weather data (solar is 1985-2015 inclusive, wind is 1980-2019)
year_start = 2015
year_end = 2015
# 1 is hourly, 3 is 3-hourly
frequency = 3
# Fixed load in MW
load = 1
# https://github.com/ERGO-Code/HiGHS
solver_name = "highs"
cts = ["DE"]
Prepare data#
[8]:
assumptions = pd.DataFrame(
columns=["FOM", "discount rate", "efficiency", "investment", "lifetime"],
index=["default", "onshore wind", "utility solar PV", "gas CCGT", "lignite"],
)
assumptions.at["default", "FOM"] = 3.0
assumptions.at["default", "discount rate"] = 0.03
assumptions.at["default", "lifetime"] = 25
assumptions.at["onshore wind", "investment"] = 2e6
assumptions.at["utility solar PV", "investment"] = 10e5
assumptions.at["gas CCGT", "investment"] = 7e5
assumptions.at["gas CCGT", "efficiency"] = 0.6
assumptions.at["lignite", "investment"] = 15e5
assumptions.at["lignite", "efficiency"] = 0.3
# fill defaults
assumptions = assumptions.fillna(
{
"FOM": assumptions.at["default", "FOM"],
"discount rate": assumptions.at["default", "discount rate"],
"lifetime": assumptions.at["default", "lifetime"],
}
)
def annuity(lifetime, rate):
if rate == 0.0:
return 1 / lifetime
else:
return rate / (1.0 - 1.0 / (1.0 + rate) ** lifetime)
# annualise investment costs, add FOM
assumptions["fixed"] = [
(annuity(v["lifetime"], v["discount rate"]) + v["FOM"] / 100.0) * v["investment"]
for i, v in assumptions.iterrows()
]
assumptions
/tmp/ipykernel_4469/322212998.py:19: FutureWarning: Downcasting object dtype arrays on .fillna, .ffill, .bfill is deprecated and will change in a future version. Call result.infer_objects(copy=False) instead. To opt-in to the future behavior, set `pd.set_option('future.no_silent_downcasting', True)`
assumptions = assumptions.fillna(
[8]:
| FOM | discount rate | efficiency | investment | lifetime | fixed | |
|---|---|---|---|---|---|---|
| default | 3.0 | 0.03 | NaN | NaN | 25 | NaN |
| onshore wind | 3.0 | 0.03 | NaN | 2000000.0 | 25 | 174855.742078 |
| utility solar PV | 3.0 | 0.03 | NaN | 1000000.0 | 25 | 87427.871039 |
| gas CCGT | 3.0 | 0.03 | 0.6 | 700000.0 | 25 | 61199.509727 |
| lignite | 3.0 | 0.03 | 0.3 | 1500000.0 | 25 | 131141.806559 |
Required functions#
[9]:
# prepare base network (without stochastic optimisation)
def prepare_network(cts, gas_price):
network = pypsa.Network()
snapshots = pd.date_range(
f"{year_start}-01-01",
f"{year_end}-12-31 23:00",
freq=str(frequency) + "H",
)
network.set_snapshots(snapshots)
network.snapshot_weightings = pd.Series(float(frequency), index=network.snapshots)
for ct in cts:
network.add("Bus", ct)
network.add("Load", ct, bus=ct, p_set=load)
network.add(
"Generator",
ct + " solar",
bus=ct,
p_max_pu=solar_pu,
p_nom_extendable=True,
marginal_cost=0.01, # Small cost to prefer curtailment to destroying energy in storage, solar curtails before wind
capital_cost=assumptions.at["utility solar PV", "fixed"],
)
network.add(
"Generator",
ct + " wind",
bus=ct,
p_max_pu=wind_pu,
p_nom_extendable=True,
marginal_cost=0.02, # Small cost to prefer curtailment to destroying energy in storage, solar curtails before wind
capital_cost=assumptions.at["onshore wind", "fixed"],
)
network.add(
"Generator",
ct + " gas",
bus=ct,
p_nom_extendable=True,
efficiency=assumptions.at["gas CCGT", "efficiency"],
marginal_cost=gas_price / assumptions.at["gas CCGT", "efficiency"],
capital_cost=assumptions.at["gas CCGT", "fixed"],
)
network.add(
"Generator",
ct + " lignite",
bus=ct,
p_nom_extendable=True,
efficiency=assumptions.at["lignite", "efficiency"],
marginal_cost=150,
capital_cost=assumptions.at["gas CCGT", "fixed"],
)
return network
[10]:
# add additional operational scenarios to the base model
def prepare_stochastic_model(n):
m = n.optimize.create_model()
nonbase_scenarios = scenarios.copy()
nonbase_scenarios.remove(base_scenario)
# we only have generators in this example, which simplifies things
c = "Generator"
sns = n.snapshots
attr = "p"
active = None
column = "bus"
sign = 1
ext_i = n.get_extendable_i(c)
min_pu, max_pu = map(DataArray, get_bounds_pu(n, c, sns, ext_i, attr))
capacity = n.model[f"{c}-{nominal_attrs[c]}"]
for scenario in nonbase_scenarios:
# add extra operational variables for each non-base scenario
dispatch = m.add_variables(
coords=m["Generator-p"].coords, name=f"Generator-p-{scenario}"
)
dispatch = reindex(dispatch, c, ext_i)
# add dispatch constraints
lhs = dispatch - max_pu * capacity # instead of the tuple formulation
m.add_constraints(lhs, "<=", 0, f"{c}-ext-{attr}-upper-{scenario}", active)
lhs = dispatch - min_pu * capacity
m.add_constraints(lhs, ">=", 0, f"{c}-ext-{attr}-lower-{scenario}", active)
# add nodal balance constraints
exprs = []
expr = DataArray(sign) * m[f"{c}-{attr}-{scenario}"]
buses = n.static(c)[column].rename("Bus")
expr = expr.groupby(
buses.to_xarray()
).sum() # for linopy >=0.2, see breaking changes log
exprs.append(expr)
lhs = merge(exprs).reindex(Bus=n.buses.index)
rhs = (
(-get_as_dense(n, "Load", "p_set", sns) * n.loads.sign)
.groupby(n.loads.bus, axis=1)
.sum()
.reindex(columns=n.buses.index, fill_value=0)
)
rhs.index.name = "snapshot"
rhs = DataArray(rhs)
mask = None
m.add_constraints(lhs, "=", rhs, f"Bus-nodal_balance-{scenario}", mask=mask)
# define the new objective
objective = []
weighting = n.snapshot_weightings.objective
weighting = weighting.loc[sns]
cost = (
get_as_dense(n, c, "marginal_cost", sns)
.loc[:, lambda ds: (ds != 0).all()]
.mul(weighting, axis=0)
)
for scenario in scenarios:
cost_modified = cost.copy()
if scenario == base_scenario:
name = f"{c}-{attr}"
else:
name = f"{c}-{attr}-{scenario}"
cost_modified["DE gas"] = (
cost_modified["DE gas"]
* gas_prices[scenario]
/ gas_prices[base_scenario]
)
operation = m[name].sel({"snapshot": sns, c: cost.columns})
objective.append((operation * (probability[scenario] * cost_modified)).sum())
ext_i = n.get_extendable_i(c)
cost = n.static(c)["capital_cost"][ext_i]
objective.append((capacity * cost).sum())
m.objective = merge(objective)
[11]:
# Check that network is created correctly:
# gas_price = 30
# n = prepare_network(cts,gas_price)
First solve capacities for each scenario deterministically#
[12]:
results = None
for scenario in scenarios:
gas_price = gas_prices[scenario]
n = prepare_network(cts, gas_price)
n.optimize(solver_name=solver_name)
if results is None:
results = pd.DataFrame(columns=n.generators.index)
results.index.name = "scenario"
results.loc[scenario] = n.generators.p_nom_opt
/tmp/ipykernel_4469/1365061072.py:5: FutureWarning: 'H' is deprecated and will be removed in a future version, please use 'h' instead.
snapshots = pd.date_range(
INFO:pypsa.components:Applying weightings to all columns of `snapshot_weightings`
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[12], line 6
3 for scenario in scenarios:
4 gas_price = gas_prices[scenario]
----> 6 n = prepare_network(cts, gas_price)
8 n.optimize(solver_name=solver_name)
10 if results is None:
Cell In[9], line 19, in prepare_network(cts, gas_price)
16 network.add("Bus", ct)
17 network.add("Load", ct, bus=ct, p_set=load)
---> 19 network.add(
20 "Generator",
21 ct + " solar",
22 bus=ct,
23 p_max_pu=solar_pu,
24 p_nom_extendable=True,
25 marginal_cost=0.01, # Small cost to prefer curtailment to destroying energy in storage, solar curtails before wind
26 capital_cost=assumptions.at["utility solar PV", "fixed"],
27 )
29 network.add(
30 "Generator",
31 ct + " wind",
(...)
36 capital_cost=assumptions.at["onshore wind", "fixed"],
37 )
39 network.add(
40 "Generator",
41 ct + " gas",
(...)
46 capital_cost=assumptions.at["gas CCGT", "fixed"],
47 )
File ~/checkouts/readthedocs.org/user_builds/pypsa/envs/v0.31.0/lib/python3.12/site-packages/pypsa/components.py:986, in Network.add(self, class_name, name, suffix, overwrite, **kwargs)
984 if isinstance(v, pd.Series) and single_component:
985 if not v.index.equals(self.snapshots):
--> 986 raise ValueError(msg.format(f"Series {k}", "network snapshots"))
987 elif isinstance(v, pd.Series):
988 # Cast names index to string + suffix
989 v = v.rename(
990 index=lambda s: str(s) if str(s).endswith(suffix) else s + suffix
991 )
ValueError: Series p_max_pu has an index which does not align with the passed network snapshots.
[13]:
results
Now solve the full problem stochastically#
[14]:
gas_price = gas_prices[base_scenario]
n = prepare_network(cts, gas_price)
prepare_stochastic_model(n)
n.optimize.solve_model(solver_name=solver_name)
/tmp/ipykernel_4469/1365061072.py:5: FutureWarning: 'H' is deprecated and will be removed in a future version, please use 'h' instead.
snapshots = pd.date_range(
INFO:pypsa.components:Applying weightings to all columns of `snapshot_weightings`
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[14], line 3
1 gas_price = gas_prices[base_scenario]
----> 3 n = prepare_network(cts, gas_price)
5 prepare_stochastic_model(n)
7 n.optimize.solve_model(solver_name=solver_name)
Cell In[9], line 19, in prepare_network(cts, gas_price)
16 network.add("Bus", ct)
17 network.add("Load", ct, bus=ct, p_set=load)
---> 19 network.add(
20 "Generator",
21 ct + " solar",
22 bus=ct,
23 p_max_pu=solar_pu,
24 p_nom_extendable=True,
25 marginal_cost=0.01, # Small cost to prefer curtailment to destroying energy in storage, solar curtails before wind
26 capital_cost=assumptions.at["utility solar PV", "fixed"],
27 )
29 network.add(
30 "Generator",
31 ct + " wind",
(...)
36 capital_cost=assumptions.at["onshore wind", "fixed"],
37 )
39 network.add(
40 "Generator",
41 ct + " gas",
(...)
46 capital_cost=assumptions.at["gas CCGT", "fixed"],
47 )
File ~/checkouts/readthedocs.org/user_builds/pypsa/envs/v0.31.0/lib/python3.12/site-packages/pypsa/components.py:986, in Network.add(self, class_name, name, suffix, overwrite, **kwargs)
984 if isinstance(v, pd.Series) and single_component:
985 if not v.index.equals(self.snapshots):
--> 986 raise ValueError(msg.format(f"Series {k}", "network snapshots"))
987 elif isinstance(v, pd.Series):
988 # Cast names index to string + suffix
989 v = v.rename(
990 index=lambda s: str(s) if str(s).endswith(suffix) else s + suffix
991 )
ValueError: Series p_max_pu has an index which does not align with the passed network snapshots.
[15]:
results.loc["stochastic"] = n.generators.p_nom_opt
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[15], line 1
----> 1 results.loc["stochastic"] = n.generators.p_nom_opt
NameError: name 'n' is not defined
[16]:
results
Now test each set of capacities against realisations of the gas price#
[17]:
for scenario in scenarios:
gas_price = gas_prices[scenario]
n = prepare_network(cts, gas_price)
n.generators.p_nom_extendable = False
for capacity_scenario in results.index:
n.generators.p_nom = results.loc[capacity_scenario, n.generators.index]
print(n.generators.p_nom)
n.optimize(solver_name=solver_name)
results.at[capacity_scenario, f"gas-p-{scenario}"] = n.generators_t.p[
"DE gas"
].sum()
results.at[capacity_scenario, f"lignite-p-{scenario}"] = n.generators_t.p[
"DE lignite"
].sum()
/tmp/ipykernel_4469/1365061072.py:5: FutureWarning: 'H' is deprecated and will be removed in a future version, please use 'h' instead.
snapshots = pd.date_range(
INFO:pypsa.components:Applying weightings to all columns of `snapshot_weightings`
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[17], line 3
1 for scenario in scenarios:
2 gas_price = gas_prices[scenario]
----> 3 n = prepare_network(cts, gas_price)
4 n.generators.p_nom_extendable = False
6 for capacity_scenario in results.index:
Cell In[9], line 19, in prepare_network(cts, gas_price)
16 network.add("Bus", ct)
17 network.add("Load", ct, bus=ct, p_set=load)
---> 19 network.add(
20 "Generator",
21 ct + " solar",
22 bus=ct,
23 p_max_pu=solar_pu,
24 p_nom_extendable=True,
25 marginal_cost=0.01, # Small cost to prefer curtailment to destroying energy in storage, solar curtails before wind
26 capital_cost=assumptions.at["utility solar PV", "fixed"],
27 )
29 network.add(
30 "Generator",
31 ct + " wind",
(...)
36 capital_cost=assumptions.at["onshore wind", "fixed"],
37 )
39 network.add(
40 "Generator",
41 ct + " gas",
(...)
46 capital_cost=assumptions.at["gas CCGT", "fixed"],
47 )
File ~/checkouts/readthedocs.org/user_builds/pypsa/envs/v0.31.0/lib/python3.12/site-packages/pypsa/components.py:986, in Network.add(self, class_name, name, suffix, overwrite, **kwargs)
984 if isinstance(v, pd.Series) and single_component:
985 if not v.index.equals(self.snapshots):
--> 986 raise ValueError(msg.format(f"Series {k}", "network snapshots"))
987 elif isinstance(v, pd.Series):
988 # Cast names index to string + suffix
989 v = v.rename(
990 index=lambda s: str(s) if str(s).endswith(suffix) else s + suffix
991 )
ValueError: Series p_max_pu has an index which does not align with the passed network snapshots.
[18]:
results
[19]:
for capacity_scenario in results.index:
for g in n.generators.index:
results.at[capacity_scenario, f"{g} CC"] = (
results.at[capacity_scenario, g] * n.generators.at[g, "capital_cost"]
)
for scenario in scenarios:
results.at[capacity_scenario, f"DE gas-{scenario} MC"] = (
n.snapshot_weightings.objective.mean()
* gas_prices[scenario]
/ n.generators.at["DE gas", "efficiency"]
* results.at[capacity_scenario, f"gas-p-{scenario}"]
)
results.at[capacity_scenario, f"DE lignite-{scenario} MC"] = (
n.snapshot_weightings.objective.mean()
* n.generators.at["DE lignite", "marginal_cost"]
* results.at[capacity_scenario, f"lignite-p-{scenario}"]
)
results.at[capacity_scenario, "DE gas-mean MC"] = sum(
[
probability[scenario]
* results.at[capacity_scenario, f"DE gas-{scenario} MC"]
for scenario in scenarios
]
)
results.at[capacity_scenario, "DE lignite-mean MC"] = sum(
[
probability[scenario]
* results.at[capacity_scenario, f"DE lignite-{scenario} MC"]
for scenario in scenarios
]
)
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
Cell In[19], line 1
----> 1 for capacity_scenario in results.index:
2 for g in n.generators.index:
3 results.at[capacity_scenario, f"{g} CC"] = (
4 results.at[capacity_scenario, g] * n.generators.at[g, "capital_cost"]
5 )
AttributeError: 'NoneType' object has no attribute 'index'
[20]:
fig, axes = plt.subplots(1, len(results.index), figsize=(len(results.index) * 4, 4))
colors = {
"wind": "b",
"solar": "y",
"lignite": "black",
"gas": "brown",
"gas MC": "orange",
"lignite MC": "gray",
}
# fig.suptitle('Horizontally stacked subplots')
for i, capacity_scenario in enumerate(results.index):
ax = axes[i]
df = pd.DataFrame(index=scenarios + ["mean"])
for tech in ["solar", "wind", "gas", "lignite"]:
df[tech] = results.at[capacity_scenario, f"DE {tech} CC"]
for scenario in scenarios + ["mean"]:
df.at[scenario, "gas MC"] = results.at[
capacity_scenario, f"DE gas-{scenario} MC"
]
df.at[scenario, "lignite MC"] = results.at[
capacity_scenario, f"DE lignite-{scenario} MC"
]
df.plot(kind="bar", stacked=True, ax=ax, color=colors)
ax.set_title(f"capacity scenario {capacity_scenario}")
ax.legend(loc="upper left")
ax.set_ylim([0, 2.5e6])
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
Cell In[20], line 1
----> 1 fig, axes = plt.subplots(1, len(results.index), figsize=(len(results.index) * 4, 4))
3 colors = {
4 "wind": "b",
5 "solar": "y",
(...)
9 "lignite MC": "gray",
10 }
12 # fig.suptitle('Horizontally stacked subplots')
AttributeError: 'NoneType' object has no attribute 'index'
[21]:
fig, ax = plt.subplots(1, 1, figsize=(4, 4))
df = (
results[
[
"DE solar CC",
"DE wind CC",
"DE gas CC",
"DE lignite CC",
"DE gas-mean MC",
"DE lignite-mean MC",
]
]
.rename(columns=lambda x: x[3:-3])
.rename(columns={"gas-mean": "gas MC", "lignite-mean": "lignite MC"})
)
df.plot(kind="bar", stacked=True, ax=ax, color=colors)
ax.set_xlabel("capacity scenario")
ax.set_title("means of results")
ax.set_ylim([0, 2e6])
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
Cell In[21], line 4
1 fig, ax = plt.subplots(1, 1, figsize=(4, 4))
3 df = (
----> 4 results[
5 [
6 "DE solar CC",
7 "DE wind CC",
8 "DE gas CC",
9 "DE lignite CC",
10 "DE gas-mean MC",
11 "DE lignite-mean MC",
12 ]
13 ]
14 .rename(columns=lambda x: x[3:-3])
15 .rename(columns={"gas-mean": "gas MC", "lignite-mean": "lignite MC"})
16 )
18 df.plot(kind="bar", stacked=True, ax=ax, color=colors)
20 ax.set_xlabel("capacity scenario")
TypeError: 'NoneType' object is not subscriptable
Analysis of a Stochastic Solution#
The Expected costs of ignoring uncertainty (ECIU)#
in some literature also defined as the Value of Stochastic Solution (VSS). Can be used interchangeably.
The natural question to ask is how much difference it really makes to the quality of the decisions reached if I use a stochastic problem instead of a deterministic problem?
The ECIU measures the value of using a stochastic model (or the expected costs of ignoring uncertainty when using a deterministic model).
[22]:
portfolios = pd.DataFrame()
costs = pd.Series()
Define the naive problem (usually – the expected value problem (EV))#
[23]:
# can be anything (e.g., the 'med' scenario). A texbook way is to take expected value of uncertain parameter.
naive_scenario = sum(pd.Series(gas_prices) * pd.Series(probability))
naive_scenario
# naive_scenario = gas_prices["med"]
[23]:
67.0
solve naive problem (deterministic)#
[24]:
scenario = "naive" # naive problem (in literature often EVP for Expected Value Problem, if the naive assumption is the expected value)
gas_price = naive_scenario
n = prepare_network(cts, gas_price)
n.optimize(solver_name=solver_name)
portfolios[scenario] = n.generators.p_nom_opt
costs[scenario] = n.objective
/tmp/ipykernel_4469/1365061072.py:5: FutureWarning: 'H' is deprecated and will be removed in a future version, please use 'h' instead.
snapshots = pd.date_range(
INFO:pypsa.components:Applying weightings to all columns of `snapshot_weightings`
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[24], line 4
1 scenario = "naive" # naive problem (in literature often EVP for Expected Value Problem, if the naive assumption is the expected value)
2 gas_price = naive_scenario
----> 4 n = prepare_network(cts, gas_price)
6 n.optimize(solver_name=solver_name)
8 portfolios[scenario] = n.generators.p_nom_opt
Cell In[9], line 19, in prepare_network(cts, gas_price)
16 network.add("Bus", ct)
17 network.add("Load", ct, bus=ct, p_set=load)
---> 19 network.add(
20 "Generator",
21 ct + " solar",
22 bus=ct,
23 p_max_pu=solar_pu,
24 p_nom_extendable=True,
25 marginal_cost=0.01, # Small cost to prefer curtailment to destroying energy in storage, solar curtails before wind
26 capital_cost=assumptions.at["utility solar PV", "fixed"],
27 )
29 network.add(
30 "Generator",
31 ct + " wind",
(...)
36 capital_cost=assumptions.at["onshore wind", "fixed"],
37 )
39 network.add(
40 "Generator",
41 ct + " gas",
(...)
46 capital_cost=assumptions.at["gas CCGT", "fixed"],
47 )
File ~/checkouts/readthedocs.org/user_builds/pypsa/envs/v0.31.0/lib/python3.12/site-packages/pypsa/components.py:986, in Network.add(self, class_name, name, suffix, overwrite, **kwargs)
984 if isinstance(v, pd.Series) and single_component:
985 if not v.index.equals(self.snapshots):
--> 986 raise ValueError(msg.format(f"Series {k}", "network snapshots"))
987 elif isinstance(v, pd.Series):
988 # Cast names index to string + suffix
989 v = v.rename(
990 index=lambda s: str(s) if str(s).endswith(suffix) else s + suffix
991 )
ValueError: Series p_max_pu has an index which does not align with the passed network snapshots.
[25]:
# pd.set_option("display.precision", 10)
portfolios
# costs
[25]:
solve stochastic problem#
[26]:
scenario = "SP" # SP for Stochastic Problem
gas_price = gas_prices[base_scenario]
n = prepare_network(cts, gas_price)
prepare_stochastic_model(n)
n.optimize.solve_model(solver_name=solver_name)
portfolios[scenario] = n.generators.p_nom_opt
costs[scenario] = n.objective
/tmp/ipykernel_4469/1365061072.py:5: FutureWarning: 'H' is deprecated and will be removed in a future version, please use 'h' instead.
snapshots = pd.date_range(
INFO:pypsa.components:Applying weightings to all columns of `snapshot_weightings`
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[26], line 4
1 scenario = "SP" # SP for Stochastic Problem
2 gas_price = gas_prices[base_scenario]
----> 4 n = prepare_network(cts, gas_price)
5 prepare_stochastic_model(n)
7 n.optimize.solve_model(solver_name=solver_name)
Cell In[9], line 19, in prepare_network(cts, gas_price)
16 network.add("Bus", ct)
17 network.add("Load", ct, bus=ct, p_set=load)
---> 19 network.add(
20 "Generator",
21 ct + " solar",
22 bus=ct,
23 p_max_pu=solar_pu,
24 p_nom_extendable=True,
25 marginal_cost=0.01, # Small cost to prefer curtailment to destroying energy in storage, solar curtails before wind
26 capital_cost=assumptions.at["utility solar PV", "fixed"],
27 )
29 network.add(
30 "Generator",
31 ct + " wind",
(...)
36 capital_cost=assumptions.at["onshore wind", "fixed"],
37 )
39 network.add(
40 "Generator",
41 ct + " gas",
(...)
46 capital_cost=assumptions.at["gas CCGT", "fixed"],
47 )
File ~/checkouts/readthedocs.org/user_builds/pypsa/envs/v0.31.0/lib/python3.12/site-packages/pypsa/components.py:986, in Network.add(self, class_name, name, suffix, overwrite, **kwargs)
984 if isinstance(v, pd.Series) and single_component:
985 if not v.index.equals(self.snapshots):
--> 986 raise ValueError(msg.format(f"Series {k}", "network snapshots"))
987 elif isinstance(v, pd.Series):
988 # Cast names index to string + suffix
989 v = v.rename(
990 index=lambda s: str(s) if str(s).endswith(suffix) else s + suffix
991 )
ValueError: Series p_max_pu has an index which does not align with the passed network snapshots.
[27]:
portfolios
[27]:
Solve stochastic problem constrained by the naive solution#
[28]:
scenario = "SP-constrained"
gas_price = gas_prices[base_scenario]
n = prepare_network(cts, gas_price)
prepare_stochastic_model(n)
n.generators.p_nom_extendable = False
n.generators.p_nom = portfolios.loc[n.generators.index, "naive"]
# n.generators.T
n.optimize.solve_model(solver_name=solver_name)
/tmp/ipykernel_4469/1365061072.py:5: FutureWarning: 'H' is deprecated and will be removed in a future version, please use 'h' instead.
snapshots = pd.date_range(
INFO:pypsa.components:Applying weightings to all columns of `snapshot_weightings`
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[28], line 4
1 scenario = "SP-constrained"
3 gas_price = gas_prices[base_scenario]
----> 4 n = prepare_network(cts, gas_price)
5 prepare_stochastic_model(n)
7 n.generators.p_nom_extendable = False
Cell In[9], line 19, in prepare_network(cts, gas_price)
16 network.add("Bus", ct)
17 network.add("Load", ct, bus=ct, p_set=load)
---> 19 network.add(
20 "Generator",
21 ct + " solar",
22 bus=ct,
23 p_max_pu=solar_pu,
24 p_nom_extendable=True,
25 marginal_cost=0.01, # Small cost to prefer curtailment to destroying energy in storage, solar curtails before wind
26 capital_cost=assumptions.at["utility solar PV", "fixed"],
27 )
29 network.add(
30 "Generator",
31 ct + " wind",
(...)
36 capital_cost=assumptions.at["onshore wind", "fixed"],
37 )
39 network.add(
40 "Generator",
41 ct + " gas",
(...)
46 capital_cost=assumptions.at["gas CCGT", "fixed"],
47 )
File ~/checkouts/readthedocs.org/user_builds/pypsa/envs/v0.31.0/lib/python3.12/site-packages/pypsa/components.py:986, in Network.add(self, class_name, name, suffix, overwrite, **kwargs)
984 if isinstance(v, pd.Series) and single_component:
985 if not v.index.equals(self.snapshots):
--> 986 raise ValueError(msg.format(f"Series {k}", "network snapshots"))
987 elif isinstance(v, pd.Series):
988 # Cast names index to string + suffix
989 v = v.rename(
990 index=lambda s: str(s) if str(s).endswith(suffix) else s + suffix
991 )
ValueError: Series p_max_pu has an index which does not align with the passed network snapshots.
[29]:
# don't forget to add the capital costs of the (fixed) generators portfolio
c = "Generator"
ext_i = portfolios["naive"].index
cost = n.static(c)["capital_cost"][ext_i]
cost_of_portfolio = (n.generators.p_nom * cost).sum()
n.objective += cost_of_portfolio
n.objective
---------------------------------------------------------------------------
KeyError Traceback (most recent call last)
Cell In[29], line 3
1 # don't forget to add the capital costs of the (fixed) generators portfolio
2 c = "Generator"
----> 3 ext_i = portfolios["naive"].index
4 cost = n.static(c)["capital_cost"][ext_i]
5 cost_of_portfolio = (n.generators.p_nom * cost).sum()
File ~/checkouts/readthedocs.org/user_builds/pypsa/envs/v0.31.0/lib/python3.12/site-packages/pandas/core/frame.py:4102, in DataFrame.__getitem__(self, key)
4100 if self.columns.nlevels > 1:
4101 return self._getitem_multilevel(key)
-> 4102 indexer = self.columns.get_loc(key)
4103 if is_integer(indexer):
4104 indexer = [indexer]
File ~/checkouts/readthedocs.org/user_builds/pypsa/envs/v0.31.0/lib/python3.12/site-packages/pandas/core/indexes/range.py:417, in RangeIndex.get_loc(self, key)
415 raise KeyError(key) from err
416 if isinstance(key, Hashable):
--> 417 raise KeyError(key)
418 self._check_indexing_error(key)
419 raise KeyError(key)
KeyError: 'naive'
[30]:
portfolios[scenario] = (
n.generators.p_nom
) # just a fixed copy of naive problem's solution
costs[scenario] = (
n.objective
) # must be >= than the stochastic solution's costs, because you do dispatch with the suboptimal first-stage decisions
costs
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[30], line 2
1 portfolios[scenario] = (
----> 2 n.generators.p_nom
3 ) # just a fixed copy of naive problem's solution
4 costs[scenario] = (
5 n.objective
6 ) # must be >= than the stochastic solution's costs, because you do dispatch with the suboptimal first-stage decisions
8 costs
NameError: name 'n' is not defined
Compute ECIU#
[31]:
# ECIU (or VSS) in M euro
eciu = (costs["SP-constrained"] - costs["SP"]) / 1e6
# ECIU in % of stochastic solution
eciu_pp = eciu / (costs["SP"] / 1e6) * 100
print(
f"ECIU: {round(eciu, 3)} Meuro \nwhich is {round(eciu_pp)}% of stochastic solution's costs"
)
---------------------------------------------------------------------------
KeyError Traceback (most recent call last)
Cell In[31], line 2
1 # ECIU (or VSS) in M euro
----> 2 eciu = (costs["SP-constrained"] - costs["SP"]) / 1e6
3 # ECIU in % of stochastic solution
4 eciu_pp = eciu / (costs["SP"] / 1e6) * 100
File ~/checkouts/readthedocs.org/user_builds/pypsa/envs/v0.31.0/lib/python3.12/site-packages/pandas/core/series.py:1121, in Series.__getitem__(self, key)
1118 return self._values[key]
1120 elif key_is_scalar:
-> 1121 return self._get_value(key)
1123 # Convert generator to list before going through hashable part
1124 # (We will iterate through the generator there to check for slices)
1125 if is_iterator(key):
File ~/checkouts/readthedocs.org/user_builds/pypsa/envs/v0.31.0/lib/python3.12/site-packages/pandas/core/series.py:1237, in Series._get_value(self, label, takeable)
1234 return self._values[label]
1236 # Similar to Index.get_value, but we do not fall back to positional
-> 1237 loc = self.index.get_loc(label)
1239 if is_integer(loc):
1240 return self._values[loc]
File ~/checkouts/readthedocs.org/user_builds/pypsa/envs/v0.31.0/lib/python3.12/site-packages/pandas/core/indexes/range.py:417, in RangeIndex.get_loc(self, key)
415 raise KeyError(key) from err
416 if isinstance(key, Hashable):
--> 417 raise KeyError(key)
418 self._check_indexing_error(key)
419 raise KeyError(key)
KeyError: 'SP-constrained'
The Expected Value of Perfect Information (EVPI)#
If system planner knew at the first stage which scenario will play out, it could optimize an expansion plan (i.e. that results in lower cost) for that scenario.
The expected value (and the corresponding mathematical problem) of such solution is denoted in the literature as „wait-and-see” solution (or wait-and-see (WS) problem).
The difference between the (probability-weighted) wait-and-see solutions and the here-and-now (stochastic) solution represents the added value of information about the future (i.e., the expected profit).
modelling perspective: How much the expected costs could be reduced if system planner in the first stage knew exactly which scenario would happen?
economic perspective: An upper bound to the amount that should be paid for improved forecasts.
[32]:
portfolios = pd.DataFrame()
costs = pd.Series()
Solve Wait-and-See problems#
where Wait-and-See (WS) is a standard textbook name for individual determinic problem (i.e. running a single scenario).
[33]:
for scenario in scenarios:
gas_price = gas_prices[scenario]
n = prepare_network(cts, gas_price)
n.optimize(solver_name=solver_name)
if results is None:
results = pd.DataFrame(columns=n.generators.index)
results.index.name = "scenario"
portfolios[scenario] = n.generators.p_nom_opt
costs[scenario] = n.objective
/tmp/ipykernel_4469/1365061072.py:5: FutureWarning: 'H' is deprecated and will be removed in a future version, please use 'h' instead.
snapshots = pd.date_range(
INFO:pypsa.components:Applying weightings to all columns of `snapshot_weightings`
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[33], line 3
1 for scenario in scenarios:
2 gas_price = gas_prices[scenario]
----> 3 n = prepare_network(cts, gas_price)
5 n.optimize(solver_name=solver_name)
7 if results is None:
Cell In[9], line 19, in prepare_network(cts, gas_price)
16 network.add("Bus", ct)
17 network.add("Load", ct, bus=ct, p_set=load)
---> 19 network.add(
20 "Generator",
21 ct + " solar",
22 bus=ct,
23 p_max_pu=solar_pu,
24 p_nom_extendable=True,
25 marginal_cost=0.01, # Small cost to prefer curtailment to destroying energy in storage, solar curtails before wind
26 capital_cost=assumptions.at["utility solar PV", "fixed"],
27 )
29 network.add(
30 "Generator",
31 ct + " wind",
(...)
36 capital_cost=assumptions.at["onshore wind", "fixed"],
37 )
39 network.add(
40 "Generator",
41 ct + " gas",
(...)
46 capital_cost=assumptions.at["gas CCGT", "fixed"],
47 )
File ~/checkouts/readthedocs.org/user_builds/pypsa/envs/v0.31.0/lib/python3.12/site-packages/pypsa/components.py:986, in Network.add(self, class_name, name, suffix, overwrite, **kwargs)
984 if isinstance(v, pd.Series) and single_component:
985 if not v.index.equals(self.snapshots):
--> 986 raise ValueError(msg.format(f"Series {k}", "network snapshots"))
987 elif isinstance(v, pd.Series):
988 # Cast names index to string + suffix
989 v = v.rename(
990 index=lambda s: str(s) if str(s).endswith(suffix) else s + suffix
991 )
ValueError: Series p_max_pu has an index which does not align with the passed network snapshots.
compute the expected value of wait-and-see scenario costs#
[34]:
ws = sum(costs * pd.Series(probability))
solve stochastic problem#
[35]:
scenario = "SP" # SP for Stochastic Problem
gas_price = gas_prices[base_scenario]
n = prepare_network(cts, gas_price)
prepare_stochastic_model(n)
n.optimize.solve_model(solver_name=solver_name)
portfolios[scenario] = n.generators.p_nom_opt
costs[scenario] = n.objective
/tmp/ipykernel_4469/1365061072.py:5: FutureWarning: 'H' is deprecated and will be removed in a future version, please use 'h' instead.
snapshots = pd.date_range(
INFO:pypsa.components:Applying weightings to all columns of `snapshot_weightings`
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[35], line 4
1 scenario = "SP" # SP for Stochastic Problem
2 gas_price = gas_prices[base_scenario]
----> 4 n = prepare_network(cts, gas_price)
5 prepare_stochastic_model(n)
7 n.optimize.solve_model(solver_name=solver_name)
Cell In[9], line 19, in prepare_network(cts, gas_price)
16 network.add("Bus", ct)
17 network.add("Load", ct, bus=ct, p_set=load)
---> 19 network.add(
20 "Generator",
21 ct + " solar",
22 bus=ct,
23 p_max_pu=solar_pu,
24 p_nom_extendable=True,
25 marginal_cost=0.01, # Small cost to prefer curtailment to destroying energy in storage, solar curtails before wind
26 capital_cost=assumptions.at["utility solar PV", "fixed"],
27 )
29 network.add(
30 "Generator",
31 ct + " wind",
(...)
36 capital_cost=assumptions.at["onshore wind", "fixed"],
37 )
39 network.add(
40 "Generator",
41 ct + " gas",
(...)
46 capital_cost=assumptions.at["gas CCGT", "fixed"],
47 )
File ~/checkouts/readthedocs.org/user_builds/pypsa/envs/v0.31.0/lib/python3.12/site-packages/pypsa/components.py:986, in Network.add(self, class_name, name, suffix, overwrite, **kwargs)
984 if isinstance(v, pd.Series) and single_component:
985 if not v.index.equals(self.snapshots):
--> 986 raise ValueError(msg.format(f"Series {k}", "network snapshots"))
987 elif isinstance(v, pd.Series):
988 # Cast names index to string + suffix
989 v = v.rename(
990 index=lambda s: str(s) if str(s).endswith(suffix) else s + suffix
991 )
ValueError: Series p_max_pu has an index which does not align with the passed network snapshots.
Compute EVPI#
[36]:
# EVPI in M euro
evpi = (
costs["SP"] - ws
) / 1e6 # must be >=0 because improved information cannot make the decision maker worse
# ECIU in % of stochastic solution
evpi_pp = evpi / (costs["SP"] / 1e6) * 100
print(
f"EVPI: {round(evpi, 3)} Meuro \nwhich is {round(evpi_pp)}% of stochastic solution's costs"
)
---------------------------------------------------------------------------
KeyError Traceback (most recent call last)
Cell In[36], line 3
1 # EVPI in M euro
2 evpi = (
----> 3 costs["SP"] - ws
4 ) / 1e6 # must be >=0 because improved information cannot make the decision maker worse
5 # ECIU in % of stochastic solution
6 evpi_pp = evpi / (costs["SP"] / 1e6) * 100
File ~/checkouts/readthedocs.org/user_builds/pypsa/envs/v0.31.0/lib/python3.12/site-packages/pandas/core/series.py:1121, in Series.__getitem__(self, key)
1118 return self._values[key]
1120 elif key_is_scalar:
-> 1121 return self._get_value(key)
1123 # Convert generator to list before going through hashable part
1124 # (We will iterate through the generator there to check for slices)
1125 if is_iterator(key):
File ~/checkouts/readthedocs.org/user_builds/pypsa/envs/v0.31.0/lib/python3.12/site-packages/pandas/core/series.py:1237, in Series._get_value(self, label, takeable)
1234 return self._values[label]
1236 # Similar to Index.get_value, but we do not fall back to positional
-> 1237 loc = self.index.get_loc(label)
1239 if is_integer(loc):
1240 return self._values[loc]
File ~/checkouts/readthedocs.org/user_builds/pypsa/envs/v0.31.0/lib/python3.12/site-packages/pandas/core/indexes/range.py:417, in RangeIndex.get_loc(self, key)
415 raise KeyError(key) from err
416 if isinstance(key, Hashable):
--> 417 raise KeyError(key)
418 self._check_indexing_error(key)
419 raise KeyError(key)
KeyError: 'SP'
Comparing the ECIU and EVPI metrics#
ECIU: an investment decision is made when uncertainty is ignored. The ECIU is the additional expected cost of assuming that future is certain.
EVPI: an investment decision is made after uncertainty is removed. The EVPI is the expected cost of being uncertain about the future.