Unsafe Events Detection in Smart Water Meter
Infrastructure via Noise-Resilient Learning*Ayanfeoluwa Oluyomi† , Sahar Abedzadeh‡ , Shameek Bhattacharjee‡ , and Sajal K. Das†
† Dept. of Computer Science, Missouri University of Science and Technology, Rolla, MO 69401, USA‡ Dept. of Computer Science, Western Michigan University, Kalamazoo, MI 49008, USA
E-mail: {aoonzb, sdas}@mst.edu, {sahar.abedzadeh, shameek.bhattacharjee}@wmich.edu
Abstract—Residential smart water meters (SWMs) collect realtime water consumption data, enabling automated billing andpeak period forecasting. The presence of unsafe events is typicallydetected via deviations from the benign profile of water usage.
However, profiling the benign behavior is non-trivial for largescale SWM networks because once deployed, the collected data
already contain those events, biasing the benign profile. To
address this challenge, we propose a real-time data-driven unsafe
event detection framework for city-scale SWM networks that
automatically learns the profile of benign behavior of water
usage. Specifically, we first propose an optimal clustering of
SWMs based on the recognition of residential similarity water
usage to divide the SWM network infrastructure into clusters.
Then we propose a mathematical invariant based on the absolutedifference between two generalized means – one with positive andthe other with negative order. Next, we propose a robust thresholdlearning approach utilizing a modified Hampel loss function thatlearns the robust detection thresholds even in the presence of
unsafe events. Finally, we validated our proposed frameworkusing a dataset of 1,099 SWMs over 2.5 years. Results showthat our model detects unsafe events in the test set, even whilelearning from the training data with unlabeled unsafe events.
Index Terms—Resilient Machine Learning, Anomaly Detection,
Smart Water Distribution, Smart Living CPS

events in a city-scale SWM network is an open problem. This
motivates our innovative work.
Since a typical city has hundreds of thousands of SWMs,
running a per-meter level anomaly detection technique at theutility head-end will require significant computational andmemory resources, which may not be scalable for large cities.
Moreover, a human operator cannot possibly monitor thehealth of each SWM. This necessitates a lightweight wide-area
anomaly detection technique at a coarse-grained resolution(cluster or zone-level) in real-time without manual supervision.

Fig. 1: Water usage data flow
I. I NTRODUCTION

A. Unsafe Events and Impacts in Consumer SWM Networks

A water distribution network is a complex network of pipes
and joints designed to ensure safe and dependable treatedwater delivery from the source tanks to the end consumers.
However, the reliability can be jeopardized due to various
unsafe events, affecting the safety and consistency of watersupply. Installing smart water meters (SWMs) on individual
households enables visibility of the entire water distribution
network and supports key functions, such as automated billingand usage tracking. Specifically, the SWMs measure hourly
water consumption for every house, and relay the readings toa data concentrator receiving the measurements from SWMs
via wireless connection. Subsequently, such data are forwarded
to the utility cloud (see Figure 1).
Motivation: Unsafe events in the context of this paper
include SWM malfunctions and/or consequent negative values,
large-scale radio frequency (RF) connectivity issues due to
network device malfunctions, SWM software issues, water
theft, or actual negative flow (i.e., water flowing to the opposite
direction). Automated and real-time detection of such unsafe

Let us first enumerate possible unsafe events in the SWM
networks and their impacts on the utilities and customers.
Water Leaks: This event shows up as higher water consumption readings from a group of smart meters. Whilehigh water usage can also be due to legitimate events (e.g.,
irrigation, gardening, or use of multiple water outlets), theease of detecting leaks depends on the severity of the damage
and the underlying causes, such as aging household piping,sabotage of pipes connecting a group of households, andnatural disasters. Therefore, real-time detection of water leaks
is highly important to the public water utilities. This is because
undetected water leaks not only pose a financial burden on
the consumers, resulting in higher bills, but also a reputation
loss for the utility. Moreover, unaccounted water usage (due towater theft or leaks) leads to revenue loss and hence financial
sustainability. The infrastructure (mainly pipes) will also be
progressively damaged if the leak is not detected or handled
early enough.
Backward Flow of Water: Backward flow or backflow occurs when the water pressure drops suddenly or completely
fails, thus allowing contaminated water, liquids, or gases to*A major part of the work was done when A. Oluyomi visited S.
Bhattacharjee at Western Michigan University in 2023.

flow back into the clean water supply. Negative water pressure
can be caused by a sudden change in the flow velocity, which
can be due to the closing or opening of valves, starting or
stopping a pump, or changes in the pipe geometry. Such an
event shows up as negative water consumption from various
smart water meters. Heavy rainfall can cause flooding, which
can form puddles around the sprinkler devices on lawns/farms.
During a pressure reversal, the contaminated water containing
fertilizer, pesticides, etc. can be drawn back into the cleanwater supply line [8], leading to a public health concern. The
dataset documentation in [4] provides further explanations fornegative readings (e.g., organized water theft and SWM faults),
which need to be detected promptly.
Sustained Missing Data: Our analysis of the SWM dataset
revealed that many SWMs did not send water consumption
data for several hours; then a single large consumption reading was sent corresponding to the total consumption during
the outage. Such events can be due to SWM malfunctions,intermittent RF connectivity, and third-party SWM software
issues, among others. Long periods of missing data impact
the water utility operations by impeding the accuracy of datadriven processes (e.g., forecasting, billing, water conservation
techniques), the service quality, and operational efficiency.
B. Aims, Objectives, and Challenges
To detect unsafe events in SWMs, we used hourly-recorded
water consumption data from 1,099 houses [3]. Since different
houses exhibit distinct water consumption patterns, our firstobjective is to design a clustering algorithm that can precisely
identify and group houses based on the inter-dependence ofsuch patterns. Water consumption in a household depends on
multiple factors, such as the number of occupants, the natureof their employment, and individual preferences, which makethe clustering method non-trivial. However, once such clusters
are identified, we run an instance of the anomaly detection
technique on the aggregate raw data from the entire cluster,
ensuring the scalability of unsafe event detection. Once therelevant “cluster” originating the anomaly is identified, theutility can opportunistically invoke a more costly per-meter
anomaly detection technique or a manual inspection only for
those houses belonging to that identified cluster.
The most important technical challenge is the presence of
unsafe events in the training dataset collected from actual
deployments as this complicates the accurate characterizationof the benign profile of any smart living cyber-physical system(CPS) such as the SWM networks. Hence, the objective is toautomatically learn the benign profile of the operation whilelearning the thresholds of benign operation characterizing theevent detection criterion, even with the inclusion of unsafe
events in the training dataset. Our research contributions aim
to address these unique objectives and challenges.
C. Contributions of This Paper
This paper proposes a novel time-series event detection
framework for large-scale water distribution networks. The
major contributions are summarized below.

