Piecewise linear basin profile approx in allocation

[SouthEndMusic]

Fixes #2060.

An additional complexity here is that because we change the way we represent basins in the allocation problems, level demand has to be reformulated as well. The previous approach was to have a 'buffer' attached to basins with a level demand which acted as a source when the level in the basin was too high, and as a UserDemand when the level in the basin was too low. This source or demand was computed with data from the physical layer.

This PR introduces the level of each basin explicitly as a variable in the optimization problems. Therefore we can optimize for the level directly. The error for other demand nodes are still formulated in terms of flow though, so the error for level demand should be formulated that way too, so all terms are treated the same way in the loss function. Therefore we formulate the error term for the level demand as

$$ \left| \frac{L_\text{end} - l_\text{target}}{\Delta t_\text{allocation}} \right|, $$

where $L_\text{end}$ is the level of a particular basin at the end of the allocation timestep.

One thing to properly think trough here is what the behavior of a basin with a level demand should be when optimizing for a different priority than the priority of the level demand. I would say intuitively: for priorities higher than the level demand, all storage of the basin should be available as a source, and for the priority of the level demand and lower priorities only the disk of water above the maximum tolerated level should be available, if applicable.

[SouthEndMusic]

Here's a bit of theoretical background to the goal of this PR.

For each basin, a table is given with levels and associated areas: $h_i, A_i \in \mathbb{R}, i = 1, \ldots n$. This data is sorted by level and the areas are strictly positive and non-decreasing. To get a continuous profile we interpolate this data linearly:

$$ A(h) = A_i + (A_{i+1} - A_i)\frac{h - h_i}{h_{i+1} - h_i}, \quad h \in [h_i, h_{i+1}]. $$

Then to get the volume (storage) of water at some level $\ell$ we integrate from the bottom, given by the first and lowest level in the table:

$$ V(\ell) = \int_{h_1}^\ell A(h)\text{d}h, $$

giving a piecewise quadratic storage - level relationship. This is the relationship we want to approximate linearly. In the code we actually use the inverse of this relationship (storage_to_level) which is more useful in the physical layer, which I implemented upstream and is described here. In this PR I evaluate this relationship in an arbitrary `n_samples_per_interval = 5 equi-spaced points in each interval between data points.

The level can also increase above the last level $h_n$, for which we extend the linear area and quadratic storage increase from the last data interval in the physical layer. It seems however to do this in the allocation problem that we have to make some educated guess for the maximum level that is going to be attained at model initialization, but @visr says that in practice this is probably doable.

[jarsarasty]

Hi @SouthEndMusic,

Originally, the objective function of this problem is quadratic, and this pull request introduces binary decision variables, so your problem becomes a mixed-integer quadratic one. Unfortunately, there are no very good open source solvers for this type of discrete nonlinear optimization problem. In particular, HiGHS cannot solve QP models where some of the variables must take integer values. You could try other free solvers like bonmin, but the solution will probably not be very fast.

Other alternatives, as discussed in our meeting last week, would be 1) to convert the objective function into a linear one (this is probably a no regret action), 2) to try to exploit the problem structure to avoid the use of binary variables (this will require a bit more investigation, and is not guaranteed to be achievable), or 3) try to solve the problem in stages, e.g. by first solving a simplified version with linear storage level tables, and then using that initial solution to solve a more detailed model with the piecewise storage curve representation (this staged approach will require some effort to set up).

I would suggest trying first action 1), implemented as an option, so that the user can choose between the quadratic or linear objective. We could discuss this approach further and refine what has been suggested here.

[visr]

It makes sense to try to use a linear objective function. I'd prefer to not support multiple objective functions unless there is a good reason for it, to keep the maintenance lower.

[SouthEndMusic]

@jarsarasty thanks for your reply. I'll start implementing the linear objective function you suggested.

[SouthEndMusic]

The current implementation of the storage level relationship is like this:

Say for basin $i$ we have the storage data and level data $s_{i, j}, h_{i, j}$ for $j = 1, \ldots, n$. Then we define:

• The Boolean variables $B_{i,j}$ for $j = 1, \ldots, n - 1$ ;
• The auxiliary variables $A_{i,j} \ge 0$ for $j = 1, \ldots n$;
• The storage $S_i$ after the time step;
• The level $H_i$.

Then we define the constraints:

• Exactly one binary variable is $1$:

$$ \sum_{j=1}^{n - 1} B_{i,j} = 1; $$

• The auxiliary variables yield a convex combination:

$$ \sum_{j=1}^{n - 1} A_{i,j} = 1; $$

• The storage (only one of these terms contributes) (quadratic expression in the variables):

$$ S_i = \sum_{j=1}^{n-1} B_{i,j} (s_{i,j}A_{i,j} + s_{i,j+1}A_{i,j+1}); $$

• The level (only one of these terms contributes) (quadratic expression in the variables):

$$ H_i = \sum_{j=1}^{n-1} B_{i,j} (h_{i,j}A_{i,j} + h_{i,j+1}A_{i,j+1}). $$

Now that I write this out, I doesn't make that much sense to me anymore. It seems that I only need 1 or 2 auxiliary variables to define the linear combination of the successive storage and level data. But then still I end up with the product of the binary variable and auxiliary variable.

Asking Claude about this, it seems that what I should do is simply have

$$ S_i = \sum_{j=1}^n A_{i,j}s_{i,j}, \quad H_i = \sum_{j=1}^n A_{i,j}h_{i,j} $$

and then make sure that only 2 consecutive auxiliary variables are non-zero by

$$ A_{i,j} \le B_{i,j-1} + B_{i, j} \quad \text{ for } j = 2, \ldots, n-1 $$

and $A_{i,1} \le B_{i,1}$, $A_{i,n} \le B_{i, n-1}$. This avoids products of variables.

[jarsarasty]

Hi @SouthEndMusic,

Indeed, the product of decision variables in the optimization model should be avoided because it introduces nonlinearity, which can significantly complicate the problem and make it harder to solve.

[SouthEndMusic]

@jarsarasty this package looks very interesting for constructing linear approximations:

https://github.com/LICO-labs/PiecewiseLinApprox.jl