{'fraction_non_alphanumeric': 0.04108729239088451, 'fraction_numerical': 0.013735998455001931, 'mean_word_length': 4.388968388187849, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Recursive Gaussian Process over graphs for Integrating Multi-timescale Measurements in Low-Observable Distribution Systems

The transition to a smarter grid is empowered by enhanced sensor deployments and smart metering infrastructure in the distribution system. Measurements from these sensors and meters can be used for many applications, including distribution system state estimation (DSSE). However, these measurements are typically sampled at different rates and could be intermittent due to losses during the aggregation process. These multi time-scale measurements should be reconciled in real-time to perform accurate grid monitoring. This paper tackles this problem by formulating a recursive multi-task Gaussian process (RGP-G) approach that sequentially aggregates sensor measurements. Specifically, we formulate a recursive multi-task GP with and without network connectivity information to reconcile the multi time-scale measurements in distribution systems. The proposed framework is capable of aggregating the multi-time scale measurements batch-wise or in real-time. Following the aggregation of the multi time-scale measurements, the spatial states of the consistent time-series are estimated using matrix completion based DSSE approach. Simulation results on IEEE 37 and IEEE 123 bus test systems illustrate the efficiency of the proposed methods from the standpoint of both multi time-scale data aggregation and DSSE.


I. INTRODUCTION
Distribution system state estimation (DSSE) techniques infer the system states based on the network model and available measurements. The distribution system typically has limited number of measurement devices to monitor the medium and low-voltage feeders, rendering the system unobservable [1], [2]. The lack of measurement data hinders the development and use of DSSE. In recent years, the installation of different measurement sensors has increased significantly. For example, smart meters are being deployed in large numbers at the secondary side of the distribution systems. They are typically sampled at 15-min intervals and used for consumer billing purposes. The load composition of a primary feeder can be calculated according to the energy consumption of all the customers served by the feeder [3]. The aggregated smart meter measurements at the primary side are critical measurements for increasing the data redundancy in the distribution system [11], [12] performs imputations using all the measurements at once. That is, the approach proposed in [11], [12] involves batch processing and cannot be used to perform imputations in real-time as measurements arrive. Furthermore, these methods do not exploit the graphical structure of the grid.
Recently, sparsity-aware DSSE approaches are proposed to address the issue of low observability at the grid edge [13]- [17]. The compressive sensing-based approach estimates the states that are sparse in a linear transformation basis [13]. Matrix completion based DSSE approach exploits the sparsity of spatial states by suitable low-rank approximation [18], [19]. Tensor completion fills the missing elements in a tensor by exploiting the spatio-temporal correlation of the measurements [20]. A comparative analysis of these sparsity-based approaches along with their robust formulations was proposed in [21]. Authors in [22] use PMU and SCADA measurements for DSSE. This approach performs DSSE by incorporating a subset of these measurements available at time t along with the predicted SCADA measurements obtained using the information from the previous state estimates. It suffers from large measurement redundancy requirements (around 1.7), which makes it impractical for low-observable distribution systems. Furthermore, [22] does not consider any missing measurements scenario that could occur while aggregating measurements over finite bandwidth communication networks. A load evolution model for the slow-rate measurements is proposed in [23] for performing day-ahead forecasting. This approach relies on recursive Kalman filter (KF) updates for dynamic DSSE. However, KF typically needs the Hessian inverse computations at every step, which can be computationally burdensome. Also, the approach is not demonstrated for unbalanced systems. A first-order prediction-correction approach using PMU and smart meter data is performed in [24]. The main limitation of the approaches in [23] and [24] approaches is that they assume smart meter measurements are available at all load bus (i.e., system is fully observable). This paper proposes a recursive multi-task Gaussian process approach that sequentially aggregates multi-time scale measurements depending on the network connectivity information. It addresses multiple limitations of the state-of-the-art approaches. For example, the proposed approach is: -(1) Flexible to incorporate various heterogeneous measurements for unbalanced systems; (2) Does not require power measurements at all the load buses. The considered measurement dataset represents an unobservable condition that is used for both imputation and DSSE; (3) Effective even with missing data within each measurement time-series, and (4) Since inversion computations are only required at the initial time step, the proposed approach is computationally efficient.

C. Contributions
The main contributions of this paper are summarized below: • We propose a novel approach that imputes the heterogeneous measurements sequentially at any desired time resolution using recursive multi-task GP with or without topology information. • The proposed approach involves sequential measurement processing and can work with intermittent measurements.
Unlike the approach in [5], the proposed method is computationally efficient and is flexible to allow for both batch-wise and real time processing of measurements. • Finally, we leverage the graphical structure of the network in the recursive multi-task Gaussian process approach. We prove that exploiting the graph structure of the distribution system leads to a decrease in the variance of the imputed measurements. • Simulation results are carried out for the IEEE 37 and IEEE 123 bus test systems to verify the efficacy of the proposed approach. Relative to the linear interpolation approach [5], the RGP-G approach offers nearly 80% improvement in error performance while reconciling the multi time-scale measurements. We further estimate the spatial states of the consistent time-series measurements using the matrix completion based DSSE proposed in [18]. It can be inferred that accurate reconstruction of states is achieved even at 50% FAD (fraction of available data, which reflects the number of available measurements in the system).
II. BACKGROUND Consider a distribution system which can be perceived as a graph G = (V, E) where V ∈ R M are the nodes and E denotes the edges. The adjacency matrix A is defined as, In the distribution system, sensors are placed at a subset of M nodes. These sensors measure power injections or voltages at different locations in the network at different sampling rates. Hence, the main goal is to first reconcile these multi timeseries measurements at the desired time scale and estimate the states. To do so, we sequentially process the measurements using a recursive multi-task Gaussian process-aided state estimation approach. This approach exploits the spatial and temporal correlations of the sequential measurements located on the graph G. We propose to leverage the graphical structure of the grid for better imputation of the unevenly sampled measurements. Here, we consider the topology of the primary distribution system. Also, the imputation of the multi time-scale measurements is performed at the primary feeder. Therefore, the proposed approach does not rely on the model information of secondary side of the network. Based on the available network connectivity information, we propose two approaches, as shown in Fig.1. The inputs in all these approaches are the unevenly sampled time-series measurements, and output is the coherent set of measurements along with their variance. Conventional full GP and multitask Recursive GP (RGP) process the unevenly sampled measurements without utilizing any graph structure. Recursive processing of heterogeneous measurements using graph information in real-time and batch mode is performed using RGP-G interpolation and RGP-G prediction, respectively. We will review the relevant concepts of graph signal processing before introducing the proposed approaches. The graph Laplacian matrix L for graph G is defined as L = D − A where D is the diagonal degree matrix whose ith diagonal element is given by the sum of the elements in the ith row of A. The observations y = [y(1), ..., y(M )] ∈ R M represents a signal on graph G. The signal variation of y over graph G is measured as, The Laplacian quadratic form y Ly denotes the smoothness of the y. Suppose we want to recover the smooth signal y d from noisy observation y = y d + w, over the graph G. In order to recover this signal, an optimization problem is formulated as, where α ≥ 0. The global solution is, Here, I M is the identity matrix. The optimal solution y * d can be seen as the graph filtering of y using the graph filter S = (I M + αL) −1 [25] [26]. This graph filter will be used in the proposed RGP-G approach for inducing the graph structure of the distribution system. In the next section, we will review the conventional full GP approach.

A. Full GP Approach
Consider a distribution system with M buses and d types of sensor tasks. Here, the sensor tasks refer to different sensor measurements available, e.g., aggregated active and reactive power injections or voltage magnitudes at the primary feeder. We consider the availability of measurements for T time instances. Let the measurements corresponding to time instant x t be y t ∈ R dM . The measurements y t are obtained by concatenating measurements from different sensor locations i.e., y t = [y 1 t , ..., y d t ] . Each y d t is drawn from a noisy process as, where, f d t ∈ R M and σ 2 is the noise variance. The entries in y d t are zero at the locations where there are no sensor measurements.
The GP prior function f d t associated with d th sensor task at time t has distribution given as, The where ⊗ denotes the Kronecker product, I c ∈ R d×d is an identity matrix between different sensor tasks in the distribution grid. For instance, a distribution grid may have active power (P), reactive power injections (Q), and voltage (V) measurements at the primary feeder. Thus, there are three sensors tasks, and I c has a size of 3 × 3. In some distribution systems, there are only P and Q measurements available. The voltage measurements are available only at the substation. Therefore, in this case, the size of I c is 2 × 2. The kernel matrix K represents the temporal covariance functions within this sensor task. There are different kernel choices, with one of the most popular being RBF (radial basis function) kernel [27] corresponding to: where hyperparameters l and σ 2 s are the length-scale and signal variance respectively. The lengthscale of the kernel function controls the smoothness of the GP function [27].
If all the measurements upto time T are represented asỹ = vec(y 1 , ..., y T ), the distribution ofỹ using (4) and (6) is given as,ỹ = N (0, ((I c ⊗ I) ⊗ K) + σ 2 I) where K is the kernel matrix defined for all time instances t = 1, ..., T . Here, the i th and j th entry of K is given as The main goal of the Gaussian process-based imputation process is to infer the unknown test values y * corresponding to the time x * given the measurementsỹ at time x and the modeled GP prior function f (·). The measurementsỹ and the test values y * are jointly Gaussian whose distribution is given as, where, the matrices A, D and F corresponds to, Here, The conditional distribution of the test values y * givenỹ is a Gaussian distribution [27] with mean and covariance, m * = D A −1ỹ (13) and Algorithm 1 summarizes the full-GP approach. (13) and (14) involves inverting the matrix A for all the x time instances which is computationally expensive. The full GP approach suffers from the following drawbacks: • The GP prior function and the corresponding measurementsỹ as defined in (8) is a simple multi-task Gaussian process with an independent kernel function among the different measurements obtained at M nodes. • This approach performs training using all the measurements in the batch, and thus the training is performed off-line in a batch mode.
is the size of x, d is the total number of sensor tasks, and M are the nodes. The inversion of the matrix A is the key contributor to this complexity. To overcome these challenges, we propose a recursive GP approach that sequentially processes the measurements corresponding to each x t by using the knowledge of graphical structure of the distribution grid. In the next section, we will formulate the RGP-G approach when the network connectivity information is known. Then, we will develop the RGP-G method when the graph information is unknown.

Algorithm 1 Full GP Approach
Input: Aggregated Active and Reactive power injection measurements at load bus y t corresponding to time x t ,ỹ = vec(y 1 , ..., y T ), Kernel choice and hyper-parameters associated to kernel function 1: Calculate the kernel matrix K that exploits the temporal correlation using any kernel function (e.g., RBF kernel). 2: Calculate the matrix I c ∈ R d×d . 3: Calculate the matrix A, D and F by means of (10), (11) and (12) respectively. 4: Perform imputation at time x * by means of the mean m * and covariance matrix C * using (13) and (14). 5: return m * , C * .

III. PROPOSED APPROACH
This section presents the formulation for recursively imputing the multi time-scale measurements with and without topology information.

A. RGP-G Approach
One of the challenges in the full GP approach is the need to to use the complete vector x ∈ R T . To overcome this challenge, we aim to use the basis vectors x = [x 1 , x 2 , ..., x n ] where n T . We perform all the calculations on the basis vectors x ∈ R n which are fixed in number and locations. The function f = f (x) is the GP function corresponding to the basis vectors x. When the network connectivity information is known, we can construct the graph filter matrix S as defined in (3). In order to exploit the spatial correlation and induce the graph information, the observations in the GP function (4) and (5) are modified as, where, f d t is defined in (5). Using (5) and (15), we obtain the distribution of y d t as, The distribution of y t obtained by concatenating y d t is given as, Here, the identity matrix I c given in (6) is replaced by the kernel matrix K c that represents correlation among different sensor tasks.
The main aim of this section is to recursively update the mean and covariance of the multi-task prior function f as the measurements arrive at time t = 1, ..., T by incorporating the graph structure of the grid. We assume that the hyperparameters of the kernel function are known apriori using historicalbased data. The proposed RGP-G approach has the flexibility of performing both interpolation and prediction described as, 1) RGP-G Interpolation-This approach operates over a set time frame (24 hours as an example). The GP function is updated recursively at those time instances where the measurements are obtained. Once the GP function in that batch is updated, we perform imputation at the finest time resolution. Here, the finest time resolution refers to the narrowest time resolution between the different measurement sources. 2) RGP-G Prediction-This approach reconciles the multitime scale measurements in real-time. The imputation at the desired time resolution is performed by predicting the GP function until the subsequent measurements is observed. Here, the prediction is performed at the finest time resolution. 1) RGP-G Interpolation: We assume that the GP prior function f at time t = 0 has an initial distribution, with mean µ f g,0 = 0 and covariance C f g,0 defined as, Here, K = K(x, x). The initial covariance of the GP prior function exploits the spatial correlation using graph structure and the temporal correlation between the time instances. The measurements y t ∈ R dM arrive sequentially at time t = 1, .., T . The goal is to calculate the posterior distribution at time t, where y 1:t = (y 1 , ..., y t ), by combining the new measurements y t with the distribution, The desired posterior distribution is expanded according to [28], where, f t is the GP function at time t and c t is the normalization constant.
Calculation of the posterior is performed in two steps: a) Inference: In this step, we infer the joint prior p(f , f t |y 1:t−1 ) using the measurements received upto time t − 1. Here, the matrices A g , D g , F g calculated at time x t be defined as, Here, the subscript g refers to notations related to recursive GP with graphs approach. The goal is to calculate the joint prior p(f , f t |y 1:t−1 ) using the information from the prior p(f |y 1:t−1 ). This can be achieved using the chain rule as, The first term p(f t |f ) follows from the assumption that f t is conditionally independent of the past measurements y 1:t−1 given f . As any finite representation of a GP is Gaussian, the joint prior is also Gaussian. Hence, the conditional distribution p(f t |f ) is Gaussian and calculated by Gaussian identities given as, Using Gaussian identities and Woodbury formula, the solution to (26) is a joint Gaussian p(f , f t |y 1:t−1 ) = N (q, Q) with mean and covariance defined as, and, b) Update: This step updates the joint prior f with new measurements y t arriving at time t. The function f t is updated by Kalman filter update step which yields, where, Here, µ p g,t and C p g,t are obtained from (27) and (32), respectively. Here, G t = C p g,t (C p g,t + σ 2 I) −1 is the Kalman gain. The update is performed at time t where measurements are available. For instance, the aggregated smart meter at the load buses provides measurements at intervals t = 1, 16, 31, ..., T . At t = 16, assume that few of the measurements y t are missing due to communication bottleneck. Therefore, at those locations, we do not update but use the predicted measurements µ p g,t .
The posterior function f has Gaussian distribution The function f is sequentially updated with the observations y t until t = T . Once updated, the imputation of y * at time x * is performed using the following steps, and where the matrices D g , A g and B g are evaluated for time x * . The complete RGP-G Interpolation approach is illustrated in Algorithm 2.

