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)

See: https://www.renewables.ninja/

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_4603/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.df(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.df(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_4603/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`
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['DE'], dtype='object', name='Bus')
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['DE'], dtype='object', name='Bus')
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 5/5 [00:00<00:00, 33.49it/s]
Writing continuous variables.: 100%|██████████| 2/2 [00:00<00:00, 69.33it/s]
INFO:linopy.io: Writing time: 0.21s
INFO:linopy.solvers:Log file at /tmp/highs.log
INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 11684 primals, 26284 duals
Objective: 6.45e+05
Solver model: available
Solver message: optimal

Running HiGHS 1.7.2 (git hash: 184e327): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e-04, 1e+00]
  Cost   [3e-02, 2e+05]
  Bound  [0e+00, 0e+00]
  RHS    [1e+00, 1e+00]
Presolving model
13239 rows, 10323 cols, 30957 nonzeros  0s
13239 rows, 10323 cols, 30957 nonzeros  0s
Presolve : Reductions: rows 13239(-13045); columns 10323(-1361); elements 30957(-14406)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     0.0000000000e+00 Pr: 2920(34348) 0s
       7405     6.4519950973e+05 Pr: 0(0) 0s
Solving the original LP from the solution after postsolve
Model   status      : Optimal
Simplex   iterations: 7405
Objective value     :  6.4519950973e+05
HiGHS run time      :          0.12
Writing the solution to /tmp/linopy-solve-kwnhh5fi.sol

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-ext-p-lower, Generator-ext-p-upper were not assigned to the network.
/tmp/ipykernel_4603/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`
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['DE'], dtype='object', name='Bus')
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['DE'], dtype='object', name='Bus')
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 5/5 [00:00<00:00, 33.85it/s]
Writing continuous variables.: 100%|██████████| 2/2 [00:00<00:00, 68.03it/s]
INFO:linopy.io: Writing time: 0.2s
INFO:linopy.solvers:Log file at /tmp/highs.log

Running HiGHS 1.7.2 (git hash: 184e327): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e-04, 1e+00]
  Cost   [3e-02, 2e+05]
  Bound  [0e+00, 0e+00]
  RHS    [1e+00, 1e+00]
Presolving model
13239 rows, 10323 cols, 30957 nonzeros  0s
13239 rows, 10323 cols, 30957 nonzeros  0s
Presolve : Reductions: rows 13239(-13045); columns 10323(-1361); elements 30957(-14406)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     0.0000000000e+00 Pr: 2920(34348) 0s
       7363     9.7601841927e+05 Pr: 0(0) 0s
Solving the original LP from the solution after postsolve
Model   status      : Optimal
Simplex   iterations: 7363
Objective value     :  9.7601841927e+05
HiGHS run time      :          0.23

INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 11684 primals, 26284 duals
Objective: 9.76e+05
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.
/tmp/ipykernel_4603/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`

Writing the solution to /tmp/linopy-solve-ghb2q46z.sol

WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['DE'], dtype='object', name='Bus')
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['DE'], dtype='object', name='Bus')
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 5/5 [00:00<00:00, 33.53it/s]
Writing continuous variables.: 100%|██████████| 2/2 [00:00<00:00, 67.85it/s]
INFO:linopy.io: Writing time: 0.21s
INFO:linopy.solvers:Log file at /tmp/highs.log

Running HiGHS 1.7.2 (git hash: 184e327): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:

INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 11684 primals, 26284 duals
Objective: 1.12e+06
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.

  Matrix [1e-04, 1e+00]
  Cost   [3e-02, 2e+05]
  Bound  [0e+00, 0e+00]
  RHS    [1e+00, 1e+00]
Presolving model
13239 rows, 10323 cols, 30957 nonzeros  0s
13239 rows, 10323 cols, 30957 nonzeros  0s
Presolve : Reductions: rows 13239(-13045); columns 10323(-1361); elements 30957(-14406)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     0.0000000000e+00 Pr: 2920(34892) 0s
       7275     1.1182753980e+06 Pr: 0(0) 0s
Solving the original LP from the solution after postsolve
Model   status      : Optimal
Simplex   iterations: 7275
Objective value     :  1.1182753980e+06
HiGHS run time      :          0.29
Writing the solution to /tmp/linopy-solve-kewls50j.sol

[13]:

results

[13]:

Generator	DE solar	DE wind	DE gas	DE lignite
scenario
low	-0.000000	-0.000000	1.000000	-0.000000
med	1.481495	1.383751	0.990452	-0.000000
high	1.631694	1.924190	-0.000000	0.986723

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_4603/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`
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['DE'], dtype='object', name='Bus')
/tmp/ipykernel_4603/15070267.py:44: FutureWarning: DataFrame.groupby with axis=1 is deprecated. Do `frame.T.groupby(...)` without axis instead.
  .groupby(n.loads.bus, axis=1)
