Milad Leyli-Abadi
email: milad.leyli-abadi@ifsttar.fr

Allou Samé
email: allou.same@ifsttar.fr

Latifa Oukhellou
email: latifa.oukhellou@ifsttar.fr

Nicolas Cheifetz
email: nicolas.cheifetz@veolia.com

Pierre Mandel
email: pierre.mandel@veolia.com

Milad Leyli-Abadi;

Cédric Féliers
email: cedric.feliers@veolia.com

Olivier Chesneau

Mixture of Joint Nonhomogeneous Markov Chains to Cluster and Model Water Consumption Behavior Sequences

Keywords: CCS Concepts:, Mathematics of computing → Markov networks, Time series analysis, Cluster analysis, • Computing methodologies → Unsupervised learning, Mixture modeling non-homogeneous Markov models, categorical time series, clustering, forecasting, water consumption behavior

The emergence of smart meters has fostered the collection of massive data that support a better understanding of consumer behaviors and better management of water resources and networks. The main focus of this paper is to analyze consumption behavior over time; thus, we first identify the main weekly consumption patterns. This approach allows each meter to be represented by a categorical series, where each category corresponds to a weekly consumption behavior. By considering the resulting consumption behavior sequences, we propose a new methodology based on a mixture of nonhomogeneous Markov models to cluster these categorical time series. Using this method, the meters are described by the Markovian dynamics of their cluster. The latent variable that controls cluster membership is estimated alongside the parameters of the Markov model using a novel classification expectation maximization (CEM) algorithm. A specific entropy measure is formulated to evaluate the quality of the estimated partition by considering the joint Markovian dynamics. The proposed clustering model can also be used to predict future consumption behaviors within each cluster. Numerical experiments using real water consumption data provided by a water utility in France and gathered over nineteen months are conducted to evaluate the performance of the proposed approach in terms of both clustering and prediction. The results demonstrate the effectiveness of the proposed method.

INTRODUCTION

The issue of effective water resource management becomes increasingly important as populations grow and the climate changes. The problem is more acute for dense urban areas that involve complex water network architectures. One significant performance leverage about management and distribution of drinking water is meter readings taken over a consumption period. By analyzing the collected information, water utilities in collaboration with national authorities can establish different types of strategies, e.g., price strategies [START_REF] Charles | The impact of price on residential water demand and its relation to system design and price structure[END_REF][START_REF] Randolph | Attitudes to conservation and water consumption[END_REF], drought strategies [START_REF] Douglas S Kenney | Residential water demand management: lessons from Aurora, Colorado1[END_REF][START_REF] Sandra | Entering an era of water scarcity: the challenges ahead[END_REF] and strategies to promote conservation [START_REF] Sheila | Comparing price and nonprice approaches to urban water conservation[END_REF]. The objective of these strategies is to encourage consumers to use water more responsibly. To achieve this goal, the collected consumption data must be acquired with a fine granularity. Recently, the increasing use of smart meters in the framework of smart water grids has provided water utilities with large amounts of telemetered consumption data at very low sampling rates (i.e., from fifteen seconds interval to one hour) compared to the traditional monthly sampling rate usually used for billing purposes. These finer spatial and temporal sampling resolutions make it possible to develop advanced data analytics in this domain. The data with a high longitudinal frequency are used for water end-use disaggregation [START_REF] Beal | A novel mixed method smart metering approach to reconciling differences between perceived and actual residential end use water consumption[END_REF] while the lower frequency data can be utilized for short-term consumption analysis [START_REF] Bougadis | Short-term municipal water demand forecasting[END_REF][START_REF] Gagliardi | A probabilistic short-term water demand forecasting model based on the Markov Chain[END_REF][START_REF] Jain | Short-term water demand forecast modelling at IIT Kanpur using artificial neural networks[END_REF] and other applications, such as detecting water leaks [START_REF] Berglund | Successive linear approximation methods for leak detection in water distribution systems[END_REF][START_REF] Tracy C Britton | Smart metering: enabler for rapid and effective post meter leakage identification and water loss management[END_REF].

One primary issue that has been addressed by several studies in the domain of water data analysis is related to clustering, in which data are reduced to an identified set of clusters to facilitate further analysis and enhance decision-making capabilities. Consumers who share common consumption behaviors are assigned to the same cluster. By partitioning the consumers into different clusters, water companies can monitor and manage the demand more efficiently. For example, one way to distinguish consumers is based on their water consumption level [START_REF] Aksela | Demand estimation with automated meter reading in a distribution network[END_REF][START_REF] Sa Mckenna | Water demand pattern classification from smart meter data[END_REF]. Another possibility is to investigate the consumer consumption behavior over different time periods [START_REF] Labeeuw | Residential electrical load model based on mixture model clustering and Markov models[END_REF]. Then, each cluster will consist of consumers who share the same consumption habits regardless of their water consumption levels.

Various exogenous factors can also influence consumption behaviors, including environmental and water conservation attitudes [START_REF] Willis | Quantifying the influence of environmental and water conservation attitudes on household end use water consumption[END_REF], climatic variables [START_REF] House-Peters | Effects of Urban Spatial Structure, Sociodemographics, and Climate on Residential Water Consumption in Hillsboro, Oregon 1[END_REF][START_REF] Shuang | Forecasting daily urban water demand: a case study of Melbourne[END_REF], socio-economic and demographic information [START_REF] Domene | Urbanisation and water consumption: Influencing factors in the metropolitan region of Barcelona[END_REF], calendar events [START_REF] Zhou | Forecasting operational demand for an urban water supply zone[END_REF], and so on. These factors should also be considered by water utilities when analyzing consumers' behavior. Several related works have studied the influence of such covariates on water consumption using statistical approaches and feature selection methods [START_REF] Cardell-Oliver | Discovering water use activities for smart metering[END_REF][START_REF] Guyon | An introduction to variable and feature selection[END_REF].

Keeping that information in mind, the main focus of this paper is on clustering of categorical sequences, which represent the evolution of consumer behavior over time. Note that in this study, clustering is performed by considering both the consumption behavior dynamics and other available exogenous factors. For this purpose, we propose a clustering approach based on a mixture of nonhomogeneous Markov models that uses contextual factors as model inputs. Aside from clustering, the proposed model can be used to perform prediction tasks. Indeed, future consumption behaviors within the identified clusters can be predicted by using the estimated parameters of each of the nonhomogeneous (input-dependent) Markov models in the mixture. We introduce several graphical tools to facilitate the interpretation of the results.

The rest of this paper is organized as follows. Section II presents the related works. Section III describes the core components of the proposed methodology for clustering and forecasting water consumption behaviors, presents the details of the estimation algorithm, and describes the analyzed real-world dataset. In Section IV, we present the experimental results using graphical and numerical tools in two separate sub-sections: the first subsection is related to the clustering results and the second subsection provides an evaluation of the forecasting performance of the proposed method.

RELATED WORK

Analyzing user behavior using Markov models has been widely studied in application domains such as energy, transport and web traffic analysis. In [START_REF] Cadez | Model-based clustering and visualization of navigation patterns on a web site[END_REF], the authors analyzed the navigation patterns of website users. The proposed approach partitioned users into clusters in which the users shared similar navigation paths through the site. To consider dynamic user behavior, the authors used a mixture of first-order Markov models. Each state of the Markov model corresponded to a page category that a user could visit. The authors showed that some higher order dependencies can also be captured using the mixture mechanism. Finally, the visualization results confirm the appropriateness of using Markov mixture models to cluster categorical time series. In contrast to this model, our proposed approach assumes the model to be input-dependent.

The consumption behavior of energy consumers has been analyzed in several different research papers [START_REF] Kwac | Household energy consumption segmentation using hourly data[END_REF][START_REF] Wang | Clustering of electricity consumption behavior dynamics toward big data applications[END_REF]. In [START_REF] Wang | Clustering of electricity consumption behavior dynamics toward big data applications[END_REF], symbolic aggregate approximation was used in conjunction with the Markov models to reduce the scale of the data and to model the dynamic behavior of consumers using transition matrices. To partition consumers in different groups, a density-based clustering technique (CFSFDP) was used that takes as input the dissimilarity matrix obtained by the Kullback-Liebler distance between each pair of transition matrices. The consumption behavior of an adjacent time period is analyzed for different groups throughout an entropy measure. However, this work did not study the influence of exogenous factors on consumption behaviors. A methodology based on an encoding system with a preprocessed load shape dictionary was proposed in [START_REF] Kwac | Household energy consumption segmentation using hourly data[END_REF]. The load shape information of energy consumers was used to classify households with regard to extracted features such as an "entropy of shape" code which measures consumption variability. This study derived five sample programs and policy-relevant energy lifestyle segmentation strategies.

