[style=theoremstyle]definitionHydraulic Dynamics
[style=theoremstyle]definitionqQuality Dynamics
[style=theoremstyle]definitionhcDecoupled Hydraulic Control
[style=theoremstyle]definitionqhc1Quality-Aware Hydraulic Control — Rank-
Informed [style=theoremstyle]definitionqhc2Quality-Aware Hydraulic Control —
Energy-Driven [style=theoremstyle]definitionqhc3Quality-Aware Hydraulic
Control — #3

# Quality-Aware Hydraulic Control in Drinking Water Networks via
Controllability Proxies

Salma M. Elsheri$\displaystyle\text{f}^{\dagger,\P}$, Mohamad H.
Kazm$\displaystyle\text{a}^{\dagger}$, and Ahmad F.
Tah$\displaystyle\text{a}^{\dagger,\ast\ast}$ †Department of Civil and
Environmental Engineering, Vanderbilt University, Nashville, TN, USA. Emails:
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>appointment: Department of
Irrigation and Hydraulics Engineering, Faculty of Engineering, Cairo
University.∗∗Corresponding author.This work is supported by the National
Science Foundation under Grant 2151392.

###### Abstract

The operation of water distribution networks is a complex procedure aimed at
efficiently delivering consumers with adequate water quantity while ensuring
its safe quality. An added challenge is the dependency of the water quality
dynamics on the system’s hydraulics, which influences the performance of the
water quality controller. Prior research has addressed either solving the
optimum operational hydraulic setting problem or regulating the water quality
dynamics as separate problems. Additionally, there have been efforts to couple
these two problems and solve one compact problem resulting in trade-offs
between the contradictory objectives. In contrast, this paper examines the
dependency and influence from a control-theoretic standpoint. More
specifically, we explore the influence of accountability for water quality
controllability improvement when addressing the pump scheduling problem. We
examine its effects on the cumulative cost of the interconnected systems as
well as the subsequent performance of the water quality controller. To achieve
this, we develop a framework that incorporates different controllability
metrics within the operational hydraulic optimization problem; its aim is
attaining an adequate level of water quality control across the system. We
assess the aforementioned aspects’ performance on various scaled networks with
a wide range of numerical scenarios.

###### Index Terms:

Optimal pump schedule, hydraulics and water quality regulation and control,
controllability metrics, model predictive control.

##  I Introduction and Literature Review

The real-time management and operation of water distribution networks (WDNs)
have been one of the most researched topics in the field of water systems. The
objective is to fulfill consumers’ and end-users’ needs and deliver clean
water in a cost-effective manner while meeting water quality (WQ) mandates.
Taking into account network-wide operational flows, pressures, and WQ,
achieving the aforementioned objective involves addressing multiple
conflicting objectives. These objectives include minimizing the costs
associated with pumping and disinfectant usage, meeting water demand,
maintaining adequate levels, and satisfying quality requirements.

For decision-makers (i.e., water utilities), these conflicting objectives
formulate challenges to be addressed utilizing control algorithms that are
built and applied in WDNs. The regulation and control of WDNs primarily
revolve around two aspects: quantity and quality. For the quantity aspect, the
energy required to operate pumps constitutes the main component of the
operational cost to be minimized while achieving the target head levels and
flows [1]. As for WQ, the objective of the control problem is to maintain
desired disinfectant levels all over the WDN with minimum injections at
treatment plants and booster stations [2].

The vintage approach to study WQ dynamics starts with first running hydraulic
simulations that generate schedules of pumps with their resulting heads and
flows through a WDN. This is then followed by investigating quality via
booster station control. This is due to the fact that hydraulic control has a
slower time-constant. This paper pursues a new approach: jointly investigating
quality-quantity control via a unified approach based on systems science and
controllability metrics of quality dynamics.

Controllability, in this context, refers to the ability to effectively steer,
regulate, and maintain disinfectant levels within the network to consistently
meet the established water health standards. In short, the paper attempts to
answer this research questions: Can quality-aware pump control significantly
improve water quality dynamics in WDN? What are the quality controllability
metrics that can be appended to a hydraulic control problem? When is it
meaningful to integrate the time-scales of quality and quantity?

To the best of our knowledge, this is the first attempt to explore this
particular topic from a system and control theoretic perspectives. The next
section delves into the corresponding literature.

###  I-A Literature Review

The literature on the regulation and control of WDNs to deliver clean water to
the end-users is rich and briefly summarized next. The literature on this
topic is divided into (i) studies focusing on determining the optimal
operational settings for pumps and/or control valves to attain the desired
water flows and levels specified by network topology and characteristics, and
consumers demands; (ii) studies that cover the regulation of WQ dynamics to
ensure meeting the standard disinfectant (i.e., chlorine) residuals throughout
the network while minimizing the source and booster stations’ injections; and
(iii) studies that propose integrating/coupling controlling water quantity and
quality, therefore considering the interdependencies between them. In this
section, we survey these areas and afterward, we end this section by
highlighting the research gaps that drive the contributions of this work.

Hydraulic Control. The operational control of WDNs depends on several factors
including network topology, demand cycle, tank dynamics, head loss model,
pumps type, valves type, etc. In addition, the resulting control problem is
inherently nonlinear and nonconvex in nature. Many studies have covered some
or many of the aforementioned factors and handled the nonlinearity and
nonconvexity under different frameworks. Studies [1, 3] perform system
linearization and apply linear programming to obtain optimal pump scheduling.
Study [4] determines the optimal scheduling by the means of relaxation and
linear programming branch and bound. Similarly, study [5] relaxes the
hydraulic constraints to second-order cone constraints with penalty terms. On
the other hand, studies [6, 7, 8] solve nonlinear nonconvex problem using
variations of methods including genetic algorithm and mixed-integer nonlinear
programming. Lastly, study [9] applies geometric programming-based model
predictive control (MPC) algorithms which turn the problem into a convex
continuous optimization problem. Many of these studies have compared their
results with built-in tools in hydraulics solvers (e.g., EPANET’s built-in
rule-based control [10])—refer to [11] for detailed review and analysis of the
literature on this topic.

Quality Control. A plethora of optimization-based approaches have been used to
solve the WQ control problem including: linear programming with the objective
of minimizing chlorine injections [12], mixed-integer linear programming with
the allocation of booster stations as a decision variable [13], and genetic
algorithm with a constraint on the formulation of disinfectant by-products
[14]. Conversely to these studies, study [15] proposes applying MPC to an
explicit representation of the single-species WQ model that guarantees
network-wide control. More on that approach, study [16] utilizes different
techniques to implement real-time MPC to nonlinear multi-species WQ dynamics—a
framework that covers controlling chlorine levels in WDNs under abnormal
conditions including contamination events. It is worth noting that these
studies rely on the assumption that systems’ hydraulics are pre-computed.

Joint Quality-Quantity Control. Several studies have investigated integrating
both the quantity and quality control problems by implicitly and/or explicitly
incorporating one or more quality control aspects within the quantity control
framework, or by turning them into one augmented formulation. In [17], the
authors use nonlinear programming to solve the pump operation problem that
accounts for disinfectant’s residuals in the constraints. On the other hand,
study [18] formulates a dual quality-quantity optimization problem with a
single augmented objective that concatenates minimizing the energy cost and
maximizing system’s protection by maximizing the injected chlorine dose. This
study utilizes a genetic algorithm to solve the optimal problem that is based
on the two conflicting objectives. The posed problem is optimally solved while
conveying the existence of trade-offs within the solution. Similarly, authors
in [19] propose applying a genetic MPC algorithm to the coupled control
problem and compare the results with real data records of a specific network,
which shows cost reduction. The authors in [20] utilize a nonlinear MPC
integrated optimizer. The control procedure proposed is divided into two
levels: an upper-level controller responsible for determining optimized pump
schedules while satisfying constraints on chlorine residuals, and a lower-
level controller that computes optimized chlorine injections. More recently,
study [21] solves a two-objectives pump scheduling problem by means of goal
programming. The first objective focuses on implicitly achieving the required
chlorine residuals by minimizing the active dynamics in storage components,
while the other objective aims towards minimizing the energy cost.

Many of the studies cited earlier on this topic solve a fully coupled quantity
and quality control problem by concatenating the objectives of each problem
into a single-objective or multi-objective formulation. This approach results
in trade-offs between conflicting objectives. However, a comprehensive
analysis of these trade-offs has not been conducted. On the other hand,
incorporating a constraint on WQ, whether implicitly or explicitly, into the
system’s operational scheduling control problem does not necessarily guarantee
the achievement of a certain level of controllability by booster stations or
reachability of the desired final states. In other words, the notion of
optimally designing a pumping schedule while attaining a certain level of WQ
controllability, based on closed-form formulation of the system’s hydraulics
and quality, has not been attempted nor investigated—a gap to be addressed by
this paper.

###  I-B Paper Contributions

This paper’s main objective is to investigate the feasibility and
applicability of maintaining a certain level of controllability for the WQ
model while computing the optimal hydraulic setting. This implies solving an
augmented water networks operational control problem that accounts for
enhancing WQ from a control-theoretic perspective. In this context, we refer
to the control problem that focuses solely on hydraulics as the decoupled
problem, while the joint quality-quantity problem as the coupled problem. The
corresponding paper contributions are as follows.

  * •

Water quality systems are inherently complex and mostly not fully
controllable. That is, we investigate the effect of changing hydraulics on the
WQ controllability. This is measured by employing different quantitative
metrics which allows us to judge the WQ system controllability from different
energy-related perspectives. Eventually, we judge the applicability and
validity of these metrics on the case-oriented application under our focus.

  * •

We formulate an augmented operational pump scheduling control problem—the
coupled problem—in a way that conserves a certain level of WQ controllability.
This level is pre-determined depending on the investigation results performed
as outlined in the previous point.

  * •

A wide comparison between the decoupled and coupled control problems through
various numerical case studies is performed.

The paper organization is as follows. We start by presenting the hydraulic and
water quality models adopted in Section II. Then in Section III, we delineate
several controllability metrics—and assess their applicability and validity to
our application. Next, we formulate the pumps scheduling optimization problem
with and without the integration of the WQ controllability preservation in
Section IV. We compare between the decoupled and coupled problems through
different case studies in Section V. Lastly, conclusion, paper’s limitations,
and future work are discussed in Section VI.

##  II Hydraulics and Water Quality Models

We model WDNs by a directed graph
$\displaystyle\mathcal{G}=(\mathcal{N},\mathcal{L})$. The set
$\displaystyle\mathcal{N}$ defines the nodes and is partitioned as
$\displaystyle\mathcal{N}=\mathcal{J}\cup\mathcal{T}\cup\mathcal{R}$ where
sets $\displaystyle\mathcal{J}$, $\displaystyle\mathcal{T}$, and
$\displaystyle\mathcal{R}$ are collections of junctions, tanks, and
reservoirs. Let
$\displaystyle\mathcal{L}\subseteq\mathcal{N}\times\mathcal{N}$ be the set of
links, and define it as,
$\displaystyle\mathcal{L}=\mathcal{P}\cup\mathcal{M}\cup\mathcal{V}$, where
sets $\displaystyle\mathcal{P}$, $\displaystyle\mathcal{M}$, and
$\displaystyle\mathcal{V}$ represent the collection of pipes, pumps, and
valves. In each network, $\displaystyle n_{\mathrm{L}}$ and $\displaystyle
n_{\mathrm{N}}$ represent the numbers of links and nodes. Specifically,
network’s nodes include numbers of $\displaystyle n_{\mathrm{R}}$ reservoirs,
$\displaystyle n_{\mathrm{J}}$ junctions, and $\displaystyle n_{\mathrm{TK}}$
tanks. The total number of links is the summation of $\displaystyle
n_{\mathrm{M}}$, $\displaystyle n_{\mathrm{V}},$ and $\displaystyle
n_{\mathrm{P}}$, representing the numbers of pumps, valves, and pipes,
respectively.