Algorithm 2 RGP-G Interpolation Approach
Input: Basis vector x, Distribution system graph laplacian L ∈ R M ×M , total time instants T , α, S, K, K c , y t 1: Initialization: µ f g,0 = 0 and C f g,0 = (K c ⊗ S 2 ) ⊗ K. Calculate the gain matrix J g,t according to (29). 4: Calculate mean µ p g,t using (27) and covariance matrix C p g,t using (32).

5:
Calculate the gain matrixG t according to (38). 6: Set y t = µ p g,t at locations where y t is missing. Calculate mean µ f g,t by means of (36) and covariance matrix C f g,t by means of (37). 7: end for 8: Perform imputation at time x * by means of mean m * g (39) and covariance matrix C * g (40). 9: return m * g , C * g .
2) RGP-G Prediction Approach: The RGP-G interpolation approach operates over a set time frame. However, it is critical to reconcile the measurements as and when they arrive. The RGP-G prediction approach achieves the reconciliation in real-time by performing the step ahead prediction of the GP function using the knowledge of the past measurements. These predictions are the imputed values at the narrowest time resolution.
The complete algorithm of the proposed sequential prediction over graphs is summarized in Algorithm 3. Fig.2 shows the proposed RGP-G prediction approach. We demonstrate this approach by illustrating an example. At time t = 0, we initialize the GP function as given in (18). At time t = 1, we receive measurements y 1 . The GP function is updated using these measurements y 1 by means of (36) and (37). The updation of the mean and covariance of the GP function f are denoted by µ f g,t and C f g,t . If any measurements at time t = 1 are missing, they are predicted using (27) and (32). From time t = 2 onwards, we perform step-ahead prediction of the GP function using (41) and (42) until the subsequent measurements are observed. We perform prediction for all x(t * ) satisfying x(t * ) > x(t) using the knowledge of function f updated at the previous time instant. The predicted mean m * g and their associated variances C * g are given as, and Algorithm 3 RGP-G prediction Approach Input: Basis vector x, Distribution graph L ∈ R M ×M , total time instants T , α, y t , S, K, K c 1: Set µ f g,0 = 0 and C f g,0 = (K c ⊗ S 2 ) ⊗ K. 2: Calculate the gain matrix J g,t according to (29). 3: Calculate mean µ p g,t using (27) and covariance matrix C p g,t using (32). 4: Calculate the gain matrixĜ according to (38). 5: Set y t = µ p g,t at locations where y t is missing. Calculate mean µ f g,t by means of (36) and covariance matrix C f t by means of (37). 6: Predict for all x * which satify x(t) < x(t * ) with mean m * g and variance C * g using (41) and (42) respectively. 7: return m * g , C * g .

