\section{Introduction}
A central tool in combinatorial and geometric group theory is the notion of the Cayley graph of a group $\Gamma$, which is a graph that encodes the structure of the group, see Section~\ref{Sec:1} for all the relevant definitions.
This construction is enlightening in two ways: on one hand it allows to see groups as geometric objects, while on the other hand it endows some graphs with an algebraic structure.
In both direction, it permits a better understanding of the objects under consideration.
It is therefore natural to ask which graphs can be realized as a Cayley graph of some group.
This question has been answered in 1958 by Sabidussi~\cite{MR0097068} who showed that a connected graph $G$ is isomorphic to a Cayley graph if and only if $\Aut(G)$ acts freely and transitively on the vertices of $G$.

It follows from Sabidussi's result that Cayley graphs are rather sparse among graphs, and thus too restrictive to study general graphs, even if we restrict ourself to the specific case of connected regular degree graphs. A more general notion of a ``graph of algebraic origin'' is hence needed.
A suitable generalization of Cayley graphs, is the concept of a Schreier graph associated to a pair $(\Gamma,\Lambda)$ where $\Lambda$ is a subgroup of $\Gamma$ (or equivalently associated to a transitive action of $\Gamma$ on some set $X$). In this context, Cayley graphs correspond to pairs $(\Gamma,\{1\})$ (equivalently to free and transitive actions).
Schreier graphs are always connected and of regular degree and one can thus wonder which connected graphs of regular degree are isomorphic to a Schreier graph.

The first step in this direction was made by Gross in his 1977 article~\cite{MR0450121}, where he proved that every finite, connected regular graph of even degree is isomorphic to a Schreier graph.
The result for infinite connected regular graphs of even degree follows by a compacity argument and was probably one of this well-known folklore result for long.
A written proof of the infinite version can be found in~\cite{MR2893544}, but it has since be rediscovered many times.
On the other hand, it is easy to find examples, even finite ones, of connected regular graphs of odd degree that are not isomorphic to a Schreier graph, see the left of Figure~\ref{Figure:Simple3reg}.
It hence remains to classify which connected regular graphs of odd degree are isomorphic to a Schreier graph.
This question turns out to be equivalent to the existence of a perfect matching, an interesting subject in its own right, and also to the existence of a covering map to some particular graph, see Section~\ref{Sec:2}.
An interesting result in this direction is the fact that all transitive graphs (with at most countably many vertices and edges) are isomorphic to a Schreier graph.
This fact was first proved by Godsil and Royle for finite graphs~\cite{GodsilRoyle}, then extended to locally finite graphs by the author~\cite{LeemannThese}, and finally to countable graphs by Hamann and Wendland in the appendix of~\cite{2020arXiv200706432G}.

Even if some connected regular graphs of odd degree are not isomorphic to a Schreier graph, one can still ask if they roughly ``look like'' a Schreier graph.
One possible interpretation of the above is to ask for the existence of a finite degree covering from a Schreier graph to $G$.
The aim of this note is to provide a proof of the following result, which concludes the above discussion.
\begin{proposition}\label{Prop:Schreier}
Let $G$ be a $d$- regular connected graph.
Then either $G$ is isomorphic to a Schreier graph or $G$ has a double-cover $H$ which is isomorphic to a Schreier graph.
\end{proposition}
While we were not able to find a reference to the above result in the literature, we do not claim any priority on it.
In fact, this note was inspired by the following remark, which can be found at the end of Section 7 of~\cite{MR2087782}: ``In fact up to the cover of degree $2$ any regular graph can be realized as a Schreier graph~\cite{MR1308046}''.
However,~\cite{MR1308046} seems to treat only the case of graphs of even degree.

\paragraph{}
Section~\ref{Sec:1} contains all the definitions as well as a discussion on the unusual concept of degenerated loop. The short Section~\ref{Sec:2} contains more materials on coverings and perfect matchings, as well as the proof of Proposition~\ref{Prop:Schreier}.

\paragraph{Acknowledgements}
The author is thankful to A. Georgakopoulos, R. Grigorchuk and M. de la Salle for comments on a previous version of this note.
\section{Preliminary definitions and results}\label{Sec:1}
Since one of our motivations is the study of Schreier graphs, we will look at graphs in a broad sense.
We allow multi-edges and loops as well as the less usual degenerated loops (also called half-edges, dangling edges or open edges).

Informally, a degenerated loop is the index $2$ quotient of a path of length~$1$, see Figure~\ref{Figure:DoubleCovers}.
Heuristicaly it can either be seen as an edge with only one end, or as a loop that add only one to the degree.
Allowing degenerated loops is necessary to treat the case of Schreier graphs of groups with generators of order~$2$.
\begin{definition}
A \emph{graph} is a $5$-tuple $(V,\vec E,\iota,\tau,\inv)$ where $V$ and $\vec E$ are sets (of \emph{vertices} and of \emph{arcs}) and $\iota,\tau\colon \vec E\to V$ and $\inv\colon \vec E\to\vec E$ are applications (\emph{initial vertex}, \emph{terminal vertex} and \emph{inverse arc}) such that $\iota\cdot\inv=\tau$, $\tau\cdot\inv=\iota$ and $\inv$ is an involution.
\end{definition}
There is a natural equivalence relation on arcs defined by $\vec e\sim\vec f$ if and only if $\vec e=\vec f$ or $\vec e=\inv(\vec f)$.
The quotient set $\vec E/\sim$ is denoted by $E$ and is the set of \emph{edges}.

A vertex $v$ and an edge $e=\{\vec e,\inv(\vec e)\}$ are \emph{adjacent} if $v=\iota(\vec e)$ or $v=\tau(\vec e)$, in which case we also say the $v$ is an \emph{end} of $e$.
An edge that is adjacent to only one vertex is called a \emph{loop}.
A loop $e=\{\vec e,\inv(\vec e)\}$ is \emph{degenerated} (or an \emph{half-edge}) if $\vec e=\inv(\vec e)$, that is if $\vec e$ is a fixed point of $\inv$.

\begin{remark}
Graph theorists usually forbid the existence of degenerated loops, but these special kinds of edges naturally arise in the context of Schreier graphs.
Indeed, when allowing degenerated loops, there is a bijection between the subgroups of $G$ (up to conjugacy) and graphs that are covered by $\Cayley{\Gamma}{S}$, see~\cite{MR3463202} for a proof.
\end{remark}

The \emph{degree} of a vertex is the number of arcs $\vec e$ such that $v=\iota(\vec e)$.
In particular, a non degenerated loop adds $2$ to the degree of its adjacent vertex, while a degenerated loop adds only $1$.
All graphs under consideration will be \emph{locally finite}, that is the degree of any vertex is finite.
Classical notions like \emph{paths} and \emph{connected components} are defined in a natural way.
In particular, a \emph{bipartite graph} has no loop.
Finally, for a vertex $v$ we define $\Star(v)\coloneqq\{\vec e\in \vec E\ |\ \iota(\vec e)=v\}$.

On pictures, we will represent degenerated loops  by dotted line and other edges (including non-degenerated loops) by plain lines.
\begin{figure}[htbp]
\centering
\begin{tikzpicture}[node distance=1cm,every node/.style=mynode]
\node (00){};
\node[fill=black] (01)[right of=00]{};
\node[fill=black] (10)[below of=00]{};
\node (11)[right of=10]{};
\node (v0)[below =1.5cm of 10]{};
\node[fill=black] (v1)[below =1.5cm of 11]{};
\graph[edge={thick}]{(00)--(01)--(11)--(10)--(00)};
\graph[edge={thick,bend left}]{(v0)--(v1)--(v0)};
\draw[->>,thick] (0.5,-1.3) -- (0.5,-2.1);
\begin{scope}[xshift=3.5cm]
\node (0){};
\node (1)[below of=0]{};
\node (v)[below =1.5cm of 1]{};
\graph[edge={thick, bend left}]{(0)--(1)--(0)};
\path	(v) edge[myloop below](v);
\draw[->>,thick] (0,-1.3) -- (0,-2.1);
\begin{scope}[xshift=2.5cm]
\node (0){};
\node (1)[below of=0]{};
\node (v)[below =1.5cm of 1]{};
\graph[edge={thick}]{(0)--(1)};
\path	(v) edge[myloop below, densely dotted,thick](v);
\draw[->>,thick] (0,-1.3) -- (0,-2.1);
\end{scope}
\end{scope}
\end{tikzpicture}
\caption{Double covers of some graphs; the middle and the right ones are canonical double covers.
The leftmost one is isomorphic to the double cover $\Cayley{\Z/4\Z}{\{\pm1\}}\twoheadrightarrow\Schreier{\Z/4\Z}{\Z/2\Z}{\{\pm1\}}$, the middle one to $\Schreier{\Z/4\Z}{\Z/2\Z}{\{\pm1\}}\twoheadrightarrow\Schreier{\Z/4\Z}{\Z/4\Z}{\{\pm1\}}$ and the rightmost one to $\Cayley{\Z/2\Z}{\{\pm1\}}\twoheadrightarrow\Schreier{\Z/2\Z}{\Z/2\Z}{\{\pm1\}}$.}
\label{Figure:DoubleCovers}
\end{figure}
A \emph{morphism} of graph from $(V,\vec E)$ to $(W,\vec F)$ is a map $\varphi\colon V\sqcup\vec E\to W\sqcup\vec F$ which is compatible with the graph structure.
For any vertex $v\in V$, this induces a map $\varphi_v\colon\Star(v)\to\Star\bigl(\varphi(v)\bigr)$.
A morphism $\varphi\colon G\to H$ between two graphs is a \emph{covering} if all the induced map $\varphi_v$ are bijections.
If moreover $\lvert\varphi^{-1}(v)\rvert=2$\footnote{The quantity $\lvert\varphi^{-1}(v)\rvert$ depends only on the connected component of $v$.} for every vertex, then it is called a \emph{double cover}.
It is an easy exercise to show that every $2d$-regular graph without degenerated loops covers a graph with $1$ vertex and $d$ loops.

Given a graph $G=(V,\vec E)$, its \emph{canonical double cover} is the tensor product $G\otimes K_2$ where $K_2$ has two vertices and an edge between them.
That is, the vertex set of $G\otimes K_2$ is $V\times \{0,1\}$, and an arc $\vec e$ from $v$ to $w$ is lifted to an arc from $(v,0)$ to $(w,1)$ and to an arc from $(v,1)$ to $(w,0)$.
In particular, a non-degenerated loop adjacent to $v$ is lifted to a pair of parallel edges between $(v,0)$ and $(v,1)$, while a degenerated loop adjacent to $v$ is lifted to a single edge between $(v,0)$ and $(v,1)$.
See Figure~\ref{Figure:DoubleCovers}.
A simple verification shows that $G\otimes K_2$ is a bipartite graph, which is a double cover of $G$. Moreover, $G\otimes K_2$ is connected if and only if $G$ is connected and not bipartite.

