Capacity expansion with InvestmentBlock¶

This notebook describes an example on how to use InvestmentBlock solver from SMS++ tools that allows using Benders Decomposition for capacity expansion, where investment decisions are located in the master problem and dispatch ones in the lower problem. To do so, the PyPSA model is converted to an InvestmentBlock and then solved using the appropriate solver.

This notebook relies on two libraries:

• pySMSpp: it aims to provide an abstract python interface for SMS++ models
• pypsa2smspp: it aims to provide a converter from PyPSA to SMS++ models exploiting pySMSpp

In this notebook, the following steps are performed:

1. Create a simple PyPSA model
2. Define the configuration to use InvestmentBlock
3. Convert the PyPSA model to SMS++ using pypsa2smspp and optimize the model
4. Verify the equivalence of objective function
5. Analyze the PyPSA model created with SMS++ and pypsa2smspp and verify optimal investment decisions

1. Creation of the PyPSA model¶

We create a simple PyPSA model arbitrary number of buses, few generators, storages and loads.

Main assumptions are:

• The network is purely radial in the form bus1 -> bus2 -> bus3 -> ... busN
• A load is added to each bus of the network
• A generator is added to the first node

In [1]:

import pypsa
import pandas as pd
import pypsa2smspp

import pypsa import pandas as pd import pypsa2smspp

The following code creates the desired pypsa model.¶

In [2]:

n = pypsa.Network()
n.set_snapshots(pd.date_range("2024-01-01T00:00", "2024-01-01T23:00", freq="h"))

# Add carriers
n.add("Carrier", "AC")

# Add buses
n_buses = 2
for b in range(n_buses):
    n.add("Bus", f"bus{b}", carrier="AC")

# Add lines in a radial topology using bidirectional links
n_lines = n_buses - 1
for l in range(n_lines):
    n.add(
        "Link",
        f"line{l}",
        bus0=f"bus{l}",
        bus1=f"bus{l+1}",
        length=1,
        capital_cost=1000,
        p_min_pu=-1,
        p_nom_extendable=True,
    )

# Add a load to each bus
n_loads = n_buses
for l in range(n_loads):
    n.add("Load", f"load{l}", bus=f"bus{l}", p_set=pd.Series(100, index=n.snapshots))

# Add a generator to the first bus
n.add(
    "Generator",
    "gen0",
    bus="bus0",
    p_nom_extendable=True,
    capital_cost=1000,
    marginal_cost=1,
)

# Add a slack unit to each node: needed when using investmentblock
for b in range(n_buses):
    n.add(
        "Generator",
        f"slack{b}",
        bus=f"bus{b}",
        carrier="slack",
        marginal_cost=1000,
        p_nom=1e6,
    )

n = pypsa.Network() n.set_snapshots(pd.date_range("2024-01-01T00:00", "2024-01-01T23:00", freq="h")) # Add carriers n.add("Carrier", "AC") # Add buses n_buses = 2 for b in range(n_buses): n.add("Bus", f"bus{b}", carrier="AC") # Add lines in a radial topology using bidirectional links n_lines = n_buses - 1 for l in range(n_lines): n.add( "Link", f"line{l}", bus0=f"bus{l}", bus1=f"bus{l+1}", length=1, capital_cost=1000, p_min_pu=-1, p_nom_extendable=True, ) # Add a load to each bus n_loads = n_buses for l in range(n_loads): n.add("Load", f"load{l}", bus=f"bus{l}", p_set=pd.Series(100, index=n.snapshots)) # Add a generator to the first bus n.add( "Generator", "gen0", bus="bus0", p_nom_extendable=True, capital_cost=1000, marginal_cost=1, ) # Add a slack unit to each node: needed when using investmentblock for b in range(n_buses): n.add( "Generator", f"slack{b}", bus=f"bus{b}", carrier="slack", marginal_cost=1000, p_nom=1e6, )

In [3]:

n.optimize(solver_name="highs")

n.optimize(solver_name="highs")

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

INFO:linopy.model: Solve problem using Highs solver

INFO:linopy.io: Writing time: 0.04s

INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 98 primals, 242 duals
Objective: 3.05e+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, Generator-ext-p-lower, Generator-ext-p-upper, Link-ext-p-lower, Link-ext-p-upper were not assigned to the network.

Running HiGHS 1.11.0 (git hash: 364c83a): Copyright (c) 2025 HiGHS under MIT licence terms
LP   linopy-problem-2wgduix1 has 242 rows; 98 cols; 386 nonzeros
Coefficient ranges:
  Matrix [1e+00, 1e+00]
  Cost   [1e+00, 1e+03]
  Bound  [0e+00, 0e+00]
  RHS    [1e+02, 1e+06]
Presolving model
96 rows, 74 cols, 216 nonzeros  0s
Dependent equations search running on 24 equations with time limit of 1000.00s
Dependent equations search removed 0 rows and 0 nonzeros in 0.00s (limit = 1000.00s)
96 rows, 74 cols, 216 nonzeros  0s
Presolve : Reductions: rows 96(-146); columns 74(-24); elements 216(-170)
Solving the presolved LP
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0    -1.0231599545e+01 Pr: 48(7200) 0s
         28     3.0480000000e+05 Pr: 0(0) 0s