/tmp/ipykernel_4603/15070267.py:44: FutureWarning: DataFrame.groupby with axis=1 is deprecated. Do `frame.T.groupby(...)` without axis instead.
  .groupby(n.loads.bus, axis=1)
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 11/11 [00:00<00:00, 26.06it/s]
Writing continuous variables.: 100%|██████████| 4/4 [00:00<00:00, 48.03it/s]
INFO:linopy.io: Writing time: 0.57s
INFO:linopy.solvers:Log file at /tmp/highs.log

Running HiGHS 1.7.2 (git hash: 184e327): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e-04, 1e+00]
  Cost   [9e-03, 2e+05]
  Bound  [0e+00, 0e+00]
  RHS    [1e+00, 1e+00]
Presolving model
39717 rows, 30961 cols, 92871 nonzeros  0s
39717 rows, 30961 cols, 92871 nonzeros  0s
Presolve : Reductions: rows 39717(-39127); columns 30961(-4083); elements 92871(-43210)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     0.0000000000e+00 Pr: 8760(49664) 0s
      26078     9.5022247221e+05 Pr: 0(0) 4s
Solving the original LP from the solution after postsolve
Model   status      : Optimal
Simplex   iterations: 26078
Objective value     :  9.5022247221e+05
HiGHS run time      :          3.77
Writing the solution to /tmp/linopy-solve-w4vwt96s.sol

INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 35044 primals, 78844 duals
Objective: 9.50e+05
Solver model: available
Solver message: optimal

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

[14]:

('ok', 'optimal')

[15]:

results.loc["stochastic"] = n.generators.p_nom_opt

[16]:

results

[16]:

Generator	DE solar	DE wind	DE gas	DE lignite
scenario
low	-0.000000	-0.000000	1.000000	-0.000000
med	1.481495	1.383751	0.990452	-0.000000
high	1.631694	1.924190	-0.000000	0.986723
stochastic	1.436681	1.325307	0.901397	0.089458

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_4603/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`
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['DE'], dtype='object', name='Bus')
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['DE'], dtype='object', name='Bus')

Generator
DE solar     -0.0
DE wind      -0.0
DE gas        1.0
DE lignite   -0.0
Name: p_nom, dtype: float64

INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 3/3 [00:00<00:00, 26.17it/s]
Writing continuous variables.: 100%|██████████| 1/1 [00:00<00:00, 35.92it/s]
INFO:linopy.io: Writing time: 0.17s
INFO:linopy.solvers:Log file at /tmp/highs.log
INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 11680 primals, 26280 duals
Objective: 5.84e+05
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.
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['DE'], dtype='object', name='Bus')

Running HiGHS 1.7.2 (git hash: 184e327): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e+00, 1e+00]
  Cost   [3e-02, 4e+02]
  Bound  [0e+00, 0e+00]
  RHS    [1e+00, 1e+00]
Presolving model
0 rows, 0 cols, 0 nonzeros  0s
0 rows, 0 cols, 0 nonzeros  0s
Presolve : Reductions: rows 0(-26280); columns 0(-11680); elements 0(-35040) - Reduced to empty
Solving the original LP from the solution after postsolve
Model   status      : Optimal
Objective value     :  5.8400000000e+05
HiGHS run time      :          0.03
Writing the solution to /tmp/linopy-solve-qrqzo6qr.sol
Generator
DE solar      1.481495
DE wind       1.383751
DE gas        0.990452
DE lignite   -0.000000
Name: p_nom, dtype: float64

WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['DE'], dtype='object', name='Bus')
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 3/3 [00:00<00:00, 25.37it/s]
Writing continuous variables.: 100%|██████████| 1/1 [00:00<00:00, 35.05it/s]
INFO:linopy.io: Writing time: 0.18s
INFO:linopy.solvers:Log file at /tmp/highs.log
INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 11680 primals, 26280 duals
Objective: 3.11e+05
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.
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['DE'], dtype='object', name='Bus')

Running HiGHS 1.7.2 (git hash: 184e327): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e+00, 1e+00]
  Cost   [3e-02, 4e+02]
  Bound  [0e+00, 0e+00]
  RHS    [1e-04, 1e+00]
Presolving model
1549 rows, 4553 cols, 4547 nonzeros  0s
179 rows, 537 cols, 537 nonzeros  0s
Presolve : Reductions: rows 179(-26101); columns 537(-11143); elements 537(-34503)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     3.1083378625e+05 Pr: 179(179) 0s
        179     3.1084113430e+05 Pr: 0(0) 0s
Solving the original LP from the solution after postsolve
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
        179     3.1084104057e+05 Pr: 2(2) 0s
        183     3.1084113430e+05 Pr: 0(0) 0s
Model   status      : Optimal
Simplex   iterations: 183
Objective value     :  3.1084113430e+05
HiGHS run time      :          0.05
Writing the solution to /tmp/linopy-solve-rgk1ipmq.sol
Generator
DE solar      1.631694
DE wind       1.924190
DE gas       -0.000000
DE lignite    0.986723
Name: p_nom, dtype: float64

WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['DE'], dtype='object', name='Bus')
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 3/3 [00:00<00:00, 25.85it/s]
Writing continuous variables.: 100%|██████████| 1/1 [00:00<00:00, 36.19it/s]
INFO:linopy.io: Writing time: 0.17s
INFO:linopy.solvers:Log file at /tmp/highs.log
INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 11680 primals, 26280 duals
Objective: 5.79e+05
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.
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['DE'], dtype='object', name='Bus')

Running HiGHS 1.7.2 (git hash: 184e327): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e+00, 1e+00]
  Cost   [3e-02, 4e+02]
  Bound  [0e+00, 0e+00]
  RHS    [2e-04, 2e+00]
Presolving model
1540 rows, 4358 cols, 4344 nonzeros  0s
210 rows, 630 cols, 630 nonzeros  0s
Presolve : Reductions: rows 210(-26070); columns 630(-11050); elements 630(-34410)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     5.7876881154e+05 Pr: 210(210) 0s
        210     5.7877726371e+05 Pr: 0(0) 0s
Solving the original LP from the solution after postsolve
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
        210     5.7877717504e+05 Pr: 2(2) 0s
        214     5.7877726371e+05 Pr: 0(0) 0s
Model   status      : Optimal
Simplex   iterations: 214
Objective value     :  5.7877726371e+05
HiGHS run time      :          0.05
Writing the solution to /tmp/linopy-solve-ndi6ahyy.sol
Generator
DE solar      1.436681
DE wind       1.325307
DE gas        0.901397
DE lignite    0.089458
Name: p_nom, dtype: float64

WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['DE'], dtype='object', name='Bus')
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 3/3 [00:00<00:00, 25.36it/s]
Writing continuous variables.: 100%|██████████| 1/1 [00:00<00:00, 36.00it/s]
INFO:linopy.io: Writing time: 0.17s
INFO:linopy.solvers:Log file at /tmp/highs.log
INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 11680 primals, 26280 duals
Objective: 3.22e+05
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.

Running HiGHS 1.7.2 (git hash: 184e327): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e+00, 1e+00]
  Cost   [3e-02, 4e+02]
  Bound  [0e+00, 0e+00]
  RHS    [1e-04, 1e+00]
Presolving model
2877 rows, 10062 cols, 10059 nonzeros  0s
292 rows, 1130 cols, 1130 nonzeros  0s
Presolve : Reductions: rows 292(-25988); columns 1130(-10550); elements 1130(-33910)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     3.2038318115e+05 Pr: 292(291.978) 0s
        292     3.2158205576e+05 Pr: 0(0) 0s
Solving the original LP from the solution after postsolve
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
        292     3.2104119973e+05 Pr: 3(3) 0s
        301     3.2158205576e+05 Pr: 0(0) 0s
Model   status      : Optimal
Simplex   iterations: 301
Objective value     :  3.2158205576e+05
HiGHS run time      :          0.05
Writing the solution to /tmp/linopy-solve-qcsijual.sol

/tmp/ipykernel_4603/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`
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['DE'], dtype='object', name='Bus')
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['DE'], dtype='object', name='Bus')

Generator
DE solar     -0.0
DE wind      -0.0
DE gas        1.0
DE lignite   -0.0
Name: p_nom, dtype: float64

INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 3/3 [00:00<00:00, 25.93it/s]
Writing continuous variables.: 100%|██████████| 1/1 [00:00<00:00, 36.31it/s]
INFO:linopy.io: Writing time: 0.17s
INFO:linopy.solvers:Log file at /tmp/highs.log
INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 11680 primals, 26280 duals
Objective: 1.02e+06
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.
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['DE'], dtype='object', name='Bus')
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['DE'], dtype='object', name='Bus')

Running HiGHS 1.7.2 (git hash: 184e327): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e+00, 1e+00]
  Cost   [3e-02, 4e+02]
  Bound  [0e+00, 0e+00]
  RHS    [1e+00, 1e+00]
Presolving model
0 rows, 0 cols, 0 nonzeros  0s
0 rows, 0 cols, 0 nonzeros  0s
Presolve : Reductions: rows 0(-26280); columns 0(-11680); elements 0(-35040) - Reduced to empty
Solving the original LP from the solution after postsolve
Model   status      : Optimal
Objective value     :  1.0220000000e+06
HiGHS run time      :          0.03
Writing the solution to /tmp/linopy-solve-x6jmt5go.sol
Generator
DE solar      1.481495
DE wind       1.383751
DE gas        0.990452
DE lignite   -0.000000
Name: p_nom, dtype: float64

INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 3/3 [00:00<00:00, 25.65it/s]
Writing continuous variables.: 100%|██████████| 1/1 [00:00<00:00, 35.49it/s]
INFO:linopy.io: Writing time: 0.17s
INFO:linopy.solvers:Log file at /tmp/highs.log
INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 11680 primals, 26280 duals
Objective: 5.44e+05
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.
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['DE'], dtype='object', name='Bus')

Running HiGHS 1.7.2 (git hash: 184e327): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e+00, 1e+00]
  Cost   [3e-02, 4e+02]
  Bound  [0e+00, 0e+00]
  RHS    [1e-04, 1e+00]
Presolving model
1549 rows, 4553 cols, 4547 nonzeros  0s
179 rows, 537 cols, 537 nonzeros  0s
Presolve : Reductions: rows 179(-26101); columns 537(-11143); elements 537(-34503)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     5.4391505674e+05 Pr: 179(179) 0s
        179     5.4392240479e+05 Pr: 0(0) 0s
Solving the original LP from the solution after postsolve
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
        179     5.4392231105e+05 Pr: 2(2) 0s
        183     5.4392240479e+05 Pr: 0(0) 0s
Model   status      : Optimal
Simplex   iterations: 183
Objective value     :  5.4392240479e+05
HiGHS run time      :          0.06
Writing the solution to /tmp/linopy-solve-sz2k2hy1.sol
Generator
DE solar      1.631694
DE wind       1.924190
DE gas       -0.000000
DE lignite    0.986723
Name: p_nom, dtype: float64

WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['DE'], dtype='object', name='Bus')
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 3/3 [00:00<00:00, 25.23it/s]
Writing continuous variables.: 100%|██████████| 1/1 [00:00<00:00, 36.00it/s]
INFO:linopy.io: Writing time: 0.17s
INFO:linopy.solvers:Log file at /tmp/highs.log
INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 11680 primals, 26280 duals
Objective: 5.79e+05
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.
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['DE'], dtype='object', name='Bus')

Running HiGHS 1.7.2 (git hash: 184e327): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e+00, 1e+00]
  Cost   [3e-02, 4e+02]
  Bound  [0e+00, 0e+00]
  RHS    [2e-04, 2e+00]
Presolving model
1540 rows, 4358 cols, 4344 nonzeros  0s
210 rows, 630 cols, 630 nonzeros  0s
Presolve : Reductions: rows 210(-26070); columns 630(-11050); elements 630(-34410)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     5.7876881154e+05 Pr: 210(210) 0s
        210     5.7877726371e+05 Pr: 0(0) 0s
Solving the original LP from the solution after postsolve
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
        210     5.7877717504e+05 Pr: 2(2) 0s
        214     5.7877726371e+05 Pr: 0(0) 0s
Model   status      : Optimal
Simplex   iterations: 214
Objective value     :  5.7877726371e+05
HiGHS run time      :          0.05
Writing the solution to /tmp/linopy-solve-kummo058.sol
Generator
DE solar      1.436681
DE wind       1.325307
DE gas        0.901397
DE lignite    0.089458
Name: p_nom, dtype: float64

WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['DE'], dtype='object', name='Bus')
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 3/3 [00:00<00:00, 25.37it/s]
Writing continuous variables.: 100%|██████████| 1/1 [00:00<00:00, 36.15it/s]
INFO:linopy.io: Writing time: 0.17s
INFO:linopy.solvers:Log file at /tmp/highs.log
INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 11680 primals, 26280 duals
Objective: 5.59e+05
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.

Running HiGHS 1.7.2 (git hash: 184e327): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e+00, 1e+00]
  Cost   [3e-02, 4e+02]
  Bound  [0e+00, 0e+00]
  RHS    [1e-04, 1e+00]
Presolving model
2877 rows, 10062 cols, 10059 nonzeros  0s
292 rows, 1130 cols, 1130 nonzeros  0s
Presolve : Reductions: rows 292(-25988); columns 1130(-10550); elements 1130(-33910)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     5.5726757927e+05 Pr: 292(291.978) 0s
        292     5.5935626817e+05 Pr: 0(0) 0s
Solving the original LP from the solution after postsolve
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
        292     5.5840978342e+05 Pr: 3(3) 0s
        301     5.5935626817e+05 Pr: 0(0) 0s
Model   status      : Optimal
Simplex   iterations: 301
Objective value     :  5.5935626817e+05
HiGHS run time      :          0.05
Writing the solution to /tmp/linopy-solve-xxiitj3h.sol