\begin{definition}
A \emph{matching} (also called a \emph{$1$-factor}) of a graph $G$ is a subgraph $M$ of $G$ such that every vertex of $M$ has degree $1$ in $M$.
A \emph{perfect matching} is a matching of $G$ that contains all vertices of $G$.
A graph is \emph{matchable} if it admits a perfect matching.
Two matchings $M_1$ and $M_2$ are \emph{orthogonal} if they share no edge.
\end{definition}
Observe that not every regular graphs of odd degree are matchable, even among graphs without degenerated loop.
See Figure~\ref{Figure:Simple3reg}.
However, every vertex-transitive graph of odd degree is matchable,~\cite{GodsilRoyle,LeemannThese}.
\begin{figure}[htbp]
\centering
\begin{tikzpicture}[every node/.style={mynode},scale=0.8]
\node(A) {};
\foreach \k in {0, 1, 2}{
	\begin{scope}[rotate=\k*120]
	\node[fill=white] (B\k) at (1,0){};
	\node (C\k) at (2,1){};
	\node (D\k) at (2,-1){};
	\node (E\k) at (2,0){};
	\node (F\k) at (3,0){};
	\end{scope}
	\graph[edge={thick}]{
	(A)--(B\k)--(C\k)--(F\k)--(D\k)--(B\k);
	(C\k)--(E\k)--(D\k);
	(E\k)--(F\k)
	};
}
\graph[edge={thick}]{(A)--(B0)};
\node (b0) at (1,0){};
\begin{scope}[xshift=7cm]
\node(A) {};
\foreach \k in {0, 1, 2}{
	\node (B\k) at (\k*120:1){};
	\node (C\k) at (\k*120-20:2){};
	\node (D\k) at (\k*120+20:2){};
	\graph[edge=thick]{
	(B\k)--(C\k)--(D\k)--(B\k)};
	\path (B\k) edge[loop, out=\k*120+160, in = \k*120+80,min distance =6mm, densely dotted,thick](B\k);
	\path (C\k) edge[loop, out=\k*120-20, in = \k*120-100,min distance =6mm, thick](C\k);
	\path (D\k) edge[loop, out=\k*120+100, in = \k*120+20,min distance =6mm,thick](D\k);
}
\graph[edge=thick]{
(B0)--(A)--(B1);
(B2)--(A);
};
\path (A) edge[loop, out=100, in = 20,min distance =6mm, densely dotted, thick](A);\end{scope}
\end{tikzpicture}
\caption{A simple $3$-regular graph (on the left) and a $4$-regular graph (on the right); both of them being not isomorphic to a Schreier graph.}
\label{Figure:Simple3reg}
\end{figure}

We finally recall the definition of Schreier graphs.
\begin{definition}
Let $\Gamma$ be a group with a symmetric generating set $S$ not containing the identity and let $\Lambda$ be a subgroup of $\Gamma$.
The corresponding (right) \emph{Schreier graph} $\Schreier{\Gamma}{\Lambda}{S}$ is the graph with vertex set the right cosets $\{\Lambda g\ |\ g\in\Gamma\}$ and with an arc from $\Lambda g$ to $\Lambda h$ for every $s$ in $S$ such that $\Lambda gs=\Lambda h$.
If $X$ is a set with a transitive right action $X\curvearrowleft \Gamma$, the corresponding (right) \emph{orbital graph} $\SchreierOr{\Gamma}{X}{S}$ is the graph with vertex set $X$ with an arc from $x$ to $y$ for every $s$ in $S$ such that $x.s=y$.
\end{definition}

These two definitions are the two faces of the same coin. Indeed, we have $\SchreierOr{\Gamma}{X}{S}=\Schreier{\Gamma}{\Stab_G(x)}{S}$ for every $x\in X$ and $\Schreier{\Gamma}{\Lambda}{S}=\SchreierOr{\Gamma}{\Lambda \backslash\Gamma}{S}$.
See Figure~\ref{Figure:DoubleCovers} for some examples.
Schreier graphs are generalizations of the so-called \emph{Cayley graphs}, where $\Cayley{\Gamma}{S}\coloneqq\Schreier{\Gamma}{\{1\}}{S}$.

In the graph $\Schreier{G}{H}{S}$ there is a loop at the vertex $Hg$ if and only if there exists $s$ in $S\cap g^{-1}Hg$.
This loop is degenerated if and only if $s$ is of order $2$.
If $G$ is a (finitely generated) group, then it admits a (finite) generating system without element of order $2$ if and only if it is not generalized dihedral~\cite{MR0498225}.
Moreover, even for non-generalized dihedral group, it is sometimes natural to look at generating set containing elements of order $2$.
This is the main reason why we allow degenerated loops in our definition of a graph.

\section{Coverings and matchability}\label{Sec:2}
First of all we have the following, which relates covers and perfect matchings:
\begin{lemma}\label{lemma:CoversMatch}
Let $G$ be a $(2d+1)$-regular graph without degenerated loop.
Then $G$ is matchable if and only if it covers a graph with one vertex, one degenerated loop and $d$ non-degenerated loops.
\end{lemma}
\begin{proof}
If $\varphi\colon G\to R$ is a covering and $R$ has a degenerated loop $e$, then $\varphi^{-1}(e)$ is a perfect matching of $G$.

On the other hand, let $M$ be a perfect matching of $G$ and $G'$ the graph obtained by removing all edges of $M$ from $G$.
Then $G'$ is $2d$-regular without degenerated loop and hence covers $R'$ a graph with one vertex and $d$ non-degenerated loops.
Let $R$ be the graph obtained by adding a degenerated loop to $R'$.
The covering $G'\twoheadrightarrow R'$ naturally extends to $G\twoheadrightarrow R$ by sending all edges of $M$ to the degenerated loop of $R$.
\end{proof}
As an easy generalization we obtain
\begin{corollary}
Let $G$ be a $(2d+n)$-regular graph without degenerated loop.
Then $G$ admits $n$ pairwise orthogonal perfect matchings if and only if it covers a graph with one vertex, $n$ degenerated loop and $d$ non-degenerated loops.
\end{corollary}

By a theorem of Gross~\cite{MR0450121}, every finite $2d$-regular connected graph without degenerated loop is isomorphic to a Schreier graph of the free group~$F_d$.
This result extends by compacity to infinite $2d$-regular connected graph without degenerated loop, see~\cite{MR2893544} for a proof.
On the other hand, Gross' Theorem does not carry over graphs of odd regular degree, see Figure~\ref{Figure:Simple3reg}.
However, every connected transitive graph of degree $2d+1$ and without degenerated loop is isomorphic to a Schreier graph of $F_d*(\Z/2\Z)$,~\cite{GodsilRoyle,LeemannThese}.
For general connected regular graphs, the situation is more complicated.
\begin{lemma}\label{Lemma:Isomorphism}
Let $G$ be a connected $(2d+1)$-regular graph without degenerated loop.
Then $G$ is isomorphic to a Schreier graph if and only if it is matchable.
In this case it is isomorphic to a Schreier graph of $F_d*(\Z/2\Z)$, where $F_d$ is the free group of rank $d$.
\end{lemma}
A proof of this fact can be found in~\cite{LeemannThese}, but the result itself was probably already general knowledge at that time and essentially (at least for finite graphs) follows from~\cite{MR0450121}.
A variation of Lemma~\ref{Lemma:Isomorphism} gives us
\begin{proposition}[\cite{LeemannThese}]\label{Prop:Matchings}
Let $G$ be a connected $(2d+n)$-regular graph without degenerated loop.
Then $G$ is isomorphic to a Schreier graph of $F_d*(\Z/2\Z)^{*n}$ if and only if it has $n$ perfect matchings that are pairwise orthogonal.
\end{proposition}
In order to prove Proposition~\ref{Prop:Schreier}, we will need the following fact, which follows (for example) from K\"onig's duality theorem for infinite bipartite graphs,~\cite{MR734985}.
It also follows, by compacity, from the similar statement for finite graphs.
\begin{lemma}
Every $d$-regular bipartite graph is matchable.
\end{lemma}
We now prove a more precise version of Proposition~\ref{Prop:Schreier}.
\begin{lemma}\label{Prop:Main}
Let $G$ be a $d$-regular connected graph.
If $G$ is bipartite, then it is isomorphic to a Schreier graph of $(\Z/2\Z)^{*d}$.
Otherwise, $G\otimes K_2$ is isomorphic to a Schreier graph of $(\Z/2\Z)^{*d}$.
\end{lemma}
\begin{proof}
If $G=(V,\vec E)$ is bipartite let $K\coloneqq G$, while if $G$ is not bipartite take $K$ to be its canonical double cover.

Now, the graph $K$ is bipartite and $d$-regular.
Let $M_1$ be a perfect matching of $K$ and let $K_1$ be the graph obtained from $K$ by removing the edges in $M_1$.
Then $K_1$ is bipartite of degree $d-1$ and hence has a perfect matching $M_2$ as soon as $d-1>0$.
It is trivial that $M_1$ and $M_2$ are orthogonal.
By repeating this we obtain $d$ pairwise orthogonal perfect matchings. By Proposition~\ref{Prop:Matchings} this implies that $K$ is isomorphic to a Schreier graph of $(\Z/2\Z)^{*d}$.
\end{proof}
\begin{corollary}
Let $G$ be a connected regular graph.
\begin{enumerate}
\item If $G$ has degree $2d$, then either $G$ (if it is bipartite) or $G\otimes K_2$ (otherwise) is isomorphic to a Schreier graph of $F_d$,
\item If $G$ has degree $2d+1$, then either $G$ (if it is bipartite) or $G\otimes K_2$ (otherwise) is isomorphic to a Schreier graph of $F_d*(\Z/2\Z)$.
\end{enumerate}
\end{corollary}
\begin{proof}
This follows from Lemma~\ref{Prop:Main} and the fact that a Schreier graph of $(\Z/2\Z)^{*2d}$ (respectively of $(\Z/2\Z)^{*(2d+1)}$) without degenerated loop is isomorphic to a Schreier graph of $F_d$ (respectively of $F_d*(\Z/2\Z)$),~\cite{LeemannThese}.
\end{proof}
Observe that a Schreier graph is necessarily connected and of regular degree.
On the other hand, there are examples of connected regular graphs of odd degree and without degenerated loop that are not isomorphic to a Schreier graph, see the left of Figure~\ref{Figure:Simple3reg}.
There are also examples of connected regular graphs of even degree (but with at least one degenerated loop) that are not isomorphic to a Schreier graph, see the right of Figure~\ref{Figure:Simple3reg}.
Finally, it isn't possible to ensure that $G$ has a double cover isomorphic to a Schreier graph.
Indeed, if $T$ is a regular tree, then its only double cover is the disjoint union of two copies of $T$, which is not connected and hence not isomorphic to a Schreier graph.
Hence, Proposition~\ref{Prop:Schreier} is the ``best possible'' for locally finite graphs.
\providecommand{\noopsort}[1]{} \def$'${$'$}


\section{Introduksyon}

Sa mga nagdaang taon, mabilis ang pag-unlad ng Natural Language Processing (NLP) dahil sa kasaganahan ng data dulot ng adopsyon ng internet sa buong mundo at sa pagiging accessible ng mas murang computing power. Nagkaron ng pag-unlad sa mga aplikasyon ng NLP tulad ng Machine Translation, Text Categorization, Text Classification, at iba pa.

Karamihan sa mga tagumpay ng NLP ay para sa mga mainstream na wika tulad ng Ingles at iba pang wika na merong accessible na malalaking text corpora at annotated texts. Ang mga wikang wala o limitado ang access sa malaking text corpora o annotated texts ay tinatawag na low-resource languages.