(1) We design a clustering algorithm that partitions theresidential smart water meters into smaller clusters using similarities in their consumption behavior. Creating such clusters,
where the data in a cluster are highly positively correlated witheach other, maximizes the invariant of the anomaly detection
metric characterizing the benign behavior. Our approach uses
mutual information that measures how much water meters areconnected to each other owing to their complex, non-linear
relationships. A graph-theoretic technique is used to partition
the houses into clusters with strong correlations. The clustering
approach promotes decentralized cluster-wise implementation,
enabling straightforward identification of the cluster where an
unsafe event has occurred. Subsequently, each smart water
meter data within that cluster can be inspected individually.
(2) We propose an unsafe event detection model that leverages the absolute difference between the negative and positive
orders of the generalized mean. Our findings demonstratethat this invariant provides optimal sensitivity to changes
induced by unsafe events while remaining insensitive to benign
fluctuations. The cyclo-stationarity of the human behavior inwater consumption ensures the effectiveness of this approach.
(3) We develop a robust method to learn the thresholds ofbenign (normal) operations even when unsafe events occur
during the training process. This method uses the Hampel
three-part estimator to automatically learn unbiased thresholds
that characterize normal behavior, even when unsafe eventsaffect the invariant’s structure during training. This contributes
to the model’s ability to distinguish normal behaviors fromunsafe events in the test and deployment stage.
(4) We validate the proposed framework using the SWM
dataset from Alicante, Spain. Experimental results demonstratethat our framework effectively detects unsafe events. The
performance is evaluated by comparing its decisions withmanually labeled water usage data, serving as ground truth.
D. Design Rationale and Theoretical Intuition
Before proceeding further, we outline the design rationale
and theoretical intuition to shed light on the essential properties required for our framework to be effective in any CPS.
(a) The framework is relevant to smart living CPS applications, where the data from the physical world is readily
affected by human behaviors, not just by the laws of physics.
The clustering method combined with the invariant helps
abstract this randomness into a time-invariant latent space.
(b) Our invariant design is such that it ameliorates theneed to manually remove the negative and zero values before
learning the benign profile essential in the existing time series
invariant using Pythagorean Means [12]. Additionally, our proposed invariant contains properties that cause deviations underunsafe events due to the asymmetry of convexity properties in
the generalized means with a negative order, as opposed to thegeneralized means computed with a positive order.
(c) The robust threshold learning involves a ‘modified’Hampel loss function (a composite loss function) rather thanthe common practice of a regular loss function, which gives
completely different fitness values depending on the error. This

explains why the learnt threshold is accurate in presence of unlabeled, unsafe events. Typically the inputs to a loss function ina regression-type learning problem are directly the regression
errors. However, since the reduction of false alarms is theprincipal design challenge in time-series anomaly detection,
we add more weights to the regression errors on moderately
outlying points rather than in-lying points before passing itto the Hampel loss function that allows the threshold learner
to be not too conservative. This explains good performance infalse alarms, unlike the usual outcome where the use of robust
loss functions comes at the cost of increased false positives.
(d) Our proposed framework will have less success if apositive covariance structure is absent (affecting the invariant
stability) under feasible spatio-temporal granularities, and if awide sense cyclo-stationarity is not observed in the extracted
time-series invariant (affecting the safe margin generation).
The rest of the paper is organized as follows. The clustering
algorithm is proposed in Section II along with the model ofwater consumption behavior. Section III describes the anomaly
detection framework while IV presents the experimental results. Section V offers conclusions.
II. P ROPOSED C LUSTERING A LGORITHM
This section explores the need for clustering and detailits design, implementation, approach, and execution process.
The clustering is done on a subset of data that are handcleaned, and this process is done only once. Once the clusters
are established, our invariant detection metric is calculated
separately per cluster using raw uncleaned data as collectedfrom the actual SWM from the wild. This reflects the reality ofthe problem of training to learn the profile of benign behavior
in the presence of unlabeled unsafe events without human supervision. Then the thresholds of benign behavior are learned
separately for each cluster. Next, in the test set, the sameinvariant is applied and compared with the thresholds learned
during the training stage that specifies benign operation. If theinvariant in the test violates the learned threshold, the proposed
framework declares an unsafe event, else it infers benign.
A. Understanding Water Consumption Patterns
We used the Hellenic Data Service [3] containing hourly
water usage data from 1,099 residential houses over 2 years
and 5 months (January 2015 – May 2017). We used a 60-20-20
split for the training, cross-validation, and testing sets.
Fig. 2 shows that the mean water usage for corresponding
hours of the entire training set is highly unpredictable. Hence,
popular methods such as the cumulative sum (CUSUM), which
requires prior knowledge of an in-control mean serving asa benign baseline, fail to solve this problem (see the proof
in Sec. IV-D). As illustrated in Fig. 2, each line represents
the average daily usage trend for an individual house (only
30 houses are displayed for brevity) over the hours of aday. Furthermore, it is important to highlight that the actual
mean usage within a specific hour of different days exhibits
substantial variations. Consequently, these fluctuations significantly influence the traditional metrics like the moving
averages and simple means. Hence, the need for anomaly

detection metrics that consider these variations, while profilingthe benign behavior and establishing a detection criterion.

Fig. 2: Mean of water usage for corresponding hours
B. Residency Similarity Recognition (RSR) Algorithm
We observed a notable correlated increase in water usage,
particularly during the early hours of the day across mosthouses, indicating the period when most household occupants
get ready for the day. Similarly, the correlated drops in usage
during the night time shows correlated behaviors in general.
1) Need for Clustering: Even under benign conditions,human behaviors are complex and non-linear. We need to
identify subgroups where the correlation is high and treat themas single points of observation for visibility of the entire SWM
network. This is why the clustering makes sense. The reason
behind considering correlation as the basis for clustering is that
this is a required property to provide stability to the invariant
metric inspired by our work [1]. We showed that the ratio of
Pythagorean means (e.g., harmonic mean to arithmetic mean
ratio) gives a stationary time series if the data comes from a
network of highly correlated random variables.
The strength of the positive correlation between water
usage of houses becomes low when considering a randomly
selected section of the city or the entire city. Low correlation
adversely affects the effectiveness of the anomaly detection
model because the derived invariant is not stable under benign conditions. Hence, clustering smart meters based on thetemporal information led to an intelligent clustering algorithm,
called the Residency Similarity Recognition (RSR) algorithm.
The primary objective of the RSR approach is to partition
the houses into highly correlated/similar clusters. We use
mutual information (MI) to capture similarity between houses
due to non-linear dependence between the data coming fromdifferent meters. The RSR also enhances the robustness ofthe detection model by improving the invariant’s stationarity(see Section III-A) and facilitates decentralized monitoring andlocalization of unsafe events to specific clusters.
2) Measuring Dependence between Meters: We use mutual information [6] to determine the dependencies among
SWM meters. The MI quantifies the knowledge gained about
one variable given another variable’s value. Unlike popular
correlation techniques such as Pearson correlation [17], theMI captures both non-linear and linear dependencies between
variables, enabling a better understanding of smart living CPS,
such as SWM networks. We calculate the MI as:
Z Zp(wi , wj )I(Wi ; Wj ) =p(wi , wj ) log
dwi dwj (1)p(wi )p(wj )where wi and wj are water usage for houses i and j overan hourly time window and I(Wi ; Wj ) is the measure of
MI between two house water usage wi and wj . The MI
between any two variables is symmetric. After measuring the
MI between every pair of houses, we use a graph-theoretic
approach to obtain the correlated clusters. An MI close tozero indicates independence between variables while higher
MI values indicate stronger mutual dependency.
3) Formulation as a Graph Problem: Assuming (theoretically) every house may have some dependencies with otherhouses, the initial formulation is an undirected complete graph,
G = (V, E), where v ∈ V is the set of vertices (i.e., SWMs)
and the e ∈ E is the set of edges denoting the dependency
between two vertices. The MI value, i.e., we = I(Wi ; Wj )between any pair of vertices is the edge weight denoting thestrength of that dependency.
The goal is to find a set of correlated clusters such that themembers within the cluster high dependency. Therefore, weassign a minimum acceptable MI value wcut such that only
edges with we > wcut are feasible, and the rest of the edges
are deleted. Note that wcut is a crucial hyperparameter, since a
low wcut results in a graph with more edges with low weights.
Thus houses having low dependencies/correlations with otherhouses will be clustered together, which will weaken theinvariant of the detection metric for that cluster. A weakinvariant under benign conditions is not conducive to thedetection of unsafe events or successful profiling of the benign
behavior of a cluster. Conversely, a high wcut will remove
a lot of edges with considerable dependency strengths. This
will cause a significant number of SWMs to be not included
in any of the clusters. Section II-C discusses the method
employed to obtain the optimal wcut along with N , whichis the maximum allowable number of houses for a single
cluster. Next, we define some key terms and then describe
the clustering algorithm.
A Cluster ck is the k-th cluster containing a set of vertices
with high similarity of water usage. The clusters are formed byincrementally adding vertices one at a time so as to maximize
the edge weight we of a cluster. While forming a cluster, thecurrent collection of vertices is the candidate cluster c. We
need the following definitions inspired from [5] to understandthe clustering.
Volume is defined as the sum of all edge weights we inside
a candidate cluster c, where each edge e = (vi , vj ) ∈ E has
both vi and vj in c, Thus,X
vol(c) =we(2)
vi ,vj ∈c