B. Multi-task Recursive GP (RGP)
When the network topology is unknown, the graph structure information by means of the graph filter matrix is not possible. In this case, we aim to use the recursive multi-task GP without graphs approach (RGP). Here, we initialize the prior function f = f (X) at time t = 0 as, with mean µ f 0 = 0 and covariance C f 0 defined as, We define the matrices A, D and F as, Instead of the matrices (22), (23) and (24) defined for RGP-G approach, we use the matrices defined in (45), (46) and (47). The sequential inference and update step remains the same for both approaches. In this approach, we can perform both the interpolation and prediction similar to the RGP-G approach.
In case of RGP Interpolation approach, the predicted mean and covariance is denoted by m * and C * given as, The RGP approach used for interpolation is summarized in Algorithm 4. The computational complexity associated with Algorithm 4 Multi-task RGP Approach Input: Basis vector x, total time instants T , α, sequential measurements y t ∈ R dM , K, K c 1: Initialization: µ f 0 = 0 and C f 0 = (K c ⊗ I) ⊗ K. Calculate the gain matrix J t according to (29) using the matrix A and D from (45) and (46) respectively . 4: Calculate mean µ p t using (27) and covariance matrix C p t using (32).

5:
Calculate the gain matrixG t according to (38). 6: Set y t = µ p t at locations where y t is missing. Calculate mean µ f t by means of (36) and covariance matrix C f t by means of (37). 7: end for 8: Perform imputation at time x * by means of mean m * (48) and covariance matrix C * (49) . 9: return m * , C * . Algorithm 2, 3 and 4 for n basis vectors and dM number of observations at step t is O(dM n 2 ). This complexity is driven by the gain matrix calculation in (29). Use of recursive GP significantly reduces the computational complexity as compared to the full GP approach.
We next show the use of graph information in RGP-G approach reduces the uncertainty of the posterior distribution when the measurements are recursively processed at time t = 1, ..., T .
Theorem 1: The variance of the estimator of f using the RGP-G (Algorithm 2 and 3) of the distribution p(f ) is less than the variance of the estimator of f using RGP at time t = 0 i.e., tr where C f 0 is defined in (44) and C f g,0 is defined in (19). Proof. In order to prove this theorem, we need to show that the trace of ∆C 0 = C f 0 − C f g,0 is nonnegative. The initial covariance matrix C f g,0 at time t = 0 is defined with a graph filter S given in (3). The Laplacian matrix L used in this graph filter has an eigen-decomposition corresponding to, where, U G = diag(U (1), U (2), .., U (M )) and V denote the diagonal eigenvalue matrix and the associated eigenvectors respectively. Every eigenvalue of the laplacian matrix L is non-negative [29]. We need to prove, The kernel matrix K and K c is positive semidefinite by construction. Hence, we need to prove, Let U (s) denote the smallest non-zero eigenvalue of L, we then need to prove (1 + αU (s)) −2 ≤ 1, For α > 0, the smallest non-zero eigenvalue of L i.e., U (s) > 0. Hence, tr(C f 0 − C f g,0 ) > 0. Theorem 2: The posterior covariance matrix of the estimator of function f for RGP-G approach at time t i.e., C * g in (40) is smaller than the posterior covariance matrix for the RGP based estimator of f given as C * in (49) evaluated without graph information ∀t = 1, ..., T .
Proof: We need to prove, tr(C * g ) < tr(C * ) for all t = 1, ..., T , where T )J t * and C * g is defined in (40). From Theorem 1, we have proved that ∆C 0 = tr(C f 0 − C f g,0 ) > 0. The matrix B for RGP approach can be defined as, Similarly, we have B g for RGP-G approach as defined in (28). Thus, we have, As tr(C f 0 −C f g,0 ) > 0 and the Schur complement of a positive definite matrix is also positive definite. Therefore, we have Similarly, Similarly, tr(C f t − C f g,t ) > 0 ∀t = 1, ..., T . The posterior distribution of the imputations tr(C * − C * g ) > 0.
Next, we discuss the robustness of the graph filter to topology uncertainties.

