Distributed PV generation estimation using multi-rate and event-driven Kalman kriging filter

: The ever-growing penetration of cost-effective photovoltaic (PV) panels within the distribution grid requires a robust and efficient method for PV system monitoring. Especially, the geographical proximity of PV panels can play an important role in lowering the dimension of measurements required for full system observability. Furthermore, the direct impact of variable cloud formation and uncertain propagation necessitates the development and validation of a spatiotemporal model. Accordingly, this study presents the modelling and validation of the spatiotemporal variability of solar power indices at 1 minute resolution for the scale of a residential neighbourhood. The spatiotemporal model is then applied to a Multi-Rate and Event-DRIven Kalman Kriging (MREDRIKK) filter to dynamically estimate behind-the-meter PV generation. The Kriging step exploits spatial correlations to estimate PV power output at locations from where measurements are unobserved. The multi-rate feature of the MREDRIKK filter enables the sampling of measurements at a rate much lower than the temporal dynamics of the associated states. A comprehensive study is undertaken to investigate the effect of multi-rate and event-driven measurement updates on the performance of the MREDRIKK filter. In addition, the superior performance of MREDRIKK filter is represented as compared to the persistence method irrespective of the observation size.


Nomenclature
number of sensors S set of sensor locations s k longitude-latitude pair for the sensor k T sampling interval t discrete time index x t value of time series at t τ amount of discrete time lag a t, τ autoregressive coefficient at t for lag τ w t Gaussian random noise at t with zero mean and variance σ w, t 2 p t order of autoregressive model at t expectation operator F t state transition matrix at t 1 p t p t -dimensional vector of ones C w, t noise covariance matrix at t C x, t state covariance matrix at t O set of observed sensors with cardinality m U set of unobserved sensors H t m × p t -dimensional observation matrix at t v t measurement noise vector at t with zero mean and diagonal covariance matrix C v, t D decimation factor ℤ + set of non-zero positive integers δ t (s r ) event detection threshold for location s r at t θ event detection sensitivity factor C r, k spatial covariance between sensor locations s r and s k α r, k Kriging weight for unobserved location s k with respect to observed location s r γ(d) semivariogram at Haversine distance d SPI(s r , t) solar power index for location s r at t l frame frame length C cluster index

Introduction
The popularity of photovoltaic (PV) systems is increasing among residential electricity customers as installed costs continue to decrease [1]. However, PV generation can become highly variable during times of intermittent cloud activity. Enhanced information about uncertain generation is crucial to maintain a reliable and stable smart grid as it enables the utilities of enhanced and realtime control of distributed PV systems [2]. However, most of the distributed PV generation monitoring is administered by the respective PV developers, whose profit-seeking data service is cost prohibitive. As a consequence, utilities have to rely on very little information about how much generation they are receiving from distributed PV generators. In addition, the data can become corrupted or lost due to meter failures and communication errors. In the smart monitoring context, the unwillingness of residential customers to share their PV readings is another reason that results in incomplete generation profiles. As a consequence, accurate estimation of unobserved data, in addition to cleansing corrupted profiles, is of paramount importance for performing grid monitoring and management tasks. Hence, a reduced order estimation technique will allow for the accurate knowledge about the instantaneous aggregated distributed PV power output at any time, using measurements from a small subset of the total number of generators. In this regard, we utilise the recursive estimation structure of the Kalman filter [3] due to it's proven and widespread applications among others in Fuzzy logic-based system modelling [4,5], dynamic neural networks [6] and biomedical signal processing [7]. In the context of observing distributed PV generation, the high resolution time series data can be modelled in the form of equal resolution state dynamics, while the measurement can be obtained at a much lower rate through decimation-a down-sampling method widely used in multi-rate signal processing [8]. The measurement updates to the Kalman filter can be 'event-driven', i.e. the recursive estimator receives the measurement only when it deviates from previous measurements beyond a pre-specified threshold [9,10]. Consequently, PV measurements can be recorded at low sampling rates and the temporal dimension of the data is reduced. Furthermore, A key characteristic of distributed PV is its high level of spatial correlation due to geographical proximity. This feature can be exploited by lowering the spatial dimension of measurements required for full system observability. In particular, the concept of Kriging utilises the spatial correlation in estimating PV generation at unobserved locations [11]. Thus, measurements are collected from a few nearby locations and used to estimate generation profile of all the distributed PV systems in geographic proximity.
In this paper, we utilise such a Multi-Rate and Event-DRIven Kalman Kriging (MREDRIKK) filter that can recursively estimate distributed and behind-the-meter solar generation in scalable spacetime. From the perspective of PV systems, we consider the time series of solar power index (SPI) [12] as the measure and to develop the necessary state-space form. The Oahu Solar Database (described in Section 3) is used to develop a data set of daily SPI time series. The distinctive weather patterns for the entire data set have been identified with the help of k-means clustering. Approximately 50% of this SPI data is used to fit the weatherdependent state-space and spatial model for the MREDRIKK filter while the rest is used to investigate the filter performance.