To analyze the behavior of public transport users and predict their future trajectories, the authors of [START_REF] Yang | Unraveling traveler mobility patterns and predicting user behavior in the Shenzhen metro system[END_REF] proposed clustering users into three categories (regular, variable and irregular) based on a spatiotemporal entropy measure and the k-means algorithm. They investigated Markov models and their hidden extensions to predict individuals' future movements within different predefined time periods throughout a day (morning, afternoon, evening, and night). The latency time (the amount of time spent at a place after exiting a metro station) and users' current locations were considered as inputs to the model. The results of a comparison indicated that the hidden Markov models were able to achieve better prediction accuracy. This model differs from ours because we used a model-based algorithm for clustering that considers the temporal dynamics of consumption behaviors. In this way, context-based clusters can be identified rather than performing grouping based on a regularity measure.

As discussed above, few papers have studied user behavior while considering the joint Markovian dynamics. In most of the developed methods consumption patterns are clustered regardless of exogenous factors. Furthermore, to capture the consumption behavior dynamics, clustering and modeling are typically performed in two separate stages. To fill these gaps, this paper presents an integrated methodology based on a mixture of nonhomogeneous Markov models that merges clustering and modeling into a single uniform platform. In this platform, each transition between an adjacent pair of consumption behaviors (called states) is assumed to be dependent on some exogenous variables.

The contributions of this study can be summarized as follows:

• The categorical water consumption series extracted from a real-world dataset are classified using a clustering approach based on a mixture of nonhomogeneous Markov models; • The consumption behavior dynamics related to each identified cluster are visualized using appropriate graphical representations and a quantitative metric; • Future consumption behaviors are predicted for each cluster.

METHODOLOGY

The proposed methodology for clustering and forecasting the consumption behaviors consists of the different phases shown in Figure 1. The consumption profiles are normalized in a preprocessing step. The typical consumption patterns are extracted from the normalized consumption profiles, and the sequences of these patterns are used to construct a new categorical database. In this way, each meter in this categorical dataset is represented by the evolution of its users' consumption behaviors (categorical states). Next, the consumption behavior dynamics are modeled using a mixture of nonhomogeneous Markov models, leading to the estimation of a latent variable that encodes the class membership. Finally, future consumption behaviors are predicted within each cluster. These phases are described in the following subsections.

Data and preprocessing

The dataset analyzed in this paper stems from 2,000 meters located in suburban areas of Paris, France. The dataset covers a 19-month consumption period: from January 26th, 2015 to July 24th, 2016. Data collection is performed hourly using automatic meter reading (AMR) [START_REF] Kimbel | Automatic meter reading data communication system[END_REF], and the meters mainly reflect both individual and collective housing. The data collected from the smart meters exhibit a significant number of missing values due to meter malfunctions or erroneous readings. A pretreatment stage is applied to address these missing consumption values. However, because this paper focuses on clustering and prediction of water consumption time series, we do not discuss the specific details of the preprocessing stage.

Depending on the desired forecasting time horizon (short-term, medium-term or long-term), either daily or weekly profiles can be investigated. In this paper, we investigate weekly profiles, since the aim is to be able to forecast the medium-term consumption behaviors (weekly behaviors); however, the proposed approach could also be performed by considering daily profiles. In a similar work [START_REF] Leyli-Abadi | Predictive Classification of Water Consumption Time Series using Non-homogeneous Markov Models[END_REF], we have used the daily profiles to forecast the short-term consumption behaviors. To analyze the evolution of consumption behaviors of consumers over time, the weekly demand profiles are log-normalized. Subsequently, all the profiles are centered and reduced. The normalization has two main purposes: it scales down the consumption data for all the consumers so they present less variability, and it causes them to approach a normal distribution. These normalized data are used as inputs for the next stage.

Apart from the consumption data, this study also considers different exogenous variables such as climatic factors and calendar events. The climatic variables are temperature and precipitation level, while the calendar variables are public and school holidays.

Data discretization

To be able to analyze the massive data collected using smart meters, one solution is to reduce the data size. To accomplish this, the normalized weekly consumption curves are discretized using a functional clustering algorithm [START_REF] Samé | Décomposition et classification de données fonctionnelles pour l'analyse de la consommation d'eau[END_REF]. The Bayesian information criterion (BIC) [START_REF] Schwarz | Estimating the dimension of a model[END_REF] is used to select the number of states. As the result of clustering, we identified eight main weekly consumption behaviors or states (K = 8), as shown in Figure 2, where each cluster center (the solid lines) is visualized with its standard deviation (brighter colors) on the left side, and the corresponding superimposed daily patterns are shown on the right. Looking at the superimposed daily patterns, the differences between consumption behaviors on weekdays (solid lines) and weekends (dashed lines) can clearly be observed. Note that for states 1-5, the morning peaks occur later in time on weekends. State 6 is characterized by a symmetrical consumption behavior on Monday and Friday with respect to Saturday and Sunday. A single peak during weekdays followed by a more consistent consumption during weekends characterizes the state 7, which can be associated with consumption behavior in a commercial activity zone. Finally, state 8 shows higher morning peaks during weekends and probably represents a consumption behavior related to holidays.

The resulted categorical dataset encoding the weekly consumption curves of each residence is shown in Figure 3. In this figure, each of the identified states (main weekly consumption behaviors) is shown by the corresponding color for each sequence. (s it ) i=1, ...,n;t =1, ...,T : dataset s it state associated to the meter i at day t (s it ∈ {1, . . . , 5}) s i (s it ) t =1, ...,T : consumption sequence for a meter i s t (s it ) i=1, ...,n : data of day t w i label associated to the time series s i w (w i ) i=1, ...,n : set of labels associated to meters s i e (e it ) i=1, ...,n;t =1, ...,T : set of input variables e it input associated to the meter i at day t ϕ set of all the parameters of the model ϕ д set of the parameters for the meters of cluster д K number of states G number of time series clusters

Modeling

Model definition.

Using a set of categorical consumption time series s = {s 1 , . . . , s i , . . . , s n }, where s i is measured over T different weeks (s i = {s i1 , . . . , s iT }), the proposed model assumes that each series s i is issued from one of G groups. A clustering random variable w = (w 1 , . . . , w n ) is associated with s where each realization of this random variable w i takes a value in д ∈ {1, . . . , G}.

In each group д, the time series are modeled by a first order nonhomogeneous Markov model from the input data or entries e [START_REF] Leyli-Abadi | Predictive Classification of Water Consumption Time Series using Non-homogeneous Markov Models[END_REF], [START_REF] Bengio | Input-output HMMs for sequence processing[END_REF], [START_REF] Bengio | Markovian models for sequential data[END_REF]. This model has the following property:

P(s it | s i1 , . . . , s it -1 , e i1 , . . . , e it ) = P(s it | s it -1 , e it ). (1)
The current state s it depends only on the previous state s it -1 and the input vector e it . A corresponding graphical representation of the proposed mixture model is shown in Figure 4. In this setting, we have considered that each state s it corresponds to one of the eight identified main weekly consumption behaviors. This model is specified by the following conditional mixture density:

f (s i |e i ; ϕ) = G д=1 p д f д (s i |e i ; ϕ д ), (2)
where p д is the proportion of the cluster д that satisfies G д=1 p д = 1, f д (.) is the distribution associated to the group д, and ϕ д is its set of parameters. The distribution of a time series s i from a group д is given by

f д (s i |e i ; ϕ д ) = P(s i1 | w i = д, e i1 ; ϕ д ) T t =2 P(s it | w i = д, s it -1 , e it ; ϕ д ), (3) with P
(s i1 = ℓ | w i = д, e i1 ) = π дℓ (e i1 ; α д ) = exp (α ⊤ дℓ e i1 ) K ℓ ′ =1 exp (α ⊤ ℓ ′ e i1 ) , (4)
P(s it = k | s it -1 = ℓ, w i = д, e it ) = π дℓk (e it ; β дℓ ) = exp (β ⊤ дℓk e it ) K k ′ =1 exp (β ⊤ дℓk ′ e it ) , (5)
where the π дℓ s designate the initial input-dependent probabilities whose parameter vector is α д = (α дℓ ) ℓ=1, ...,K , and the π дℓk s designate the input-dependent transition probabilities of partition д whose set of parameters are β дℓ = (β дℓk ) ℓ=1, ...,K . In the following, we call the proposed mixture of joint nonhomogeneous Markov models "MixJNMM. "