C. Robustness to topology uncertainties
The distribution grid topology can be unreliable and incorrectly estimated [30], [31]. Thus, it would be desirable for our GP predictions to be robust against the uncertainties in topology. The following discussion presents the conditions under which the proposed graph filter will be stable to perturbations in the graph topology, as proved in [32]. The graph filter is said to be stable against a perturbation if the perturbation does not lead to large changes in the filter output. The graph Laplacian matrix L is also known as the graph shift operator (GSO) which can be decomposed as L = VU G V , where U G and V denote the eigenvalues and the eigenvectors of the Laplacian matrix, respectively. Let the Laplacian matrix of the perturbed graph G p be L p . The magnitude of the error matrix is defined as E 2 = L − L p 2 . In [32], a graph filter g is said to be linearly stable for any GSO if, for any GSO ∆ and ∆ p , we have the following conditions satisfied, i.e., where, C is a positive constant. The graph filter defined in (3) can also be written as, where, y, V l is the graph Fourier transform of y on the vertices of the graph. Equivalently, y d * =ĥ(L)y, wherê h(L) = 1 1+αU can be viewed as low-pass filter. The low-pass filter is said to be linearly stable for perturbations in the graph as it satisfies (54). This property was proven in [32] and given as, Complete analysis of the impact of topology uncertainties on graph filter stability and eventually the RGP-G approach performance will be pursued as part of our future work.