Literature review
There are two important aspects in the design of a scalable, recursive spatiotemporal PV generation estimation routine: (a) preprocessing of raw time series and (b) development of the necessary state-space form. In this regard, the following prior works on solar irradiance, generation and other weather parameters are important to highlight.
The raw time series of irradiance can be normalised with respect to the clear sky value [13]. However, authors in [14] take additive model approach and decompose the raw measurement into: (a) large scale variability (the trend) and (b) small scale variability (the error). The latter is assumed to be an unobserved state process, which is modelled as a pth order autoregressive (AR(p)) model. The fitted model is used for the Kalman filtering and Kriging to predict weekly and yearly solar radiation. However, the work needs to be expanded to assess Kriging performance under varying dimension of observation. Authors in [15] used Fourier series to detrend the global horizontal irradiance (GHI) data (assuming an additive model); the seasonal component is based on diurnal patterns. The AR (2) and Lucheroni models are then combined to produce 1 h ahead forecasts. The combination is based on the positivity test of the difference of deseasoned data and does not involve dynamic estimation using Kalman filter or spatial estimation through Kriging. In [16], solar irradiance is assumed to be an additive model to predict hourly or subhourly PV power. The deterministic component is approximated through a sinusoidal equation. Power spectral analysis is then used to find the AR moving average model parameters, and then the state-space model for Kalman filtering. Data for one day is used for model fitting and tested on another day of similar generation pattern. This work, however, demands more comprehensive effort in capturing seasonal PV generation patterns. Authors in [17] use a fusion approach for the regressor and Kalman filter prediction outputs, i.e. short-term 10-20 min horizon nowcasting predictions. One month of solar irradiance data are used to obtain regressor input feature vectors and ramp-down event forecasting through a support vector machine. The ramp-down event forecasting gives the adjusting factor within the state transition matrix, which is used in the Kalman filter. However, the fusion alternatives are data dependent and seem inconclusive. In [18], a third order polynomial is used to define the state-space model, although the selection of an identity matrix to define state transition lacks sufficient justification. Unlike the recursive approach of Kalman filters, ratio methods are used in [19] to estimate 1 min resolution PV power given the PV power recorded at a nearby PV system.
In forming the observation space, a historical data-driven approach is investigated to model the observation matrix and state vector. In this regard, Gaussian weight kernels [20], spline backfitted kernels [21], and also a probabilistic approach [22] are reported in capturing the mean dynamics of states over a large volume of measurements. Authors in [20], apply Gaussian weight kernels to fit observation model and assume 'random walk' state evolution to make Kriged forecast of quarterly mean rainfall and annual mean temperature. However, as with [14], the work needs to be expanded to assess Kriging performance under varying dimension of observation. In [21], space-time and time-space models are developed to fit 1 min resolution GHI data where the spline back-fitted kernels are used to find the model coefficients. The time-dependent state vectors obtained can then be used to train the system evolution model through vector AR analysis [23]. This work, like [16], focuses on day-specific modelling and demands for more comprehensive effort in capturing seasonal irradiance patterns. The pair copula-based versatile probability model is used in [22] to present 15 min resolution PV generation for probabilisitc load flow study. However, the consistency in spatial resolution for PV generation data and the power system test case needs to be addressed.
In summary, the prior methods involve the usage of standard time series models and machine learning to forecast solar irradiance, solar generation, temperature and rainfall. An additive modelling approach is evident in all of these works with the absence of any normalisation efforts. Furthermore, only one work involves both the Kalman filtering and Kriging, while rest of them either use Kalman filter or Kriging or none of them. Finally, none of the prior works addresses the issue of capturing weatherdependent variability in temporal evolution as well as effort in temporal down-scaling for measurement access. On the other hand, it is important to design a real-time, weather-dependent and scalable spatiotemporal estimation platform that not only supports limited communication bandwidth through lesser data aggregation but also provides better visibility to the spatially distributed grid contribution of behind-the-meter PV systems. In this regard, we summarise the contribution of proposed method • Develop a generic procedure of utilising available information to train and implement a scalable spatiotemporal estimation of distributed PV generation. The procedure is validated using publicly available and highly variable irradiance database. • Incorporate weather-dependent solar generation patterns within the dynamic estimation algorithm. • Utilise the multi-rate and event-driven mechanism to collect accessible and available PV power measurements. Exploit spatial correlation to estimate PV generation from unobserved locations. This approach reduces utilities' dependence on PV developers' data, preserve user privacy, and can also provide expanded monitoring of ever-growing distributed PVs without additional installation of PV production meters. • Illustrate estimation performance sensitivity to various parameters. This will serve as guidelines for tuning-parameter adjustment in achieving desired estimation performance.
The rest of the study is organised as follows: Section 2 describes the steps in designing the MREDRIKK filter and provides the algorithm. The solar irradiance database used in this research is described in Section 3. The steps in fitting the state-space and spatial model are described in Section 4. The estimation performance is discussed in Section 5. Finally, the findings and future work are summarised in Section 6.

