Stochastic Optimization with pypsa2smspp¶
This notebook showcases an example on how to solve a PyPSA StochasticNetwork with SMS++.
To this goal, the following steps are performed:
- Create a small PyPSA network
- Optimizes it using PyPSA's built-in optimization capabilities to establish a baseline for comparison
- Convert it to SMS++ same instance through the
pypsa2smspptransformation pipeline - Optimizes the SMS++ model using the Two Stage Stochastic Block (ttsb_solver) solver from SMS++ tools
- Compare the results obtained from both optimizations
This notebook will use the Two Stage Stochastic Block (ttsb_solver) solver from SMS++ tools to optimize the model. To do so, pypsa2smspp will convert the PyPSA model to an TwoStageStochastic Block and optimize it using the default settings.
Note: it requires smspp-project version 0.6.0 or higher.
In [1]:
Copied!
from pathlib import Path
import numpy as np
import pandas as pd
import pypsa
from pypsa2smspp.transformation import Transformation
from pypsa2smspp.network_correction import add_slack_unit
from pathlib import Path
import numpy as np
import pandas as pd
import pypsa
from pypsa2smspp.transformation import Transformation
from pypsa2smspp.network_correction import add_slack_unit
Settings¶
The example uses three equiprobable scenarios. Both demand and renewable availability are scenario-dependent.
In [2]:
Copied!
CASE_NAME = "stochastic_uc_example"
SOLVER_NAME = "gurobi"
SOLVER_OPTIONS = {
"Threads": 32,
"Method": 2,
"Crossover": 0,
"Seed": 123,
"AggFill": 0,
"PreDual": 0,
}
SCENARIO_PROBABILITIES = {
"low": 1 / 3,
"med": 1 / 3,
"high": 1 / 3,
}
STOCHASTIC_PARAMETERS = ["demand", "renewables"]
CASE_NAME = "stochastic_uc_example"
SOLVER_NAME = "gurobi"
SOLVER_OPTIONS = {
"Threads": 32,
"Method": 2,
"Crossover": 0,
"Seed": 123,
"AggFill": 0,
"PreDual": 0,
}
SCENARIO_PROBABILITIES = {
"low": 1 / 3,
"med": 1 / 3,
"high": 1 / 3,
}
STOCHASTIC_PARAMETERS = ["demand", "renewables"]
Build the deterministic base network¶
The network has two AC buses, one transmission line, two loads, two renewable generators, and one committable gas generator.
The gas generator makes this a unit-commitment example. Demand profiles are generated directly in the notebook.
In [3]:
Copied!
snapshots = pd.date_range("2025-01-01 00:00", periods=24, freq="h")
n = pypsa.Network()
n.set_snapshots(snapshots)
# Keep snapshot weights explicit.
n.snapshot_weightings.loc[:, :] = 1.0
# Carriers
n.add("Carrier", "AC")
n.add("Carrier", "solar")
n.add("Carrier", "onwind")
n.add("Carrier", "slack")
# Buses
n.add("Bus", "IT0 0", carrier="AC", v_nom=380)
n.add("Bus", "IT0 1", carrier="AC", v_nom=380)
# Transmission line
n.add(
"Line",
"Line 0-1",
bus0="IT0 0",
bus1="IT0 1",
carrier="AC",
x=0.1,
r=0.01,
s_nom=700,
)
# Loads
n.add("Load", "load_IT0_0", bus="IT0 0", carrier="AC")
n.add("Load", "load_IT0_1", bus="IT0 1", carrier="AC")
# Renewable generators
n.add(
"Generator",
"solar_IT0_0",
bus="IT0 0",
carrier="solar",
p_nom=45,
p_min_pu=0,
p_max_pu=1,
marginal_cost=0.01,
p_nom_extendable=True,
capital_cost=55
)
n.add(
"Generator",
"wind_IT0_1",
bus="IT0 1",
carrier="onwind",
p_nom=35,
p_min_pu=0,
p_max_pu=1,
marginal_cost=0.015,
p_nom_extendable=True,
capital_cost=70
)
# Deterministic base demand.
hours = np.arange(len(snapshots))
load_0 = 220 + 60 * np.sin(2 * np.pi * (hours - 7) / 24)
load_1 = 260 + 70 * np.sin(2 * np.pi * (hours - 8) / 24)
load_0 = np.maximum(load_0, 150)
load_1 = np.maximum(load_1, 180)
n.loads_t.p_set = pd.DataFrame(
{
"load_IT0_0": load_0,
"load_IT0_1": load_1,
},
index=snapshots,
)
# Deterministic renewable availability.
solar_profile = np.maximum(0, np.sin(np.pi * (hours - 6) / 12))
wind_profile = 0.45 + 0.20 * np.sin(2 * np.pi * (hours + 3) / 24)
wind_profile = np.clip(wind_profile, 0.10, 0.85)
n.generators_t.p_max_pu = pd.DataFrame(
{
"solar_IT0_0": solar_profile,
"wind_IT0_1": wind_profile,
},
index=snapshots,
)
# Add high-cost slack units to guarantee feasibility in all scenarios.
for b in range(len(n.buses)):
n.add(
"Generator",
f"slack_{b}",
bus=f"IT0 {b}",
carrier="slack",
marginal_cost=1000,
p_nom=1e6,
)
n
snapshots = pd.date_range("2025-01-01 00:00", periods=24, freq="h")
n = pypsa.Network()
n.set_snapshots(snapshots)
# Keep snapshot weights explicit.
n.snapshot_weightings.loc[:, :] = 1.0
# Carriers
n.add("Carrier", "AC")
n.add("Carrier", "solar")
n.add("Carrier", "onwind")
n.add("Carrier", "slack")
# Buses
n.add("Bus", "IT0 0", carrier="AC", v_nom=380)
n.add("Bus", "IT0 1", carrier="AC", v_nom=380)
# Transmission line
n.add(
"Line",
"Line 0-1",
bus0="IT0 0",
bus1="IT0 1",
carrier="AC",
x=0.1,
r=0.01,
s_nom=700,
)
# Loads
n.add("Load", "load_IT0_0", bus="IT0 0", carrier="AC")
n.add("Load", "load_IT0_1", bus="IT0 1", carrier="AC")
# Renewable generators
n.add(
"Generator",
"solar_IT0_0",
bus="IT0 0",
carrier="solar",
p_nom=45,
p_min_pu=0,
p_max_pu=1,
marginal_cost=0.01,
p_nom_extendable=True,
capital_cost=55
)
n.add(
"Generator",
"wind_IT0_1",
bus="IT0 1",
carrier="onwind",
p_nom=35,
p_min_pu=0,
p_max_pu=1,
marginal_cost=0.015,
p_nom_extendable=True,
capital_cost=70
)
# Deterministic base demand.
hours = np.arange(len(snapshots))
load_0 = 220 + 60 * np.sin(2 * np.pi * (hours - 7) / 24)
load_1 = 260 + 70 * np.sin(2 * np.pi * (hours - 8) / 24)
load_0 = np.maximum(load_0, 150)
load_1 = np.maximum(load_1, 180)
n.loads_t.p_set = pd.DataFrame(
{
"load_IT0_0": load_0,
"load_IT0_1": load_1,
},
index=snapshots,
)
# Deterministic renewable availability.
solar_profile = np.maximum(0, np.sin(np.pi * (hours - 6) / 12))
wind_profile = 0.45 + 0.20 * np.sin(2 * np.pi * (hours + 3) / 24)
wind_profile = np.clip(wind_profile, 0.10, 0.85)
n.generators_t.p_max_pu = pd.DataFrame(
{
"solar_IT0_0": solar_profile,
"wind_IT0_1": wind_profile,
},
index=snapshots,
)
# Add high-cost slack units to guarantee feasibility in all scenarios.
for b in range(len(n.buses)):
n.add(
"Generator",
f"slack_{b}",
bus=f"IT0 {b}",
carrier="slack",
marginal_cost=1000,
p_nom=1e6,
)
n
Out[3]:
PyPSA Network 'Unnamed Network' ------------------------------- Components: - Bus: 2 - Carrier: 4 - Generator: 4 - Line: 1 - Load: 2 Snapshots: 24
Convert the network to a stochastic PyPSA network¶
The deterministic profiles are copied before calling set_scenarios. Then each scenario receives its own demand and renewable availability.
In [4]:
Copied!
base_load = n.loads_t.p_set.copy()
base_p_max_pu = n.generators_t.p_max_pu.copy()
n.set_scenarios(SCENARIO_PROBABILITIES)
# Scenario-dependent demand.
demand_multiplier = {
"low": 0.85,
"med": 1.00,
"high": 1.20,
}
# Scenario-dependent renewable availability.
renewable_multiplier = {
"low": 0.70,
"med": 0.85,
"high": 1.00,
}
for scenario in SCENARIO_PROBABILITIES:
n.loads_t.p_set[scenario] = base_load * demand_multiplier[scenario]
n.generators_t.p_max_pu[scenario] = (base_p_max_pu * renewable_multiplier[scenario]).clip(upper=1.0)
n
base_load = n.loads_t.p_set.copy()
base_p_max_pu = n.generators_t.p_max_pu.copy()
n.set_scenarios(SCENARIO_PROBABILITIES)
# Scenario-dependent demand.
demand_multiplier = {
"low": 0.85,
"med": 1.00,
"high": 1.20,
}
# Scenario-dependent renewable availability.
renewable_multiplier = {
"low": 0.70,
"med": 0.85,
"high": 1.00,
}
for scenario in SCENARIO_PROBABILITIES:
n.loads_t.p_set[scenario] = base_load * demand_multiplier[scenario]
n.generators_t.p_max_pu[scenario] = (base_p_max_pu * renewable_multiplier[scenario]).clip(upper=1.0)
n
Out[4]:
Stochastic PyPSA Network 'Unnamed Network' ------------------------------------------ Components: - Bus: 6 - Carrier: 12 - Generator: 12 - Line: 3 - Load: 6 Snapshots: 24 Scenarios: 3
Solve the stochastic unit-commitment problem with PyPSA¶
This is the PyPSA reference solution. The solve uses normal unit commitment, not linearized unit commitment.
In [5]:
Copied!
n_pypsa = n.copy()
n_pypsa.optimize(
solver_name=SOLVER_NAME,
solver_options=SOLVER_OPTIONS,
)
obj_pypsa = n_pypsa.objective + n_pypsa.objective_constant
statistics_pypsa = n_pypsa.statistics()
print(f"PyPSA objective: {obj_pypsa}")
statistics_pypsa
n_pypsa = n.copy()
n_pypsa.optimize(
solver_name=SOLVER_NAME,
solver_options=SOLVER_OPTIONS,
)
obj_pypsa = n_pypsa.objective + n_pypsa.objective_constant
statistics_pypsa = n_pypsa.statistics()
print(f"PyPSA objective: {obj_pypsa}")
statistics_pypsa
/tmp/ipykernel_1940/3437769620.py:3: FutureWarning: The default value of `include_objective_constant` will change from True to False in version 2.0. Set `include_objective_constant` explicitly to suppress this warning. Using False improves LP numerical conditioning by not including the objective constant as a variable.
n_pypsa.optimize(
WARNING:linopy.expressions:Constant RHS contains dimensions {'scenario'} not present in the expression, which might lead to inefficiencies. Consider collapsing the dimensions by taking min/max.
INFO:linopy.model: Solve problem using Gurobi solver
INFO:linopy.model:Solver options: - Threads: 32 - Method: 2 - Crossover: 0 - Seed: 123 - AggFill: 0 - PreDual: 0
INFO:linopy.io: Writing time: 0.04s
Restricted license - for non-production use only - expires 2027-11-29
INFO:gurobipy:Restricted license - for non-production use only - expires 2027-11-29
Read LP format model from file /tmp/linopy-problem-mdwnmv70.lp
INFO:gurobipy:Read LP format model from file /tmp/linopy-problem-mdwnmv70.lp
Reading time = 0.00 seconds
INFO:gurobipy:Reading time = 0.00 seconds
obj: 870 rows, 365 columns, 1263 nonzeros
INFO:gurobipy:obj: 870 rows, 365 columns, 1263 nonzeros
Set parameter Threads to value 32
INFO:gurobipy:Set parameter Threads to value 32
Set parameter Method to value 2
INFO:gurobipy:Set parameter Method to value 2
Set parameter Crossover to value 0
INFO:gurobipy:Set parameter Crossover to value 0
Set parameter Seed to value 123
INFO:gurobipy:Set parameter Seed to value 123
Set parameter AggFill to value 0
INFO:gurobipy:Set parameter AggFill to value 0
Set parameter PreDual to value 0
INFO:gurobipy:Set parameter PreDual to value 0
Gurobi Optimizer version 13.0.2 build v13.0.2rc1 (linux64 - "Ubuntu 24.04 LTS")
INFO:gurobipy:Gurobi Optimizer version 13.0.2 build v13.0.2rc1 (linux64 - "Ubuntu 24.04 LTS")
INFO:gurobipy:
CPU model: AMD EPYC 7R13 Processor, instruction set [SSE2|AVX|AVX2]
INFO:gurobipy:CPU model: AMD EPYC 7R13 Processor, instruction set [SSE2|AVX|AVX2]
Thread count: 1 physical cores, 2 logical processors, using up to 32 threads
INFO:gurobipy:Thread count: 1 physical cores, 2 logical processors, using up to 32 threads
INFO:gurobipy:
Warning: Thread count (32) is larger than processor count (2)
INFO:gurobipy:Warning: Thread count (32) is larger than processor count (2)
Reduce the value of the Threads parameter to improve performance
INFO:gurobipy: Reduce the value of the Threads parameter to improve performance
INFO:gurobipy:
INFO:gurobipy:
Non-default parameters:
INFO:gurobipy:Non-default parameters:
Method 2
INFO:gurobipy:Method 2
Crossover 0
INFO:gurobipy:Crossover 0
AggFill 0
INFO:gurobipy:AggFill 0
PreDual 0
INFO:gurobipy:PreDual 0
Seed 123
INFO:gurobipy:Seed 123
Threads 32
INFO:gurobipy:Threads 32
INFO:gurobipy:
Optimize a model with 870 rows, 365 columns and 1263 nonzeros (Min)
INFO:gurobipy:Optimize a model with 870 rows, 365 columns and 1263 nonzeros (Min)
Model fingerprint: 0xf9a8ddb8
INFO:gurobipy:Model fingerprint: 0xf9a8ddb8
Model has 293 linear objective coefficients
INFO:gurobipy:Model has 293 linear objective coefficients
Coefficient statistics:
INFO:gurobipy:Coefficient statistics:
Matrix range [2e-01, 1e+00]
INFO:gurobipy: Matrix range [2e-01, 1e+00]
Objective range [3e-03, 3e+02]
INFO:gurobipy: Objective range [3e-03, 3e+02]
Bounds range [5e+03, 5e+03]
INFO:gurobipy: Bounds range [5e+03, 5e+03]
RHS range [1e+02, 1e+06]
INFO:gurobipy: RHS range [1e+02, 1e+06]
INFO:gurobipy:
Presolve removed 814 rows and 325 columns
INFO:gurobipy:Presolve removed 814 rows and 325 columns
Presolve time: 0.01s
INFO:gurobipy:Presolve time: 0.01s
Presolved: 56 rows, 40 columns, 112 nonzeros
INFO:gurobipy:Presolved: 56 rows, 40 columns, 112 nonzeros
Ordering time: 0.00s
INFO:gurobipy:Ordering time: 0.00s
INFO:gurobipy:
Barrier statistics:
INFO:gurobipy:Barrier statistics:
AA' NZ : 7.040e+02
INFO:gurobipy: AA' NZ : 7.040e+02
Factor NZ : 1.244e+03
INFO:gurobipy: Factor NZ : 1.244e+03
Factor Ops : 3.372e+04 (less than 1 second per iteration)
INFO:gurobipy: Factor Ops : 3.372e+04 (less than 1 second per iteration)
Threads : 1
INFO:gurobipy: Threads : 1
INFO:gurobipy:
Objective Residual
INFO:gurobipy: Objective Residual
Iter Primal Dual Primal Dual Compl Time
INFO:gurobipy:Iter Primal Dual Primal Dual Compl Time
0 2.53752135e+05 1.07113175e+05 7.33e+02 0.00e+00 7.31e+03 0s
INFO:gurobipy: 0 2.53752135e+05 1.07113175e+05 7.33e+02 0.00e+00 7.31e+03 0s
1 1.75850787e+05 1.16578764e+05 1.72e+01 2.09e-15 6.03e+02 0s
INFO:gurobipy: 1 1.75850787e+05 1.16578764e+05 1.72e+01 2.09e-15 6.03e+02 0s
2 1.63451597e+05 1.49865716e+05 5.71e+00 5.64e-15 1.33e+02 0s
INFO:gurobipy: 2 1.63451597e+05 1.49865716e+05 5.71e+00 5.64e-15 1.33e+02 0s
3 1.65031279e+05 1.60833816e+05 7.99e-03 3.20e-15 3.53e+01 0s
INFO:gurobipy: 3 1.65031279e+05 1.60833816e+05 7.99e-03 3.20e-15 3.53e+01 0s
4 1.63024582e+05 1.62774026e+05 1.05e-05 5.68e-14 2.11e+00 0s
INFO:gurobipy: 4 1.63024582e+05 1.62774026e+05 1.05e-05 5.68e-14 2.11e+00 0s
5 1.62957085e+05 1.62946627e+05 1.59e-08 6.10e-15 8.79e-02 0s
INFO:gurobipy: 5 1.62957085e+05 1.62946627e+05 1.59e-08 6.10e-15 8.79e-02 0s
6 1.62951411e+05 1.62951078e+05 3.52e-11 1.14e-13 2.80e-03 0s
INFO:gurobipy: 6 1.62951411e+05 1.62951078e+05 3.52e-11 1.14e-13 2.80e-03 0s
7 1.62951255e+05 1.62951255e+05 2.76e-11 7.11e-15 3.78e-06 0s
INFO:gurobipy: 7 1.62951255e+05 1.62951255e+05 2.76e-11 7.11e-15 3.78e-06 0s
8 1.62951255e+05 1.62951255e+05 5.65e-11 1.63e-14 3.78e-09 0s
INFO:gurobipy: 8 1.62951255e+05 1.62951255e+05 5.65e-11 1.63e-14 3.78e-09 0s
INFO:gurobipy:
Barrier solved model in 8 iterations and 0.03 seconds (0.00 work units)
INFO:gurobipy:Barrier solved model in 8 iterations and 0.03 seconds (0.00 work units)
Optimal objective 1.62951255e+05
INFO:gurobipy:Optimal objective 1.62951255e+05
INFO:gurobipy:
INFO:linopy.constants: Optimization successful: Status: ok Termination condition: optimal Solution: 365 primals, 870 duals Objective: 1.63e+05 Solver model: available Solver message: 2
INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Generator-ext-p-lower, Generator-ext-p-upper, Line-fix-s-lower, Line-fix-s-upper were not assigned to the network.
PyPSA objective: 167876.2550517915
Out[5]:
| Optimal Capacity | Installed Capacity | Supply | Withdrawal | Energy Balance | Transmission | Capacity Factor | Curtailment | Capital Expenditure | Operational Expenditure | Revenue | Market Value | |||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Generator | high | onwind | 2.087682e+03 | 35.0 | 10844.32020 | 0.00000 | 10844.32020 | 0.00000 | 0.216435 | 1.170264e+04 | 146137.71251 | 162.66480 | 54.22160 | 0.005000 |
| slack | 2.000000e+06 | 2000000.0 | 0.00000 | 0.00000 | 0.00000 | 0.00000 | 0.000000 | 4.800000e+07 | 0.00000 | 0.00000 | 0.00000 | NaN | ||
| solar | 3.922823e+02 | 45.0 | 2979.67980 | 0.00000 | 2979.67980 | 0.00000 | 0.316490 | 5.000000e-05 | 21575.52618 | 29.79680 | 14.89840 | 0.005000 | ||
| low | onwind | 2.087682e+03 | 35.0 | 7706.22414 | 0.00000 | 7706.22414 | 0.00000 | 0.153803 | 8.076649e+03 | 146137.71251 | 115.59336 | 146176.24363 | 18.968595 | |
| slack | 2.000000e+06 | 2000000.0 | 0.00000 | 0.00000 | 0.00000 | 0.00000 | 0.000000 | 4.800000e+07 | 0.00000 | 0.00000 | 0.00000 | NaN | ||
| solar | 3.922823e+02 | 45.0 | 2085.77586 | 0.00000 | 2085.77586 | 0.00000 | 0.221543 | 4.000000e-05 | 21575.52618 | 20.85776 | 21573.29142 | 10.343053 | ||
| med | onwind | 2.087682e+03 | 35.0 | 8987.27216 | 0.00000 | 8987.27216 | 0.00000 | 0.179371 | 1.017764e+04 | 146137.71251 | 134.80908 | 44.93636 | 0.005000 | |
| slack | 2.000000e+06 | 2000000.0 | 0.00000 | 0.00000 | 0.00000 | 0.00000 | 0.000000 | 4.800000e+07 | 0.00000 | 0.00000 | 0.00000 | NaN | ||
| solar | 3.922823e+02 | 45.0 | 2532.72784 | 0.00000 | 2532.72784 | 0.00000 | 0.269016 | 4.000000e-05 | 21575.52618 | 25.32728 | 12.66364 | 0.005000 | ||
| Line | high | AC | 7.000000e+02 | 700.0 | 3541.52771 | 185.20751 | 3356.32020 | -3356.32020 | 0.221829 | 0.000000e+00 | 0.00000 | 0.00000 | 16.78160 | 0.004503 |
| low | AC | 7.000000e+02 | 700.0 | 2518.16640 | 115.94226 | 2402.22414 | -2402.22414 | 0.156792 | 0.000000e+00 | 0.00000 | 0.00000 | 53315.33461 | 20.240370 | |
| med | AC | 7.000000e+02 | 700.0 | 2932.10452 | 184.83236 | 2747.27216 | -2747.27216 | 0.185532 | 0.000000e+00 | 0.00000 | 0.00000 | 13.73636 | 0.004407 | |
| Load | high | AC | 0.000000e+00 | 0.0 | 0.00000 | 13824.00000 | -13824.00000 | 0.00000 | NaN | 0.000000e+00 | 0.00000 | 0.00000 | -69.12000 | -0.005000 |
| low | AC | 0.000000e+00 | 0.0 | 0.00000 | 9792.00000 | -9792.00000 | 0.00000 | NaN | 0.000000e+00 | 0.00000 | 0.00000 | -167749.53505 | -17.131284 | |
| med | AC | 0.000000e+00 | 0.0 | 0.00000 | 11520.00000 | -11520.00000 | 0.00000 | NaN | 0.000000e+00 | 0.00000 | 0.00000 | -57.60000 | -0.005000 |
Build and solve the SMS++ model¶
The transformation is configured for a Two-Stage Stochastic Block (tssb) with stochastic demand and renewable availability.
In [6]:
Copied!
n_smspp = n.copy()
transformation = Transformation(
name=CASE_NAME,
configfile="TSSBlock/TSSBSCfg_grb.txt",
capacity_expansion_ucblock=True,
stochastic_parameters={
"stochastic_type": "tssb",
"parameters": STOCHASTIC_PARAMETERS,
},
)
smspp_network = transformation.create_model(n_smspp, verbose=True)
smspp_network.print_tree()
n_smspp = n.copy()
transformation = Transformation(
name=CASE_NAME,
configfile="TSSBlock/TSSBSCfg_grb.txt",
capacity_expansion_ucblock=True,
stochastic_parameters={
"stochastic_type": "tssb",
"parameters": STOCHASTIC_PARAMETERS,
},
)
smspp_network = transformation.create_model(n_smspp, verbose=True)
smspp_network.print_tree()
[START] consistency_check [FAIL ] consistency_check (0.000s) -> Unsupported stochastic parameters: ['renewables']. Supported values are ['demand', 'hydro_inflow', 'renewable_marginal_cost', 'renewable_maxpower', 'thermal_marginal_cost', 'thermal_marginal_cost_quadratic', 'thermal_stand_by_cost'].
--------------------------------------------------------------------------- ValueError Traceback (most recent call last) Cell In[6], line 13 9 "parameters": STOCHASTIC_PARAMETERS, 10 }, 11 ) 12 ---> 13 smspp_network = transformation.create_model(n_smspp, verbose=True) 14 15 smspp_network.print_tree() File ~/checkouts/readthedocs.org/user_builds/pypsa2smspp/checkouts/latest/pypsa2smspp/transformation.py:299, in Transformation.create_model(self, n, verbose) 296 n_direct = get_base_scenario_network(n) 298 with step(self.timer, "consistency_check", verbose=verbose): --> 299 self.consistency_check(n) 301 with step(self.timer, "direct", verbose=verbose): 302 self.direct(n_direct) File ~/checkouts/readthedocs.org/user_builds/pypsa2smspp/checkouts/latest/pypsa2smspp/transformation.py:2036, in Transformation.consistency_check(self, n) 2033 raise TypeError("capacity_expansion_ucblock must be a boolean.") 2035 # ---- Describe high-level problem structure ---- -> 2036 self.problem_structure = describe_problem_structure( 2037 n, 2038 capacity_expansion_ucblock=self.capacity_expansion_ucblock, 2039 stochastic_parameters=self.stochastic_parameters, 2040 ) 2042 # ---- Minimal stochastic consistency checks ---- 2043 if self.problem_structure["is_stochastic"]: File ~/checkouts/readthedocs.org/user_builds/pypsa2smspp/checkouts/latest/pypsa2smspp/stochastic_utils.py:253, in describe_problem_structure(n, capacity_expansion_ucblock, stochastic_parameters) 250 is_stochastic = is_stochastic_network(n) 251 scenario_names = get_scenario_names(n) --> 253 sp = _normalize_stochastic_parameters(stochastic_parameters) 255 stochastic_type = sp["stochastic_type"] if is_stochastic else None 256 stochastic_parameters_list = list(sp["parameters"]) if is_stochastic else [] File ~/checkouts/readthedocs.org/user_builds/pypsa2smspp/checkouts/latest/pypsa2smspp/stochastic_utils.py:232, in _normalize_stochastic_parameters(stochastic_parameters) 229 invalid = sorted(set(parameters) - valid_parameters) 231 if invalid: --> 232 raise ValueError( 233 f"Unsupported stochastic parameters: {invalid}. " 234 f"Supported values are {sorted(valid_parameters)}." 235 ) 237 return { 238 "stochastic_type": stochastic_type, 239 "parameters": parameters, 240 } ValueError: Unsupported stochastic parameters: ['renewables']. Supported values are ['demand', 'hydro_inflow', 'renewable_marginal_cost', 'renewable_maxpower', 'thermal_marginal_cost', 'thermal_marginal_cost_quadratic', 'thermal_stand_by_cost'].
In [7]:
Copied!
transformation.optimize(verbose=True)
print(f"SMS++ objective: {transformation.result.objective_value}")
transformation.optimize(verbose=True)
print(f"SMS++ objective: {transformation.result.objective_value}")
--------------------------------------------------------------------------- ValueError Traceback (most recent call last) Cell In[7], line 1 ----> 1 transformation.optimize(verbose=True) 2 3 print(f"SMS++ objective: {transformation.result.objective_value}") File ~/checkouts/readthedocs.org/user_builds/pypsa2smspp/checkouts/latest/pypsa2smspp/transformation.py:314, in Transformation.optimize(self, verbose) 312 def optimize(self, verbose: bool = True): 313 if self.sms_network is None: --> 314 raise ValueError("Model must be created before optimization. Call create_model(n) first.") 315 with step(self.timer, "optimize", verbose=verbose, extra={"mode": "auto"}): 316 return self._optimize() ValueError: Model must be created before optimization. Call create_model(n) first.
Retrieve the SMS++ solution back into PyPSA¶
In [8]:
Copied!
n_smspp = transformation.retrieve_solution(n_smspp, verbose=True)
obj_smspp = transformation.result.objective_value
statistics_smspp = n_smspp.statistics()
relative_error_pct = (obj_smspp - obj_pypsa) / obj_pypsa * 100
print(f"PyPSA objective: {obj_pypsa}")
print(f"SMS++ objective: {obj_smspp}")
print(f"Relative error SMS++ vs PyPSA: {relative_error_pct:.6f}%")
statistics_smspp
n_smspp = transformation.retrieve_solution(n_smspp, verbose=True)
obj_smspp = transformation.result.objective_value
statistics_smspp = n_smspp.statistics()
relative_error_pct = (obj_smspp - obj_pypsa) / obj_pypsa * 100
print(f"PyPSA objective: {obj_pypsa}")
print(f"SMS++ objective: {obj_smspp}")
print(f"Relative error SMS++ vs PyPSA: {relative_error_pct:.6f}%")
statistics_smspp
--------------------------------------------------------------------------- ValueError Traceback (most recent call last) Cell In[8], line 1 ----> 1 n_smspp = transformation.retrieve_solution(n_smspp, verbose=True) 2 3 obj_smspp = transformation.result.objective_value 4 statistics_smspp = n_smspp.statistics() File ~/checkouts/readthedocs.org/user_builds/pypsa2smspp/checkouts/latest/pypsa2smspp/transformation.py:320, in Transformation.retrieve_solution(self, n, verbose) 318 def retrieve_solution(self, n, verbose: bool = True): 319 if self.result is None: --> 320 raise ValueError("Optimization must be run before retrieving solution. Call optimize() first.") 321 with step(self.timer, "parse_solution_to_unitblocks", verbose=verbose): 322 self.parse_solution_to_unitblocks(self.result.solution, n) ValueError: Optimization must be run before retrieving solution. Call optimize() first.
Print the summary¶
In [9]:
Copied!
summary = pd.DataFrame(
[
{
"case": CASE_NAME,
"objective_pypsa": obj_pypsa,
"objective_smspp": obj_smspp,
"relative_error_pct": relative_error_pct,
"n_snapshots": len(n.snapshots),
"n_scenarios": len(SCENARIO_PROBABILITIES),
"stochastic_parameters": ", ".join(STOCHASTIC_PARAMETERS),
}
]
)
summary
summary = pd.DataFrame(
[
{
"case": CASE_NAME,
"objective_pypsa": obj_pypsa,
"objective_smspp": obj_smspp,
"relative_error_pct": relative_error_pct,
"n_snapshots": len(n.snapshots),
"n_scenarios": len(SCENARIO_PROBABILITIES),
"stochastic_parameters": ", ".join(STOCHASTIC_PARAMETERS),
}
]
)
summary
--------------------------------------------------------------------------- NameError Traceback (most recent call last) Cell In[9], line 6 2 [ 3 { 4 "case": CASE_NAME, 5 "objective_pypsa": obj_pypsa, ----> 6 "objective_smspp": obj_smspp, 7 "relative_error_pct": relative_error_pct, 8 "n_snapshots": len(n.snapshots), 9 "n_scenarios": len(SCENARIO_PROBABILITIES), NameError: name 'obj_smspp' is not defined