[breaking] Configuration templates for the BundleSolver 2.0

[dmeoli]

This PR brings the configuration templates of pySMSpp to BundleSolver 2.0, which is now in the develop of SMS++, and adds one template per way of solving a two-stage stochastic Block (TSSB) or a multi-stage one (MSSB) with the Solvers of SMS++. The six templates that use the BundleSolver lose the parameters of the master of the 1.0 (intMPName, intMPlvl, intQPmp1, intQPmp2, intOSImp1-3, dblCtOff), since any of them makes a ComputeConfig of the 2.0 fail to load, and they name the Solver of the new master in MPBCfg.txt, one per folder (HiGHS by default, Gurobi in MPBCfg_grb.txt). In exchange, with the SMS++ of conda-forge, which still ships the 1.0 until the next release, the templates that use the BundleSolver do not load (they stop at "Invalid int parameter name intMPStbl"). To say which template fits which network we measured all of them against PyPSA with Gurobi on three families of networks written by pypsa2smspp, and the same networks, in smaller sizes, are now in the test batteries of SMS++.

The structure of this description is as follows. We first list the new templates, and then describe the forms in which pypsa2smspp writes a network and the setting of each Solver, since the comparison depends on both. Then we report the measurements family by family, the nested forms (a decomposition solving the subproblems of another one) included, and draw from them the choice of the form and of the Solver for a class of networks. Finally, we list the batteries that check the same networks, the limits of what was measured, and the commands to reproduce it.

Templates

All the new templates are in TSSBlock/.

template	Solver	what it gives	tool option
TSSBSCfg-IP.txt	:MILPSolver on the deterministic equivalent	the optimum, to the gap asked of Gurobi
TSSBSCfg-LP.txt	the same on its continuous relaxation	a lower bound, in seconds
TSSBSCfg-PIPS.txt	PIPS-IPM++ on the same linear program	the same bound, spread over one process per group of blocks	mpirun -np <n>
TSSBSCfg-LD-IP.txt	Lagrangian dual of the scenarios, each scenario a MILP	a lower bound, the strongest of the three duals (on the thermal family, the optimum to 6 digits)
TSSBSCfg-LDLD.txt	a Lagrangian dual per scenario inside the one over the scenarios	the bound of the convex hull of each unit, weaker than the one above, on instances the MILP does not close
TSSBSCfg-LDrec.txt	recursive Lagrangian dual, scenarios and units at once	the same bound as LD+LD, usually the fastest of the three on long horizons
*-par.txt	the same three duals with a ParallelBundleSolver, 8 threads	the same bounds to all the printed digits, usually in less time
TSSBSCfg-PPH.txt	PrimalProximalHeur	a feasible solution and its distance from the bound
TSSBSCfg-BDS.txt	BendersDecompositionSolver, one cut per scenario, 8 threads (BDSCfg.txt)	the optimum, exactly	-k
TSSBSCfg-BDS-LD.txt	the same, each subproblem solved by the recursive Lagrangian dual	the optimum where the subproblems are linear programs; otherwise the bound of the dual and a feasible solution, recovered at the design of the master (strRecoveryBSC in BDSCfg-LD.txt)	-k
TSSBSCfg-BDS-CVX-LD.txt	the same with the master solved by a BundleSolver (BDSMCfg-CVX.txt)	idem	-k
TSSBSCfg-BDS-CVX-int.txt	BDS with the master solved by a BundleSolver that keeps an integer design integer (intIntVars, BDSMCfg-CVX-int.txt)	the optimum with the integer design, the same as TSSBSCfg-BDS.txt	-k
TSSBSCfg-IB.txt	InvestmentBlock over the UCBlock, the TSSB or the MSSB	the optimum of the problem with a continuous design	smspp_investmentblock_solver, with -B InnerBCfg-LD.txt the inner Block solved by the recursive Lagrangian dual