= n i=1 log G д=1 p д P(s i1 | w i = д, e i1 ; ϕ д ) T t =2 P(s it | s it -1 , w i = д, e it ; ϕ д ) , (6)
where ϕ = (p д , α д , β д ) is the complete set of model parameters, where β д = (β дℓk ) ℓk . The values of P(s i1 | w i = д, e i1 ; ϕ д ) and P(s it | s it -1 , w i = д, e it ; ϕ д ) are given by Equations ( 4) and ( 5), respectively. This criterion cannot be maximized directly. However, one can use the expectationmaximization (EM) algorithm [START_REF] Arthur P Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] to maximize it. Considering the large global data size, we adopted the classification EM (CEM) algorithm [START_REF] Celeux | A classification EM algorithm for clustering and two stochastic versions[END_REF] for this study because it converges faster than the EM algorithm. Using CEM, both the latent random variable w and the other parameters of the model are estimated by maximizing the following complete data likelihood: P(s, w | e; ϕ) = P(w) P(s | w, e; ϕ)

= n i=1 P(w i ) P(s i | w i , e i ; ϕ w i ). (7)
Developing 7, we get the following complete log-likelihood:

L c (ϕ, w) = log P(s, w | e; ϕ) (8)
= G д=1 n i=1 w iд log p д + n i=1 K ℓ=1 w iд s i1ℓ log π дℓ (e i1 ; α д ) + n i=1 T t =2 K k =1 K ℓ=1 w iд s it k s i(t -1)ℓ log π дℓk (e it ; β дℓ ) ,
for which we have used the following indicator variables to specify the cluster memberships and designate the transitions of the Markov model:

• w iд = 1 if w i = д and w iд = 0 otherwise; • s i1ℓ = 1 if s i1 = ℓ and s i1ℓ = 0 otherwise; • s it k = 1 if s it = k and s it k = 0 otherwise; • s it k s i(t -1)ℓ = 1 if s it = k and s i(t -1) = ℓ simultaneously and s it k s i(t -1)ℓ = 0 otherwise.
To promote the robustness of the parameter estimation, we add a regularization term to the likelihood function of Equation [START_REF] Tracy C Britton | Smart metering: enabler for rapid and effective post meter leakage identification and water loss management[END_REF], which leads to the following optimization problem:

arg max (ϕ,w) (L c (ϕ, w) - ξ 2 ∥ϕ ∥ 2 2 ), (9)
where ∥.∥ 2 is the L 2 norm and ξ is a hyperparameter that controls the importance of the regularization term.

In the following section, we introduce the proposed CEM algorithm to estimate the parameters of the mixture model.

Proposed CEM algorithm. The CEM algorithm [START_REF] Celeux | A classification EM algorithm for clustering and two stochastic versions[END_REF] used here is a variant of EM that incorporates a classification step between the E-and M-steps of the EM algorithm. Starting from an initial partition w (0) , the q-th iteration of CEM for the proposed mixture approach is defined as follows:

-E-step. For i = 1, . . . , n and д = 1, . . . , G, compute the current posterior probabilities τ (q) д (s i ) that s i belongs to cluster д as follows:

τ (q) д (s i ) = p (q) д f д (s i ; ϕ (q) д ) G д ′ =1 p (q) д ′ f д ′ (s i ; ϕ (q) д ′ ) , (10)
given the current parameter estimates ϕ (q) .

-C-step. The partition w (q+1) is defined by assigning each observation s i to the cluster that provides the maximum current posterior probability τ

(q) д (s i ), (1 ≤ д ≤ G) as follows: w (q+1) i = arg max 1≤д ≤G τ (q) д (s i ). (11)
-M-step. Compute the maximum likelihood estimates (p

q+1 д , ϕ (q+1) д
), which leads to

p (q+1) д = n i=1 w (q) iд n , ∀ д = {1, . . . , G}. (12)
Furthermore, the parameters of the model (α

(q+1) д , β (q+1)
дℓ ) can be computed by maximizing the quantity:

L 1 (α д ) + K ℓ=1 L 2, ℓ (β дℓ ), (13)
where

L 1 (α (q+1) д ) = K ℓ=1 n i=1 w iд s i1ℓ log π дℓ (e i1 ; α (q+1) д ), (14)
and

L 2, ℓ (β (q+1) дℓ ) = T t =2 K k =1 n i=1 w iд s it -1ℓ s it k log π дℓk (e it ; β (q+1) дℓ ). (15)
Consequently, the problem ( 9) can be solved by the following G × (K + 1) separate maximization problems: arg max

β дℓ L 2, ℓ (β дℓ ) - ξ 2 ∥β дℓ ∥ 2 2 , ∀ д, ℓ arg max α д L 1 (α д ) - ξ 2 ∥α д ∥ 2 2 , ∀ д = 1, . . . , G. (16)
In this paper, we opted for the Newton algorithm [START_REF] Roos | Theory and algorithms for linear optimization -an interior point approach[END_REF], known in this situation as the iteratively reweighted least squares (IRLS) algorithm [START_REF] Paul | Robust regression using iteratively reweighted least-squares[END_REF]. It converges to the optimal solution at a sufficient speed.

Model selection. In the proposed method, the best model is defined as the model with the optimal value of the number of time series clusters G. To select the best model, we use the following integrated classification likelihood (ICL) criterion [START_REF] Biernacki | Assessing a mixture model for clustering with the integrated completed likelihood[END_REF] which is essentially the ordinary BIC [START_REF] Schwarz | Estimating the dimension of a model[END_REF] penalized by the entropyi,д τ д (s i ) log τ д (s i ):

ICL(G) = L c ( φ) - ϑ (G) 2 log(n × T ), ( 17
)
where φ is the maximum likelihood estimate of the parameter vector ϕ, and ϑ (G) is the number of free parameters in the model, which is given by

ϑ (G) = G × (K -1)m α д + K × (K -1)m (β дℓ ) ℓ=1, . . ., K + G -1 (p д ) д=1, . . .,G , ( 18
)
where G is the number of clusters, K is the number of states of the variables s it , and m is the dimension of the input vector e it . The ICL values are computed for different numbers of parameters G. Finally, the model with the minimum ICL value is selected.

3.3.3 Forecasting. After the parameters have been estimated from the training sequence (s 1 , . . . , s T ), the following one-step-ahead forecast permits future state prediction within each group д:

s iT +1 = arg max k P(s iT +1 = k | s iT , w i = д, e iT +1 ) = arg max k π дs iT k (e iT +1 ; β s iT ). (19)
The pseudocode in 1 summarizes the proposed MixJNMM approach.

д , β (0) д ), w (0) ; fix a threshold ϵ > 0; set q ← 0.

while

L (q+1) c (ϕ) -L (q) c (ϕ) > ϵ do E-Step: forall i, д do compute τ (q)
д (s i ) using Eq. ( 10) end C-Step:

forall i do compute w (q+1) i
using Eq. ( 11) end M-Step:

for д = 1, . . . , G do compute p (q+1) д
using Eq. ( 12) for ℓ = 1, . . . , K do compute the parameters ϕ

(q+1) д = (α (q+1) д , β (q+1)
дℓ ) by maximizing the likelihood (13)

end q ← q + 1 end ϕ = ( p (q) д , α (q) д , β (q)
д ) and w = w (q) . end ACM Trans. Intell. Syst. Technol., Vol. 1, No. 1, Article 1. Publication date: January 2019.

EXPERIMENTAL RESULTS

We applied the proposed methodology to the previously discussed categorical dataset with the objectives of learning the heterogeneous structure of the consumption behaviors time series and predicting future consumption behaviors with respect to the identified structure. To this end, we used two-thirds of the consumption period (12 months of consumption during 2015) and all the available exogenous covariates to learn the mixture parameters, including the meter clustering. Then, we used the remaining data (corresponding to 7 months of consumption during 2016) to test the forecasting capabilities of the proposed method.

The exogenous variables used are the following climatic factors and calendar events (school holidays or summer vacations):

• T t -1 is the temperature (°C) at the previous timestep;

• P t -1 is the precipitation (in mm) at the previous timestep;

• C comprises the calendar events (school holidays or summer vacations) encoded using a binary variable at the current timestep; • Y t -1 is the hourly consumption level at the previous timestep.