In our paper, we have two models (i.e., hydraulic and water quality models)
with different numbers of states and representations. Next, we succinctly list
the governing equations for both models and their final representations. It is
worth mentioning that the hydraulic time-step $\displaystyle\Delta
t_{\mathrm{H}}$ is different than the water quality one $\displaystyle\Delta
t_{\mathrm{WQ}}$. The hydraulic time-step is taken to be within an hourly
scale to reflect the patterned demand, while the water quality time-step is
chosen between minutes and seconds to allow a stable accurate numerical
simulation [22]. Then, variable $\displaystyle t$ represents a specific time
in the simulation period $\displaystyle[0,T_{s}]$ and it is updated
incrementally by $\displaystyle\Delta t_{\mathrm{WQ}}$ within each
$\displaystyle\Delta t_{\mathrm{H}}$ reaching the end of the simulation period
(i.e., $\displaystyle t=T_{s}$).

###  II-A Modeling Hydraulics

We apply the principles of conservation of mass and energy to obtain the
amount of water flowing in each network link and the head at each node. For
all the network components, we give brief description of the equations that
model these principles in the next section. In the latter section and based on
these equations, we formulate a compact state-space representation of the
network hydraulics.

####  II-A1 Water Network Components Modeling

For each of the following network elements, we model the hydraulics variables
at/through this element depending on their characteristics and their
connection to other elements.

– Reservoirs: We follow the valid assumption that reservoirs are infinite
source of water with fixed head [23]. Thence, the head at Reservoir
$\displaystyle i$ is calculated as $\displaystyle h_{i}^{\mathrm{R}}(t+\Delta
t_{\mathrm{H}})=h_{i}^{\mathrm{R}}(t)$.

– Tanks: We consider tanks with constant cross section along its height.
Change in the head at the tank depends on its cross section’s area and the
algebraic difference between the inflows and outflows. The head at Tank
$\displaystyle i$ is described as

$h_{i}^{\mathrm{TK}}(t+\Delta
t_{\mathrm{H}})=h_{i}^{\mathrm{TK}}(t)+\frac{\Delta
t_{\mathrm{H}}}{A_{i}^{\mathrm{TK}}}\Big{(}\sum_{j\in
L_{\mathrm{in}}}{q_{\mathrm{in}}^{j}(t)}-\sum_{k\in
L_{\mathrm{out}}}{q_{\mathrm{out}}^{k}(t)}\Big{)},$ (1)

where $\displaystyle A_{i}^{\mathrm{TK}}$ is the tank cross section’s
area;$\displaystyle j$ and $\displaystyle k$ are the counters for total
$\displaystyle L_{\mathrm{in}}$ links flowing into the node and $\displaystyle
L_{\mathrm{out}}$ links extracting flow from the node; and $\displaystyle
q_{\mathrm{in}}^{j}(t)$ and $\displaystyle q_{\mathrm{out}}^{k}(t)$ are the
inflows and outflows from these links connected to the node.

– Junctions: The conservation of mass law at junctions is expressed as

$\sum_{j\in L_{\mathrm{in}}}{q_{\mathrm{in}}^{j}(t)}-\sum_{k\in
L_{\mathrm{out}}}{q_{\mathrm{out}}^{k}(t)}=q^{\mathrm{D}_{\mathrm{J}}}_{i}(t),$
(2)

where $\displaystyle q^{\mathrm{D}_{\mathrm{J}}}_{i}(t)$ is the consumers
demand withdrawn from this junction.

– Pipes: As water flows in a pipe, it starts with the head of the upstream
node and reaches the end-node head by applying the conservation of energy. The
difference between the two heads is due to the losses caused by friction along
the pipe’s length and localized minor losses (e.g., at bends, fittings, etc.).
In our study, we neglect the minor losses as they are relatively small in
comparison to the friction losses in water networks [24]. That is, the change
in head through Pipe $\displaystyle i$ connecting and flowing water between
node $\displaystyle j$ and node $\displaystyle k$ is expressed in Eq.
(3)—these nodes can be junctions, tanks, or reservoirs.

$\leavevmode\nobreak\ \Delta
h^{\mathrm{P}}_{i}(t)=h_{j}(t)-h_{k}(t)=r_{i}q^{\mathrm{P}}_{i}(t)|q^{\mathrm{P}}_{i}(t)|^{\mu-1},$
(3)

where $\displaystyle q^{\mathrm{P}}_{i}$ is the pipe flow; $\displaystyle
r_{i}$ is the pipe resistance coefficient, which is a function of pipe size,
length, and material; and $\displaystyle\mu$ is the constant flow exponent.
These parameters’ values depend on the chosen head loss formula. In our study
we use the Hazen-Williams equation; $\displaystyle\mu=1.852$ [25].

– Pumps: As an active component with variable speeds, pumps can provide the
system with different values of head gain according to its operating speed and
the corresponding head-flow relationship [25]. For Pump $\displaystyle i$
adding energy to water flowing from node $\displaystyle j$ to node
$\displaystyle k$, the head gain is calculated as

$\leavevmode\nobreak\ \Delta
h^{\mathrm{M}}_{i}(t)=h_{j}(t)-h_{k}(t)=-s^{2}_{i}(h_{i}^{0}-\alpha_{i}(s^{-1}_{i}q^{\mathrm{M}}_{i}(t))^{\nu_{i}}),$
(4)

where $\displaystyle s_{i}$ is the pump relative speed varying between 0 and
$\displaystyle s_{i}^{\max}$; $\displaystyle h_{i}^{0}$ is the shutoff head;
$\displaystyle q^{\mathrm{M}}_{i}$ is the pump flow; and
$\displaystyle\alpha_{i}$ and $\displaystyle\nu_{i}$ are pump characteristics
coefficients. Note that the head gain is strictly negative, as the pump
provides the water with more energy, causing the head at the delivery node to
exceed that at the suction node, and no back-flow is allowed.

– Valves: In our model, we consider pumps to be the only controller of the
system. Therefore, we formulate and solve pump operation problem to obtain the
optimal pumping schedule that fulfills the water flow and head constraints.
That being the case, valves in our networks are considered as on-off valves
and they have two states; fully closed or fully open. The knowledge of valve
state is predetermined along the simulation period. In the case of fully
closed, the two nodes connected by the valve are considered decoupled. For the
other case of fully open valve, it is treated as a straight pipe section with
minor losses [26] that can be expressed as in Eq. (5) for Valve $\displaystyle
i$ connecting nodes $\displaystyle j$ and $\displaystyle k$.

$\leavevmode\nobreak\ \Delta
h^{\mathrm{V}}_{i}(t)=h_{j}(t)-h_{k}(t)=m_{i}q^{\mathrm{V}}_{i}(t)|q^{\mathrm{V}}_{i}(t)|,$
(5)

where $\displaystyle q^{\mathrm{V}}_{i}$ is the flow through the valve and
$\displaystyle m_{i}$ is the minor losses coefficient that depends on the
valve type (e.g., ball, butterfly, gate, etc.) and its cross-sectional area.

####  II-A2 Hydraulics in State-Space Form

The hydraulics model explained in Section II-A1 can be written in a form of
nonlinear difference algebraic equations (NLDAE) as expressed in (6). In these
equations, we collect the system state variables and inputs in vectors of
appropriate dimensions as follows: heads at tanks in
$\displaystyle\bm{\mathrm{w}}$; heads at junctions at
$\displaystyle\bm{\mathrm{l}}$; flows at pipes, pumps, and valves in
$\displaystyle\bm{\mathrm{z}}$; and relative speed of pumps
$\displaystyle\bm{\mathrm{s}}$. The vector
$\displaystyle\bm{\Omega}\in\mathbb{R}^{n_{\mathrm{J}}}$ encapsulates the
junctions demands which are considered predetermined in our study.

[Hyd-NLDAE]

$\displaystyle\displaystyle\bm{\mathrm{w}}(t+\Delta t_{\mathrm{H}})$
$\displaystyle\displaystyle={\bm{A}_{\mathrm{H}}}\bm{\mathrm{w}}(t)+{\bm{B}_{\mathrm{H}}}\bm{\mathrm{z}}(t),\leavevmode\nobreak\
$ (6a) $\displaystyle\displaystyle\bm{0}_{n_{\mathrm{J}}}$
$\displaystyle\displaystyle=\bm{E}_{z}\bm{\mathrm{z}}(t)+\bm{E}_{\Omega}{\bm{\Omega}}(t),\leavevmode\nobreak\
$ (6b)
$\displaystyle\displaystyle\bm{0}_{n_{\mathrm{P}}+n_{\mathrm{M}}+n_{\mathrm{V}}}$
$\displaystyle\displaystyle=\bm{E}_{w}\bm{\mathrm{w}}(t)+\bm{E}_{l}\bm{\mathrm{l}}(t)+{\bm{\Psi}}(\bm{\mathrm{z}},\bm{\mathrm{s}}),\leavevmode\nobreak\
$ (6c)

where $\displaystyle\\{\bm{A},\bm{B},\bm{E}\\}_{\bullet}$ are constant
matrices that depend on the network’s topology, components’ characteristics,
and hydraulic parameters. Additionally, $\displaystyle\bm{\Psi}(\cdot)$
gathers the nonlinear terms in (3), (4), and (5), and
$\displaystyle\bm{0}_{n}$ is a zero-vector of size $\displaystyle n$.

###  II-B Water Network Quality Modeling

This model traces the disinfectant concentrations throughout the network’s
components. The evolution of chlorine concentrations follows the conservation
of chemical’s mass, transport, and single-species reaction and decay models.
In each component, we represent the chlorine concentration as $\displaystyle
c$ with a superscript of the component symbol. Additionally, we accentuate the
hydraulics variables (e.g., velocities, flows, volumes, etc.) by coloring them
in violet whenever they appear in the WQ model.

####  II-B1 Conservation of Chlorine Mass in Nodes

For reservoirs, tanks, and junctions, the principles of conservation of mass
are applied. Reservoirs are considered a continuous source of chlorine with
constant concentrations over time; for a Reservoir $\displaystyle i$,
$\displaystyle c^{\mathrm{R}}_{i}(t+\Delta
t_{\mathrm{WQ}})=c^{\mathrm{R}}_{i}(t)$. On the other hand, junctions and
tanks are assumed to have complete immediate mixing at place with no storage
for junction and changing volume with time for tanks[27, 28]. Accordingly, if
node $\displaystyle i$ is defined as a Junction, chlorine concentration at
this node is calculated as

$\leavevmode\nobreak\ c_{i}^{\mathrm{J}}(t)=\frac{\sum_{j\in
L_{\mathrm{in}}}{\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}q_{\mathrm{in}}^{j}(t)}c_{\mathrm{in}}^{j}(t)+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}q^{\mathrm{B}}_{i}(t)}c^{\mathrm{B}}_{i}(t)}{{\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}q^{\mathrm{D}_{\mathrm{J}}}_{i}(t)}+\sum_{k\in
L_{\mathrm{out}}}{\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}q_{\mathrm{out}}^{k}(t)}}.$
(7)

While, if it is defined as a Tank, $\displaystyle c_{i}^{\mathrm{TK}}$ is
calculated as follows

$\small\begin{split}&{\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}V_{i}^{\mathrm{TK}}(t+\Delta
t_{\mathrm{WQ}})}c_{i}^{\mathrm{TK}}(t+\Delta
t_{\mathrm{WQ}})={\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}V_{i}^{\mathrm{TK}}(t)}c_{i}^{\mathrm{TK}}(t)\\\
&+\sum_{j\in
L_{\mathrm{in}}}{\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}q^{j}_{\mathrm{in}}(t)}c^{j}_{\mathrm{in}}(t)\Delta
t_{\mathrm{WQ}}+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}V^{\mathrm{B}}_{i}(t+\Delta
t_{\mathrm{WQ}})}c^{\mathrm{B}}_{i}(t+\Delta t_{\mathrm{WQ}})\\\ &-\sum_{k\in
L_{\mathrm{out}}}{\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}q^{k}_{\mathrm{out}}(t)}c_{i}^{\mathrm{TK}}(t)\Delta
t_{\mathrm{WQ}}+R^{\mathrm{TK}}(c_{i}^{\mathrm{TK}}(t)){\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}V_{i}^{\mathrm{TK}}(t)}\Delta
t_{\mathrm{WQ}},\end{split}$ (8)