Ang pangangailangan ng malalaking text corpora o datasets para mapakinabangan ang mga benepisyo ng NLP ay isang balakid na pumipigil o nagpapabagal ng adopsyon ng teknolohiya sa mga low-resource languages. Sa papel na ito, pagtutuunan ng pansin ang paggamit ng Transfer Learning (TL) upang malagpasan ang balakid ng low-resource language tulad ng wikang Filipino.

Sa mga nakaraang taon, nanaig ang mga transformer-based transfer learning techniques pagdating sa pagiging epektibo sa mga iba't ibang aplikasyon ng NLP sa low-resource languages tulad ng Fake News Detection in Filipino \cite{FAKENEWS}, Named Entity Recognition in Vietnamese \cite{PHOBERT}, at iba pa,

Bagamat ang transformer-based models ay napatunayang epektibo sa iba't ibang aplikasyon ng NLP, ito ay hindi ganoon ka-accessible sa kasalukuyang panahon dahil mataas ang kailangang compute power at memory upang ma-train ito. Habang ang RNN-based models naman kagaya ng AWD-LSTM \cite{MERITY} ay mas higit na accessible o mura i-train kumpara sa transformer-based models. Ang accessibility ng AWD-LSTM ay mahalaga para mapalawak ang adopsyon ng NLP sa Pilipinas at mas mapakinabangan ng nakararami ang benepisyo nito.

Ang papel na ito ay may tatlong kontribusyon: 1) Maglabas ng pre-trained AWD-LSTM language model (LM) sa wikang Filipino upang maging tuntungan sa pagbuo ng mga NLP applications sa wikang Filipino\footnote{https://github.com/danjohnvelasco/Filipino-ULMFiT}. 2) Mag benchmark ng AWD-LSTM sa Hate Speech classification task at ipakita na kayang makipagsabayan nito sa mga transformer-based models. 3) Pangatlo, suriin ang performance ng AWD-LSTM sa low-resource setting gamit ang degradation test at ikumpara ito sa mga transformer-based models.

\section{Background}
\subsection{Natural Language Processing}
Ang Natural Language Processing (NLP) ay isang subfield ng linguistics, computer science, at artificial intelligence na nauukol sa pag proseso at pag-unawa ng natural na wika \cite{WIKI:NLP}. Ang ilan sa mga aplikasyon ng NLP ay ang email spam filters (Text Classification), pag-unawa ng nais sabihin tulad ng mga smart assistants (Language Understanding), pagsasalin ng isang wika sa iba pang wika (Machine Translation), mag predict ng susunod na salita base sa mga naunang salita (Language Modelling), at marami pang iba. Dahil sa kaunlaran sa kasaganahan sa datos at pagiging accessible ng malakas na compute power, nabuhay muli ang machine learning approach. Sa maikling salita, ang machine learning approach ay gumagamit ng malaking datos na ginagamit ng isang computer algorithm upang matutunan ang mga patterns ng datos na ito. Dahil dito, naging epektibo siyang approach sa mga komplikadong problema dahil hindi na kailangan direktang i-program ang mga rules para malutas ang isang problema.

\subsection{Transfer Learning}
Notorious ang machine learning approach sa pangangailangan nito ng sobrang laking datos para mapakinabangan. Ang Transfer Learning (TL) ay isang area ng research na concerned sa problemang ito \cite{WIKI:TL}. Sa maikling salita, ang TL ay ang pag retain o pagpapanatili ng mga natutunan ng isang model sa isang gawain at paggamit o "transfer" ng mga natutunan nito sa iba pero may kaugnayan na gawain. Halimbawa, ang mga natutunan ng isang model sa pag detect ng muka ng tao ay maaring gamitin bilang tuntungan para sa pag-aaral ng model na matutunan kung ang muka ng tao ay galit, masaya, at iba pang facial expressions \cite{LI}.

\section{Metodolohiya}
Ang metodolohiya ay ang mga sumusunod 1) Mag train ng model gamit ang ULMFiT \cite{ULMFIT}, isang RNN-based transer learning method. 2) Sukatin ang kahusayan ng model sa text classification task gamit ang Hate Speech Dataset \cite{CABASAG}. Gamit ang parehong dataset, magsasagawa ng degradation test upang masuri ang performance ng AWD-LSTM sa low-resource setting.\cite{CRUZCHENG:20}. Makikita ang kabuuang proseso ng ULMFiT sa Figure \ref{fig:ulmfit_summary}.

Ang API na ginamit sa papel na ito ay fastai v2.0.13\footnote{https://pypi.org/project/fastai/2.0.13/}. Kapag hindi binanggit ang isang partikular na configuration sa model, ibig sabihin ang default setting lang ang ginamit. Para mas mapabilis ang training, gumamit ng mixed precision training \cite{MICIKEVICIUS}. Ang GPU na ginamit sa training ay Tesla T4. Sa buong proseso ng model training, ang learning rate schedule na 1cycle policy \cite{LESLIESMITH} ang ginamit. Ang buong code ay available sa public repository\footnote{https://github.com/danjohnvelasco/Filipino-ULMFiT}.

\begin{figure*}[t]
\includegraphics[width=16cm]{ulmfit}
\caption{ULMFiT approach summary. Adapted from fastai documentation\footnote{https://docs.fast.ai/tutorial.text}. Figure was modified to match the dataset used in this paper.}
\label{fig:ulmfit_summary}
\centering
\end{figure*}

\subsection{ULMFiT}
Ang ULMFiT o Universal Language Model Fine-tuning \cite{ULMFIT} ay isang epektibong transfer learning technique na gumagamit ng language model na natuto sa malaking unlabeled text corpora at gagamitin ito bilang tuntungan sa iba pang gawain. Napatunayan na epektibo ang approach na ito kahit na maliit lang ang target text corpus o datos na gagamitin para sa partikular na gawain. Ang technique na ito ay may tatlong hakbang: 1) LM Pretraining Phase o pag pretrain ng AWD-LSTM language model sa isang malaking unlabeled text corpus. 2) LM Fine-tuning Phase o ang paggamit ng pre-trained LM bilang tuntungan at sanayin o i-train pa ng husto ang model sa target text corpus. 3) Text Classifier Fine-tuning Phase o ang pagsanay ng LM fine-tuned model sa text classification task.