Although other demographic and socioeconomic information could improve the richness of the analysis, such data were not available at the time of writing. For clustering purposes, all the introduced exogenous variables are used as model inputs. However, to examine the influences of these factors on the forecasting capability of the proposed method, the input vector e can contain one or more of these variables. The following two sections describe the results of these stages.

Learning and clustering

The graphical representation of the proposed method (see Figure 4) shows that each meter issues from a conditional mixture density function. To be able to compute this density function for the meters in Figure 5 (a), the number of mixture components should be known beforehand. After this parameter is estimated, the learning procedure can be initiated, and the consumption behavior dynamics can be captured by considering the exogenous factors. Note that we set the hyperparameter ξ to 10 -8 . The following sections examine these steps in more detail by providing the corresponding results.

Number of components.

The number of components G of the mixture must be set beforehand, as is the case for most unsupervised clustering algorithms. To estimate G for the proposed model, we computed the ICL criterion introduced earlier in the modeling section for a varying number of components G ∈ {1, . . . , 30}. The corresponding results are shown in Figure 6. A better model is obtained when the associated ICL value is minimum. In Figure 6, the minimum value of the ICL criterion is highlighted in red and corresponds to G = 13. In this figure, the meters (rows) are reorganized to reflect their cluster membership. The identified clusters are separated using horizontal red lines. Note that some clusters are made up of meters with more homogeneous consumption behaviors (e.g., clusters five and eight). To analyze the obtained clusters in more depth, the next section introduces some graphical tools and quantitative criteria.

4.1.3

Clustering sensitivity with respect to the number of states. In this section, we have conducted a sensitivity analysis of the model quality with respect to the number of states. This analysis is performed using simulated data which are constituted from three main Markovian dynamics (G = 3) and eight states (K = 8). We retrospectively varied the number of states from 2 to 20 increased by steps of 2, an operation that leads to the creation of 10 categorical datasets. Thereafter, using the proposed method (MixJNMM), we have classified the sequences in each of these datasets by fixing G to 3. Finally, we have computed the accuracy between the clustering results obtained on each dataset and the ground truth clusters of the simulated data (see Figure 7). In our case, the accuracy is estimated through the rate of the correctly classified sequences. As can be seen in Figure 7, the performances do not fluctuate significantly from K = 8. This can be explained by the fact that, increasing the number of states (more than 8) may result to the emergence of redundant states, which does not improve the clustering result. 4.1.4 Analysis of the clustering results. One way to represent the obtained clusters is to plot the proportion of states per cluster during the consumption time period and look for significant changes in consumption behavior. In Figure 8, the state proportions are illustrated for each cluster along with the school holidays and summer vacations, which are surrounded with red and cyan dashed rectangles, respectively. The distinctions between the proportions of the consumption behaviors at different time periods is immediately apparent. As shown in Figure 8, the state proportions change significantly during summer vacations for clusters C 1 and C 2 . This change is highlighted by a significant emergence of the eighth state, which has a later morning peak compared with states 1 and 2 (both of which are prevalent in these clusters).

Taking a different case, the meters of cluster C 4 exhibit a double behavior over time-the significant emergence of state 6 after the month of June 2015. Then, this state persists until the next summer. The symmetrical behavior of state 6 could be induced by some scheduled task (e.g., a scheduled maintenance operation) undertaken starting from this date.

The more homogeneous clusters (i.e., C 5 and C 8 ) exhibit a stable consumption behavior over time that does not change significantly. Cluster 8 consists mainly of the state 7, which is characterized by a one-peak pattern during weekdays and more constant consumption during weekends. The meters in this cluster may be associated with a commercial activity zone.

To summarize global cluster behavior, we used the estimated Markov model transition matrices. Figure 9 represents the transition matrices associated with three of the thirteen identified clusters. These are averaged transition matrices over time. As in the case of the time-variant Markov model, the transition probabilities between states belonging to any pair of consecutive timestamps are held in a unique transition matrix. In these matrices, the transition probabilities are encoded using a range of [0, 1]-from white (the lowest probability) to black (the highest probability)-by K ℓ=1 π дℓk = 1. The rows of these matrices indicate the states at time t -1, whereas the columns indicate the states at time t. In each matrix, the predominant states are surrounded by orange boxes.

By observing the general behaviors of these matrices, we can immediately perceive the difference between homogeneous and less-homogeneous clusters. In homogeneous clusters, the transition probabilities accumulate in a specific zone (e.g., cluster 8), whereas in the less-homogeneous clusters, the transition probabilities are shared between a higher number of states (e.g., clusters 2 and 9).

Considering the transition matrix of cluster 2 (Figure 9 (a)), we can see that there is a high transition probability from all the states towards state 1, which is one of the predominant states in this cluster. State 8, which is another predominant state in cluster 2 during summer vacations exhibits a high probability of staying in the same state. To summarize, most of the consumers associated with this cluster change their consumption behavior at the onset of summer vacations and then maintain that same behavior until the summer vacation ends.

Considering the more homogeneous cluster 8 (Figure 9 (b)), most of the transitions are towards state 7, which is the only predominant state in this cluster. The meters associated with this cluster tend not to change their consumption behavior over time, which acts as a confirmation that they belong to a spatial zone with commercial activity.

In addition to graphical representations, using a quantitative metric that measures the variability within clusters can also be helpful. For this purpose, we propose an entropy measure based on Markovian dynamics. For each observation i belonging to a cluster д, the entropy is determined as follows:

H (s i , д) = - t ℓ k P (д) k ℓ (e it ) log P (д) k ℓ (e it ), (20) where P (д)
k ℓ (e it ) indicates an input-dependent transition probability from state k to state ℓ for meter i belonging to group д. To compute the entropy at the cluster level, the mean entropy over the meters belonging to each cluster are computed as follows:

H (s д ) = 1 n д i ∈д H (s i , д), (21)
where n д is the number of meters in cluster д. As a result, the entropy measure is independent of the cluster size. Higher entropy values for a cluster д correspond to higher variability in the consumer consumption behaviors associated with that cluster. The entropy distribution of the meters with respect to the clusters of Figure 5 (b) is shown in Figure 10 through boxplots. In addition, the water consumption (in cubic meters) associated with the identified clusters is superimposed over the boxplots and visualized with respect to the right vertical axis. To facilitate the interpretation, we divide the entropy distribution values into the following three regions: (i) low variability (entropy < 2), (ii) medium variability (2 ≤ entropy ≤ 5), and (iii) high variability (entropy > 5). Fig. 10. Entropy distribution of meters per clusters using boxplots and mean water consumption (in cubic meters) superimposed using blue filled curves and plotted with respect to the right-hand y-axis Looking at the entropy distribution in Figure 10, clusters 1 and 2 are characterized by medium variability in their consumption behaviors, and their water consumption is much than that of the other clusters. This confirms the fact that residential housing accounts for a large amounts of water consumption and their behaviors do not change much. In contrast, cluster 4, which encompasses only a few meters, exhibits medium-to-high variability in consumption behavior and has the lowest consumption volume. This follows the interpretation given previously when discussing Figure 8. We can observe that a very low entropy distribution exists for the homogeneous clusters 5 and 8, which consist mainly of one state; however, cluster 8 includes a higher number of meters and, consequently, a higher consumption volume.

Another interesting pattern that can be observed regarding the entropy is the relationship between cluster size, water consumption volume per cluster, and the variability inside the clusters. In most cases, high variability is associated with the clusters that include fewer meters, yet at the same time, present a lower consumption volume. These clusters present irregular behaviors over time and can be the subject of further analysis.

Analysis of the influence of the context variables.

To conduct this analysis, we relied on the estimated coefficients of the model to provide more details concerning the influence of the exogenous factors on consumption behaviors. To this end, we have extracted the estimated coefficients β corresponding to climatic variables (i.e., temperature and precipitation) for the most significant transitions of two different clusters (i.e., C 2 and C 9 ). First, let us recall that, within a cluster д, the transition probability from a state ℓ to a state k is defined by

P(s it = k | s it -1 = ℓ, w i = д, e it ) = π дℓk (e it ; β дℓ ) = exp (β ⊤ дℓk e it ) K k ′ =1 exp (β ⊤ дℓk ′ e it )
.

The above mentioned coefficients for clusters 2 and 9 are shown respectively in the tables 1 and 2. In these tables, each value corresponds to the coefficient of the above logistic regression model, for a climatic variable, given transitions from the states located in the rows to the states located in the columns. The most significant values are highlighted in bold face in these tables. Cluster 2 is dominated by the state 8 during the summer period and by state 1 during the remaining period (see figure 5 (b)). Regarding cluster 9, it is constituted mainly by states 1, 5 and 8.