MREDRIKK filter design
In this section, we describe the state-space model formation from a generic AR time series model and introduce the multi-rate and event-driven approach for low-dimensional capturing of temporal measurements. Subsequently, the spatial Kriging will be discussed and the MREDRIKK algorithm will be described for low dimensional spatiotemporal estimation of distributed PV generation.

System model
We consider a scalar continuous-time phenomena observed by N sensors collocated with geographically distributed PV systems at S = s 1 Here μ t is an intercept term reflecting the mean of the time series [24], Defining the p t -dimensional vectors, w t = w t , 0, …, 0 and we have the following model of system dynamics Here and C w, t denotes the covariance matrix of w t -the first diagonal being σ w, t 2 while the rest of entries are zero. We obtain the state covariance matrix C x, t as the solution of the Lyapunov equation Now, let O denote the set of m active sensors from which measurements are being collected, i.e. O ⊂ S and cardinality of O is m. Therefore, the set of unobserved sensors, U = S∖O. Now, the collective observation model for the active sensors here the m × p t observation matrix H t (O) consists of the set of observation coefficients {h t (s r ) r ∈ O} in the first column while rest of the entries are zero. In the same way, the matrix H t (U) is defined for the set of unobserved sensors. The measurement noise covariance matrix C v, t is diagonal in nature indicating the uncorrelated noise variance for the spatial measurements. The parameters of the state-space model just described are time dependent and hence the subscript for discrete time 't' is included. Such dependence comes from the time varying weather pattern that constitutes itself by switching from one pattern to another after certain time period. In this research, the weather patterns are captured by a set of disjoint clusters. The frame length of the time series forming the clusters is determined based on the training data sets' performance in fitting to the AR model. The detailed procedure of finding the number of clusters, frame length and cluster specific model parameter estimation is described in Section 4. Measurements following the observation model (3) are used to estimate the state in the minimum mean squared error (MMSE) sense for a distributed PV system described by the above model. Accordingly, the state estimate at discrete time t is The corresponding error covariance matrix is We use this estimate and the spatial covariance information to estimate unobserved PV systems through the Kriging step (Section 2.3).

Multi-rate and event driven approach
The temporal down-scaling is achieved by incorporating the multirate and event driven approach in collecting measurements from active sensors (i.e. y t (O)). In this regard, we utilise a D-fold decimation [8] to define the following indicator function The event-driven approach to collecting measurements is based upon a pre-determined disparity threshold δ [10] and is prioritised over the D-fold decimation [8]. As a consequence, we introduce a binary factor ϕ(t) in order to reflect the combined (multi-rate and event-driven) process of measurement collection therefore, ϕ(t) indicates the measurement collection status for an active sensor. Here, 0 < θ ≤ 1 is a multiplying factor weighing the sensitivity of the event detection mechanism. The detection becomes more sensitive as the value of θ decreases. Such design is motivated by the fact that the stability of a dynamic estimator is not affected by the choice of event-triggering threshold δ [25]. Now, during the sub-sampling interval [iD, (i + 1)D); i ∈ ℤ + , measurement from an active sensor is collected only if there is any significant event, otherwise it is estimated using the Kalman filter and the most recent measurement (i.e. collected at t = iD). This approach is carried out in two steps. First, at any discrete time t(iD ≤ t < (i + 1)D; i ∈ ℤ + ), the most recent measurement from active sensors is recorded as which is then used with the MMSE estimate x t t in order to correct measurements for the active sensor It is worthwhile to mention that the concept of multi-rate and event-driven approach was first introduced in [9] to estimate the range map for a low-altitude helicopter with significant computational reduction. However, the persistence based measurement update and correction proposed in (8) and (9) is not considered. Consequently, the estimation error covariance matrix and Kalman gain based estimation correction take place only when an actual measurement is collected. On the contrary, the proposed approach in (8) and (9) allow the occurrence of such recursive steps independent of measurement collection status inside the MREDRIKK algorithm (discussed in Section 2.4).