Solving the original LP from the solution after postsolve
Model name          : linopy-problem-2wgduix1
Model status        : Optimal
Simplex   iterations: 28
Objective value     :  3.0480000000e+05
P-D objective error :  0.0000000000e+00
HiGHS run time      :          0.00

Out[3]:

('ok', 'optimal')

2. Define the configuration to use InvestmentBlock¶

In [4]:

# Create the transformation object and run the transformation
# The `capacity_expansion_ucblock` option is set to `False` to enable capacity expansion using InvestmentBlockSolver instead of UnitCommitmentBlockSolver
tran = pypsa2smspp.Transformation(
    capacity_expansion_ucblock=False,
)

# Create the transformation object and run the transformation # The `capacity_expansion_ucblock` option is set to `False` to enable capacity expansion using InvestmentBlockSolver instead of UnitCommitmentBlockSolver tran = pypsa2smspp.Transformation( capacity_expansion_ucblock=False, )

3. Convert the PyPSA model to SMS++ and optimize it¶

In [5]:

n_smspp = tran.run(n, verbose=False)  # create the model and optimizes it in one line

# If you wish to create the model and optimize it in separate steps, you can do so as follows:
# tran.create_model(n)  # create the model
# tran.optimize(solver_name="highs")  # optimize the model
# n_smspp = tran.retrieve_solution(n)  # retrieve the solution and store it in a new Network object

n_smspp = tran.run(n, verbose=False) # create the model and optimizes it in one line # If you wish to create the model and optimize it in separate steps, you can do so as follows: # tran.create_model(n) # create the model # tran.optimize(solver_name="highs") # optimize the model # n_smspp = tran.retrieve_solution(n) # retrieve the solution and store it in a new Network object

Executing command:
InvestmentBlock_test temp_test_case.nc -S BSPar.txt -c /home/docs/checkouts/readthedocs.org/user_builds/pypsa2smspp/conda/latest/lib/python3.13/site-packages/pysmspp/data/configs/InvestmentBlock/ -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_test_case.nc -v 1

temp_test_case.nc is a Block file

{1-0-0-0.0001} t = 1.00e-10 ~ D*_1( z* ) = 0.00e+00 ~ Sigma = 0.00e+00
            Fi undefined
    Lambda1 = [ 0.00e+00 0.00e+00 ]
    UB[ 0 ] = INF, LB[ 0 ] = -INF
            Fi[ 0 ]: UB = 4.8000000000e+06, LB = 4.8000000000e+06 [0.0014] 
            New subgradient for Fi[ 0 ] ~ Alfa1 = -4.80e+06 ~ gd = -0.00e+00 stored in 0 (0)
            [0.0016] Fi1 = 4.8000000000e+06
            Fi1 defined ==> SS 

{1-1-1-0.0020} t = 1.00e+00 ~ D*_1( z* ) = 2.87e+08 ~ Sigma = 0.00e+00
            Fi = 4.8000000000e+06 ~ eU = 5.99e+03
    Lambda1 = [ 2.40e+04 0.00e+00 ]
    UB[ 0 ] = INF, LB[ 0 ] = -5.7004857600e+08
            Fi[ 0 ]: UB = 2.6378400000e+07, LB = 2.6378400000e+07 [0.0005] 
            New subgradient for Fi[ 0 ] ~ Alfa1 = 2.40e+06 ~ gd = 2.40e+07 stored in 1 (1)
            [0.0024] Fi1 = 2.6378400000e+07

 NS[0]: DFi = 2.16e+07 ~ Lw1(2.6378400000e+07) >= LwTrgt(-3.9759400320e+08)

{1-2-2-0.0035} t = 1.00e+00 ~ D*_1( z* ) = 1.38e+08 ~ Sigma = 1.19e+06
            Fi = 4.8000000000e+06 ~ eU = 2.87e+03
    Lambda1 = [ 1.15e+04 1.19e+04 ]
    UB[ 0 ] = INF, LB[ 0 ] = -2.7208387931e+08
            Fi[ 0 ]: UB = 2.3483217388e+07, LB = 2.3483217388e+07 [0.0006] 
            New subgradient for Fi[ 0 ] ~ Alfa1 = 4.80e+06 ~ gd = 2.35e+07 stored in 2 (2)
            [0.0033] Fi1 = 2.3483217388e+07

 NS[1]: DFi = 1.87e+07 ~ Lw1(2.3483217388e+07) >= LwTrgt(-1.8901871551e+08)

{1-3-3-0.0050} t = 1.00e+00 ~ D*_1( z* ) = 2.23e+04 ~ Sigma = 4.47e+06
            Fi = 4.8000000000e+06 ~ eU = 1.40e+00
    Lambda1 = [ 1.88e+02 9.60e+01 ]
    UB[ 0 ] = INF, LB[ 0 ] = 2.8894493288e+05
            Fi[ 0 ]: UB = 5.7308986577e+05, LB = 5.7308986577e+05 [0.0007] 
            New subgradient for Fi[ 0 ] ~ Alfa1 = 0.00e+00 ~ gd = -4.23e+06 stored in 3 (3)
            [0.0043] Fi1 = 5.7308986577e+05

 SS[0]: DFi = -4.23e+06 ~ Up1(5.7308986577e+05) <= UpTrgt(4.3488944933e+06)