/tmp/ipykernel_4603/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`
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['DE'], dtype='object', name='Bus')
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['DE'], dtype='object', name='Bus')

Generator
DE solar     -0.0
DE wind      -0.0
DE gas        1.0
DE lignite   -0.0
Name: p_nom, dtype: float64

INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 3/3 [00:00<00:00, 26.13it/s]
Writing continuous variables.: 100%|██████████| 1/1 [00:00<00:00, 35.98it/s]
INFO:linopy.io: Writing time: 0.17s
INFO:linopy.solvers:Log file at /tmp/highs.log
INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 11680 primals, 26280 duals
Objective: 1.46e+06
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.
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['DE'], dtype='object', name='Bus')
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['DE'], dtype='object', name='Bus')

Running HiGHS 1.7.2 (git hash: 184e327): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e+00, 1e+00]
  Cost   [3e-02, 5e+02]
  Bound  [0e+00, 0e+00]
  RHS    [1e+00, 1e+00]
Presolving model
0 rows, 0 cols, 0 nonzeros  0s
0 rows, 0 cols, 0 nonzeros  0s
Presolve : Reductions: rows 0(-26280); columns 0(-11680); elements 0(-35040) - Reduced to empty
Solving the original LP from the solution after postsolve
Model   status      : Optimal
Objective value     :  1.4600000000e+06
HiGHS run time      :          0.03
Writing the solution to /tmp/linopy-solve-evz6n31m.sol
Generator
DE solar      1.481495
DE wind       1.383751
DE gas        0.990452
DE lignite   -0.000000
Name: p_nom, dtype: float64

INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 3/3 [00:00<00:00, 25.37it/s]
Writing continuous variables.: 100%|██████████| 1/1 [00:00<00:00, 36.11it/s]
INFO:linopy.io: Writing time: 0.17s
INFO:linopy.solvers:Log file at /tmp/highs.log
INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 11680 primals, 26280 duals
Objective: 7.77e+05
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.
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['DE'], dtype='object', name='Bus')

Running HiGHS 1.7.2 (git hash: 184e327): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e+00, 1e+00]
  Cost   [3e-02, 5e+02]
  Bound  [0e+00, 0e+00]
  RHS    [1e-04, 1e+00]
Presolving model
1549 rows, 4553 cols, 4547 nonzeros  0s
179 rows, 537 cols, 537 nonzeros  0s
Presolve : Reductions: rows 179(-26101); columns 537(-11143); elements 537(-34503)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     7.7699632723e+05 Pr: 179(179) 0s
        179     7.7700367527e+05 Pr: 0(0) 0s
Solving the original LP from the solution after postsolve
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
        179     7.7700358154e+05 Pr: 2(2) 0s
        183     7.7700367527e+05 Pr: 0(0) 0s
Model   status      : Optimal
Simplex   iterations: 183
Objective value     :  7.7700367527e+05
HiGHS run time      :          0.05
Writing the solution to /tmp/linopy-solve-9ox0t593.sol
Generator
DE solar      1.631694
DE wind       1.924190
DE gas       -0.000000
DE lignite    0.986723
Name: p_nom, dtype: float64

WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['DE'], dtype='object', name='Bus')
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 3/3 [00:00<00:00, 25.72it/s]
Writing continuous variables.: 100%|██████████| 1/1 [00:00<00:00, 35.79it/s]
INFO:linopy.io: Writing time: 0.17s
INFO:linopy.solvers:Log file at /tmp/highs.log
INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 11680 primals, 26280 duals
Objective: 5.79e+05
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.
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['DE'], dtype='object', name='Bus')

Running HiGHS 1.7.2 (git hash: 184e327): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e+00, 1e+00]
  Cost   [3e-02, 5e+02]
  Bound  [0e+00, 0e+00]
  RHS    [2e-04, 2e+00]
Presolving model
1540 rows, 4358 cols, 4344 nonzeros  0s
210 rows, 630 cols, 630 nonzeros  0s
Presolve : Reductions: rows 210(-26070); columns 630(-11050); elements 630(-34410)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     5.7876881154e+05 Pr: 210(210) 0s
        210     5.7877726371e+05 Pr: 0(0) 0s
Solving the original LP from the solution after postsolve
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
        210     5.7877717504e+05 Pr: 2(2) 0s
        214     5.7877726371e+05 Pr: 0(0) 0s
Model   status      : Optimal
Simplex   iterations: 214
Objective value     :  5.7877726371e+05
HiGHS run time      :          0.05
Writing the solution to /tmp/linopy-solve-sa_xzg7x.sol
Generator
DE solar      1.436681
DE wind       1.325307
DE gas        0.901397
DE lignite    0.089458
Name: p_nom, dtype: float64

WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['DE'], dtype='object', name='Bus')
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 3/3 [00:00<00:00, 25.58it/s]
Writing continuous variables.: 100%|██████████| 1/1 [00:00<00:00, 36.26it/s]
INFO:linopy.io: Writing time: 0.17s
INFO:linopy.solvers:Log file at /tmp/highs.log
INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 11680 primals, 26280 duals
Objective: 7.86e+05
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.

Running HiGHS 1.7.2 (git hash: 184e327): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e+00, 1e+00]
  Cost   [3e-02, 5e+02]
  Bound  [0e+00, 0e+00]
  RHS    [1e-04, 1e+00]
Presolving model
2877 rows, 7529 cols, 7526 nonzeros  0s
292 rows, 1130 cols, 1130 nonzeros  0s
Presolve : Reductions: rows 292(-25988); columns 1130(-10550); elements 1130(-33910)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     7.8331602855e+05 Pr: 292(291.978) 0s
        292     7.8599792695e+05 Pr: 0(0) 0s
Solving the original LP from the solution after postsolve
Model   status      : Optimal
Simplex   iterations: 292
Objective value     :  7.8599792695e+05
HiGHS run time      :          0.04
Writing the solution to /tmp/linopy-solve-0sturp8t.sol

[18]:

results

[18]:

Generator	DE solar	DE wind	DE gas	DE lignite	gas-p-low	lignite-p-low	gas-p-med	lignite-p-med	gas-p-high	lignite-p-high
scenario
low	-0.000000	-0.000000	1.000000	-0.000000	2920.000000	0.000000	2920.000000	0.000000	2920.000000	0.000000
med	1.481495	1.383751	0.990452	-0.000000	1553.875137	0.000000	1553.875137	0.000000	1553.875137	0.000000
high	1.631694	1.924190	-0.000000	0.986723	0.000000	1285.992264	0.000000	1285.992264	0.000000	1285.992264
stochastic	1.436681	1.325307	0.901397	0.089458	1585.161416	9.968200	1585.161416	9.968200	1362.510344	232.619272

[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
        ]
    )

[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])

../_images/examples_stochastic-problem_29_0.png

[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])

[21]:

(0.0, 2000000.0)

../_images/examples_stochastic-problem_30_1.png

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_4603/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`
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['DE'], dtype='object', name='Bus')
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['DE'], dtype='object', name='Bus')
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 5/5 [00:00<00:00, 33.51it/s]
Writing continuous variables.: 100%|██████████| 2/2 [00:00<00:00, 69.33it/s]
INFO:linopy.io: Writing time: 0.21s
INFO:linopy.solvers:Log file at /tmp/highs.log

