Changes for the BundleSolver 2.0

[dmeoli]

The BundleSolver 2.0 it is written for is now in the develop of SMS++ (#59, on which it was stacked, is merged).

A modular extendable asset. An asset with p_nom_mod > 0 is built by PyPSA in an integer number of modules, and the UCBlock form now writes the design of such an IntermittentUnitBlock as that number, reading the capacity back as the number of modules times their size; the InvestmentBlock form does the same, with the cost and the bounds of the design per module and the asset marked Integer, and the BundleSolver of InvestmentBlock/BSPar-int.txt keeping its design integer. It needs UCBlock 87edf11f, and for the InvestmentBlock form InvestmentBlock e43b146, BundleSolver e1c624f and SMS++ 8f41958. The case is mod_1n_2c_2g, whose continuous optimum builds 2.36 modules and whose modular one builds 3, so that rounding would give the wrong answer.

The susceptance of the lines. The lines were written with LineSusceptance 0, i.e., as a transport model without Kirchhoff's voltage law, which went unnoticed on networks without cycles of AC lines: on the 20-node networks of power-IT-noUC (30 lines, 12 cycles) SMS++ was 0.08% below PyPSA. They now get 1/x_pu_eff (and 0 where the reactance is 0, i.e., a DC line), and SMS++ gives the value of PyPSA to about 1e-5.

The batteries of PyPSA-Eur. merge_links never merged their charger and discharger: the discharger has no capital cost, and the extendable assets without a capital cost are fixed at a large capacity before the merge, which only merges an extendable pair. The links that the merge absorbs are now left extendable; merge_links also accepts a preset or a list of presets, as its docstring says, and the merged links are read back also in the stochastic case. @AlessandroPampado99, the merged link has p_min_pu = -efficiency of the discharger, so that it discharges at most efficiency * p_nom on the electric side, while the constraint of PyPSA-Eur (p_nom of the charger = efficiency * p_nom of the discharger) lets it discharge p_nom there: is the factor intended? On power-IT-noUC the battery constraint is not binding, so the values are the same either way.

The bounds of the extendable assets. Gone from this PR, as agreed below: an asset PyPSA gives no p_nom_max keeps an infinite MaxCapacityDesign, which the conversion already writes, and the vector carries a finite entry only where PyPSA has one.

The templates of the tools for the BundleSolver 2.0 are in SPSUnipi/pySMSpp#116, and the generators of the SMS++ instances in SPSUnipi/pypsa-eur-instances#16.

[dmeoli]

The bound this PR gives an extendable asset with no bound of its own is the right answer measured from the Lagrangian side as well, and there is one sentinel left that bound_extendable_assets() does not reach.

What the sentinels do to the Lagrangian dual. The four sector-coupled TSSB instances of batch-pypsa were written when the conversion replaced an infinite e_nom_max with a finite number, and carried ConverterMaxCapacityDesign = 1e10 on nine converters and BatteryMaxCapacityDesign = 1e9 on nine batteries. On the one whose stochastic datum is the demand, the LagrangianDualSolver ends in error while a :MILPSolver solves the deterministic equivalent in 0.7 s with the right value. With the log on, the bundle dies at iteration 136 of 160, i.e., a step away from the optimum, with Fi = 4.8187786914e+09 against a reference of 4.8226936e+09: the master stops being solvable, the bundle deletes its own items one at a time for sixty-odd rounds and ends with Bundle::FormD: unrecoverable MP failure, the multipliers run away and the value it reports is -1.07e+22.

The cause is the cap, and it reaches the master through the subgradients. In the Lagrangian relaxation of the non-anticipativity rows, the design of a converter that the multiplier makes profitable goes straight to its bound, so the residual of the relaxed row carries the bound itself: 262 of the 311 linearizations of that solve have an infinity norm of exactly 1.00e+10, while the ones on the physical scale are at 1e4 to 1e5. The multipliers are ordinary, 1.3e+04; it is their product with those subgradients that reaches 1e+12 while the function is worth 1e9, so that a row of the master of the bundle is a difference of terms of 1e12 that has to come out at 1e8, and Gurobi answers GRB_NUMERIC near the optimum.

Two counter-proofs, on the same tree and with the master configured exactly as it was when it failed:

• rewriting in that instance the caps >= 1e8 to 1e6, and nothing else, it solves in 21.4 s with the exact value instead of dying at 37.6 s;
• regenerating the four instances with the conversion as it stands, where an uncapped e_nom_max stays infinite, batch-pypsa passes its eight runs with no KO. The objective value of each of the four is the one it was to the tenth digit, the caps having never been binding; the Lagrangian dual is slower on them, 62.5 s instead of 22.9 on the complete one, which is the price of the instance being the one the network describes.

The sentinel that is left. preprocess_zero_capital_cost_extendable_generators() and preprocess_zero_capital_cost_extendable_lines_links() freeze an extendable asset whose capital cost is zero at fixed_capacity = 1e9, and that lands as MinCapacityDesign = MaxCapacityDesign = 1e9 in the instances this branch writes as well. Where I looked it does no harm, the design being fixed and the residual of its row therefore identically zero, but it is a billion MW among the coefficients of every row that asset appears in, and it is the same number in the same place as the ones this PR removes. The reason it is there, that an asset which expands for free and without a finite bound is ill-posed, is met by the physical bound just as well, so fixed_capacity could read it off SMSPP_DESIGN_BOUNDS instead of being a constant.

[AlessandroPampado99]

The bounds of the extendable assets of the instances. The instances were written with the extendable assets capped at 1e7 (1e8 on the links, which become the converters of the stores). The caps are not a property of the networks: they were there because with an unbounded design some Lagrangian subproblem is unbounded, and the bundle in the released SMS++ cannot remove from its master problem the feasibility linearizations such a subproblem produces. They are also what kills the master of the BundleSolver 2.0: the design is bang-bang, the value of a component becomes the bound times the investment cost, and the master is left with coefficients of 1e9 against a quadratic term of 10, on which Gurobi turns to quad precision and ends in a numerical error. bound_extendable_assets(n, mode, safety) in network_correction.py gives an extendable asset with no bound one the demand of the network pays for, and the two generators call it, reading SMSPP_DESIGN_BOUNDS:

Well, this assumption should be wrong... We are helping the solution but some optima issues could arise with physical

[AlessandroPampado99]

@dmeoli Could you test PyPSA–SMS++ for the new p_nom_mod conversion?

For example, with 50 MW modules, p_nom_min = 70 MW and p_nom_max = 180 MW, the design should allow only 2 or 3 modules. The test should check the investment cost per module, the dispatch limit including p_max_pu, and the conversion of the solution back to MW. Ideally, choose a case where the continuous relaxation would build a non-integer number of modules, so that the test actually checks integrality, and compare both the capacity and objective with PyPSA.

Also, since emit_modular_tssb.py explicitly relaxes modularity for investment_outside, that form should be checked against the corresponding continuous PyPSA model rather than the integer reference.

[dmeoli]

Could you test PyPSA–SMS++ for the new p_nom_mod conversion?

Done in 1f76a36, test/test_modular.py. The network has one bus, a flat load of 100 MW, a gas unit at 100 €/MWh and an extendable wind generator with modules of 50 MW, p_nom_min = 70, p_nom_max = 180, capital_cost = 800 and a p_max_pu that varies over 24 snapshots, so that the design may take 2 or 3 modules. The continuous model builds 118.15 MW, i.e., 2.36 modules, while the modular one builds 150 MW, i.e., 3 modules, hence rounding the continuous optimum would give 2 modules and the wrong answer; the first test checks that this is so.

The second test reads back from the netCDF file, with no SMS++, what the conversion writes: MaxCapacityDesign = -3 (at most 3 modules, the sign making the design integer), MinCapacityDesign = 2, InvestmentCost = 40000 (the cost of a module) and MaxPower = 50 * p_max_pu. The third one solves the UCBlock form and compares both the objective and the capacity that comes back, the number of modules times their size, with PyPSA:

form               SMS++                          PyPSA
UCBlock            199968.696135, 150 MW          199968.696135, 150 MW (modular)
InvestmentBlock    194318.95788, 118.146 MW       194318.95788, 118.146 MW (continuous)

The test found a defect in UCBlock and not in the conversion: IntermittentUnitBlock::check_data_consistency() refused a MinCapacityDesign above 1 when MaxCapacityDesign is negative, as if the design were binary, while the design is generated as an integer in {MinCapacityDesign, ..., |MaxCapacityDesign|}. It is fixed in UCBlock develop (87edf11f); since the released SMS++ still has the old check, the UCBlock comparison is skipped where the solver on PATH refuses the file, in the same way as the pollutant budgets.

Also, since emit_modular_tssb.py explicitly relaxes modularity for investment_outside, that form should be checked against the corresponding continuous PyPSA model rather than the integer reference.

This is the fourth test: the InvestmentBlock form of the same modular network writes a continuous capacity in [70, 180] at 800 €/MW, and its optimum is compared with the continuous PyPSA model (the second row of the table). The test uses the default configuration, i.e., the one of the released SMS++; the values above come from a build of the BundleSolver 2.0 with BSPar_2.0.txt, and I have not run the test against the released SMS++.

[dmeoli]

Well, this assumption should be wrong... We are helping the solution but some optima issues could arise with physical

You are right, and it is not only a risk: the physical bound cuts the optimum of one of the test networks. The sentence in the docstring and in the description, that no asset can be used beyond what the demand absorbs, confuses the capacity of an asset with what it produces: a generator produces its capacity times p_max_pu, so an intermittent one needs more than the peak of the load just to cover it, and more again when a storage unit takes what it produces in excess.

I solved with PyPSA every test network with an extendable asset twice, once with no bound (mode="none") and once with the physical one, and looked for the assets whose unbounded optimum lies beyond the bound:

network                        unbounded      physical       asset at the bound
3n_3c_2gext_1lext_1mlext       4.862560e+06   1.090674e+07   wind: 47674 MW against 37336 MW
2n_1c_1gext_1bext_2lext        2.991550e+10   2.991550e+10   H2 Fuel Cell at the bound, same objective
the other 9                    same           same           none

On 3n_3c_2gext_1lext_1mlext the largest p_max_pu of the wind generator is 0.541, so it needs 1.85 times the peak to cover it alone, and the optimum builds 2.5 times the peak because the batteries store the excess. The load of the test networks is drawn at random, so the numbers change from one run to the next, but the bound was active on all 3 draws I tried. The instances of the batches show the same: when the physical bound replaced the 1e7 caps, the reference of smspp_3n_3c_2gext_1lext_1mlext went from 4.682891301e+06 to 1.495190064e+07, i.e., the reference is now the optimum of a network that has been cut. SMS++ and PyPSA still solve the same network there, so the comparison is sound, but it is not the network that was written.

No bound taken from the demand alone can be proved not to cut, and the one that can, i.e., capital_cost * p_nom_max at most the cost of a feasible point, is useless here: with the slack units that cost reaches 7e11, and the bounds that come out (1e9 to 1e10 MW) are the same caps that broke the master.

What I propose is to check the bound after solving, instead of trusting it: the generators solve every instance with PyPSA anyway to write its reference, so an asset whose design ends at its bound gets the bound doubled and the instance is solved again, until no design is at its bound. For a linear network this is exact, since a bound that is not active at the optimum of a convex problem has a zero multiplier and does not cut it; for a network with committable units it is the same check, not a proof. The docstring and the description then say what the bound actually is, a starting point verified on the solution. If you agree, I will push it here and regenerate the instances whose bound was active.

[AlessandroPampado99]

Well, this assumption should be wrong... We are helping the solution but some optima issues could arise with physical

You are right, and it is not only a risk: the physical bound cuts the optimum of one of the test networks. The sentence in the docstring and in the description, that no asset can be used beyond what the demand absorbs, confuses the capacity of an asset with what it produces: a generator produces its capacity times p_max_pu, so an intermittent one needs more than the peak of the load just to cover it, and more again when a storage unit takes what it produces in excess.

I solved with PyPSA every test network with an extendable asset twice, once with no bound (mode="none") and once with the physical one, and looked for the assets whose unbounded optimum lies beyond the bound:

network                        unbounded      physical       asset at the bound
3n_3c_2gext_1lext_1mlext       4.862560e+06   1.090674e+07   wind: 47674 MW against 37336 MW
2n_1c_1gext_1bext_2lext        2.991550e+10   2.991550e+10   H2 Fuel Cell at the bound, same objective
the other 9                    same           same           none

On 3n_3c_2gext_1lext_1mlext the largest p_max_pu of the wind generator is 0.541, so it needs 1.85 times the peak to cover it alone, and the optimum builds 2.5 times the peak because the batteries store the excess. The load of the test networks is drawn at random, so the numbers change from one run to the next, but the bound was active on all 3 draws I tried. The instances of the batches show the same: when the physical bound replaced the 1e7 caps, the reference of smspp_3n_3c_2gext_1lext_1mlext went from 4.682891301e+06 to 1.495190064e+07, i.e., the reference is now the optimum of a network that has been cut. SMS++ and PyPSA still solve the same network there, so the comparison is sound, but it is not the network that was written.

No bound taken from the demand alone can be proved not to cut, and the one that can, i.e., capital_cost * p_nom_max at most the cost of a feasible point, is useless here: with the slack units that cost reaches 7e11, and the bounds that come out (1e9 to 1e10 MW) are the same caps that broke the master.

What I propose is to check the bound after solving, instead of trusting it: the generators solve every instance with PyPSA anyway to write its reference, so an asset whose design ends at its bound gets the bound doubled and the instance is solved again, until no design is at its bound. For a linear network this is exact, since a bound that is not active at the optimum of a convex problem has a zero multiplier and does not cut it; for a network with committable units it is the same check, not a proof. The docstring and the description then say what the bound actually is, a starting point verified on the solution. If you agree, I will push it here and regenerate the instances whose bound was active.

@davide-f

[AlessandroPampado99]

Configs should not stay here as said in #58

[dmeoli]

Done in d310a94: the configurations are no longer in pypsa2smspp. Those of the BundleSolver 2.0 for an InvestmentBlock are in config/smspp/InvestmentBlock-BS2 of SPSUnipi/pypsa-eur-instances#16, together with the generators that use them, and the tool templates in SPSUnipi/pySMSpp#116. The PR now carries only the conversion of the modular assets and the bounds of the extendable ones (see the description).

[davide-f]

Well, this assumption should be wrong... We are helping the solution but some optima issues could arise with physical

You are right, and it is not only a risk: the physical bound cuts the optimum of one of the test networks. The sentence in the docstring and in the description, that no asset can be used beyond what the demand absorbs, confuses the capacity of an asset with what it produces: a generator produces its capacity times p_max_pu, so an intermittent one needs more than the peak of the load just to cover it, and more again when a storage unit takes what it produces in excess.
I solved with PyPSA every test network with an extendable asset twice, once with no bound (mode="none") and once with the physical one, and looked for the assets whose unbounded optimum lies beyond the bound:

network                        unbounded      physical       asset at the bound
3n_3c_2gext_1lext_1mlext       4.862560e+06   1.090674e+07   wind: 47674 MW against 37336 MW
2n_1c_1gext_1bext_2lext        2.991550e+10   2.991550e+10   H2 Fuel Cell at the bound, same objective
the other 9                    same           same           none

On 3n_3c_2gext_1lext_1mlext the largest p_max_pu of the wind generator is 0.541, so it needs 1.85 times the peak to cover it alone, and the optimum builds 2.5 times the peak because the batteries store the excess. The load of the test networks is drawn at random, so the numbers change from one run to the next, but the bound was active on all 3 draws I tried. The instances of the batches show the same: when the physical bound replaced the 1e7 caps, the reference of smspp_3n_3c_2gext_1lext_1mlext went from 4.682891301e+06 to 1.495190064e+07, i.e., the reference is now the optimum of a network that has been cut. SMS++ and PyPSA still solve the same network there, so the comparison is sound, but it is not the network that was written.
No bound taken from the demand alone can be proved not to cut, and the one that can, i.e., capital_cost * p_nom_max at most the cost of a feasible point, is useless here: with the slack units that cost reaches 7e11, and the bounds that come out (1e9 to 1e10 MW) are the same caps that broke the master.
What I propose is to check the bound after solving, instead of trusting it: the generators solve every instance with PyPSA anyway to write its reference, so an asset whose design ends at its bound gets the bound doubled and the instance is solved again, until no design is at its bound. For a linear network this is exact, since a bound that is not active at the optimum of a convex problem has a zero multiplier and does not cut it; for a network with committable units it is the same check, not a proof. The docstring and the description then say what the bound actually is, a starting point verified on the solution. If you agree, I will push it here and regenerate the instances whose bound was active.

@davide-f

An initial comment: @dmeoli / @AlessandroPampado99 let's keep the description always coincise and readable by a human. Details can be added later but it is very hard to review the PRs otherwise.

If I understand correctly, here are are exploring the inclusion of bounds to investment capacities.

Preliminary question: without bounds the lagrangian cannot work? I guess so, I'm asking for confirmation

The heuristic must be light. For sure, initial pypsa optimizations to calibrate the bounds are a no-go.
I'd definitely prefer a heuristic like "1e8/1e9" and add a warning in such cases. We can have other heuristics, but they must be very light both in terms of methodology and code to be maintainable and understandable.

[dmeoli]

Short answers, and I am dropping the calibration runs.

Does the Lagrangian need a bound? No. What hurts is not the missing bound, it is a bound written as a big finite number, because the design of an asset the multiplier makes attractive goes straight to it and the bound lands in the subgradient: on network_small_fewsectors__snap100__scen3_demand, 262 of the 311 linearizations had an infinity norm of exactly 1e10, the master row b + <g,L> cancelled three orders of magnitude, Gurobi answered GRB_NUMERIC near the optimum and the bundle emptied itself and died. The same run with those caps rewritten to 1e6 closes in 21 s at the exact value. The other way round, on the small test networks, setting only the sentinels (1e7 to 1e10) to infinity changes no optimum at all.

So 1e8/1e9 is the shape that hurt, and no bound at all is safer than a big one. What I propose, and what is also the least code, is to write no MaxCapacityDesign where PyPSA gives no p_nom_max, and keep the bound only where PyPSA has a real one. No PyPSA run to calibrate anything, no doubling loop, and the warning goes where an asset has no bound at all.

Worth adding that the BundleSolver 2.0 this PR targets scales the rows of the master by default and scales the cuts inside each hard component, which is the cure for rows that differ by orders of magnitude, so a large sentinel is much less dangerous there than it was on the 1.0. That makes the choice a matter of what is simple rather than of what survives.

Description shortened.

[davide-f]

I like the proposal: do not write the max where pypsa does not specify it.
Note that some assets may define it and some may not, so the vector maxcapacity is needed but some entries may be infinite.
We need to support this case.

[dmeoli]

Already supported, and now that is all the PR does about bounds: MaxCapacityDesign is written per asset, finite where PyPSA has a p_nom_max and infinite where it has none, so a mixed Block is fine, and an infinite bound is safe here because a design with a cost cannot run away (an extendable asset with zero capital cost is made non-extendable upstream).

Dropped in a6e9c4f: bound_extendable_assets, verify_extendable_bounds and their helpers, 174 lines, so network_correction.py is back to what main has. I will take them out of the generators of SPSUnipi/pypsa-eur-instances#16 too.

[davide-f]

@dmeoli have you added support for modular capacity expansion in InvestmentBlock?

[AlessandroPampado99]

docs/readthedocs.org:pypsa2smspp — Read the Docs build failed!
This fails, then we can merge @dmeoli (after Davide's comments)

[dmeoli]

@davide-f yes, in both forms (this comment is edited: at first the InvestmentBlock form relaxed the modules).

With capacity_expansion_ucblock: true, an extendable generator with p_nom_mod > 0 (i.e., an IntermittentUnitBlock; the batteries are not covered) gets an integer design, i.e., the number of its modules (MaxCapacityDesign = -floor(p_nom_max / p_nom_mod), the negative sign being how UCBlock declares an integer design), and the capacity is read back as the number of modules times their size. This is a6e9c4f, tested by the case mod_1n_2c_2g, whose continuous optimum builds 2.36 modules and whose modular one builds 3.

In the InvestmentBlock form the design of such an asset is now the number of its modules as well: the InvestmentBlock reads a netCDF variable Integer (InvestmentBlock e43b146), and its BundleSolver keeps an integer design integer with intIntVars (BundleSolver e1c624f, SMS++ 8f41958), i.e., the stabilized cutting-plane method of van Ackooij, Frangioni and de Oliveira, whose master is a mixed-integer program. pypsa2smspp writes the cost and the bounds per module, marks the asset Integer, and picks the template InvestmentBlock/BSPar-int.txt of SPSUnipi/pySMSpp#116; on mod_1n_2c_2g the InvestmentBlock gives 199968.69614 as PyPSA (relative difference 2.4e-11), with 3 modules of wind.

@AlessandroPampado99 the Read the Docs failure is not in the docs: the build of 24/09 started while the branch was being rewritten and stopped on fatal: reference is not a tree: a6e9c4f. I am restarting it.

[AlessandroPampado99]

@AlessandroPampado99 the Read the Docs failure is not in the docs: the build of 24/09 started while the branch was being rewritten and stopped on fatal: reference is not a tree: a6e9c4f. I am restarting it.

@dmeoli If it's ok I'll merge it :)