Kriging step
The concept of Kriging is about the utilisation of spatial covariance to make linear MMSE estimates at locations from where no measurements are recorded [11,26,27]. As a consequence, Kriging facilitates the spatial downscaling as well as provides estimates for any malfunctioning sensor's measurement. In this research, we establish the simple Kriging of unobserved sensors based on the active sensors' measurement update. Let C r, k denote the spatial covariance between sensor locations s r and s k . Now, given the unobserved sensor location s k , the Kriging weights {α r, k } are obtained by solving the following equation The spatial covariances, especially those from unobserved sensor locations are obtained from a standard semivariogram model fitted to historical data. In this research, we fit the exponential model to the semivariogram between locations s r and s k [28] Here γ(d) = 0.5 × variance y(s r ) ∼ y(s k ) and d is the Haversine distance between locations s r and s k in the geographic coordinate system. The choice of exponential fitting is motivated by the exponentially decaying spatial variability of solar generation [29].
Once, the semivariogram model is fitted, the spatial covariance is then calculated as, C r, k = c 0 + c 1 − γ(d) and used in (10) to calculate the Kriging weights. Therefore, the estimate for the unobserved sensor k Here A = [α : , 1 , …, α : , N − m ] ⊤ . In the next subsection, we constitute the spatiotemporal dynamic estimation algorithm that entails the concepts described so far.

MREDRIKK algorithm
We incorporate the spatial (Kriging) and temporal (multi-rate and event-driven) downscaling approach within the recursive structure of Kalman filtering and develop Algorithm 1 (Fig. 1).
Step 1 initialies the filtering, and steps 2-7 represent the Kalman filter. The multi-rate and event-driven measurement collection is represented in step 6.
Step 8 represents active sensor measurement correction if necessary based on step 6. Finally, Kriging is performed in step 9 for the unobserved locations. x is a p tdimensional state vector. Correspondingly, the error covariance matrix M, state covariance matrix C x , state transition matrix F, and the process noise covariance matrix C w are all of dimension p t × p t . These matrices are defined in Section 2.1 and associated weatherdependent coefficients are found according to the method described in . The relation between known data, and the models and algorithm described in this section are illustrated in Fig. 2. To summarise, the proposed algorithm adds novelty to the recursive structure of Kalman filter by integrating: (i) time varying weather pattern through distinctive set of parameters for state-space and observation models, (ii) temporal downscaling through an upgraded multi-rate and event-driven approach and (iii) spatial downscaling through Kriging. In the following section, we describe the database used to fit various parameters of the MREDRIKK algorithm and to analyse filter performance.

