An allocation meta-optimization

[SouthEndMusic]

The idea

#3133 got me thinking about how the allocation optimization problems can be simplified to increase performance when needed.

Looking at that issue, we see that we can use all the equality constraints to reduce the problem to only a small number of unknowns, namely all the flows that are controlled by allocation. The remaining complexity comes from the bounds, where there are 2 categories:

• Bounds for goal programming (relatively small amount)
• Others: non-negative storage, storage upper bounds, etc.

So here is my idea:

• Get rid of all non-goal programming bounds (but keep the ones on the controlled flows);
• Investigate whether HiGHS is then able to reduce the problem to only the allocation controlled flows. If not, we can do our own problem transformation;
• Enforce the non-goalprogramming bounds with an isoutofdomain-type mechanism as it exists in the physical layer: check the bounds after the whole allocation solve, if any of the bounds are not met, run allocation again with a smaller timestep.

The performance gain here really depends on how cheap the new formulation is and how often allocation timesteps get rejected because of bound violations.

@jarsarasty I'm really interested in your take on this!

The problem transformation

The equations in #3133 yield the following linear system in $\Delta S$:

$$ \left[I_n - \frac{\Delta t}{2}M\left[I_m + CS_\text{p}^\top S_\text{p}\right]H\right] \Delta S = \Delta t \left[f + M\left[\frac{1}{2}S_\text{p}^\top S_\text{p}\hat{q} + \frac{1}{2}v + \frac{1}{2}CS_\text{p}^\top S_\text{p}v\right]\right] + \Delta t M\left[\left[I_m + \frac{1}{2}C \right]S_\text{c}^\top S_\text{c}Q + \left[ I_m + \frac{1}{2}C \left[S_\text{UD, in}^\top + S_\text{UD, out}^\top \text{diag}\left(\overline{\rho}\right) \right]\right]S_\text{UD, in}Q \right]. $$

Here:

$$ \begin{align} C &=& \frac{\partial q}{\partial c}\frac{\partial c}{\partial q} \\ H &=&\frac{\partial q}{\partial s} - \frac{\Delta t}{2}\frac{\partial q}{\partial i}S_\text{PID}\text{diag}\left(\frac{1}{\hat{a}}\right) + \frac{\partial q}{\partial c}\frac{\partial c}{\partial s} \\ v &=& \left[I_m - C\right]\hat{q} + \frac{\partial q}{\partial i}\left[\int_{t_n}^{t_n + \Delta t}r(\tau)\text{d}\tau - \Delta t S_\text{PID}\hat{h}\right] \\ \overline{\rho} &=& \frac{1}{\Delta t}\int_{t_n}^{t_n + \Delta t}\rho(\tau)\text{d}\tau \end{align} $$

This linear system is similar to the ones discussed #3113, #3125 and #3127, but note that the second term in the right hand side is expressed in terms of $S_\text{c}Q$ and $S_\text{UD, in}Q$, precisely the flows to be optimized by the allocation algorithm. Now we can solve this linear system to express the storage changes in terms of the to be optimized flows. This has 2 goals:

• Some storage changes are needed in the optimization problems; for level demand and for putting bounds so that the to be optimized flows do not cause negative storage in their directly upstream Basin;
• After the optimization we want to know $\Delta S$ for output and bounds checking (also the physical flows which can be directly computed with the expression for $S_\text{p}Q$ in Systematic linearization for allocation #3133).

For the route optimization objective we can also express the cost in terms of the to be optimized flows.

[jarsarasty]

Hi @SouthEndMusic,

Nice idea! Your reformulation reduces the number of variables and keeps the bounds as inequalities on the controlled flows, without iterative retries to restore violated bounds.

It may be worth having a look at the typical model size and time performance of the reformulation you proposed, against the current formulation.

If your reformulation is fast enough, we can avoid using an isoutofdomain-type mechanism, which could be risky for allocation: just dropping the non-goal-programming bounds would enlarge the feasible region, so the LP can choose flows that violate them (e.g. draining a basin below zero within a step); a smaller time step shrinks the violation but doesn't guarantee the same flows as the bounded problem, and retries add solver calls with no guaranteed speedup.

Separately, it's worth checking presolve/setup time vs. simplex time per solver call, since this overhead applies regardless of which formulation we use. If presolve/setup dominates, I would try disabling presolve (or tuning its options) and instead carrying the optimal basis between consecutive goals, which could work as a warm start for the simplex iterations. There is currently a warm start from the physical layer, but it's set once per time step, and presolve may discard it anyway.