Running HiGHS 1.7.2 (git hash: 184e327): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e-04, 1e+00]
  Cost   [3e-02, 2e+05]
  Bound  [0e+00, 0e+00]
  RHS    [1e+00, 1e+00]
Presolving model
13239 rows, 10323 cols, 30957 nonzeros  0s
13239 rows, 10323 cols, 30957 nonzeros  0s
Presolve : Reductions: rows 13239(-13045); columns 10323(-1361); elements 30957(-14406)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     0.0000000000e+00 Pr: 2920(34348) 0s
       7362     9.5241319722e+05 Pr: 0(0) 0s
Solving the original LP from the solution after postsolve
Model   status      : Optimal
Simplex   iterations: 7362
Objective value     :  9.5241319722e+05
HiGHS run time      :          0.21

INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 11684 primals, 26284 duals
Objective: 9.52e+05
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.

Writing the solution to /tmp/linopy-solve-bamlimij.sol

[25]:

# pd.set_option("display.precision", 10)
portfolios
# costs

[25]:

	naive
Generator
DE solar	1.439943
DE wind	1.319087
DE gas	0.990898
DE lignite	-0.000000

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_4603/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`
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['DE'], dtype='object', name='Bus')
/tmp/ipykernel_4603/15070267.py:44: FutureWarning: DataFrame.groupby with axis=1 is deprecated. Do `frame.T.groupby(...)` without axis instead.
  .groupby(n.loads.bus, axis=1)
/tmp/ipykernel_4603/15070267.py:44: FutureWarning: DataFrame.groupby with axis=1 is deprecated. Do `frame.T.groupby(...)` without axis instead.
  .groupby(n.loads.bus, axis=1)
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 11/11 [00:00<00:00, 25.82it/s]
Writing continuous variables.: 100%|██████████| 4/4 [00:00<00:00, 47.99it/s]
INFO:linopy.io: Writing time: 0.58s
INFO:linopy.solvers:Log file at /tmp/highs.log

Running HiGHS 1.7.2 (git hash: 184e327): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e-04, 1e+00]
  Cost   [9e-03, 2e+05]
  Bound  [0e+00, 0e+00]
  RHS    [1e+00, 1e+00]
Presolving model
39717 rows, 30961 cols, 92871 nonzeros  0s
39717 rows, 30961 cols, 92871 nonzeros  0s
Presolve : Reductions: rows 39717(-39127); columns 30961(-4083); elements 92871(-43210)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     0.0000000000e+00 Pr: 8760(49664) 0s
      26078     9.5022247221e+05 Pr: 0(0) 4s
Solving the original LP from the solution after postsolve
Model   status      : Optimal
Simplex   iterations: 26078
Objective value     :  9.5022247221e+05
HiGHS run time      :          3.74

INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 35044 primals, 78844 duals
Objective: 9.50e+05
Solver model: available
Solver message: optimal

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

Writing the solution to /tmp/linopy-solve-wtvj6dao.sol

[27]:

portfolios

[27]:

	naive	SP
Generator
DE solar	1.439943	1.436681
DE wind	1.319087	1.325307
DE gas	0.990898	0.901397
DE lignite	-0.000000	0.089458

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_4603/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`
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['DE'], dtype='object', name='Bus')
/tmp/ipykernel_4603/15070267.py:44: FutureWarning: DataFrame.groupby with axis=1 is deprecated. Do `frame.T.groupby(...)` without axis instead.
  .groupby(n.loads.bus, axis=1)