Oahu solar database
The Oahu Database consists of daily GHI data collected from 18 irradiance sensors [30]. The database covers one and a half year of observation (March 2010 through October 2011). Fig. 3 shows the geographic distribution of these sensors in the vicinity of Kalaeloa airport, Oahu, Hawaii -one of the very high solar variability zones [31]. The sensor nearest to the geo-centroid of this layout is the 'AP01', marked in Fig. 3. All sensors except 'AP02' measure GHI at 1 s resolution. However, the 'AP02' sensor collects GHI, temperature, pressure and wind speed at 1 min resolution. Thus, for consistency, the 1 s resolution data set is averaged to 1 min making the sampling interval, T = 1 min. The maximum separation is 0.7 miles between 'DH01' and 'AP06' indicating a circular area of ∼0.4 square miles. The number of residences covered by this would be around 1000 and 500 in a dense urban and suburban neighbourhood, respectively [32]. The 'REST 2' model [33] based half hourly clear sky GHI is then retrieved at a nearby location from the NSRDB database [34]. With the help of cubic spline interpolation, a 1 min resolution clear sky GHI is then obtained for each day. Since, the goal is to estimate PV generation, we use these spatial time series of irradiance and weather data in synthesising the PV inverter AC power output profile at the specified locations of N = 18 sensors. The longitude-latitude pair of these locations form the set S. At each location, the AC power output from a 2.5 kW inverter (with inverter loading ratio of 1.48 [35]) is calculated by feeding the Oahu Database into the PVLIB toolbox [36]. For 500 residences this would represent a 1.25 MW distributed PV in the dense suburban neighbourhood. Such power profiles obtained are then used to calculate the spatial time series of SPI SPI(s r , t) = P(s r , t) P cs (s r , t) ; s r ∈ S Here, SPI stands for the SPI -a ratio of actual P and clear sky AC power P cs [12], and is illustrated in Fig. 4 for the sensor 'AP01'. The formation of daily SPI time series is repeated for seven months from April 2010 to October 2010. For each day, the time horizon is kept between 7:07 AM and 5:06 PM totaling 600 time stamps at 1 min resolution. The particular time horizon is chosen so that it cuts off any spurious SPI value caused by the diminishing clear sky AC power in the denominator during sunrise and sunset periods. The aggregate SPI series thus obtained is named as the training data set and used in the next section to fit the system model. It consists of 214 realisations (i.e. 214 days of seven months) of SPI time series at 18 locations. Similarly the test data set is formed for April 2011 through October 2011, which is used in Section 5 to investigate the performance of the fitted model and MREDRIKK algorithm. The procedure of creating the training and test database from known information is summarised in Fig. 2, which also illustrates how these database are being used in the upcoming sections for model fitting, algorithm implementation and performance assessment. It should be noted that, the proposed approach leverages two tuning parameters, which can be adjusted for desired performance while only three types of information are required as input to implement the algorithm. No additional inputs are required for internal data processing to fit the associated statespace model.

Fitting of state-space model
From the training data set, we pick the 214 realisations of SPI(s AP01 , t) to estimate various parameters of the state-space model. This choice is motivated by the closest vicinity of 'AP01' to the geo-centroid (Fig. 3). Each of these highly variable and nonstationary SPI time series is then broken into (600/l frame ) time frames for a given frame length, l frame min. Thus, a total of (214 × 600/l frame ) time frames are used to perform cluster analysis. The clustering is performed to identify distinctive weather patterns as represented by different clusters. In this regard, the squared Euclidean distance based k-means + + algorithm is chosen to perform unsupervised clustering of the entire set of time frames just formed. Such choice is motivated by its proven application in weather-dependent solar irradiance classification [37]. Now, due to the random initialisation, the algorithm is run multiple times and the clusters with the lowest within-cluster sums of point-tocentroid distances are identified [38]. For a given number of clusters being sought, the performance of k-means clustering is analysed by calculating the average silhouette coefficient (ASC)a measure of a time frame's similarity to its own cluster as compared to other clusters. The value of ASC lies in between −1 and + 1 and closer to the later indicates higher quality of clustering performed [39]. The entire process of k-means clustering, and ASC calculation are repeated for different number of clusters. Furthermore, the frame length l frame min is varied to create new set of time frames and fed to the clustering and ASC calculation process. It has been found that the ASC value becomes the highest for partitioning into two clusters, irrespective of the frame length, l frame min. This indicates the existence of two possible weather patterns for the data set under consideration -a direct consequence of continuous intermittent cloud movement over such a high solar variability zone [31]. Consequently, we continue with two clusters to fit various models in the following subsections.

Fitting of system dynamics
We denote the set of parameters to be estimated for the dynamic model (1) as β = μ t , a t, 1 , a t, 2 , a t, 3 , …, a t, p t − 1 , a t, p t . The partitioning into two clusters is repeated by varying the length of time frames and the corresponding ASC is shown in Fig. 5a. For each cluster and frame length, the ARFIT toolbox [40] is used to fit an AR model to the training data set while capturing the approximate 95% confidence intervals for the estimated parameters [24]. The set of intervals is denoted by β 95 and is normalised with respect to estimated parameters to calculate the metric, γ = 100 × β 95 / β. Fig. 5b shows the maximum and minimum values of γ for each cluster. The narrowest range of confidence interval is observed for both clusters when the frame length is set to 12 min. Based on these observations, l frame is set to 12 min and both the training and test data set are broken into 12 min time frames and two clusters are obtained for each data set. Intuitively, the weather switches between two possible patterns after every 12 min and the subscript t in the system model and Algorithm 2.4 captures this time varying behaviour. The clustered 12 min long training data set is used to fit the models (1), (2) and the cluster labels of the test data set are used as weather forecast while investigating the performance of proposed algorithm in Section 5.