IV. MATRIX COMPLETION BASED DSSE
While not the primary focus or contribution of our work, we provide a brief summary of the matrix completion based DSSE proposed in [18] for the sake of completeness. Unlike [18], [21] that assumes time synchronized subset of measurements, here the GP based reconciled measurements are used within the Matrix completion (MC) based DSSE. MC based DSSE estimate the spatial states of the network (i.e., the voltage phasors and power injections of all the buses at a single instant of time) by exploiting the sparsity of raw measurements. Specifically, matrix completion aims to estimate the complete matrix X from an incomplete and noisy observation matrix by suitable low rank approximation.
The consistent multi time-scale measurements are limited to specific spatial locations in the network where measurements are aggregated. Assume that the measurements at the slack bus are known. Thus, we use the measurements at the nonslack buses to construct a data matrix. Let m denote the set of phases at all the non-slack buses. The noisy matrix Z is constructed such that each row represents a phase and each column represents the measurement associated with the phase of each bus. For each b ∈ m, each row of the matrix Z ∈ R m×n , n = 5 is structured as, where, P b and Q b represent the active power and reactive power injections at each phase of non-slack bus b respectively. The terms (v b ) and (v b ) represent the real and imaginary parts of voltage phasors at each phase of non-slack buses respectively. Let Ω ⊆ {1, ..., m} × {1, ..., n} describe the known entries in Z. The observation matrix P Ω (Z) is represented as, The matrix completion formulation (59) recovers the complete low-rank matrix, aŝ Here, the nuclear norm X * is the sum of the singular values of the matrix X. (60) and (61) captures the linearized power-flow constraint relating voltage phasors v and voltage magnitude |v| to the power measurements as given in [33]. More details about the matrix completion based DSSE can be found in [18], [21]. Additionally, the impact of uncertain topology on matrix completion based DSSE is considered in [8], [21]. Additionally, a more comprehensive integrated robustness analysis of graph filter and matrix completion based DSSE will be considered in our future work.