Looking at Table 1 (a), we can notice a high positive coefficient value for the transitions from state 1 to states 5 and 8. An increase in temperature therefore results in an increase of such transitions (occurring mainly when entering the summer period). Looking at the Table 1 (b), regarding the precipitation variable coefficients, we can notice a high positive value for the transition from the state 8 to the state 1, following which an increase in rainfall may lead to an increase of such transitions (occurring mainly when entering the autumn period). Additionally, we can notice the high negative values for transitions from states 1 and 8 to state 2, according to which a low precipitation rate may lead to an increase of such transitions.

Table 2 represents the transition coefficients of cluster 9 corresponding to the above-mentioned climatic variables. Looking at the Table 2 (a), we can notice a high positive value for the transition from the state 8 to itself and a high negative value for the transition from the state 5 to the state 1. This would signify that high temperatures lead to an increasing number of transitions from state 8 to itself and low temperatures lead to an increasing number of transitions from state 5 to state 1. Looking at Table 2 (b), we can notice a high positive coefficient value for the transition form state 8 to state 5. It means that a high precipitation rate leads to an increasing number of such transitions for this cluster.

Forecasting

Using the mixture parameters and the input-dependent Markov model transition matrices estimated during the learning process, future consumption behaviors can be predicted within each cluster.

As mentioned earlier, we used the first 7 months of data from 2016 to test the performance of the proposed method. Here, we provide the details of the experimental setup used to evaluate the performance of the proposed method in terms of the prediction error.

Evaluated methods.

The performance of the proposed method in terms of forecasting accuracy is evaluated with respect to four other methods:

• The homogeneous Markov model (MM) is a time-independent model and assumes that a single transition matrix summarizes the evolution of consumption behavior dynamics for all the categorical sequences; • The mixture of homogeneous Markov models (MixMM) fits a homogeneous Markov model within each cluster to model the dynamics of the consumption behaviors; • The joint nonhomogeneous Markov model (JNMM), unlike the MM, considers the temporal regularity of states; • The K-means+JNMM is a combination of the k-means algorithm and the nonhomogeneous Markov model. The initial time series (with an hourly frequency) are clustered in the first step by the k-means algorithm; then, the categorical series within each cluster are predicted using the nonhomogeneous Markov model.

Evaluation criteria.

To compare the forecasting results, we compute different criteria for each configuration as follows:

• The adjusted Rand index (ARI) [START_REF] Hubert | Comparing partitions[END_REF] is a measure of agreement between two clustering results and is a multi-class criterion; • Accuracy is the percentage of correctly predicted daily consumption behaviors; • Precision is defined for each category by the rate of correctly predicted instances among the instances affected by this category; • Recall is defined for each category by the rate of correctly predicted instances among the instances of that category; • the F-measure represents the harmonic mean of precision and recall. 3 summarizes the comparison results between the aforementioned methods when using different combinations of input variables. Note that no input variable exists for the methods based on the homogeneous Markov model (i.e., MM and MixMM) because it is inherently time-invariant. The evaluation results obtained by each criterion appear as bars inside each associated table cell. All the evaluated criteria can take values within [0, 1], where 1 indicates the best performance. The bar lengths are designed to be proportional to the values of the metrics and to facilitate the interpretation of the results. The best performances are highlighted with red bars. and consumption level). Figure 8 illustrated how calendar events can affect consumer consumption behavior over time.

Comparisons. Table

The proposed MixJNMM method obtains the best prediction results. The proposed method may give better results due to the presence of heterogeneous latent structure among the different categorical time series. The mixture property of the proposed method allows for the proper identification of this structure. Its Markovian property provides us with the ability to capture the joint dynamics (state-state and state-input) of the consumption behaviors.

Figure 11 provides greater detail concerning the forecasting results when using the proposed method. In this figure, the prediction errors are measured for each identified cluster and are expressed as percentages. We can see that the prediction error changes significantly between clusters. For the clusters that exhibit homogeneous consumption behaviors over time (i.e., clusters 5 and 8), the prediction error remains very low. In contrast, as the consumption behavior variability increases inside the clusters, the prediction error also increases.

DISCUSSION AND CONCLUSION

In this paper, we propose a novel methodology based on a mixture of Markov models for analyzing consumer behavior. Using the methodology proposed in this study, categorical time series representing sequences of consumption behaviors were classified into homogeneous groups while considering the joint Markovian dynamics. Subsequently, we predicted the future consumption behaviors of each identified cluster using a set of input variables and cluster-specific estimated parameters. The numerical and graphical evaluations on a real-world dataset demonstrated the effectiveness of the proposed methodology for clustering and forecasting of water consumption time series. In this article, we have analyzed weekly consumption behaviors, however, the proposed method could potentially be extended to treat daily behaviors and can also be applied to other application domains involving categorical time series.

The proposed method allows clustering of water consumers and forecasting of their future consumption habits. This method can help water utilities in different ways to better manage the water network. By clustering the consumers in different groups, inside which they share similar consumption habits, it can accelerate the decision making process and forecasting the future consumption habits could be necessary in some critical conditions such as a drought period for better distribution of this valuable resource. It provides also water utilities with the ability of targeting a group of consumers based on their consumption behavior dynamics.

As an insight for future works, in order to improve the forecasting accuracy of the proposed method, the transition dynamics can be learned on a set of homogeneous segments in each cluster. It requires the model modification, following which, each segment should be composed of adjacent time instants in order to respect the time dependence constraint of the proposed method. Additionally, spatial localization of smart meters may be used in future works to expand the analysis of consumption behaviors and to assign spatial implications to clusters. The proposed method may also be adapted to detect consumption behavior changes and conduct change analyses. "

Fig. 1 .

 1 Fig. 1. Scheme of the proposed methodology, which consists of 4 steps: (1) preprocessing and normalization, (2) data discretization per week, (3) modeling consumption behavior dynamics, (4) clustering of categorical time series and forecasting of future consumption behaviors.

Fig. 2 .

 2 Fig. 2. Clusters centers or main weekly consumption behaviors, called states (solid lines), with their standard deviations (brighter colors) on the left and the corresponding superimposed daily profiles on the right. The plots on the right-hand side clearly reveal how consumption behavior on weekdays (solid lines) differ from those on weekends (dashed lines).

Fig. 3 .

 3 Fig. 3. Categorical dataset encoding the weekly consumption behaviors of 2,000 smart meters

Fig. 4 .

 4 Fig. 4. A graphical representation of the proposed mixed joint nonhomogeneous Markov model (MixJNMM) with input variables

ALGORITHM 1 :

 1 MixJNMMData: (s 1 , . . . , s n ), (e 1 , . . . , e n ). Result: Parameters ϕ; Partition variable w. initialization: G; ϕ(0)

Fig. 5 .

 5 Fig. 5. Clustering process: (a) the consumption behavior sequences of 2,000 consumers and (b) the 13 clusters (separated by horizontal red lines) obtained by applying the proposed method.

Fig. 6 .

 6 Fig. 6. ICL computed using an interval of 1-30 clusters; estimating the number of clusters.

Fig. 7 .

 7 Fig. 7. Sensitivity of clustering result for different number of states K.

Fig. 8 .

 8 Fig. 8. State proportions over months per cluster. The states are visualized using different color codes associated with the signatures identified in Figure 2. The red dashed rectangles indicate school holidays and the cyan dashed rectangles indicate summer holidays.

Fig. 9 .

 9 Fig. 9. Average transition matrices associated with three of the thirteen identified clusters. Darker colors indicate higher transition probabilities. In each matrix, the predominant states are surrounded by orange boxes. The rows of the matrices indicate the states at time t -1, and the columns indicate the states at time t.

Fig. 11 .

 11 Fig.11. Prediction error with respect to the clusters using the proposed MixJNMM and the following inputs: temperature, precipitation and calendar events

Table 1 .

 1 Estimated coefficients of the most significant transition probabilities of cluster 2 for temperature

	and precipitation factors
	test test test test (a) Temperature coefficients			(b) Precipitation coefficients
	s t -1	s t	1	2	5	8	s t -1	s t	1	2	5	8
	1		0.17 0.19 0.30 0.28	1		2.13 -3.63 0.68 0.41
	8		-0.12 -0.06 -0.01 0.12	8		3.90 -4.77 2.30 -1.37