The integer design is also available in an InvestmentBlock, whose netCDF variable Integer marks the assets whose design is integer (e.g., the number of modules of a modular asset, which pypsa2smspp writes from p_nom_mod in SPSUnipi/pypsa2smspp#62): InvestmentBlock/BSPar-int.txt is the BundleSolver of BSPar.txt keeping it integer, i.e., the stabilized cutting-plane method of van Ackooij, Frangioni and de Oliveira, with a trust region of radius 1 and a master that is a mixed-integer linear program (MPBCfg-IP.txt), which HiGHS solves; the BDS template above uses the same master. On mod_t168_s5_b1 the BDS template gives 8.26022e7, as TSSBSCfg-BDS.txt and PyPSA (82602180.128), while with the continuous design it gives 8.25888e7; the InvestmentBlock gives the value of PyPSA on mod_1n_2c_2g of pypsa2smspp.

All the thermal templates need -B InnerBCfg.txt, and -k is the option of smspp_tssb_solver (develop of SMS++) that builds the Benders form. As for anything else a Solver does, a feasible solution is asked for in the configuration: the PrimalProximalHeur of TSSBSCfg-PPH.txt and the BDS of TSSBSCfg-BDS-LD.txt and of TSSBSCfg-BDS-CVX-LD.txt recover it with the BlockSolverConfig in their strRecoveryBSC (here BSCfg1-IP.txt, i.e., each scenario a MILP at a fixed design), and report its value as their upper bound. For BDS the fixed design is that of the master, which with a bundle master is the best point of the bundle.

Forms and setting of the Solvers

pypsa2smspp writes the netCDF file in one of several forms, selected by a few of its knobs. capacity_expansion_ucblock puts the design in the units (the form that BDS and the MILP take), and investment_outside states it once in an InvestmentBlock above them. Also, stochastic_type is tssb or mssb (with tree grouping the scenarios), and p_nom_mod on a generator makes its design integer. We set the Solvers compared below as follows.

• PyPSA writes the model with linopy and passes it to Gurobi with its defaults, MIPGap=1e-6 apart, i.e., the concurrent method at the root (primal and dual simplex racing the barrier, crossover on) and then its branch and bound.
• MILP SMS++ is the GRBMILPSolver on the model written by the Block, with the same defaults and the same gap (dblRelAcc). Its option intCutSepPar 7 gives Gurobi a callback, and with it PreCrush and LazyConstraints, which keep more columns after the presolve; in our runs they make it faster: without the callback (intCutSepPar 0), and on the same unloaded machine, the same runs take 1.4 to 5.4 times longer (129.3 against 93.0 s on u80_t96_s5, 465.9 against 87.4 s on s200_b10), and on s200_b20 and u80_t168_s3 they do not end within the hour.
• The ThermalUnitBlock has 7 formulations of the commitment of a unit (the first value of TUBCfg.txt: 0 = 3bin, 1 = T, 2 = pt, 3 = DP, 4 = SU, 5 = SD, 6 = SUSD). Our templates use 3bin, which in these runs is usually the faster of the two that Gurobi closes. Against T it takes 12.1 against 16.7 s on u80_t48_s3, 92.6 against 122.8 s on u80_t96_s5 and 2229.0 against more than 3600 s on u80_t168_s3, while it is 6% slower on u40_t96_s10 (78.1 against 73.5 s), with the same optimum. In principle the continuous relaxation of 3bin is weaker than that of T; however, on these instances the two give the same value (5.39874e6 on u10_t24_s3), and the root relaxation of Gurobi brings both to the bound of LD+LD (cf. the table of values below). Of the other five, only DP describes the convex hull of the feasible set of a unit, and hence only its continuous relaxation is in general as strong as the Lagrangian dual. Instead, pt, SU, SD and SUSD approximate it with different trade-offs between size and bound (on u10_t24_s3 the four happen to give the value of DP, 5.40782e6, which need not hold on other instances). Their models are much larger (15.8 million rows for DP against 110 thousand for T on u80_t48_s3): Gurobi does not close DP, SU, SD and SUSD in 1800 s, and pt takes 529.7 work units against the 20.3 of T.
• In LD+LD, in the recursive LD and in PPH each unit is solved by the dynamic programming of ThermalUnitExtDPSolver (TUBSCfg-DP.txt); ThermalUnitDPSolver is as fast within the noise of the machine, and takes twice the memory. In all of them the bundle is that of BundleSolver 2.0, sequential or a ParallelBundleSolver with 8 threads (-par), and its master is solved by Gurobi in all the measurements below (MPBCfg_grb.txt).
• BDS uses one optimality cut per scenario with a MILP master, which here is the fastest of its configurations (a single aggregated cut costs 13% to 150% more). Its convex regime and the Pareto-optimal (Papadakos), deepest (Hosseini-Turner) and static (Brandenberg-Stursberg) cuts are options of the Solver. With the static cut BDS makes the fewest rounds, but its separation problem is degenerate; choosing its vertex by the perturbation of Sherali and Lunday (int_BDSlv_Pareto 2, dbl_BDSlv_ParetoMu 0.1) takes it from 42-73 to 24-56 rounds, and from 120-201 s to 93-105 s, on the four capacity-expansion instances of PyPSA, which is still about 10 times the multi-cut standard.
• InvestmentBlock is a BundleSolver on the design over a UCBlock, a TSSB or an MSSB, whose inner Block is solved in full by Gurobi, as a linear program, at each evaluation (BSCfg-IB.txt).
• PIPS-IPM++ solves linear programs only, with one MPI rank per group of blocks. Its first stage is the Variable of the root of the Block: on the extensive form of a TSSB it is empty, and the non-anticipativity becomes linking Constraint, while on the Benders form of -k it is the design, which is the case it is made for. In fact, each Modification rebuilds its model, and therefore it is best suited to a large LP solved once rather than inside a decomposition. On s50_b10c the setting took it from 477 s to 94 s, at 46 to 51 iterations: OMP_NUM_THREADS=1 per rank (without it each rank takes all the cores, and the run is 3 times slower), MUMPS as the linear solver (PARDISO of MKL hangs at the root with 25 or more ranks), 25 ranks and the form of -k. Whether it converges at all depends on the scaler: SCALER equilibrium, that of PIPSCfg.txt, converges on all the instances here, and it is 7 to 14% faster than Curtis-Reid where both converge, whereas Curtis-Reid does not converge on the relaxation of a unit commitment (it stops at its limit of 300 interior point iterations without a bound).
• On its LPs PyPSA-Eur gives Gurobi other options (barrier without crossover, BarConvTol=1e-5, PreDual=0, AggFill=0, TSSBSCfg_grb.txt), which may pay, at least on the modular family: with them, on s200_b10c PyPSA takes 85.3 s and MILP SMS++ 50.7 s, and the optimum is 8e-6 higher in relative terms, since the barrier stops without a crossover. They are not used in the tables below.

The parameters one most often changes are all in the templates: dblRelAcc is the relative gap of a :MILPSolver (the MIPGap of Gurobi) and dblMaxTime its time limit, and intMaxThread the threads of BDS and those of a ParallelBundleSolver. Also, strRecoveryBSC is the BlockSolverConfig of the recovery of a feasible solution, in the PrimalProximalHeur and in BDS, intRecoveryThreads the threads of the recovery of the PrimalProximalHeur, and the first value of TUBCfg.txt the formulation of a thermal unit.

Computational results

We ran everything on one machine with 2 AMD EPYC 9654 (192 cores) and 1.5 TB, with Gurobi 13.0.2 in PyPSA 1.2.4 and 13.0.1 in SMS++, and with the BundleSolver 2.0 of 20/09 (2b77538). The InvestmentBlock was measured with that of 25/09 (86ba410), and the nested forms with the develop of 27/09 (umbrella e6ba4025). We made the runs one at a time, with a limit of 3600 s each (1800 s for the nested forms); the load average of the machine stayed under 60 in all of them, and each run records the load when it ends. A time is that of the complete process in seconds, which includes the reading of the network or of the netCDF file but not the conversion by pypsa2smspp, and "> 3600" means that the run did not end within the limit. Where a cell holds two numbers separated by "/", the second is the peak memory in GiB (for PIPS-IPM++, that of its largest process, the others being about as large). PyPSA and MILP SMS++ set in Gurobi the same relative gap, 1e-6, and reach the same optimum to all the digits Gurobi prints. In the tables, "PyPSA" is Gurobi on the model of PyPSA, and "MILP SMS++" and "LP SMS++ (Gurobi)" are the GRBMILPSolver on the model of SMS++ and on its continuous relaxation, or on the model itself where it is a linear program. "LP SMS++ (PIPS)" is PIPS-IPM++ on the latter, and the other columns are the Solvers of the templates.

Thermal family

gen_thermal_tssb.py generates a network with 1 bus (units x periods x scenarios in the name). Since PyPSA cannot optimize a network that is at once stochastic and committable, extensive_thermal_tssb.py writes its deterministic equivalent for PyPSA, i.e., one copy of the network per scenario, with the costs of each copy weighted by the probability of its scenario and the design tied across the copies. This is the same two-stage problem that SMS++ solves, and indeed the two reach the same optimum; however, PyPSA does not offer this form, and hence the column "PyPSA" of this family is Gurobi on the extensive form written with PyPSA rather than a path of PyPSA. We use this family mainly for the Lagrangian duals and the PPH, which have no counterpart on an extensive form. On this family the Solvers do not compute the same thing: the MILPs give the optimum, the three Lagrangian Duals (LD) lower bounds of different strength, the LPs the bound of the continuous relaxation, and the PrimalProximalHeur (PPH) feasible solutions. Hence we report the values first, and we compare the times only among the Solvers that deliver the same value. Values are given to the 6 digits SMS++ prints; the optimum is that of Gurobi at a relative gap of 1e-6, the same for PyPSA and MILP SMS++.

instance	LP of PyPSA	LP of SMS++	LD+LD and recursive LD	LD	optimum
u80_t48_s3	1.04606e7	1.04614e7	1.04698e7	1.04711e7	1.04711e7
u80_t96_s5	1.98643e7	1.98646e7	1.98738e7	1.98789e7	1.98789e7
u40_t96_s10	1.97684e7	1.97687e7	1.97759e7	1.97877e7	1.97877e7
u80_t168_s3	2.85920e7	2.85920e7	2.86008e7	2.86133e7	2.86138e7
u80_t336_s3	5.12560e7	5.12561e7	5.12639e7	(> 3600 s)	(> 3600 s)

LD+LD and the recursive LD give the same bound, i.e., that of the relaxation in which each unit is replaced by the convex hull of its feasible set. The LD of the scenarios keeps the units integer inside a scenario, and hence gives a stronger bound, which on these instances is the optimum to 6 digits (0.0016% below it on u80_t168_s3). In turn, the LP of PyPSA is only slightly weaker than that of SMS++, since the two models write the constraints of the units somewhat differently, and both are 0.07 to 0.1% below the optimum. Note that the root relaxation of Gurobi (after its presolve, before its cuts) is already at the bound of LD+LD on the model of SMS++ (1.046984e7 on u80_t48_s3, 2.860082e7 on u80_t168_s3), which matters for the next table.

Times are grouped by the value delivered, and the column "PyPSA to the LD+LD bound" is Gurobi on the model of PyPSA, stopped as soon as its bound reaches that of LD+LD (BestBdStop); the peak memory of the same runs follows in the second table.

instance	PyPSA	MILP SMS++	LD	LD+LD	recursive LD	PyPSA to the LD+LD bound	LP SMS++ (Gurobi)	LP SMS++ (PIPS)
u80_t48_s3	27.1	12.3	6.6	3.4	1.8	12.0	1.4	28.0
u80_t96_s5	158.0	93.4	30.0	24.4	17.6	25.1	5.6	220.5
u40_t96_s10	142.3	78.1	75.5	38.5	43.6	21.9	4.3	82.3
u80_t168_s3	2839.5	2286.5	1654.7	85.1	27.7	21.0	7.0	431.1
u80_t336_s3	> 3600	> 3600	> 3600	586.7	123.7	40.1	15.7	2002.2

Peak memory in GiB of the same runs.

instance	PyPSA	MILP SMS++	LD	LD+LD	recursive LD	PyPSA to the LD+LD bound	LP SMS++ (Gurobi)	LP SMS++ (PIPS)
u80_t48_s3	3.1	1.2	0.7	0.2	0.2	1.2	0.4	0.2
u80_t96_s5	6.4	4.4	2.2	0.4	0.4	2.1	1.0	0.5
u40_t96_s10	14.8	1.9	1.2	0.4	0.5	2.1	1.0	0.3
u80_t168_s3	28.9	38.9	9.9	0.5	0.4	2.2	1.1	0.9
u80_t336_s3	51.0	60.7	21.2	0.8	0.8	3.8	2.1	2.2

Among the Solvers that give the optimum, the LD of the scenarios is 1.03 to 3.1 times faster than MILP SMS++ and takes 1.6 to 3.9 times less memory, while MILP SMS++ is 1.2 to 2.2 times faster than PyPSA. The bound of LD+LD and of the recursive LD, instead, is reached sooner by Gurobi on 3 instances out of 5 (u40_t96_s10, u80_t168_s3 and u80_t336_s3, in 21.0 to 40.1 s against 27.7 to 586.7 s), which the root relaxation noted above explains, at least in part. Hence, on this family the two duals are not a faster way to that bound; they win on memory (0.2 to 0.8 GiB against 1.2 to 3.8 GiB), and they are the starting point of the feasible solutions below. With a ParallelBundleSolver of 8 threads (-par) the same bounds come out, to all the printed digits, in 3.5, 9.1, 47.4, 780.1 and > 3600 s for LD, in 1.4, 6.4, 11.0, 52.4 and 311.3 s for LD+LD, and in 1.5, 13.0, 35.2, 20.8 and 102.2 s for the recursive LD. Threads pay where the components are few and large (2.1 to 3.8 times on LD+LD), and much less where they are many and small (1.2 to 1.4 times on the recursive LD, whose master is the bottleneck). PIPS-IPM++, with one rank per scenario and the form of -k, gives the value of Gurobi on the same relaxation to 6 digits, but 20 to 130 times more slowly.

For the feasible solutions, PPH starts from the primal solution of its LD and fixes the design of each scenario to its mean over the scenarios (rounded where it is integer). It then solves each scenario by itself with the BlockSolverConfig of its strRecoveryBSC (here a MILP, intRecoveryThreads scenarios at a time), and it reports the sum of their values, i.e., the value of a feasible solution, as its upper bound.

instance	PyPSA	MILP SMS++	PPH
u80_t48_s3	27.1 / 3.1	12.3 / 1.2	7.2 / 1.2
u80_t96_s5	158.0 / 6.4	93.4 / 4.4	33.4 / 4.5
u40_t96_s10	142.3 / 14.8	78.1 / 1.9	43.6 / 3.3
u80_t168_s3	2839.5 / 28.9	2286.5 / 38.9	150.8 / 8.3
u80_t336_s3	> 3600 / 51.0	> 3600 / 60.7	1809.8 / 41.1

On the first three instances PPH gives the optimum to the digits Gurobi prints. On u80_t168_s3 it is 0.0005% above the optimum (2.86139e7 against 2.8613751e7), in 150.8 s against the 2286.5 s of MILP SMS++ and the 2839.5 s of PyPSA. On u80_t336_s3, which neither MILP closes in an hour, it gives 5.12783e7, i.e., 0.028% above the bound of LD+LD; there most of the time goes into the recovery, since each scenario is a MILP of 80 units over two weeks. It is fair to remark that some runs on u80_t336_s3 were interrupted by the license service of Gurobi, which allows 2 sessions at a time. The table reports the runs that were not, which give the same bound to all digits.

Modular family, integer design

gen_modular_tssb.py builds modules of solar and wind in the first stage, over 168 periods (scenarios x buses in the name). BDS runs on the Benders form that -k builds, with 1, 8 or 32 threads (intMaxThread), and it reaches the optimum exactly in the same number of rounds with any number of threads; "PyPSA, 1 thread" is PyPSA with Gurobi restricted to one thread.

instance	PyPSA	MILP SMS++ (Gurobi)	PyPSA, 1 thread	BDS, 1 thread	BDS, 8 threads	BDS, 32 threads
s50_b10	665.7 / 47.8	846.8 / 52.1	> 3600 / 3.5	51.6 / 1.4	8.9 / 1.4	4.6 / 1.5
s100_b10	59.4 / 5.6	38.8 / 4.5	> 3600 / 3.4	101.0 / 2.7	17.0 / 2.8	8.4 / 2.8
s200_b10	134.7 / 10.0	83.7 / 8.9	> 3600 / 6.3	193.0 / 5.3	33.2 / 5.4	15.3 / 5.6
s100_b20	185.8 / 10.4	189.5 / 9.0	> 3600 / 6.1	272.0 / 5.2	43.6 / 5.5	19.8 / 5.9
s200_b20	> 3600 / 18.0	478.6 / 18.0	> 3600 / 12.1	540.2 / 10.4	84.4 / 10.7	36.7 / 11.3

Since the rounds of BDS do not depend on the threads, we see 11 to 15 times less time going from 1 to 32 threads; on s200_b20 BDS reaches the optimum in 36.7 s, where PyPSA does not close in an hour with all the threads of the machine. With one thread each, i.e., with the same resources, PyPSA does not close any of the five in an hour at the same gap of 1e-6, while BDS takes 52 to 540 s. The exception is s50_b10, the one cell where SMS++ is the slower of the two: both PyPSA and MILP SMS++ find the optimum at the root in about 30 s, and they spend the rest closing the last millionth of the gap (43347 nodes).

Modular family, continuous design, TSSB and MSSB

With the modules turned off (c in the name) the same networks are linear programs, so that all the Solvers of SMS++ for a stochastic investment apply, PIPS-IPM++ included. Each is written as a TSSB and as an MSSB whose scenarios are grouped into 5 or 10 outer realizations, with the same deterministic equivalent. BDS is the generic Benders decomposition on the leaves, with 8 threads, and the InvestmentBlock is the ad hoc one, which solves the TSSB or the MSSB in full at each evaluation; PyPSA has no MSSB, hence its cells are empty there.

instance	form	PyPSA	LP SMS++ (Gurobi)	LP SMS++ (PIPS)	BDS, 8 threads	InvestmentBlock
s50_b10c	TSSB	26.1 / 2.6	13.3 / 2.1	44.9 / 0.7	8.2 / 1.4	302.9 / 2.6
s50_b10c	MSSB		32.8 / 2.2	46.9 / 0.8	8.6 / 1.5	340.9 / 2.8
s100_b10c	TSSB	50.4 / 4.7	30.5 / 4.3	85.6 / 1.3	16.6 / 2.8	1081.3 / 3.8
s100_b10c	MSSB		98.3 / 4.4	84.8 / 1.5	16.8 / 2.9	1099.9 / 4.0
s200_b10c	TSSB	102.9 / 8.1	60.7 / 8.6	169.0 / 2.6	32.7 / 5.4	> 3600 / 8.5
s200_b10c	MSSB		60.1 / 8.6	170.8 / 2.7	34.5 / 5.6	> 3600 / 8.5

BDS is the fastest of the four on all the instances here. On the tree and on the flat problem it makes the same rounds, and both give the same optimum; also the LP takes the same time on both (10.5 to 52.1 s of Gurobi), the difference in the column coming from the reading of the tree. PIPS-IPM++ takes about 3 times the time of Gurobi on the same linear program (3.4 at 50 scenarios, 2.8 at 100 and 200 on the TSSB), with its memory spread over 25 or 50 processes of 0.7 to 2.6 GiB each. At each evaluation the InvestmentBlock solves the stochastic Block (about 15 s on s50, 50 s on s100 and 190 s on s200). It closes s50 and s100 in 18 to 21 iterations, i.e., in 37 to 65 times the time of BDS, the tree costing it no more than the flat problem, and on s200 the hour ends with its value 0.14% (TSSB) and 0.02% (MSSB) above the optimum. Its stopping test is relative to the norm of the first full subgradient (intWZNorm 10 in TSSBSCfg-IB.txt), since as an absolute one the threshold on the aggregated subgradient is below what it reaches on a design in MW, and the bundle would go on long after the optimum.

Deterministic capacity expansion

One scenario of the modular family over a longer horizon (continuous design, periods and buses in the name), where the InvestmentBlock wraps a UCBlock. We compare PyPSA, the LP of the UCBlock with its design in the units, PIPS-IPM++ on the same LP (whose blocks are the units), BDS on the Benders form with one subproblem, and the InvestmentBlock over the UCBlock.

instance	PyPSA	LP SMS++ (Gurobi)	LP SMS++ (PIPS)	BDS, 8 threads	InvestmentBlock
t720_b10c	6.2 / 0.5	1.1 / 0.2	292.5 / 2.2	14.4 / 0.2	2.8 / 0.2
t2016_b10c	8.2 / 0.9	2.5 / 0.6	> 3600	74.1 / 0.3	7.5 / 0.5
t2016_b20c	14.1 / 1.4	5.7 / 1.0	> 3600	167.4 / 0.6	12.6 / 0.9

Here the ranking is the opposite of that on the stochastic instances above. With one subproblem BDS has no evaluations to run side by side, and it is the slowest of the four; the InvestmentBlock, which closes in 8 to 10 iterations, is instead as fast as PyPSA or faster. PIPS-IPM++ converges on the smallest of the three only, since the blocks of a UCBlock are its units, i.e., many small ones instead of the few large ones its decomposition is made for.

Nested forms

A decomposition can solve the subproblems of another one. The nested forms that we found to make sense here are BDS whose subproblems are solved by the recursive LD (TSSBSCfg-BDS-LD.txt), the same with the master solved by a bundle (TSSBSCfg-BDS-CVX-LD.txt, i.e., the convex regime), and the InvestmentBlock whose inner Block is solved by the recursive LD (-B InnerBCfg-LD.txt instead of InnerBCfg.txt). For comparison, we also give the rows of BDS and of the InvestmentBlock with Gurobi inside.

instance	form	value	time	memory
u80_t48_s3 (thermal)	BDS, recursive LD in the subproblems	1.04698e7	1.1	0.2
t48_s10_b2c (modular, continuous)	BDS	3.68235e7	0.17	0.06
	BDS, recursive LD in the subproblems	3.68235e7	152.2	0.2
	BDS with a bundle master, recursive LD in the subproblems	3.68235e7	75.9	0.2
	InvestmentBlock	3.68250e7	1.1	0.06
	InvestmentBlock, recursive LD inside	3.68235e7	277.7	0.4
s50_b10c (modular, continuous)	BDS	6.19693e7	8.4	1.4
	BDS, recursive LD in the subproblems		> 1800	4.3
	BDS with a bundle master, recursive LD in the subproblems		> 1800	3.1
	InvestmentBlock	6.19698e7	224.7	1.7
	InvestmentBlock, recursive LD inside		> 1800	152

On the thermal instance, whose subproblems are unit commitments, BDS with the recursive LD gives the bound of LD+LD in about a second, since the units are solved by their dynamic programming; this bound (1.04698e7) is above the 1.04614e7 that BDS gives with the relaxed recourse (cf. the limits below). With strRecoveryBSC it also recovers a feasible solution at the design of its master, each scenario solved as a MILP, with either master (the bundle one at the best point of its bundle). On u80_t48_s3 both give the optimum (1.0471130e7), i.e., 0.012% above the bound, as PPH does on the extensive form. On the modular family, whose subproblems are linear programs, the results are negative: the nested forms reach the optimum of t48_s10_b2c, but they are 450 to 900 times slower than BDS and about 250 times slower than the InvestmentBlock with Gurobi inside, and none of them closes s50_b10c in 1800 s. This is most likely because the dual works through thousands of small components where Gurobi solves one linear program per scenario. With the dual inside, the InvestmentBlock does not even complete its first evaluation on s50_b10c, and it takes 152 GiB, which we have not investigated yet. Also, BDS with a bundle master and Gurobi in the subproblems reaches the optimum of s50_b10c in 76.4 s, against the 8.2 s of the MILP master.

Choice of the form and of the Solver

From the measurements above one obtains, for each class of networks, the form and the Solver to use; the ranges are those of the tables.

class of network	form	Solver	what it gives	against Gurobi (in PyPSA or in SMS++)
committable thermal units, a few scenarios, each scenario a MILP Gurobi closes	TSSB	LD of the scenarios (TSSBSCfg-LD-IP.txt)	a bound equal to the optimum to 6 digits	1.03 to 3.1 times faster and 1.6 to 3.9 times less memory than MILP SMS++ to the optimum
committable thermal units, long horizons or scenarios Gurobi does not close	TSSB	LD+LD or recursive LD for the bound, PPH (TSSBSCfg-PPH.txt) for a solution	the bound of the convex hull of each unit, and a feasible solution with its gap	a solution 0.0005% above the optimum in 151 s, where MILP SMS++ takes 2287 s and PyPSA 2840 s (u80_t168_s3); Gurobi alone reaches the same bound at its root, sooner on 3 instances out of 5
integer first stage (modules with p_nom_mod), a linear program per scenario	TSSB, Benders form (-k)	BDS (TSSBSCfg-BDS.txt), threads in intMaxThread	the optimum	at 1 thread 52 to 540 s, where PyPSA at 1 thread does not close in 3600 s; at 32 threads 4.6 to 36.7 s, 7 to more than 98 times faster than PyPSA with all the threads
continuous first stage, a linear program per scenario	TSSB or MSSB, Benders form	BDS	the optimum	1.6 to 1.9 times faster than LP SMS++ (Gurobi) on the TSSB, and PIPS-IPM++ about 3 times slower than Gurobi
continuous design, one scenario	UCBlock with the design in the units	LP SMS++ (Gurobi)	the optimum	2.5 to 5.6 times faster than PyPSA; the InvestmentBlock over the UCBlock is 1.1 to 2.2 times faster than PyPSA, BDS 2.3 to 12 times slower than PyPSA

Over a TSSB or an MSSB the InvestmentBlock is not in the table, since it solves the stochastic Block in full at each evaluation; its place is rather the deterministic case of the last row, and possibly the cases in which the inner Block is not a linear program (e.g., an SDDPBlock). Similarly, the nested forms may be a choice when the subproblem is not a linear program that Gurobi closes (e.g., a unit commitment), whereas on linear subproblems they are usually much slower than the plain forms.

The networks in the batteries of SMS++

Small instances of the three families are in the batteries of tests (branch develop), each in the battery of the module that reads its form, with the reference value of PyPSA; the data are in the nc4 archives of 2026-09-26 of UCBlock and InvestmentBlock in the Package Registry, which their CMake downloads. Each battery cross-checks all the Solvers of one BlockSolverConfig on the same instance, and a Solver that gives a bound is declared as such (-R of the tester for a relaxation, -E inf for a heuristic, whose interval has to contain the optimum).

battery	instance and form	Solver cross-checked
TwoStageStochasticBlock/batches/batch-pypsa	thermal u10_t24_s3, TSSB	MILP, PrimalProximalHeur, recursive LD (BSPar-2S-LD-IP.txt); on the Benders form, BDS whose subproblems are solved by the recursive LD, with its feasible solution recovered at the design of the master, with a MILP master (BSPar-BDS-2S-LD.txt) and with a bundle master (BSPar-BDS-2S-CVX-LD.txt)
TwoStageStochasticBlock/batches/batch-pypsa-modular	modular t48_s10_b2 and t168_s5_b1 (integer design), t48_s10_b2c (continuous), as TSSB and as MSSB	MILP and BDS on the Benders form (BSPar-BDS-2S-IP.txt); on the continuous ones also BDS whose subproblems are solved by the recursive LD, with a MILP master (BSPar-BDS-2S-LD.txt) and with a bundle master (BSPar-BDS-2S-CVX-LD.txt), and the recursive LD on the extensive form (BSPar-2S-LDrec-IP.txt)
MultiStageStochasticBlock/batches/batch-pypsa	modular t48_s10_b2c, MSSB	MILP
InvestmentBlock/batches/batch-stochastic	modular t48_s10_b2c, InvestmentBlock over the TSSB and over the MSSB	InvestmentBlock, with the inner Block solved by Gurobi and by the recursive LD (InnerBCfg-LD.txt)
InvestmentBlock/batches/batch-pypsa	modular t48 b2c, one scenario, InvestmentBlock over the UCBlock	InvestmentBlock, with the inner UCBlock solved by Gurobi and by the recursive LD (InnerBCfg-LD.txt)
UCBlock/batches/batch-pypsa	the same network, design in the UCBlock	MILP, LD, PrimalProximalHeur, 3 formulations of the units

Each configuration has one Lagrangian dual, the recursive one, which relaxes the most at once. Several LagrangianDualSolvers attached to the same Block (a PrimalProximalHeur included) give each the bound it gives when it is the only one, in any order.

Limits and open questions

Our measurements are limited to the Solvers that apply to these networks as pypsa2smspp writes them. BDS with an integer recourse, e.g., a TSSB whose scenarios are unit commitments, gives a bound instead of the optimum, since its subproblems are solved as their continuous relaxation: on the thermal family it gives the LP value (1.04614e7 on u80_t48_s3, 5.12561e7 on u80_t336_s3, in 1.1 and 12.8 s), which is weaker than the Lagrangian bound. Its relaxed binaries have bounds and no OneVarConstraint; hence BDSSCfg.txt sets intThrowReducedCostException 0, without which the subproblem throws. A Lagrangian dual over an MSSB does not run yet, since relaxing the non-anticipativity of the outer stage gives each LagBFunction a TSSB as its inner Block, and a LagBFunction requires an Objective of that Block, which TSSB and MSSB do not have. Other Solvers are left out because they answer a different question or give nothing new here: (i) ScenarioReductionSolver solves the problem on K representative scenarios, and its measure would be the error against the full set; (ii) SDDP reads an SDDPBlock, which the converter does not write yet; (iii) BranchAndXSolver works on knapsack problems only, for now; (iv) FrankWolfeSolver is on these problems the Dantzig-Wolfe side of the Lagrangian dual, and hence it gives the same bound. Also, intDoEasy, which keeps the network as an "easy" component of the inner master, matters on networks with many buses, while the thermal family has one. Of course, CPLEX, SCIP and HiGHS can replace Gurobi wherever it appears in the templates, by changing one line of the configuration; we used Gurobi because PyPSA-Eur uses it, so that the two sides are comparable.

Two questions remain open. The first is the memory of the InvestmentBlock with the recursive LD inside on s50_b10c (152 GiB before its first evaluation), which we have not explained yet and intend to investigate. The second is whether the nested forms pay on networks whose subproblems are unit commitments of realistic size, which Gurobi does not close. The thermal family here is small enough for Gurobi, and hence it only shows that they give the right bound and solution.

Reproduction

<cfg> is the folder pysmspp/data/configs/TSSBlock/ of this PR, with MPBCfg_grb.txt copied over MPBCfg.txt (the master of the measurements is Gurobi, while the default of the templates is HiGHS). Apart from gen_thermal_tssb.py and emit_thermal_tssb.py, which are in test/references of pypsa2smspp, the generators are in scripts/smspp_instances/references of SPSUnipi/pypsa-eur-instances#16.

# thermal (the formulation of the ThermalUnitBlock is the first value of TUBCfg.txt in <cfg>)
python gen_thermal_tssb.py --units 80 --snapshots 96 --scenarios 5 --buses 1
python emit_thermal_tssb.py tuc_u80_t96_s5_b1
python extensive_thermal_tssb.py tuc_u80_t96_s5_b1 gurobi MIPGap=1e-6     # or: relax
smspp_tssb_solver -v2 -c <cfg> -B InnerBCfg.txt -S <template> smspp_tuc_u80_t96_s5_b1.nc4

# modular, integer and continuous (c), TSSB and MSSB
python gen_modular_tssb.py --snapshots 168 --scenarios 200 --buses 10
python emit_modular_tssb.py mod_t168_s200_b10 ucblock
python emit_modular_tssb.py mod_t168_s200_b10c mssb_ucblock 10
python emit_modular_tssb.py mod_t168_s200_b10c mssb_investment_outside 10
python solve_modular_tssb.py mod_t168_s200_b10 gurobi MIPGap=1e-6
smspp_tssb_solver -v2 -c <cfg> -S TSSBSCfg-IP.txt smspp_mod_t168_s200_b10_ucblock.nc4
smspp_tssb_solver -v2 -k -c <cfg> -S TSSBSCfg-BDS.txt smspp_mod_t168_s200_b10_ucblock.nc4

# nested forms: BDS whose subproblems are solved by the recursive LD, with a MILP or a bundle master, and the InvestmentBlock with it inside
smspp_tssb_solver -v2 -k -c <cfg> -B InnerBCfg.txt -S TSSBSCfg-BDS-LD.txt smspp_mod_t168_s50_b10c_ucblock.nc4
smspp_tssb_solver -v2 -k -c <cfg> -B InnerBCfg.txt -S TSSBSCfg-BDS-LD.txt smspp_tuc_u80_t48_s3_b1.nc4
smspp_tssb_solver -v2 -k -c <cfg> -B InnerBCfg.txt -S TSSBSCfg-BDS-CVX-LD.txt smspp_tuc_u80_t48_s3_b1.nc4
smspp_tssb_solver -v2 -k -c <cfg> -B InnerBCfg.txt -S TSSBSCfg-BDS-CVX-LD.txt smspp_mod_t168_s50_b10c_ucblock.nc4
smspp_investmentblock_solver -v2 -c <cfg> -B InnerBCfg-LD.txt -S TSSBSCfg-IB.txt smspp_mod_t168_s50_b10c_investment_outside.nc4
smspp_investmentblock_solver -v2 -c <cfg> -B InnerBCfg.txt -S TSSBSCfg-IB.txt smspp_mod_t168_s200_b10c_o10_mssb_investment_outside.nc4

# deterministic: one scenario, the design in the UCBlock or in an InvestmentBlock over it
python gen_modular_tssb.py --snapshots 2016 --scenarios 1 --buses 10
python emit_modular_tssb.py mod_t2016_s1_b10c det_ucblock
python emit_modular_tssb.py mod_t2016_s1_b10c det_investment
python emit_modular_tssb.py mod_t2016_s1_b10c ucblock
smspp_ucblock_solver -v2 -c <cfg> -S BSCfg1-IP.txt smspp_mod_t2016_s1_b10c_det_ucblock.nc4
smspp_investmentblock_solver -v2 -c <cfg> -B InnerBCfg.txt -S TSSBSCfg-IB.txt smspp_mod_t2016_s1_b10c_det_investment.nc4
smspp_tssb_solver -v2 -k -c <cfg> -S TSSBSCfg-BDS.txt smspp_mod_t2016_s1_b10c_ucblock.nc4

# PIPS-IPM++, one process per group of blocks, one thread each
OMP_NUM_THREADS=1 mpirun -np 10 -x OMP_NUM_THREADS smspp_tssb_solver -v2 [-k] -c <cfg> -S TSSBSCfg-PIPS.txt smspp_mod_t168_s200_b10c_ucblock.nc4

[davide-f]

@dmeoli this will be needed in the next release correct? it is not yet in develop or yes?

[dmeoli]

Right: they need the BundleSolver 2.0, which is not in the develop of SMS++ yet, and that is why this one is a draft. The only red job is the one that installs SMS++ from conda-forge, which stops on Invalid int parameter name intMPStbl, i.e., the released BundleSolver is still the 1.0. It goes green with the release that carries the 2.0, and I keep the PR in draft until then.

[AlessandroPampado99]

@dmeoli CI does not pass, could you trigger manually Tests / Test package with latest SMSpp (pull_request) from your branch?

[AlessandroPampado99]

"pypsa2smspp writes the netCDF file in one of several forms, selected by a few of its knobs. capacity_expansion_ucblock puts the design in the units (the form that BDS and the MILP take), investment_outside states it once in an InvestmentBlock above them, and design_cost_outside leaves it in the units but takes its cost out."

design_cost_outside is depcrecated and should not be used SPSUnipi/pypsa2smspp#65

Modular family, continuous design, TSSB and MSSB section has MILP SMS++ in table. Is it a mistake? It should be an LP

[dmeoli]

@AlessandroPampado99 thanks, both points are fixed in the description: design_cost_outside is gone, and the column of the continuous instances is now "LP SMS++ (Gurobi)", as in the deterministic table, since those are linear programs.

About the CI: the job with the latest SMS++ runs only on schedule or workflow_dispatch, so I triggered it on the branch of my fork: https://github.com/dmeoli/pySMSpp/actions/runs/36426995339. It was red on main as well (the scheduled runs fail since 22/09), for reasons that are not in the templates:

• macOS: the CMake of four modules extracted their data with tar --warning=no-unknown-keyword, an option of GNU tar that the tar of macOS does not have; they now use cmake -E tar, on the develop of SMS++.
• Ubuntu: add_hub_to_ucblock counted the HydroUnitBlock as one electrical generator, while its 3 arcs are 3 generators, which UCBlock now checks when it reads the Block; and test_optimize_svmblock_formulations used GRBMILPSolver, which the CI does not build. Both are fixed in this PR (a7cc5c3): the conftest counts the arcs, and the SVM test skips a Solver that is not in the factory, as it already did for LIBSVM.
• Windows: vcpkg stops on liblemon, because lemon.cs.elte.hu now serves a regenerated archive of the same revision with a different hash; I opened [liblemon] Update archive hash microsoft/vcpkg#54163, and once it is merged the baseline of SMS++ moves to it.

The job with SMS++ from conda-forge stays red until the release that carries BundleSolver 2.0.

[dmeoli]

The run with the latest SMS++ is green on macOS and Ubuntu now (https://github.com/dmeoli/pySMSpp/actions/runs/36426995339); Windows still stops on the hash of liblemon, and goes green once microsoft/vcpkg#54163 is merged.