/tmp/ipykernel_4603/15070267.py:44: FutureWarning: DataFrame.groupby with axis=1 is deprecated. Do `frame.T.groupby(...)` without axis instead.
  .groupby(n.loads.bus, axis=1)
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 11/11 [00:00<00:00, 25.98it/s]
Writing continuous variables.: 100%|██████████| 4/4 [00:00<00:00, 48.24it/s]
INFO:linopy.io: Writing time: 0.57s
INFO:linopy.solvers:Log file at /tmp/highs.log

Running HiGHS 1.7.2 (git hash: 184e327): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e-04, 1e+00]
  Cost   [9e-03, 2e+05]
  Bound  [0e+00, 0e+00]
  RHS    [1e+00, 1e+00]
Presolving model
39717 rows, 30961 cols, 92871 nonzeros  0s
39717 rows, 30961 cols, 92871 nonzeros  0s
Presolve : Reductions: rows 39717(-39127); columns 30961(-4083); elements 92871(-43210)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     0.0000000000e+00 Pr: 8760(49664) 0s
      26078     9.5022247221e+05 Pr: 0(0) 4s
Solving the original LP from the solution after postsolve
Model   status      : Optimal
Simplex   iterations: 26078
Objective value     :  9.5022247221e+05
HiGHS run time      :          3.75
Writing the solution to /tmp/linopy-solve-47dmnotx.sol

INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 35044 primals, 78844 duals
Objective: 9.50e+05
Solver model: available
Solver message: optimal

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

[28]:

('ok', 'optimal')

[29]:

# don't forget to add the capital costs of the (fixed) generators portfolio
c = "Generator"
ext_i = portfolios["naive"].index
cost = n.df(c)["capital_cost"][ext_i]
cost_of_portfolio = (n.generators.p_nom * cost).sum()
n.objective += cost_of_portfolio
n.objective

[29]:

1367406.1058110292

[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

[30]:

naive             9.524132e+05
SP                9.502225e+05
SP-constrained    1.367406e+06
dtype: float64

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"
)

ECIU: 0.417 Meuro
which is 44% of stochastic solution's costs

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_4603/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`
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['DE'], dtype='object', name='Bus')
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['DE'], dtype='object', name='Bus')
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 5/5 [00:00<00:00, 34.00it/s]
Writing continuous variables.: 100%|██████████| 2/2 [00:00<00:00, 69.78it/s]
INFO:linopy.io: Writing time: 0.2s
INFO:linopy.solvers:Log file at /tmp/highs.log

Running HiGHS 1.7.2 (git hash: 184e327): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e-04, 1e+00]
  Cost   [3e-02, 2e+05]
  Bound  [0e+00, 0e+00]
  RHS    [1e+00, 1e+00]
Presolving model
13239 rows, 10323 cols, 30957 nonzeros  0s
13239 rows, 10323 cols, 30957 nonzeros  0s
Presolve : Reductions: rows 13239(-13045); columns 10323(-1361); elements 30957(-14406)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     0.0000000000e+00 Pr: 2920(34348) 0s
       7405     6.4519950973e+05 Pr: 0(0) 0s
Solving the original LP from the solution after postsolve
Model   status      : Optimal
Simplex   iterations: 7405
Objective value     :  6.4519950973e+05
HiGHS run time      :          0.13
Writing the solution to /tmp/linopy-solve-t18xncqm.sol

INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 11684 primals, 26284 duals
Objective: 6.45e+05
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.
/tmp/ipykernel_4603/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`
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['DE'], dtype='object', name='Bus')
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['DE'], dtype='object', name='Bus')
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 5/5 [00:00<00:00, 33.38it/s]
Writing continuous variables.: 100%|██████████| 2/2 [00:00<00:00, 69.71it/s]
INFO:linopy.io: Writing time: 0.21s
INFO:linopy.solvers:Log file at /tmp/highs.log

Running HiGHS 1.7.2 (git hash: 184e327): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e-04, 1e+00]
  Cost   [3e-02, 2e+05]
  Bound  [0e+00, 0e+00]
  RHS    [1e+00, 1e+00]