where $\displaystyle c_{\mathrm{in}}^{j}(t)$ and $\displaystyle
c_{\mathrm{out}}^{k}(t)$ are the concentrations in the inflow and outflow
solute; $\displaystyle q^{\mathrm{B}}_{i}(t)$ and $\displaystyle
V^{\mathrm{B}}_{i}(t+\Delta t)$ are the flow and the volume of chlorine
injected to the node with concentration $\displaystyle c^{\mathrm{B}}_{i}(t)$
by booster station if located; and $\displaystyle
R^{\mathrm{TK}}(c^{\mathrm{P}}_{i}(x,t))$ is the decay and reaction expression
in tanks (refer to Section II-B3).

####  II-B2 Chlorine Transport and Reaction Model in Links

In the water quality model, pumps and valves are considered links with no
actual length and accordingly and accordingly there are no changes in the
chemical concentration from the upstream node. That is, for Pump
$\displaystyle i$ and Valve $\displaystyle k$ place after Junction
$\displaystyle j$, concentrations are expressed as $\displaystyle
c_{i}^{\mathrm{M}}(t+\Delta t_{\mathrm{WQ}})=c_{j}^{\mathrm{J}}(t+\Delta
t_{\mathrm{WQ}}),$ and $\displaystyle c_{k}^{\mathrm{V}}(t+\Delta
t_{\mathrm{WQ}})=c_{j}^{\mathrm{J}}(t+\Delta t_{\mathrm{WQ}}).$

Nonetheless, the transport and reaction model in pipes is simulated by the
one-dimension advection-reaction (1-D AR) partial differential equation (PDE),
which for Pipe $\displaystyle i$ is expressed as

$\partial_{t}c_{i}^{\mathrm{P}}=-{\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}v_{i}(t)}\partial_{x}c_{i}^{\mathrm{P}}+R^{\mathrm{P}}(c_{i}^{\mathrm{P}}(x,t)),$
(9)

where $\displaystyle c^{\mathrm{P}}_{i}(x,t)$ is the concentration in pipe at
location $\displaystyle x$ along its length and time $\displaystyle t$;
$\displaystyle v_{i}(t)$ is the mean flow velocity; and $\displaystyle
R^{\mathrm{P}}(c^{\mathrm{P}}_{i}(x,t))$ is the decay reaction expression in
pipes (more explanation is given in Section II-B3). Eq. (9) is discretized
over a fixed spatio-tamporal grid using the Explicit Upwind scheme—Eulerian
Finite-Difference based method [29, 30]. This scheme is conditionally stable
by satisfying the Courant-Friedrichs-Lewy condition (CFL). This condition puts
limits on the Courant number
$\displaystyle{\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}{\lambda}_{i}(t)}=\frac{{\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}v_{i}(t)}\Delta
t}{\Delta x_{i}}$ to be $\displaystyle
0<{\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}{\lambda}_{i}(t)}\leq
1$, for a Pipe $\displaystyle i$.

Subsequently, the Pipe $\displaystyle i$ with length $\displaystyle L_{i}$ is
split into a number of segments $\displaystyle
n_{s_{i}}=\Big{\lfloor}\frac{L_{i}}{{\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}v_{i}(t)}\Delta
t}\Big{\rfloor}$ of length $\displaystyle\Delta x_{i}=\frac{L_{i}}{s_{i}}$.
The chemical concentrations at the first segment and any other Segment
$\displaystyle s$ along the pipe’s length are calculated as follows with the
assumption that Junction $\displaystyle j$ is upstream this pipe

$\begin{bmatrix}c^{\mathrm{P}}_{i}(1,t+\Delta t_{\mathrm{WQ}})\\\
c^{\mathrm{P}}_{i}(s,t+\Delta
t_{\mathrm{WQ}})\end{bmatrix}=(1-{\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}{\lambda}_{i}(t)})\begin{bmatrix}c^{\mathrm{P}}_{i}(1,t)\\\
c^{\mathrm{P}}_{i}(s,t)\end{bmatrix}\\\
+{\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}{\lambda}_{i}(t)}\begin{bmatrix}c^{\mathrm{J}}_{j}(t)\\\
c^{\mathrm{P}}_{i}(s-1,t)\end{bmatrix}+\Delta
t_{\mathrm{WQ}}\begin{bmatrix}R^{\mathrm{P}}(c^{\mathrm{P}}_{i}(1,t))\\\
R^{\mathrm{P}}(c^{\mathrm{P}}_{i}(s,t))\end{bmatrix}.$ (10)

####  II-B3 Single-species Decay Model

The single-species decay model is a first-order model where the chlorine
concentrations are decaying due to wall reaction dynamics in pipes and bulk
reaction dynamics in both; pipes and tanks. Henceforward, the chlorine decay
reaction rates for Pipe $\displaystyle i$ and Tank $\displaystyle j$ are
$\displaystyle
k_{i}^{\mathrm{P}}=k_{b}+\frac{2k_{w}k_{f}}{r_{\mathrm{P}_{i}}(k_{w}+k_{f})},\,\,\,\,\,k_{j}^{\mathrm{TK}}=k_{b}$,
where $\displaystyle k_{b}$ is the bulk reaction rate constant; $\displaystyle
k_{w}$ is the wall reaction rate constant; $\displaystyle k_{f}$ is the mass
transfer coefficient between the bulk flow and the pipe wall; and
$\displaystyle r_{\mathrm{P}_{i}}$ is the pipe radius. Eventually, the decay
and reaction expressions for Segment $\displaystyle s$ of Pipe $\displaystyle
i$ and Tank $\displaystyle j$ are

$R^{\mathrm{P}}(c^{\mathrm{P}}_{i}(s,t))=-k_{i}^{\mathrm{P}}c^{\mathrm{P}}_{i}(s,t),\;\;R^{\mathrm{TK}}(c_{i}^{\mathrm{TK}}(t))=-k_{i}^{\mathrm{TK}}c^{\mathrm{TK}}_{i}(t).$
(11)

####  II-B4 Water Quality State-space Representation

The water quality model described in Sections II-B1, II-B2, and II-B3 can be
formulated as the following linear difference equations (LDE): [WQ-LDE]

$\displaystyle\displaystyle\leavevmode\nobreak\ \bm{\mathrm{x}}(t+\Delta
t_{\mathrm{WQ}})$
$\displaystyle\displaystyle={\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}\bm{A}_{\mathrm{WQ}}(t)}\bm{\mathrm{x}}(t)+{\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}\bm{B}_{\mathrm{WQ}}(t)}\bm{\mathrm{u}}(t),$
(12) $\displaystyle\displaystyle\bm{\mathrm{y}}(t)$
$\displaystyle\displaystyle=\bm{C}_{\mathrm{WQ}}(t)\bm{\mathrm{x}}(t),$ (13)

where vector
$\displaystyle\bm{\mathrm{x}}(t):=\\{\bm{\mathrm{c}}^{\mathrm{R}}(t),\;\bm{\mathrm{c}}^{\mathrm{J}}(t),\;\bm{\mathrm{c}}^{\mathrm{TK}}(t),\;\bm{\mathrm{c}}^{\mathrm{P}}(t),\;\bm{\mathrm{c}}^{\mathrm{M}}(t),\;$
$\displaystyle\bm{\mathrm{c}}^{\mathrm{V}}(t)\\}\in\mathbb{R}^{n_{x}}$ depicts
the concentrations of chlorine in the entire network and the total number of
states $\displaystyle
n_{x}=n_{\mathrm{R}}+n_{\mathrm{J}}+n_{\mathrm{TK}}+\sum_{i=1}^{n_{\mathrm{P}}}n_{s_{i}}+n_{\mathrm{M}}+n_{\mathrm{V}}$;
vector $\displaystyle\bm{\mathrm{u}}(t)\in\mathbb{R}^{n_{u}}$ represents the
dosages of injected chlorine; vector
$\displaystyle\bm{\mathrm{y}}(t)\in\mathbb{R}^{n_{y}}$ denotes the sensor
measurements of chlorine concentrations at specific locations in the network.
The state-space matrices
$\displaystyle\\{\bm{A}_{\mathrm{WQ}},\bm{B}_{\mathrm{WQ}},\bm{C}_{\mathrm{WQ}}\\}$
are all time-varying matrices that depend on the network topology and
parameters, hydraulic parameters, decay rate coefficients for the
disinfectant, and booster stations and sensors locations. It is customary to
assume that these matrices evolve at a slower pace than the states
$\displaystyle\bm{\mathrm{x}}(t)$ and control inputs
$\displaystyle\bm{\mathrm{u}}(t)$.

##  III Water Quality Controllability Proxies

In this section we introduce the notions of controllability for linear water
quality dynamics. We survey controllability measures that allow us to
qualitatively and quantitatively evaluate the controllability of the developed
water quality system (12). These metrics are assayed based on their
applicability to the water quality model and their suitability to the optimal
pump scheduling problem formulation. We consider the notion of water quality
system controllability and relate it to the hydraulic pump scheduling problem
formulated in Section IV.

From a control-theoretic perspective, controllability is the ability to steer
a system from initial states $\displaystyle\bm{x}_{o}:=\bm{x}_{t(0)}$ to
$\displaystyle\bm{x}_{s}:=\bm{x}_{T_{s}}$ by some input
$\displaystyle\bm{u}(t)$ [31]. That is, the goal is to be able to steer
complex dynamical systems to a desired state or trajectory. Specifically, for
water quality control we want to maintain chlorine concentrations within
certain levels.

######  Definition 1.

A linear system [e.g., the water quality system expressed as (12)] is
controllable if for any finite time interval $\displaystyle[0,T_{s}]$ and for
any initial state $\displaystyle\bm{\mathrm{x}}_{o}\in\mathbb{R}^{n_{x}}$, the
initial state $\displaystyle\bm{\mathrm{x}}_{o}$ can be steered or driven to a
final state $\displaystyle\bm{\mathrm{x}}_{T_{s}}\in\mathbb{R}^{n_{x}}$ for
some input $\displaystyle\bm{\mathrm{u}}(t)$ under the specified time
interval.

That being said, the dynamic linear system (12) is said to be controllable if
only if the controllability matrix for $\displaystyle
N_{s}=\frac{T_{s}}{\Delta t_{\mathrm{WQ}}}$ time-steps given as

$\leavevmode\nobreak\
\mathcal{C}_{N_{s}}:=\\{\begin{matrix}\bm{B}_{\mathrm{WQ}},\;\;\bm{A}_{\mathrm{WQ}}\bm{B}_{\mathrm{WQ}},\;\;\bm{A}_{\mathrm{WQ}}^{2}\bm{B}_{\mathrm{WQ}},\end{matrix}\\\
\begin{matrix}\dots,\;\;\bm{A}_{\mathrm{WQ}}^{N_{s}-1}\bm{B}_{\mathrm{WQ}}\end{matrix}\\}\in\mathbb{R}^{n_{x}\times
N_{s}n_{u}},$ (14)

is full row rank, i.e, $\displaystyle\mathrm{rank}(\mathcal{C}_{N_{s}})=n_{x}$
[32], without loss of generality as we assume that $\displaystyle
N_{s}n_{u}>n_{x}$. This is known as Kalman’s rank condition [31]. However,
matrix rank is a generic property that might lead to similar values depending
on the relations between the variables; it therefore is informative in a
qualitative sense but fails to indicate how controllable the system is under
many cases and various scenarios.

For the WQ dynamics (12), full row rank of $\displaystyle\mathcal{C}_{N_{s}}$
seldom occurs—this is due to the complexity and the high dimensionality of the
system. With that in mind, it is more practical to consider quantitative
measures of controllability which, unlike the aforementioned rank metric, are
able to reflect the difficulty in controlling the WQ system.

To that end, the notion of control energy
$\displaystyle\mathcal{E}(\bm{A}_{\mathrm{WQ}},\bm{B}_{\mathrm{WQ}},N_{s},\bm{x}_{o}\rightarrow\bm{x}_{T_{s}})$
is introduced to quantify the energy needed to steer the system from
$\displaystyle\bm{x}_{o}$ to $\displaystyle\bm{x}_{T_{s}}$ [33]. Ideally, we
want to minimize the energy required to control the system. The concept of
energy-related control depends on the application. For the case of WQ control,
it is related to the amount of chlorine needed to be injected into the system
to keep a desired chlorine level at the network’s components.