Sumcut is the sum of weights of the edges that connect cto the rest of the graph. Each edge e = (vi , vi′ ) ∈ E has one
vertex vi inside c and the other vertex vi′ outside c. Thus,
X
sumcut (c) =wei,i′(3)
vi ,vi′ |vi ∈c∧vi′ ∈c/4) Description of Clustering Algorithm: Algorithm 1 returns a set C of clusters. The algorithm starts by verifying that

Algorithm 1: Algorithm for Clustering Houses
Data: Edge-list (ei,j = {Vi , Vj , wei,j })Result: Clusters C = c1 , . . . , ck
1 while edgelist ̸= ∅ do
2
Select random ei,j such that ei,j ∈ edgelist;3c ← vi , vj ;4vol(c) = wei,jP;
5sumcut (c) = vi ,v′ |vi ∈c∧v′ ∈cwei,i′ ;ii/6edgelist\ = {(vi , vj )};7
while vol(c) ≤ sumcut (c) and count(c) ≤ N do
8Neighbors(c) ← {ei,i′ | (vi , vi′ )};′9// vi ∈ c, and/c
P vi ∈
10
sumcut (c) = e ′ ∈Neighbors(c) wei,i′ ;i,i11// Sum of edge weights of vertices
that connect c to the rest of thegraph
12vi′ +1 = max[wei,i′ suchthat ei,i′ ∈ Neighbors(c)];13c ← vi′ +1 ;14
vol(c)+ = wei,i′ +1 ;15
edgelist\ = vi′ +1 ;end
16
return ck ;17
remove all (v, e) ∈ c from edgelist;end
18 return C = c1 , . . . , ck ;Fig. 3: A cluster formulation with RSR
the edgelist is not empty (line 1). Using Fig. 3 as an example,
all edges in Fig. 3a have we > wcut . The initial edge v1 , v2 israndomly selected from the edge list to establish the process ofbuilding the first cluster that contains vertices v1 and v2 (lines
2 and 3). Then the algorithm calculates the candidate cluster’svol(c) and sumcut (c) (lines 4 and 5). In Fig. 3b, vol(c) is
we of edge v1 v2 , while sumcut (c) is the sum of we for edges
v1 v4 , v1 v5 , v2 v5 , and v2 v3 . Vertices v1 , v2 are removed fromthe edgelist to prevent duplication (line 6).
The algorithm proceeds iteratively from lines 7 to 15 whilevol(c) ≤ sumcut (c) and the number of vertices, count(c)within the candidate cluster c does not exceed N . Here,N eighbors(c) is the set of all vertices that connect thecandidate cluster to the rest of the graph (line 8). Then, we
calculate the sum of we of all edges that connect the candidate
cluster to the rest of the graph, i.e., the sum of all we in
N eighbors(c) (line 10).
Next, the algorithm selects the edge with the highest wefrom N eighbors(c) (line 12), and includes the adjacent vertex
vi′ +1 (currently outside the candidate cluster) into that candidate cluster (line 13). Then, we of edge (vi , vi′ +1 ) is added
to vol(c) (line 14). This is illustrated in Fig. 3b and Fig. 3c,where vertex v5 is included in the cluster since the edge v1 v5
has the highest we . Then vol(c) is calculated as the sum of
we for the edges v1 v5 and v1 v2 .
The sumcut (c) is updated as the sum of we for edges v1 v4 ,v2 v3 , and v5 v6 . Subsequently, all rows where vertex vi′ +1
and all its adjacent vertices are removed from the edgelist to
prevent duplication (line 15). Once any of the two conditions
on line 7 is violated, the algorithm returns cluster ck alongwith the vertices formed in the process. For the example inFig. 3c, vertices v1 , v2 , v5 form the first cluster since the sumof we of v1 v2 and v1 v5 is greater than the sum of v1 v4 ,v2 v3 , and v5 v6 . The vertices and adjacent vertices are removed
from the edgelist (line 17) and the iterative process continues.
Once the edge list becomes empty, the set C consisting of
all c1 , ..., ck clusters formed is returned (line 18) and thealgorithm terminates.
5) Using the RSR Algorithm with Cleaned Data: The RSR
algorithm aims to create clusters that contain houses withhighly correlated water usage patterns. We employ a handcleaned version of the dataset and apply the RSR algorithm
to create a one-time cluster formation. In a real-life scenario,
the clustering only has to be done only once in a few years
to consider water patterns that might have changed over theyears. After the clusters are formed, the raw uncleaned datafrom the corresponding SWMs within each cluster is usedin the rest of this paper. The use of an uncleaned versionof the data preserves all the challenges in the problem of‘automatically’ learning the profile of benign behavior in thepresence of unsafe events from a real smart living CPS. The
benign profile and its thresholds, if learned accurately, can beused for detecting unsafe events in the deployment/test set. The
process of cleaning the data is explained in Appendix VI-A.
C. Hyperparameter Learning for the RSR Algorithm
We derive the optimal wcut and N by finding which∗
combination of the candidate wcut
and N ∗ maximizes thestrength of stationarity in the time-series data for all theclusters. To quantify the strength of stationarity, we use an
Augmented Dickey-Fuller (ADF) [11] test with the timeseries data ad(tϕ )ck = |HM (tϕ )ck − AM (tϕ )ck | (see Appendix VI-B for details). Here ad(tϕ )ck is obtained using aninvariant proposed in our earlier work [2], which however fails
under uncleaned SWM data containing unsafe events, thereby
reducing the detection accuracy.
Now, let ∆ytck be the ADF test statistics for a cluster ck
and ∆avg (C) the mean of ∆ytck produced for a given set C.
Furthermore, let Υ(C ′ ) be the number of houses that do notbelong to any cluster ck (after C has been obtained) becauseof their low edge weight to adjacent houses.
By jointly minimizing a weighted sum, RSRmin =w∆ ∆avg (C) + wΥ Υ(C ′ ), we obtain the optimal values of
wcut and N (see Algorithm 2). Here, it is crucial to note thatw∆ > wΥ , since we want more importance to be placed on∆avg (C). This choice reflects the importance of stationarity
over the number of houses covered in our optimization process.

Algorithm 2: RSR Hyperparameter Learning
Result: N , wcut∗
for N ∗ , wcut
do∗∗for w∆ , wΥdo
2
Obtain set C with Algorithm 1;c3∆avg (C) = mean(∆yt k );4
Υ(C ′ ) = number of houses not clustered;
5
RSRmin ← w∆ .∆avg (C) + wΥ .Υ(C ′ )endend
6 N, wcut , w∆ , wΥ = arg minN ∗ ,w∗ ,w∗ ,w∗ (RSRmin);
cut
1

∆

Υ

∗
For each combination of wcut
and N ∗ , a unique set ofclusters C = c1 , ..., ck , is produced, and we calculate the
RSRmin for that set C. Finally, we identify the global minima
of RSRmin to determine the optimal values of wcut = 0.061,
N = 80, w∆ = 0.99, and wΥ = 0.01. The optimization∗
process involves varying wcut
from 0.001 to 0.205 with a step∗
size of 0.005, and N from 50 to 100 with a step size of 5 toobtain the optimal values for wcut and N .
Clustering with Optimal Hyperparameters: Using the optimal values, we found 14 clusters, out of which the last 3
clusters had only two houses, making the event detection
trivial. Hence, we report the performance on the remaining
11 clusters. We provide cluster-wise test statistics results inAppendix VI-C. We also presented a comparison of the performance of RSR algorithm with other conventional algorithms
in Appendix VI-D.