Presolving model
13239 rows, 10323 cols, 30957 nonzeros  0s
13239 rows, 10323 cols, 30957 nonzeros  0s
Presolve : Reductions: rows 13239(-13045); columns 10323(-1361); elements 30957(-14406)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     0.0000000000e+00 Pr: 2920(34348) 0s
       7363     9.7601841927e+05 Pr: 0(0) 0s
Solving the original LP from the solution after postsolve
Model   status      : Optimal
Simplex   iterations: 7363
Objective value     :  9.7601841927e+05
HiGHS run time      :          0.23

INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 11684 primals, 26284 duals
Objective: 9.76e+05
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.

Writing the solution to /tmp/linopy-solve-u_h1r9an.sol

/tmp/ipykernel_4603/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`
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['DE'], dtype='object', name='Bus')
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['DE'], dtype='object', name='Bus')
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 5/5 [00:00<00:00, 33.60it/s]
Writing continuous variables.: 100%|██████████| 2/2 [00:00<00:00, 69.99it/s]
INFO:linopy.io: Writing time: 0.2s
INFO:linopy.solvers:Log file at /tmp/highs.log

Running HiGHS 1.7.2 (git hash: 184e327): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:

INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 11684 primals, 26284 duals
Objective: 1.12e+06
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.

  Matrix [1e-04, 1e+00]
  Cost   [3e-02, 2e+05]
  Bound  [0e+00, 0e+00]
  RHS    [1e+00, 1e+00]
Presolving model
13239 rows, 10323 cols, 30957 nonzeros  0s
13239 rows, 10323 cols, 30957 nonzeros  0s
Presolve : Reductions: rows 13239(-13045); columns 10323(-1361); elements 30957(-14406)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     0.0000000000e+00 Pr: 2920(34892) 0s
       7275     1.1182753980e+06 Pr: 0(0) 0s
Solving the original LP from the solution after postsolve
Model   status      : Optimal
Simplex   iterations: 7275
Objective value     :  1.1182753980e+06
HiGHS run time      :          0.30
Writing the solution to /tmp/linopy-solve-hu7bg_t6.sol

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_4603/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`
WARNING:pypsa.consistency:The following buses have carriers which are not defined:
Index(['DE'], dtype='object', name='Bus')
/tmp/ipykernel_4603/15070267.py:44: FutureWarning: DataFrame.groupby with axis=1 is deprecated. Do `frame.T.groupby(...)` without axis instead.
  .groupby(n.loads.bus, axis=1)
/tmp/ipykernel_4603/15070267.py:44: FutureWarning: DataFrame.groupby with axis=1 is deprecated. Do `frame.T.groupby(...)` without axis instead.
  .groupby(n.loads.bus, axis=1)
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.io:Writing objective.
Writing constraints.: 100%|██████████| 11/11 [00:00<00:00, 25.43it/s]
Writing continuous variables.: 100%|██████████| 4/4 [00:00<00:00, 48.13it/s]
INFO:linopy.io: Writing time: 0.58s
INFO:linopy.solvers:Log file at /tmp/highs.log

Running HiGHS 1.7.2 (git hash: 184e327): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e-04, 1e+00]
  Cost   [9e-03, 2e+05]
  Bound  [0e+00, 0e+00]
  RHS    [1e+00, 1e+00]
Presolving model
39717 rows, 30961 cols, 92871 nonzeros  0s
39717 rows, 30961 cols, 92871 nonzeros  0s
Presolve : Reductions: rows 39717(-39127); columns 30961(-4083); elements 92871(-43210)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0     0.0000000000e+00 Pr: 8760(49664) 0s
      26078     9.5022247221e+05 Pr: 0(0) 4s
Solving the original LP from the solution after postsolve
Model   status      : Optimal
Simplex   iterations: 26078
Objective value     :  9.5022247221e+05
HiGHS run time      :          3.79
Writing the solution to /tmp/linopy-solve-2sy9p95z.sol

INFO:linopy.constants: Optimization successful:
Status: ok
Termination condition: optimal
Solution: 35044 primals, 78844 duals
Objective: 9.50e+05
Solver model: available
Solver message: optimal

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

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"
)

EVPI: 0.064 Meuro
which is 7% of stochastic solution's costs

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.

Stochastic optimisation example

Contents

Stochastic optimisation example#

Required data#

Dependencies#

Retrieve PV & Wind data#

Major settings#

Prepare data#

Required functions#

First solve capacities for each scenario deterministically#

Now solve the full problem stochastically#

Now test each set of capacities against realisations of the gas price#

Analysis of a Stochastic Solution#

The Expected costs of ignoring uncertainty (ECIU)#

Define the naive problem (usually – the expected value problem (EV))#

solve naive problem (deterministic)#

solve stochastic problem#

Solve stochastic problem constrained by the naive solution#

Compute ECIU#

The Expected Value of Perfect Information (EVPI)#

Solve Wait-and-See problems#

compute the expected value of wait-and-see scenario costs#

solve stochastic problem#

Compute EVPI#

Comparing the ECIU and EVPI metrics#