Metrics related to the input energy required are based on the controllability
Gramian
$\displaystyle\bm{W}_{c}(\bm{A}_{\mathrm{WQ}},\bm{B}_{\mathrm{WQ}},N_{s}):=\bm{W}_{c}\in\mathbb{R}^{n_{x}}$
that is defined for $\displaystyle N_{s}$ sum of matrices pair
$\displaystyle\bm{A}_{\mathrm{WQ}}$ and $\displaystyle\bm{B}_{\mathrm{WQ}}$ as

$\leavevmode\nobreak\
\bm{W}_{c}:=\sum_{\tau=0}^{N_{s}-1}\bm{A}_{\mathrm{WQ}}^{\tau}\bm{B}_{\mathrm{WQ}}\bm{B}_{\mathrm{WQ}}^{\top}(\bm{A}_{\mathrm{WQ}}^{\top})^{\tau}=\mathcal{C}_{N_{s}}\mathcal{C}_{N_{s}}^{\top},$
(15)

where the controllability Gramian $\displaystyle\bm{W}_{c}$, that is a
positive semidefinite metrics, provides an energy-related quantification of
controllability such that,
$\displaystyle\mathcal{E}\propto\mathrm{trace}(\bm{W}_{c})^{-1}$. We note here
that $\displaystyle\bm{W}_{c}$ is non-singular if the system is controllable
after time $\displaystyle{T_{s}}$, otherwise it is uncontrollable. In the
literature [34, 35] a myriad of measures exist; these measures provide a
scalar energy-related quantification of the controllability Gramian
$\displaystyle\bm{W}_{c}$. These measures include:
$\displaystyle\mathrm{log}\mathrm{det}(\bm{W}_{c})$,
$\displaystyle\mathrm{trace}(\bm{W}_{c})$,
$\displaystyle\mathrm{rank}(\bm{W}_{c})$, and minimum eigenvalue
$\displaystyle\lambda_{\min}(\bm{W}_{c})$. A discussion on the aforementioned
measures is given as follows.

  1. 1.

The $\displaystyle\mathrm{log}\mathrm{det}(\bm{W}_{c})$ metric is proportional
to the volumetric measure of the ellipsoid enclosing the set of states that
can be reached with at most a unit control energy input.

  2. 2.

The $\displaystyle\mathrm{trace}(\bm{W}_{c})$ metric is inversely related to
the average controllability energy in all directions of the state-space.

  3. 3.

The $\displaystyle\mathrm{rank}(\bm{W}_{c})$ metric quantifies the size of the
controllable subspace.

  4. 4.

The $\displaystyle\lambda_{\min}(\bm{W}_{c})$ metric is inversely related to
the control energy in the most difficult control direction. The smallest
eigenvalue quantifies that the worst-case direction that requires the largest
amount of control energy.

It is worth mentioning that the notion of controllable subspace is equivalent
to the notion of reachable subspace, this is related to the representation of
the Gramian and its associated metrics [36]. A system is said to be reachable
in a specific state space if the subspace of all the reachable states from an
initial state $\displaystyle\bm{\mathrm{x}}_{o}$ is equal to the whole state
space. Attaining this propriety is important while controlling complex water
quality dynamics. That is, our goal is to preserve and maintain a certain
level of energy, within a controllable and reachable subspace.

For an uncontrollable WQ system with total number of states $\displaystyle
n_{x}$, let the rank of the controllability matrix/Gramian be $\displaystyle
k<n_{x}$. Then there exists a nonsingular matrix
$\displaystyle\bm{T}\in\mathbb{R}^{n_{x}\times n_{x}}$ such that

$\leavevmode\nobreak\
\begin{split}&\bar{\bm{A}}_{\mathrm{WQ}}=\bm{T}\bm{A}_{\mathrm{WQ}}\bm{T}^{-1}=\begin{bmatrix}\bar{\bm{A}}_{\mathrm{WQ},11}&\bar{\bm{A}}_{\mathrm{WQ},12}\\\
\bm{0}&\bar{\bm{A}}_{\mathrm{WQ},22}\end{bmatrix},\\\
&\bar{\bm{B}}_{\mathrm{WQ}}=\bm{T}\bm{B}_{\mathrm{WQ}}=\begin{bmatrix}{\bar{\bm{B}}_{\mathrm{WQ},1}}\\\
\bm{0}\end{bmatrix},\end{split}$ (16)

where
$\displaystyle\bar{\bm{A}}_{\mathrm{WQ},11},\;\bar{\bm{A}}_{\mathrm{WQ},12}$
and $\displaystyle\bar{\bm{A}}_{\mathrm{WQ},22}$ have dimensions of
$\displaystyle k\times k,\;k\times(n_{x}-k)$ and
$\displaystyle(n_{x}-k)\times(n_{x}-k)$, and
$\displaystyle\bar{\bm{B}}_{\mathrm{WQ},1}$ has $\displaystyle k$ rows.
Namely, $\displaystyle\bar{\bm{A}}_{\mathrm{WQ},11}$ and
$\displaystyle\bar{\bm{B}}_{\mathrm{WQ},1}$ define a controllable subspace.
Defining the controllable subspace is instrumental in measuring metrics like
$\displaystyle\mathrm{trace}$, $\displaystyle\mathrm{logdet}$, and
$\displaystyle\lambda_{\min}$ when the system is not full rank. Readers are
referred to [37] for the theorem pertaining to the development of reachable
subspace decomposition.

To that end, we investigate the use of the aforementioned metrics on the real-
time control and pump scheduling optimization problem to achieve the stated
goal of preserving energy within specific subspaces of the system.

##  IV Hydraulics Control

In this section, we build the decoupled and coupled pump scheduling problems.

###  IV-A Decoupled Pump Scheduling Problem

We first follow the approximation approach for the system components proposed
in [38] to formulate the pump control problem. The system components include
the head loss in pipes, head loss in valves, head gain in pumps, and pump
power consumption. For the head losses in pipes and valves, we apply piecewise
linear approximations that transform (3) and (5) to multiple linear
constraints.

Specifically, each pipe/valve’s head loss curve is segmented into linear
segments, which are determined by connecting points that can be calculated
offline as a pre-optimization step. For the curve and piecewise linearization
presented in Fig. 1a as an example, four points are located and the segmented
lines are connected. Using the equations for these lines, three constraints
are added to the optimization problem for this specific pipe. When aiming for
closer fits (i.e., employing more segmentation), a drawback emerges in the
context of large network models, which is that scalability is negatively
affected. For formulations of $\displaystyle N_{\mathrm{PW}}$ pieces, the
constraints for Pipe $\displaystyle i$ connecting and flowing water between
node $\displaystyle j$ and node $\displaystyle k$, are as follows

$\displaystyle\displaystyle
h_{j}(t)-h_{k}(t)-\sum_{n=1}^{N_{\mathrm{PW}}}\tilde{m}_{n}\zeta_{n}\zeta_{n}(t)-\sum_{n=1}^{N_{\mathrm{PW}}}\tilde{b}\omega_{n}\omega_{n}(t)=0,\leavevmode\nobreak\
$ (17a) $\displaystyle\displaystyle
q^{\mathrm{P}}_{i}(t)-\sum_{n=1}^{N_{\mathrm{PW}}}\zeta_{n}\zeta_{n}(t)=0,\leavevmode\nobreak\
$ (17b)
$\displaystyle\displaystyle\sum_{n=1}^{N_{\mathrm{PW}}}\omega_{n}\omega_{n}(t)=1,\leavevmode\nobreak\
$ (17c)
$\displaystyle\displaystyle\begin{cases}-\zeta_{n}\zeta_{n}(t)+q_{n,\min}\omega_{n}\omega_{n}(t)\leq
0,\\\ \zeta_{n}\zeta_{n}(t)-q_{n,\max}\omega_{n}\omega_{n}(t)\leq
0,\end{cases}\leavevmode\nobreak\ $ (17d)

where $\displaystyle\tilde{m}_{n}$ and $\displaystyle\tilde{b}$ are the
$\displaystyle n$ segment’s line slope and intercept, while $\displaystyle
q_{n,\min}$ and $\displaystyle q_{n,\max}$ are the flow boundary limits. In
addition, for the same segment $\displaystyle n$, $\displaystyle\zeta_{n}(t)$
and the binary $\displaystyle\omega_{n}(t)$ are decision variables to enable
falling within the right segment range. Constraint (17a) represents the
linearized head loss through the pipe within this segment, (17b) defines the
pipe segments flow equality constraint, (17c) allows the segment selection,
and (17d) are the boundary constraints for each segment. By adding these
equality and inequality constraints to the optimization problem, new integer
decision variables are introduced; mixed integer programming.

Moreover, in contrast to study [38] we use variable speed pumps, making the
problem more general but a bit more complex. This complexity emerges as a
result of the pump curve variation based on the selected pump speed; see Fig
1b. Therefore, we follow the proposed methodology by [39, 40], making minor
modifications to transform the decision variables into the flow through the
pump and pump speed, all while ensuring convexity. This transformation is
achieved by applying the affinity rules that govern the relationship between
pump shaft speed, discharge, and head gain. These rules relate the discharge
and head gain to the shaft speed via parabolic relations. Specifically, pump
discharge exhibits a linear relationship with pump speed, whereas head gain
depends on the square of the speed. As a pre-optimization step, both the pump
performance curve and the corresponding power consumption that we aim to
minimize are approximated to formulate two convex expressions to be integrated
in the optimization problem.

First, we approximate the characteristic curve of Pump $\displaystyle i$ from
Eq. (4) to the following

$\leavevmode\nobreak\ \Delta
h_{i}^{M}\big{|}_{\mathrm{App}}=\beta_{1}\Big{(}q^{\mathrm{M}}_{i}\Big{)}^{2}+\beta_{2}q^{\mathrm{M}}_{i}+\beta_{3}\Big{(}s^{\mathrm{M}}_{i}\Big{)}^{2}+\beta_{4},$
(18)

where $\displaystyle\beta_{1},\beta_{2},\beta_{3}$ and
$\displaystyle\beta_{4}$ are coefficient calculated by minimizing the error
between $\displaystyle\Delta h_{i}^{M}$ in Eq. (4) and Eq. (18) with
$\displaystyle\beta_{1},\beta_{3}\geq 0$ to ensure convexity.

(a)  (b)  (c)

Figure 1: (a) Linear, (b) variable-speed pump curve, and (c) the resultant
power consumption ($\displaystyle h_{0}=393.7,\alpha=3.7\times
10^{-6},\nu=2.59$).

Furthermore, alongside the system dynamics and their approximations addressed
in the preceding sections of the paper, we also account for the physical
constraints on the head levels and flow values among the various network
components—expressed in Eq. (19). Pump speed vector
$\displaystyle\bm{\mathrm{s}}(t)$ is constrained to be between
$\displaystyle\bm{0}$ and $\displaystyle\bm{\mathrm{s}}_{\max}$, where the
zero values indicate that the pump is off. These considerations can all be
written as box constraints

$\leavevmode\nobreak\
\begin{split}\bm{\mathrm{w}}(t)\in[\bm{\mathrm{w}}_{\min},&\bm{\mathrm{w}}_{\max}],\;\;\bm{\mathrm{l}}(t)\in[\bm{\mathrm{l}}_{\min},\bm{\mathrm{l}}_{\max}],\\\
\bm{\mathrm{z}}(t)\in[\bm{\mathrm{z}}_{\min},&\bm{\mathrm{z}}_{\max}],\;\;\bm{\mathrm{s}}(t)\in[\bm{0},\bm{\mathrm{s}}_{\max}].\end{split}$
(19)

Lastly, the objective function for this problem enforces minimizing the cost
of power consumption by pumps. This objective function is expressed as

$\Pi(t)=\varphi_{\mathrm{EL}}\sum_{i}^{n_{\mathrm{M}}}\frac{\rho_{\mathrm{W}}g}{\eta_{i}(t)}\Delta
h_{i}^{\mathrm{M}}(t)q_{i}^{\mathrm{M}}(t),$ (20)

where $\displaystyle\varphi_{\mathrm{EL}}$ is the price of each kW of
electricity per hour ($/kWh), $\displaystyle\rho_{\mathrm{W}}$ is water
density, $\displaystyle g$ is the gravity acceleration,
$\displaystyle\eta_{i}(t)$ is the efficiency of Pump $\displaystyle i$ under a
head gain of $\displaystyle\Delta h_{i}^{\mathrm{M}}(t)$ and flow of
$\displaystyle q_{i}^{\mathrm{M}}(t)$.