III. U NSAFE E VENT D ETECTION M ETHOD
In this section, we first present our novel invariant design,
followed by a method that learns a robust detection thresholdthat can successfully detect unsafe events in the test set, andfinally the details of optimizing the hyperparameters of theinvariant design and the detection threshold learning method.
A. Characterizing Benign Invariant
As outlined in section II, the clustering process produces
clusters that maximize the correlation between houses included
within the cluster. In this section, the uncleaned raw datacollected from the wild are used for each house.
Consider a cluster ck comprising n(ck ) houses. Next, we
compute an invariant metric ADck (tϕ ) for each cluster ck at
every time slot tϕ , where ϕ represents the hour of the day
such that ϕ ∈ {1, 24}. For the rest of this section, we remove
the symbol ck from all notations/equations, since the detection
model runs per cluster ck in the same way.
Let wi (tϕ ) be the water usage recorded by the i-th SWM atthe time slot t corresponding to the ϕ hour of the day. Thenthe Generalized Mean (GM) value from all wi (tϕ ) from nmeters in the cluster can be calculated as: ! 1nGM (p) (tϕ ) =p
1 X
wi (tϕ )
n i=1

p

(4)

where the parameter p ∈ R denotes the order of the GM
statistic. Instead of GM, we propose a new AD(tϕ ) invariant,defined as the absolute difference between the two generalized−+
mean values GM p (tϕ ) and GM p (tϕ ), written as:−+

AD(tϕ ) = |GM p (tϕ ) − GM p (tϕ )|(5)−such that the first term GM p (tϕ ) has a strictly negative orderparameter such that p− ∈ R− \ −1, while the other term+
GM p (tϕ ) strictly has a positive order of p, where p+ ∈ R+ .
This equation can be expanded into the following:AD(tϕ ) = 


 1
 1
n 
n 
p− p−
p+ p+
1 X
1 X


−
wi (tϕ )wi (tϕ )
n i=1

(6)

n i=1

The orders p+ and p− is chosen by hyperparameter tuning
such that it magnifies deviations during unsafe events. The−1
GM p i.e., HM , do not exist under zero and negative values,and hence is not a feasible value of p− . An explanation is givenin the appendix VI-E, validating the observed deviation in theinvariant whenever an unsafe event occurs.
The optimization approach employed to determine p+ and−p is outlined in section III-E. To emphasize the importanceof the clustering algorithm and demonstrate the robustness of
the invariant under benign conditions, we illustrate a scenario
in Fig. 4a displaying the behavior of AD(tϕ ) for a cluster
where an unsafe event has taken place in 93rd and 101st hour.
It is clear from Fig. 4a, that our invariant shows a deviation
matching the 93rd and 101st hour. In contrast, in the absenceof clustering and previously proposed invariants (See fig. 4b),no considerable deviation is present.

(a) AD(tϕ ) with our clusteringobtaining a range that wraps around each hour as seen in
Fig. 5.
1) Expected Value of Invariant per Slot: Given tϕ as aunique hour slot in a day, the expected value for the invariant
E[AD(tϕ )] is defined as the average of the ADtϕ values
at each corresponding hour tϕ (i.e. same ϕ) over the entire
historical training data, mathematically written as:Ph(ϕ) AD(tϕ )(7)
E[AD(tϕ )] =h(ϕ)(b) ad(tϕ ) without clustering

