Native multi objective optimization

[SouthEndMusic]

Fixes #2233 (To be precise: the multi-objective optimization is now handled by MultiObjectiveAlgorithms.jl and not by HiGHS, this was advised here.)
Fixes #2323 (it has been made a goal anyway, the reduction factors are maximized)
Fixes #2507
Fixes #2399

In the following we assume lexicographic multi-objective optimization.

UserDemand and FlowDemand

Say we have a subnetwork with UserDemand nodes $UD$ and FlowDemand nodes $FD$ and each node $v \in UD \cup FD$ has demands for priorities $P_v$. Then for every one of these demand nodes for every priority for which they have a demand (at some point in time) we define a special flow decision variable which denotes how much is allocated to that node for that priority:

$$ 0 \le F^p_v \le d_v^p, \quad \forall v \in (UD \cup FD),\; p \in P_v. $$

As indicated, each of these allocated flow variables is bounded between 0 and the corresponding demand (NOTE: in this text we assume that the demands are already expressed in the scaled flow unit). For UserDemand nodes we have that the sum of these allocated flow variables over the priorities is equal to the total flow into that node:

$$ F_{(b_\text{upstream}, v)} = \sum_{p \in P_v} F^p_v, \quad \forall v \in UD. $$

Note that this means that we enforce that only allocated flow can enter a UserDemand node, and so there can not be more flow into the UserDemand node than the total demand of that node.
For FlowDemand we do allow that there is more flow through the node with the flow demand than the total demand (or even negative flow). To allow this freedom, we introduce a special 'demand priority 0' variable

$$ 0 \le F^0_v \le \frac{\text{MAX\_ABS\_FLOW}}{\text{scaling.flow}}, \quad \forall v \in FD. $$

We only let $F^0_v$ account for surplus flow, because if we allow that variable to be negative it can be used to allocate to demands. To make sure we still account for the possibility that the flow through the node with the flow demand is negative, we allow the allocated flow variable with the earliest priority to be negative, and account for this when writing the results.

With this we formulate the constraints

$$ F_{(b_\text{upstream}, v)} = \sum_{p \in P_v \cup \{0\}} F^p_v, \quad \forall v \in FD. $$

We define relative error variables per demand node per priority for which they have a demand:

$$ E^p_v \ge 0, \quad \forall v \in (UD \cup FD),\; p \in P_v. $$

The constraints on these relative error terms are similar as before:

$$ d_v^p \cdot E_v^p \ge d_v^p - F_v^p, \quad \forall v \in (UD \cup FD),\; p \in P_v. $$

With these relative error terms we formulate an objective per demand priority for which any of these nodes have a demand ($P^\text{flow} = \bigcup_{v \in UD \cup FD} P_v$):

$$ \min \sum_{v \in UD \cup FD \;: \; p \in P_v} d_v^p \cdot E_v^p, \quad \forall p \in P^\text{flow}. $$

Note the multiplication by the demand here, which makes this effectively a minimization in the absolute error. There is a reason we formulate relative in stead of absolute errors: to optimize for a fair distribution of water. To that end we can formulate per demand priority a global relative allocation error:

$$ G^p = \frac{\sum_{v \in UD \cup FD \;: \; p \in P_v} d_v^p \cdot E_v^p}{\sum_{v \in UD \cup FD \;: \; p \in P_v} d^p_v}, \quad \forall p \in P^\text{flow}, $$

where in practice we multiply both sides by the denominator to prevent division by zero problems. Also note that the enumerator is precisely the expression for the aforementioned objective.
Now in a subsequent objective, we want to penalize relative errors being larger than $G^p$. So we define new errors

$$ \overline{E}_v^p \ge 0, \quad \forall v \in (UD \cup FD),\; p \in P_v, $$

on which we define the constraints

$$ \overline{E}_v^p \ge E_v^p - G^p, \quad \forall v \in (UD \cup FD),\; p \in P_v $$

and with which we formulate the objectives

$$ \min \sum_{v \in UD \cup FD \;: \; p \in P_v} \overline{E}_v^p, \quad \forall p \in P^\text{flow}. $$

Of course the order of objectives is important. For each demand priority of the flow type described here, first the weighted sum of the errors $E_v^p$ must be minimized, then the sum of the errors $\overline{E}_v^p$.

Edge cases