Fitting of observation model
The observation model (3) is fitted to the clustered training data set through least square fitting. For each cluster C, a concatenated time series is formed as a row vector using all the 12 min long cluster members. We denote its transpose as ℛ(s AP01 , C). At the second step, the spatial observation from the remaining 17 sensors are broken into 12 min frames and labelled with the cluster of concurrent 12 min frame of SPI(s AP01 , t). Thus, the set of 12 min SPI series labelled as cluster C is denoted by ℬ(s r , C) for the sensor located at s r . The concatenation of all series in ℬ(s r , C) results in a row vector, whose transpose is denoted as b(s r , C). Therefore, the least square estimate of the observation coefficients for cluster C is given by [41] h^(s r , C) = b(s r , C) ⊤ ℛ(s AP01 , C) ℛ(s AP01 , C) ⊤ ℛ(s AP01 , C) The measurement noise variance is then modelled as the error variance of this fitting With the prior knowledge about the set of active sensors, the observation matrices H t (O), H t (U) and noise covariance matrix C v, t are then formed using the fitted coefficients

Calculating threshold δ
We utilise the training data set to calculate the disparity threshold for each location s r and cluster C as follows: The threshold values are then chosen according to the cluster type for any given 12 min frame. It should be noted that the threshold calculation in (16) is rather comprehensive as compared to previous work of determining the event detection threshold in [10].

Semivariogram fitting
The sample semivariograms for each cluster C are calculated using the series set {b(s r , C), s r ∈ S}. In fitting the exponential model (11), we vary the sum (c 0 + c 1 ) and use the calculated semivariograms to find the lowest fitting error. Fig. 6 shows the fitting performance for both clusters and suggests the following choice for (c 0 + c 1 ) c 0 + c 1 = 0.12; Cluster1 0.08; Cluster2.
Using these values, we obtain the exponential fitting of sample semivariograms of each clusters, which is illustrated in Fig. 7. Table 1 shows the fitted coefficients rounded to two decimal places. These coefficients are then used to calculate the spatial covariance matrix for simple Kriging expressed in (12).

MREDRIKK performance analysis
The fitted state-space models are applied to the MREDRIKK filter according to the cluster sequence of the test data set (April 2011 through October 2011). As mentioned in Section 4.1, for any time frame ( j − 1)l frame + 1 ≤ t ≤ jl frame ; j ∈ ℤ + , the associated model parameters tagged as Cluster 1(or 2) are selected, if the jth cluster is 1(or 2). As time evolves, this parameter selection for Algorithm 1 is repeated every 12 min -the value of l frame set for the database being used. With such formation, Fig. 8 illustrates the estimation and Kriging of SPIs for one of the days within the test data set, at locations of 'DH07' and 'DH11', respectively. In this example, 'DH07' represents one of the nine active sensors and 'DH11' The above metric represents performance of estimation when r ∈ O and Kriging when r ∈ U. The spatial average of MAE is then obtained for each day while varying the sets O and U. The cardinality of O is chosen from {3,6,9,12,15}, while U is set as (17 − O ). In the following subsections, we utilise the cumulative probability distributions (CDF) of these daily average MAEs of the entire test database to showcase the effect of active sensor volume O , decimation factor D and event detection sensitivity factor θ.
The empirical CDFs are obtained by creating 64-bin histograms of the MAEs.

Effect of active sensor volume
The estimation and Kriging performance for varying size of active sensor set is shown in Fig. 9a. For this investigation, D and θ is set to 5 and 1, respectively. It is observed that the estimation and Kriging performance improve as more sensors are active while less sensors are needed to be Kriged. The increase in the number of active sensors leads to a better observability of state in (3) and hence better estimation for the active sensors can be achieved.
Consequently, Kriging has to be done for fewer sensors from improved estimates, resulting in better Kriging performance.