Fig. 4: Effect of cluster level correlation on the invariance.
However, the AD(tϕ ) is also vulnerable to occasional spikes
under benign conditions; which reduces detection confidence.
Hence, we implement a two-tier approach that smoothens
the invariant into a stateless residual space, that provides
the ability to improve detection accuracy while simultaneously reducing the false alarms. We generate safe margins
around AD(tϕ ) (Sec. III-B), resulting in stateless residuals(Sec. III-C), which are used to learn the appropriate robust
threshold in Sec. III-D.
B. Safe Margin Generation
Stateless residuals measure the difference between the expected range of the invariant value at a time window t (calledsafe margin) and the actual invariant value at the same time
window. We devise a method for determining the safe margin
through temporal reasoning. This entails calculating certain
indicators of central tendency and variability of the invariant
within a specific temporal context. This idea is based onwhere h(ϕ) is the total number of data entries corresponding to
the ϕ-th hour of any given day over the historical training set.
In our training dataset h(ϕ) is 522, and since ϕ corresponds toeach hour of the day, we end up with 24 entries of E[AD(tϕ )],one corresponding to each hour.
2) Safe Range of the Invariant: Now, we determine a safe
margin around each E[AD(tϕ )] (location parameter) by multiplying a scalar factor ϵ to the median absolute deviation of all
samples in AD(tϕ ) (denoted as M AD(AD(t)) in the training
set regardless of the hour ϕ. We use MAD as a measureof dispersion around the E[AD(tϕ )] given its robustness tooutliers when compared to the standard deviations. The ϵprovisions the trade-off that allows for benign fluctuations in
the metric to be tolerated. The resulting safe margin consists ofboth the upper margin Γhigh (tϕ ) and lower margin Γlow (tϕ )corresponding to any given hour ϕ of the day:Γhigh (tϕ ) = E[AD(tϕ )] + ϵ.M AD(AD(t))(8)Γlow (tϕ ) = E[AD(tϕ )] − ϵ.M AD(AD(t))(9)

Finding the optimal ϵ is critical for minimizing false alarms(FA) while concurrently controlling missed detection (MD).
This hyperparameter is determined in Section III-E. A large
ϵ results in a wider range of safe margins, leading to more
missed detections (MDs). Conversely, a smaller ϵ yields anarrower safe margin, resulting in an increase in false alarms.
Figure 5a shows the AD(tϕ ) and the safe margins for a test
set cluster after obtaining the safe margins from the training
set. We only show t for 8 days (192 time slots).

(a) AD(tϕ ) samples shown continu- (b) Stateless Residual shown continously for 8 days with the upper and uously for 8 days with the optimal
lower safe margins
learnt detection standard limits

Fig. 5: Snapshot of the Detection model.
C. Formation of Stateless Residuals
As defined in Eqn. 10, we compute a signed residual distance RU C(t) (Residual Under Curve (RUC)) by subtracting
the observed sample AD(tϕ ) at time t for a specific hour
ϕ from the corresponding safe margins for that hour i.e. theΓhigh (tϕ ) and Γlow (tϕ ). The ϕ was omitted from subsequent
equations because we are not particular about the specific hourof the day once the RU C(t) has been obtained (See Fig. 5(b)).
Mathematically,AD(tϕ ) − Γhigh (tϕ ),RU C(t) = AD(tϕ ) − Γlow (tϕ ),0,

if AD(tϕ ) > Γhigh (tϕ );if AD(tϕ ) < Γlow (tϕ );otherwise(10)

Specifically, Eq. 10 can theoretically produce zero (0), positiveor negative values. Formally, the set of positive, negative, andzero stateless residuals is denoted as RU C + (t), RU C − (t),and RU C 0 (t) respectively.
D. Learning Resilient Thresholds of Stateless Residuals
The RU C(t) is obtained from raw uncleaned data that
already contains the unsafe events. Unless painstaking manual
labeling is done to identify and remove RU C(t) samples corresponding to the unsafe events/errors, the thresholds identified
by typical regression or even max(residual) [13], will create
biased thresholds that are unable to detect most events.
Therefore, we need to design a robust learning mechanismthat automatically learns a correct event detection thresholdwithout manual cleaning of abnormal events. Our method
combines weighted regression with Hampel’s three-part redescending M-estimator. Formally, to learn the upper threshold- the inputs to the regression problem are the RU C + (t)values, the candidate upper threshold hypothesis space isparameterized by τ + , and the optimal output learned thresholdis τmax . Similarly, to learn the lower threshold, the inputs are
RU C − (t), the candidate threshold hypothesis space is τ − ,and the optimal lower threshold is τmin . Since the method is
the same for learning both thresholds, henceforth we use just
RU C(t) and τ notations for brevity.
The generic regression error is denoted as r(t) = RU C(t)−τ . Now we make a design change by not providing r(t) as
direct inputs to the Hampel loss function. Instead, we defines(t) = r(t).λ, where s(t) is a weighted regression error suchthat, λ = w2 ∀ r(t) ≥ 0 and λ = w1 ∀ r(t) < 0 and
w1 < w2 . As seen in Algorithm 3, r(t) ≥ 0 is checked, thenr(t).w2 is checked with the different a, b, and c values suchthat less importance are placed on the bigger outliers, while
the weights provide the trade-off that reduces false alarms.
Hampel’s Three-part Loss Function (HTF): is one of therobust estimators for outlier resilient learning in the broad classof M-estimators. The outliers in this context are the presenceof RU C(t) samples triggered by unsafe events already present
in the training data that is used to learn the profile of benign
operation. For weighted regression error s(t), the HTF loss
function is defined as:

ρ(s(t)) = (s(t))2

 2


2

a|s(t)| − a22a(c−b)
ab − a2 +

2




2
a(c−b)ab − a2 +2

if |s(t)| ≤ a
if a < |s(t)| ≤ b1−

c−|s(t)|c−b

2 , if b < |s(t)| ≤ c
if |s(t)| > c(11)

the a, b, and c are different threshold values of the regression
error such that 0 < a ≤ b < c < ∞, such that thecorresponding loss function value depends on the magnitudeof s(t). The HTF has a special property, such that the firstorder partial derivative of the estimator can be customized togradually descend to a zero. This unique characteristic ensuresthat moderate and large outliers are entirely disregarded whilesmaller outliers are not entirely disregarded, thus, significantly
enhancing the efficiency of the re-descending M-estimator bybalancing legitimate benign fluctuations and large RU C(t)spikes that are due to unsafe events [7].
Why HTF for Training Data: High values of RU C(t) areless frequent (highly unsafe events) compared to moderatehigh values to low values signifying more benign conditions.
Thus, high s(t) should contribute less to the threshold. Thus,
the range |s(t)| ≤ a is assigned the most importance (first partof Eqn 11). All the other parts of Eqn. 11) are such that theygive decreasing importance with an increase in the s(t). The
next task is to determine the appropriate a, b, and c values inthe HTF function.
Obtaining a, b and c in HTF: The mechanism to find a, b,
and c is data-specific unlike the other parts of the paper.
The K-means clustering algorithm was employed to categorize
the RU C(t) values into distinct clusters. K-means’ ability
to cluster without preconceived assumptions about the data’s
distribution made it an ideal choice for this purpose. The
optimal number of k-means clusters was determined to beeight, as it effectively separated residual data into 8 clusters.
Since the HTF requires four distinct regions, we need to
combine the clusters into 4 segments demarcated by 3 values
a, b, c. To ensure the robustness of the borders, we computed
the MAD of two preceding and two succeeding clusters.
Obtaining Thresholds τmax and τmin : After determining
the components for HTF, we use algorithm 3 to calculate
the noise resilient upper limit τmax which searches among
RU C + ≥ 0. To obtain the lower limit τmin , we searchedamong RU C < 0, and fig. 5b, shows the learned thresholds.
E. Hyperparameter Optimization
Next, we show how to find the optimal hyperparameters p+ ,p , ϵ for the robust threshold τmax and τmin . As discussed
in Appendix VI-E, when an unsafe event occurs, a significant
deviation in the invariant can be guaranteed with the propertiesof the Schur concave and Schur convex. Therefore, we varied
p+ and p− from 1 to 5 and -2 to -6 respectively with astep size of 0.5, omitting p− = −1 as it corresponds toHarmonic Mean that does not work with negative readings(a result of backflow events). We also varied the ϵ from 0
to 4 with a step size of 0.05. We generated different setsof RU C(t) and learned the threshold τmax and τmin . To
obtain the optimal combination of these hyperparameters, weformulated the optimization problem as an error minimization
problem aimed at minimizing the number of False Alarms(FA) and Missed Detection (MD) for each combination of
the parameters such that it minimizes wf a .F A + wmd .M D,
such that wf a = 0.6 and wmd = 0.4 to emphasize the baserate prior between benign and unsafe events. We assign moreweight to FA as occurrences of unsafe events are infrequent−

Algorithm 3: Robust Learning of τmaxck+

