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", "renewable_maxpower"]
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", "renewable_maxpower"]
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_2981/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.05s
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-d0mfq3zy.lp
INFO:gurobipy:Read LP format model from file /tmp/linopy-problem-d0mfq3zy.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: gurobi Runtime: 0.03s Dual bound: 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.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.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 [ OK ] consistency_check (0.000s) [START] direct [ OK ] direct (0.144s) [START] prepare_tssb_interface
[ OK ] prepare_tssb_interface (0.126s) [START] convert_to_blocks
[ OK ] convert_to_blocks (0.208s)
SMSNetwork
└── Block_0 [TwoStageStochasticBlock]
├── DiscreteScenarioSet [DiscreteScenarioSet]
├── StaticAbstractPath [AbstractPath]
└── StochasticBlock [StochasticBlock]
├── AbstractPath [AbstractPath]
└── Block [UCBlock]
├── UnitBlock_0 [IntermittentUnitBlock]
├── UnitBlock_1 [IntermittentUnitBlock]
├── UnitBlock_2 [SlackUnitBlock]
└── UnitBlock_3 [SlackUnitBlock]
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}")
[START] optimize
Executing command:
tssb_solver temp_stochastic_uc_example.nc -S TSSBSCfg.txt -c /home/docs/checkouts/readthedocs.org/user_builds/pypsa2smspp/conda/latest/lib/python3.14/site-packages/pysmspp/data/configs/TSSBlock/ -p /home/docs/checkouts/readthedocs.org/user_builds/pypsa2smspp/checkouts/latest/docs/examples/output/ -O /home/docs/checkouts/readthedocs.org/user_builds/pypsa2smspp/checkouts/latest/docs/examples/output/solution_stochastic_uc_example.nc
temp_stochastic_uc_example.nc is a Block file
Running HiGHS 1.14.0 (git hash: n/a): Copyright (c) 2026 under MIT licence terms
WARNING: LP matrix packed vector contains 183 |value| in [0, 1.22465e-16] less than or equal to 1e-09: ignored
Solver: HiGHSMILPSolver
LP has 724 rows; 582 cols; 1409 nonzeros
Coefficient ranges:
Matrix [2e-01, 1e+00]
Cost [3e-03, 3e+02]
Bound [7e+02, 1e+09]
RHS [1e+02, 4e+02]
WARNING: Problem has some excessively large column bounds
WARNING: Consider scaling the bounds by 1e-3, or setting the user_bound_scale option to -10
Presolving model
249 rows, 323 cols, 603 nonzeros 0s
138 rows, 173 cols, 342 nonzeros 0s
Dependent equations search running on 33 equations with time limit of 1000.00s
Dependent equations search removed 0 rows and 0 nonzeros in 0.00s (limit = 1000.00s)
138 rows, 140 cols, 309 nonzeros 0s
Presolve reductions: rows 138(-586); columns 140(-442); nonzeros 309(-1100)
Solving the presolved LP
Using dual simplex solver
Iteration Objective Infeasibilities num(sum)
0 0.0000000000e+00 Ph1: 0(0) 0.0s
71 1.6787625505e+05 Pr: 0(0) 0.0s
Performed postsolve
Solving the original LP from the solution after postsolve
Model status : Optimal
Simplex iterations: 71
Objective value : 1.6787625505e+05
P-D objective error : 8.6682133115e-17
HiGHS run time : 0.00
Elapsed time: 0.00439385 s
Status = 10 (Success)
Upper bound = 167876
Lower bound = 167876
Peak CPU Usage: 0.00 % Peak Memory Usage: 3.05 MB Total Time: 0.20 seconds [ OK ] optimize (0.250s) SMS++ objective: 167876.0
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
[START] parse_solution_to_unitblocks [split] No merged DCNetworkBlock entries found; nothing to do. [ OK ] parse_solution_to_unitblocks (0.000s) [START] inverse_transformation [ OK ] inverse_transformation (0.077s)
PyPSA objective: 167876.2550517915 SMS++ objective: 167876.0 Relative error SMS++ vs PyPSA: -0.000152%
Out[8]:
| 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.32015 | 0.00000 | 10844.32015 | 0.00000 | 0.216435 | 1.170264e+04 | 146137.71251 | 162.66480 | 0.0 | NaN |
| slack | 2.000000e+06 | 2000000.0 | 0.00000 | 0.00000 | 0.00000 | 0.00000 | 0.000000 | 4.800000e+07 | 0.00000 | 0.00000 | 0.0 | NaN | ||
| solar | 3.922823e+02 | 45.0 | 2979.67985 | 0.00000 | 2979.67985 | 0.00000 | 0.316490 | 0.000000e+00 | 21575.52618 | 29.79680 | 0.0 | NaN | ||
| low | onwind | 2.087682e+03 | 35.0 | 7706.22411 | 0.00000 | 7706.22411 | 0.00000 | 0.153803 | 8.076649e+03 | 146137.71251 | 115.59336 | 0.0 | NaN | |
| slack | 2.000000e+06 | 2000000.0 | 0.00000 | 0.00000 | 0.00000 | 0.00000 | 0.000000 | 4.800000e+07 | 0.00000 | 0.00000 | 0.0 | NaN | ||
| solar | 3.922823e+02 | 45.0 | 2085.77589 | 0.00000 | 2085.77589 | 0.00000 | 0.221543 | 0.000000e+00 | 21575.52618 | 20.85776 | 0.0 | NaN | ||
| med | onwind | 2.087682e+03 | 35.0 | 8987.27213 | 0.00000 | 8987.27213 | 0.00000 | 0.179371 | 1.017765e+04 | 146137.71251 | 134.80908 | 0.0 | NaN | |
| slack | 2.000000e+06 | 2000000.0 | 0.00000 | 0.00000 | 0.00000 | 0.00000 | 0.000000 | 4.800000e+07 | 0.00000 | 0.00000 | 0.0 | NaN | ||
| solar | 3.922823e+02 | 45.0 | 2532.72787 | 0.00000 | 2532.72787 | 0.00000 | 0.269016 | 0.000000e+00 | 21575.52618 | 25.32728 | 0.0 | NaN | ||
| Line | high | AC | 7.000000e+02 | 700.0 | 3541.52767 | 185.20752 | 3356.32015 | -3356.32015 | 0.221829 | 0.000000e+00 | 0.00000 | 0.00000 | 0.0 | NaN |
| low | AC | 7.000000e+02 | 700.0 | 2518.16638 | 115.94227 | 2402.22411 | -2402.22411 | 0.156792 | 0.000000e+00 | 0.00000 | 0.00000 | 0.0 | NaN | |
| med | AC | 7.000000e+02 | 700.0 | 2932.10450 | 184.83237 | 2747.27213 | -2747.27213 | 0.185532 | 0.000000e+00 | 0.00000 | 0.00000 | 0.0 | NaN | |
| Load | high | AC | 0.000000e+00 | 0.0 | 0.00000 | 13824.00000 | -13824.00000 | 0.00000 | NaN | 0.000000e+00 | 0.00000 | 0.00000 | 0.0 | NaN |
| low | AC | 0.000000e+00 | 0.0 | 0.00000 | 9792.00000 | -9792.00000 | 0.00000 | NaN | 0.000000e+00 | 0.00000 | 0.00000 | 0.0 | NaN | |
| med | AC | 0.000000e+00 | 0.0 | 0.00000 | 11520.00000 | -11520.00000 | 0.00000 | NaN | 0.000000e+00 | 0.00000 | 0.00000 | 0.0 | NaN |
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
Out[9]:
| case | objective_pypsa | objective_smspp | relative_error_pct | n_snapshots | n_scenarios | stochastic_parameters | |
|---|---|---|---|---|---|---|---|
| 0 | stochastic_uc_example | 167876.255052 | 167876.0 | -0.000152 | 24 | 3 | demand, renewable_maxpower |