# Economic Predictive Control with Periodic Horizon for Water Distribution
Networks

Mirhan Ürkmez  Carsten Kallesøe  Jan Dimon Bendtsen  John Leth  Aalborg
University, Fredrik Bajers Vej 7c, DK-9220 Aalborg, Denmark
(e-mail<EMAIL_ADDRESS>
###### Abstract

This paper deals with the control of pumps in large-scale water distribution
networks with the aim of minimizing economic costs while satisfying
operational constraints. Finding a control algorithm in combination with a
model that can be applied in real-time is a challenging problem due to the
nonlinearities presented by the pipes and the network sizes. We propose a
predictive control algorithm with a periodic horizon. The method provides a
way for the economic operation of large water networks with a small linear
model. Economic Predictive control with a periodic horizon and a terminal
state constraint is constructed to keep the state trajectories close to an
optimal periodic trajectory. Barrier terms are also included in the cost
function to prevent constraint violations. The proposed method is tested on
the EPANET implementation of the water network of a medium size Danish town
(Randers) and shown to perform as intended under varying conditions.

###### keywords:

Water distribution networks, Pump Scheduling, Predictive control, Periodic
horizon, Economic model predictive control

††thanks: This work is funded by Independent Research Fund Denmark (DFF). We
acknowledge Verdo company, Peter Nordahn, and Steffen Schmidt for providing us
with the EPANET model and the network information.

##  1 Introduction

Water distribution networks (WDNs) deliver drinkable water from water sources
to consumers using elements such as pumps, pipes, tanks etc. About $7\%-8\%$
of the world’s energy is used for water production and distribution (Sharif et
al., 2019). Water pumps account for a significant part of the energy required
for water distribution with their percentage ranging from $90\%$ to $95\%$ of
the total (Abdelsalam and Gabbar, 2021). There have been many works to
schedule the operation of the pumps in WDNs with proper methods so as to
reduce energy costs. However, pump scheduling is not an easy task because of
the nonlinearities governing the network elements. The problem gets
complicated with increasing network size. Also, there are constraints to be
satisfied such as limits on tank levels.

In the literature, WDNs with both constant and variable speed pumps are
studied extensively. The control input is turning on and off the pump for
constant speed pumps. In Lindell Ormsbee (2009), a constant-speed pump
scheduling problem is posed as an optimization problem in which the decision
variables are the operation times of the pumps and the objective is the energy
cost. After observing that the optimal solution would be not running the pumps
at all without the constraints, the authors try to find the solution closest
to the origin that also complies with the constraints. The proposed way to
find such a solution is using a Genetic Algorithm (GA). In Bagirov et al.
(2013), the Hooke-Jeeves method is used for finding optimal pump operating
times for a similar problem. Then, network simulation algorithms are used to
check if the constraints are satisfied. In Castro-Gama et al. (2017), binary
decision variables are used to represent the opening and the closing of each
pump. The feasibility of the solution found with GA is checked with EPANET, a
WDN modeling software, and a high cost is assigned to the infeasible
solutions. The number of open pumps is also taken as the input to the system
in some works, e.g., Wang et al. (2021). The problem is then solved using
mixed-integer nonlinear programming. In Berkel et al. (2018), a network in
which pressure zones are connected via constant speed pumps is considered.
Each pressure zone is treated as a subsystem and distributed model predictive
control (DMPC) is applied.

The flow rate of the pumps should be determined for networks with variable
speed pumps. In Pour et al. (2019), Linear Parameter Varying (LPV) system
modeling is used to replace the nonlinear part of the network, and an Economic
Model Predictive Control (EMPC) is applied on top of the LPV system to find
the optimal flow rates. In Kallesøe et al. (2017), a network structure with an
elevated reservoir is considered. Available data is used for the
identification of a reduced system model. Then, EMPC is applied to the model.
In the EMPC formulation, node pressures are not constrained. It is assumed
that the pressures would be in the accepted range because there is an elevated
reservoir. Relaxation of the original problem into a simpler one is commonly
used because of the large network sizes. The relaxation is generally achieved
by approximating the nonlinear pipe equations with some sort of linear
equations or inequalities. In Baunsgaard et al. (2016), pipe equations are
linearized around an operating point, and model predictive control (MPC) is
applied. In Wang et al. (2018), an EMPC is applied to a network where the
nonlinear pipe equations are relaxed into a set of linear inequalities. Before
simplifying the system model, the network structure is also simplified in
Fiedler et al. (2020). A hierarchical clustering method is used to represent
the original network with a smaller one which originally had 378 junctions. A
system model is derived from the simplified structure using a Deep Neural
Network (DNN) structure. Lagrangian relaxation is used to approximate the
original problem in Ghaddar et al. (2015).

In this paper, a way for optimal pump scheduling of large-scale WDNs is
presented. To control the pumps, a linear model of the system is derived.
Then, a predictive control method with a periodic horizon is constructed.
Barrier functions are utilized to prevent constraint violation due to the
model-plant mismatch. With the introduction of the periodic horizon and the
terminal state constraint, the chance of finding a feasible solution is
increased by keeping trajectories close to an optimal periodic trajectory. The
method is applied to a medium-sized Danish town’s network (Randers).

The outline of the rest of the paper is as follows. The network model is given
in Section 2. The proposed control method is explained in Section 3. The
experimental results are presented in Section 4. The paper is concluded with
final remarks in Section 5.

##  2 Network Model

A typical water distribution network consists of pipes, pumps, tanks, junction
nodes and reservoirs. Water in the network flows from high hydraulic head to
low head. Hydraulic head is a measure of the fluid pressure and is equal to
the height of a fluid in a static column at a point.

Hydraulic head loss occurring in a pipe can be approximated by the Hazen-
Williams Equation as

$\Delta h=h_{1}-h_{2}=Kq\lvert q\rvert^{0.852}$ (1)

where $K$ is the pipe resistance that depends on the physical features of a
pipe such as diameter and length, $q$ is the flow rate, and $h_{1}$ and
$h_{2}$ are the heads at the two ends of the pipe.

At each node $j$, the mass conservation law is satisfied. It can be expressed
as

$\sum_{i\in\mathpzc{N}_{j}}q_{ij}=d_{j}$ (2)

where $q_{ij}$ is the flow entering the node $j$ from node $i$ and $d_{j}$ is
the demand at node $j$, which is the water requested by the user at node $j$.
The symbol $\mathpzc{N}_{j}$ denotes the set of neighbor nodes of node $j$.
Note that $q_{ij}$ is positive if the flow is from node $i$ to the neighbor
node $j$ and negative vice versa.

Tanks are storage elements that provide water to the users. In the network,
tanks are elevated so that water can be pressurized enough to be delivered to
the consumers. The change in the water level of a tank is dependent on the
flow coming from neighbor nodes and can be written for the tank $j$ as

$A_{j}\dot{{h}}_{j}=\sum_{i\in\mathpzc{N}_{j}}q_{ij}$ (3)

where $A_{j}$ is the cross-sectional area, $h_{j}$ is the level of the tank.
Tank levels change according to the flow passing through the pipes connected
to the tanks. Those flows are determined by a set of pipe head loss equations
(1), and mass balance equations (2) throughout the whole network. As Equation
(1) is nonlinear, flow through pipes connected to the tanks are nonlinear
functions $f_{i}$ of the demand at each node, tank levels, and the amount of
water coming from the pumps. Explicit forms of those nonlinear functions could
be derived if the vector $d=[d_{1},d_{2}...]^{T}$ containing the demands of
all the nodes is known, which is not possible unless demand data for all nodes
are available. In our work, we assume that the total demand of the zones that
are supplied by the pumps can be estimated through available data with time
series analysis methods, but not require $d$ vector to be known. Since $f_{i}$
functions can not be found without $d$ vector, we approximate them using
linear models and write tank level change equations as

$\dot{{h}}(t)=Ah(t)+B_{1}u(t)+B_{2}d_{a}(t)$ (4)

where $h(t)\in\mathbb{R}^{n}$ includes tank levels, $A\in\mathbb{R}^{n\times
n}$, $B_{1}\in\mathbb{R}^{n\times m}$, $B_{2}\in\mathbb{R}^{n\times 1}$ are
constant system matrices and $d_{a}(t)$ is the aggregated demand of controlled
zone at time $t$, $u(t)\in\mathbb{R}^{m}$ is the input containing pump flows.
The reason we chose a linear model is to increase the chance of finding a
feasible solution for the controller which is posed as an optimization problem
in the next section. Although capturing the full dynamics of a large-scale
network is not possible with a linear model, the proposed control method is
designed to compensate for model inaccuracies and we have observed that it was
enough to control the system while satisfying the constraints.

##  3 Periodic Horizon Control

In this section, a predictive control algorithm for pump scheduling is
presented to minimize the economical costs. The aim is to run the pumps when
the electricity price is low and let tanks deliver water when the price is
high while also satisfying input and output constraints. The problem at time
$t$ is posed as

$\displaystyle\min_{u_{0}^{t},u_{1}^{t}\cdots
u_{N(t)-1}^{t}}\sum_{j=0}^{N(t)-1}J(h_{j}^{t},u_{j}^{t})$ (5a) $\displaystyle
h_{j}^{t}=A_{d}h_{j-1}^{t}+B_{d1}u_{j-1}^{t}+B_{d2}d_{a}(j-1)$ (5b)
$\displaystyle h_{0}^{t}=h(t)$ (5c) $\displaystyle
u_{j}^{t}\in\mathcal{U}\subseteq\mathbb{R}^{m}$ (5d) $\displaystyle
h_{j}^{t}\in\mathcal{H}\subseteq\mathbb{R}^{n}$ (5e) $\displaystyle
h_{N(t)}^{t}\in\mathcal{H}_{tf}\subseteq\mathbb{R}^{n}$ (5f)

where $J(h_{j}^{t},u_{j}^{t})$ is the economic cost function,
$h^{t}=[h_{1}^{t}\cdots h_{N(t)}^{t}]\in\mathbb{R}^{n\times N(t)}$ is the
predicted future states, $u_{j}^{t}$ is the input vector, $N(t)$ is the
prediction horizon, $\mathcal{U}\subseteq\mathbb{R}^{m}$ and
$\mathcal{H}\subseteq\mathbb{R}^{n}$ denotes the input and state constraints
respectively and $\mathcal{H}_{tf}\subseteq\mathbb{R}^{n}$ is the terminal
state set. The continuous system (4) is discretized and (5b) is obtained. The
optimization problem (5) is solved at every time step separated by
$\Delta_{t}$ and the first term $u_{0}^{t}$ of the optimal input sequence
$\mathbf{u^{t}}=[u_{0}^{t}\cdots u_{N(t)-1}^{t}]\in\mathbb{R}^{m\times N(t)}$
is applied to the system.

Input constraints come from the physical limitations and working principles of
the pumps. A pump can not provide water in the opposite direction and it can
deliver a maximum amount of water per unit of time. These conditions are
expressed as

$\mathcal{U}=\\{[u_{1},\cdots u_{m}]\in\mathbb{R}^{m}\mid 0\leq
u_{1}\leq\overline{u}_{1},\cdots 0\leq u_{m}\leq\overline{u}_{m}\\}$ (6)

where $\overline{u}_{1}\cdots\overline{u}_{m}$ are upper flow limits. Tank
levels are also constrained so that there is always enough water in the tanks
in case of an emergency and there is no overflow of water. The set
$\mathcal{H}$ can be defined as

$\mathcal{H}=\\{[h_{1},\cdots h_{2}]\in\mathbb{R}^{n}\mid\tilde{h}_{1}\leq
h_{1}\leq\overline{h}_{1},\cdots\tilde{h}_{n}\leq
h_{n}\leq\overline{h}_{n}\\}$ (7)

The cost function includes the electricity costs of the pumps. The power
provided to the network by the pump $i$ is equal to
${q}_{pi}(p^{out}_{i}-p^{in}_{i})$, where ${q}_{pi}$ is the pump flow,
$p_{out}^{i}$ and $p_{in}^{i}$ are the outlet and inlet pressures of the pump
$i$. The inlet pressures $p^{in}=[p^{in}_{1},p^{in}_{2}]$ are the pressures of
the related reservoirs and are assumed to be constant. The outlet pressures
$p^{out}=[p^{out}_{1},p^{out}_{2}]$ are given as the output of the linear
model

$p^{out}(t)=A_{p}h(t)+B_{p}u(t)$ (8)

where $A_{p}$ and $B_{p}$ are found using system identification on data
generated by the EPANET model. Electricity cost at time $t$ is then found by
multiplying total power consumption $u(t)^{T}(p^{out}(t)-p^{in}(t))$ with the
electricity price $c(t)$.

We acknowledge a certain degree of model-plant mismatch by using a linear
model (4) to represent the whole network. This causes actual states $h(t)$ to
be different than the predicted states $h^{t}$. We know that the predicted
states satisfy the state constraints (7) since they are the solution to the
optimization problem 5, but the actual states might violate them. To ensure
the satisfaction of the state constraints with the model-plant mismatch, we
introduce new terms to the cost function. First, we rewrite state constraints
(7) as

$C_{i}(h)\leq 0,\quad i=0,1,\cdots 2\times n-1$ (9)

where $C_{0}(h)=\tilde{h}_{1}-h_{1}$ and the rest of the $C_{i}$ functions are
chosen in a similar manner. The cost function terms are then defined as

$J_{hi}(h)=e^{a_{i}(C_{i}(h)+b_{i})}\quad i=0,1,\cdots 2\times n-1$ (10)

where $a_{i},b_{i}\in\mathbb{R}_{>0}$. This can be seen as an exponential
barrier function. The parameters $a_{i},b_{i}$ determine a dangerous region
close to the boundaries of the state constraints where cost function $J_{hi}$
attains high values. The predicted optimal state trajectories $h^{t}$ do not
enter the dangerous region if possible because of the high cost values in the
dangerous region. Then, the actual states $h(t)$ do not violate the state
constraints (7) assuming the difference between the predicted state and the
actual state is small. If the state trajectory enters one of the dangerous
regions at any step due to the model-plant mismatch, then the cost function
will try to drive the trajectory out of the region.

The overall cost function includes both the electricity expense term and the
constraint barrier functions and it can be expressed as

$J(h(t),u(t))=c(t)u(t)^{T}(p^{out}(t)-p^{in}(t))+\sum_{i=0}^{2\times
n-1}J_{hi}(h(t))$ (11)

Both electricity price $c(t)$ and total water demand $d_{a}(t)$ signals can be
viewed as consisting of a periodic signal with a period of 1 day and a
relatively small deviation signal. This can be leveraged to find a feasible
controller. Suppose a sequence of inputs can be found for some initial tank
levels such that levels after 1 day are equal to initial levels. In that case,
the problem after 1 day is the same as in the beginning assuming deviation
signals of the electricity price and the demand are zero, hence they are
periodic. Then, the input sequence from the previous day could be applied and
produce the same path for tank levels. Taking into account the deviation
signals and supposing that a solution exists such that levels after 1 day are
close to initial levels, the input sequence from the previous day could be a
good point of start to search for a feasible solution if the map from the
initial conditions and the demand profile to the optimal input sequences is
continuous. Therefore, we choose a terminal state constraint for the end of
each day to increase the chance of finding a feasible solution.

$\Delta_{t}$$N(t)\Delta_{t}$$N(t+\Delta_{t})\Delta_{t}$$h(t)$
$h(t+\Delta_{t})$ $h_{1}^{t}$
$u_{0}^{t}$$u_{0}^{t+\Delta_{t}}$$h^{t+\Delta_{t}}$ $h^{t}$
$\mathcal{B}_{r}(h_{T_{day}/\Delta_{t}}^{*})$

Figure 1: Predicted state trajectories $h^{t},h^{t+\Delta_{t}}$ at times
$t,t+\Delta_{t}$. Sampling time $\Delta_{t}$, prediction horizons
$N(t),N(t+\Delta_{t})$ and the applied inputs $u_{0}^{t},u_{0}^{t+\Delta_{t}}$
are shown. The true state $h(t+\Delta_{t})$ and the predicted state
$h_{1}^{t}$ are indicated to emphasize the deviation from the prediction. The
terminal set $\mathcal{B}_{r}(h_{T_{day}/\Delta_{t}}^{*})$ is also
illustrated.

Now, the remaining problem is to decide which tank levels should the
trajectories turn back to at the end of each day. We define the optimal
periodic trajectory of the system as the solution of

$\displaystyle(\mathbf{u^{*}},\mathbf{h^{*}})=\operatorname*{arg\,min}_{{u_{i}},{h_{i}}}\sum_{i=0}^{(T_{day}/\Delta_{t})-1}J({h_{i}},{u_{i}})$
(12a)
$\displaystyle{h_{i}}=A_{d}{h_{i-1}}+B_{d1}{u_{i-1}}+B_{d2}d_{a}^{*}(i-1)$
(12b) $\displaystyle{u_{i}}\in\mathcal{U}\subseteq\mathbb{R}^{m}$ (12c)
$\displaystyle{h_{i}}\in\mathcal{H}\subseteq\mathbb{R}^{n}$ (12d)
$\displaystyle{h_{0}}=h_{T_{day}/\Delta_{t}}$ (12e)

where $T_{day}$ is the duration of a whole day, $d_{a}^{*}$ is the average
daily demand profile obtained from the past measurements. The resulting state
trajectory $\mathbf{h^{*}}=[h^{*}_{0}\cdots
h^{*}_{T_{day}/\Delta_{t}}]\in\mathbb{R}^{n\times(T_{day}/\Delta_{t}+1)}$ is
the optimal periodic trajectory because of the constraint (12e). The terminal
set $\mathcal{H}_{tf}$ and the prediction horizon $N(t)$ is chosen to make
tank levels at the end of each day close to $h_{T_{day}/\Delta_{t}}^{*}$. At
any time $t$, $t+N(t)\Delta_{t}$ should be equal to the end of the day.
$\mathcal{H}_{tf}$ and $N(t)$ could be written as

$\displaystyle\mathcal{{H}}_{tf}=\mathcal{B}_{r}(h_{T_{day}/\Delta_{t}}^{*})$
(13a) $\displaystyle N(t)=(T_{day}-t\bmod T_{day})/\Delta_{t}$ (13b)

where $\mathcal{B}_{r}(h_{T_{day}/\Delta_{t}}^{*})$ is the open ball centered
at $h_{T_{day}/\Delta_{t}}^{*}$ with radius $r$. Note that $N(t)$ changes so
that $t+N(t)\Delta_{t}$ is the end of the day for all $t$. With these
definitions, the condition (13a) will translate to tank levels at the end of
the day being close to the final point in optimal periodic trajectory
$h_{T_{day}/\Delta_{t}}^{*}$ as shown in Figure 1. Therefore, not only chance
of finding a feasible solution is increased but also the solutions are kept
around the optimal periodic trajectory $\mathbf{h^{*}}$. If the problem (5)
becomes infeasible at any time step $t$, we apply the second term of the input
sequence from the previous step $u_{1}^{t-\Delta_{t}}$. The reason behind this
choice is as follows: If we apply the optimal control input
$u_{0}^{t-\Delta_{t}}$ to the network model (4) at time $t-\Delta_{t}$, then
the optimal sequence in the next time step will be
$\mathbf{u^{t}}=[u_{1}^{t-\Delta_{t}}\cdots
u_{N(t-\Delta_{t})-1}^{t-\Delta_{t}}]$ following Bellman’s principle of
optimality. Then, at time $t$, $u_{1}^{t-\Delta_{t}}$ will be applied to the
system as calculated at $t-\Delta_{t}$. Assuming the model-plant mismatch is
small enough, $u_{1}^{t-\Delta_{t}}$ is still a good input candidate if the
problem is infeasible at time $t$.

##  4 Application

High ZoneLow Zone

Figure 2: Water Distribution Network of Randers. The pumping stations to be
controlled are shown in red. Tanks are shown with a ’T’ shaped symbol in
yellow.

The presented method is applied to WDN of Randers, a Danish city, which is
shown in Figure 2. The network consists of 4549 nodes and 4905 links
connecting them. There are 8 pumping stations in the network, 6 of which are
shown in the figure whereas the other 2 are stationed where tanks are placed.
The goal is to derive the schedules for 2 of the pumping stations while other
pumps are already working according to some predetermined strategies. The
stations to be controlled are shown in red in the figure. Their task is to
deliver water mostly to the High Zone (HZ) and Low Zone (LZ). However,
connections exist between HZ-LZ and the rest of the city, so we can not think
of the system as composed of isolated networks entirely. There are also 3
tanks in the HZ. While 2 of them are directly connected via pipes, the third
one stands alone as shown in the figure.

Figure 3: Structure of the WDN.

The overall structure of the Randers WDN with tanks and pumps to be controlled
are given in Figure 3. There are 3 water tanks in the network, 2 of which have
been connected with a pipe directly. The tank level changes can be written by
applying the mass conservation law (3) to the tanks in Figure 3 as

$\displaystyle A_{1}\dot{{h}}_{1}$ $\displaystyle=q_{1down}+q_{1up}+q_{inter}$
(14a) $\displaystyle=f_{1}(h_{1},h_{2},h_{3},q_{p1},q_{p2},d),$ $\displaystyle
A_{2}\dot{{h}}_{2}$ $\displaystyle=q_{2down}+q_{2up}-q_{inter}$ (14b)
$\displaystyle=f_{2}(h_{1},h_{2},h_{3},q_{p1},q_{p2},d),$ $\displaystyle
A_{3}\dot{{h}}_{3}$
$\displaystyle=q_{3}=f_{3}(h_{1},h_{2},h_{3},q_{p1},q_{p2},d),$ (14c)

where $d$ is the vector containing the demands of all the nodes,
$q_{p1},q_{p2}$ are the pump flows, $A_{1},A_{2},A_{3}$ are the cross
sectional areas of the tanks and $f_{1},f_{2},f_{3}$ are nonlinear flow
functions. Water levels at the two connected tanks are almost equal
$h_{1}\thickapprox h_{2}$ all the time since the pipe connecting respective
tanks is big enough to oppose the water flows coming from neighbor nodes. That
enables us to consider $h_{1},h_{2}$ together as

$(A_{1}+A_{2})\dot{{h}}_{1,2}\thickapprox
q_{1down}+q_{2down}+q_{up}=f_{1}+f_{2}.$ (15)

We have used the EPANET model of the network to generate the data required for
approximating $f_{1}+f_{2}$ and $f_{3}$. The model is simulated with various
tank level initial conditions and flow rates of 2 pumping stations to be
controlled. The control laws for the remaining pumping stations are already
defined in the EPANET model. Then, the linear model (4) is fitted to
simulation data using least squares. The state variables for the model are
$h(t)=[{h}_{1,2}(t),{h}_{3}(t)]\in\mathbb{R}^{2}$ and the inputs are
$u(t)=[{q}_{p1}(t),{q}_{p2}(t)]\in\mathbb{R}^{2}$. The total demand of High
and Low Zone is used as aggregated demand $d_{a}$ in the model since mainly
those areas are supplied by the controlled pumps.

###  4.1 Simulation Results

The proposed control method is tested on EPANET model of Randers water
network. Epanet-Matlab toolkit Eliades et al. (2016) is used to set the flow
of the 2 pumps at each time step and simulate the network. The remaining pumps
are controlled with rule-based control laws that are previously defined on
EPANET.

The parameters of exponential barrier functions $J_{hi}$ are chosen as
$a_{i}=80,$ $b_{i}=0.3$ for all $i$. It is assumed that the electricity prices
are known in advance during the test. Tank levels $h_{1},h_{2}$ have a maximum
value of 3m while $h_{3}$ has 2.8m. Tanks are required to be at least half
full. Maximum pump flows are set to 100. Sampling time $\Delta_{t}$ is set to
1 hour in the experiments, so the control input is recalculated at each hour.
We assume that total demand $d_{a}(t)$ of HZ and LZ can be estimated up to 1
day from available data. Although we do not have historical data on the
demand, we imitate this behaviour by using a slightly perturbed version of the
real demand used in EPANET simulation during MPC calculations. The
perturbations are adapted from a real demand data set of a small Danish
facility. Normalized difference between the average demand and the demand of a
random day in data set is added to EPANET demand to replicate estimated
demand. In each experiment a different day from the data set is used, so the
assumed estimated demand is different each time.

The simulation results when the presented method is applied to the EPANET
model are given in Figure 4. The initial tank levels are equal to
$h_{T_{day}/\Delta_{t}}^{*}$ in the simulation. The top plot shows the
evolution of tank levels along with the upper and lower thresholds. It is seen
that the thresholds are not violated and moreover tank levels are not getting
too close to them, which was the idea behind exponential barrier functions.
Both the real demand and the assumed estimated demand of HZ and LZ are in the
figure below. Total applied pump flows and electricity prices are in the
following figures. The expected result is pump flows being higher when
electricity prices are low, and lower when they are high, which seems to be
the case as can be seen in the plot. Pump flows drop significantly when prices
are at the peak and they reach their highest value at the end of the day when
prices are low. A more aggressive controller can be obtained by picking a
smaller $b_{i}$ value for barrier functions at the expense of risking
constraint violation.

(a)

(b)

(c)

(d)

Figure 4: Sample simulation. (a) evolution of tank levels through 1 day with
upper and lower level thresholds; (b) real total demand of HZ and LZ used in
EPANET simulation and the demand used in MPC calculations; (c) total flow
provided by the 2 pumps; (d) electricity price.

In Figure 5, the tank level simulation results and control inputs for
different initial conditions and different assumed estimated demands are
given. The electricity price profile is the same as before. It is seen that
the algorithm is able to control the network on various cases while satisfying
the constraints. Regardless of initial tank levels, the pumping profiles have
a similar profile: high pump flows close to midnight and in the middle of the
day. The only exception is the bottom plot. In the beginning, prices are low
but pump flows are not high. This can be attributed to water levels
$h_{1},h_{2}$ being close to the upper thresholds and water demand being low
in the beginning.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 5: Tank levels and pump flows for different initial conditions

The assumption that the optimal input sequences $U(t)$ would not diverge a lot
from the one found in previous step $U(t-\Delta_{t})$ is the reason we apply
$u_{1}^{t-\Delta_{t}}$ at time $t$ if the problem (5) is infeasible at time
$t$. This assumption is tested with initial conditions
$h_{1,2,3}=h_{T_{day}/\Delta_{t}}^{*}$. In figure 6, total pump flow
$[1,1]^{T}u_{i}^{t},i=0\cdots N(t)-1$ of the found optimal input sequences
$U(t),t=0,\Delta_{t}\cdots T_{day}-\Delta_{t}$, except when the problem were
infeasible, are given. It can be seen that $u_{1}^{t-\Delta_{t}}$ is close to
the $u_{0}^{t}$ for all $t$, which shows that our assumption is valid at least
for this experiment.

Figure 6: Evolution of found input sequences $U(t)$ through 1 day. It can be
seen that the solutions remain close to the initial optimal sequence $U(0)$.

Finally, the ability of the algorithm to decrease economic costs is tested
with various initial conditions. For each case, a demand follower pumping
strategy is used as a benchmark. The flow of the 2 pumps is adjusted with
trial and error for each demand follower such that the total flow of the 2
pumps is equal to water demand at each time step and tank levels satisfy the
terminal constraint (13a). The demand follower is a natural candidate to be a
benchmark method since providing as much water as demand is an intuitive idea
and the constraints in (5) can be satisfied with manual adjustments of pump
flows. The economic costs are presented relatively in Table 1 As it is seen,
the proposed algorithm saves between $40\%$ and $45\%$ of the cost with
different demand profiles.

Proposed Method | Demand Follower
---|---
0.5967 | 1
0.5745 | 1
0.5826 | 1
0.5558 | 1
Table 1: Relative economic costs of the proposed method and demand follower
strategy for various demand profiles

##  5 Conclusion

We have presented a predictive control algorithm with a periodic horizon for
WDNs. The aim is to minimize the economic cost and satisfy the operational
constraints. A linear model is used to represent Randers WDN to increase the
chance of finding a solution to the problem (5) at expense of a model-plant
mismatch. Periodic horizon is introduced to the predictive control formulation
to keep the resulting state trajectories around the optimal periodic
trajectory. Barrier functions are used to prevent constraint violation since
there is a model-plant mismatch.

The presented algorithm is tested on Randers WDN using EPANET. It is shown in
various situations that the method is able to find an economic solution where
pump flows are adjusted according to electricity prices. Also, it is shown
that the system trajectories do not enter dangerous zones introduced by
barrier functions as long as the predicted demand and the actual demand are
somewhat close.

As future work, we plan to work on theoretical guarantees of the existence of
solutions to the proposed method. Also, the robustness of periodic horizon
control of periodical systems with barrier functions will be investigated.

## References

  * Abdelsalam and Gabbar (2021) Abdelsalam, A.A. and Gabbar, H.A. (2021).  Energy saving and management of water pumping networks.  _Heliyon_ , 7(8), e07820.  https://doi.org/10.1016/j.heliyon.2021.e07820.
  * Bagirov et al. (2013) Bagirov, A.M., Barton, A., Mala-Jetmarova, H., Nuaimat, A.A., Ahmed, S.T., Sultanova, N., and Yearwood, J. (2013).  An algorithm for minimization of pumping costs in water distribution systems using a novel approach to pump scheduling.  _Math. Comput. Model._ , 57, 873–886.
  * Baunsgaard et al. (2016) Baunsgaard, K.M.H., Ravn, O., Kallesøe, C.S., and Poulsen, N.K. (2016).  Mpc control of water supply networks.  _2016 European Control Conference (ECC)_ , 1770–1775.
  * Berkel et al. (2018) Berkel, F., Caba, S., Bleich, J., and Liu, S. (2018).  A modeling and distributed mpc approach for water distribution networks.  _Control Engineering Practice_.
  * Castro-Gama et al. (2017) Castro-Gama, M.E., Pan, Q., Lanfranchi, E.A., Jonoski, A., and Solomatine, D.P. (2017).  Pump scheduling for a large water distribution network. milan, italy.  _Procedia Engineering_ , 186, 436–443.
  * Eliades et al. (2016) Eliades, D.G., Kyriakou, M., Vrachimis, S., and Polycarpou, M.M. (2016).  Epanet-matlab toolkit: An open-source software for interfacing epanet with matlab.  In _Proc. 14th International Conference on Computing and Control for the Water Industry (CCWI)_ , 8. The Netherlands.  10.5281/zenodo.831493.
  * Fiedler et al. (2020) Fiedler, F., Cominola, A., and Lucia, S. (2020).  Economic nonlinear predictive control of water distribution networks based on surrogate modeling and automatic clustering.  _IFAC-PapersOnLine_ , 53, 16636–16643.
  * Ghaddar et al. (2015) Ghaddar, B., Naoum-Sawaya, J., Kishimoto, A., Taheri, N., and Eck, B. (2015).  A lagrangian decomposition approach for the pump scheduling problem in water networks.  _Eur. J. Oper. Res._ , 241, 490–501.
  * Kallesøe et al. (2017) Kallesøe, C.S., Jensen, T.N., and Bendtsen, J.D. (2017).  Plug-and-play model predictive control for water supply networks with storage.  _IFAC-PapersOnLine_ , 50, 6582–6587.
  * Lindell Ormsbee (2009) Lindell Ormsbee, Srini Lingireddy, D.C. (2009).  Optimal pump scheduling for water distribution systems.  URL http://www.uky.edu/WDST/PDFs/[73.3]%20Ormsbee%20Optimal%20Pump%20Scheduling%20Paper.pdf.
  * Pour et al. (2019) Pour, F.K., Puig, V., and Cembraño, G. (2019).  Economic mpc-lpv control for the operational management of water distribution networks.  _IFAC-PapersOnLine_.
  * Sharif et al. (2019) Sharif, N., Haider, H., Farahat, A., Hewage, K., and Sadiq, R. (2019).  Water energy nexus for water distribution systems: A literature review.  _Environmental Reviews_ , 27.  10.1139/er-2018-0106.
  * Wang et al. (2018) Wang, Y., Alamo, T., Puig, V., and Cembraño, G. (2018).  Economic model predictive control with nonlinear constraint relaxation for the operational management of water distribution networks.  _Energies_ , 11, 991.
  * Wang et al. (2021) Wang, Y., Yok, K.T., Wu, W., Simpson, A.R., Weyer, E., and Manzie, C. (2021).  Minimizing pumping energy cost in real-time operations of water distribution systems using economic model predictive control.  _ArXiv_ , abs/2010.07477.