Table 2 .

 2 Estimated coefficients of the most significant transition probabilities of cluster 9 for temperature

	and precipitation factors
	test test test test (a) Temperature coefficients	(b) Precipitation coefficients
	s t -1	s t	1	5	8	s t -1	s t	1	5	8
	1		-0.11 -0.03 -0.01	1		0.46 1.75 0.44
	5		-0.17 -0.06 0.03	5		1.61 2.42 -1.49
	8		-0.03 0.08 0.14	8		0.46 3.89 -0.69

Table 3 .

 3 Comparison of the models in terms of the supervised metrics during the 26 test weeks: MM: homogeneous Markov model, MixMM: mixture of homogeneous Markov models, JNMM: joint nonhomogeneous Markov model, k-means+JNMM: joint nonhomogeneous Markov model within the clusters identified by the k-means algorithm, MixJNMM: mixture of nonhomogeneous Markov models, Y: meter water consumption, T: temperature, P: precipitation, C: calendar informationBy looking at this table, we can see that the methods depending on the homogeneous Markov chains (MM and MixMM) achieve low forecasting precision levels. In contrast, the models that take exogenous factors as input (JNMM, K-means+JNMM and MixJNMM) obtain their best results when the input vector e includes all the available covariates (temperature, precipitation, calendar events

	Models	Inputs (e it )	ARI	Accuracy Recall	Precision F-
							Measure
	MM		0.36	0.65	0.67	0.67	0.67
	MixMM		0.37	0.65	0.68	0.68	0.68
		(T t -1 , P t -1 )	0.57	0.66	0.67	0.67	0.67
		(Y t -1 )	0.59	0.71	0.72	0.73	0.72
	JNMM	(Y t -1 , C)	0.62 0.61	0.71 0.73	0.72 0.72	0.73 0.72	0.73 0.72
			0.63	0.74	0.73	0.74	0.73
		(T t -1 , P t -1 )	0.58	0.69	0.67	0.69	0.68
	K-means (Y t -1 )	0.61	0.74	0.74	0.72	0.73
	+	(Y t -1 ,T t -1 , P t -1 )	0.63	0.75	0.75	0.76	0.75
	JNMM	(Y t -1 , C)	0.62	0.74	0.75	0.74	0.74
			0.65	0.76	0.76	0.75	0.76
		(T t -1 , P t -1 )	0.59	0.71	0.73	0.71	0.72
		(Y t -1 )	0.64	0.79	0.77	0.78	0.77
	MixJNMM	(Y t -1 ,T t -1 , P t -1 )	0.67	0.80	0.81	0.79	0.80
		(Y t -1 , C)	0.64	0.79	0.78	0.76	0.77
		(Y t -1 , C,T t -1 , P t -1 )	0.70	0.82	0.81	0.80	0.80

(Y t -1 ,T t -1 , P t -1 ) (Y t -1 , C,T t -1 , P t -1 ) (Y t -1 , C,T t -1 , P t -1 )

ACM Trans. Intell. Syst. Technol., Vol. 1, No. 1, Article 1. Publication date: January 2019.

s 1,1 s 1,2 s 1,T s 1,T +1 s n,1 s n,2 s n,T s n,T +1 e 1,1 e 1,2 e 1,T e1,T +1 e n,1 e n,2 e n,T en,T +1 . . . . . . . . .

ACKNOWLEDGMENTS

This work was conducted within the framework of the project "Analyse de données massives relevées à distance sur les compteurs d'eau" involving Veolia Eau d'Ile de France, the Syndicat des Eaux d'Ile de France (SEDIF) and the French institute of science and technology for transport, development and networks (IFSTTAR).


Barnabé Potel
email: barnabe.potel@g2elab.grenoble-inp.fr

Vincent Debusschere
email: vincent.debusschere@g2elab.grendoble-inp.fr

Florent Cadoux
email: florent.cadoux@g2elab.grenoble-inp.fr

Leticia De

Alvaro Garcia

Under-frequency Load Shedding Schemes Characteristics and Performance Criteria

Keywords: Under-frequency Load Shedding (UFLS), Defense plan, Power system, Stability

Our work focuses on the European Network Code on Emergency and Restoration (NC-ER) which introduced an harmonization requirement of Under-frequency Load Shedding (UFLS) schemes. This new requirement implies that in the medium term, the UFLS scheme of all European countries will have to evolve. It stipulates an acceptable range for the main factors that define an UFLS scheme, namely: the number of load shedding steps, the percentage of load shed in each step, and the accuracy of the frequency measurement that is implemented in the protection relays. In this context, the contribution of this paper is twofold: first, we define a new performance criterion for UFLS schemes; and second, we apply this criterion to two opposite archetypal schemes, in order to highlight the pros and cons of both approaches.

I. INTRODUCTION

An Under-frequency Load Shedding (UFLS) scheme is designed to prevent a power system from a black-out due to a large power imbalance. Most current used schemes in the European grid follow the following traditional design: a predefined amount of load is shed automatically whenever frequency decreases below certain predefined thresholds [START_REF] Delfino | Implementation and comparison of different under frequency loadshedding schemes[END_REF].

The Network Code on Emergency and Restoration (NC-ER) [START_REF] Entso-E | Network Code on Emergency and Restoration[END_REF] is a legally binding document which compels ENTSO-E members to update their UFLS schemes in order to harmonize technical operations in the European power grid. When a large grid event which involves load shedding occurs, a balanced participation of each country seems to be the fairest and safest response. As a general rule, UFLS schemes that are currently implemented by European TSOs, and that were inherited from the past, do not comply with the latest NC-ER and will have to be updated.

Two archetypal examples of UFLS schemes are the following. In what could be called the "current" scheme, only a few UFLS steps are used by each TSO, the associated frequency threshold being different from one TSO to another, the amount of load shed in each step is large and frequency measurement is relatively inaccurate due to the use of protective relays designed decades ago. Conversely, in what could be called the "future" situation, many UFLS steps containing only relatively

The work reported in the paper has been developed in the framework of the Enedis Industrial Chair On Smart Grids research program. little load would be used, and the accuracy of frequency measurement would be improved. When implementing the NC-ER, national TSOs will have to choose the parameters of their UFLS scheme inside the allowed range and to make sure that the new settings actually improve system stability with respect to the current situation; this raises the question of evaluating the performance of a given UFLS scheme.

As an example, "Scheme A", defined by Table I, illustrates the type of scheme that is currently used by TSOs. It is composed of four UFLS steps. Such a traditional scheme typically allows a relatively large inaccuracy in frequency measurement; hence this scheme is not compliant with the guidelines described in the NC-ER (the new UFLS scheme shall be designed with six steps between 48 and 49 Hz, each step should shed between 5 and 10 % for a total of 45 % shed load). In that context, the objective of this work is to compare the performance of a typical current UFLS scheme with the performance of a newer NC-ER-compliant scheme. To this end, we will first introduce a new performance metric for UFLS schemes, and then apply this metric to perform our comparison.

A prospect of this UFLS design, not studied in this paper, is the future challenge due to the steady increase of distributed generation: in this context, disconnecting an entire MV feeder (both load and generation at once) by triggering an UFLS relays in the primary substation may prove too brutal. As a consequence, new solutions will probably be required in order to preserve the efficiency of the load shedding mechanism [3]- [START_REF] Xie | Distributed underfrequency load shedding using a multi-agent system[END_REF].

II. UFLS MECHANISM

A. Principle

Different kinds of events may cause the frequency to drop to the point where automatic load shedding triggers. For example, large generators may trip, or transmission lines may disconnect leading to the grid separation into subareas with potentially large power imbalances. Then frequency decreases according to the swing equation [START_REF] Delfino | Implementation and comparison of different under frequency loadshedding schemes[END_REF]. The primary reserve of frequency regulation may not be sufficient to contain the frequency drop in such a severe situation; whenever power consumption exceeds power generation so much that the imbalance cannot be dealt with by quickly ramping up generators, some consumers must be immediately disconnected from the grid to contain the contingency. As a consequence, whenever frequency reaches the first threshold of the UFLS scheme (which is commonly set to 49 Hz), a fraction of the consumption is shed. If frequency keeps dropping, additional load is shed at lower frequency thresholds. Under-frequency load shedding is thus the last action to avoid an extensive black-out.

B. UFLS scheme design

The steps of an UFLS Scheme are defined by their frequency threshold and the percentage of load to shed. After enough load is shed, frequency eventually goes up and the frequency drop is contained. The amount of load shed should thus be sufficient, but not excessive, since an overly aggressive UFLS scheme would actually be counter-productive. Indeed, if too much load is shed, the under-frequency situation may be followed by an over-frequency situation that could lead to the disconnection of generators [START_REF] Entso-E | Network Code on Requirements for Grid Connection Applicable to all Generators[END_REF]. In such a case, the UFLS mechanism would overshoot and actually worsen the power generation deficit.

Fig. 1 shows an example. The simulation was based on equations that will be presented below in Section III. An imbalance of 20 % between generation and consumption is simulated and three schemes triggering at 49 Hz are compared. The first one is shedding 19 % of the consumption, the second one 20 % and the last one 22 %. The final value of frequency is reported in the legend of the graph on Fig. 1; this value may be determined by a longer simulation or simply by using equations ( 2) and (3) below. By contrast with the UFLS Scheme A, Table II shows one of the many possible schemes that comply with the NC-ER. Note that this scheme contains six steps, two more than Scheme A.

III. POWER GRID MODELING FOR UFLS STUDIES A. Grid model

The grid is modeled by a single synchronous machine which provides the generation, inertia and primary frequency regulation. The secondary frequency regulation is not taken into account because its time constant is too large to be considered.

Swing equation: the frequency evolution is ruled by the swing equation:

df dt = ∆p 2H (1)
f frequency in p.u. ∆p generation and consumption difference in p.u.

H inertia

Primary frequency regulation: The frequency regulation is represented by a proportional gain g d and saturates at 49.8 Hz and 50.2 Hz [START_REF] Entso-E | Network Code on Load Frequency Control and Reserves[END_REF], followed by a first-order low-pass filter with time constant τ . The symbol s denotes Laplace-transform.

g = g 0 1 + g d 1 -f 1 (1 + τ s) (2)
g 0 initial generation (after triggering event) g d generation droop τ time constant for generator response Consumption: We modeled the load as depending linearly from frequency [START_REF] Concordia | Load representation in power system stability studies[END_REF].

p = p 0 1 + k pf f -1 (3)
p 0 initial consumption considering the load shedding k pf frequency-sensitivity of consumption Simulations were conducted using the following values:

H 6 seconds g d 5 %/Hz k pf 1 %/Hz
A sensitivity study with respect to these parameters is provided in the next section.

B. Relays

On November 4th, 2006, an incident led to triggering the first step of the UFLS scheme in several European countries. It revealed that a part of the relays contained in the first step did not trigger because of frequency measurement inaccuracy [START_REF]Rapport d'enquête de la commission de régulation de l'énergie sur la panne d'électricité du samedi 4 novembre[END_REF].

This behavior is added in our model, considering that the load shed at a step follows a normal distribution depending on the frequency. The standard deviation of this distribution is directly linked to the frequency measurement accuracy. The tripping delay of our "inaccurate relays" is considered as depending on the accumulated time spent in the vicinity of the frequency threshold: the longer frequency remains at a certain value (frequency threshold ±3 standard deviation) the more relays trigger. See Fig. 2 for an example of this behavior.

C. Premature loss of generation caused by low frequency

During an under-frequency situation, generator tripping is generally assumed to begin at 47.5 Hz. In fact, the event that occurred on November 4th, 2006 in Europe and the countrywide black-out that occurred on September 28th, 2003 in Italy showed that generators may disconnect before frequency actually reaches the 47.5 Hz limit, and sometimes even before it reaches 49 Hz. We do not implement these situations explicitly in our simulation, since they may simply be considered as a larger initial loss of generation.

Moreover, most of the generators which tripped during the event of 2006 were meeting grid requirements (tripping at those frequencies was their expected response); but these requirements have evolved since then, and generators, nowadays, should not trip before reaching 47.5 Hz. Nevertheless, the physical design of generators (that depends on choices made by their manufacturers) may be such that the actual range of frequency supported by the machine is smaller than the range of frequency implemented in the protection relays (which depends on grid requirements). Thus the possibility of losing generators before reaching 47.5 Hz or 51 Hz may still exist.

IV. COMPARISON OF VARIOUS UFLS SCHEMES

A. Robustness of an UFLS scheme with respect to the triggering event

As developed in the introduction, under-frequency situations may be caused by various contingencies. In particular, they may arise due to the loss of transmission lines causing the network to split into subareas. It is important to notice that the parameters of these sub-networks cannot be determined in advance; hence, no assumption should be made regarding (A) the amount of imbalance that the UFLS scheme will have to deal with, nor (B) the characteristics (such as inertia) of the network after the event, nor (C) the location of the resources (such as MV-feeders or entire HV/MV substations) that the UFLS will be able to rely on.

The last point (C) leads to the relatively straightforward conclusion that UFLS resources should be spread relatively evenly in the network. It would probably not be a good idea, for instance, to locate all the resources of the first UFLS step in the North of a country, and all the resources of the second UFLS step in the South.

The first and second points (A) and (B) are however less obvious. Ensuring that an UFLS scheme will be robust with respect to the total initial imbalance and the system inertia calls for detailed analysis, which will be the object of the following sections.

B. Performance criterion

The performance metric used by ENTSO-E to support the new requirements is [START_REF] Entso-E | Technical background for the Low Frequency Demand Disconnection requirements[END_REF] the steady-state value reached by the frequency after the load-shedding : the scheme is considered as efficient if this value is close to the initial value of 50 Hz for 8 different imbalances (1, 5, 10, 20, 30, 40, 50 and 60%).

Taking a different perspective at the matter, we argue that it would also be reasonable to measure the performance of an UFLS scheme by the maximum frequency which can be reached after load shedding, for any imbalance lower than the total amount of load in the UFLS scheme. By comparison with the performance criterion of ENTSO-E [START_REF] Entso-E | Technical background for the Low Frequency Demand Disconnection requirements[END_REF], this means two things:

• first, we are checking the value of frequency for the whole range of possible imbalances, and not only for a few predefined values of the imbalance. • Second, we are checking the value of frequency at any time after load shedding occurs, not only its steady-state value. The rationale of our performance metric is the following: UFLS is an automatic and fast response to re-establish power balance, it should thus be judged by its ability to maintain system stability and not by its ability to return frequency to its nominal value. If the steady-state frequency value after load shedding stabilizes to a low (yet above 47.5 Hz) value, this is not detrimental to the power system: the point is that a full black-out was avoided, and TSOs will then be able to take actions to get the frequency back to its initial value, by starting up additional generation units or manually shedding some more load if necessary. As a consequence, a given UFLS plan should not be penalized based on the fact that it leads to an under-frequency situation. On the contrary, over frequency situations should be avoided because generators as a whole will disconnect much faster in case of over-frequency. Avoiding over-frequency situations, at any time (not only in steady-state) and for any reasonable value of the initial imbalance, is thus the aim of our performance criterion.

The value of the maximum frequency after load shedding depends on relays characteristics and grid parameters such as inertia. If the grid splits into subareas, each subarea may have its own characteristics which are not predictable. We thus need to carry out a sensitivity study with respect to the three following important grid parameters: system inertia, droop constant and frequency-sensitivity of consumption.

The principle of our work is to time-simulate various imbalances and to plot the maximum value of frequency that is reached after load shedding for each scenario. Fig. 3 shows how the maximum frequency changes when loss of generation is increased. Two schemes are compared: one is representative of "current" schemes, and the other is representative of "future" schemes that will fulfill the NC-ER requirements (we used the schemes shown in Table I and II).

Curves start from 3 % (as the primary reserve is around 1 % and the consumption decrease is around 1 % too for triggering the first UFLS step at 49Hz). These curves look like saw-teeth. Their number is directly linked to the number of steps for each scheme: six for Scheme B and four for Scheme A. It is clear that taking a limited numbers of imbalances does not allow to have an overview of the performances of an UFLS scheme.

Local minima are observed on Fig. 3 whenever the loss of generation reaches a value such that a new step of the UFLS scheme must be activated; in such a case, the maximum frequency after shedding (on the y-axis) increases quickly with the amount of generation lost (on the x-axis). Indeed, losing only slightly more generation will lead the frequency to reach a slightly lower value and will potentially cause a Fig. 3. Frequency progression for different shed percentages large number of additional UFLS relays to trigger. In the extreme situation where all relays would trigger exactly at the same frequency, that is to say, if the frequency measurement inaccuracy was not taken into account, the curves depicted on Fig. 3 would go up vertically after each local minimum; local maxima would be reached immediately after local minima; and whenever a minimum would be attained, one would read the frequency threshold of the corresponding UFLS step on the yaxis. However, because of frequency measurement inaccuracy, the slope of the curve after local minima is not vertical, local maxima are reached about halfway between two UFLS steps (not right after a step), and minima are reached for frequency values that are slightly higher than the theoretical frequency thresholds of the UFLS scheme.

Fig. 3 shows that the "current" scheme, where the thresholds of UFLS relays are distributed with a standard deviation of 50 mHz [START_REF]Rapport d'enquête de la commission de régulation de l'énergie sur la panne d'électricité du samedi 4 novembre[END_REF], never exceeds the nominal frequency of 50 Hz: as a consequence, there is no risk of over-frequency situation. In comparison, the "future" scheme, where the standard deviation is only 10 mHz (representing 30 mHz precision according to [START_REF] Entso-E | [END_REF]), drives frequency to values well above 50 Hz. The newer scheme benefits from having more steps of smaller size than the "current" scheme; but the fact that frequency is measured more accurately favors the occurrence of over-frequency situations. The main differences between these schemes are shown in Table III.

The results of Table III may be generalized, but the curves from Fig. 3 tremendously depend on grid and relays parameters (inertia, accuracy of frequency measurements, etc). If those parameters vary, the curves will keep approximately the same shape but may be shifted upwards along the frequency axis. Indeed, they are lower-bounded by the characteristics of their UFLS schemes (recall that the minima on Fig. 3 derive from the UFLS steps in a relatively straightforward manner), which prevents them from being shifted downwards; but they may reach high values of frequency, potentially above Maximum frequency contained due to the significant inaccuracy reaches higher values Last frequency step 47.5 Hz -risk of generators tripping for low frequencies is higher before shedding the total load 48 Hz 50 Hz. We will thus now turn to studying how the maximum frequency of Fig. 3 changes when the main parameters of the simulation are changed for each UFLS scheme.

C. Sensitivity of the max frequency reached post-shedding Fig. 3 shows the maximum frequency reached after loadshedding for various imbalances. In this section, we will focus on the maximum frequency reached for a given scheme, when the initial loss of generation varies. For instance, Fig. 3 shows that for the value of inertia H, droop constant g d and load sensitivity k pf that are indicated in Section III, the maximum frequency that could be reached for UFLS Scheme A is about 49.50 Hz, and it is about 50.25 Hz for Scheme B. Now, how do these two values change when parameters H, g d and k pf vary? Note that the initial imbalance that leads to the maximum frequency may not be the same for different values of these three parameters. In the following figures, we will vary one (and only one at a time) of the three parameters under consideration, and plot the maximum postload-shedding frequency with respect to this parameter. The other two parameters are kept constants; they are set to the values indicated in Section III, namely: H = 6 seconds, g d = 5 %/Hz, k pf = 1 %/Hz. Note that our model does not capture the phenomenon of generators tripping due to excessive frequency; as a consequence, the following graphs display unrealistic high values of frequency (that would be prevented in practice by the fact that generators would trip). In our studies, these values indicate that we detected a "danger zone": values of parameters H, g d and k pf that could lead to an overshoot of the UFLS scheme under consideration.

1) System inertia: power system inertia, H, is an energy stock provided by rotating machines in the network. The larger this energy stock, the smaller the time derivative of frequency (for a given imbalance); see Equation [START_REF] Delfino | Implementation and comparison of different under frequency loadshedding schemes[END_REF].