Effect of event detection
The effect of event detection sensitivity is shown in Fig. 9c by varying θ while keeping the decimation factor and active sensor volume at D = 15 and O = 9, respectively. It can be observed that both the estimation and Kriging performance improve as the event detection becomes more sensitive (i.e. θ is reduced). According to (7), a smaller value of θ causes the detection of greater number of SPI significant events by overriding decimation in measurement update. As a result, the estimation for the active sensors improves by capturing more frequent measurements leading to further betterment in Kriging. However, a careful observation of Algorithm 1 reveals that unlike Kriging (Step 9), the estimation for active sensors (Step 8) have direct impact from the measurement collection status ϕ(t). This characteristic can be observed in Fig. 9c as the improvement in estimation for active sensors is more significant than that for Kriging. It is expected that, the number of SPI events detected as significant would increase when the event detection mechanism becomes more sensitive. However, the observability of active sensors can mitigate such detection based measurement augmentation. The box-plots in Fig. 10 show that this phenomena where the number of significant events reduces with the increase of active sensors. This is a direct consequence of the state observability, which is enhanced as more sensors are actively  • A zero value for KSTS in the last column of Table 2 indicates consistent estimation performance when at least six sensors are active. This is a direct consequence of achieving the best estimation performance (i.e. average MAE reaches zero) when θ ≤ 0.1. Fig. 9c Fig. 9c about more significant variation in estimation performance compared to Kriging. • While the decimation has similar impact on both estimation and Kriging, the rate of performance degradation drops with the increase of decimation.
To summarise, the estimation and Kriging performance vary significantly with respect to active sensor volume, decimation as well as event detection threshold, which are validated through the KSTS values.

MREDRIKK versus persistent estimation
In this section, we compare the performance of MREDRIKK filter with the persistent estimation of distributed PV generation. The concept of 'persistence' is widely used to assess performance of solar forecasting methods [12], which motivates the derivation of persistent estimation in this research. This applies both to the active sensor measurement collection as well as estimation for unobserved sensors. For the active sensors, the persistence implies no change in the measurement during the decimation interval unless there is a significant event, and hence no measurement correction is required (in contrast to step 8 of Algorithm 1). Thus, from (8), the persistence in active sensor measurement is given by For the unobserved sensors, the persistence implies that the leastsquare fitted observation coefficients are sufficient to capture distributed PV generation from the intermediary states (in contrast to step 9 of Algorithm 1). Accordingly, (12) is modified as follows to give persistent estimation at unobserved sensors Now, we obtain persistent estimation algorithm by replacing steps 8 and 9 of Algorithm 1 by (19) and (20), respectively. Next the performance of MREDRIKK filter is compared with persistent estimation. In this regard, the global average MAE is obtained from the daily spatial MAEs (pre-calculated using (18)) for the entire test database. The global performance in estimation and Kriging thus obtained are shown in Fig. 11 for different size of active sensor set. It can be observed, that the MMSE estimation performs better when the filter receives measurement from at least nine sensors. On the other hand, the Kriging performs better irrespective of the number of sensors when the event detection threshold is set to θ = 0.01.

Conclusion
A MREDRIKK filter is utilised for dynamic estimation of distributed PV generation in this research. The scalable estimation is performed in temporal domain through the multi-rate and eventdriven measurement update, whereas the Kriging enables a scalable spatial estimation. A time series of SPI is used to define the behind-the-meter PV states for a set of irradiance sensors located in the proximity corresponding to a high PV penetrated residential neighbourhood. The weather-dependent state-space and semivariogram model have been fitted to publicly available high resolution irradiance database and the performance of the MREDRIKK filter has been thoroughly investigated by varying (i) the size of active sensor set, (ii) decimation factor and (iii) event detection sensitivity. The MREDRIKK filter shows improvement in both estimation and Kriging with the increase of active sensors while the Kriging is found to be less sensitive to event detection mechanism as compared to the decimation. The enhanced observability from increased number of active sensors is further illustrated through the reduction in significant SPI events irrespective of the sensitivity of event detection. The proposed filter would not only reduce utilities' dependence on PV developers provided costly data service, or avoid incentives to access residential customers' PV readings, but would also enable the best utilisation of available information in maintaining voltage stability  and reliable grid. In future work, we plan to further investigate and compare the performance of MREDRIKK filter with a much larger spatial irradiance database of comparable resolution and variability. Furthermore, the data-driven design of the MREDRIKK filter for distributed PV generation estimation will be extended to a distribution feeder scale with the incorporation of standard decentralised optimisation approaches.