Mauricio Gonzalez
email: mauricio.gonzalez@outlook.fr

Antoine Girard
email: antoine.girard@centralesupelec.fr

Data-Driven Modeling of the Temporal Evolution of Breakers' States in the French Electrical Transmission Grid

Keywords: Data-driven modeling, Hidden Markov Models, EM algorithms, Electrical transmission grid

In electrical transmission grids, it is common to observe the states of circuit breakers. While they are known at irregular times, system modeling and grid state estimation are of the highest importance to ensure secure operations. This paper proposes a richer method to estimate the grid state over its reference configurations based on the temporal evolution of its breakers' states. The first contribution consists in developing a general multi-observation continuous-time finite-state Hidden Markov Model with filter-based parameter estimation to infer the hidden state (e.g., the grid reference configuration) handling multiple observed processes with irregular "jump" times (e.g., the breakers' states). As a second contribution, we build a numerical scheme with no discretization error adapted to all state jumps generated by the observed processes. Finally, we apply our model to simulated and real data to illustrate the approach's performance. The available data consists of historical records of breakers' states during the electrical transmission grid operated normally. For this real-data-driven application, we also present a clustering approach to identify the set of grid reference configurations.

Introduction

An electrical transmission grid is an interconnected network that permits the electrical energy movements from producers, e.g., nuclear plants, to electrical substations. Central components of these networks are circuit breakers, which are electrical switches designed to interrupt or continue the electrical energy flow.

Hence, the grid configuration is determined by all breakers' states, which can be used to adapt the grid to diverse operating conditions. However, while the breakers' states are known at all irregular times, the current grid configuration is assumed to be unknown. As the grid can be prone to failures or malicious attacks, system modeling is highly important to ensure secure operations [START_REF] Wang | Cyber security in the smart grid: Survey and challenges[END_REF][START_REF] Sun | Cyber security of a power grid: State-of-the-art[END_REF]. This paper proposes an adhoc Hidden Markov Model (HMM) [START_REF] Elliott | Hidden Markov models: estimation and control[END_REF] to estimate the state temporal evolution of the grid over its reference configurations that maximize the likelihood of observing the different breakers' states.