{1-4-4-0.0064} t = 1.50e+00 ~ D*_1( z* ) = 3.12e+01 ~ Sigma = 2.72e+05
            Fi = 5.7308986577e+05 ~ eU = 4.80e-01
    Lambda1 = [ 2.00e+02 9.60e+01 ]
    UB[ 0 ] = INF, LB[ 0 ] = 3.0079615631e+05
            Fi[ 0 ]: UB = 3.9679231262e+05, LB = 3.9679231262e+05 [0.0006] 
            New subgradient for Fi[ 0 ] ~ Alfa1 = 1.88e+05 ~ gd = 1.19e+04 stored in 4 (4)
            [0.0049] Fi1 = 3.9679231262e+05

 SS[0]: DFi = -1.76e+05 ~ Up1(3.9679231262e+05) <= UpTrgt(5.4586049482e+05)

{1-5-5-0.0076} t = 2.25e+00 ~ D*_1( z* ) = 1.58e+00 ~ Sigma = 9.20e+04
            Fi = 3.9679231262e+05 ~ eU = 2.32e-01
    Lambda1 = [ 2.00e+02 1.00e+02 ]
    UB[ 0 ] = INF, LB[ 0 ] = 3.0480000000e+05
            Fi[ 0 ]: UB = 3.0480000000e+05, LB = 3.0480000000e+05 [0.0009] 
            New subgradient for Fi[ 0 ] ~ Alfa1 = 9.47e-09 ~ gd = -9.20e+04 is copy of 4 (4)
            [0.0061] Fi1 = 3.0480000000e+05

 SS[0]: DFi = -9.20e+04 ~ Up1(3.0480000000e+05) <= UpTrgt(3.8759308136e+05)

{1-6-5-0.0094} t = 3.38e+00 ~ D*_1( z* ) = 2.07e-25 ~ Sigma = -9.47e-09
            Fi = 3.0480000000e+05 ~ eU = -3.11e-14 ~ stop (optimal)

Call 1: 7 ~ 6 ~ 0.0097 ~ 0.0061 -> optimal ~ Fi* = 3.0480000000e+05
Solver status: 10
Peak CPU Usage: 118.90 %
Peak Memory Usage: 51.66 MB
Total Time: 0.31 seconds

[split] No merged DCNetworkBlock entries found; nothing to do.

Objective value

In [6]:

n_smspp.objective

n_smspp.objective

Out[6]:

304800.0

4. Verify the equivalence of SMS++ and PyPSA results¶

In [7]:

pypsa_obj = n.objective
smspp_obj = n_smspp.objective
print("SMS++ obj         : %.6f" % smspp_obj)
print("PyPSA obj         : %.6f" % pypsa_obj)
print("Error SMS++ - PyPSA [%%]: %.5f" % (100*(smspp_obj - pypsa_obj)/pypsa_obj))

pypsa_obj = n.objective smspp_obj = n_smspp.objective print("SMS++ obj : %.6f" % smspp_obj) print("PyPSA obj : %.6f" % pypsa_obj) print("Error SMS++ - PyPSA [%%]: %.5f" % (100*(smspp_obj - pypsa_obj)/pypsa_obj))

SMS++ obj         : 304800.000000
PyPSA obj         : 304800.000000
Error SMS++ - PyPSA [%]: 0.00000

4. Analyze the PyPSA model created with SMS++ and pypsa2smspp¶

Retrieve PyPSA statistics

In [8]:

n_stat = n_smspp.statistics.optimal_capacity()
n_stat

n_stat = n_smspp.statistics.optimal_capacity() n_stat

Out[8]:

component  carrier
Generator  -              200.0
           slack      2000000.0
Link       AC             100.0
dtype: float64

Retrieve SMS++ statistics

In [9]:

n_parsed_stat = n_smspp.statistics.optimal_capacity()
n_parsed_stat

n_parsed_stat = n_smspp.statistics.optimal_capacity() n_parsed_stat

Out[9]:

component  carrier
Generator  -              200.0
           slack      2000000.0
Link       AC             100.0
dtype: float64

Calculate difference between the two

In [10]:

error_stat = n_stat - n_parsed_stat

merged_stat = pd.concat([n_stat, n_parsed_stat, error_stat], axis=1)
merged_stat.columns = ["pypsa", "smspp", "difference"]
merged_stat

error_stat = n_stat - n_parsed_stat merged_stat = pd.concat([n_stat, n_parsed_stat, error_stat], axis=1) merged_stat.columns = ["pypsa", "smspp", "difference"] merged_stat

Out[10]:

		pypsa	smspp	difference
component	carrier
Generator	-	200.0	200.0	0.0
	slack	2000000.0	2000000.0	0.0
Link	AC	100.0	100.0	0.0