V. SIMULATION RESULTS
The efficacy of the proposed approach is verified on the three-phase unbalanced IEEE 37 bus [34], and IEEE 123 bus test system [35]. An aggregated 24-hr load profile at the primary nodes consists of a mixture of load profiles, i.e., industrial and commercial load profiles obtained from [36], and residential loads obtained from [37]. Reactive power profiles are obtained by assuming a power factor of 0.9 lagging. Other profiles at different nodes were obtained by adding a random noise term and a sinusoidal wave of random amplitude spanning the 24-hr period. By utilizing this data, the voltage profile at all nodes is obtained by running load flow. The aggregated smart meter data are averaged over 15-min intervals while the voltage magnitude measurements are sampled at a 1-min interval. Thus, we have considered two sensor types for the case study. We have assumed RBF kernel for all the GP-based approaches. The imputation is performed for the aggregated smart meter data at a 1-min interval. We compare the performance of the proposed RGP-G Interpolation against the linear interpolation approach [5], RGP (Algorithm 3), and full GP (Algorithm 1) approach. The RGP-G prediction approach is compared with [38]. Algorithms 2, 3, and 4 are initialized using their respective mean and co-variance function associated with the GP function at time t = 0. Here, the hyper-parameters associated with the GP function can be obtained by either training the proposed approaches using historical data or using cross-validation techniques. The hyper-parameters involved in the proposed approach are θ = [l, σ 2 s , σ 2 ], where l, σ 2 s , σ 2 , are defined in (7). We have used the grid search method guided by a five-fold cross-validation technique to obtain the hyper-parameters for our problem. In the cross-fold validation technique, one fold of the measurement set is retained as a validation set, and the other folds as a training set. Each time a different set is chosen as the validation set, and this procedure is repeated five times. We select a finite set of reasonable hyper-parameter values to perform a grid search. The performance of each combination is evaluated through cross-validation on the training set. This approach evaluates the MAPE for each possible combination of hyperparameter values and chooses the set that minimizes the error on the validation set. More details on the grid-searchbased cross-fold validation technique for Gaussian process hyperparameter tuning can be found in [27]. Another approach is to consider the historical data for hyperparameter tuning. The historical measurements of multi time-scale measurements can be used to obtain the hyper-parameters by maximizing the log marginal likelihood of the historical time-series data. The log-likelihood can be computed in closed form as given in [27]. It is important to note that the proposed approach does not require any extra training set for imputation. The parameter α for the RGP-G approach is set to 0.05. There are three cases by which we illustrate the performance of the multi-task RGP-G approach.
1) Case 1: In this case, we consider the measurement noise as mean 0 with standard deviation equal to 1% of the actual values. Fig. 3 shows the performance of the RGP-G interpolation approach at 0% missing measurement case for an IEEE 37 bus test system. Here, the time-series is the active power injection at node 11 of phase A. The RGP-G interpolation approach recursively updates the GP function in the 24-hr batch and later performs imputation at 1-min interval. The 95% confidence interval indicates the uncertainty bounds associated with the imputed measurements. The ideal case, i.e., 0% missing measurements dataset, has no missing measurements, but the dataset consists of a subset of the total measurements, representing an unobservable condition. For instance, let us assume that the aggregated AMI measurements in this dataset are available at 15-min intervals. Therefore, if we consider a 24-hr duration to perform imputations every minute, we have only 96 measurements per AMI sensor out of the total 1440-time instances. 2) Case 2: In this case, we perform the reconciliation and state estimation using the noisy time-series measurements corrupted by 0 mean and 10% standard deviation for a total of four hours duration. In this case, we fix the number of sensors and their locations corresponding to a particular FAD. The meters are placed randomly in the network. Fig.4 shows the meter deployment for IEEE 37 bus test system for 50% FAD. Thus, there are no timeseries measurements at the nodes where the sensors are absent. We then perform imputation using these incomplete measurements. Fig.5 shows the comparison of the RGP-G approach with linear interpolation. It can be seen that the former approach provides smoother imputation than the latter one. We compare the performance of all the five approaches using the mean absolute percentage error (MAPE) metric. Tables I and II tabulate their performances for IEEE 37 and IEEE 123 bus test systems, respectively. It can be inferred that the performance of the RGP-G interpolation is superior to the other four approaches. In the linear interpolation approach, each time-series data is imputed individually without exploiting any spatio-temporal property of the data. In contrast, the full GP approach (Algorithm 1) exploits temporal relationships for imputation. The GP function update is performed using all the measurements in batch mode at once, which is computationally expensive. The proposed approach can impute at any missing measurement level and the uncertainty as well as MAPE will increase with the increase in percentage of missing measurements. Uncertain imputed measurements affects the state estimation process. The knowledge of the uncertainty in imputations is used to guide the DSSE process using a Bayesian MC framework as proposed in [12]. Analytical bounds on the estimation error of matrix completion approach in the presence of missing and noisy measurements are derived in [39]. Since the primary goal of our work is to introduce the novel RGP-G and RGP algorithms as viable options to deal with multi time-scale measurements, derivations of error bounds based on [39] will be pursued as a part of the future work.   The proposed recursive GP approaches assumes that the measurement data is corrupted by Gaussian noise as seen from (15). Also, the optimization formulation in (2) assumes Gaussian noise. These are common assumptions used in many prior efforts [25]- [27]. The proposed approach in its current form can be applied to non-Gaussian noise but will not be optimal. Table III shows the performance of the RGP-G approach with Laplacian noise with 0 mean with standard deviation equal to 5%, 10% of the actual power values. As seen from Table III, it can be inferred that the performance of the recursive GP approach under non-Gaussian noise scenarios is not optimal. Alternately, under a non-Gaussian measurement noise scenario, a warped GP approach [40] can be used. In this approach, the observations are transformed into a latent space such that the transformed data has Gaussian noise and will be better modeled by the GP. Developing a multi-task recursive GP framework using the warped GP will be pursued as part of our future work. The consistent time-series measurements are further used to estimate the states using the matrix completion-based DSSE approach discussed in section IV. While performing the matrix completion for a particular FAD, the corresponding entries in the measurements matrix are zero, with no sensor measurements. Table IV shows the absolute errors and relative error reductions for RGP-G and linear interpolation methods. It can be deduced that RGP-G based technique significantly reduces error at all FADs. For example, the error in estimating reactive power using GP-based imputed time-series is reduced by 46% at 90% FAD compared to the linearly interpolated time-series. The reduction in voltage state estimation error is more modest mainly due to the robustness of matrix completion based DSSE that includes the topology information. 3) Case 3: In this case, we compare the RGP-G prediction aided matrix completion with [38]. The method in [38] uses a data collation method to reconcile heterogeneous measurements and a Kalman filter method to perform DSSE. The data collation consists of an exponential moving average method to extrapolate the slow-rate measurements. Table V shows the prediction errors for different percentages of missing temporal data. Our proposed approach has several advantages over [38]. Firstly, the DSSE in [38] requires the measurement redundancy (ratio of number of measurements to the total states) higher than two. Thus, unlike our proposed approach, the method in [38] fails in low-observable conditions. Secondly, it can be observed that the imputation error in [38] is higher than the proposed RGP-G prediction approach as seen from Table V.
We now discuss the scalability of the proposed approach for large test systems.