I think this formulation is quite robust against edge cases. When a demand $d_v^p$ is zero the constraints on the corresponding allocated flow variable reduce to $F_v^p = 0$. The corresponding relative error variable $E_v^p$ is neither in any constraint nor in the objective function, is that bad or is that variable then filtered out in the presolve? When all demands of a particular demand priority $p$ are zero, $G_p$ is not in any constraint, will that also be handled by the presolve?

LevelDemand

Say we have basins with LevelDemand $B$ where each basin $b \in B$ has a minimum level $h_{b,\min}^p$ for demand priorities $P^\min_b$ and a maximum level $h_{b,\max}^p$ for demand priorities $P_b^\max$. Furthermore each basin has a level to storage function $s_b(h)$. For all basins $b$ in the subnetwork (not just the ones with a level demand) we also have

• The known initial level at the beginning of the allocation timestep $h_b^\text{init}$
• The decision variable of the change in storage over the allocation timestep $\Delta S_b$

For basins with a level demand $b \in B$ we additionally have

• The lower absolute storage error

$$ E_{b, \text{lower}}^p \ge 0 \quad \forall b \in B, \; p \in P^\text{min}_b $$

• The incoming storage demand

$$ d^p_{b, +} = s\left(h_{b, \min}^p\right) - s_b\left(\text{clamp}\left(h_b^\text{init}, h_{b,\min}^{\text{prev}(p)}, h_{b,\min}^p \right)\right) \quad \forall b \in B, \; p \in P^\text{min}_b $$