\subsubsection{Language Model Pre-training}
Gumagamit ito ng AWD-LSTM \cite{MERITY} at sasanayin ito sa language modelling task na kung saan base sa isang pagkasunod sunod na mga salita, magbibigay ka ng prediksyon kung anong salita ang may pinakamataas na probability na kasunod nito. Sa LM pre-training phase, kailangan ng isang malaking unlabeled text corpora upang matutunan ng model ang mga patterns at structure ng wikang ginagamit sa text corpora. Ang ideal na text corpora ay dapat malaki, diverse, at nacacapture ang mga general properties ng isang wika. Ang ginamit na training data ay ang WikiText-TL-39 \cite{CRUZCHENG:19} na mula sa mga artikulo sa Tagalog Wikipedia \footnote{https://tl.wikipedia.org/wiki/Unang\_Pahina}. Pinagisa ang ang train, valid, and test set at randomly na kinuha ang 10\% ng data bilang validation set at ang natirang 90\% ay ginamit bilang training set. Ang text data ay dumaan sa preprocessing\footnote{Para sa partikular na preprocessing rules, bisitahin ang https://docs.fast.ai/text.core\#Preprocessing-rules} bago gamitin sa training. Kinukuha lamang dito ang 60,000 na salita na pinaka madalas makita sa datos. Ang model ay sinanay ng 20 epochs na may learning rate na 1e-2, batch size na 128, at dropout multiplier na 0.5. Ang buong proseso ng training ay tumagal ng 26 hours.

\subsubsection{Language Model Fine-tuning}
Gamit ang pre-trained model mula sa unang phase, mas sasanayin pa ang model sa target text corpus upang maka adapt ang model sa wika at kung ano mang patterns at vocabulary ang meron ito. Ang target text corpus na ginamit sa papel na ito ay ang Hate Speech Dataset \cite{CABASAG}. Ang dataset ay hinati sa train, validation, at test set. Sa pag fine-tune ng model, sinanay muna yung last layer ng model for 1 epoch na may learning rate na 4e-2. Pagkatapos nito, ang lahat naman ng layers ng model ay sasanayin for 7 epochs na may learning rate na 4e-3.

\subsubsection{Text Classifier Fine-tuning}
Gamit ang fine-tuned LM mula sa pangalawang phase, nag append ng karagdagang layers para sa text classification task \cite{ULMFIT}. Ang model ay muling sinanay sa target text corpus na Hate Speech dataset pero hindi na sa language modelling task kundi sa text classification task na. Dito, kasama na ang labels sa text (0 = not hate, 1 = hate). Ang dropout multiplier ay 0.3, weight decay ay 0.1, at momentum ay (0.8,0.7,0.6). Ginamit ang fine-tuning techniques na gradual unfreezing at discriminative learning rates \cite{ULMFIT}. Makikita sa Table \ref{tbl:finetune-table} ang buong proseso ng fine-tuning.

\begin{table}
  \includegraphics[width=\linewidth]{finetune-table}
  \caption{Set ng hyperparameters na ginamit sa fine-tuning na may gradual unfreezing and discriminative learning rates. Ang lr ay learning rate at lr = 5e-2.}
  \label{tbl:finetune-table}
\end{table}

\subsection{Degradation Test}
Ang degradation test \cite{CRUZCHENG:20} ay isang paraan ng pagsukat ng resillience ng model sa performance degradation kapag binawasan ang training samples. Ang performance degradation ay nirereport bilang percentage drop ng metric ng isang task. Ito ay:

\[Degradation\% = \frac{Metric_{full}-Metric_{reduced}}{Metric_{full}}\]

na kung saan ang full ay ang performance sa kapag nag train sa buong training set at reduced ay ang performance kapag nag train sa reduced training set.

\begin{table*}[t]
\centering
\includegraphics[width=16cm]{degradation-table}
\caption{Degradation test results sa Hate Speech Dataset}
\label{tbl:degradation-table}
\end{table*}

\begin{table}
  \includegraphics[width=\linewidth]{weakpoints}
  \caption{Halimbawa ng mga top losses ng model. Ang target ang tamang label at ang prediction ay ang prediction ng model.}
  \label{tbl:weakpoints}
\end{table}

Ang model na may mabagal na degradation ay mas epektibo sa low-resource setting. Gamit ang Hate Speech dataset, hahatiin sa ang proseso sa tatlong setup. Sa unang setup, gagamitin ang buong training set o 10k samples. Sa pangalawang setup ay 50\% split o 5k training samples nalang at ang panghuli ay 10\% split o 1k training samples. Ang buong proseso ng training sa bawat setup ay kapareho lang ng nabanggit sa Section 3. Limang beses itong uulitin sa bawat setup at kukunin ang average test loss at accuracy nito sa bawat setup. Makikita sa Table \ref{tbl:degradation-table} ang resulta ng degradation test.

\section{Resulta at Diskusyon}
\subsection{Resulta ng Fine-tuning}
Matapos ang fine-tuning sa Hate Speech dataset, ang AWD-LSTM ay naka score ng 76.84\% accuracy sa test set ng Hate Speech dataset. Ito ay lamang ng 2.08\% sa best model sa baseline \cite{CRUZCHENG:20} ngunit marginal improvement lamang ito.

\subsection{Resulta ng Degradation Test}
Ang AWD-LSTM ay may nagkaron ng performance drop sa accuracy na 4.01\% sa 5k split na may degradation na 5.26\%. At sa 1k split naman, sobrang laki ng binaba nito sa accuracy na bumaba ng 8.6\% at may degradation na 11.28\%.

Mula sa mga resultang ito, makikita na mas epektibo lang ng kaunti ang AWD-LSTM sa Hate Speech classification task kumpara sa BERT models kapag ginamit ang 10k training samples. Ngunit ang AWD-LSTM ay mas mababa ang performance pagdating sa mga pagkakataon na maliit lang ang datos tulad ng sa 5k at 1k training samples. Ang initial baseline ng BERT \cite{CRUZCHENG:20} ay may degradation na 3.28\% on average sa 5k split. Ang worst model sa baseline o ang pinakamabilis ang degradation ay ang DistilBERT \cite{DISTILBERT} na may degradation na 4.34\%  sa 5k split. Samantala ang AWD-LSTM naman ay may 5.26\% na degradation sa 5k split. Makikita dito na mas epektibo ang transformer-based models kagaya ng BERT at DistilBERT pagdating sa mga low-resource tasks.

Hindi ito nakakagulat dahil ang transformer-based models ay nakadisenyo talaga para makakuha ng mas malalim na patterns mula sa data. Ang advantage lang ng AWD-LSTM ay mas mabilis itong i-train lalo na kapag from scratch gagawin ang pre-trained language model at kaya itong gawin sa isang GPU lang hindi kagaya ng transformer-based models na kailangan gumamit ng Tensor Processing Unit (TPU).

\subsection{Weak Points ng Model}
Makikita sa Table \ref{tbl:weakpoints} ang halimbawa ng mga kaso ng top losses o ang mga prediksyon ng model na may mataas na confidence level pero mali pala. Ang paraan ng pagsusulat na ginamit sa case 1 ay high-context kung saan ang kahulugan nito ay nakadepende sa context ng paguusap at hindi sa literal na kahulugan ng mga salitang ginamit. Isa pang maaring dahilan ay ang paggamit ng *facepalm* o mga expression na hindi ganoon kadalas gamitin. Dahil hindi ito madalas gamitin, ibig sabihin hindi ito madalas makita sa dataset. At dahil hindi ito madalas makita sa dataset, hindi matututunan ng lubusan ng model ang patterns at relationships ng salita o expression na ito.

Sa case 2, ang salitang "disappointing" ang may pinaka malaking impluwensiya sa prediction ng model. Ang hinala ay dahil bihirang salita lang ito sa dataset at madalas gamitin ang salitang ito sa negatibong konteksto. Napatunayan ang hinala dahil meron lang tatlo na occurences sa training set ang salitang "disappointing" at lahat ito ay naka label as "hate". Ito ang dahilan kung bakit mataas ang association ng salitang "disappointing" sa label na "hate" at naging malakas ang impluwensiya nito sa pag predict sa isang sentence as hate.

In general, mahina ang performance ng model kapag ang mga salitang ginamit ay hindi common. Para matugunan ang problema na ito, kinakailangan ng mas malaki at mas diverse na dataset. Dahil kapag mas malaki at mas representative ng real world ang dataset, mataas ang chance na mataas ang general performance ng isang model.

\section{Kongklusyon}
Kahit na mas malalim ang patterns na kayang makuha ng transformer-based models, kaya parin makipagsabayan ng AWD-LSTM sa simpleng text classification task ngunit ang performance nito ay naghihingalo kapag ang training data ay maliit. Naipakita sa papel na ito na Nananaig padin ang transformer-based models kagaya ng BERT pagdating sa low-resource tasks. Nirerekomenda namin na gumamit ng transformer-based models kapag merong available na pre-trained language model sa iyong target na wika. At kapag from scratch gagawin ang pre-trained language model at kung issue ang oras at high compute power, mas mainam na gumamit ng AWD-LSTM dahil higit na mas mabilis ang pre-training step nito kumpara sa transformer-based models. Kung nais mong gumawa ng aplikasyon na kailangan gamitan ng NLP, nirerekomenda namin na magsimula muna sa mas murang option na AWD-LSTM at tignan kung sapat na ang performance nito para magawa ang isang task. Kapag hindi pa sapat ang performance nito, saka lamang subukan ang transformer-based models.


\section{Introduction}
Ongoing efforts in the development of online measurement sensors and data acquisition \cite{Nardo2018} let new techniques emerge in real-time control of water distribution systems (WDSs).
The need for smart water networks and economic water supply is motivated mostly by the rapidly growing and urbanizing human population (\cite{urb2018} and \cite{pop2019}).
Ever-growing water networks challenge operators in achieving cost-effective, yet satisfactory operation.
According to the comprehensive literature review of \citet{Mala-Jetmarova2018}, operation optimization is the third most active topic among water supply-related applications.

Operation optimization is often treated as a pump scheduling problem, where dynamic programming (DP) methods are widely utilized because they are able to provide the global optimum if the number of variables is low.
The iterative DP method \cite{Zessler1989} is able to cope with the curse of dimensionality in large WDSs but its application is restricted to convex objective functions.
A more recent work by \citet{Bene2013} introduces a method based on discrete dynamic programming that finds the optimal pump scheduling for medium-sized water networks.
However, the multiple restrictions on the problem definition shows that modern variants of dynamic programming are computationally still too expensive.

\citet{Savic2018} found that many state-of-the-art (SOTA) optimization techniques in WDS operation optimization involves some kind of metaheuristics.
\citet{Ritter2020} presents a multiobjective direct policy search technique where different heuristic approaches share the individuals of a population among different metaheuristics.
This way the different beneficial properties are combined from different metaheuristics but incorrectly chosen hyperprameters of the distinct methods can alter the whole optimization process.

Moreover, these methods require a high amount of hydraulic evaluation which makes the use of surrogate models a reasonable approach.
Surrogate models are approximation models that are computationally cheaper than the corresponding -- physically valid -- numerical simulation.
The input of a surrogate model is usually the same as for the simulation (e.g. boundary and initial conditions) but its output is restricted to just a handful of variables.

Surrogate models with high complexity (e.g. Gaussian process regression, artificial neural networks, k-nearest neighbors) are found to be misleading for non-convex objective functions as these smooth the response surface that leads to suboptimal solutions \cite{Simpson2004}.
In concordance, \citet{Razavi2012} shows for water network simulations that substitutes with less complexity are more robust in cases where an unseen part of the explanatory variable space has to be discovered.
\citet{Machell2010} introduces a hybrid approach where, the governing equations of the hydraulic simulation are supported with an extended number of boundary conditions that reduces the computation cost of the hydraulic model.
The additional information is extracted from online sensory data including pressure and flow rate measures.

Optimization convergence can be enhanced by improving the optimum search itself to reduce the time of the overall process.
A promising branch of these efforts is about finding the optimal parameters of the optimization algorithm.
These parameters are called hyperparameters, and methods that optimize these hyperparameters during the single optimization process are referred to as hyperheuristic methods \cite{Savic2018}.
\citet{MccLymont2013} present a Markov chain-based approximator that is superior in finding optimal hyperheuristics compared to Simple Random and TSRoulWheel methods.

Though heuristics are robust in finding the global optima, the computational burden is still too high for real-life WDSs.
Model predictive control (\cite{Pascual2011} and \cite{Wang2017}) overcomes the computational issue by using a two-layer optimization approach where the layers predict optimal ON-OFF schedules for pumps and valves in different time horizons.
The method is able to run large WDSs in real-time, the capabilities are presented on the example of the water network of Barcelona.
Iterative extended lexicographic programming (iELGP), introduced by \citet{Abdallah2017} and \citet{Abdallah2019} copes even with continuous decision variables like the set-point of variable speed pumps (VSPs) and chlorine dosing machines.
With iELGP, \citet{Abdallah2019} achieved SOTA performance on the C-Town benchmark problem wherein both fixed-speed and variable speed pumps were presented.

Water distribution system operation also benefits from recent advances in machine learning.
Artificial neural networks are used in a wide variety of applications from demand forecast \cite{Guo2018} through attack detection \cite{Taormina2018} to pump scheduling \cite{Bhattacharya2003}.

\citet{Bhattacharya2003} uses an ANN-based reinforcement learning agent as an optimal controller for switching on or off pumps in a water network.
The input of the ANN is water levels in reservoirs besides the water demand in the WDS, while its output is the decision about a pump to operate it or not.
The neural network is trained with reinforcement learning (RL) in a manner the RL agent is encouraged to replicate the moves of the master controller algorithm.

\citet{Mahootchi2007} introduces an agent, based on reinforcement learning for the control of one reservoir in a water distribution network and finds that the RL approach achieves comparable results to stochastic dynamic programming in the benchmark case.

A sophisticated optimization technique is introduced by \citet{Castelletti2013} for multiobjective tasks, where the importance of the standalone objectives is unknown or uncertain at the beginning of the optimization.
Recombining these objectives yields a different optimal solution, hence a movement on the Pareto frontier.
\citet{Castelletti2013} approximate the optimal policy with fitted q-iteration (FQI) during an optimization run but also add the weights of the standalone objectives to the state space of the environment used in FQI.
The approximation of the optimal policy is constructed as the function of the importance weights, and the frontier can be drawn from one single run with sampling the weight space as well.

All the pump scheduling methods above rely on the simulation of the water distribution system that limits their usability when the size or the complexity of the water network is too large.
The most resource-intensive method is finding optimal pump speeds for a time step since the descriptive variables are continuous and the evaluation of every trial needs to run the hydraulic simulation as well.

A reinforcement learning agent for real-time control for variable speed pumps is presented in this paper.
The pump scheduling problem is omitted and only the optimum speed search for a given demand-distribution is realized in a smart water network.
A significant novelty of the method is that the hydraulic simulation of the WDS is not needed for the controller because it relies only on measurement data.
Moreover, the training is carried out in such a way that the agent is able to develop different policies compared to the reference, so that the agent can outperform the master algorithm it was trained with.
If online pressure data are not available, the hydraulic simulation of the water distribution system can be used as a surrogate of the real WDS, like in the case of existing optimization methods.
In such cases, the agent-based control is able to outperform current optimization techniques in the number of hydraulic model evaluation.
As the proposed technique is able to cope with multiple objectives with a simple weighted sum of objectives, it can serve as the optimal pump speed search algorithm in hybrid pump schedule optimization solutions.

The paper is organized as follows. First, the simplified problem is introduced including the reference methods. The details of reinforcement learning methodology are described thereafter. Then, the settings of the case study, including the water distribution networks and the training settings are presented. The performance of the agent is evaluated on the WDSs and the conclusions of the training are drawn. Finally, the key findings and possible applications are mentioned.

\section{Methodology for pump speed control with reinforcement learning}
\subsection{Problem definition and reference solutions}
The present study aims at providing optimal set-points for variable speed pumps in a WDS according to given nodal pressures.
Hence, the conventional pump scheduling problem is simplified and the problem is formulated as follows.

Instantaneous operation point optimization is carried out neglecting temporal changes in the demand.
The competing algorithms have to find the optimal pump speeds in terms of hydraulic efficiency with respect to randomly generated demands.
Besides, the inflow and outflow into and from the tanks should be minimized, and pressure (head) should be kept in a specific range in all of the nodes.
Customer satisfaction implies a lower pressure limit, while the mechanical strength of water fittings and other elements of the water network implies an upper limit as well.

These objectives are formulated into a simple weighted sum of objectives.
The meaning of this combined value is two-fold, as it serves as an objective value for the wide-spread optimization algorithms and is interpreted as the value of the state in the Markovian decision process described later.
The rating of the final state is the same regardless, of whether it is achieved by the reference optimization methods or the reinforcement learning-based agent.

The calculation depends on the observation of the environment that is the ensemble of the state variables and some other quantities that are not involved in the decision-making but are measurable.
The algorithm of the calculation is presented in Alg.~\ref{alg:state-value}.
The input of the algorithm are the state variables $\matr{S}$, the set of the pump efficiencies $\matr{E}_\mathrm{pump}$, the volumetric flow rate $\matr{F}_\mathrm{pump}$ delivered by the pumps and the signed volumetric flow rate $\matr{F}_\mathrm{tank}$ of the tanks.

\begin{algorithm}
\caption{Calculating the state value $v(\matr{S})$}
\label{alg:state-value}
\begin{algorithmic}[1]
    \Function{CompStateValue}{$S, E_\mathrm{pump}, F_\mathrm{pump}, F_\mathrm{tank}$}
        \State $\matr{H}=\{h_{1},\dots,h_{N}\} \gets \Call{ExtractNodalHeads}{\matr{S}}$
        \State $n_\mathrm{wrong} \gets 0$
        \ForAll {$h \in \matr{H}$}
            \If{$h < h_\mathrm{min} \OR h > h_\mathrm{max}$}
                \State $n_\mathrm{wrong} \gets n_\mathrm{wrong}+1$
            \EndIf
        \EndFor
        \State $c_\mathrm{tot} \gets \sum{\matr{F}_\mathrm{pump}}+\sum{\matr{F}_\mathrm{tank}}$
        \State $\eta_\mathrm{tot} \gets \prod{\matr{E}_\mathrm{pump}}$
        \State $r_\mathrm{feed} \gets \frac{c_\mathrm{tot}}{c_\mathrm{tot}+\sum_{i=1}^{{\#}\matr{F}_\mathrm{tank}}{\vert{f_i}\vert}}$
        \State $r_\mathrm{satisfaction} \gets 1-\frac{n_\mathrm{wrong}}{\#\matr{H}}$
        \State $r_\mathrm{eff} \gets \frac{\eta_\mathrm{tot}}{\eta_\mathrm{limit}}$
        \State $v \gets w_\mathrm{satisfaction} \cdot r_\mathrm{satisfaction} + w_\mathrm{eff} \cdot r_\mathrm{eff} + w_\mathrm{feed} \cdot r_\mathrm{feed}$
        \State \Return $v$
    \EndFunction
\end{algorithmic}
\end{algorithm}

The state value $v$ is a weighted combination of $3$ objectives.
\begin{enumerate}
     \item $r_\mathrm{satisfaction}$ represents the satisfaction of the consumers that is characterized by the number of the insufficiently supplied nodes $n_\mathrm{wrong}$ wherein the head is out of the acceptable range bounded by $h_\mathrm{min}$ and $h_\mathrm{max}$. The number of problematic nodes is divided by the number of all nodes in the water distribution system, and the ratio is subtracted from $1$. The outcome is the basis of the first weighted objective, the satisfaction ratio.
     \item $r_\mathrm{eff}$ represents the efficiency of the pumps, and it is just the product of the standalone pumps feeding the network divided by the product of the theoretical peak efficiencies of all of the pumps.
     \item $r_\mathrm{feed}$ represents the feed ratio of the water network that is solved by the pumps supplying water from the wells or other sources, on the one hand and by the tanks and reservoirs, on the other hand.
\end{enumerate}
In the current experiment, using the water reserve is considered expensive, and using the reserve is punished by a lower objective value, which is a simplification of the problem.
However, this methodology is easily changeable by adding the target inflows and outflows to the state value and comparing the actual ones to those.
In this way, the agent would be forced to learn to maintain the water flux of the tanks but this is a harder problem to solve.
For the time being, this is ignored by forcing the agent to minimize the water flux of the tanks.
Thus, the feed ratio is the amount of total consumption divided by the sum of the total consumption and the water flux of the tanks.

The standalone objectives are also presented in Eq.~\ref{eq:objectives}, where the notations are the same as for Alg.~\ref{alg:state-value} except for $f_{\mathrm{tank}}$ that denotes the sum of the absolute value of inflow/outflow for every tank.

\begin{equation}\label{eq:objectives}
    \begin{split}
        r_{\mathrm{satisfaction}} &= 1 - \frac{n_{\mathrm{wrong}}}{n_{\mathrm{nodes}}},\\
        r_{\mathrm{eff}} &= \frac{\eta_{\mathrm{tot}}}{\eta_{\mathrm{limit}}},\\
        r_{\mathrm{feed}} &= \frac{c_{\mathrm{tot}}}{c_{\mathrm{tot}}+f_{\mathrm{tot, tank}}}.
    \end{split}
\end{equation}

The above formulation of the standalone objectives leads to $3$ measures, whose maximum values are $1.0$ all.
These are then weighted separately by constant factors.
The coefficients in the presented experiments are chosen arbitrarily $\nicefrac{8}{16}$, $\nicefrac{5}{16}$ and $\nicefrac{3}{16}$ as the ratio of the satisfaction, the efficiency and the feed, respectively.

The task was solved with different optimization techniques to provide reference: fixed step-size random search (FSSRS), particle swarm optimization (PSO), differential evolution (DE), and Nelder-Mead method.
PSO, DE and Nelder-Mead method are able to solve operation point optimization of VSPs, see the work of \citet{Baltar2008}, \citet{Odan2015} and \citet{Celi2017}, respectively.

A suboptimal optimization method was also applied on the task that is referred to as one-shot random trial (One-shot RT).
This method sets pump speeds randomly inside the feasibility domain and it serves as a baseline: the reinforcement learning agent has to be able to outperform this approach if it can interpret any information on the connection of nodal pressures and optimal pump speed settings.

\subsection{Dueling deep Q-network agent as a real-time controller for water networks}
The drawback of the reference methods is the incapability to accomplish real-time control because of the necessity of computationally expensive simulation of the hydraulics.
In this paper, a dueling deep q-network agent is introduced to optimize the speeds of the pumps in a water distribution system.
The agent relies on simulation data in the training phase only; thus, the controlling is feasible real-time, regardless of the size of the water network and the number of pumps controlled.

The proposed technique is based on Markov decision processes (MDPs) that are key building blocks of reinforcement learning.
In the terminology of MDPs an \textit{agent} works in an \textit{environment} by taking \textit{actions} according to a  \textit{policy}.
The environment is observable to the agent through its \textit{state} that changes according to the actions.
The instantaneous quality of a state is measured by the \textit{reward} that is signaled to the agent.
Rewards of a decision sequence are combined in the \textit{value function}.
The value function measures the overall quality of the starting state in the decision sequence according to the given policy.
The goal of the learning is to find a policy with which the agent can govern the environment to the most advantageous state from any initial state.

In the present study, the agent is a dueling deep q-network and it interacts with the environment that is based on the water network.
The decision-making sequence and most of its concepts are depicted in Fig.~\ref{fig:rl-wds} (originating from \citet{Sutton2018}) where $t$ and $t+1$ refer to two consecutive steps (also referred to as time steps) in the sequence.

\begin{figure}
    \centering
    \def3cm{3cm}
    \def1cm{1cm}
    \def.3cm{.3cm}
    \resizebox{\linewidth}{!}{
    \begin{tikzpicture}
        \node (agent) [agent] {
            {Dueling deep q-network.}
            };
        \node (env) [env, yshift=-{3cm}] {
            {Manipulation and observation of the WDS, reward calculation of the state}
            };
        \draw [thick, dotted] ($(env.west)+(-1cm,-1cm)$) -- ($(env.west)+(-1cm,1cm)$);
        \draw [arrow] (agent.east)
            -- ($(env.east)+(1.5cm,3cm)$)
            -- node[anchor=south, rotate=90, align=center, text width=3cm] {tuning the speed of a pump} ++(0,-{3cm})
            -- (env.east);
        \draw [arrow] ($(env.west)+(0,+.3cm)$)
            -- node[anchor=south, align=center] {$t\!+\!1$} ++($(-1cm,0)$);
        \draw [arrow] ($(env.west)+(0,-.3cm)$)
            -- ++($(-1cm,0)$);
        \draw [arrow] ($(env.west)+(-1cm,+.3cm)$)
            -- node[anchor=south, align=center] {$t$} ++(-.751cm,0)
            -- node[anchor=south, rotate=90, align=center] {reward} ++($(0,3cm-2*.3cm)$)
            -- ($(agent.west)+(0,-.3cm)$);
        \draw [arrow] ($(env.west)+(-1cm,-.3cm)$)
            -- ++(-1.51cm,0)
            -- node[anchor=south, rotate=90, text width=3cm, align=center] {nodal pressures \& pump speeds} ++($(0,3cm+2*.3cm)$)
            -- ($(agent.west)+(0,+.3cm)$);
    \end{tikzpicture}
    }
    \caption{The agent-environment ecosystem tailored to the controlling of the pumps in a  water distribution system.}
    \label{fig:rl-wds}
\end{figure}

According to Fig.~\ref{fig:rl-wds}, the communication between the environment and the agent during the training is as follows.
\begin{enumerate}
    \item The environment provides its state variables -- the ensemble of the nodal pressures and the pump speed ratios -- to the agent in a concatenated vector. Speed ratios are the actual pump speed divided by the nominal pump speed.
    \item The agent maps the state values to an action according to its current policy and executes it in the environment. The action can be the increase or decrease of one of the pump speed ratios or keeping the speed ratios as they were set in the previous step. The latter is referred to as \textit{siesta}.
    \item The new hydraulic properties are provided by the hydraulic simulator, and the value of the new state is calculated. Besides, the reward is calculated from the state value and some other properties of the environment.
    \item The environment propagates the new state and the corresponding reward to the agent to induce the next move.
    \item At last, (1) the preceding state, (2) the preceding action taken, (3) the corresponding reward, and (4) the new state is stored as a training sample wherewith the action-value estimator is updated later.
\end{enumerate}

This cyclic process goes on until the environment reaches a terminal state, which can be caused by two distinct events.
On the one hand, the agent is allowed to run until a fixed number of steps to get to the most advantageous state.
The environment keeps track of the number of steps taken, and it terminates the episode when the agent reaches the limit.
On the other hand, the agent can decide to do nothing instead of changing one of the speed ratios.
If the agent adheres to a state for a predefined number of consecutive steps, the episode is terminated as well.

Seeking to achieve the optimal policy leads to the grading of different policies.
A policy is regarded better than another one if the value function is equal or higher for all states by taking actions under that one.
The value function can be formulated to depend on different variables, the related form is the action-value function, that is for policy $\pi$:
\begin{equation}\label{eq:action-value_fun}
    q_{\pi}(s, a) \doteq \mathbb{E}_{\pi} \left[ \sum_{k=t+1}^{T}{\gamma^{k-t-1} R_{k}} \biggr\rvert S_{t}=s, A_{t}=a \right] \;.
\end{equation}

The value is represented as an expected value of the sum of rewards $R_k$ in the forthcoming sequence started from state $S_t$ and taking action $A_t$ under policy $\pi$.
The importance of future rewards can be attenuated with the discount-factor $\gamma$, such that $\gamma \in (0, 1]$.

If the agent runs under the optimal policy, its actions maximize the reward in every state-transition; thus, the optimum of the action-value function is:
\begin{equation}\label{eq:max_action-value}
    q_{\pi^{*}}(s, a) \doteq \mathbb{E}_{\pi^{*}} \left[ R+\gamma\max_{a'}q^{*}(s', a') \biggr\rvert S_{t}=s, A_{t}=a \right] \;,
\end{equation}
where the forthcoming action and state is denoted by a prime and $q^{*}$ refers to action-values under the optimal policy.

Eq.~\ref{eq:max_action-value} can be solved in different ways, for the relevant techniques see \citet{Sutton2018}.
SOTA solutions use artificial neural networks as value function estimators to cope with the complexity of the equation.
ANNs are good candidates for the task, because they are universal function approximators as shown by \citet{Cybenko1989}.

A variant of deep q-networks (DQN, introduced by \citet{Mnih2013}) is chosen in the present study to approximate the action-value function.
The term deep refers here to the fact that the neural network architecture has two or more hidden layers.
The agent can be trained by minimizing the loss
\begin{equation}\label{eq:dqn}
    L_{i}(\theta_{i}) = \mathbb{E}_{\pi} \left[ (r+\gamma\max_{a'}q(s',a',\theta_{i-1}) - q(s,a,\theta_{i}))^2 \right]\;,
\end{equation}
where $\theta_{i}$ corresponds to the weights in the ${i}^{th}$ update of the deep q-network (a deep neural network in this case) and $L_{i}$ is the mean squared error defined between the $(i-1)^{th}$ and $i^{th}$ approximation of the value function.
A significant advantage of Eq.~\ref{eq:dqn} over Eq.~\ref{eq:max_action-value} is that it does not have to evaluate a whole decision sequence to get a new assumption on the value of the updated state; hence, it can be trained off-policy.

Eq.~\ref{eq:dqn} is optimized with stochastic gradient descent method as it is a general, yet effective optimization algorithm for training ANNs \cite{bottou-bousquet-2008}.
Moreover, the last layer of the deep q-network is forked: one layer approximates the value of the current state, while the other layer estimates the advantage in the value by taking the action.
This is intended to speed up training convergence and the technique is called dueling deep q-networks, as introduced by \citet{Wang2015}.

\begin{algorithm}[ht]
\caption{Step procedure of the agent}
\label{alg:step}
\begin{algorithmic}[1]
    \Procedure{Step}{$a$}
        \State $n_\mathrm{siesta} \gets 0$
        \State $n_\mathrm{steps} \gets n_\mathrm{steps}+1$
        \State $terminal \gets n_\mathrm{steps} = N$
        \State{$r \gets penalty$}
        \If{$a\not=siesta$}
            \State{$\matr{PS} \gets \Call{UpdatePumpSpeeds}{a}$}
            \State{$\delta \gets \lVert{\matr{PS}-\matr{PS}_\mathrm{ref}}\rVert$}
            \If{$any(\matr{PS}<ps_\mathrm{lo}) \OR any(\matr{PS}>ps_\mathrm{hi})$}
                \State{$\matr{PS} \gets \Call{RevertAction}{a}$}
            \Else
                \If{$\delta<\delta^*$}
                    \State{$r \gets \delta \cdot constant$}
                \EndIf
                \State{$\matr{S} \gets \Call{UpdateState}{a}$}
                \State{$\matr{E}_\mathrm{pump}, \matr{F}_\mathrm{pump}, \matr{F}_\mathrm{tank} \gets \Call{Observe}{\matr{S}}$}
                \State{$v \gets \Call{CompStateValue}{\matr{S}, \matr{E}_\mathrm{pump}, \matr{F}_\mathrm{pump}, \matr{F}_\mathrm{tank}}$}
                \State{$\delta^* \gets \delta$}
            \EndIf
        \Else
            \State{$n_\mathrm{siesta} \gets n_\mathrm{siesta}+1$}
            \If{$n_\mathrm{siesta}<3$}
                \If{$1-{v}/{v_\mathrm{ref}}<\epsilon$}
                    \State{$r \gets n_\mathrm{siesta} \cdot constant$}
                \EndIf
            \Else
                \If{$1-{v}/{v_\mathrm{ref}}<\epsilon$}
                    \State{$r \gets bonus$}
                \EndIf
                \State{$terminal \gets \texttt{True}$}
            \EndIf
        \EndIf
        \State \Return $r, \matr{S}, terminal$
    \EndProcedure
\end{algorithmic}
\end{algorithm}

A crucial component of the reinforcement learning setting is the rewarding of the agent as it influences the behavior of the trained agent.
Alg.~\ref{alg:step} summarizes the step procedure with the rewarding system of the present study.

The reward is based upon the state value but it is not a direct match, as other goals have to be reached beside optimality, as shown in Alg.~\ref{alg:state-value}.
As the control is realized stepwise, the agent has to be taught to set the optimal speeds in as few steps as possible.
The basic idea behind the presented step function, including the reward calculation is as follows.
At the beginning of every new scenario, i.e., after the generation of a new demand-distribution, the optimal setting of the pump speeds is obtained by the Nelder-Mead method.
This speed setting $\matr{PS}_\mathrm{ref}$ and the corresponding state-value $v_\mathrm{ref}$ is considered as a reference, guiding the agent softly in the search space.
After the agent intervenes, the reward is based upon the closeness $\delta$ of the reference solution and the agent's solution in the search space, defined by the feasible pump speed ratios.
At this stage, the agent gets a reward $r$ if it gets closer to the reference solution on a monotonous trajectory, i.e., it does not make steps increasing the distance $\delta$ between the solutions.
The agent can also choose to make a \textit{siesta}; in this case, the reward is computed differently.
In siestas, the actual value of the state is computed according to Alg.~\ref{alg:state-value} and compared to the reference value.
If the current state value is higher than the $98\%$ percent of the reference value (i.e., the value of $1-v/v_\mathrm{ref}$ is at most $\epsilon$), the agent is rewarded with a slightly higher score.
The immediate reward achieved in the first siesta is doubled if the agent chooses to stay in this beneficial state for the second time.
Finally, the episode is terminated after the third consecutive siesta, and the agent gets an extra reward that is significantly higher than the rewards that are achievable with the active moves.
In contrast, the agent gets a negative reward whenever it recedes from the reference speed setting with an active action or if it stays in a suboptimal position, i.e., if it makes a siesta in a state, with a value below $98\%$ of the reference.
A penalty is given even when the agent tries to set pump speeds lower or higher as the feasibility limits $ps_\mathrm{lo}$ and $ps_\mathrm{hi}$, respectively.
With this concept, the agent is incited to move along the optimal position continuously and terminate the episode there.
The $2\%$ tolerance in the comparison of the state values is needed because of the different variable space from where the solution is acquired with the reference method and with the reinforcement learning agent: the Nelder-Mead method chooses candidate solutions from a continuous variable space while the dueling deep q-network handles a discrete action space; thus only a restricted set of pump speeds is available compared to Nelder-Mead.
Although the discretization can be arbitrarily fine, the dimensionality of the q-function -- hence, the length of the training of its estimator -- depends on the size of the discrete parameter space.

The soft guidance means in the above process that the agent can achieve high rewards even far from the reference solution if it has a similar value.
This is encouraged during training by taking random actions in a part of the whole process, thanks to the $\varepsilon$-greedy policy.

After training, the agent relies only on the observation of the environment to take its next action; thus, the calculation of the reward -- including the computationally expensive reference calculation -- is completely omitted.

\section{Case study}
\subsection{Reference water distribution systems}
The proposed controller is trained and tested on the modified Anytown and D-Town (referred to as C-Town in some literature) water distribution systems that are common among the WDS design and management research community.
These networks serve as benchmarks for water network design methods wherewith the networks can be modified and extended to comply with increased water needs and expanding consumption areas.
The original topology of both WDSs is accessible online from the web page of the Centre of Water Systems research group at the \citet{cwsExeter}.
Moreover, the numerical models of the WDSs are also available in an \verb+EPANET+-compatible format.
The modified versions are referred to as Anytown-mod and D-Town-mod in the following.

\begin{figure}[ht]
    \centering
    \includegraphics[width=\linewidth,keepaspectratio]{./figures/anytown.pdf}
    \caption{Topology of the Anytown water distribution system. Network elements plotted: nodes, pipes, pumps and tanks.}
    \label{fig:anytown_topo}
\end{figure}
\begin{figure}
    \centering
    \includegraphics[width=.9\linewidth, keepaspectratio]{./figures/dtown.pdf}
    \caption{Topology of the D-Town water distribution system. Network elements plotted: nodes, pipes, pumps and tanks.}
    \label{fig:dtown_topo}
\end{figure}

The topology of the networks is depicted in Fig.~\ref{fig:anytown_topo} and in Fig.~\ref{fig:dtown_topo} for the Anytown and D-Town test cases, that are the same for Anytown-mod and D-Town-mod, respectively.
The networks are slightly modified as described below.
The primary purpose of Anytown-mod is to provide a simple, easily solvable problem that serves as a proof of concept: it has a global optimum that can be found with conventional methods easily.
If the agent performs poorly on Anytown-mod, then there is no chance that it can solve more complex tasks.
D-Town-mod is a more challenging problem whose solution is not trivial because of the multi-dimensional parameter space of the optimization loss function.
If the agent performs well on D-Town-mod, then the RL-based agent with the proposed reward function is a good candidate to serve as a real-time controller in real water networks.
The main properties of the two water distribution systems are listed in Table~\ref{table:benchmark_prop}.
Both water networks share their topology with its original version published by the \citet{cwsExeter}; the differences are summarized below.

\begin{table}
    \caption{Main properties of the water distribution systems}
    \label{table:benchmark_prop}
    \centering
    \small
    \renewcommand{\arraystretch}{1.25}
    \begin{tabular}{l l l}
    \hline\hline
    \multicolumn{1}{c}{Attribute of the water network} &
    \multicolumn{1}{c}{Anytown-mod} &
    \multicolumn{1}{c}{D-Town-mod} \\
    \hline
    Number of nodes & $22$ & $399$ \\
    Number of pipes & $41$ & $443$ \\
    Number of pump stations & $1$ & $5$\\
    \hline\hline
    \end{tabular}
    \normalsize
\end{table}

Anytown-mod is almost identical to Anytown as only the number of pumps in the single pumping station is reduced to two from three, and both of the pumps operate with the same continuously variable speed.
The performance curves of the pumps are the same as in the original network.

In the case of D-Town-mod, the control logic is removed, and all pumps operate in the water distribution network.
The performance curves are kept for the head, and fictitious curves are added for the efficiency.

An environment is created that serves as an interface between the agent and the water distribution system.
The purposes of the environment are the interpretation of the commands of the agent, reward calculation and re-initialization of the water network.
The state of the water network can be initialized by randomizing the demand map and by randomizing the pump speeds or resetting them to the original value.
The original demands in the task definition of the water network are stored, and the total sum of the original nodal consumption is calculated and stored at the beginning of training.

Demand variations are substituted with randomized demand distributions over the water distribution system.
This can be done because of training the agent with an off-policy method.
It is not necessary to show real-life demand variations to the agent in the training phase as it is trained off-policy.

The randomization of the nodal consumption is done in three steps to be able to mimic the total demand variation during day and night and the nodal demand variation as well.
First, new demand is generated in every node as the product of the original demand and a random value from a truncated normal distribution.
Second, new total demand is generated as the product of the original total demand and a random value from a uniform distribution.
Third, the nodal consumptions are rescaled such that the sum is equal to the one generated in the second step.
The process is summarized in Alg.~\ref{alg:consumption_init} where the set of original nodal consumptions and the randomized nodal consumptions are denoted by $\matr{C}_0$ and $\matr{C}$, respectively.

\begin{algorithm}
\caption{Initializing nodal consumptions}
\label{alg:consumption_init}
\begin{algorithmic}[1]
    \Function{InitNodalConsumption}{$C_{0}$}
        \State{$\matr{C} \gets \matr{C}_{0}$}
        \State{$c_{\mathrm{tot}} \gets \Call{RndUniform}{{\scriptscriptstyle\texttt{lo}=0.3, \texttt{hi}=1.1}} \cdot \sum{\matr{C}_{0}}$}
        \State{$c_{tot^*} \gets 0$}
        \ForAll{$c \in \matr{C}$}
            \State{$c \gets \Call{RndTruncNormal}{{\scriptscriptstyle\mathrm{mean}=1.0, \mathrm{stddev}=1.0, \texttt{lo}=0.7, \texttt{hi}=1.3}} \cdot c$}
            \State{$c_{tot^*} \gets c_{tot^*} + c$}
        \EndFor
        \ForAll{$c \in \matr{C}$}
            \State{$c \gets c \cdot \frac{c_{tot}}{c_{tot^*}}$}
        \EndFor
        \State{\Return{$\matr{C}$}}
    \EndFunction
\end{algorithmic}
\end{algorithm}

This randomization method allows control on the demand map both in relative and in absolute terms.
The nodal consumption randomized one-by-one mimics the relative variation of the demand map.
Meanwhile, the overall change in consumption through the day-night cycle is imitated by the variation of the total demand.
I.e., the total demand of the area covered by the water distribution network is considered first; then, the total value is partitioned among the nodes.

Both distributions are chosen arbitrarily and their probability density functions are depicted in Fig.~\ref{fig:proba_dists}.
The lower and upper bounds of the uniform distribution are $0.3$ and $1.1$, respectively, while the normal distribution has a mean $1.0$ and standard deviation $1.0$ with tails cut off below $0.7$ and above $1.3$.
The shape of the truncated normal distribution is similar to the uniform distribution with this parameter set, although multipliers close to one have slightly higher probabilities during the demand randomization than remote values.

\begin{figure}
    \centering
    \includegraphics[width=.8\linewidth, keepaspectratio]{./figures/proba_densities.pdf}
    \caption{Probability distributions for the three-stage demand randomization.}
    \label{fig:proba_dists}
\end{figure}

The pumps feeding the water network are considered as variable-speed pumps driven by e.g., frequency converters, in case of both water networks.
The absolute pump speed is not considered; the ratio compared to the nominal pump speed is recorded instead.
In the current environment, the speeds can be lowered to half of the maximum speed as water transportation is infeasible at lower speeds due to insufficient pump heads relative to the geodetic height of the nodes.
The cause of the lower speed limit is valid for both WDSs.

Finally, the well from where the pumps feed the network is considered as a reservoir with a water surface large enough to ensure a stationary water level during the speed-adjustment of the pumps.

The presented results are obtained with the default hyperparameters of optimization algorithms.
In Fig.~\ref{fig:opti_results}, the average performance of the different algorithms is presented on the test set in the case of Anytown-mod and D-Town-mod.
Moreover, the box plot of the difference calculated for every test scenario is depicted in Fig.~\ref{fig:opti_diffs} where the results achieved by the differential evolution are chosen as a reference.

\begin{figure}[ht]
    \centering
    \includegraphics[width=\linewidth,keepaspectratio]{./figures/opti_values.pdf}
    \caption{Results of the optimization on the test set with different techniques. The higher state values are the better.}
    \label{fig:opti_results}
\end{figure}

\begin{figure}
    \centering
    \includegraphics[width=\linewidth,keepaspectratio]{./figures/rel_opti_values.pdf}
    \caption{Box plots and swarm plots of the performance of the optimization algorithms compared to the differential evolution method. $1.00$ means the same performance that is achieved with DE.}
    \label{fig:opti_diffs}
\end{figure}

The similarity of the results suggests that each technique is able to find the empirical global optima, except for the one-shot random trial, which is reasonable.

The single pump station in Anytown-mod means that the pump speed optimization problem has a one-dimensional parameter space; hence, evolutionary algorithms and related techniques are inefficient in finding the solution as the crossover-like operations are meaningless on a one-dimensional candidate solutions.
Thus, only the performance achieved by the Nelder-Mead method was considered in the evaluation of scenarios from the Anytown-mod case.

\subsection{Training procedure of the agent}
\begin{table*}
    \caption{Possible hyperparameters for training}
    \label{table:ho}
    \centering
    \small
    \renewcommand{\arraystretch}{1.25}
    \begin{tabularx}{\textwidth}{>{\raggedright\arraybackslash\hsize=.25\hsize\linewidth=\hsize}X>{\centering\arraybackslash\hsize=.375\hsize\linewidth=\hsize}X>{\centering\arraybackslash\hsize=.375\hsize\linewidth=\hsize}X}
    \hline\hline
    \multicolumn{1}{c}{Hyperparameter} &
    \multicolumn{1}{c}{Anytown-mod} &
    \multicolumn{1}{c}{D-Town-mod} \\
    \hline
    Learning rate & \multicolumn{2}{c}{$\{10^{-4}, 10^{-3}, 10^{-2}\}$} \\
    Discount-factor & \multicolumn{2}{c}{$\{0.9, 0.99, 0.999, 0.9999\}$}  \\
    Batch size & $\{8, 16, 32, 64\}$ & $ \{16, 32, 64, 128\}$ \\
    \makecell[l]{\# of neurons \\ in the hidden layers} & \makecell{$\{[24]$, $[24, 6]$, \\ $[48,12]$, $[24, 12, 6]$, \\ $[48, 32, 12]\}$} & \makecell{$\{[300]$, $[300, 10]$, \\ $[128, 12]$, $[256, 32, 12]$, \\ $[256, 128, 12]\}$} \\
    \hline\hline
    \end{tabularx}
    \normalsize
\end{table*}

Two separate agents are trained during the experiment: one to solve the Anytown-mod problem and one to solve the D-Town-mod problem, respectively.

The agents were trained on the water networks separately, thus, during a training session, only one of the water distribution systems is utilized in the environment, and the trained agent will perform well only on the specific water network that was used for training.
The agents trained with different hyperparameters; the common part is just the underlying algorithm and the fact that a densely connected neural network is used as an action-value estimator.

The training will be described next.
At the beginning of every episode, the demands are randomized according to Alg.~\ref{alg:consumption_init}, the pump speed ratios are randomized according to a uniform distribution between the feasibility limits.
After every $\nicefrac{1}{25}th$ of the total number of training steps the training is suspended, and the performance of the agent is tested on $100$ different scenarios.
These scenarios are generated beforehand the training and they serve as validation data.
Additional $50$ test scenarios are generated, these scenarios are solved by the agent only after the training is completed.

A hyperparameter optimization was carried out prior to the training with a tree-structured parzen estimator (TPE) \cite{NIPS2011_4443}.
The critical hyperparameters were identified according to the previous experience of the authors.
Table~\ref{table:ho} summarizes possible hyperparameter values from where the optimal values were selected by the TPE.
Additional information on the sensitivity on the parameters and on the hyperparameter optimization is not presented due to length considerations.

The input layers of the neural networks have exactly as many neurons, as nodes in the water network as the nodal pressures divided by the highest shut-off head of the pumps are handed over the action-value estimator.
The activation functions are rectified linear units (ReLU) that are introduced in \citet{Nair2010}.
ReLU activation can handle input data outside the normalized range when the data are positive-valued, as is the case with nodal pressures under normal circumstances.
Nodal pressures are handed over to the neural network in vectors in every scenario, i.e., the $k^\mathrm{th}$ element of the pressure vector belongs to the same node under all circumstances.
The training is carried out with the $\varepsilon$-greedy policy starting with an exploration factor of $0.95$ and linearly decreased to $0.0$ at the final episode.
$\varepsilon$-greedy policy means that the agent makes random moves with a probability of the exploration factor, i.e., acting off-policy.
This behavior ensures the possibility to find better policies than the reference one, i.e., the agent is able to achieve high rewards even without mimicking the reference strategy.
The discount-factor affects the anticipatory behavior of the agent, the higher the value, the more prescient the agent is.
A scenario is $40$ and $200$ timesteps long in the case of Anytown-mod and D-Town-mod, respectively.
Moreover, the agents make moves respectively $50000$ and $1000000$ times in the smaller and the larger water distribution network before the final test.
The state-action-reward-new state data are saved to a so-called replay memory from where the training algorithm takes samples to fill the training batch for every update of the action-value estimator.
The hyperparameters of both agents are summarized in Table~\ref{table:hyperparams}.

\begin{table}
    \caption{Hyperparameters of the agents}
    \label{table:hyperparams}
    \centering
    \small
    \renewcommand{\arraystretch}{1.25}
    \begin{tabular}{l c c}
    \hline\hline
    \multicolumn{1}{c}{Hyperparameter} &
    \multicolumn{1}{c}{Anytown-mod} &
    \multicolumn{1}{c}{D-Town-mod} \\
    \hline
    Learning rate & \multicolumn{2}{c}{$0.0001$} \\
    Discount-factor & $0.99$ & $0.9$\\
    Batch size & $8$ & $64$ \\
    \# of neurons in hidden layers & $[48, 32, 12]$ & $[256, 128, 12]$ \\
    Size of replay memory & $25000$ & $350000$ \\
    Type of activation function & \multicolumn{2}{c}{rectified linear unit} \\
    Initial exploration factor & \multicolumn{2}{c}{$0.95$} \\
    Final exploration factor & \multicolumn{2}{c}{$0.00$} \\
    Maximum episode length & $40$ & $200$ \\
    Speed ratio increment & \multicolumn{2}{c}{$0.05$} \\
    \# of initialization steps & $1000$ & $10000$ \\
    \hline\hline
    \end{tabular}
    \normalsize
\end{table}

The hydraulic backend system is ensured by the numerical model of the water distribution systems, and these models are evaluated in \verb+EPANET+, which is a popular, open-source hydraulic simulator introduced by \citet{Rossmann1993}.
As the main parts of the computer code are written in \verb+Python+, only the core hydraulic calculations are computed by \verb+EPANET+, modification of the problem definition, and querying the results are conducted by the \verb+EPYNET+ interface that is a \verb+Python+ wrapper for \verb+EPANET+.
\verb+EPYNET+ is developed by IT experts at the Dutch drinking water company Vitens and is publicly available by \citet{epynet}.
The dueling deep q-network agent is a pre-packaged one from the reinforcement learning library \verb+Stable Baselines+ \cite{stableBL}.
The authors have developed the environment with which the agent interacts in \verb+Python+, and its interface is compatible with \verb+OpenAI gym+ for the seamless pairing of pre-packaged agents with the same standard interface.
Hyperparameter optimization was carried out with the \verb+Optuna+ package \cite{Akiba2019} for \verb+Python+.

\section{Evaluation and Results}
The history of the rewards and episode rewards is presented first for the D-Town-mod case in Fig.~\ref{fig:trn_rwd_dtown} and Fig.~\ref{fig:trn_ep_rwd_dtown}, respectively.
The history of the rewards and the episode rewards are similar in the Anytown-mod case; thus, those are not presented in the paper due to length considerations.
As the parameters of the action-value estimator are initialized randomly, the agent has no prior knowledge of the problem.
However, some initialization steps are taken by the agent to collect at least as many training samples with which the replay memory can be filled.
The underlying model -- the artificial neural network -- is updated first only after the initialization steps; thus, the training history is not available for the preceding steps.

\begin{figure}
    \centering
    \includegraphics[width=\linewidth, keepaspectratio]{./figures/trn_reward_Dtown.pdf}
    \caption{Reward during the training at D-Town-mod.}
    \label{fig:trn_rwd_dtown}
\end{figure}

Because of the random initialization, different training sessions can produce different results.
Moreover, the exploration capability of the agent introduces some stochastic behavior in the test performance of the agent through training.
Thus, $4$ training sessions are presented for each problem, and the results are smoothed by exponential moving averaging at $\alpha=0.3$.
As the agent starts taking steps according to a random policy, the achieved rewards are low at the beginning of the training.
The randomness of acting is maintained with a decreasing magnitude as the training goes on; thus, the rewards should increase as the action-value estimator is updated with more and more training samples.
Both the immediate and episode rewards saturate at the end of the training indicating that the action-value estimator (the deep q-network) cannot be developed further.
However, the ability of the agent to solve the task cannot be judged by only the training rewards because if the reward is improperly specified, the optimum-looking behavior can be related to a completely false policy from the viewpoint of the original objective.
Evolving false policies is a common problem in optimization tasks as most of them are exposed to the curse of Goodhart's law, see e.g. the work of \citet{Manheim2018} and \citet{Lehman2018}.

\begin{figure}
    \centering
    \includegraphics[width=\linewidth, keepaspectratio]{./figures/trn_episode_reward_Dtown.pdf}
    \caption{Episode reward during the training at D-Town-mod.}
    \label{fig:trn_ep_rwd_dtown}
\end{figure}

Thus, two additional metrics are examined in the agent's performance.
The most important aspect of the DDQN agent's behavior is whether it can reach comparable results to the reference solutions.
This property is examined on the test scenarios by looking at the terminal state value of the episodes denoted by the reference state value of the same scenario.
This is depicted in Fig.~\ref{fig:val_rat_anytown} for the Anytown-mod case.
The value ratio rises above $99\%$ in the first quarter of the training, and it stays there from then on in $3$ of the $4$ training sessions.
It implies that learning the optimal strategy in the sense of reaching high state value is relatively easy in the Anytown-mod case.

\begin{figure}
    \centering
    \includegraphics[width=\linewidth,keepaspectratio]{./figures/value_ratio_Anytown.pdf}
    \caption{Average of the final state value achieved by the agent on the test set during the training at Anytown-mod.}
    \label{fig:val_rat_anytown}
\end{figure}

However, the value ratio does not tell about the number of steps required to reach the terminal state.
Thus, the average length of the test scenarios are depicted in Fig.~\ref{fig:ep_len_anytown} for the Anytown-mod case.
The average number of steps taken before the termination is $6$, which means, that the agent takes at most $3$ active steps on average before terminating an episode with $3$ consecutive siestas.
As a reminder: there is only $1$ pump station in this environment whose speed ratio can be varied between $0.7$ and $1.1$ with increments of $0.05$.
Moreover, the initial speed ratio is randomized uniformly between the feasibility limits; thus a full sweep can be done in $8$ to $11$ steps with these settings (the initialization of the state is not counted as a step) and then the termination could be signaled with $3$ additional siesta steps.
Thus, $3$ active steps in the Anytown-mod means that the agent starts to move in the right direction.
The fact that both the high rewards and the shortest trajectory are achieved at the beginning of the training means that this amount of training steps are unnecessarily high for the Anytown-mod problem.

\begin{figure}[ht]
    \centering
    \includegraphics[width=\linewidth,keepaspectratio]{./figures/episode_length_Anytown.pdf}
    \caption{Average length of an episode in the test set during the training at Anytown-mod.}
    \label{fig:ep_len_anytown}
\end{figure}

The same metrics, i.e., average value ratio and average episode length over the test scenarios during the training are plotted from $4$ different runs in Fig.~\ref{fig:val_rat_dtown} and Fig.~\ref{fig:ep_len_dtown} achieved with the D-Town-mod environment, respectively.
According to Fig.~\ref{fig:val_rat_dtown}, DDQN also performs well in terms of value ratio in the more challenging environment.
However, the absolute performance is slightly lower than in the previous case: the agent achieves at most $98\%$ of the reference state value on average.
The average episode length differs from the previous case, where the shortest trajectories were found soon after the beginning of the training.
In this case, the agent exhausts all its possibilities to move around the parameter space until it learns that choosing a shorter trajectory leads to a higher reward.

\begin{figure}
    \centering
    \includegraphics[width=\linewidth,keepaspectratio]{./figures/episode_length_Dtown.pdf}
    \caption{Average length of an episode in the test set during the training at D-Town-mod.}
    \label{fig:ep_len_dtown}
\end{figure}

\begin{figure}[ht]
    \centering
    \includegraphics[width=\linewidth,keepaspectratio]{./figures/value_ratio_Dtown.pdf}
    \caption{Average of the final state value achieved by the agent on the test set during the training at D-town-mod.}
    \label{fig:val_rat_dtown}
\end{figure}

Finally, Fig.~\ref{fig:val_rat_dtownos} and Fig.~\ref{fig:ep_len_dtownos} represents the same performance metrics for the D-Town-mod environment but with the reference solution acquired with one-shot random trial.
So, in this case, the agent was guided by a suboptimal policy, but it can choose moves off-policy thanks to the fairly high exploration factor at the beginning of the training.
Owing to this, the DDQN agent performs reasonably well even with a poorly chosen reference algorithm as it outperforms that.

\begin{figure}
    \centering
    \includegraphics[width=\linewidth,keepaspectratio]{./figures/value_ratio_DtownOS.pdf}
    \caption{Average of the final state value achieved by the agent on the test set during the training at D-town-mod used one-shot random trial as a reference.}
    \label{fig:val_rat_dtownos}
\end{figure}
\begin{figure}
    \centering
    \includegraphics[width=\linewidth,keepaspectratio]{./figures/episode_length_DtownOS.pdf}
    \caption{Average length of an episode in the test set during the training at D-Town-mod used one-shot random trial as a reference.}
    \label{fig:ep_len_dtownos}
\end{figure}

The average of the value ratios achieved on the test scenarios are summarized in Table~\ref{table:tst_results}.
The performance on the test set is similar to the the performance on the validation set in the end of the training that means the agent is able to generalize among different demand distributions.

\begin{table*}
    \caption{Agent performance on test scenarios}
    \label{table:tst_results}
    \centering
    \small
    \renewcommand{\arraystretch}{1.25}
    \begin{tabularx}{\textwidth}{>{\centering\arraybackslash\hsize=.10\hsize\linewidth=\hsize}X>{\centering\arraybackslash\hsize=.3\hsize\linewidth=\hsize}X>{\centering\arraybackslash\hsize=.3\hsize\linewidth=\hsize}X>{\centering\arraybackslash\hsize=.3\hsize\linewidth=\hsize}X}
    \hline\hline
    \multicolumn{1}{c}{run} &
    \multicolumn{1}{c}{Anytown-mod} &
    \multicolumn{1}{c}{D-Town-mod} &
    \multicolumn{1}{c}{D-Town-mod with one-shot RT} \\
    \hline
    1 & $0.991$ & $0.987$ & $1.04$ \\
    2 & $0.992$ & $0.987$ & $1.04$ \\
    3 & $0.991$ & $0.979$ & $1.04$ \\
    4 & $0.991$ & $0.988$ & $1.04$ \\
    \hline\hline
    \end{tabularx}
    \normalsize
\end{table*}

\section{Conclusions}
A methodology is presented for instantaneous pump control in water distribution systems with reinforcement learning.
The proposed technique is novel in the sense that it relies entirely on live measurement data in the decision-making process.
This property makes real-time optimal control of the pumps possible in a smart water network.

The underlying algorithm of the agent is a dueling deep q-network (DDQN).
As DDQN handles discrete action spaces only, the pump speed setting process is done by changing the speeds with fixed-size increments.

The performance of the agent is measured on two water networks with significantly different number of nodes, pipes and pump stations.
Metrics are defined as value ratio and episode length.
The former measures the value of the terminal state compared to the reference found by the Nelder-Mead method, while the latter shows whether the agent uses the shortest path to the optimum or is it just randomly seeking over the parameter space.
The novel approach is found to be useful, as the agent can achieve comparable results to widespread optimization algorithms even with the restrictions mentioned above.
In the more challenging case, the agent gets close to the reference solution with a difference in value ratio smaller than $2\%$.
As the agent gets an extra bounty if it reaches a value ratio higher than $98\%$, this performance ought to be upgraded by fine-tuning the rewarding function.
Another critical aspect of the agent-based pump control is whether the agent can do its task in a limited number of time steps.
The proposed rewarding algorithm is found to be appropriate to incite the agent to use the shortest possible trajectories towards the optimal state.
In other words, the agent was able to solve the more challenging task in around $20$ timesteps ($n_\mathrm{DQN}$) compared to the $70$ evaluations ($3.5 n_\mathrm{DQN}$) of the Nelder-Mead method and the few hundred evaluations ($\approx\!10 n_\mathrm{DQN}$) of the evolutionary methods.
Besides the number of evaluations, the classic methods rely on the hydraulic model of the water network, but the proposed technique is able to rely only on measurement data after training.

Moreover, the exploratory behavior during the training provides an opportunity to outperform the reference algorithm that is used for the training guide.
This phenomenon is also observed with the presented technique in the case when suboptimal guidance is used.

The proposed technique is a novel operation point optimizer for variable speed pumps in the presented state that depends only on measurement data during control.
It is also suitable for hybrid pump scheduling methods to assign optimal set-points to VSPs in real-time.
As the deep q-network can be changed to handle temporal processes, the presented technique can be further developed to a real-time pump scheduling method for water distribution networks.

\subsection{Data Availability Statement}
\begin{itemize}
\item Some or all data, models, or code generated or used during the study are available in a repository online in accordance with funder data retention policies (Reinforcement learning algorithms implementation \cite{stableBL}, Python API for EPANET \cite{epynet}, original numerical models of Anytown and D-Town \cite{cwsExeter}, the code repository of the presented research including everything to reproduce the results \cite{rl-wds}.)
\end{itemize}

\section{Acknowledgements}
 The research presented in this paper has been supported by the BME-Artificial Intelligence FIKP grant of Ministry of Human Resources (BME FIKP-MI/SC) and by the National Research, Development and Innovation Fund (TUDFO/51757/2019-ITM), Thematic Excellence Program. Bálint Gyires-Tóth is grateful for the financial support of the Doctoral Research Scholarship of Ministry of Human Resources (ÚNKP-19-4-BME-189) in the scope of New National Excellence Program and of János Bolyai Research Scholarship of the Hungarian Academy of Sciences.

\bibliographystyle{elsarticle-num-names}