Systematic linearization for allocation

[SouthEndMusic]

The tooling suggested in #3113 for evaluating the Jacobian can also be used for linearizations of the physics for allocation in a systematic way. Not only will this make this part of the allocation code simpler, it is also a good opportunity to:

• Let allocation take PID control into account;
• Let allocation take ContinuousControl into account (Let allocation take ContinuousControl into account #2322);
• Use the higher order trapezoid integration for physical flows.

To see how, let:

• $\Delta S$ be the unknown storages change over the timestep $\Delta t$;
• $Q$ be all unknown flows at the end of the timestep;
• $S_\text{p}$ be the matrix selecting the unknown physical flows at the end of the timestep (that includes PID and continuously controlled flows);
• $S_\text{c}$ be the matrix selecting the unknown allocation controlled flows (structures) at the end of the timestep;
• $S_\text{UD, in}, S_\text{UD, out}$ be the matrices selecting the UserDemand inflow and outflow respectively at the end of the timestep;
• $S_\text{PID}$ be the matrix selecting the PID controlled Basins;
• $P$ be the permutation that sorts $Q$ into sections of physical flows, allocation controlled flows and UserDemand inflows for notational convenience:

$$ P = \begin{bmatrix} S_\text{p} \\ S_\text{c} \\ S_\text{UD, in} \\ S_\text{UD, out} \end{bmatrix}; $$

• $f$ be the known forcing term onto the Basins;
• $\Delta I$ be the unknown change in the PID integral terms over the timestep;
• $\Delta C$ be the unknown change in the continuous control compound variables over the timestep;
• $M$ be the signed incidence matrix of the model graph;
• $q$ be the flow function;
• $\hat{q}$ be the instantaneous flows now;
• $\hat{h}$ be the Basin levels now;
• $\hat{a}$ be the Basin areas now;
• $r(t)$ be the PID control setpoint;
• $\rho(t)$ be the UserDemand return fraction;
• c be the continuous control compound variable function.

Then we have the following equality constraint sets:

• The UserDemand return flow (averaged return fraction over timestep): $S_\text{UD, out} Q = \left(S_\text{UD, in}Q\right) \odot \frac{1}{\Delta t}\int_{t_n}^{t_n + \Delta t} \rho(\tau)\text{d}\tau$.
• The water balance, using trapezoid integration for the physical flows and reverse fill integration for the allocation controlled flows and UserDemand inflows:

$$ \begin{align} \Delta S &=& \Delta t \left[ M\left[ P^\top \begin{bmatrix} I_\text{p} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{bmatrix} P\frac{1}{2}\left(\hat{q} + Q\right) + P^\top \begin{bmatrix} 0 & 0 & 0 & 0 \\ 0 & I_\text{c} & 0 & 0 \\ 0 & 0 & I_\text{UD} & 0 \\ 0 & 0 & 0 & I_\text{UD} \end{bmatrix} PQ \right] + f \right] \\ &=& \Delta t \left[ MP^\top \left( \begin{bmatrix} \frac{1}{2}I_\text{p} & 0 & 0 & 0 \\ 0 & I_\text{c} & 0 & 0 \\ 0 & 0 & I_\text{UD} & 0 \\ 0 & 0 & \text{diag}\left(\frac{1}{\Delta t}\int_{t_n}^{t_n + \Delta t} \rho(\tau)\text{d}\tau\right) & 0 \end{bmatrix} PQ + \begin{bmatrix} \frac{1}{2}S_\text{p}\hat{q} \\ 0 \\ 0 \\ 0 \end{bmatrix} \right) + f \right] \end{align} $$

• PID Integral terms (exactly integrated setpoint): $\Delta I = \int_{t_n}^{t_n + \Delta t}r(\tau)\text{d}\tau - \Delta t S_{PID}\left(\hat{h} + \frac{\Delta S}{2\hat{a}}\right)$;
• The continuous control compound variables: $\Delta C = \frac{\partial c}{\partial s}\Delta s + \frac{\partial c}{\partial q}\left(Q - \hat{q}\right)$;

We linearize $q$ to obtain for the physical flows:

$$ S_\text{p}Q = S_\text{p}\left[\hat{q} + \frac{\partial q}{\partial s}\Delta S + \frac{\partial q}{\partial i}\Delta I + \frac{\partial q}{\partial c} \Delta C\right]. $$

The partial derivatives needed here are exactly the ones computed for the Jacobian in #3113.

Note: This will not work when continuous control (or PID control) acts across subnetworks, then we would have have to assume fixed integral terms and continuous control compound variables with lower accuracy (or should this throw an error, @gijsber).

[gijsber]

At this moment I would say it is fair, to that continupous control can not work across subnetworks (and thus throw an error). By the time we see that this limits the modeller too much, we probably have to assume fixed integral terms and continuous control compound variables with lower accuracy