Data: [RU C (t), [τ ], , a, b, c, w1, w2]
Result: τmax
for τ + do
if r(t) ≥ 0 thenif (|r(t)| · w2) ≤a then
cost ← 12 · (r(t) · w2)2
end
else if a < (|r(t)| · w2) ≤ b thencost ← (a · (|r(t)| · w2)) − (0.5 · (a2 ))endelse if b < (|r(t)| · w2) ≤ c thencost ←
c−(|r(t)|·w2) 2(a · b) − (0.5 · (a2 )) + (a · c−b) )2 ) · (1 − (c−bendelse if (|r(t)| · w2) > c thencost ← (a · b) − (0.5 · (a2 )) + (a · c−b
2 )
end
end
else
if (|r(t)| · w1) ≤ a then
cost ← 0.5 · (|r(t)| · w1)2
end
else if a < (|r(t)| · w1) ≤ b thencost ← (a · (|r(t)| · w1)) − (0.5 · (a2 ))endelse if b < (|r(t)| · w1) ≤ c thencost ←
c−(|r(t)|·w1) 2) )(a · b) − (0.5 · (a2 )) + (a · c−b2 ) · (1 − (c−bendelse if (|r(t)| · w1) > c thencost ← (a · b) − (0.5 · (a2 )) + (a · c−b
2 )
end
end
end
τmax = arg minτ + (cost) ;compared to benign. The wf a and wmd in our framework
should be modified accordingly for other systems. Algorithm 4
formalizes the hyperparameter learning.
Algorithm 4: Model Hyperparameter Learning
Data: Data clusters
Result: p+ , p− , ϵfor each cluster ck do
for p+ , p− do−+

Calculate AD(tϕ ) = |GM p (tϕ ) − GM p (tϕ )|;for ϵ doGenerate Γhigh (tϕ ) and Γlow (tϕ ) ;Generate new RUC(t);Robust learning for τmax and τmin ;error ← wf a .F A + wmd .M D;
end
end
−
p+
∗ , p∗ , ϵ∗ = arg minp+ ,p− ,ϵ (error);end

F. Detection Criterion in Test Set
This section discusses the process of arriving at the predicted events using the two-tier detection approach.
1) Obtaining RUC(t) in Test Set: Let AD(tc (ϕ) be theinstantaneous current value of our invariant obtained from theobserved water usage in cluster ck at time tϕ in the test set.
Now we use the 24 distinct values of Γhigh (tϕ ) (fromEqn. 8) and Γlow (tϕ ) (from Eqn. 9) known from the training
phase, and compare it with the current AD(tc (ϕ) sample
using the Eqn. 10 by replacing AD(tϕ ) with AD(tc (ϕ), which
produces the current RU C(tc ) sample.

2) Detection Decision: We flag an unsafe event if the RUC
gotten from the safe margins (tier 1) is not within the lower
and upper threshold τmin and τmax respectively. This is given
in equation 12 and shownwith figure 5b.
(
RU C(tc ) :∈ [τmin , τmax ], No Events;̸∈ [τmin , τmax ], Unsafe Event.

(12)3) Predicted Events (Z(t)): We obtain the predicted eventat time window t, denoted as Z(t) calculated as the difference
between RU C(tc ) and the standard limits τmax and τmin .
ccRU C(t ) − τmax , if RU C(t ) > τmax ;ccZ(t) = RU C(t ) − τmin , if RU C(t ) < τmin ;0,otherwise(13)4) Ground Truth Formation: We manually identified theunsafe events described in section I-A by rigorously studying
the data and referring to available documentation for bothcleaned and uncleaned versions of the data. It should be notedthat the different events explained in Section I-A have two
effects - which is either a large positive reading (+ve Events)
or a negative reading (-ve Events) from one or more SWMs.
For sustained missing data, we observed that a huge positive
reading always follows the missing data, thus, the effect of thisis seen as a moderately large positive number (+ve Events).
To generate a ground truth events dataset Dgt from the raw
data Draw we applied the following rules/conditions to label
the ground truth of events:
2 if ∃ i , wi (tϕ ) > θhζ(t) = 3 if ∃ i , wi (tϕ ) < θl
0 otherwise(14)

where θh and θl are a certain high water reading threshold
and a certain low water reading threshold respectively, andζ(t) is the time-stamped event label. We now discuss how wedetermined θh and θl for the SWM application.
Determining θl : Practically speaking, water usage should
be strictly positive, thus, we postulate that any water usage
significantly less than 0 is an unsafe event, hence θl = −1.
Determining θh : Determining the appropriate water usage
per hour is a complex task due to varying water usage per hour,
day, week, and season, as well as household occupants andbuilding types. Hence to obtain θh , we employ Chebyshev’s
Inequality on cleaned water usage data (the cleaning method is
explained in appendix VI-A). We used a clean dataset because
we can obtain the appropriate standard deviation (σ) and mean(µ) of the actual water consumption unlike the uncleaned datathat contains large positive numbers (up to 13 million liters),
these values influence the σ which in turn leads to taking larger
step sizes to obtain θh .
Chebyshev’s Inequality is represented as P (|X − µ| ≥gσ) ≤ g12 , where, X in our context is the water usage, g isscalar multiplicand of the standard deviation σ. This inequality
provides a reliable means of setting the upper bounds on the
proportion of data that can deviate significantly from the mean.
With this, we calculate θh = µ + gσ, where µ and σ are themean and standard deviation of the entire cleaned dataset.

The θh is common for all clusters for the following reason:
When thresholds are calculated on a per-cluster basis, the RSR
algorithm (section II), cluster-specific θh leads to significant
impacts on µ and σ, causing the threshold to shift either tothe far right or far left of the distribution. High water usage
clusters can cause the threshold to shift to the far right, leading
to the categorization of unsafe events as benign. Conversely,
low water usage clusters can cause the threshold to shift to theleft, labeling benign water usage as unsafe events. To addressthese challenges, the entire dataset is used to determine thethreshold, hence introducing noise and ensuring the thresholdvalue, θh , is more robust in dealing with these two extremes.
We obtain the value of g = 30 while calculating θh , whichtranslated into having only 0.11% of the data outside of thebound. From this, we obtained θh = 740.7. We validated θhwith real-world statistics [9] of water usage per hour and we
observe that θh is higher than the expected water usage perhour. Thus, if the water usage in an hour t from house i isgreater than 740.6529 liters, we consider it an unsafe event.

Fig. 6: Ground Truth vs Z(t) for test data
Illustration: Figure 6 shows a comparison of predicted
events Z(t) (lower plot) with the ground truth ζ(t) (upperplot) where levels 2 and 3 correspond to events that lead topositive and negative outliers respectively. Thus, using Z(t)helps us understand what kind of event might have happened.
G. Performance Evaluation Metrics
Apart from false alarm and missed detection counts andrates per cluster, we report a less used but NIST recommended performance metric for cyber-physical systems, called
Expected
Time Between False Alarms, defined as: E(Tf a ) =PnF A
i=1 TEnF A , where nF A is the total number of false alarms,
and TE is the duration between two consecutive false alarms.
This metric provides valuable insights into the occurrence andspacing between false alarms [13].
IV. E XPERIMENTAL E VALUATION
A. Experimental Setup
We used residential household water usage data collectedby Hellenic Data Service [3]. The dataset contains water
usage from 1,099 houses from Jan 2015 - May 2017. We
formed a training set of 521 days (Jan 2015 - June 2016),
a cross-validation set of 174 days (May 2016 - Nov 2016)
for hyperparameter selection, and a test set of 174 days (Nov.

2016 - May 2017) for performance evaluation, all containing
uncleaned raw data.
As per Section II-B, we used the clustering approach todivide the houses into 11 clusters. For each cluster, we
identified optimal hyperparameters using the cross-validation
set with Algorithm 4. We fit the optimal hyperparameters
into our model with the training data to extract optimal the
RU C ck (t) set per cluster, which is then fed as input tock
ck
Algorithm 3 to obtain optimal thresholds τmaxand τminfor each cluster ck . The optimal hyperparameters and thethresholds are then applied to the test set. We calculated theinvariant ADck (tcϕ ) for each cluster by using the optimal p+
and p− for that cluster (Sec III-A). Then we use the safe
margin gotten during the training using the optimal ϵ andcompare the invariant ADck (tcϕ ) against the safe margins foreach cluster to obtain the RU C(tc ). We compare RU C(tc )against the robust thresholds using Eqn. 12 to represent thedecision taken by our method on unsafe events. For the sakeof visualization, we use the Eqn. 13 to compare the predictedunsafe events against the ground truth events.
B. Overall Performance with Learnt Hyperparameters
To assess the performance of our anomaly detection approach, we applied the learned p− , p+ , ϵ, τmax , and τmin to eachof the 11 clusters. The results areP
presented in table I, where
ck represents the cluster name,
ζ(t) is the total numberofgroundtruth
unsafeeventsgottenfrom the raw dataset,PZ(t) is the total number of unsafe events our model wasable to predict, the missed detection and false alarms aredenoted as MD and FA respectively.
TABLE I: Experimental Results for the 11 Clusters
ck
1
2
3
4
5
6
7
8
9
10
11
Total

P

ζ(t)
43
45
72
45
37
83
37
213
45
51
43
704

P

Z(t)
38
35
70
43
35
72
32
213
40
46
41
665

MD
5
0
2
2
2
11
5
0
5
5
2
39

FA
6
5
8
7
2
6
3
6
4
1
4
52

E(Tf a )(Hours)
443.0
693.4
320.75231.571307.0504.83824.67543.67773.254175.0713.5

From Table I, we observe that over a testing time of 4,175
hours across all clusters- 94.5% of the unsafe events weredetected, while the overall false alarm rate was 1.4% andthe missed detection rate was 5.5%.
P Table I also
P provides thecluster-wise results for MD, FA,
Z(t) andζ(t).
Furthermore, we give the average expected time between
consecutive false alarms. The worst E(Tf a ) is 231.57 hours
for cluster 4, implying a FA once every 9 days, while cluster
10 has the highest E(Tf a ) with a FA expected only once
every 174 days. The average E(Tf a ) across all clusters in
957.3 hours. In summary, our anomaly detection approach
exhibits accuracy in the identification of unsafe events in SWM
infrastructure based on learned benign behavioral models,
while maintaining a low incidence of FA. While the false
alarm rate is low, more research is needed to further reduce
the frequency of false alarms.

To further understand the significant improvement provided
by Algorithm 3, our results are 15 times better in terms of false
alarm rate and 4.5 better in missed detection rate compared tolinear regression used with our proposed invariant.
V. C ONCLUSION(a) ROC curve(b) Mean E(Tf a )Fig. 7: Evaluating Sensitivity of Performance
C. Performance Sensitivity Analysis
Here we show sensitivity analysis using ϵ as a free parameter, examining its impact on true positive rate (TPR), falsepositive rate (FPR), and E(Tf a ).
Figure 7a presents the average Receiver Operating Characteristic (ROC) curve across the 11 clusters, serving as a robust
indicator of our framework’s performance. Notably, at a 98%
TPR, the FPR remains impressively low, standing at just 0.037.
This low FPR is crucial to Cyber-Physical Systems (CPS).
This is because of the infrequent occurrences of events, thecost of false alarms is high. Specifically, Figure 7b provides
insight into the expected time gap between two successive
false alarms across varying values of ϵ.
D. Comparison with Existing Works
Now, we compare our proposed approach with four existing
methods for anomaly detection: (i) CUSUM control [14], [15];(ii) α-Winsorization [16]; (iii) the closest theoretically similar
work proposed by us [12] that uses the idea of correlation
clustering and extracts a Pythagorean mean ratio invariant(from time series) and Cauchy Lorentz loss function for arobust threshold learning algorithm under noisy training data;(iv) our invariant AD(tϕ ) when combined with simple linear
regression (i.e., square loss) to show the design benefit of themodified Hampel loss function.

Fig. 8: Comparison with Previous Methods
Fig. 8 demonstrates that our proposed approach (marked
in blue bars) outperforms other methods in terms of anomaly
detection rate, missed detection, and false alarm rates.
Specifically, in terms of detection rate, our proposed method
is two times better than [12] and approximately three times
better than CUSUM and Winsorized Statistics. In terms offalse alarm rate, our method exhibits only 1.4% while theframework in [12] has 8%. Thus, the proposed method is morethan four times better compared to the next best approach.

We propose a framework for real-time detection of unsafe
events in SWM infrastructure using a time-series approach.
First, a clustering algorithm groups SWMs with similar water
usage and then utilizes the absolute difference between thenegative and positive order of the generalized means to create
an invariant that takes higher values under unsafe events andlower values under benign conditions. We proposed a modified
application of Hampel three three-part loss function as a noiseresilient technique to learn the bounds of normal behavior
in the presence of unsafe events in the training set. The
results show that the proposed framework allows decentralized
detection of unsafe events with high accuracy and low falsepositive rates, and outperforms various well-known methods.
Acknowledgements: The work is supported by NSF grants
CNS-2030611, CNS-2030624, and DGE-1433659.
R EFERENCES[1] S. Bhattacharjee and S. K. Das, “Detection and Forensics against Stealthy Data
Falsification in Smart Metering Infrastructure,” IEEE Transactions on Dependable
and Secure Computing, 18(1): 356–371, Jan 2021.
[2] S.Bhattacharjee, P. Madhavarapu, S. Silvestri, S.K. Das, “Attack Context Embedded Data Driven Trust Diagnostics in Smart Metering Infrastructure”, ACM Trans.
on Security and Privacy, 2021.
[3] Helix, “Smart Water Meter Consumption Time Series - Datasets - HELIX,”
Hellenic Data Service, 2020. https://data.hellenicdataservice.gr/dataset/78776f38a58b-4a2a-a8f9-85b964fe5c95 (accessed Apr. 03, 2022).
[4] Daiad, “Trials Evaluation and Social Experiment Results,” European
Commission’s 7th Framework Programme, 2017. daiad.eu/wp-contentuploads/
201711D7.3 Trials Evaluation v1.0.pdf. (accessed Apr. 03, 2022).
[5] J. Islam et al., ”Anomaly based Incident Detection in Large Scale Smart Transportation Systems,” ACM/IEEE 13th International Conference on Cyber-Physical
Systems (ICCPS), Milano, Italy, 2022, pp. 215-224.
[6] A. Kraskov, H. Stögbauer, and P. Grassberger, “Estimating Mutual Information,”
Phys Rev E Stat Phys Plasmas Fluids Relat Interdiscip Topics, 69(6), p. 16, 2004.
[7] D. F. Andrews, P. J. Bickel, F. R. Hampel, P. J. Huber, W. H. Rogers, and J. W.
Tukey, Robust Estimates of Location: Survey and Advances, Princeton University
Press, Princeton, NJ, USA, 1972.
[8] N. X. Truong and T. H. Thiep, “Selection of Hydrological and Hydraulic Models
Applied in Urban Drainage,” South Asian Research Journal of Engineering andTechnology, vol. 3, no. 6, 2021.
[9] A. Gutierrez-Escolar, A. Castillo-Martinez, J. M. Gomez-Pulido, J. M. GutierrezMartinez, and E. Garcia-Lopez, “A New System for Households in Spain to
Evaluate and Reduce Their Water Consumption,” Water, 6(1): 181–195, Jan 2014.
[10] B. J. Frey and D. Dueck, “Clustering by Passing Messages between Data Points,”
Science, 315(5814): 972–976, Feb 2007.
[11] R. Mushtaq, “Augmented Dickey Fuller Test,” SSRN Electronic Journal, Aug 2011.
[12] P. Roy, S. Bhattacharjee, S. Abedzadeh, and S. K. Das, “Noise Resilient Learning
for Attack Detection in Smart Grid PMU Infrastructure” IEEE Trans Dependable
and Secure Computing, to appear, 2024.
[13] D. I. Urbina et al., “Limiting the Impact of Stealthy Attacks on Industrial
Control Systems,” Proceedings of ACM SIGSAC Conference on Computer andCommunications Security, 2016, doi: 10.1145/2976749.
[14] “6.3.2.3. Cusum Control Charts.” Accessed: Jan. 31, 2024. [Online]. Available:
https://www.itl.nist.gov/div898/handbook/pmc/section3/pmc323.htm[15] V. R. Palleti, V. K. Mishra, C. M. Ahmed, and A. Mathur, “Can Replay Attacks
Designed to Steal Water from Water Distribution Systems Remain Undetected?,”ACM Transactions on Cyber-Physical Systems, 5(1): 1–19, 2020.
[16] R. Wilcox, “Trimming and Winsorization,” Encyclopedia of Biostatistics, Feb.
2005.
[17] P. Schober and L. A. Schwarte, “Correlation coefficients: Appropriate use andinterpretation,” Anesth Analg, vol. 126, no. 5, pp. 1763–1768, May 2018.

VI. A PPENDIX
A. Data Cleaning for Clustering
Let us discuss the data cleaning process required for forming
the cluster with the RSR algorithm in section II-B5.
To clean a household’s negative values at a time t, weanalyzed adjacent values, selected the positive value closest
in magnitude to the absolute value of the negative number,and summed them up. We consider t′ of the selected value as
the correctivet′ . If the sum is negative, we attribute it to theoriginal hour and set the correctivet′ to zero, else, the originalhour becomes zero and the correctivet′ retains the sum. This
also was able to handle some large positive numbers. The
rationale behind this approach lies in the assumption that theoccurrence of negative values can be addressed by summingthem up with adjacent positive values. This assumption is
also backed up by further study into the data where it wasobserved that some instances of large negative readings weremirrored (or corrected) by approximately equal large instancesof positive values, which we attributed to potential meter
malfunctions. Summing these values refines the dataset toaccurately depict the actual water usage.
If two or more positive values exist within an hour, we
obtained and assigned their sum to that specific hour. Thisis because SWM could potentially capture and send readings
more frequently to the data concentrator than what is expected.
For missing data over multiple hours, we identified thepotential cause as network malfunction. We addressed this by
averaging the next available value over the past tm +1 missing
values and recorded this average as the reading for the current
hour and the tm consecutive missing hours. This was able
to handle some large positive numbers. This is because we
assume that each hour preceding the spike contains some datarather than being entirely missing.
For applications where location or geographical data ispresent, the clustering could be made more accurate.
B. Details on Augmented Dickey Fuller (ADF) Test
In this section, we present the step-by-step processes wefollowed to test for the goodness of the RSR algorithm informing clusters presented in section II that is crucial toachieving invariance in metrics.
Time-series: [2] proposed an invariant that works well on
cleaned data based on the absolute difference between the
Harmonic Mean (HM) and the Arithmetic mean (AM) at time
tϕ . We use this to obtain the time series for the ADF, which
is denoted by ad(tϕ )ck such thatad(tϕ )ck = |HM (tϕ )ck − AM (tϕ )ck |and variance, remain constant over time. A more negative (i.e.,higher in magnitude) test statistic indicates stronger evidence
for stationarity [11]. This is defined as:∆yt = a + βt + γyt−1 + δ1 ∆yt−1 + δ2 ∆yt−2 +δ3 ∆yt−3 + . . . + δq−1 ∆yt−(q−1) + et

(16)

Where ∆yt is the differenced time series; a is a constant; β
is the time trend coefficient; γ is the coefficient of the lagged
level of the time series; δ1 , . . . , δq−1 are the lagged differences
coefficients of the time series; and et is the error term.
In the framework of the ADF test, the null hypothesis (H0 )and alternative hypothesis (H1 ) is formulated as follows:
H0 : The time series data possesses a unit root
H1 : The time series data is stationary
The presence of a unit root in a time series data analysis indicates a persistent stochastic trend, making it non-stationary,
its absence however signifies stationarity. We embraced thisconcept in our study. To use this test, once a cluster isformed, we aggregate the water usage of each house basedon equation 15 and then subject ad(tϕ )ck to the ADF test
(equation 16) to assess the stationarity which would indicate
the efficacy of the RSR algorithm. We provided our hypothesis
based on this and explained it with two clauses.
H0 : The time series data possesses a unit root
H1 : RSR is effective; time series data is stationary
Clause A (H1 ): The ADF test can determine if a cluster
is stationary given the expected cyclostationary behavior ofad(tϕ )ck for corresponding hours of each day, thus, indicating
the RSR’s effectiveness of clustering houses with similar water
usage. The concept is depicted in fig. 9a. The ad(tϕ )ck for 3
clusters is shown for 6 days. We observe the stationarity of
the ad(tϕ )ck indicating the efficacy of the RSR in clustering
houses with similar water usage patterns, thus, the visual
similarity for each corresponding hours of each day, i.e., thead(tϕ )ck at t = 00 : 00 on day 1 is similar to the ad(tϕ )ck at
t = 00 : 00 from day 2 to day 6.
Clause B (H0 ): Conversely, if the algorithm is ineffective, itwill group houses with dissimilar water usage patterns leading
to a unit root, thus, diverging the ad(tϕ )ck significantly,
leading to non-stationarity when subjected to the ADF test.

(15)

where Eq. 15 is calculated per cluster ck . Also, the potencyof this combination was confirmed by preliminary results. Thead(tϕ )ck is then fed as the time series data into the ADF toassess the stability.
Augmented Dickey-Fuller Test (ADF): The ADF test is astatistical method used to determine if the time series is stationary, indicating that its statistical properties, such as mean(a) Stationary ad(tϕ )ck with RSR(b) Test statistics for all clusters

Fig. 9: RSR results for wcut = 0.061 and N = 80

+C. Results of Test Statistics for Each Cluster
After obtaining the optimal parameter for wcut and N
in section II-C, we inserted these parameters into the RSR
algorithm to calculate ad(tϕ )ck . We present a comprehensive
breakdown of the ∆ykck per cluster in Fig. 9b. We observethat the RSR initially clustered strongly correlated houses;
however, as the iterations progressed, these correlations weakened, resulting in higher test statistics. Despite this, all formedclusters passed the stationarity test. The stationarity test isdetermined by subjecting the test statistics to a critical value,
and if the test statistic is greater than the critical value, thenull hypothesis cannot be rejected, suggesting that the timeseries is non-stationary.
D. Comparison of the RSR Algorithm
We compared the RSR algorithm with time series K-means(TSKM), Agglomerative clustering (Agg), and Density-BasedSpatial Clustering of Applications with Noise (DBSCAN). We
provide the results in Table II where we observe that DBScan
and RSR perform better than TSKM and Agg in terms of thepercentage of clusters passing the stationarity test. We obtained
lower test statistics of −17.48 for the proposed RSR whencompared to DBSCAN which is −11.03, indicating that RSR
gives more stable invariants compared to DBSCAN even withthe same percentage of stationary clusters (those that pass thestationarity test).
TABLE II: Comparison with Other Clustering Algorithms
Algorithm
Number of Clusters Formed
Percentage of stationary cluster (%)

TSKM
8
62.5

Agg
18
77.8

DBSCAN
1
100

RSR
11
100

E. Deviation Reasoning
Here, we give a theoretical explanation of why our proposedAD(tϕ ) in section III deviates under unsafe events.
In mathematics, a Schur − convex function, alternatively
referred to as an S-convex or isotonic, is a type of function
denoted as f : Rd → R. This mathematical concept has afundamental property: for any pair of vectors, namely, x andy belonging to the d-dimensional real space Rd , if it holds
true that x is majorized by y or x is more spread out than y,then this function satisfies the condition f (x) ≤ f (y). This
property implies that the function assigns smaller values tovectors that are ”dominated” by other vectors.
In contrast, f is Schur − concave, which means that any
pair of vectors such that x is majorized by y, then this function
satisfies the condition f (x) ≥ f (y). Let’s start interpretingthe behavior of the generalized mean,GM p (x1 , x2 , . . . , xn ) in
Eqn. 4, in terms of Schur concavity and Schur convexity:−
If p < 0: the generalized mean GM p ≡ GM p is Schurconcave. In terms of the power mean formula, since smaller
values are given more weight, an unsafe event causing a large
decrease (e.g., negative water flow) in one or more numbers in−the series (x1 , x2 , . . . , xn ), causes the resulting GM p value
to decrease sharply, and the magnitude of the decrease iscontrolled by the value of p− .

If p ≥ 1: the generalized mean GM p ≡ GM p is Schurconvex. In terms of the power mean formula, since smaller
values are given more weight, an unsafe event causing a large
increase (e.g., leakage, network malfunction) in one or morenumbers in the series (x1 , x2 , . . . , xn ), causes the resulting+
GM p value to increase sharply, and the magnitude of theincrease is controlled by the value of p+ .
Since our invariant definition contains a combination of both+−
GM p and GM p such that−+AD(tϕ ) = |GM p (tϕ ) − GM p (tϕ )|regardless of the nature of the unsafe event. Given we are+−taking the absolute difference between GM p and GM p , ourinvariant will increase under unsafe events following properties
listed in this work, as long as p ̸= −1

Fig. 10: The Absolute Difference (AD(tϕ ))As illustrated in Fig. 10, there are discernible spikes withina specific time window. These spikes correspond to exceptionally high values of AD(tϕ ) in Equation 5. This phenomenon−
arises when there is a substantial separation between GM pp+p−p+
and GM . In such cases, either GM
or GM
assumesa significantly elevated value.
−
When GM p exhibits a very large value, it implies thatwithin the given time window, numerous instances of water
usage data fall within the range of very small to moderatevalues, while larger values are relatively less prominent. Con+
versely, when GM p gets a high value, it signifies that within
the same time window, very large to moderate water usage
values are given priority over smaller values.
Thus, the detection of highly significant small-scale water
usage or water drainage events can be accomplished throughthe analysis of GM p , where p < 1. Similarly, the identificationof substantial large-scale water usage events can be achievedby scrutinizing GM p , where p > 1.