A. Scalability Analysis
As discussed in section III, the computational complexity associated with Algorithm 2, 3 and 4 for M node system with d different sensor data streams and n time instants will be of the order of O(dM n 2 ). Hence, while the approach can be used for M = 8500, the complexity grow as M increases. To address the scalability issue, a distributed implementation of the proposed multi-task recursive GP approach is possible. We perform the distributed implementation on the 11,000node feeder proposed in [41]. The 11,000-node test feeder is constructed by connecting an IEEE 8,500-node test feeder and an EPRI Ckt7 test feeder at the substation. To perform the distributed recursive GP approach, we partition the 11,000 node network into four areas. We assume that smart meter measurements are available at 30% of nodes in each area. The load profiles assigned to each node consist of industrial, residential, and commercial load profiles and and scaled according to their base loads provided in [41]. Reactive power profiles are obtained by assuming a power factor randomly varying between 0.9 and 0.95 lagging. The smart meter measurements are averaged over 15-minute intervals with measurement noise as mean 0 and standard deviation equal to 1% of the actual power values. We perform the multi-task RGP Interpolation approach in each area for 4 hours. The MAPE for each area is tabulated in Table VI.

VI. CONCLUSION AND FUTURE WORK
This paper proposes a recursive Gaussian process with graphs for effectively aggregating heterogeneous intermittent time-series data and using it to estimate the distribution system states in low observability conditions. The proposed approach leverages the graphical structure of the network for accurately imputing the multi time-scale measurements. It has the flexibility to perform imputations in batch mode or realtime mode. Superior imputation performance of the active and reactive power time-series measurements are obtained with the proposed approach. Further, state estimation in IEEE 37 and the IEEE 123 bus system reveals that the power and voltage states are recovered with high fidelity.
Our proposed novel approach has significant strengths and certain limitations which will be addressed as part of our future efforts. These limitations include: 1) The proposed approach is sensitive to outliers in the measurement data. Hence, we aim to develop a robust Gaussian process framework against outliers as a part of our future work. 2) The performance of the proposed approach depends on the hyper-parameter values of the GP function. Hence, future work will involve recursively learning the hyperparameters as the GP functions are updated.
3) The proposed GP function assumes that measurement data is distributed as multivariate Gaussian. Developing a multi-task recursive GP approach for non-Gaussian noise will be pursued as a part of our future work.