As illustrated in Fig. 1c, the power consumption of pumps exhibits a nonlinear
relationship with respect to the flow rate, and this relationship shifts when
the pump speed is altered. For a specific pump speed, this objective function
can be approximated to a convex second-order function [40]. However, when the
actual pump speed deviates from the speed at which the curve is approximated,
the optimized pump power consumption varies substantially from the actual
power consumption—a limitation of [40] to work in a vacancy of the
approximation region. To overcome this issue, we approximate the objective
function to be function in the pump speed and flow. This approach allows us to
avoid constructing an approximate formulation reliant on head gain with
inherent error, thereby preventing an increase in inaccuracies. We define this
function for Pump $\displaystyle i$ to be

$\leavevmode\nobreak\
\begin{split}\Pi_{\mathrm{App}}(t)&=\theta_{1}+\theta_{2}q^{\mathrm{M}}_{i}+\theta_{3}\Big{(}q^{\mathrm{M}}_{i}\Big{)}^{2}+\theta_{4}s^{\mathrm{M}}_{i}\\\
&+\theta_{5}\Big{(}s^{\mathrm{M}}_{i}\Big{)}^{2}+\theta_{6}q^{\mathrm{M}}_{i}s^{\mathrm{M}}_{i}.\end{split}$
(21)

The vector $\displaystyle\bm{\theta}:=\\{\theta_{i}|i\in\\{1,\cdots,6\\}\\}$
collects the coefficients of the approximate second-order power consumption
function. These coefficients are derived by solving a straightforward
optimization problem aimed at ensuring the convexity of this function. This
function is convex under the condition that its Hessian is positive
semidefinite. The Hessian is defined to be
$\displaystyle\mathbfcal{H}=\begin{bmatrix}2\theta_{3}&\theta_{6}\\\
\theta_{6}&2\theta_{5}\end{bmatrix}$. That is, we obtain these coefficients by
solving the following optimization problem after pre-calculating the power
consumption for different $\displaystyle N_{\mathrm{M,OP}}$ operating points
in the $\displaystyle(q^{\mathrm{M}}_{i},s^{\mathrm{M}}_{i})$ domain.

$\displaystyle\displaystyle\underset{\bm{\theta}}{\mbox{minimize}}\;\;$
$\displaystyle\displaystyle\sum_{j=1}^{N_{\mathrm{M,OP}}}(\Pi_{\mathrm{App}}(q^{\mathrm{M}}_{i,j},s^{\mathrm{M}}_{i,j})-\Pi(q^{\mathrm{M}}_{i,j},s^{\mathrm{M}}_{i,j}))^{2}$
(22a) $\displaystyle\displaystyle\begin{split}\leavevmode\nobreak\
\mbox{subject to}\;\;&\mathbfcal{H}\succeq 0.\end{split}$ (22b)

Eventually, the decoupled pump scheduling optimization problem is formulated
to minimize the pump power consumption while being subjected to the system
dynamics and functional constraints. The modifications made to the pipe and
pump dynamics to ensure a convex formulation lead to the elimination of the
nonlinear formulation in (6); Eq. (6c). Henceforward, we incorporate these
modification in addition to the linear expressions in (6), referring to this
combination as Hyd-LDAE. The decision variables are collected in vector
$\displaystyle\bm{\Upsilon}(t):=\\{\bm{\mathrm{w}}(t),\bm{\mathrm{l}}(t),\bm{\mathrm{z}}(t),\bm{\zeta}(t),\bm{\omega}(t),\bm{\mathrm{s}}(t)\\}$.
Therefore, the final optimization problem is represented as mixed-integer
quadratic constrained quadratic problem (MIQCQP) for each hydraulic time-step:

$\displaystyle\displaystyle\underset{\bm{\Upsilon}(t)}{\mbox{minimize}}\hskip
28.45274pt$
$\displaystyle\displaystyle\Pi_{\mathrm{App}}(\bm{\Upsilon}(t))\Delta
t_{\mathrm{H}}$ (23a)
$\displaystyle\displaystyle\begin{split}\leavevmode\nobreak\ \mbox{subject
to}\hskip 28.45274pt&\mbox{{Hyd-
LDAE}}\;\eqref{equ:HydSSa},\eqref{equ:HydSSb},\\\
&\eqref{equ:PipeConstraints},\;\eqref{equ:AppPumpCurve},\;\eqref{equ:HydPhysConstraints}.\end{split}$
(23b)

###  IV-B Coupled WQ-Aware Pump Control

In this paper, we include the discussed water quality controllability Gramian
(WQ-CG) metrics in the previously developed pump scheduling problem to
formulate the coupled WQ controllability-guided pump scheduling optimization
problem. The dependency of those metrics calculations on the hydraulics
settings is direct, yet complex and lead to highly nonlinear
expressions—Appendix A showcases the raised concern for the case of the Three-
node network. In addition, multiple factors should be taken into consideration
that reflect the general and specific characteristics of the operation of
WDNs. These factors impose logical and physical constraints on the choice and
purpose of controllability metric integrated in our problem. For instance,
during the initial simulation time with zero initial chlorine concentrations
throughout the network, we aim for a higher controllability Gramian rank along
with high energy. This approach also applies to branched networks with
numerous dead-ends to ensure chlorine concentrations remain within standard
limits. Another scenario arises when we need to store water with sufficient
chlorine concentrations in tanks during off-peak demand periods, ready for
distribution when the tanks are in demand, supplying either specific network
sections or the entire network.

Formulating and solving a problem that takes into account these factors and
achieves the desired level of WQ controllability presents several challenges.
These challenges primarily arise from the high nonlinearity and complexity
involved in defining these metrics. Additionally, the formulation of the WQ
state-space matrices, and consequently, its controllability Gramian, depends
on factors such as flow directions in each pipe and the number of segments
into which it is discretized. Yet, these factors are to be determined by
solving the problem itself. In response to these challenges, we propose the
approach we have developed to address these issues effectively.

First, we overcome the flow directions issue while formulating the
$\displaystyle\bm{A}_{\mathrm{WQ}}$ matrix by building it for both cases and
use the introduced binary variables $\displaystyle\omega(t)$ for each pipe.
These variables define which pipe piecewise-linearization segment is chosen
and accordingly the flow direction. For each element of the matrix depending
on the flow direction for a specific pipe, it is multiplied by the summation
of half of the $\displaystyle\omega(t)$ variables representing the
corresponding direction.

In addition, as explained in Section II, the number of segments defined for
each pipe depends on the water velocities (i.e., hydraulics in the system and
decision variable of our problem). This number of segments defines the water
quality model dimensions for the simulation window and accordingly the
dimensions of the $\displaystyle\bm{A}_{\mathrm{WQ}}$ and
$\displaystyle\bm{B}_{\mathrm{WQ}}$ matrices. Yet, the hydraulics variables
are to be determined by the problem through which we aim to account for the
water quality controllability. To that end, we define the water quality system
dimensions offline as a prior-control step that preserves WQ model stability
and hydraulics applicable scenarios. We randomly generate pump speeds between
0 and 1 and solve the systems’ flow and heads according to variety of demands
patterns. Then, we calculate the number of segments needed for each pipe to
ensure a fulfilled stability condition. Finally, we define that number to be
the minimum for each pipe out of all the scenarios. This approach guarantees
that after solving the pump scheduling problem and obtaining actual
operational hydraulic setting, the water quality model is stable and has been
well-represented.

Last, given that our optimization problem is tailored for the specific purpose
of enhancing dynamics within WDNs, we simplify the utilization of WQ
controllability metrics within this context. This simplification is based on
the characteristics of the dynamics inherent to these systems. These
distinctive dynamics include: (i) booster stations are located on nodes only,
(ii) as each pipe is discretized into number of segments for WQ dynamics
simulation then ensuring controllability should be over the whole length of
it, and (iii) in many scenarios, we aim for higher controllability coverage
and/or energy to reach specific junctions/dead-ends that serve large areas
and/or tanks scheduled for on-demand operation during various simulation
intervals. To that end, we adopt the concept of target controllability [41]
given in Definition 2. Target controllability allows us to choose the desired
target nodes and accordingly eliminates the large dimentionality issue
associate with the WQ representation. In this case, the metrics are applied to
the targeted controllability Gramian
$\displaystyle\bm{W}_{\mathcal{T}}=C_{\mathcal{T}}\bm{W}_{c}C_{\mathcal{T}}^{\top}$.

######  Definition 2.

A discrete linear system is said to be target controllable with respect to the
target set $\displaystyle\mathcal{T}\subseteq\mathcal{G}$;
$\displaystyle|\mathcal{T}|=n_{y_{f}}$ over time $\displaystyle[0,T_{s}]$, if
for any final output $\displaystyle y_{f}(t)=C_{\mathcal{T}}x(t),$
$\displaystyle y_{f}\in\mathbb{R}^{n_{y_{f}}}$ and
$\displaystyle\bm{\mathrm{x}}_{o}\in\mathbb{R}^{n_{x}}$, the initial state
$\displaystyle\bm{\mathrm{x}}_{o}$ can be steered or driven to final state of
the target nodes as $\displaystyle y_{f}$ for some input
$\displaystyle\bm{\mathrm{u}}(t)$ under the specified time interval. The
output matrix $\displaystyle C_{\mathcal{T}}$ identifies the set of target
nodes $\displaystyle\mathcal{T}$.

By integrating constraints on the target controllability Gramian to Eq. (23),
the coupled optimization problem is formulated in (24) as nonconvex nonlinear
problem.

$\displaystyle\displaystyle\underset{\bm{\Upsilon}(t)}{\mbox{minimize}}\;\;\;$
$\displaystyle\displaystyle\Pi_{\mathrm{App}}(\bm{\Upsilon}(t))\Delta
t_{\mathrm{H}}-\Theta_{1}{\rm{trace}}(\bm{W}_{\mathcal{T}}(\bm{\Upsilon}(t)))$
$\displaystyle\displaystyle-\Theta_{2}{\lambda_{\min}}(\bm{W}_{\mathcal{T}}(\bm{\Upsilon}(t))),$
(24a) subject to $\displaystyle\displaystyle\mbox{{Hyd-
LDAE}}\;\eqref{equ:HydSSa},\eqref{equ:HydSSb},\eqref{equ:PipeConstraints},\eqref{equ:AppPumpCurve},\eqref{equ:HydPhysConstraints},$
$\displaystyle\displaystyle{\rm{rank}}(\bm{W}_{\mathcal{T}}(\bm{\Upsilon}(t)))\geq
n_{y_{f}}-1,$ (24b)

where $\displaystyle\Theta_{1}$ and $\displaystyle\Theta_{2}$ are scaling
factors.

We choose the target nodes to be the one listed in point (iii) along with the
first and last segment of the connecting pipes to assure the validation of the
flow directions. However, while employing the targeted controllability Gramian
helps reduce the number of constraints, integrating the metrics results in
nonlinear formulations. Nonlinear solvers can address these issues, but they
demand significantly more runtime, which escalates exponentially as network
size increases. Additionally, the presence of the "$\displaystyle\rm{rank}$"
constraint imposes limitations on this formulation, as many solvers do not
support constraints related to rank. That is, we employ the following
simplifications to accelerate the computational process while achieving the
desired output. Some of these simplification strategies are suitable to
relatively small systems, while others are tailored for larger networks. In
addition, some preparation steps can be taken for all sizes.

Figure 2: Proposed Quality-Aware Hydraulic Control Framework Flowchart.

First, we approximate the Gramian by eliminating the denominator in all values
of $\displaystyle\bm{A}_{\mathrm{WQ}}$ and $\displaystyle\bm{B}_{\mathrm{WQ}}$
to only result in polynomial expressions in the Gramian and change its
notation to $\displaystyle\bm{\widetilde{W}}_{\mathcal{T}}$. This helps the
distinguish process between the columns/rows for the $\displaystyle\rm{rank}$
determination and also the increase in the $\displaystyle\rm{trace}$ and
$\displaystyle\lambda_{\min}$. Second, depending on the scenario under
consideration, the optimization problem (24) can be modified to be rank-
oriented to achieve desired controllability coverage or/and energy-oriented to
achieve desired controllability energy. To solve the $\displaystyle\rm{rank}$
constraint issue, we employ the approach detailed in [42], which involves
applying a convex relaxation of the rank constraint by using a nuclear norm
penalty and specifying the required rank—formulated as Eq. (25). On the other
hand, the energy-oriented problem (Eq. (26)) is formulated to avoid the
transformation needed to maximize the $\displaystyle\lambda_{\min}$ by
maximizing the $\displaystyle\rm{trace}$ of the target controllability Gramian
built for specific direction.