which is the volume of water needed to get to the minimum level of the current demand priority from the minimum level of the previous demand priority (or the starting level if that's higher). If there is no priority in $P_b^\min$ lower than the current priority $p$ we say $h^{\text{prev}(p)}_{b, \min}$ is equal to the Basin bottom.

• The upper absolute storage error

$$ E_{b, \text{upper}}^p \ge 0 \quad \forall b \in B, \; p \in P_b^\max $$

• The outgoing storage demand:

$$ d^p_{b, -} = s_b\left(\text{clamp}\left(h_b^\text{init}, h_{b, \max}^p, h_{b, \max}^{\text{prev}(p)}\right)\right) - s\left(h_{b, \max}^p\right) \quad \forall b \in B, \; p \in P^\text{max}_b $$

which is the volume of water needed to get from the maximum level of the current demand priority to the maximum level of the previous demand priority (or the starting level if that's lower). If there is no priority in $P_b^\max$ lower than the current priority $p$ we say $h^{\text{prev}(p)}_{b, \max} := \infty$.

We define the following constraints on the errors:

$$ E^p_{b, \text{lower}} \ge s\left(h^p_{b, \min}\right) -\left(s\left(h_b^\text{init}\right) + \Delta S_b\right) \quad \forall b \in S, \; p \in P^\min_b $$

$$ E^p_{b, \text{upper}} \ge \left(s\left(h_b^\text{init}\right) + \Delta S_b\right) - s\left(h^p_{b, \max}\right) \quad \forall b \in S, \; p \in P^\max_b. $$

Note that unlike with UserDemand and LevelDemand, the allocated amounts are not present as variables in the problem. The allocated amounts have to be derived from $\Delta S_b$. Also note that the demands, which depend on the state of the physical layer and the start of the allocation timestep, do not appear in the LP problem. These are only computed for output.

Then given all this we define the objective for minimizing the absolute error sum per priority for which there is at least one level demand ($P_\text{lower}^\text{level} = \bigcup_{b \in B} P^\min_b$, $P_\text{upper}^\text{level} = \bigcup_{b \in B} P^\max_b$, $P^\text{level} = P_\text{lower}^\text{level} \cup P_\text{upper}^\text{level}$):

$$ \min \sum_{b \in B \; : \; p\in P^\min_b} E^p_{b, \text{lower}} + \sum_{b \in B \; : \; p\in P^\max_b} E^p_{b, \text{upper}}, \qquad \forall p \in P^\text{level} $$

Now, just as in the previous section, we want to add a subsequent objective for fair distribution. Given the Basin area $A_b$ (which is assumed to be constant over the allocation timestep), we can define the average lower and upper level error

$$ \Delta H^p_\text{lower} = \frac{\sum_{b \in B \; : \; p\in P^\min_b} E^p_{b, \text{lower}}}{\sum_{b \in B \; : \; p\in P^\min_b} A_b}, \qquad \forall p \in P_\text{lower}^\text{level}; $$

$$ \Delta H^p_\text{upper} = \frac{\sum_{b \in B \; : \; p\in P^\max_b} E^p_{b, \text{upper}}}{\sum_{b \in B \; : \; p\in P^\max_b} A_b}, \qquad \forall p \in P_\text{upper}^\text{level}; $$

Again we define new errors:

$$ \overline{E}_{b, \text{lower}}^p \ge 0, \quad \forall b \in B, \; p \in P_b^\min; $$

$$ \overline{E}_{b, \text{upper}}^p \ge 0, \quad \forall b \in B, \; p \in P_b^\max. $$

On these errors we define the following constraints:

$$ \overline{E}^p_{b, \text{lower}} \ge \frac{E^p_{b, \text{lower}}}{A_b} - \Delta H^p_\text{lower}, \quad \forall b \in B, \;p \in P_b^\text{min}; $$

$$ \overline{E}^p_{b, \text{upper}} \ge \frac{E^p_{b, \text{upper}}}{A_b} - \Delta H^p_\text{upper}, \quad \forall b \in B, \;p \in P_b^\text{max}, $$

and the following objectives:

$$ \min \sum_{b \in B \; : \; p \in P_b^\min} \overline{E}^p_{b, \text{lower}} + \sum_{b \in B \; : \; p \in P_b^\max} \overline{E}^p_{b, \text{upper}}, \qquad \forall p \in P^\text{level}. $$

@jarsarasty

[SouthEndMusic]

An overview of changes I make in this PR (code deletions not included):

parameter.jl

• A restructure of the objectives data structure. There is now a AllocationObjectives structure, which has as fields a vector of all objective expressions in the right order to pass to the optimizer, and a vector of metadata objects (the kind of objective, priority if applicable, objective expressions). The idea is that the objective expressions can be updated via the metadata objects because these point to the same expressions as in the vector of all expressions. There are more objective expressions than metadata objects because demand objectives are split into 2 objectives (the wording here is a bit vague, but I hope you understand by what I explain above).
• There's now a NodeIDsInSubnetwork which is a new field of the allocation models. It has fields which tell you of relevant node types which nodes of that type are in the subnetwork. This is nicer than obtaining those from the indices of the JuMP variables and constraints, which is what I did before.
• Basin and structs for connector node types now know which node provides them with a demand with the level_demand_id or flow_demand_id field.
• Quite a few fields were removed from LevelDemand. This is because this data doesn't have to be cached anymore because with the refactor they are computed and then used all in one place.

read.jl

• Populate the level_demand_id and flow_demand_id fields

allocation_util.jl

• I created a custom iterator DemandPriorityIterator which gives the demand priorities for which a particular node has a demand based on the has_demand_priority field of the demand node structs.
• I created the function get_external_demand_id which for a given node ID gives the LevelDemand or FlowDemand node that defines the demand of the input node. The method that takes the graph as input should only be used in initialization because this method uses a relatively expensive dictionary lookup.

allocation_init.jl

• I started using the fact that variable and constraint containers can have multiple indices. For example user_demand_allocated has a node ID index and a demand priority index. Here and in many other places I use the DemandPriorityIterator to only define variables and constraints for demand priorities for which the particular node has a demand (so it's not simply a Cartesian product, see e.g. here the part Sets can depend upon previous indices).
• For the UserDemand and FlowDemand nodes I now define an allocated amount variable per demand priority for which there is a demand. For UserDemand, these sum up to the total inflow of the UserDemand node. FlowDemand is very similar, but there is a flow_demand_extra variable which allows for the flow to be more than the total demand, and earliest priority demand priority allocation is allowed to be negative.
• As mentioned, for each demand priority there are 2 objectives; one for how much can be allocated in total and one for fair distribution of that total. I wasn't inspired to give these informative names, so I just refer to them as 'first' and 'second'.
• Objective initialization was moved to the end of the AllocationModel constructor instead of partly at the start and partly in the add_*! functions. In the function add_demand_objectives! also the variables and constraints for the 'second' objectives are defined, which rely on variables for the 'first' objectives already existing.
• I removed the piecewise linear function for the low storage factor of each Basin, so now all binary variables are removed. There is still a low storage factor variable bounded in [0, 1] for each Basin (initialized in add_basin!), but its value is unconstrained until the sum of the low storage factors is maximized in the low_storage_factor objective type (see @enumx AllocationObjectiveType).
• I split the Basin forcing into a positive and a negative term, see add_conservation!. These are no longer JuMP variables by themselves; the positive term is a constant and the negative term is a coefficient of the low storage factor in the volume conservation constraints.

Linearizations

• I suggested before to formulate the physics in terms of the storage change over the allocation timestep instead of the (known) storage at the start of the timestep and the (unknown) storage at the end of the allocation timestep. That's what I did, and it makes the linearizations even simpler, because the level dependencies can just be formulated as the partial derivative w.r.t. the level times the storage change divided by the area without any rhs contribution
• In line with implicit Euler, I evaluate levels and flows using the end time of the allocation timestep in stead of the start time (see t_after in linearize_connector_node!).
• Linearization with respect to the level is only needed when that level comes from a Basin. I apply this in add_linearized_connector_node! and linearize_connector_node!.

allocation_optim.jl

• Infeasiblity analysis is a bit different: we now cannot say at which objective the infeasibility was introduced. However, it's highly unlikely that the infeasibility is introduced during the lexicographic algorithm, because the constraints that are added there are known to be satisfiable. I also found that writing the problem to a .lp file with a vector-valued objective is not supported, so I have to change the objective (and the optimizer with in) to a scalar one before doing so.
• Goal programming resetting is no longer necessary (it might be worth checking this to be sure though: do the constraints added by MultiObjectiveAlgorithms.jl remain in the model?)
• Demands are set according to the descriptions above
• The warm start takes the current instantaneous flow rates from the physical layer as initial guess for the flows in allocation

[SouthEndMusic]

I have a problem with LevelDemand and FlowDemand that I don't know how to solve. Let's look at FlowDemand as an example:

Say I have a connector node with a FlowDemand $d \ge 0$ for a single priority. Let $F$ be the flow through that node in the allocation LP. Let $0 \le A \le d$ be the amount allocated to the connector node. Then

$$ F = A + F_\text{extra}, $$

where $F_\text{extra}$ is flow through the node that is either negative (which can happen for resistance nodes) or the flow trough the connector node that is more than the demand. Because we have these 2 possibilities, we cannot put tight bounds on $F_\text{extra}$. The problem that arises is that $F_\text{extra}$ can be used to allocate; say $A = d$, $F = 0$ and $F_\text{extra} = - d$, then the demand is satisfied despite no net flow going through the connector node.

For FlowDemand we could circumvent this by only allowing FlowDemand to connect to Pump and Outlet, so that we do not have to take the negative flow scenario into account and we can restrict $F_\text{extra} \ge 0$. For LevelDemand however the problem is not fixable in a similar way.

Annoyingly, binary variables can be used to fixed this, while I just got rid of them all. The binary variable basically says: if I'm in a positive flow scenario the extra flow cannot be negative, and vice versa.

The drain_surplus model got me into this problem, since it made a test fail.

@gijsber

[SouthEndMusic]

There is a way to fix this for FlowDemand that somewhat ties in to our discussions of today @gijsber. In the example above, we can set $F_\text{extra} \ge 0$, so that it only signifies flow over the total demand. Then for the first priority allocated value we remove the lower bound of $0$, so that negative flow ends up there. We could decide to clip that to $0$ in the output.

For LevelDemand this stil doesn't work because demands can be either positive or negative. My preference to solve this would be to get rid of the maximum level in LevelDemand, because as discussed earlier, in practice excess storage in reservoirs is drained by physical means and not optimized for, and in most of the discussions today with @Fati-Mon we only spoke about minimum levels. If in the future there comes a request for maximum levels this could maybe be achieved by splitting the optimization for the maximum level and the minimum level per priority in different goals.

Edit: I think actually it might be fixable for LevelDemand too because it is known a priori whether the Basin is in a demand or surplus situation.

[SouthEndMusic]

OK, the dust has finally settled on this PR. The PR description and the overview of changes are up to date. I think the end result is really nice:

• The code is simplified quite a bit
• It uses a more standard approach to multi-objective optimization
• I'm quite happy with the fair distribution approach which has been a difficult aspect from the start
• The performance of models that use allocation has increased greatly, at least in part because there are no longer any binary variables

[jarsarasty]

Hi @SouthEndMusic
Regarding your question:

The corresponding relative error variable $E_v^p$ is neither in any constraint nor in the objective function, is that bad or is that variable then filtered out in the presolve? When all demands of a particular demand priority $p$ are zero, $G_p$ is not in any constraint, will that also be handled by the presolve?

Yes, we expecte these variables to be filtered out by presolve.

[jarsarasty]

This implementation should significantly improve performance. Further testing is recommended to evaluate the effectiveness of this new approximation approach that does not use binary decision variables. Excellent job, @SouthEndMusic!

[verheem]

Nice work, this was a huge chunk. I have only some small remarks

All comments were addressed