The literature has generally investigated system state estimation in the discrete-time framework under the (most recurrent) assumption of partially observing the state in Gaussian noise [START_REF] Mo | Detecting integrity attacks on scada systems[END_REF][START_REF] Mo | Secure control against replay attacks[END_REF][START_REF] Mishra | Secure state estimation against sensor attacks in the presence of noise[END_REF][START_REF] Amin | Safe and secure networked control systems under denial-of-service attacks[END_REF][START_REF] Befekadu | Risk-sensitive control under markov modulated denial-of-service (dos) attack strategies[END_REF][START_REF] Forti | Joint attack detection and secure state estimation of cyber-physical systems[END_REF][START_REF] Mo | False data injection attacks against state estimation in wireless sensor networks[END_REF][START_REF] Guo | Optimal linear cyber-attack on remote state estimation[END_REF][START_REF] Shi | Causality countermeasures for anomaly detection in cyber-physical systems[END_REF][START_REF] Shi | On finite-state stochastic modeling and secure estimation of cyber-physical systems[END_REF][START_REF] Shi | Event-based state estimation of linear dynamic systems with unknown exogenous inputs[END_REF]. However, applying this standard modeling is not appropriate for our problem. In fact, we need to consider two specific characteristics of our framework. First, breakers' states switch between the modes "off /on," which may happen at any time. Hence, our model must handle breakers' state changes at any time, which is a feature of the continuous-time modeling. Second, we observe several breakers. Hence, we need to extend standard HMMs [START_REF] Elliott | Hidden Markov models: estimation and control[END_REF][START_REF] Aggoun | Finite-dimensional models for hidden Markov chains[END_REF][START_REF] Elliott | Discrete-time expectation maximization algorithms for Markov-modulated poisson processes[END_REF] to handle multiple observed stochastic processes with irregular "jump" times. We then propose a general multi-observation HMM with parameters fitted (that define the transition rate matrices of the hidden and observed state processes) based on the standing assumption that the finite-state spaces of all stochastic processes are known. Our continuous-time finite-state HMM is then applied to our framework as a running example. Precisely, the HMM drives the grid and consists of a hidden Markov process (the reference configuration) and multiple observed Markov processes (the breakers' states). This approach may be potentially applicable in a system monitoring scheme. Some works focus on different assumed-types of malicious attacks to design system monitoring on electrical grids. For instance, replay attack is considered in [START_REF] Mo | Detecting integrity attacks on scada systems[END_REF][START_REF] Mo | Secure control against replay attacks[END_REF][START_REF] Mishra | Secure state estimation against sensor attacks in the presence of noise[END_REF], in which an attacker hijack sensors, observe and record the outputs. Denial-of-Service (DoS) attack models were assumed in [START_REF] Amin | Safe and secure networked control systems under denial-of-service attacks[END_REF][START_REF] Befekadu | Risk-sensitive control under markov modulated denial-of-service (dos) attack strategies[END_REF], in which an optimal control problem under security constraints is solved. False data injection was considered in [START_REF] Forti | Joint attack detection and secure state estimation of cyber-physical systems[END_REF][START_REF] Mo | False data injection attacks against state estimation in wireless sensor networks[END_REF][START_REF] Guo | Optimal linear cyber-attack on remote state estimation[END_REF][START_REF] Liu | False data injection attacks against state estimation in electric power grids[END_REF], in which some attacker can inject error measurements in the state estimation. For a study considering all these types of attacks, see, e.g., [START_REF] Shi | Causality countermeasures for anomaly detection in cyber-physical systems[END_REF]. However, it is difficult to know a priori the attack type on a system, and the assumptions made by these works may not be close to reality. In [START_REF] Shi | Event-based state estimation of linear dynamic systems with unknown exogenous inputs[END_REF][START_REF] Liu | False data injection attacks against state estimation in electric power grids[END_REF][START_REF] Fawzi | Secure state-estimation for dynamical systems under active adversaries[END_REF], it had been analyzed the case when an arbitrary error or unknown parameter is additively injected on the state model or measurement to represent, e.g., malicious attacks. However, in these works, it is common to assume that all the model parameters are known to estimate the system state. That is not the case in our problem because we need to estimate the parameters. While our HMM is later intended to be included in online monitoring algorithms to detect abnormal behaviors (e.g., as in [START_REF] Shi | On finite-state stochastic modeling and secure estimation of cyber-physical systems[END_REF]), the current paper mainly focuses on modeling and state and parameter estimation problems.

We provide an iteratively filter-based Expectation-Maximization (EM) approach [START_REF] Wonham | Some applications of stochastic differential equations to optimal nonlinear filtering[END_REF][START_REF] Zakai | On the optimal filtering of diffusion processes[END_REF][START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF][START_REF] Mclachlan | The EM algorithm and extensions[END_REF] to estimate the model parameters and the hidden state (e.g., the reference configuration). This approach aims to maximize a log-likelihood function over parameter space. While our parameter estimation method is close as obtained in previously cited papers and [START_REF] Elliott | Hidden Markov models: estimation and control[END_REF], we adapt the filter-based approach to handle multi-observed processes. Each one has its transition rate matrix (i.e., parameters), and no average is considered over them as in [START_REF] Korolkiewicz | A hidden markov model of credit quality[END_REF]. In detail, we suppose first that all state processes belong in a probability space representing the "real world". Then, we use a change of probability measure technique (Girsanov's Theorem, see, e.g., [START_REF] Elliott | Hidden Markov models: estimation and control[END_REF][START_REF] Brémaud | Point processes and queues: martingale dynamics[END_REF]) to define a new probability measure representing a "fictitious world". In this new space, filters for estimating the hidden state and the model's parameters are easy to obtain. They are linear Stochastic Differential Equations (SDEs) modulated by counting processes. Then, instead of using a classical Euler-Maruyama discretization (with small-time step) for all the SDEs obtained (e.g., as in [START_REF] Elliott | Discrete-time expectation maximization algorithms for Markov-modulated poisson processes[END_REF][START_REF] James | Time discretization of continuous-time filters and smoothers for HMM parameter estimation[END_REF][START_REF] Krishnamurthy | Finite-dimensional filters for passive tracking of Markov jump linear systems[END_REF]), we present a strong scheme with no discretization error for numerical purposes. This scheme adapts to all state change times generated by the temporal evolution of all observed processes.

The proposed modeling approach is finally confronted with available real data provided by France's transmission system operator (RTE), consisting of Boolean temporal sequences describing a set of breakers' states (off /on). The data has been collected during the normal behavior of the network over a given period.

The reference configurations are also obtained from this data using a version of the well-known K-means method [START_REF] Bock | Clustering methods: a history of K-means algorithms[END_REF]. We then identify the normal behavior of the French electrical grid, which could be embedded in a monitoring and detection algorithm in the future.

The rest of this paper is organized as follows. In Section 2, we introduce the general finite-state continuous-time HMM and dynamics. Section 3 provides a filtering approach to estimate the hidden state based on all observed processes. We also briefly recall the EM algorithm and compute a filter-based EM algorithm for all model parameters. In Section 4, we present a numerical method for all filter estimates. We also show how to obtain the initial estimation of all parameters. Section 5 shows the numerical results in a simulated scenario and the application of the breakers' states in the French electrical transmission grid, which illustrates a real-world application of the HMM considered in this paper. The final Section 6, gives concluding remarks.

The Modeling

Let us assume that the grid has finitely many reference configurations h 1 , ..., h N known a priori, N ∈ N * . Each reference configuration is a vector of length P ∈ N * (the number of breakers in the grid), and each of its components represents some information of a breaker (e.g., modes, states, etc.). Suppose that a state process H := {H t } t≥0 represents the evolution over time of the grid between the reference configurations h 1 , ..., h N , where H t denotes the unknown reference configuration at time t ≥ 0. Then, H constitutes a hidden process, i.e., it is not directly observable. Suppose that we observe each breaker state over time represented by the state process K p := {K p t } t≥0 . Then, K p t denotes the available information of the breaker p = 1, ..., P at time t ≥ 0 and constitutes the observable information of the grid. We aim in this paper to estimate the hidden evolution of the grid H in some optimal way based on the temporal evolution of K p , p = 1, ..., P .

First, we fix a complete probability space (Ω, F, P), where P is the probability measure of the "realworld", and we denote by E the expectation operator under P. We suppose that all state processes are continuous-time finite-state Markov Chains (MCs) defined on the common probability space (Ω, F, P). It is also assumed that (almost) all sample functions are right-continuous with left limits.

To be more general, we work with general hidden and observable processes. The case study considered of the electrical grid is a running example while explaining the general proposed approach.

The State Processes

Consider that the state space of the hidden process H = {H t } t≥0 is the finite set:

H := h 1 , h 2 , ..., h N ⊆ P p=1 M p , (1)
where, for each p = 1, ..., P ,

M p := m p 1 , m p 2 , ... , m p Mp ⊂ R R (2)
represents a finite set, with M p , R ∈ N * .

Example (Electrical grid). The sets H and M p can represent the set of known reference configurations of the grid and the set of different modes or states that the breaker p can eventually take, respectively (resp.).

For instance, if we consider M p = 2, then M p = m p 1 , m p 2 . Thus, if e.g., we assume R = 1, then M p represents the breaker states "off " and "on", i.e., we can let m p 1 = 0 and m p 2 = 1, resp.; and H is therefore a set of binary vectors. On the contrary, if we assume R = 2, then M p represents breaker modes with information in two dimensions. For example, a state mode is a vector in R 2 representing first the number of jumps between beaker's state values ("off " and "on") and second the time spent at each beaker's state value ("off " and "on"), both on average in one hour.

In order to simplify the modeling of the hidden process, we introduce a transformation over H. Let

1 n : H → {0, 1} the indicator function defined for each n = 1, ..., N by 1 n (h m ) = 1 if n = m, and
1 n (h m ) = 0 otherwise. Then, the vector (1 1 , 1 2 , ..., 1 N
) is a bijection from H to the set of unit vectors

X := {e 1 , e 2 , ..., e N } ⊂ R N
, where e n ∈ X denotes a vector in R N with unity in the n-th position and zero elsewhere, n = 1, ..., N . Thus, without loss of generality (w.l.o.g.), we shall consider an underlying state process X := {X t } t≥0 with state space X and defined by

X t := 1 1 (H t ) , 1 2 (H t ) , ... , 1 N (H t ) . (3)
Note that at any time t ≥ 0, just one component of X t is one and the others are all zero. In this way, X t can be represented as X t = N n=1 X t , e n e n , where

H t = N n=1 X t , e n h n . (4)
In this way, each state h n ∈ H is identified with an unit vector e n ∈ X, n = 1, ..., N . Thus, instead of inferring the hidden evolution of H, we can estimate w.l.o.g. the temporal evolution of the underlying state process X.

Similarly, the set of different state values of the observable process K p = {K p t } t≥0 , p = 1, ..., P , is defined by K p := 0, 1 , ..., K p -1 ⊂ N, where K p ∈ N * . Note in particular that we can consider in eq. ( 2), R = 1, M p = K p and M p = K p . For instance, if K p = 2, then K p takes the binary values 0 and 1. This set is identified with the set of unit vectors

Y p := f p 1 , f p 2 , ..., f p Kp ⊂ R Kp ,
where f p k ∈ Y p denotes the vector in R Kp with unity in the k-th position and zero elsewhere, k = 1, ..., K p .

Thus, w.l.o.g., the temporal evolution of the observable process K p can be represented by a state process

Y p := {Y p t } t≥0 . When Y p t is in the state f p k , that means that the observable process K p t is in state k -1 ∈ K p at time t ≥ 0.
Example (Electrical grid). If we consider K p t as the observable state process of the breaker p at time t ≥ 0, then Y p t being in the state f p k means that breaker p is in state k -1 ∈ K p at time t ≥ 0. If K p = 2 for each p = 1, ..., P , then all breakers take binary values (representing the values "off " and "on"), and

f p 1 = (1, 0)
and f p 2 = (0, 1) represent the state values 0 ∈ K p and 1 ∈ K p , resp. In addition, from eq. ( 2), if we consider R = 1, M p = K p = 2 and M p = K p , then M p is also the set of values "off " and "on" and H is a set of binary vectors.

The Dynamic of the Hidden Process

Since X = {X t } t≥0 is a Markov chain by assumption, we shall suppose that X has a transition rate

matrix A = (a ij ) ∈ R N ×N
, where

a ij := d dt P X t = e i | X 0 = e j t=0
represents the transition probability rate of X from state e j ∈ X to the state e i ∈ X, for each i, j = 1, ..., N , i = j. In addition, the transpose of A belongs to the Q-matrix class3 . Thus, defining p i,t := P[X t = e i ],

i = 1, ..., N , t ≥ 0, the probability distribution vector p t := (p 1,t , p 2,t , ..., p N,t ) satisfies the forward equation

dp t /dt = Ap t .
On the other hand, X is adapted to the (complete) right-continuous increasing family of the natural σ-fields generated by himself, i.e., to the natural filtration F t := σ(X s ; s ≤ t) ⊂ F. Then, the process [START_REF] Aggoun | Measure theory and filtering: Introduction and applications[END_REF]Lemma 2.6.18]. The semi-Martingale representation of X is therefore:

V t := X t -X 0 - t 0 AX s ds is a (F t , P)-Martingale
X t = X 0 + t 0 AX s ds + V t . (5)

The Dynamic of the Observed Processes

The observed state process

Y p = {Y p t } t≥0 is directly related with H = {H t } t≥0 . First, since X = {X t } t≥0
takes unit vectors in X ⊂ R N , then we can express any matrix of real-valued functions with finite range, let say C p : X → R Kp×Kp in function of X t as:

C p (X t ) = N n=1 C p (e n ) X t , e n . (6)
c p kl (e n ) = d dt P Y p t = f p k | Y p 0 = f p l , X 0 = e n t=0 = d dt P Y p t = f p k | Y p 0 = f p l , H 0 = h n t=0
by eq. ( 4), because X 0 , e n = 1 and X 0 , e n = 0, for each n = 1, ..., N , n = n

= d dt P Y p t = f p k | Y p 0 = f p l , H p 0 = m p m t=0
, where the last equality holds by the local assumption over the vectors h n ∈ H. In such a way, for each

m p m ∈ M p , m = 1, ..., M p ,
c p,m kl := d dt P Y p t = f p k | Y p 0 = f p l , H p 0 = m p m t=0 , (7)
for each k, l = 1, ..., K p , k = l. In addition, based on the eq. ( 6), we express C p (X t ) in function of the number of elements in M p , by:

C p (X t ) = Mp m=1 C p m n∈I p m X t , e n , (8)
where C p m ∈ R Kp×Kp is the matrix with components of the eq. ( 7), and I p m ⊆ I := {1, ..., N } is the subset of indices n = 1, ..., N for which the p-th position of the vector h n ∈ H is m p m ∈ M p , i.e.,

I p m := n ∈ I | proj p (h n ) = m p m , h n ∈ H , m p m ∈ M p .
For convenience, we associate with each I p m a diagonal matrix diag p m ∈ R N ×N defined by:

diag p m := diag 1 {1∈I p m } , 1 {2∈I p m } , ..., 1 {N ∈I p m } , (9)
where

1
Y p t = Y p 0 + t 0 C p (X s )Y p s ds + W p t . (10)
We denote by Y t := σ(Y p s ; s ≤ t, p = 1, ..., P ) the corresponding right-continuous complete filtration generated by all observed processes Y p , p = 1, ..., P .

Summary

In summary, H t represents a hidden MC (e.g., the temporal evolution of the electrical grid) that takes values in the known set H of eq. ( 1) (e.g., a set of reference configurations in the grid). Each element in H is a vector constructed from the sets M p of eq. ( 2) (e.g., the set of different modes or states of the breakers p = 1, ..., P ). Instead of estimating H t over time, we estimate w.l.o.g. the hidden underlying process X t , whose state space is X of unit vectors of R N . This quantities are related by the eq. ( 4). The definition of X t is given by eq. ( 3) and its semi-Martingale representation is given in eq. [START_REF] Mo | Secure control against replay attacks[END_REF]. The evolution at time t ≥ 0 of the observable process p = 1, ..., P is represented by K p t (e.g., the available information of the breaker p in the grid). This process takes numerical values in the set K p (e.g., the set of binary values representing the values "off" and "on" of breaker p). In the same way, we identify this set with the set Y p of unit vectors of R Kp and we work with the observable underlying process Y p t . This is the observable information to estimate X t , and then H t . The semi-Martingale representation of Y p t is given by eq. [START_REF] Mo | False data injection attacks against state estimation in wireless sensor networks[END_REF].

Example (Electrical grid). For instance, suppose that there are P = 4 breakers in the grid, the beaker's states are 0 and 1, and that the reference configurations are h 1 = (1, 0, 1, 0), h 2 = (0, 1, 0, 1), h 3 = (0, 0, 0, 0), and h 4 = (1, 1, 1, 0). Also suppose that the temporal evolution of each breaker state is as in Figure 1. Then, under these observations over time, we want to estimate the temporal evolution of the grid represented by H, but equivalently, using the process X as it is shown in Figure 2.

Finite-Dimensional Filter for the Temporal Evolution of the Grid

In this section, we provide a filtering approach to estimate the hidden process from all observed state processes. To do that, we also need filter estimates for the parameters involved in our model, i.e., the matrices A = (a ij ) and C p m = c p,m kl , for each p = 1, ..., P and m = 1, ..., M p , that define resp. the transition rates of X and Y p .

Moving to a "Fictitious World"

We are looking for an estimation of the hidden MC over time by using the filtration Y t = σ(Y p s ; s ≤ t, p = 1, ..., P ) generated by all observed processes Y p . This is done through the estimation of X t at each time t ≥ 0. The filtered estimate of X t under P is the expectation operator E over X t given Y t . In fact, since X t is defined as an indicator function in eq. ( 3), the filtered estimate is a conditional probability

distribution, i.e., E[X t | Y t ] = P[X t = e 1 | Y t ], ..., P[X t = e N | Y t ] .
It can be shown that an explicit equation for E[X t | Y t ] can be obtained, but it will be nonlinear. In contrast, by using some change of the probability measure P, we can obtain filtered estimate that will be is linear, as it will be shown below.

To

obtain E[X t | Y t ]
we can use a simple Bayes' rule.

Suppose that on the probability space (Ω, F, P) there is for each p = 1, ..., P a counting process N p kl,t of the number of jumps of the state process

Y p from state f p k ∈ Y p to state f p l ∈ Y p within the time interval [0, t], k, l = 1, ..., K p , k = l.
The semi-Martingale representation of N p kl,t can be obtained via the following decomposition:

N p kl,t = t 0 f p k , Y p s-f p l , dY p s = t 0 f p k , Y p s-f p l , C p (X s )Y p s ds + t 0 f p k , Y p s-f p l , dW p s , (11)
where, we have used the eq. ( 10) in differential form, and Y p t -:= lim s↑t Y p s is the left limit of the state process Y p t at t ≥ 0. Note that each N p kl,t is Y t -measurable for each t ≥ 0 and have no common jumps for indices (k , l ) = (k, l). Now, since C p (X t ) is given by ( 8), the semi-Martingale representation of N p kl,t is given by:

N p kl,t = t 0 f p k , Y p s- Mp m=1 f p l , C p m Y p s n∈I p m X s , e n ds + t 0 f p k , Y p s-f p l , dW p s = t 0 f p k , Y p s- Mp m=1 c p,m lk n∈I p m X s , e n ds + t 0 f p k , Y p s-f p l , dW p s = t 0 λ p kl,s ds + M p kl,t ,
where M p kl,t := N p kl,t -t 0 λ p kl,s ds is a (G t , P)-Martingale [START_REF] Aggoun | Finite-dimensional models for hidden Markov chains[END_REF], and λ p kl,t represents the "P-intensity" of the counting process N p kl,t , defined by:

λ p kl,t := f p k , Y p t- Mp m=1 c p,m lk n∈I p m X t , e n . (12)
Example (Electrical grid). Suppose that R = 1, M p = K p = 2 and M p = K p , then M p is a set of binary values. In this context, the process N p 12,t (resp. N p 21,t ) counts the number of jumps that Y p does from state

f p 1 = (1, 0) (resp., f p 2 = (0, 1)), to the state f p 2 = (0, 1) (resp., f p 1 = (1, 0)) in the time interval [0, t],
i.e., the representation of the number of jumps of the breaker p from state "off " to "on" (resp., from "on" to "off ").

From eq. ( 12), we find back the intuition that a higher transition rate c p,m 21 (resp.c p,m 12 ) between state "off "

to "on" (resp., from "on" to "off ") in the breaker p (which are entries in the matrix C p m ), is related to a higher intensity in the counting process N p 12,t (resp. N p 21,t ).

The idea is then to introduce a new probability measure P for a "fictitious world" from the probability measure P of the "real world" to change all intensities to one under P. This is described by means of the Radon-Nikodym derivative, see, e.g. [24, Ch. VI, Sec.2-3]. By using [24, Ch. VI, eq. (3.

3)] but for multidimensional case 4 , we define P by putting:

dP dP Gt = Λ t := exp        - P p=1 Kp k,l=1 k =l t 0 ln λ p kl,s dN p kl,s + P p=1 Kp k,l=1 k =l t 0 (λ p kl,s -1 ds        , (13)
which is a (G t , P)-martingale. Using now Ito's Lemma, see, e.g., [START_REF] Sennewald | Itô's lemma" and the bellman equation for poisson processes: An applied view[END_REF], we have:

Λ t = 1 - P p=1 Kp k,l=1 k =l t 0 Λ s-λ p kl,s -1 λ p kl,s -1 dN p kl,s -λ p kl,s ds . (14)
We also define the reverse counterpart of ( 13) by putting:

Λ t := exp        P p=1 Kp k,l=1 k =l t 0 ln λ p kl,s dN p kl,s - P p=1 Kp k,l=1 k =l t 0 (λ p kl,s -1 ds        , (15)
so that Λ t Λ t = 1. Again by Ito's Lemma, it holds:

Λ t = 1 + P p=1 Kp k,l=1 k =l t 0 Λ s-λ p kl,s -1 d N p kl,s -s . (16)
In this way, Λ t and N p kl,t -t are (G t , P)-martingale ∀ t ≥ 0. It can be also shown that, under P, the dynamic for X t is still given by ( 5

E[F t | Y t ] = E[Λ t F t | Y t ] E[Λ t | Y t ] , (17)
where E denotes the expectation operator under the probability measure P. We denote by

σ t (F t ) the expectation E[Λ t F t | Y t ]. Consequently σ t (1) = E[Λ t | Y t ].
Note that σ t (1) can be computed as the sum of the components of σ t (X t ). Indeed, since X t takes values in the space X of unit vectors of R N , then X t , 1 N = 1 for all t ≥ 0, where 1 N := N n=1 e n , and therefore

σ t (F t ) = σ t (F t X t , 1 N ) = σ t (F t X t ), 1 N .
Thus, in particular taking F t ≡ 1 we have σ t (1) = σ t (X t ), 1 N . The linear filtered estimate of X t is given in the next Proposition 3.1.

Proposition 3.1. The finite-dimensional (unnormalized) estimator for the states of X t is of the form:

σ t (X t ) = σ 0 (X 0 ) + A t 0 σ s (X s )ds - P p=1 Kp k,l=1 k =l t 0 σ s (X s-)d N p kl,s -s + P p=1 Kp k,l=1 k =l Mp m=1 t 0 f p k , Y p s-c p,m lk diag p m σ s (X s-)d N p kl,s -s , (18)
where diag p m is the diagonal matrix defined in eq. (9).

Proof.

The proof is postponed after that of the Theorem 3.3.

To obtain σ t (X t ), we need the estimation of all parameters involved in eq. ( 18), i.e., the matrices A = (a ij ) and C p m = c p,m kl , for each p = 1, ..., P and m = 1, ..., M p . This is the purpose of the next section.

Parameter Estimation

To estimate the parameters that define the transition rate matrices of X and Y p , p = 1, ..., P , we focus on the EM algorithm for continuous-time stochastic processes, see, e.g., [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF][START_REF] Campillo | MLE for partially observed diffusions: direct maximization vs. the EM algorithm[END_REF][START_REF] Dembo | Parameter estimation of partially observed continuous time stochastic processes via the EM algorithm[END_REF]. The idea is to maximize a likelihood function in an iterative form. Let {P θ , θ ∈ Θ} be a family of probability measures on the measurable space (Ω, F), all absolutely continuous with respect to the (initial) fixed probability measure P, wherein our case,

Θ := a ij , c p,m kl ; 1 ≤ i, j ≤ N , i = j , 1 ≤ k, l ≤ K p , k = l , 1 ≤ m ≤ M p , 1 ≤ p ≤ P . (19)
The log-likelihood for an estimation of a θ ∈ Θ can be defined by:

L(θ) := ln E dP θ dP Y ,
where Y ⊂ F, and then, the Maximum Likelihood Estimator (MLE) is defined by θ * ∈ arg max θ∈Θ L(θ).

In general, computing directly the MLE is challenging. The Expectation-Maximization (EM) algorithm provides an iterative approximation method starting from an initial estimation θ 0 , see Section 4.4. This algorithm is based on the following straightforward application of the well-known Jensen's inequality:

L(θ) -L( θ ) = ln E θ dP θ dP θ Y ≥ E θ ln dP θ dP θ Y =: Q θ, θ .
This gives a global minoration for the log-likelihood mapping θ → L(θ) by means of the auxiliary mapping θ → L( θ ) + Q θ, θ . At each iteration r ∈ N 0 , the EM algorithm consists of two main steps:

(1) E-step: set θ = θ r and compute Q • , θ ,

(2) M-step:

find θ r+1 ∈ arg max θ∈Θ Q θ, θ .
This algorithm can be stopped when a stopping test is satisfied, see Section 4.2. The generated sequence θ r r∈N0 gives nondecreasing values of the likelihood function, i.e., L θ r+1 > L θ r unless θ r+1 = θ r . For convergence issues, see, e.g., [START_REF] Campillo | MLE for partially observed diffusions: direct maximization vs. the EM algorithm[END_REF][START_REF] Dembo | Parameter estimation of partially observed continuous time stochastic processes via the EM algorithm[END_REF][START_REF] Wu | On the convergence properties of the EM algorithm[END_REF].

In our context, suppose our model is determined by some parameters θ ∈ Θ, i.e., we have computed already the E-step under θ. To compute the new parameters θ ∈ Θ that maximize the log-likelihood, i.e., the M-step, we have the following Theorem 3.2.

a ji = E[J ij,t | Y t ] E[O i,t | Y t ] , c p,m lk = n∈I p m E[L p,n kl,t | Y t ] n∈I p m E[S p,n k,t | Y t ]
,

where, for i, j = 1, ..., N , i = j,

J ij,t := t 0 e i , X s-e j , dX s , O i,t := t 0 e i , X s ds , (20)
and, for p = 1, ..., P , k, l = 1, ..., K p , k = l, and n = 1, ..., N , From eq. ( 17), the estimation a ij and c p,m kl can be obtained via the probability measure P by:

L p,
a ji = σ t (J ij,t X t ), 1 N σ t (O i,t X t ), 1 N , c p,m lk = n∈I p m σ t L p,n kl,t X t , 1 N n∈I p m σ t S p,n k,t X t , 1 N . ( 22
)
In this way, it is sufficient to compute the estimators σ t (J ij,t X t ), σ t (O i,t X t ), σ t S p,n k,t X t and σ t L p,n kl,t X t . Now, if we consider the process:

F t = F 0 + t 0 α(X s )ds + t 0 β(X s ), dV s (23)
then, for each i, j = 1, ..., N , i = j, p = 1, ..., P , k = 1, ..., K p , n = 1, ..., N , the processes J ij,t , O i,t and S p,n k,t are considered into F t , where F 0 ∈ R is known, and α : X → R and β : X → R N are known functions with finite range, G t -adapted and integrable for each t ≥ 0. Indeed, by using eq. ( 5), and taking

(i) F 0 = 0 ∈ R, α(X t ) = e i
, X t a ji , and β(X t ) = e i , X t e j , we obtain

F t = J ij,t , (ii) F 0 = 0 ∈ R, α(X t ) = e i , X t , and β(X t ) = 0 N ∈ R N , we obtain F t = O i,t , (iii) F 0 = 0 ∈ R, α(X t ) = f p k , Y p t e n , X t , and β(X t ) = 0 N ∈ R N , we obtain F t = S p,n k,t .
Therefore, to compute σ t (J ij,t X t ), σ t (O i,t X t ), and σ t S p,n k,t X t , we can compute once σ t (F t X t ) and restrict afterwards to the particular cases of α(X t ) and β(X t ). On the other hand, we know that σ t (F t ) = σ t (F t X t ), 1 N , so that, we make the inner product between σ t (F t X t ) and 1 N to have the estimation for all parameters, see eq. ( 22). The following Theorem 3.3 gives the linear filter estimate σ t (F t X t ). The filter estimate σ t L p,n kl,t X t is given in Theorem 3.6.

Theorem 3.3. The finite-dimensional (unnormalized) estimator for F t X t is of the form:

σ t (F t X t ) = σ 0 (F 0 X 0 ) + A t 0 σ s (F s X s )ds - P p=1 Kp k,l=1 k =l t 0 σ s (F s-X s-)d N p kl,s -s + t 0 σ s (X s α s )ds + N i,j=1 i =j t 0 σ s (β j,s -β i,s )X s , e i a ji (e j -e i )ds (24)
+ P p=1 Kp k,l=1 k =l Mp m=1 t 0 f p k , Y p s-c p,m lk diag p m σ s (F s-X s-)d N p kl,s -s ,
where α t := α(X t ) ∈ R and β t := β(X t ) ∈ R N are known functions with finite range, G t -adapted and integrable for each t ≥ 0, F t is given in eq. ( 23), and diag p m is the diagonal matrix of eq. (9).

Proof.

See Proof 2 in Appendix A.

Remark 3.4. Note that if we consider F t = F 0 = 1, α t = 0 ∈ R, and β t = 0 N ∈ R N in the eq. ( 23), then the Proposition 3.1 is a particular case of Theorem 3.3.

In this way, the filter estimates for the parameter estimation are given in the next Corollary 3.5 by taking the particular cases of α t and β t within F t , see eq. ( 23) and the cases (i), (ii) and (iii).

Corollary 3.5. The finite-dimensional (unnormalized) estimator for J ij,t X t , O i,t X t and S p,n k,t X t are resp.

of the form:

σ t (J ij,t X t ) =A t 0 σ s (J ij,s X s )ds + t 0 σ s (X s
), e i e j a ji ds -

P p=1 Kp k,l=1 k =l t 0 σ s (J ij,s-X s-)d N p kl,s -s + P p=1 Kp k,l=1 k =l Mp m=1 t 0 f p k , Y p s-c p,m lk diag p m σ s (J ij,s-X s-)d N p kl,s -s , σ t (O i,t X t ) = A t 0 σ s (O i,s X s )ds + t 0 σ s (X s ), e i e i ds - P p=1 Kp k,l=1 k =l t 0 σ s (O i,s-X s-)d N p kl,s -s + P p=1 Kp k,l=1 k =l Mp m=1 t 0 f p k , Y p s-c p,m lk diag p m σ s (O i,s-X s-)d N p kl,s -s , σ t (S p,n k,t X t ) = A t 0 σ s (S p,n k,s X s )ds + t 0 f p k , Y p s-σ s (X s ), e n e n ds - P q=1 Kq u,v=1 u =v t 0 σ s (S p,n k,s-X s-)d N q uv,s -s + P q=1 Kq u,v=1 u =v Mq m=1 t 0 f q u , Y q s-c q,m vu diag q m σ s (S p,n k,s-X s-)d N q uv,s -s ,
where diag p m is the matrix of eq. ( 9).

Finally, the filter estimate we need to complete the estimation of all parameters is σ t L p,n kl,t X t . For this, we write the semi-Martingale representation of L p,n kl,t from eq. ( 21), k = l, by:

L p,n kl,t = t 0 e n , X s-d N p kl,s -s + t 0 e n , X s-ds . ( 25
)
Theorem 3.6. The finite-dimensional (unnormalized) estimator for L p,n kl,t X t , k = l, is of the form:

σ t L p,n kl,t X t = A t 0 σ s L p,n kl,s X s ds - P q=1 Kq u,v=1 u =v t 0 σ s L p,n kl,s-X s-d N q uv,s -s + Mp m=1 t 0 f p k , Y p s-c p,m lk diag p m e n , σ s (X s-) e n dN p kl,s (26)
+ P q=1 Kq u,v=1 u =v Mq m=1 t 0 f q u , Y q s-c q,m vu diag q m σ s L p,n kl,s-X s-d N q uv,s -s ,
where diag p m is the diagonal matrix defined in eq. (9).

Proof.

See Proof 3 in Appendix A.

Numerical Methods

This section presents a general numerical method for all SDEs involved in our model, i.e., a general scheme involving each SDE of the filter estimates σ t (X t ), σ t (J ij,t X t ), σ t (O i,t X t ), σ t S p,n k,t X t and σ t L p,n kl,t X t , obtained resp. in Proposition 3.1 Corollary 3.5 and Theorem 3.6. We also show how to obtain the initial estimation of all parameters and the set H when M p = K p for each p = 1, ..., P . This is a particular case that applies in the case study considered in this paper of the electrical grid by using available data of all breakers states.

First, we write t 0 < t 1 < ... < t W +1 the increasingly ordered instances of all "jump times" of the observed processes. This is shown in Figure 3

Jump-Adapted Scheme for Filters

The scheme presented here is a strong approximation with no discretization error on the SDEs solution, in which a jump-adapted scheme is given by a superposition of all "jump times", i.e., all state change times generated by the temporal evolution of the observed processes, see Figure 3. The jump effects are then added at the correct jump times. To use this kind of method, one has to check whether the SDE concerned belongs to the particular subclass of SDEs for which the corresponding non-jump part admits an exact solution, see, e.g., [33, Ch. II]. In our case, all SDEs admit an explicit solution in the non-jump parts as we see in next.

Instead of building a scheme for each filter estimate σ t (X t ), σ t (J ij,t X t ), σ t (O i,t X t ), σ t S p,n k,t X t and σ t L p,n kl,t X t , we present a generalized scheme for the following SDEs system:

d K t = Ξ t K t dt + P p=1 Kp k,l=1 k =l Q p kl,t-K t-dN p kl,t , (27)
d G t = Υ t G t + Γ t K t dt + P p=1 Kp k,l=1 k =l Π p kl,t-G t-+ Λ p kl,t-K t-dN p kl,t , (28)
where K t , G t ∈ R N are G t -adapted and integrable for any t ≥ 0, and

Ξ t , Q p kl,t , Υ t , Γ t , Π p kl,t , Λ p kl,t ∈ R N ×N
are constant matrices between "jumps" for each k, l = 1, ..., K p , p = 1, ..., P . They are defined by Ξ t = Ξ tw , Υ t = Υ tw , Γ t = Γ tw , Q p kl,t = Q p kl,tw , Π p kl,t = Π p kl,tw , Λ p kl,t = Λ p kl,tw , for any t ∈ [t w , t w+1 ), w = 0, ..., W .

Theorem 4.1. The jump-adapted exact solution scheme for K t and G t with initial conditions K 0 , G 0 ∈ R N is written by: As a consequence, the jump-adapted scheme for each filter estimate σ t (X t ), σ t (J ij,t X t ), σ t (O i,t X t ), σ t S p,n k,t X t and σ t L p,n kl,t X t can be easily obtained from the Theorem 4.1 as follows. First, to relax the notation of all equations, we define, for each k, l = 1, ..., K p , p = 1, ..., P , the matrices Φ t and Ψ p kl,t for each t ≥ 0, by:

K t0 = K 0 , G t0 = G 0 , K t - w+1 = exp Ξ tw ∆t w+1 K tw , G t - w+1 = Id N 0 N ×N exp Υ tw Γ tw 0 N ×N Ξ tw ∆t w+1 G tw K tw , K tw+1 = K t - w+1 + P p=1 Kp k,l=1 k =l Q p kl,t - w+1 K t - w+1 ∆N p kl,tw+1 , G tw+1 = G t -
Φ t := A - P p=1 Kp k,l=1 k =l Ψ p kl,t , (29)
Ψ p kl,t := Mp m=1 f p k , Y p t c p,m lk diag p m -Id N , (30)
where Id N is the identity matrix in R N ×N , and diag p m is the diagonal matrix defined in eq. ( 9). Note that Φ t and Ψ p kl,t are constant between "jumps". That is, Φ t = Φ tw and Ψ p kl,t = Ψ p kl,tw for any t ∈ [t w , t w+1 ), w = 0, . . . , W .

Note that we can obtain the SDE of σ t (X t ) (see eq. ( 18)) from eq. ( 27) by considering K t = σ t (X t ), Ξ t = Φ t , and for each p = 1, ..., P , k, l = 1, ..., K p , Q p kl,t = Ψ p kl,t . Thus, from Theorem 4.1, the following scheme holds for the filter estimate σ t (X t ) of the grid state.

Corollary 4.2. A jump-adapted exact solution scheme for σ t (X t ) is written by: For the other filter estimates, σ t (J ij,t X t ), σ t (O i,t X t ), and σ t S p,n k,t X t , we can express the scheme for σ t (F t X t ) from the Theorem 3.3, see the cases (i), (ii) and (iii) after eq. ( 23). The following scheme holds for the filter estimate σ t (F t X t ).

σ t0 (X t0 ) = σ 0 (X 0 ) , σ tw+1 X t - w+1 = exp Φ tw ∆t w+1 σ tw (X tw ) , σ tw+1 X tw+1 = σ tw+1 X t -
Corollary 4.3. The jump-adapted exact solution scheme for σ t (F t X t ) with initial condition σ 0 (F 0 X 0 ) ∈ R N is written by: β(e i ), e j -e i a ji (e j -e i ) e i , where diag α(e 1 ), ..., α(e N ) is a diagonal matrix in R N ×N with diagonal (α(e 1 ), ..., α(e N )), and α : X → R and β : X → R N are the known functions in F t , see eq. [START_REF] Korolkiewicz | A hidden markov model of credit quality[END_REF].

σ t0 (F t0 X t0 ) = σ 0 (F 0 X 0 ) , σ tw+1 F t - w+1 X t - w+1 = Id N 0 N ×N exp Φ tw Γ tw 0 N ×N Φ tw ∆t w+1 σ tw (F tw X tw ) σ tw (X tw ) , σ tw+1 (F tw+1 X tw+1 ) = σ tw+1 F t -

Proof.

See Proof 5 in Appendix A. Now, for each u, v = 1, .., K q , u = v, q = 1, ..., P , and n = 1, ..., N , the SDE of σ t L q,n uv,t X t (see eq. ( 26)) can be obtained from eq. ( 28) by considering

G t = σ t L q,n uv,t X t , K t = σ t (X t ), Υ t = Ξ t = Φ t , Γ t = 0 N ×N ∈ R N ×N
, and for each p = 1, ..., P , k, l = 1, ..., K p , Π p kl,t = Q p kl,t = Ψ p kl,t and Λ p kl,t = Ψ q uv,t + Id N diag(e n ) for k = u, l = v, p = q, and Λ p kl,t = 0 N ×N ∈ R N ×N otherwise, where diag(e n ) is a diagonal matrix in R N ×N with diagonal e n ∈ X. Thus, from Theorem 4.1, the following scheme holds for the filter estimate

σ t L q,n uv,t X t .
Corollary 4.4. A jump-adapted exact solution scheme for σ t L q,n uv,t X t , u, v = 1, .., K q , u = v, q = 1, ..., P , and n = 1, ..., N , is written by:

σ t0 L q,n uv,t0 X t0 = 0 N , σ tw+1 L q,n uv,t - w+1 X t - w+1 = Id N 0 N ×N exp Φ tw 0 N ×N 0 N ×N Φ tw ∆t w+1 σ tw L q,n uv,tw X tw σ tw (X tw ) , σ tw+1 L q,n uv,tw+1 X tw+1 = σ tw+1 L q,n uv,t - w+1 X t - w+1 + P p=1 Kp k,l=1 k =l Ψ p kl,t - w+1 σ tw+1 L q,n uv,t - w+1 X t - w+1 ∆N p kl,tw+1 + Ψ q uv,t - w+1 + Id N diag(e n ) σ tw+1 X t - w+1 ∆N q uv,

An Stopping Criteria for the EM Algorithm

Instead of using the strict stopping criteria for the EM algorithm5 θ r+1 = θ r for some r ∈ N, we define the following stopping test for numerical purposes:

σ (r) t (X t ) -σ (r-1) t (X t ) + A (r) -A (r-1) + P p=1 Mp m=1 C p,(r) m -C p,(r-1) m σ (r-1) t (X t ) + A (r-1) + P p=1 Mp m=1 C p,(r-1) m ≤ ε , (31)
where ε > 0 is a given stopping parameter. We note that the normalization by the sum of parameter norms limits the dependency of ε to the magnitude of the parameters.

Uncertainty and Final Hidden State

Note that in Proposition 3.1 we obtain the estimator σ t (X t ) = E[Λ t X t | Y t ] and then, by applying the Bayes' rule of eq. ( 17), we obtain the filtered estimate

σ t (X t ) = E[X t | Y t ] which is the probability distribution P[X t = e 1 | Y t ], ..., P[X t = e N | Y t ] .
In order to know the exact estate of the hidden process over the states h 1 , ..., h N , we then take for each time t ≥ 0, the value:

X t = e n , n ∈ arg max n =1,...,N σ t (X t ), e n . (32)
Considering this choice, at fixed time t ≥ 0, X t , e n = 1 and X t , e n = 0 for each n = 1, ..., N , n = n, and by eq. ( 4), the estimated hidden state is therefore h n at such time t ≥ 0.

Concerning the uncertainty of the choice in eq. ( 32), we can use the estimator σ t (X t ) for each t ≥ 0 in the following way. First, for each t ≥ 0, let

t := 1 -max n=1,...,N σ t (X t ), e n (33)
be the function that represents how far the highest σ t (X t ) is from the value one. Recall that this estimator is a probability distribution. So, for a fixed t ≥ 0, if one component of σ t (X t ) is near to one, then selecting the estimated hidden state by using the eq. ( 32), is an almost-sure choice. Thus, for any time t ≥ 0, the function t represents the uncertainty on the hidden states. For the Section 5 of numerical results, we also define the mean

δ t := 1 t t 0 s ds (34)
to see how far is t of the mean δ t , in particular at the jump times. In this section, we consider the particular case when M p = K p , R = 1, M p = K p for each p = 1, ..., P .

For the EM algorithm, a θ 0 ∈ Θ must be initialized, i.e., we need to choose initial values for the matrices At first glance, we can use Theorem 3.2 empirically, i.e., we can discretize all the involved integrals in the parameter estimation (those in [START_REF] Zakai | On the optimal filtering of diffusion processes[END_REF] and eq. ( 21)), and use all the information of the state processes to estimate θ 0 .

A (0) = ( a ij ( 0
First, since we know the change of each observed process, then we know the values of each Y p t , p = 1, ..., P , at any time t ≥ 0, since it remains constant between "jumps". Second, because it is assumed that we know a priori the space H, then we can compute an empirical estimation of H, and so that, an empirical estimation of the underlying state process X t by eq. ( 3) and eq. ( 4). This can be computed through a distance measure by finding the closest state to the information vector of all observed processes at each jump time, i.e., by clustering and classification method. We show that in the following.

The joint state of all observed processes can be represented by a piecewise constant state process Y :=

{Y t } t≥0
, where Y t is the joint information of observed processes at time t ≥ 0, defined by:

Y t := P p=1 Kp k=1 f p k , Y p t (k -1)g p ∈ P p=1 K p , (35)
where for each p = 1, ..., P , k = 1, ..., K p , the k-th unit vector f p k ∈ K p , and g p ∈ R P is the p-th unit vector of R P . Note that Y t is piecewise right-continuous with left limits. In this way, if a state process Y p changes of value at time t ≥ 0, then Y t changes too. This occurs at the times t 0 < t 1 < ... < t W , see Figure 3.

Let d : R P × R P → R + a distance measure. Under the knowledge of the set H = {h 1 , ..., h N }, we can compute an empirical estimation H := { H t } t≥0 of the hidden state over h 1 , ..., h N , by:

H t ∈ arg min hn∈H n=1,...,N d(Y t , h n ) . ( 36
)
This approach is a classification procedure over the joint information of all observable processes at time t ≥ 0, which is represented by Y t in eq. ( 35). Now, with this empirical estimation H t , we can compute the empirical values of X t by means of its definition in eq. ( 3). We denote this estimation by X := { X t } t≥0 . Note that H and X are also piecewise right-continuous with left limits. Since we know now the values of X t at the state change times of Y , t 0 , ..., t W , we discretize all the involved integrals in Theorem 3.2 to estimate empirically all the parameters of our model. For each i, j = 1, ..., N , i = j, the initial estimation of A = (a ij ), i.e., the matrix

A (0) = a ij (
a ij (0) := W -1 w=0 e i , X tw e j , X tw+1 W -1 w=0 e i , X tw ∆t w+1 , c p,m lk (0) := W -1 w=0 f p k , Y p tw f p l , Y p tw+1 n∈I p m e n , X tw W -1 w=0 f p k , Y p tw n∈I p m e n , X tw ∆t w+1 . ( 37
)
On the other hand, the initial estimation σ 0 (X 0 ) of the process X at time t = 0 that we need in eq. ( 18), is given by the empirical estimation σ 0 ( X 0 ) defined for each n = 1, ..., N , by:

σ 0 ( X 0 ), e n := W -1 w=0 e n , X tw ∆t w+1 N n=1 W -1 w=0 e n , X tw ∆t w+1 . ( 38
)

A Method to Find the Reference Configurations of the Grid

We want to build the set H by using the available data which consists of the temporal evolution of the breakers states in the grid, by using Y = {Y t } t≥0 of eq. ( 35). Here, we also consider the particular case when M p = K p , R = 1, M p = K p , for each p = 1, ..., P .

To obtain the reference configurations of the grid, we construct clusters from data by partitioning it into N subsets. Each subset U 1 , ..., U N , called cluster, is represented by its representative state µ 1 , ..., µ N , resp.

To obtain optimal clusters, we use K-means method, see, e.g., [START_REF] Bock | Clustering methods: a history of K-means algorithms[END_REF]. The extension of K-means in continuous time is given by the minimization of the following cost function in the horizon time T > 0:

J = N n=1 T 0 1 {Yt∈Un} d(Y t , µ n )dt ,
where 1 {Yt∈Un} is the indicator function for sets, i.e., 1 {Yt∈Un} = 1 if Y t ∈ U n , and 1 {Yt∈Un} = 0 otherwise.

Since Y = {Y t } t≥0 is piecewise right-continuous on 0 = t 0 < t 1 < ... < t W = T , see eq. ( 35), we have

J = N n=1 W -1 w=0
1 {Yt w ∈Un} d(Y tw , µ n )∆t w+1 ,
where ∆t w+1 = t w+1 -t w . This corresponds to the classical discrete K-means approach with a weighted cost. Here, the representative state of a cluster U n , n = 1, ..., N , is:

µ n ∈ arg min η∈Un W -1 w=0
1 {Yt w ∈Un} d(Y tw , η) .
We use the traditional approach to construct the clusters by classification. For each t ≥ 0, the cluster U n , n = 1, ..., N , is obtained as the set:

U n = W -1 w=0 Y tw ∈ P p=1 K p d(Y tw , µ n ) ≤ d(Y tw , µ n ) , for each n = 1, ..., N . (39)
Note that the empirical estimation of the grid states in eq. ( 36) also means that for any time t ≥ 0, there is

a n ∈ {1, ..., N } such that H t ∈ U n .
This method generates a sequence {µ r } r∈N0 , where µ r := µ is the vector of all representative states of the clusters at iteration r ∈ N 0 . This can be initialized from the available data, e.g., randomly, heuristically, or by K-means++ approach [START_REF] Arthur | K-means++: The advantages of careful seeding[END_REF]. It should be noted that the performance of an iterative clustering algorithm may converges to numerous local minima and depends highly on initial cluster centers [START_REF] Khan | Cluster center initialization algorithm for K-means clustering[END_REF].

Finally, the reference configurations of the grid are given at the last iteration of the method, that is when µ r+1 = µ r for some r ∈ N 0 . In such a way, the reference configuration h n is µ

n for each n = 1, ..., N .

Numerical Results

In this section, we present our model's numerical results. Because we focus on the application for an electrical grid, we consider Boolean temporal sequences describing the breakers' states (off /on) of the network. First, we evaluate a simulated scenario in which the Markov state process, that represents the hidden process, is known. Recall that only the observed processes are used to infer the hidden process. After the temporal evolution of the hidden process is estimated, we compare it with the "real" state process for validation. Second, the proposed modeling is confronted with real data provided by the France's transmission system operator (RTE). In both cases, we fix the states for the observable processes as binary values, i.e., M p = K p = {0, 1}, R = 1, and M p = K p = 2 for each p = 1, ..., P , see eq. ( 2). In this way, the underlying space of the observed states is Y p = {f p 1 , f p 2 }, where f p 1 = (1, 0) and f p 2 = (0, 1).

Simulated Data

The simulated scenario is performed by the initial parameters shown in next.

Initial Parameters

We fix first the space of N = 4 hidden states, belonging to the space H considered here as:

H = (1, 1, 1, 1, 0, 1), (0, 1, 0, 0, 0, 1), (1, 0, 1, 1, 1, 1)
, (0, 1, 0, 1, 0, 0) .

Thus, the space X = {e 1 , e 2 , e 3 , e 4 } stands for the space of canonical vectors in R 4 , where e 1 = (1, 0, 0, 0), e 2 = (0, 1, 0, 0), e 3 = (0, 0, 1, 0), e 4 = (0, 0, 0, 1). The simulation of X = {X t } t≥0 and Y p = {Y 38). Under the simulated data, we obtain σ 0 ( X 0 ) = (0.307, 0.196, 0.388, 0.109). The initial state for X is therefore chosen to be X 0 = e 3 by eq. ( 32). Thus, by eq. ( 4), the initial hidden state is h 3 = (1, 0, 1, 1, 1, 1). Finally, for the stopping criteria of eq. ( 31), we choose ε = 10 -5 .

Hidden State Estimation by Clustering Method

On the left side in Figure 5, we observe the "real" temporal evolution of the hidden Markov process X.

The first estimation that we do is the empirical estimation by clustering method by using eq. ( 36). This approach is a classification procedure that takes the joint information of all observed processes at each time and computes the arg min set to know at which cluster the hidden process is. A cluster is obtained by finding the points of the observed processes' joint information that have the minimum distance to a cluster center, see eq. [START_REF] Wong | Stochastic processes in engineering systems[END_REF], where the centers are the states of eq. ( 40). The metric distance that we use is the euclidean one. The purpose of using this method is only for comparison. This paper does not aim to find the best clustering method with the appropriate distance measure. A study about this subject will be addressed in the future, based, e.g., on [START_REF] Choi | A survey of binary similarity and distance measures[END_REF].

Under the clustering method with the euclidean distance, the arg min set might not be a single point several times because there are points of the joint information of all observable processes at the same distance of different vectors of H. In such cases, the clustering method is not exact because we cannot know which state the hidden process is. This is represented with green points on the left side in Figure 5. When a cluster is active, i.e., when the hidden state can be chosen by the clustering method, it is represented on the right side in Figure 5. This is shown with a red point in the same picture when there is more than one choice.

Hidden State Estimation by HMM

Figure 6 shows the filter estimate of the hidden state over the set H of eq. ( 40). The hidden process X represents this over the canonical vectors of R 4 . On the left side in Figure 6, the filtering is done using the (empirical) initial parameter estimation of eq. [START_REF] Choi | A survey of binary similarity and distance measures[END_REF]. At first glance, the filter estimate "jumps" several times when the "real" temporal evolution of the hidden process does not. This is because the first parameter estimation is not exact. However, when the parameter estimation is computed through the Theorem 3.2, the filter estimate fits better as the number of iterations increases. This is confirmed by the Mean Squared Error (MSE) between the "real" values of X, and the filter estimate X. The MSE value obtained at the last iteration is less then 0.09, showing the accuracy of our model over the simulated scenario. The filter estimate for the last iteration is shown in Figure 6 on the right side.

Parameter Estimation

Figure 7 shows the values of the parameters of our model over the number of iterations, i.e., the values of the matrices A (r) and C p(r) m , p = 1, ..., 6, m = 1, 2, for each iteration r = 0, ..., 14. The first estimation is made by eq. ( 37) and the last one is obtained when the stopping criteria of eq. ( 31) is verified with ε = 10 -5 .

Uncertainty and Final Hidden State

Using the filter estimate of the hidden state, we can compute the function t of eq. ( 33), that represents a temporal evolution's uncertainty signal on the vectors in H. Figure 8 shows the values of this uncertainty over time for the first and the last iterations. Comparing the two pictures, we observe that there are fewer peaks at the last iteration compared to the first iteration. This is because finding the optimal parameters to fit our model, the uncertainty decreases as the number of iterations increases. In the Figure 8, we also see the mean δ t of eq. [START_REF] Arthur | K-means++: The advantages of careful seeding[END_REF]. At the last iteration when the parameters are fitted, less than 1% of the time, t is over δ t . When t is over δ t could mean that the state at which the hidden process "jumps" is known with some uncertainty. However, such uncertainty is almost instantaneous because it is not remaining in time, as reflected in Figure 8.

Conclusion and Remarks

In this paper, we have proposed a general data-driven approach for the hidden MC with several observed process. While the application is based on the breakers' states in an electrical transmission grid, we believe the model is general enough to serve other types of dynamic networks. Our framework was based on a continuous-time finite-state Hidden Markov Model (HMM) driven by multiple-observed counting processes.

The central assumption in the application was that the grid's state varies around a finite set of reference configurations. The grid's current reference configuration is unknown and constitutes the hidden state, while each breakers' state is an observable process. We have provided a filter-based expectation-maximization approach using a change of probability measure method to estimate recursively the model parameters and the hidden reference configuration of the grid. Filter estimates are also given for various processes related to the Markov state processes.

Further, we have presented a strong scheme with no discretization error for a general filter dynamic for numerical purposes. The state change effects of the breakers are then added at the correct "jump" times.

In addition, a clustering approach was also presented to identify the set of reference configurations of the grid. A future work will be to identify the "best" number of hidden states, based for example on [START_REF] Pohle | Selecting the number of states in hidden markov models: pragmatic solutions illustrated using animal movement[END_REF].

Using our theoretical results, we have then shown the performance of the framework by considering Boolean temporal sequences describing the breakers' states (off /on) of the grid. First, we have evaluated a simulated scenario, showing the advantages of the HMM approach with the proposed strong scheme.

Second, we have confronted our model with real data provided by the France's transmission system operator (RTE), showing stability in the filter estimates for the parameters and the hidden state of the grid. We finally identify the normal behavior of the French electrical grid, which could be embedded in monitoring and detection algorithm in the future.

Appendix A. Mathematical Proofs

Proof 1 [Theorem 3.2]. The idea is to change the parameters that define the "intensities" of the counting processes J ij,t and N p kl,t of X and Y p resp., i.e., to modify a ij and c p,m kl to a ij and c p,m kl resp. Next, we proceed to the maximization step. More precisely, for each p = 1, ..., P , the parameters c p,m kl , m = 1, ..., M p , define λ p kl,t of the process N p kl,t , see eq. ( 12), and a ij appears in J ij,t since its semi-Martingale decomposition is:

J ij,t = t 0 e i , X s-e j , dX s = e i , X s-a ji dt + Q ij,t ,
from eq. ( 5), where Q ij,t := J ij,t -t 0 e i , X s-a ji ds is a (G t , P)-Martingale [START_REF] Aggoun | Finite-dimensional models for hidden Markov chains[END_REF].

To estimate all the parameters, we define a new probability measure6 P = P θ , with θ ∈ Θ, see eq. ( 19), for a "fictitious world" from the probability measure P = P θ , with θ ∈ Θ, of the "real world", by putting:

d P dP Gt = Λ t :=
F t X t = F 0 X 0 + t 0 X s-dF s + t 0 F s-dX s + [F, X] t = F 0 X 0 + t 0 X s-α s ds + t 0 X s-β s , dV s + A t 0 F s-X s ds + t 0 F s-dV s + 0<s≤t
β s , ∆X s ∆X s by using eqs. ( 5) and ( 23), and since

[F, X] t = 0<s≤t ∆F s ∆X s = 0<s≤t β s , ∆X s ∆X s = F 0 X 0 + t 0 X s-α s ds + t 0 X s-β s , dV s + A t 0 F s-X s ds + t 0 F s-dV s (A.3) + N i,j=1 i =j t 0 (β j,s -β i,s )X s , e i a ji (e j -e i )ds + N i,j=1 i =j t 0 (β j,s -β i,s )X s-, e i e j , dV s (e j -e i ) since 0<s≤t β s , ∆X s ∆X s = N i,j=1 i =j t 0
(β j,s -β i,s )X s-, e i e j , dX s (e j -e i ) and using eq. ( 5) .

Taking now Ito's product rule on Λ t F t X t , we have that:

Λ t F t X t = Λ 0 F 0 X 0 + t 0 F s-X s-dΛ s + t 0 Λ s-dF s X s + [Λ, F X] t = F 0 X 0 + P p=1 Kp k,l=1 k =l t 0 Λ s-F s-X s-λ p kl,s-d N p kl,s -s - P p=1 Kp k,l=1 k =l t 0 Λ s-F s-X s-d N p kl,s -s + t 0 Λ s-X s-α s ds + t 0 Λ s-X s-β s , dV s + A t 0 Λ s-F s-X s ds + t 0 Λ s-F s-dV s + N i,j=1 i =j t 0 Λ s-(β j,s -β i,s )X s , e i a ji (e j -e i )ds + N i,j=1 i =j t 0 Λ s-(β j,s -β i,s )X s-, e i e j , dV s (e j -e i )
since Λ 0 = 1, using eq. (A.3), and because there is no jump at the sametime a.s, [Λ, F X] t = 0

= F 0 X 0 + P p=1 Kp k,l=1 k =l t 0 f p k , Y p s-Λ s-F s- Mp m=1 c p,m lk n∈I p m X s-, e n e n d N p kl,s -s - P p=1 Kp k,l=1 k =l t 0 Λ s-F s-X s-d N p kl,s -s + t 0 Λ s-X s-α s ds + t 0 Λ s-X s-β s , dV s + A t 0 Λ s-F s-X s ds + t 0 Λ s-F s-dV s + N i,j=1 i =j t 0 Λ s-(β j,s -β i,s )X s , e i a ji (e j -e i )ds + N i,j=1 i =j t 0 Λ s-(β j,s -β i,s
)X s-, e i e j , dV s (e j -e i ) by using eq. ( 12) and because X t X t , e n = X t , e n e n , ∀t ≥ 0

= F 0 X 0 - P p=1 Kp k,l=1 k =l t 0 Λ s-F s-X s-d N p kl,s -s + P p=1 Kp k,l=1 k =l t 0 f p k , Y p s-Λ s-F s- Mp m=1 c p,m lk N n=1 diag p m X s-, e n e n d N p kl,s -s + t 0 Λ s-X s-α s ds + t 0 Λ s-X s-β s , dV s + A t 0 Λ s-F s-X s ds + t 0 Λ s-F s-dV s + N i,j=1 i =j t 0 Λ s-(β j,s -β i,s )X s , e i a ji (e j -e i )ds + N i,j=1 i =j t 0 Λ s-(β j,s -β i,s )X s-, e i e j , dV s (e j -e i )
by using the definition of the diagonal matrix diag p see eq. ( 9) where we have used the eqs. ( 5) and ( 25), and since there is no jump at the same time a.s, L p,n kl , X t = 0.

= F 0 X 0 - P p=1 Kp k,l=1 k =l t 0 Λ s-F s-X s-d N p kl,s -s + t 0 Λ s-X s-β s , dV s + A t 0 Λ s-F s-X s ds + P p=1 Kp k,l=1 k =l Mp m=1 t 0 f p k , Y p s-c p,m lk diag p m Λ s-F s-X s-d N p kl,s -s + t 0 Λ s-X s-α s ds + t 0 Λ s-F s-dV s + N i,j=1 i =j t 0 Λ s-(β j,s -β i,
Taking now Ito's product rule on Λ t F t X t , and considering the eq. ( 16) and the eq. (A.4), we have: Λ s-L p,n kl,s-X s-λ q uv,s-d N q uv,s -s - f q u , Y q s-c q,m vu diag q m Λ s-L p,n kl,s-X s-d N q uv,s -s (i.e., before the next jump-time t w+1 ), and the initial condition Z tw described above. First, we multiply by the matrix Id N 0 N ×N on both sides of the eq. (A.6), where Id N is the identity matrix in R N ×N . We have:

Λ t L p,
-
G t - w+1 = Id N 0 N ×N exp Υ tw Γ tw 0 N ×N Ξ tw ∆t w+1 G tw K tw .
Second, we multiply by the matrix 0 N ×N Id N on both sides of the eq. (A.6), and we obtain: where we have used the fact that the matrix within the exponential is block triangular, see for instance [START_REF] Van Loan | Computing integrals involving the matrix exponential[END_REF].

K t -
The jump-adapted almost exact solution scheme for K t and G t can be easily obtained from eq. ( 27) and eq. ( 28) by adding the effect of jumps, resp.

Proof 5 [Corollary 4.3]. First, note that σ t (F t X t ) in eq. ( 24) can be rewritten by splitting it as: Regrouping the terms and using the definition of Φ t and Ψ t in eq. ( 29) and eq. ( 30) resp., the following SDE holds: Now, let us focus on the second term in the no-jump part of the SDE of eq. (A.7). First, because by definition α t = α(X t ) and β t = β(X t ), with α : X → R and β : X → R N known function with finite range, we can use the eq.( 6) under this functions. Thus, we have for α t :

d σ t (F t X t ) = A σ t (
d σ t (
X t α t = N n=1
α(e n ) X t , e n e n = N n=1 diag(α(e 1 ), ..., α(e N )) X t , e n e n = diag(α)X t ,

where α := (α(e 1 ), ..., α(e N )) ∈ R N , diag(α) is a diagonal matrix with α in the diagonal, and because we have used the fact that X t X t , e n = X t , e n e n for each n = 1, ..., N , and that X t = N n=1 X t , e n e n . Thus, applying σ t ( • ) and by linearity, we have σ t (X t α t ) = diag(α) σ t (X t ) . (A.8)

In the same way, using eq.( 6) under β t , and taking σ t ( • ), we have: β(e n ), e j -e i σ t (X t ), e n e n , but, because in the second term in the no-jump part of the SDE of eq. (A.7), we make inner product with e i ∈ X, i = 1, .., N , it holds σ t (β j,t -β i,t )X t , e i = β(e i ), e j -e i σ t (X t ), e i . Computing all the sum involved in eq. (A.7), we have that:

σ t (
N i,j=1 i =j σ t (β j,t -β i,t )X t , e i a ji (e j -e i ) = N i,j=1 i =j

β(e i ), e j -e i σ t (X t ), e i a ji (e j -e i ) = N i,j=1 i =j

β(e i ), e j -e i a ji (e j -e i ) e i σ t (X t ) , (A.9) since σ t (X t ), e i = e i σ t (X t ) ∈ R for any i = 1, ..., N .

Using then eq. (A.8) and eq. (A.9) into the SDE in eq. (A.7), the following holds:

  we define the transition rate matrix C p m = c p,m kl of Y p (whose transpose belongs to the Q-matrix class), by:

Figure 1 :

 1 Figure 1: Temporal evolution of breakers' states and its resp. component in the reference configuration of the grid.

Figure 2 :

 2 Figure 2: Temporal evolution of the hidden process H and its estimator H representing the grid over the reference configurations.

Theorem 3 . 2 .

 32 The estimation A = ( a ij ) of A = (a ij ), and C p m = c p,m kl of C p m = c p,m kl , for each p = 1, ..., P and m = 1, ..., M p ; are given for i = j and k = l, by:

Figure 3 :

 3 Figure 3: Representation of the jump-adapted scheme given by the jump times of all observed processes.

4. 4 .

 4 An Estimation Procedure for the Initial Parameters when M p = K p

  )) and C p(0) m = c p,m kl (0) , for each observed process p = 1, ..., P and index m = 1, ..., M p .

  0) , and of C p m = c p,m kl , i.e, the matrix C p(0) m = c p,m kl (0) , for each p = 1, ..., P and m = 1, ..., M p , are obtained as follows:

p t } t≥0 , p = 1 ,

 1 ..., 6, are obtained by the classical simulation procedure of jump chains and holding times with exponential distribution, see, e.g., [36, Section 2.6]. The Markov processes are performed under the fixed matrices of eq. (B.1) and Table B.1 in Appendix B.1.1. The sample path for X is shown on the left side in Figure 5. This simulation was stopped at 50 jump-events, giving a horizon time of two years from 2018 to 2020. The total number of jump-events of the observed states is 1964. Each simulated observed process is shown in Figure 4. Concerning the parameters of our model, we compute the initial estimation of the matrices A (0) and C p(0) m , p = 1, ..., 6, m = 1, 2, by the empirical estimation of eq. (37). These values are shown in eq. (B.2) and Table B.2 in Appendix B.1.2, resp. The initial filter estimate of X is obtained from eq. (

Figure 4 :

 4 Figure 4: Simulated temporal evolution of the observable processes.

Figure 5 :

 5 Figure 5: On the left side, a sample path of the hidden Markov state process X (blue line) representing the hidden state over time, and the estimation by clustering method (green points). At several times, the arg min set for clustering could not be a singleton. On the right side, the temporal evolution of the active clusters (blue lines) representing when the arg min set has more than one point (red points).

20

 20

Figure 6 : 2 klFigure 7 :

 627 Figure 6: Filter estimate of the hidden process X. The "real" state over time is shown by the blue dashed line. The filter is obtained by the parameter estimation at the first iteration (on the left side) and the final one (on the right side).

Figure 8 :

 8 Figure 8: Temporal evolution's uncertainty signal on the states of H at the first iteration (on the left side) and on the last iteration 14 (on the right side).

Figure 11 :

 11 Figure 11: Final filter estimate X for the temporal evolution of the grid.

Figure 12 :

 12 Figure 12: Temporal evolution's uncertainty signal on the reference configurations of the grid at the first iteration (on the left side) and on the last iteration 62 (on the right side).

0 ΛA t 0 Λ 0 Λ 0 Λ 0 Λ

 00000 s-L p,n kl,s-X s-d N q uv,s -s + s-L p,n kl,s X s ds + t s-L p,n kl,s dV s + t s-e n , X s-X s-λ p kl,s-dN p kl,s since there are common jumps, Λ, L p,n kl X t = t s-λ p kl,s--1 e n , X s-X s-

w+1= 0 N

 0 ×N Id N exp Υ tw Γ tw 0 N ×N Ξ tw ∆t w+1 G tw K tw = 0 N ×N Id N   exp Υ tw ∆t w+1 ∆tw+1 0 exp Υ tw (∆t w+1 -s) Γ tw exp Ξ tw s ds 0 N ×N exp Ξ tw ∆t w+1   G tw K tw = exp Ξ tw ∆t w+1 K tw ,

  Second, it is assumed that the transition probability rate of Y p from state f p l ∈ Y p to the state f p k ∈ Y p also depends on the local p-th position of the value that H takes, i.e., the p-th position in the vector h n ∈ H, see ) as a function of the p-th position of h n ∈ H. Indeed, let e n ∈ X. For each p = 1, ..., P and k, l = 1, ..., K p , k = l, the transition probability rate of Y p from state f p l ∈ Y p to the state

	eq. (1) and eq.(2).
	Mathematically, let proj p : H → M p the p-th projection function. We denote by H p t := proj p (H t ) the p-th projection in H t at time t ≥ 0, i.e., the value of H p t is m p m ∈ M p for some m = 1, ..., M p . In line with [15], we can relate each observed state process Y p with the underlying hidden process X by its transition rate
	matrix C p (X t ) ∈ R Kp×Kp , where the transpose of C p (X t ) belongs to the Q-matrix class. This matrix can
	be expressed as the sum of the eq. (6). In our case, we can obtain for each n = 1, ..., N an expression of
	the matrix C p (e n ) = c p kl (e n f p k ∈ Y p can be expressed by:

  ), N p kl,t are independent Poisson processes, and that they have fixed intensity one, see, e.g., [16, Lemma 1], [24, Ch. II, Theorem T6] and [28, Lemma 4.7.1] resp., mutatis mutandi.

3.2. Filter Estimate for the Grid States

The idea is to use Λ t to compute the estimator σ t (X t ) := E[X t | Y t ] by means of a version of Bayes' rule, see, e.g., [24, Ch. VI, Lemma L5]. More precisely, for any G t -adapted and integrable process F t , the filtered estimate of F t can be computed via:

  Second, O i,t is the time that the grid spends on the reference configuration h i ∈ H up

		t				t
	n kl,t :=	0	e n , X s-dN p kl,s	,	S p,n k,t :=	0	f p k , Y p s	e n , X s ds .	(21)
	Proof.						See Proof 1 in Appendix A.
	to time t ≥ 0. Third, L p,n 12,t (resp. L p,n 21,t ) increases only when the breaker p changes resp. from state "off " to
	"on" (resp. from state "on" to "off "), and simultaneously the grid is in the reference configuration h n ∈ H within the time interval [0, t]. Fourth, S p,n 1,t (resp. S p,n 2,t ) is the time that the grid spends on the reference
	configuration h

In Theorem 3.2, note that J ij,t represents a counting process of the number of jumps of X from state e i ∈ X to state e j ∈ X within [0, t], i = j, O i,t stands for the occupation time by X on the state e i ∈ X within [0, t], L p,n kl,t represents the process that increases only when Y p jumps from state f p k ∈ Y p to state f p l ∈ Y p and X is in state e n ∈ X, k = l; and S p,n k,t stands for the total time up to t ≥ 0 for which X is in state e n ∈ X and simultaneously Y p is in state

f p k ∈ Y p .

Example (Electrical grid). For instance, consider R = 1, M p = K p = 2 and M p = K p . First, J ij,t counts the number of jumps that the grid does from the reference configuration h i ∈ H to the reference configuration

h j ∈ H within [0, t]

. n ∈ H up to time t ≥ 0 and simultaneously the breaker p is in the state "on" (resp. the state "off ").

  Y p to state f p l ∈ Y p at time t = t w+1 , and σ tw+1 X t -

	where ∆t w+1 = t w+1 -t w , ∆N p kl,tw+1 = N p kl,tw+1 -N p kl,t -w+1 state f p k ∈	is defined by ∆N p kl,tw+1 = 1 if Y p t jumps from
			P	Kp
	w+1	+	p=1	k,l=1	Ψ p kl,t -w+1	σ tw+1 X t -w+1	∆N p kl,tw+1 ,
				k =l

w+1

:= lim s↑tw+1 σ tw+1 (X s ) is the values "before" the jump at time t w+1 .

  Y p to state f p l ∈ Y p at time t = t w+1 , and σ tw+1 L q,n

	where ∆t w+1 = t w+1 -t w , ∆N p kl,tw+1 = N p kl,tw+1 -N p kl,t -w+1 state f p k ∈ uv,t -is defined by ∆N p kl,tw+1 = 1 if Y p t jumps from w+1 X t -w+1 := lim s↑tw+1 σ tw+1 L q,n uv,s X s and
	w+1 σ tw+1 X t -	:= lim
		tw+1 ,

s↑tw+1

σ tw+1 (X s ) are the values "before" the jump at time t w+1 .

  expand that, since X t takes values in the space X of unit vectors of R N ,

	ln Λ t = -	P p=1	Kp k,l=1 k =l	Mp m=1	ln	c p,m lk c p,m lk	n∈I p m	L p,n kl,t -	i =j N i,j=1	ln	a ji a ji	J ij,t
					P			Kp		Mp						N
				+							c p,m lk	-c p,m lk		S p,n k,t +	a ji -a ji O i,t .
					p=1	k,l=1 k =l	m=1					n∈I p m	i,j=1 i =j
	Taking now E[ • | Y • ] above and using notation of eq. (A.1),
	P p=1 Kp k,l=1 k =l Kp + p=1 P P P p=1 p=1 k,l=1 k =l a derivatives: Kp k,l=1 k =l t 0 ln λ p kl,s-λ p kl,s-dN p kl,s -E ln d P dP Y t = -Mp m=1 ln c p,m lk c p,m lk n∈I p E L p,n kl,t | Y t -N ln i,j=1 m i =j N i,j=1 i =j t 0 ln a ji a ji Mp N Kp k,l=1 k =l t 0 λ p kl,s -λ p kl,s ds + N i,j=1 i =j t 0 X (A.1) a ji E J ij,t | Y t a ji dJ ij,s + m=1 c p,m lk -c p,m lk n∈I p m E S p,n k,t | Y t + i,j=1 i =j
	where λ p kl,t = f p k , Y p t-∂ ∂ a ji E ln m=1 Mp d P c p,m lk dP	n∈I p Y t = 0			and		∂ ∂ c p,m lk	E ln	d P dP	Y t = 0 .
	Thus, we have a ji =	P p=1 E[Oi,t|Yt] and c p,m Kp k,l=1 E[Jij,t|Yt] lk	t = 0	ln n∈I n∈I p m m p	λ p kl,s-λ p kl,s-E L p,n kl,t |Yt dN p kl,s -E S p,n k,t |Yt	N i,j=1	ln	a ji a ji	J ij,t
								P		k =l Kp	t					i =j N	(A.2)
					+	p=1	k,l=1	0	λ p kl,s -λ p kl,s ds +	i,j=1	a ji -a ji O i,t ,
										k =l						i =j
	where we have used the definition of the process O i,t in eq. (20). Now, noting that
														Mp
				λ p kl,t -λ p kl,t = f p k , Y p t		c p,m lk	-c p,m lk	X t , e n ,
													m=1		n∈I p m
	ln	λ p kl,t λ p kl,t	= ln	     	Mp m=1 Mp m=1	c p,m lk c p,m lk	n∈I p m m n∈I p	X t , e n X t , e n	     	=	Mp m=1	ln	c p,m lk lk c p,m	m n∈I p	X t , e n .
	Thus, by rearranging terms in eq. (A.2) and using the definition of L p,n kl,t and S p,n k,t see eq. (21), we have:

s , e i a ji -a ji ds , m X t , e n . The log-likelihood is, therefore:

ln Λ t = -ji -a ji E O i,t | Y t .

The unique maximum above over a ji and c p,m lk occurs at such values obtained by equaling to zero the partial .

Proof 2 [Theorem 3.3]. We proceed as in

[START_REF] Aggoun | Finite-dimensional models for hidden Markov chains[END_REF] Theorem 3.2]

but in our context. To derive an equation for σ t (F t X t ), we use first Ito's product rule on F t X t . Let ∆X s := X s -X s -. The following holds:

  s )X s , e i a ji (e j -e i )ds + (β j,s -β i,s )X s-, e i e j , dV s (e j -e i ) , where the last equality holds because X t = N n=1 X t , e n e n , ∀ t ≥ 0. Taking now E[ • | Y • ] above and interchanging expectation and integration[START_REF] Wong | Stochastic processes in engineering systems[END_REF], we get the result.Proof 3 [Theorem 3.6]. We proceed as in Theorem 3.3. To derive an equation for σ t L p,n kl,t X t , we use first Ito's product rule on L p,n kl,t X t . The following holds:L p,n kl,t X t = L p,n kl,0 X 0 +

								N	t
				i,j=1 i =j Λ s-t 0 t
			0	X s-dL p,n kl,s +	0	L p,n kl,s-dX s + L p,n kl , X t	(A.4)
		t			t					t	t
	=	0	X s e n , X s ds +	0	X s-e n , X s-d N p kl,s -s + A	0	L p,n kl,s X s ds +	0	L p,n kl,s dV s ,

  , Λ s-X s-e n dN p X t and X t e n , X t = e n , X t e n , ∀ t ≥ 0 .Taking E[ • | Y • ] above and interchanging expectation and integration[START_REF] Wong | Stochastic processes in engineering systems[END_REF], we get the result.Proof 4 [Theorem 4.1]. First, let's focus on the no-jump part of the SDE in eq. (27) and eq. (28), i.e, for t not be jump time, let say t ∈ (t w , t w+1 ), w = 0, ..., W , we considerd G t = Υ t G t + Γ t K t dt d K t = Ξ t K t dt ,which is an ordinary linear differential equation system. Recall Υ t , Γ t , Ξ t ∈ R N ×N are constant matrices matrices between jumps. Now, let2N . Thus we can redefine the ordinary differential equation system by: ×N is the zero matrix in R N ×N and Z tw is the initial condition ∀t ∈ (t w , t w+1 ). The solution of the latter differential equation with initial condition at t w is:Z t = exp Υ tw Γ tw 0 N ×N Ξ tw ∆t w+1 Z tw , (A.6)where ∆t = t -t w is the length of [t w , t]. Consider now the definition of Z t in eq. (A.5) but for t = t -

	P q=1 Mp m=1 u,v=1 Kq u =v t 0 f p 0 k , Y p t Λ s-L p,n kl,s-X s-d N q uv,s -s + A s-c p,m lk diag p m e n kl,s t 0 Λ s-L p,n kl,s X s ds + because X t λ p + kl,t = Mp m=1 f p k , Y p t c p,m lk diag p G t K t the column vector in R dZ t = Υ t Γ t 0 N ×N Ξ t Z t dt , Z tw = G tw K tw , m Z t := where 0 N w+1 t 0 Λ s-L p,n kl,s dV s (A.5)

  F t X t )dt + (F t-X t-)dN p kl,t + σ t (X t α t )dt (β j,t -β i,t )X t , e i a ji (e j -e i )dt -(F t-X t-)dN p kl,t ,wherein the integral w.r.t. Lebesgue measure dt, we are allowed to evaluate σ t (F t-X t-) as σ t (F t X t ).

			P	Kp	P	Kp
	p=1 σ s + k,l=1 k =l σ t (F t X t )dt -p=1 k,l=1 k =l N Kp σ t P	Mp	f p k , Y p t c p,m lk diag p m σ t (F t X t )dt
	i,j=1			p=1	k,l=1	m=1
	i =j				k =l
	P	Kp	Mp
	+		f p k , Y p t-c p,m lk diag p m σ t
	p=1	k,l=1	m=1
		k =l

  F t X t ) = Φ t σ t (F t X t )dt + σ t (X t α t ) + (β j,t -β i,t )X t , e i a ji (e j -e i ) dt

		N
	i,j=1 i =j kl,t-σ t (F t-X t-)dN p Ψ p kl,t , σ t + P Kp	(A.7)
	p=1	k,l=1
		k =l

  β j,t -β i,t )X t =

N n=1 (β j,t (e n ) -β i,t (e n )) σ t (X t ), e n e n = N n=1

That is, for each j = 1, ..., N , N i=1 a ij = 0 and a ij ≥ 0, ∀ i = j.

See, e.g.,[24, Ch. VI, Theorem T2] or [15, eq. (14)] for a general case.

Because, e.g., it could take several iterations to have the equality in all the parameters.

See, e.g.,[24, Ch. VI, Eq. (2.17)] or[START_REF] Aggoun | Measure theory and filtering: Introduction and applications[END_REF] Lemma 4.7.3].

the RTE-CentraleSupélec Chair.

In this section, our model is confronted with real data provided by France's transmission system operator (RTE). The available data consists of historical records of the evolution of the breakers' states over a period during which the grid operated normally.

Initial Parameters

We take from the data an electrical transmission grid with P = 6 breakers, with states off /on for each one. The observations' horizon time takes place on three months from 2016 -01 to 2019 -01, giving a total number of jump-events of 1264. The breakers' states over time are shown in Figure 9. Using the clustering method presented in Section 4.5, we choose N = 3 reference states defining the space H = {h 1 , h 2 , h 3 } by: H = (1, 1, 0, 1, 0, 1), (0, 0, 1, 0, 0, 0), (1, 1, 1, 1, 1, 1) .

(41)

Thus, the space of X = {e 1 , e 2 , e 3 } stands for the space of canonical vectors of R 3 , where e 1 = (1, 0, 0, 0), e 2 = (0, 1, 0, 0), e 3 = (0, 0, 1, 0). The initial filter estimate of the grid (represented equivalently by X) is obtained from eq. [START_REF] Pohle | Selecting the number of states in hidden markov models: pragmatic solutions illustrated using animal movement[END_REF]. Under the real data, we obtain σ 0 ( X 0 ) = (0.111, 0.785, 0.104). The initial state for X is therefore chosen to be X 0 = e 2 by eq. ( 32). Thus, by eq. ( 4), the initial state of the grid is h 2 , see eq. ( 41). The initial estimation of the matrices A (0) and C p(0) m , p = 1, ..., 6, m = 1, 2, is computed by the empirical estimation of eq. [START_REF] Choi | A survey of binary similarity and distance measures[END_REF]. Finally, for the stopping criteria of eq. ( 31), we fix ε = 10 -3 .

Grid State Estimation by Clustering Method

In Figure 10, we observe on the left side the temporal evolution of the empirical estimation of the grid by clustering method, see eq. [START_REF] Norris | Markov chains[END_REF]. At any time, all breakers states' joint information is classified with the minimum distance to the cluster centers, see eq. [START_REF] Wong | Stochastic processes in engineering systems[END_REF]. Here, the cluster centers are the reference configurations of eq. (41). However, the argmin set might not be a single point because of the nature of the dataset. In such cases, the clustering method is not exact because we cannot know at which configuration the grid is. The green dots show all the points of the arg min set. Because of the several choices, we represent the first choice of the arg min set by the blue line in the same picture. Also, we show on the right side in Figure 10 when a cluster is active, i.e., when the reference configurations can be chosen by clustering. The red dots represent when there is more than one choice. The metric distance used is the euclidean distance.

This is a crucial aspect in our estimation since the first estimation of the matrices A (0) and C p(0) m , p = 1, ..., 6, m = 1, 2, may be affected by this issue. So, by taking the first point of the arg min set at each time, we have estimated using eq. ( 37) these matrices. Thus, this estimation might not be correct and could affect the estimation by HMM.

Grid State Estimation by HMM

We show the filter estimate of the grid's temporal evolution over the reference configurations in Figure 11.

Recall that the estimation is represented (equivalently) by the process X over the canonical vectors of R 3 .

On the left side in Figure 11, the first filter estimation is done using the (empirical) initial parameter estimation by clustering method that could not be correct. So, we iterate the estimation of the parameters using Theorem 3.2. In the simulated case shown in Figure 6, it shows that the filter estimate fits better as the number of iterations increases. However, we can not compute the MSE between the real values of X and the filter estimate X, because the real grid's temporal evolution is not known. The filter estimate for the last iteration is shown in Figure 11 on the right side.

Uncertainty and Final Hidden State of the Grid

Now, we focus on the uncertainty in the grid by using t , see eq. ( 33). We show in Figure 12 the values of the uncertainty signal t that represents a temporal evolution's uncertainty signal on the reference configurations. We notice that there are fewer peaks in t at the last iteration compared to the first iteration. This is because finding the optimal parameters to fit our model, the uncertainty decreases as the number of iterations increases. Recall that when t is over δ t , the state at which the grid "jumps" is known with some uncertainty. However, such uncertainty is almost instantaneous because it does not remain in time, as reflected in Figure 12.

β(e i ), e j -e i a ji (e j -e i ) e i σ t (X t ) dt

The results follows by considering in Theorem 4.1,

β(e i ), e j -e i a ji (e j -e i ) e i .

Appendix B. Parameter Values

In this section, we present the matrices used in the Section 5 of our model's numerical results. , m = 0, 1, for each breaker p = 1, ..., 6.


Jonas H Osório

Joao B Rosolem

Fabio R Bassan

Foued Amrani

Frédéric Gérôme

Cristiano Fetah Benabid

M B Cordeiro

Cristiano M B Cordeiro

Hollow-core photonic crystal fibers as platforms for Power-over-Fiber applications

come

Introduction

The current framework in photonics identifies hollowcore photonic crystal fibers (HCPCFs) as building blocks for the next developments in the field. The recent research endeavors have allowed the achievement of technological breakthroughs (with remarkable reduction of the loss figures [START_REF] Jasion | Hollow core NANF with 0.28 dB/km attenuation in the C and L bands[END_REF][START_REF] Amrani | Low-loss single-mode hybrid-lattice hollow-core photonic-crystal fibre[END_REF][START_REF] Osório | Hollow-core fibers with ultralow loss in the ultraviolet range and sub-thermodynamic equilibrium surface-roughness[END_REF]) as well as the demonstration of a myriad of applications ranging from the development of optical sources [START_REF] Chafer | Near-and middle-ultraviolet reconfigurable Raman source using a record-low UV/visible transmission loss inhibited-coupling hollow-core fiber[END_REF] to the study of sensing opportunities [START_REF] Cordeiro | Azimuthally asymmetric tubular lattice hollow-core optical fiber[END_REF][START_REF] Guimarães | Angle-resolved hollow-core fiberbased sensing approach[END_REF].

In HCPCFs, the optical mode is guided through an empty core, and with minimal interaction with the microstructure that demarks the fiber core contour. Indeed, this property is even more salient for HCPCF guiding via the inhibited coupling (IC) mechanism, for which the power fraction of the core mode overlapping the silica surrounding can have values as low as 10 -6 [START_REF] Debord | Hollow-core fiber technology: the rising of 'gas photonics[END_REF]. The tiny overlap between the core mode and the silica microstructure unlocks the possibility of using HCPCFs to transport high-power beams at levels at which solidcore fibers cannot withstand. This capability of HCPCFs has been proven when dealing with both pulsed [START_REF] Hädrich | Scalability of components for kW-level average power few-cycle lasers[END_REF] and continuous-wave (CW) laser beams [START_REF] Zhu | Delivery of CW laser power up to 300 watts at 1080 nm by an uncooled low-loss anti-resonant hollowcore fiber[END_REF]. Remarkably, and very recently, a HCPCF has been demonstrated to be capable to transmit, for 1 km, a CW laser beam with 1kW power level at 1075 nm [START_REF] Mulvad | Kilometre-scale, kilowatt average power, single-mode laser delivery through hollow core fibre: overcoming nonlinear limits of glass fibre[END_REF].

In a different context, and via the use of solid-core fibers, power-over-fiber (PoF) systems have been widely studied in the past years [START_REF] Rosolem | Power-over-fiber applications for telecommunications and for electric utilities[END_REF][START_REF] Matsuura | Recent advancement in power-over-fiber technologies[END_REF]. PoF systems can convert the optical power emitted from a light source (and transmitted through an optical fiber) to electrical power thanks to the utilization of a photovoltaic (PV) converter. As typical PoF systems employ silica-core fibers, the utilizable optical power levels are inherently limited by silica properties. Indeed, approaches using large-core and double-clad fibers [START_REF] Li | A large-fore microstructure optical fiber for cotransmission of signal and power[END_REF][START_REF] Matsuura | 150-W power-over-fiber using doubleclad fiber[END_REF] have been reported as alternatives for increasing the power levels in silica-based PoF systems. Still, such systems remain fundamentally limited by the power levels that silica can withstand.

Thus, considering the potential of HCPCFs to deliver high-power laser beams, we here propose the utilization of this class of optical fibers for PoF applications. Hence, we employ a single-ring tubular-lattice (SR-TL) HCPCF to deliver a 976 nm laser beam onto a photovoltaic converter and activate a representative electrical circuit and display. To the best of our knowledge, this is the first application of a HCPCF in PoF. We consider that the demonstration reported herein identifies a promising path for the development of future PoF systems that can potentially lift the fundamental limitations which hinder the increase of the practicable power levels in PoF configurations based on silica-core fibers.

HCPCF structure and loss characterization

Fig. 1a shows the cross-section of the fiber we use herein. It stands for a SR-TL HCPCF whose 39 µm-diameter core (identified as a yellow circle in Fig. 1a) is demarked by eight 720 nm-thick untouching silica tubes with 20 µm diameter. The fiber transmission spectrum alternates highand low-attenuation regions, which are set by the cladding tubes' thickness, t, and the refractive index of the core and cladding tubes, n1 and n2, respectively. The high-loss regions occur at wavelengths λm = (2t/m) (n2 2 -n1 2 ) 1/2 , where m is an integer. Fig. 1b inset shows the transmission spectrum of a 3m-long section of the fiber used here and Fig. 1b main plot exhibits the fiber loss spectrum between 850 nm and 1350 nm. The fiber minimum loss figure is 1.0 dB/km at 1120 nm. At 976 nm (wavelength of the laser in PoF experiments), the fiber loss amounts to 4.0 dB/km.

Power-over-Fiber demonstration

To demonstrate the use of the above-mentioned HCPCF in a PoF configuration, light from a 976 nm laser (SkyEra TY976) has been launched into a 7 m-long HCPCF by butt-coupling its fibered output (62.5 µm core diameter) to the HCPCF -a schematic diagram of the experimental setup is presented in Fig. 2a and a picture of the buttcoupling site is shown in Fig. 2b. Optimization of the fibers' relative positions entailed a maximum transmitted power of 130 mW at the HCPCF output. Although the attained power level allows demonstrating a HCPCFbased PoF operation, there is considerable still room for improvement of the coupling efficiency between the laser output and the HCPCF, which is currently lower than 10%. Future work will, therefore, focus on improving such coupling efficiency, and, also, on using optical sources with higher power levels.

For the PoF demonstration, a silicon-based photovoltaic (PV) converter optimized for 915 to 980 nm sources (MH GoPower YCH-L200-06-FC) has been placed in front of the HCPCF end-face (Fig. 2c). Due to the small numerical aperture of the HCPCF (which, typically, has the order of 0.02 [START_REF] Delahaye | Double-clad hollow-core photonic crystal fiber for nonlinear optical imaging[END_REF]), the tip of the HCPCF has been placed 2.5 cm away from the window of the PV converter so to achieve adequate illumination on the latter. When shining a 130 mW-optical power level on the PV converter, a 1.39 mA current intensity have been measured between the terminals of the same. To illustrate the capability of using a HCPCF as a transmission channel in a PoF application, a representative electrical circuit (Texas Instruments MSP430) was connected to the PV converter. Fig. 2d shows that the light delivered by the HCPCF has been able to power up the tested display.

Conclusion

In this paper, we have proposed the employment of HCPCFs for PoF applications and described an experiment in which a 976 nm laser beam, after being transported by a 7 m-long HCPCF, has been used to activate a representative electrical circuit and display. To our knowledge, this is the first time that a HCPCF has been used to deliver light within a PoF context.

Fig. 1 .

 1 Fig. 1. HCPCF (a) cross-section and (b) loss spectrum (inset: 3m-long fiber transmission spectrum).

Fig. 2 .

 2 Fig. 2. (a) Setup diagram. Pictures of the (b) input coupling site, (c) the HCPCF output placement in front of the PV converter, and (d) the representative electrical circuit and display used in the PoF demonstration.

Acknowledgment J. H. Osório thanks Faepex/UNICAMP (grant 2499/21). J. B. Rosolem and C. M. B. Cordeiro thank CNPq (grants 303791/2021-7 and 309989/2021-3).