$\displaystyle\displaystyle\underset{\bm{\Upsilon}(t)}{\mbox{minimize}}\;\;\;$
$\displaystyle\displaystyle\Pi_{\mathrm{App}}(\bm{\Upsilon}(t))\Delta
t_{\mathrm{H}}-\Theta_{3}||\bm{\widetilde{W}}_{\mathcal{T}}(\bm{\Upsilon}(t))||_{*},$
(25a) subject to $\displaystyle\displaystyle\mbox{{Hyd-
LDAE}}\;\eqref{equ:HydSSa},\eqref{equ:HydSSb},\eqref{equ:PipeConstraints},\eqref{equ:AppPumpCurve},\eqref{equ:HydPhysConstraints},$
$\displaystyle\displaystyle
l_{r}||(\bm{\widetilde{W}}_{\mathcal{T}}(\bm{\Upsilon}(t)))||_{2}-{\rm{trace}}(\bm{\widetilde{W}}_{\mathcal{T}}(\bm{\Upsilon}(t)))\leq
0,$ (25b)

where $\displaystyle\Theta_{3}$ is scaling factor,
$\displaystyle||\bm{\widetilde{W}}_{\mathcal{T}}||_{*}$ is the nuclear norm,
$\displaystyle||\bm{\widetilde{W}}_{\mathcal{T}}||_{2}$ is the second norm,
and $\displaystyle l_{r}$ is the desired rank of the target controllability
Gramian.

$\displaystyle\displaystyle\underset{\bm{\Upsilon}(t)}{\mbox{minimize}}\;\;\;$
$\displaystyle\displaystyle\Pi_{\mathrm{App}}(\bm{\Upsilon}(t))\Delta
t_{\mathrm{H}}-\sum_{i=1}^{n}\Theta_{i}{\rm{trace}}(\bm{\widetilde{W}}_{\mathcal{T}_{i}}(\bm{\Upsilon}(t)))$
(26a) subject to $\displaystyle\displaystyle\mbox{{Hyd-
LDAE}}\;\eqref{equ:HydSSa},\eqref{equ:HydSSb},\eqref{equ:PipeConstraints},\eqref{equ:AppPumpCurve},\eqref{equ:HydPhysConstraints},$
$\displaystyle\displaystyle{\rm{trace}}(\bm{\widetilde{W}}_{\mathcal{T}_{i}}(\bm{\Upsilon}(t)))\geq
0\;\;\;\forall i=1,\dots,n.$ (26b)

The proposed framework is suitable for networks of all sizes. However, in
scenarios where multiple controllability constraints need to be satisfied for
large networks, the computational time can become demanding. That is, for
these networks the controllability Gramian can be built for a defined
important path of the network.In such cases, we propose focusing on building
the controllability Gramian for specific important network paths. Depending on
network characteristics, this path can be determined based on factors like
high demands, water transfer between reservoirs and elevated tanks, or
mainline locations before branching occurs. To that end, our approach is
applicable for the entire or parts of the network.

It is worth mentioning that the final formulation is nonlinear and nonconvex,
yet tackled by several solvers and the level of complexity and accordingly the
runtime are determined by the scenario under consideration. Additionally, in
some cases constraining the optimization problem to reach specific level of
controllability results in infeasible problem. Therefore, we put condition on
our formulation to reduce that level for such case and re-solve the problem
successively and work on more improvement in the next time-step. Eventually,
Fig. 2 summarizes the flowchart of the proposed framework.

In last, after the quality-aware pump schedules are computed for the
simulation period, the WQ control problem can be solved by the appropriate
chosen technology. In our paper, we apply the model predictive control (MPC)
algorithm adopted in [43]. For brevity, we do not include the details and the
derivation of the control problem and we refer the reader to the cited study
reaching the final formulation. Note that, the WQ control problem is
constrained by upper and lower bounds on the chlorine concentrations at all
network components. These bounds are specified by EPA regulations to be
between $\displaystyle\bm{\mathrm{x}}_{\min}=0.2$ mg/L and
$\displaystyle\bm{\mathrm{x}}_{\max}=4$ mg/L [44].

Figure 3: Networks under study and their layouts: (a) Three-node network, (b)
Net1, and (c) Richmond skeleton network.

##  V Case Studies

In our study, we investigate the validity of our proposed control algorithms
on three networks; Three-node network, Net1, and the Richmond skeleton network
[45] (see Fig. 3 illustrating their layouts). The Three-node network consists
of a reservoir, a pump, a junction, a pipe, and a tank, in addition to one
booster station located at Junction J1 and one WQ sensor located at Tank TK1.
Net1 has a reservoir, a pump, a tank, 9 junctions, and 12 pipes. Two booster
stations are positioned at Junctions 1 and 6 of Net1, and two WQ sensors are
placed at Junctions 4 and 9 within the network. Lastly, the Richmond skeleton
network is a schematic representation of the Richmond water distribution
system, which is composed of one reservoir, 7 pumps, 41 junctions, 6 tanks, 8
valves, 37 pipes, 4 booster stations, and 3 WQ sensors.

Figure 4: Flow velocity pattern in P1 of the Three-node network in comparison
to the minimum required velocity to achieve maximum reachability of the flow
along the pipe’s length within the hydraulic time-step vs. the corresponding
percentage (%) of the system’s WQ-CG rank out of the # states in the subspace.
Change in bars colors to highlight (darker shade) the windows where the WQ-CG
is not full rank.

###  V-A Hydraulics Settings vs. Water Quality Controllability

Before we showcase the developed optimal pump scheduling approach, that is
augmented with a desired level of controllability of water quality states, we
first provide an investigation towards how varying the pump schedule actually
affects the water quality dynamics. First, we run a hydraulic simulation on
the Three-node network which has a simple layout where chlorine travels along
a single available path. This simulation results in the heads at TK1 and J1
and the head loss in P1 shown in Fig. 5 and the corresponding velocity profile
of P1 as demonstrated in Fig. 4 under the flow directions indicated in Fig. 3.
Next, we construct the WQ-CG and employ the rank metric to assess the system’s
performance. It is important to note that in our analysis, we exclude
Reservoir R1 and Pump M1 from the assessment as they are upstream of the
booster station located at J1. This exclusion allows us to focus on evaluating
the metric results specifically within the subspace of interest. We present
these results as a percentage of the rank calculated within every hydraulic
time-step, relative to the total number of states within this subspace
($\displaystyle n_{\mathrm{J}}+n_{s_{1}}+n_{\mathrm{TK}}$), number of segments
$\displaystyle n_{s_{1}}$ of P1 is 100 segments.

Figure 5: Head at Tank TK1 and Junction J1 (top) and the corresponding head
loss in Pipe P1 (bottom) of the Three-node network.

In this particular scenario, the hydraulic time-step is set to 1 hour, while
the WQ time-step is 10 seconds. Considering that Pipe P1 has a length of 1,000
ft, it is necessary for a water parcel to achieve a minimum velocity of 0.278
ft/sec in order to traverse the entire length of the pipe and reach Tank TK1
within the specified hydraulic time-step. This velocity directly influences
the change of chlorine concentrations over time and space as expressed in the
advection-reaction equation (Eq. (9)). It is obvious in Fig. 4 that this
characteristic has vital influence on the WQ controllability; when the
velocities through P1 surpass this velocity boundary, the WQ-CG exhibits full
rank, indicating full controllability of the system. On the other hand, in
most cases where the velocities are lower, the system is uncontrollable with
different deviations in the rank of the WQ-CG. It is worth mentioning that
other factors are affecting the response of the system states to the inputs
(in this case, the chlorine injections at J1). For instance, the rate of
change of water volume at TK1 and the flows rates in the system’s components,
which are directly related to the head levels at the nodes and the head
loss/gain at the links—refer to Fig. 5 and Appendix A. Note that, the rank is
calculated with keeping original MATLAB’s machine epsilon which is equal to
$\displaystyle 2.2204\mathrm{e}{-16}$. Accordingly, with relative difference
in that metric less than that machine epsilon, the two elements are considered
dependent.

Figure 6: Controllability metrics, $\displaystyle\mathrm{trace}$ and
$\displaystyle\lambda_{\min}$ of the controllable subspace Gramian
$\displaystyle(\bm{W}_{c})$ of the Three-node network vs. the percentage (%)
of the WQ-CG rank out of the # states. Change in bars colors to highlight
(darker shade) the windows where the WQ-CG is not full rank.

In addition, in Fig. 6 we showcase the
$\displaystyle\mathrm{trace}(\bm{W}_{c})$ and
$\displaystyle\lambda_{\min}(\bm{W}_{c})$ metrics of the controllable
subspaces of the system over the same simulation period of 24 hours. It is
noticeable that in many cases where the WQ-CG is full rank, there is change in
the $\displaystyle\mathrm{trace}$ and $\displaystyle\lambda_{\min}$ values
which reflects varying levels of energy stored in the system and their
respective directions. Furthermore, the greater the disparity between the
$\displaystyle\mathrm{trace}$ and $\displaystyle\lambda_{\min}$ values, the
more pronounced the sparsity in the distribution of energy.

Figure 7: Change in the volume of Tank 11 and the operational pump head gain
in Net1 under three different hydraulic scenarios; Hyd#1, Hyd#2, and Hyd#3.

Secondly, we apply three different hydraulic settings on Net1 resulting in the
change in tank volume and the operational pump head gain depicted in Fig. 7.
In the Hyd#1 scenario, Tank 11 is filling throughout the whole simulation
window. Whilst, Hyd#2 and Hyd#3 scenarios have smaller windows where Tank 11
is on demand. This results in different flow directions for specific pipes,
which directly influence the water quality dynamics. In addition, it results
in different total number of states (i.e, changes the number of segments into
which each pipe is divided). All these scenarios are run over a 24 hours
simulation period with a hydraulic time-step of 1 hour and a WQ time-step of
10 seconds. The total number of states for Hyd#1 is 490, 470 for Hyd#2, and
425 for Hyd#3.

(a)  (b)

