feat: add m.dualize()

[bobbyxng]

Changes proposed in this Pull Request

• Adds Model.dualize(), a method that constructs the LP dual of a linopy model, and Model.bounds_to_constraints(), a preprocessing step that converts variable bounds to explicit constraints so they are correctly reflected in the dual.
• The dual is constructed following standard LP duality theory. For a primal minimization problem, the dual is a maximization problem with one dual variable per primal constraint. Variable bounds are converted to explicit constraints before dualization via bounds_to_constraints(), so that they appear in the constraint matrix and are correctly reflected in the dual.
• Signs follows linopy's dual convention, allowing direct comparison between m2.variables[con_name].solution and m1.constraints[con_name].dual without sign adjustments.
• Building on PR feat: add m.copy() method to create deep copy of model #623 which adds m.copy()

The motivation behind this was an automated LP dualizer (originally developed for adaptive robust optimization for energy system planning) which requires constructing an independent dual model without modifying the original. While similar implementations exist in other modelling frameworks, notably JuMP's Dualization.jl https://jump.dev/JuMP.jl/stable/packages/Dualization/, this feature was still missing in linopy.

Checklist

• Code changes are sufficiently documented; i.e. new functions contain docstrings and further explanations may be given in doc.
• Unit tests for new features were added (if applicable).
• A note for the release notes doc/release_notes.rst of the upcoming release is included.
• I consent to the release of this PR's code under the MIT license.

[bobbyxng]

Here's a successful dualization tested on a 50-node PyPSA network
Notice the small deviations between primal and dual solution, given the nature/dual degeneracy of such energy system problems.

Primal objective: 95047351841.07184
Dual objective: 95047351841.07161
Abs. gap: 0.0002288818359375
Rel. gap: 2.4080821980207514e-15

Solution Generator-p

Solution StorageUnit-state_of_charge

[brynpickering]

@bobbyxng I can't comment on the actual dualisation so I'm focussing on cleanliness and efficiency of code. I've opened another PR (#629) in which I've made some suggestions (I started writing them in comments but they became a bit too verbose).

I feel like it shouldn't be necessary to have your lookup dictionary or rule and that you should be able to store your coefficients your mapping from constraint name to dual variables. I haven't investigated it in detail but the linear expression rule is probably reasonably slow to build compared to a vectorised array operation. I tried with a dummy model with 10000 timesteps and 2 active spatial nodes and it took 20 seconds to build. I can see dualisation therefore exploding on practical models, so more vectorisation is probably necessary.

You still also need:

• unit tests
• docs

BTW, for ease of review, it would be easier if you opened this PR with the copy branch as the base, rather than master. It will automatically revert to master when copy is merged.

[FabianHofmann]

@bobbyxng what is the status here? should someone review again?

[FBumann]

@bobbyxng Im not really getting what this PR does, as im not that deep into the LP world, but maybe you could use the new features from #634 to do Variables.relax()and maybe Variables.fix()?

EDIT: Probably not, just needed if a user wants to dualize a MILP. He would call Variables.relax() BEFORE Model.dualize()

[FBumann]

Quick review:

1. Could we move bounds_to_constraints() into dual.py to not expose it to users directly?
2. We should rename bounds_to_constraints() to lift_bounds_to_constraints() and clearly state that its converting bounds to constraints AND removing the bounds from the variables. We might add an unlift_*method later if we need it.

But I'd defer the unlift_* method.

[FBumann]

@bobbyxng Ignore the pypsa models workflow and possible revert your changes. It will be fixed in #522

[bobbyxng]

Thank you all @brynpickering @FBumann @FabianHofmann!

Sorry for the delay, I wanted to make sure the dualisation table was solid (also in my ARO workflow). I am confident that the rules all apply correctly. All PyPSA example networks yield the same results in primal and dual formulation (both on the objective and on the variable level, if dual degeneracy is accounted for). I have also added plenty of unit tests.

And thanks @brynpickering for the comment on (not) using the linear expression rule, I refactored everything to use numpy operations. Runtime dropped from ~6 minutes to ~2 seconds for a typical 50-node PyPSA network. But I think all this refactoring could use a careful look by someone more experienced with linopy internals, xarray and numpy operations. The dual.py code still relies on a few lookup helpers, happy to hear suggestions on how to avoid them :)

From my side, the PR is ready for review. :)