Fig. 4 shows the evolution of the maximum frequency reached for the two UFLS schemes previously compared. If inertia is larger, the system response will be long, fewer relays will trip because of frequency measurement inaccuracy, and the final value of frequency will be lower, which is beneficial. Conversely of inertia is low, frequency will evolve fast, many relays will trigger, frequency will reach a potentially high value after load shedding.

Particularly, the effect of inertia is stronger on Scheme A because UFLS steps are larger. Thus, for low values of inertia, Scheme B appears better than Scheme A (the beneficial effect of having many steps dominates), while Scheme A appears better than Scheme B for high values of inertia (the beneficial effect of using less accurate frequency measurements dominates). 2) Generators droop: the generator droop, g d , is the primary reserve of the frequency regulation. Interestingly, the qualitative effect of this parameter is not the same for both schemes, see Fig. 5. This may be explained simply by the fact that droop control tends to drive frequency towards (rather than "away from") 50 Hz, as follows. When Scheme A is used, the post-load-shedding frequency is lower than 50 Hz. As a consequence, the effect of droop control when Schema A is used is to increase post-load-shedding frequency (hence the fact that the curve for Scheme A is nearly monotonically increasing on Fig. 5).

Conversely, when Scheme B is used, the post-load-shedding frequency is higher than 50 Hz, because more accurate relays are used and such relays tend to trigger all-at-once; so the effect of droop control is to decrease frequency (hence the fact that the curve for Scheme A is monotonically decreasing on Fig. 5).