Figure 8: For each of the three hydraulic scenarios (Hyd#1, Hyd#2, and Hyd#3)
applied on Net1, (a) the percentage (%) of the WQ-CG rank out of the # states
and (b) $\displaystyle\mathrm{trace}$ of the controllable subspace Gramian.

Fig. 8a shows the change in the rank of the controllability Gramian for each
of the scenarios. The controllability matrix/Gramian does not reach full rank
under the conditions where the booster stations are located at Junctions 1 and
6 and the water is drawn from Tank 11. In addition, the
$\displaystyle\mathrm{trace}$ metric for the controllable subspace is
calculated and the results are illustrated in Fig. 8b. These values are
affected by the direction in which the energy is distributed and the response
of the system states to the inputs. Such response is affected by the input
values in comparison to the systems states—the rate in which chlorine is
injected into the system and the flow rates in the system. Note that, the
results of the $\displaystyle\mathrm{log}\mathrm{det}$ metric for all
hydraulic simulations under consideration are equal to the values of the
$\displaystyle\mathrm{trace}$ presented in Fig. 8b.This observation implies
that either of these metrics can be effectively employed in our study. In
conclusion, each of the $\displaystyle\mathrm{rank}$,
$\displaystyle\mathrm{trace}$, and $\displaystyle\mathrm{\lambda_{\min}}$
reflects an important behavior of the WQ dynamics and can be taken into
consideration to reach the desired level of controllability over the system
tailored to the scenario under focus.

###  V-B WQ Controllability-Aware Optimal Pump Scheduling

In this section we showcase the results of solving the decoupled and coupled
pump optimal scheduling problems. The optimization problems are interfaced
using YALMIP in MATLAB R2023a and solved using Gurobi and/or BMIBNB
optimization solvers. The use of two optimization solvers is to compensate for
the difference in the underlying problems formulated for each network and each
scenario. The problems formulations differ in the their level of complexity
and nonlinearity order. That being said, Gurobi is utilized to solve problems
with lower nonlinearity order as it is a robust and fast solver capable of
handling binary decision variables. Since Gurobi is limited under a highly
nonlinear setting, the BMIBNB solver is used. BMIBNB is a global nonlinear
solver capable of handling nonconvex problems, however as with all global
solvers, the computational time is relatively more expensive in comparison to
Gurobi.

(a)  (a)

(b)  (b)

Figure 9: Optimal Tank TK1 head and Pump M1 speed obtained by solving the
decoupled and coupled pump scheduling problems of the Three-node network under
the (a) first and (b) second hydraulic scenarios.

With that in mind, first we apply the decoupled and coupled problem on the
Three-node network under different scenario. The first scenario imposes no
restrictions on the times when Tank TK1 is in an on-demand or off-demand
state. Tank TK1 has a minimum of 904 ft, maximum of 924 ft, and an initial
head level of 912 ft. Fig. 9a shows the optimal pump speed and Tank TK1 head
levels obtained by solving the decoupled and coupled problems. for this
scenario, the coupled problem is formulated to achieve higher and wider
controllability at the beginning of the simulation period by constraining it
with higher energy. This is done while filling the tank so that when it is in
the on-demand state; it can supply the system with water that has sufficient
chlorine levels. Filling the tank for this specific window is achieved by
initially operating the pump at a higher speed. However, this increase in the
pump speed is balanced later in the day when the tank supplies the network and
the pump speed is lower compared to the decoupled problem. In this scenario,
the water quality controller’s performance is effectively improved by
injecting chlorine at the start of the simulation period with fully
controllable system. This approach reduces the reliance on the water volume
stored in TK1, where chlorine concentration tends to decay over time. In
addition, when utilizing the coupled optimization approach for pump
scheduling, the total resulting cost is 8.9% higher compared to the decoupled
approach. In terms of computational efficiency, the solver completes the
decoupled problem in 3.8 seconds, while the coupled problem requires 4.7
seconds to reach a solution. The second hydraulic scenario that we adopt has
an initial tank head level at 908 ft and safe water level of 909.5 ft (Fig.
9b). The goal of maintaining a safe water level in TK1 leads to it being in
the filling state at the beginning of the simulation period. Additionally, the
coupled problem aims to higher controllability energy and coverage levels
resulting in higher velocity in Pipe P1. In addition, the flow directions for
both problems differ for the remainder of the simulation period, and the total
operational cost of pump operation for both problems is comparable. In fact,
the coupled problem’s total cost is lower by 2% in comparison to the decoupled
problem. Moreover, under zero initial chlorine concentrations at Tank TK1 and
Pipe P1, the WQ controller succeeds to achieve the set point concentration
faster by 20 minutes.

(a)  (b)

Figure 10: Optimal Tank 11 head and Pump speed obtained by solving the coupled
pump scheduling problem of the Net1 network under the (a) first and (b) second
hydraulic scenarios.

Second, we validate the performance of the coupled problem on Net1 network. It
is worth noting that the pre-computation of the target controllability Gramian
for this network requires approximately 2 seconds. Consequently, the runtime
listed below exclusively pertains to solving the optimization problem
separated for a 24-hour simulation period. As for the decoupled problem, the
time required to solve it for this network using Gurobi is 8.7 sec and 161.1
sec using BMIBNB. On the other hand, the computational time for the coupled
problem–where both solvers are used alternately during the hydraulic time-
steps based on the constraints developed for each time-step–amounts to 145.2
seconds. For the first scenario, the initial Tank 11 head is 920 ft, while the
minimum bound on the head is 910 ft. In this scenario, Junctions 2, 4, and 9
have water demands with different patterns and bases. Results from this
simulation scenario are illustrated in Fig. 10a. For the first 6 hrs of the
simulation period, Tank 11 is being filled and accordingly is included in the
target nodes while solving the rank-informed control problem. For the next
simulation window till the 20th hour, Tank 11 is on demand allowing the pump
to work on lower speeds. During this window, the coupled problem is formulated
as an energy-driven problem with the network’s junctions as the target nodes.
However, the water head in Tank 11 reaches the minimum head level which
requires activating a constraint in the optimization problem to recover the
water level in the tank for the next simulation period. During the last 4 hrs
of the simulation, the energy level required at start for some directions
resulted in an infeasible solutions and accordingly, the problem is solved
decoupled.

Another scenario for Net1, Tank 11 has initial head of 915 ft and minimum
safety head of 908 ft. In addition, the stored water volume in the tank has
sufficient chlorine levels of 1.5 mg/L to serve parts of the network, that is,
it is excluded from the target nodes. Also, the demands’ bases at the
junctions are higher than the first scenario. As shown in Fig. 10b, the pump
speed is lower for the first half of the simulation period to fulfill the
demands at the start of the network while Tank 11 is serving the rest. For the
second half, the tank reaches the minimum safety head and keeps alternating
between filling and emptying states. The coupled problem for this scenario if
formulated to achieve higher controllability energy as the tank helps with
serving the network with water that has sufficient chlorine levels and
accordingly higher network chlorine coverage, specifically at the start of the
simulation. The integration of the energy metrics in the hydraulic
optimization problem for this scenario results in balancing the injections of
the booster stations located at Junction 1 and Junction 6, in comparison to
scenarios where the booster station at the start of the network is overworked
and may run to an over capacity state.

Last, we implement our proposed quality-aware pump control framework on the
Richmond skeleton network which has a more branched layout. This network
comprises multiple elevated tanks at varying elevations, each serving distinct
areas of the network with differing and varying water demands. For this
network, we run multiple hydraulic scenario with varying target nodes
throughout the simulation window. The selection of target nodes is based on
the main pipeline leading to the tank during filling operations and the
inclusion of dead-end points covered by booster stations. Essentially, the
coupled problem for this network is designed to achieve extensive
controllability coverage in specific network segments while maintaining a high
controllability energy level within areas with low expected disinfectant
residuals. This case study serves as a validation of the proposed framework’s
performance on a complex, branched real network and provides insights into the
implementation of WQ control based on the obtained pump schedule. The average
computational time to solve the coupled optimization problem for these
scenarios is 432 seconds. Furthermore, across various scenarios involving
distinct hydraulic settings and constraints, the WQ controller demonstrates
improved and more efficient performance. This conclusion is drawn from the
following observations made during the aforementioned scenarios: (i) chlorine
is injected by booster stations near tanks that are being filled to reach
desired levels that are set to 2 mg/L; (ii) chlorine injections are balanced
among booster stations to distribute the workload. For instance, booster
station at J9 is not overworked for the whole network and alternates based on
tanks’ on- and off-demand schedules and J15 supports the network while
ensuring adequate residuals at J20; and (iii) reduction in chlorine injections
varies between 3-11% when compared to the decoupled framework for scenarios
with no firm restrictions on flow directions. However, for scenarios with such
constraints, the increase in injections does not exceed 9%, yet, dead-ends
maintain sufficient chlorine residuals.

##  VI Summary, Limitations, and Future Work

This paper provides a comprehensive analysis of how the accountability of WQ
controllability impacts the solution of the hydraulic settings optimization
problem. This accountability aims to enhance the performance of the WQ control
and regulation algorithms which is directly affected by the system hydraulics.
To accomplish this, we incorporate the goal of improving various
controllability metrics into the operational hydraulics optimization problem,
aiming to achieve a targeted level of water quality controllability throughout
the system. The performance of this approach is evaluated on three networks: a
Three-node network, Net1, and the Richmond skeleton network. These networks
vary in terms of scale and configurations. Additionally, different initial
hydraulic and quality dynamics are examined for each network.

The results demonstrate the efficacy of the proposed approach in enhancing WQ
controllability metrics, leading to a more efficient controller performance in
achieving desired chlorine levels across the network. However, the enhancement
of this performance is significantly influenced by variations in consumer
demand patterns, leading to substantial shifts in system hydraulics in certain
scenarios. Moreover, the network’s configuration has a pivotal role in
determining the operational schedule and chlorine injections, with limitations
imposed in some cases on the feasibility of chlorine reaching uncontrollable
regions due to resultant flow directions. In conclusion, this approach is a
step-forward in advancing the performance of regulating the water distribution
networks dynamics and further improvement can be achieved by including these
aspects while designing the network’s layout and the functional constraints,
which can be posed as work extensions in our group’s future studies. Other
future work encompasses the formulation of joint real-time WDNs control where
hydraulics are updated with feedback that implies the foreseen effect on both;
future hydraulic and WQ settings and dynamics not an optimization problem that
is solved every hydraulic time-step with no feedback from the future of the
impact of the control action.

It is important to acknowledge certain limitations in our study, particularly
regarding the assumption of pre-allocated booster station locations throughout
the network, despite their direct influence on water quality controllability.
In our study, we set a specific good configuration that ensures a satisfactory
level of controllability across the entire network under different scenarios
of hydraulic settings. Another limitation of our study is the reliance on
simplified models (e.g., approximate pump power, pump affinity rules, and pipe
head loss dynamics piece-wise linearization) and assumptions, which may not
capture the full complexity and variability of real-world hydraulic dynamics.

##  Appendix A Water Quality Controllability Dependency on System
Hydraulics—An Example

In this appendix, we demonstrate the dependence of water quality
controllability on system hydraulics using a straightforward example of the
Three-node network. In our example, flow directions are assumed to be as
illustrated in Fig. 3. Tank TK1 is considered filling and off-demand and a
booster station is located at Junction J1 dosing chlorine into the system at a
rate of $\displaystyle q_{1}^{\mathrm{B}}(t)$. The WQ system matrices of Eq.
(12) are as expressed in (27) for this scenario. Note that, to calculate the
concentrations at Junction J1 at time-step $\displaystyle t+\Delta
t_{\mathrm{WQ}}$ following Eq. (7), the flow rates need to be at the same
time-step. Nevertheless, the water quality time-step operates at the
seconds/minutes scale, whereas the hydraulic time-step is on an hourly scale.
That is, within the same hydraulic time-step, $\displaystyle q(t+\Delta
t_{\mathrm{WQ}})=q(t)$ for all links.

$\displaystyle\displaystyle\bm{A}_{\mathrm{WQ}}(t)$
$\displaystyle\displaystyle={\small\begin{tabular}[]{c
c}$\displaystyle\begin{bmatrix}1&0&0&0&0&0&0\\\
0&0&0&{\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}{a}_{{\mathrm{J}}}}&0&0&0\\\
0&0&{\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}{a}_{{\mathrm{TK}}}}&0&0&0&{\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}\underline{a}_{{\mathrm{TK}}}}\\\
1&0&0&0&0&0&0\\\
0&{\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}\underline{a}_{{\mathrm{P}}}}&0&0&{\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}{a}_{{\mathrm{P}}}}&0&0\\\
0&0&0&0&{\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}\underline{a}_{{\mathrm{P}}}}&{\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}{a}_{{\mathrm{P}}}}&0\\\
0&0&0&0&0&{\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}\underline{a}_{{\mathrm{P}}}}&{\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}{a}_{{\mathrm{P}}}}\end{bmatrix}$&$\displaystyle\begin{matrix}\begin{array}[]{c}${\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>${\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>${\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>${\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>\\\
${\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>\mathbb{}\end{array}\end{matrix},$\end{tabular}}$ (27b)
$\displaystyle\displaystyle\bm{B}_{\mathrm{WQ}}(t)$
$\displaystyle\displaystyle={\small\begin{tabular}[]{c
c}$\displaystyle\begin{bmatrix}0&0\\\
{\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}{a}^{\mathrm{B}}_{{\mathrm{J}}}}&0\\\
0&0\\\ 0&0\\\ 0&0\\\ 0&0\\\
0&0\end{bmatrix}$&$\displaystyle\begin{matrix}\begin{array}[]{c}${\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>${\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>${\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>${\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>\\\
${\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>\mathbb{}\end{array}\end{matrix},$\end{tabular}}$ (27d)

where

$\displaystyle\displaystyle{\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}{a}_{{\mathrm{J}}}}=\frac{{\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}q_{1}^{\mathrm{M}}(t)}}{{\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}q^{\mathrm{D}_{\mathrm{J}}}_{1}(t)}+{\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}q^{\mathrm{P}}_{1}(t)}},$
$\displaystyle\displaystyle{\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}{a}^{\mathrm{B}}_{{\mathrm{J}}}}=\frac{{\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}q_{1}^{\mathrm{B}}(t)}}{{\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}q^{\mathrm{D}_{\mathrm{J}}}_{1}(t)}+{\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}q^{\mathrm{P}}_{1}(t)}},$
$\displaystyle\displaystyle{\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}{a}_{{\mathrm{TK}}}}=\frac{(1-k_{b}){\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}V_{1}^{\mathrm{TK}}(t)}}{{\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}V_{1}^{\mathrm{TK}}(t+\Delta
t_{\mathrm{WQ}})}},$
$\displaystyle\displaystyle{\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}\underline{a}_{{\mathrm{TK}}}}=\frac{{\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}q^{\mathrm{P}}_{1}(t)\Delta
t_{\mathrm{WQ}}}}{{\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}V_{1}^{\mathrm{TK}}(t+\Delta
t_{\mathrm{WQ}})}},$
$\displaystyle\displaystyle{\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}{a}_{{\mathrm{P}}}}=1-{\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}\lambda(t)}-k_{b}\Delta
t_{\mathrm{WQ}},$
$\displaystyle\displaystyle{\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}\underline{a}_{{\mathrm{P}}}}={\color[rgb]{.5,0,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,0,.5}\lambda(t)}.$

