Some tricks to help numerical issues

[odow]

The bound projection thing

So you'll probably run into an issue that took me quite a while to debug.

If you take the solution of a node in the forward pass and fix it as the incoming state in the next node, that subproblem is not guaranteed to be feasible, even if, in theory, you have relatively complete recourse.

That's because the first solve may violate tolerances by a small amount, and then presolve in the next node deduces that the problem is infeasible. (See, e.g., https://jump.dev/JuMP.jl/stable/tutorials/getting_started/tolerances/.)

This will hit you here:

PlasmoAlgorithms.jl/lib/PlasmoBenders/src/solution.jl

Lines 172 to 175 in 6c9e489

	# Fix primal solutions
	for (j, var) in enumerate(comp_vars)
	JuMP.fix(var_copy_map[var], last_primals[j], force = true)
	end

There isn't a universal fix, but the most common cause of this is a solution violating the bounds on the state variables. A simple fix that catches almost every case is to project the state back onto its bounds:

https://github.com/odow/SDDP.jl/blob/aec55e2279c4abbd2a15ab1a430f6bb5a526976e/src/algorithm.jl#L182-L206

Resetting the optimizer

Another thing that is super(-ish) common is the following:

• You have a problem running along fine for many iterations
• Each iteration adds a cut that changes optimality but does NOT change feasibility (because your cost to go variable won't be bounded above)
• Suddenly, the solver says your problem is unbounded/infeasible/has numerical error.

The reason is that solvers may work from an existing basis. The sequence of updates can eventually lead to numerical issues, which, if the solver started from scratch, it wouldn't encounter. One solution is to use MOI.Utilities.reset_optimizer to force a reset if you detect a non-optimal status, and then bail completely only if that can't repair the solution.

Here's how we do it:

https://github.com/odow/SDDP.jl/blob/aec55e2279c4abbd2a15ab1a430f6bb5a526976e/src/algorithm.jl#L303-L400
https://github.com/odow/SDDP.jl/blob/aec55e2279c4abbd2a15ab1a430f6bb5a526976e/src/algorithm.jl#L481-L521

This feels like a massive hack, but it works surprisingly well.

[dlcole3]

@odow This is super helpful. Thank you! I have run into this where I get an infeasible model even when I am pretty confident it should be feasible, and I expect it is from these issues you note above. I will update this in a future PR.

[odow]

I have run into this where I get an infeasible model even when I am pretty confident it should be feasible

Yeah it's pretty common and super confusing to debug. @joaquimg suggested the resetting trick.