From a certain value of the droop constant, there is no effect on the maximum frequency reached after load shedding (this is particularly visible for Scheme A on Fig. 5). The droop is so important that frequency stability of the system is guaranteed even with imbalances which lead to apply load shedding.

3) Frequency sensitivity of consumption, k pf : on one hand, the primary frequency regulation is almost a proportional gain (limited between ± 0.2 Hz around 50 Hz) which allows the system to reach a power balance. On the other hand, the load tends to decrease when frequency drops, which further favors the power balance. See Fig. 6 for the effect on both schemes of this parameter. The higher k pf , the better from the viewpoint of system stability. On the contrary, if load is not sensitive to frequency (or if the sensitivity is even negative), then its effect is detrimental from the viewpoint of system stability. This behavior is the same for both UFLS schemes. Nevertheless the effect if more visible for Scheme A than for Scheme B: the reason is the same that for the inertia.

V. CONCLUSIONS

The performance criterion defined in Section IV-B allows us to have a basis to compare the two archetypal types of UFLS schemes: the scheme that has been commonly used in the last decades, corresponding to Scheme A, and the future one (yet to be implemented), corresponding to Scheme B. The curve of Fig. 3 representing maximum frequency reached post-shedding in function of the power imbalance looks like saw-teeth for both schemes. Then, whichever scheme is implemented, the frequency reached may correspond to one of the frequency thresholds of the UFLS scheme (the threshold reached will depend on the power imbalance). Thus, it is not possible to design a UFLS scheme in order to regulate the frequency. Moreover, as previously mentioned, an UFLS scheme is the last part of frequency stability and has not a regulation aim but a protection/defense goal of the electrical system.

The sensitivity study depending on grid parameters shows that inertia H and the load frequency dependence k pf have the same beneficial effect on maximum frequency reached for both schemes. Scheme A is however more sensitive than Scheme B to these two parameters, because of its inaccuracy in frequency measurement. For our reference values of H and k pf , the droop parameter g d has an opposite effect on the two schemes: when g d is increased, the frequency reached with Scheme A increases towards 50 Hz while the frequency reached with Scheme B decreases towards 50 Hz. Because of differences in accuracy of frequency measurement, the maximum postshedding frequency reached with Scheme A is relatively low, while Scheme B presents a higher risk of over-frequency situations. This explains why the frequency reached after the shedding that occurred on November 4th, 2006 in Europe was just over 49 Hz. It was due of UFLS schemes of countries involved in shedding which were of the type of the Scheme A developed in this paper.

In the future, the gradual replacement of existing UFLS relays with higher-accuracy ones should be done with care, in order to avoid unintentionally degrading the performance of the UFLS mechanism.

Fig. 1 .

 1 Fig. 1. Evolution of frequency for different percentages of load shed.

Fig. 1

 1 Fig. 1 demonstrates that our three schemes exhibit very different behaviors.

Fig. 2 .

 2 Fig. 2. Percentage of triggered relays of a step for a frequency threshold at 0 Hz and a standard deviation of 0.1 Hz

Fig. 4 .

 4 Fig. 4. Evolution of frequency for different values of system inertia, H.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Evolution of frequency for different generators droop constants, g d .

• Scheme 3 is shedding too much load and leads to an over-frequency situation. The most sensitive generators may disconnect. • Scheme 1 simply contains the frequency drop and the final frequency is lower than its initial value.• Scheme 2 not only contains the frequency drop, but also drives frequency back to its initial value. Nevertheless, a small temporary overshoot in frequency may be observed; this is due to the frequency regulation. In section IV-B below, we will elaborate on this example to devise our performance criterion for UFLS schemes.