The formulation of the water quality system matrices clearly demonstrates a
strong dependency on the system hydraulics, which extends to the water quality
controllability matrix and Gramian as well. However, when these matrices are
multiplied to calculate the water quality controllability matrix and Gramian
in (14) and (15), the nonlinearity order increases exponentially.

## References

  * [1] V. Puleo, M. Morley, G. Freni, and D. Savić, “Multi-stage Linear Programming Optimization for Pump Scheduling,” _Procedia Engineering_ , vol. 70, pp. 1378–1385, 2014.
  * [2] I. Fisher, G. Kastl, F. Shang, and A. Sathasivan, “Framework for Optimizing Chlorine and Byproduct Concentrations in Drinking Water Distribution Systems.” _Journal: American Water Works Association_ , vol. 110, no. 11, 2018.
  * [3] G. Bonvin and S. Demassey, “Extended linear formulation of the pump scheduling problem in water distribution networks,” in _International Network Optimization Conference, June 12-14, 2019, Avignon, France._ , 2019, pp. 13–18.
  * [4] G. Bonvin, S. Demassey, and A. Lodi, “Pump scheduling in drinking water distribution networks with an LP/NLP-based branch and bound,” _Optimization and Engineering_ , vol. 22, no. 3, pp. 1275–1313, Sep. 2021\.
  * [5] M. K. Singh and V. Kekatos, “Optimal Scheduling of Water Distribution Systems,” Aug. 2019.
  * [6] D. Kang, “Real-time Optimal Control of Water Distribution Systems,” _Procedia Engineering_ , vol. 70, pp. 917–923, 2014.
  * [7] Y. Wang, V. Puig, and G. Cembrano, “Non-linear economic model predictive control of water distribution networks,” _Journal of Process Control_ , vol. 56, pp. 23–34, 2017.
  * [8] Y. Wang, K. Too Yok, W. Wu, A. R. Simpson, E. Weyer, and C. Manzie, “Minimizing Pumping Energy Cost in Real-Time Operations of Water Distribution Systems Using Economic Model Predictive Control,” _Journal of Water Resources Planning and Management_ , vol. 147, no. 7, p. 04021042, Jul. 2021.
  * [9] S. Wang, A. F. Taha, N. Gatsis, and M. H. Giacomoni, “Receding Horizon Control for Drinking Water Networks: The Case for Geometric Programming,” _IEEE Transactions on Control of Network Systems_ , vol. 7, no. 3, pp. 1151–1163, Sep. 2020.
  * [10] A. Marchi, A. R. Simpson, and M. F. Lambert, “Pump Operation Optimization Using Rule-based Controls,” _Procedia Engineering_ , vol. 186, pp. 210–217, 2017.
  * [11] A. L. Reis, M. A. Lopes, A. Andrade-Campos, and C. Henggeler Antunes, “A review of operational control strategies in water supply systems for energy and cost efficiency,” _Renewable and Sustainable Energy Reviews_ , vol. 175, p. 113140, Apr. 2023.
  * [12] D. L. Boccelli, M. E. Tryby, J. G. Uber, L. A. Rossman, M. L. Zierolf, and M. M. Polycarpou, “Optimal scheduling of booster disinfection in water distribution systems,” _Journal of water resources planning and management_ , vol. 124, no. 2, pp. 99–111, 1998.
  * [13] M. E. Tryby, D. L. Boccelli, J. G. Uber, and L. A. Rossman, “Facility location model for booster disinfection of water supply networks,” _Journal of Water Resources Planning and Management_ , vol. 128, no. 5, pp. 322–333, 2002\.
  * [14] Z. Ohar and A. Ostfeld, “Optimal design and operation of booster chlorination stations layout in water distribution systems,” _Water research_ , vol. 58, pp. 209–220, 2014.
  * [15] S. Wang, A. F. Taha, and A. A. Abokifa, “How Effective is Model Predictive Control in Real-Time Water Quality Regulation? State-Space Modeling and Scalable Control,” _Water Resources Research_ , vol. 57, no. 5, May 2021.
  * [16] S. M. Elsherif, A. F. Taha, A. A. Abokifa, and L. Sela, “Comprehensive Framework for Controlling Nonlinear Multi-Species Water Quality Dynamics,” Feb. 2023.
  * [17] A. B. A. Sakarya and L. W. Mays, “Optimal Operation of Water Distribution Pumps Considering Water Quality,” _Journal of Water Resources Planning and Management_ , vol. 126, no. 4, pp. 210–220, Jul. 2000.
  * [18] A. Ostfeld and E. Salomons, “Conjunctive optimal scheduling of pumping and booster chlorine injections in water distribution systems,” _Engineering Optimization_ , vol. 38, no. 3, pp. 337–352, Apr. 2006.
  * [19] M. Drewa, M. Brdys, and A. Cimiński, “Model Predictive Control of Integrated Quantity and Quality in Drinking Water Distribution Systems,” _IFAC Proceedings Volumes_ , vol. 40, no. 5, pp. 95–100, 2007\.
  * [20] M. Xie and M. Brdys, “Nonlinear Model Predictive Control of Water Quality in Drinking Water Distribution Systems with DBPs Objectives,” vol. 9, no. 8, 2015.
  * [21] M. Abdallah and Z. Kapelan, “Fast Pump Scheduling Method for Optimum Energy Cost and Water Quality in Water Distribution Networks with Fixed and Variable Speed Pumps,” _Journal of Water Resources Planning and Management_ , vol. 145, no. 12, p. 04019055, Dec. 2019.
  * [22] A. G. Seyoum and T. T. Tanyimboh, “Integration of hydraulic and water quality modelling in distribution networks: EPANET-PMX,” _Water Resources Management_ , vol. 31, pp. 4485–4503, 2017.
  * [23] A. S. Zamzam, E. Dall’Anese, C. Zhao, J. A. Taylor, and N. D. Sidiropoulos, “Optimal Water–Power Flow-Problem: Formulation and Distributed Optimal Solution,” _IEEE Transactions on Control of Network Systems_ , vol. 6, no. 1, pp. 37–47, Mar. 2019.
  * [24] F. M. White, _Fluid Mechanics_ , 1966.
  * [25] R. K. Linsley, J. B. Franzini, and G. Tchobanoglous, “Water resources engineering,” _McGrow-Hill Book Co_ , 1997.
  * [26] D. C. Rennels, _Pipe Flow: A Practical and Comprehensive Guide_. John Wiley & Sons, 2022.
  * [27] P. F. Boulos, K. E. Lansey, and B. W. Karney, _Comprehensive Water Distribution Systems Analysis Handbook for Engineers and Planners_. American Water Works Association, 2006.
  * [28] J. B. Burkhardt, J. Szabo, S. Klosterman, J. Hall, and R. Murray, “Modeling fate and transport of arsenic in a chlorinated distribution system,” _Environmental Modelling & Software_, vol. 93, pp. 322–331, 2017.
  * [29] B. Koren, _A Robust Upwind Discretization Method for Advection, Diffusion and Source Terms_. Centrum voor Wiskunde en Informatica Amsterdam, 1993, vol. 45.
  * [30] S. M. Elsherif, S. Wang, A. F. Taha, L. Sela, M. H. Giacomoni, and A. A. Abokifa, “Control-theoretic modeling of multi-species water quality dynamics in drinking water networks: Survey, methods, and test cases,” _Annual Reviews in Control_ , 2022.
  * [31] R. E. Kalman, “Mathematical description of linear dynamical systems,” _Journal of the Society for Industrial and Applied Mathematics, Series A: Control_ , vol. 1, no. 2, pp. 152–192, 1963.
  * [32] J. P. Hespanha, _Linear Systems Theory_ , 2nd ed. Princeton: Princeton University Press, 2018.
  * [33] A. Olshevsky, “On (Non)Supermodularity of Average Control Energy,” _IEEE Transactions on Control of Network Systems_ , vol. 5, no. 3, pp. 1177–1181, Sep. 2018.
  * [34] F. Pasqualetti, S. Zampieri, and F. Bullo, “Controllability Metrics, Limitations and Algorithms for Complex Networks,” _IEEE Transactions on Control of Network Systems_ , vol. 1, no. 1, pp. 40–52, Mar. 2014\.
  * [35] T. H. Summers, F. L. Cortesi, and J. Lygeros, “On Submodularity and Controllability in Complex Dynamical Networks,” _IEEE Transactions on Control of Network Systems_ , vol. 3, no. 1, pp. 91–101, Mar. 2016\.
  * [36] P. A. Fuhrmann, “On weak and strong reachability and controllability of infinite-dimensional linear systems,” _Journal of Optimization Theory and Applications_ , vol. 9, no. 2, pp. 77–89, Feb. 1972.
  * [37] B. Datta, _Numerical Methods for Linear Control Systems_. Academic Press, 2004.
  * [38] R. Menke, E. Abraham, P. Parpas, and I. Stoianov, “Approximation of System Components for Pump Scheduling Optimisation,” _Procedia Engineering_ , vol. 119, pp. 1059–1068, 2015.
  * [39] R. Menke, E. Abraham, and I. Stoianov, “Modeling Variable Speed Pumps for Optimal Pump Scheduling,” in _World Environmental and Water Resources Congress 2016_. West Palm Beach, Florida: American Society of Civil Engineers, May 2016, pp. 199–209.
  * [40] K. Horváth, B. Esch, D. Vreeken, I. Pothof, and J. Baayen, “Convex Modeling of Pumps in Order to Optimize Their Energy Use,” _Water Resources Research_ , vol. 55, no. 3, pp. 2432–2445, Mar. 2019.
  * [41] A. Gokhale, S. M. Valli, R. Kalaimani, and R. Pasumarthy, “Optimizing controllability metrics for target controllability,” in _2021 Seventh Indian Control Conference (ICC)_. Mumbai, India: IEEE, Dec. 2021, pp. 141–146.
  * [42] V. S. Mai, D. Maity, B. Ramasubramanian, and M. C. Rotkowitz, “Convex Methods for Rank-Constrained Optimization Problems,” in _2015 Proceedings of the Conference on Control and Its Applications_ , C. Bonnet, B. Pasik-Duncan, H. Ozbay, and Q. Zhang, Eds. Philadelphia, PA: Society for Industrial and Applied Mathematics, Jul. 2015, pp. 123–130.
  * [43] S. Wang, A. F. Taha, and A. A. Abokifa, “How Effective is Model Predictive Control in Real-Time Water Quality Regulation? State-Space Modeling and Scalable Control,” _Water Resources Research_ , vol. 57, no. 5, p. e2020WR027771, 2021.
  * [44] O. C. Acrylamide, “National Primary Drinking Water Regulations,” _Kidney_ , vol. 2, no. 4-D, p. 0.07, 2009.
  * [45] L. A. Rossman, H. Woo, M. Tryby, F. Shang, R. Janke, and T. Haxton, “EPANET 2.2 User Manual Water Infrastructure Division, Center for Environmental Solutions and Emergency Response,” 2020.