Unsupervised Disaggregation of PhotoVoltaic Production from Composite Power Flow Measurements of Heterogeneous Prosumers

TL;DR

This paper develops and compares four algorithms to estimate unobserved photovoltaic generation and demand from aggregated power flow data and irradiance measurements, validated with real-world data from a heterogeneous prosumer setup.

Contribution

It introduces four novel estimation algorithms that leverage irradiance models and demand assumptions to disaggregate PV and demand in power flow measurements.

Findings

Algorithms accurately estimate PV generation in real-world data.

Using demand models improves disaggregation performance.

PV generation patterns are explainable through irradiance and installation configurations.

Abstract

We consider the problem of estimating the unobserved amount of photovoltaic (PV) generation and demand in a power distribution network starting from measurements of the aggregated power flow at the point of common coupling (PCC) and local global horizontal irradiance (GHI). The estimation principle relies on modeling the PV generation as a function of the measured GHI, enabling the identification of PV production patterns in the aggregated power flow measurements. Four estimation algorithms are proposed: the first assumes that variability in the aggregated PV generation is given by variations of PV generation, the next two use a model of the demand to improve estimation performance, and the fourth assumes that, in a certain frequency range, the aggregated power flow is dominated by PV generation dynamics. These algorithms leverage irradiance transposition models to explore several azimuth/tilt configurations and explain PV generation patterns from multiple plants with non-uniform installation characteristics. Their estimation performance is compared and validated with measurements from a real-life setup including 4 houses with rooftop PV installations and battery systems for PV self-consumption.

Tables

Table 1. TABLE I: PV and battery systems in the experimental setup

House ID	PV capacity (kWh)	Azimuth	Tilt ($°$)	Distance from House 1 (m)	Battery rating (kVA/kWh)
1	10.0	95	14	0	3/8.8
2(a)	7.2	187	36	100	3/4.4
2(b)	3.5	266	40	100	–
3	8.0	187	40	260	3/8.8
4	6.6	180	24	170	3/4.4

Table 2. TABLE II: Estimation performance statistics (%)

Statistic	A	B	C	D
nRMSE with self-consumption
min	5.22	5.2	4.0	4.3
max	20.6	17.8	10.4	8.2
mean	9.7	9.2	7.0	5.5
median	8.4	8.4	7.0	5.4
nRMSE without self-consumption
min	3.5	4.4	3.4	4.1
max	17.9	17.7	9.5	8.8
mean	6.6	7.62	5.3	4.6
median	4.6	4.8	4.7	4.4
nMAE with self-consumption
min	3.3	3.3	2.5	2.7
max	14.0	11.8	8.5	5.1
mean	6.2	5.9	4.7	3.4
median	5.3	5.3	4.4	3.4
nMAE without self-consumption
min	2.2	2.8	2.3	2.5
max	12.1	11.8	5.4	5.6
mean	4.3	4.9	3.2	2.9
median	3.0	3.0	3.0	2.7
nME with self-consumption
min	-13.9	-11.7	-8.1	-4.3
max	-1.8	-1.8	8	4.3
mean	-5.7	-5.2	2.3	-1.0
median	-4.7	-4.8	-2.8	-1.3
nME without self-consumption
min	-11.9	-1.2	-4.8	-5.0
max	1.9	-1.1	4.6	4.0
mean	-2.2	-3.8	0.8	0.5
median	-1.2	-1.7	-1.0	0.2

Table 3. TABLE III: Estimation nRMSE (%) at different levels of PV penetration with self-consumption.

Nominal PV Capacity (kWp)	PV Penetration (%)	A	B	C	D
35.3	294	5.2	5.2	4.0	4.3
17.6	146	9.4	6.9	9.7	4.6
8.8	73	18.7	10.4	11.8	5.2

Equations

$${\hat{G}}_k(\mathbf{\alpha })=\sum _{j=1}^J{\alpha }_j\cdot I_{jk}^{╲},$$

$$\displaystyle\begin{aligned} \boldsymbol{\alpha}^{o}=\underset{\boldsymbol{\alpha}\in\mathbf{R}^{J}_{+}}{\operatorname{arg}\,\operatorname{min}}\;\Bigg{\{}\sum_{k=1}^{K}\bigg{|}&\left(P_{k}-P_{k-1}\right)+\\ -&\left(\widehat{G}_{k}(\boldsymbol{\alpha})-\widehat{G}_{k-1}(\boldsymbol{\alpha})\right)\bigg{|}\Bigg{\}}\end{aligned}$$

$${\hat{G}}_k(\mathbf{\alpha })=\sum _{j=1}^J{\alpha }_j\cdot I_{jk}^{╲},$$

$${\hat{P}}_k(\mathbf{\alpha },\hat{\mathbf{L}})={\hat{L}}_k-{\hat{G}}_k(\mathbf{\alpha }),$$

$$\displaystyle\begin{aligned} \begin{bmatrix}\boldsymbol{\alpha}^{o}\\ \widehat{\boldsymbol{L}}^{o}\end{bmatrix}=\underset{{\scriptsize\begin{bmatrix}\boldsymbol{\alpha}\in\mathbf{R}^{J}_{+}\\ \widehat{\boldsymbol{L}}\in\mathbf{R}^{K}_{+}\end{bmatrix}}}{\operatorname{arg}\,\operatorname{min}}\;\Bigg{\{}&\sum_{k=1}^{K}\left(P_{k}-\widehat{P}_{k}(\boldsymbol{\alpha},\widehat{\boldsymbol{L}})\right)^{2}+\\ +\lambda&\sum_{k=1}^{K}\left|\widehat{L}_{k}-\widehat{L}_{k-1}\right|\Bigg{\}}\end{aligned}$$

$${\hat{P}}_k(\mathbf{\alpha },\hat{\mathbf{L}})={\hat{L}}_k-\sum _{j=1}^J{\alpha }_j\cdot I_{jk}^{╲},$$

$${\hat{L}}_1={\hat{L}}_2=\cdots ={\hat{L}}_c,$$

$${\hat{L}}_{c(i-1)+1}=\cdots ={\hat{L}}_{c(i-1)+c}$$

$$\left[\begin{array}{c}{\mathbf{\alpha }}^o\\ {\hat{\mathbf{L}}}^o\end{array}\right]=\underset{\left[\begin{array}{c}\mathbf{\alpha }\in {\mathbf{R}}_{+}^J\\ \hat{\mathbf{L}}\in {\mathbf{R}}_{+}^K\end{array}\right]}{\mathrm{argmin}}\left\{\sum _{k=1}^K{\left(P_k-{\hat{P}}_k(\mathbf{\alpha },\hat{\mathbf{L}})\right)}^2\right\}$$

$${\hat{P}}_k(\mathbf{\alpha },\hat{\mathbf{L}})={\hat{L}}_k-\sum _{j=1}^J{\alpha }_j\cdot I_{jk}^{╲},$$

$${\hat{L}}_{c(i-1)+1}=\cdots ={\hat{L}}_{c(i-1)+c}$$

$${\mathbf{\alpha }}^o=\underset{\mathbf{\alpha }\in {\mathbf{R}}_{+}^J}{\mathrm{argmin}}\left\{\sum _{k=1}^K\rho \left({\mathcal{P}}_k-G_k\left(\mathbf{\alpha }\right)\right)\right\}$$

$${\hat{G}}_k(\mathbf{\alpha })=\sum _{j=1}^J{\alpha }_j\cdot {\mathcal{I}}_{jk}^{╲},$$

$$I_{jk}^{╲}={\mathcal{I}}_{jk}^{╲}\left[1+\gamma (T_{\text{cell},k}-T_{\text{ref}})\right]$$

$$T_{cell,k}=T_{\text{air}}+\beta {\mathcal{I}}_{jk}^{╲},$$

$$\hat{\mathbf{P}}=\hat{\mathbf{L}}-M\mathbf{\alpha }=\left(\begin{array}{cc}{1}_{K\times K},& -M\end{array}\right)\left(\begin{array}{c}\hat{\mathbf{L}}\\ \mathbf{\alpha }\end{array}\right)=S\mathbf{x}$$

$$J$$

$$={\mathbf{x}}^T{S}^{T}S\mathbf{x}-2{\mathbf{P}}^{T}S\mathbf{x}=\frac{1}{2}{\mathbf{x}}^{T}H\mathbf{x}-{\mathbf{f}}^T\mathbf{x},$$

Full text

Nomenclature and Acronyms

PV

Photovoltaic.

GHI

Global Horizontal Irradiance.

GNI

Global Normal Irradiance.

nRMS

Normalized Root Mean Square Error.

PCC

Point of Common Coupling.

MPPT

Maximum Power Point Tracking.

$k$

Discrete time index.

$P_{k}$

Active power flow at the PCC at time interval $k$.

$\mathcal{I}^{\text{\----}}_{k}$

GHI observation at time $k$.

$j$

Index for panel tilt/azimuth configuration.

$I^{\diagdown}_{jk}$

Estimated GNI corrected for temperature for configuration $j$ at time $k$.

$\widehat{G}_{k}$

Estimated PV production at time interval $k$.

$\widehat{L}_{k}$

Estimated demand at time interval $k$.

$\alpha_{j}$

PV nominal capacity at configuration $j$.

I Introduction

Incresed levels of distributed photovoltaic (PV) generation determine higher reserve requirements at the system level and violations of voltage and line ampacity constraints in distribution systems during peak production hours [1, 2]. Technical solutions envisaged to mitigate PV generation drawbacks are curtailment strategies, control of converters active/reactive power, PV self-consumption schemes and dispatch of local power flows according to network-safe power consumption trajectories (e.g. [3, 4, 5, 6, 7, 8, 9, 10, 11]). A requirement for the implementation of those strategies is the availability of real-time production measurements from PV facilities. Incidentally, these are also useful to train data-driven prediction models (e.g. [12, 13]). However, such a precondition is not always met in real-life conditions because installations are not always monitored, and, even when they are, factors such as i) privacy concerns, ii) conflicts due to the different owners of the metering infrastructures, and iii) lack of standards for monitoring and aggregation of measurements, and their communication, play against the possibility of collecting real-time PV production measurements.

As an alternative to direct monitoring of PV systems, we consider in this paper the problem of disaggregating PV generation from the aggregated active power measurements of a group of prosumers. The estimation principle relies on modelling PV generation as a function of the global horizontal irradiance (GHI), assumed known from local measurements. Four estimation algorithms are proposed and compared: the first assumes that the variability in the aggregated power flow measurements are mostly given by variations of the PV generation, the second and third leverage a model of the demand to improve estimation performance, and the fourth assumes that there is a certain frequency range in which the aggregated power flow measurement is dominated by PV generation components. All four algorithms use a transposition model to project GHI into a number of pre-defined differently oriented tilted planes to explain production from sites with different configurations. The algorithms are designed to be unsupervised, i.e., they do not require measurements of the PV power profiles to be trained. The algorithms are tested with measurements from a real-life setup of four houses with monitored rooftop PV plants and grid-connected battery systems, enabling the testing of estimation performance even when the demand is correlated with PV generation.

Even if in the existing literature the problem of disaggregation has been extensively investigated for nonintrusive load monitoring (e.g., [14, 15]), its application to PV disaggregation was considered in [16], which develops estimators of the total PV generation using the active power profile of a nearby installation and GHI proxy measurements as explanatory signals. With respect to [16], we leverage PV transposition models to identify PV production patterns from installations with different tilt/azimuth configurations, a key factor in urban contexts where PV generation is from rooftop PV plants and tilt and azimuth configurations are dictated by roofs characteristics and might not be uniform.

The paper is organized as follows: Section II states the problem, III presents the disaggregation algorithms, IV discusses procedures and measurements for validation, V presents and discusses performance, VI summarizes key results and contributions.

II Notation and Problem Statement

II-A Configuration of the system

We consider a feeder with distributed PV production, possibly from installations with different tilt and azimuth, and demand (e.g., Fig. 1). The power injections at the single buses are not measured, however the total prosumption (PV generation + demand) is known thanks to sensing the active power flow111Reactive power is not of special interest since PV plants normally operate at unitary power factor and, more in general, it is not possible to do assumptions on the kind of reactive power control implemented. at the point of common coupling (PCC). Local GHI values are from a pyranometer (although other methods could be considered, e.g. leveraging information available from nearby monitored PV installations [17]). We do the modeling assumption that the PV installations in the area are subject to the same GHI. Local GHI measurements are known to be accurate in a range of 50 meters[18]. Therefore, the proposed algorithms are expected to perform adequately when PV plants are spread over a small area, and their performance to decrease when considering larger areas. Due to the small size of the networks that these methods target, grid transmission losses are neglected at this stage due to the short length of the cables.

II-B Notation

The active power flow measured at the PCC at the discrete time interval $k$ is denoted as $P_{k}$ ($\mathrm{kW}$). Positive flows denote consumption and vice-versa (i.e., passive sign notation). GHI measurements are denoted by $\mathcal{I}^{\text{\----}}_{k}$ in $\mathrm{kW}\text{\,}{\mathrm{m}}^{-2}$, while $I^{\diagdown}_{jk}$ ($\mathrm{kW}\text{\,}{\mathrm{m}}^{-2}$) denotes the estimated global normal irradiance (GNI) to a certain tilted plane $j$ corrected for temperature (as described in III-E), where $j=1,\dots,J$ denotes the plane tilt/azimuth configuration. We consider $J=21$ tilted planes with tilt and azimuth values equally spaced on a south-facing semi-sphere, chosen to have a reasonably representative set of potential configurations of PV installations in the northern hemisphere. GNI estimations are from the Hay-Davies transposition model [19, 20]. The quantity $\widehat{G}_{k}$ and $\widehat{L}_{k}$ ($\mathrm{kW}$) respectively denote the estimated PV production and estimated demand, which are to be determined. A practical example of the disaggregation process is described in the following paragraph.

II-C Problem statement

The problem is estimating the trajectories of the demand and total PV generation from measurements of the active power flow at the PCC and local GHI observations. This is exemplified in Fig. 2a (night time observations are omitted) which shows the inputs, intermediate results and outputs of one among the proposed algorithms. The inputs are the aggregated power flow $P_{k}$ at the PCC (top panel) and GHI (middle panel, solid fill). The middle panel plot in Fig. 2b also shows the GNI trajectories $I^{\diagdown}_{jk}$, used to explore the potential PV production from plants with various tilt/azimuth configurations, as typical in urban feeders where panels are installed according to roof characteristics. Finally, the lower panel plot in Fig. 2c shows the output of one of the proposed algorithms222Method C, parameters $T_{s}=30$ s $c=10$, mean nRMSE 5%., with the estimated demand and estimated PV generation (orange line), the latter close to the measured ground truth PV generation (solid gray fill).

III Disaggregation Algorithms

The estimated global PV generation $\widehat{G}_{k}$ at the PCC is modelled as the sum over all the tilt/azimuth configurations $j=1,\cdots,J$ of the transposed irradiance $I^{\diagdown}_{jk}$ times $J$ nonnegative coefficients $\alpha_{j}\in\mathbf{R}_{+}$:

[TABLE]

where $\boldsymbol{\alpha}=\{\alpha_{j},j=1,\cdots,J\}\in\mathbf{R}^{J}_{+}$ denotes the set of $\alpha_{j}$. Physically, $\alpha_{j}$ is the PV generation capacity installed at configuration $j$. It is measured in kWp, kilowatt peak, and denotes the amount of power produced in standard test conditions (STC, $\text{\,}\mathrm{kW}\text{\,}{\mathrm{m}}^{-2}$ GNI at $20\text{\,}\mathrm{\SIUnitSymbolCelsius}$). The temperature effect on the PV conversion efficiency is accounted for by preprocessing the input GNI time series as described in Section III-E. By modelling the PV generation as in (1), we assume that PV plants operate in the maximum power point tracking (MPPT) mode; in case the output of a PV plant is controlled (i.e. curtailed), it is likely to be monitored and its contribution can be removed from the aggregated power flow, still allowing to apply the algorithms. Partial shading effects is not explicitly modeled, even if, as discussed in the following, some of the proposed methods are robust against it.

As detailed in the following, the proposed disaggregation algorithms estimate $\alpha$ (thus PV generation $\widehat{G}_{k}$) and require a training phase. Four algorithms, denoted as Method A, B, C and D, are discussed. They attempt to estimate PV generation by exploiting different modeling principles inspired by the following empirical considerations:

• •

variations of PV generation dominate the variations of the power flow at the PPC. Method A estimates PV generation by seeking for a trajectory with variations as close as possible to those observed at the PCC. The drawback is that variations are also due to demand changes (e.g. load inrushes);

• •

the power flow at the PCC is modelled as the sum between PV generation and demand, where the latter is described by using a load model. Method B and Method C exploit two different load models, as explained later;

• •

in a certain frequency range, the dynamics of the power flow at the PCC are dominated by those of PV generation. Method D exploits the fact that demand and PV generation have different time dynamics: filtering the power flow measurements at the PCC makes possible to estimated the PV generation.

The validity of these empirical modelling considerations are tested in the results section by assessing and comparing the algorithms performance.

III-A Method A

The unknowns $\boldsymbol{\alpha}$ are determined by assuming that the variability in the observed aggregated power flow are due to variations of the aggregated PV power. This is modeled by minimizing the norm-1 of the difference between the once differentiated time series $P_{k}$ and $\widehat{G}_{k}$ while subject to the estimated total PV production model (1):

[TABLE]

subject to:

[TABLE]

The problem in (2)-(3) is linear, thus convex and tractable. The resolution of the input time series is a parameter of the Method A, and its importance is discussed in the performance assessment. Method A does not allow to model demand dynamics, e.g. load inrushes would be considered as variations of PV generation. The next two methods use a model of the demand to work around this limitation.

III-B Method B

Let $\widehat{\boldsymbol{L}}=[\widehat{L}_{1},\dots,\widehat{L}_{K}]$ be the estimated demand trajectory. The estimated active power flow at the PCC $\widehat{P}_{k}$ is now written as:

[TABLE]

i.e. the sum between the estimated total PV generation and demand, the former with negative sign because corresponding to generation. Method B attempts to determine $\widehat{L}_{k}$ and $\boldsymbol{\alpha}$ by minimizing the norm-2 of the estimation error ${P}_{k}-\widehat{P}_{k}$. However, this problem is under-determined since the $K+J$ free variables are more than the number of observations $K$. Therefore, we augment the cost function and consider the sum of the least square and norm-1 of the once differentiated $\widehat{L}_{k}$ series (i.e., a combined linear regression and trend filtering problem, as for example in [21] [22]):

[TABLE]

subject to:

[TABLE]

The cost function is the sum of a vector norm-1 and a quadratic cost function and convex if the latter term is convex. As shown in Appendix A, the convexity of the quadratic term cannot depends on the input data and can be verified a-priori.

III-C Method C

In Method B, under-determination was solved by minimizing the norm-1 of the demand trajectory $\widehat{\boldsymbol{L}}$. As an alternative, we apply here a piecewise constant model of the demand, i.e. we require the unknown sequence $\widehat{\boldsymbol{L}}$ to be piecewise constant for $c$ consecutive samples, where $c$ is a parameter, by enforcing the following $c-1$ equality constraints:

[TABLE]

for the case of the first $c$ samples. Extending to the set of $K$ measurements ($K$ multiple of $c$) yields:

[TABLE]

Modelling the demand as piecewise constant is a reasonable assumption when the length of the constant segment does not overlap with typical intra-day dynamics of the demand, i.e. for small $c$ values and densely sampled series. In other words, it is reasonable when considering short periods of time (e.g., seconds), when the persistence model of the demand has unbeaten performance, see, e.g., [23, 24]. When the demand has shorter variations than $c$ (e.g., load inrushes), the estimated demand will have the average value of the waveform and the residuals will be the estimation error. The sampling time and $c$ are design parameters: the sensitivity of the algorithm performance with respect to their values is assessed in Section V. Method C consists in minimizing the norm-2 of the estimation error $P_{k}-\widehat{P}_{k}$ subject to the estimated aggregated power flow and piecewise constant demand models:

[TABLE]

subject to:

[TABLE]

The additional constraints (11) are linear and do not impact convexity, thus the same consideration as for Method B applies.

III-D Method D

Method D splits the observed active power flow at the PCC by exploiting similarities between the signals $P_{k}$ and $\widehat{G}_{k}(\boldsymbol{\alpha})$ in a certain frequency range. This approach is motivated by having verified similarities in the spectral density of the measured aggregated power flow and measured PV generation (available from the test site) with the Welch’s periodogram method for coherence analysis [25]. Method D initially filters the input GNI $I^{\diagdown}_{jk}$ and aggregated power flow $P_{k}$ with a sixth order Butterworth band-pass filter, where the low and high cut-off frequencies are parameters that reflect the frequency range where the aggregated power profile and PV generation are similar. They are free parameters and the sensitivity of the algorithm performance to their values is investigated in Section V. Let $\mathcal{P}_{k}$ and $\mathcal{I}^{\diagdown}_{jk}$ respectively denote the above mentioned filtered version of the sequence $P_{k}$ and transposed irradiance $I^{\diagdown}_{jk}$. The vector of unknowns $\boldsymbol{\alpha}^{o}$ is computed by the following robust linear regression:

[TABLE]

subject to:

[TABLE]

where $\rho(\cdot)$ is the bisquare loss function, see [26], a nonconvex relantionship which gives less weight to extreme values in the cost function to be robust against outliers. The problem (12)-(13) is solved by applying an iterative least square approach with guaranteed covergence [27].

III-E Temperature correction

To account for the dependency between PV conversion efficiency and temperature, GNI values are corrected to reflect temperature variations from the reference value $T_{\text{ref}}$ ($25\text{\,}\mathrm{\SIUnitSymbolCelsius}$) by using the empirical model proposed in [28]:

[TABLE]

where $\mathcal{I}^{\diagdown}_{jk}$ is the original irradiance value, and $T_{\text{cell},k}$ the cell temperature at time $k$, estimated as[28]

[TABLE]

where $T_{\text{air}}$ is the air temperature, assumed known from local measurements, and $\beta=$3.78\text{\times}{10}^{-2}\text{,} and $\gamma=$-4.3\text{\times}{10}^{-3}\text{\,} are empirical and plant specific values here assigned considering typical average values333$\beta$ is the average of values for the close roof mount and open rack configurations from [28] , and $\gamma$ the average of the values for polycrystalline modules from [29]..

IV Methods for performance evaluation

The performance assessment of the disaggregation algorithms is performed by evaluating their ability to reconstruct the PV generation time series starting from the measurements of the power flow at the PCC and GHI. Let $e_{k}=G_{k}-\widehat{G}_{k}$ ($G_{k}$ and $\widehat{G}_{k}$ are the PV generation ground-truth value and estimation, respectively) be the estimation error. The performance metrics are: normalized root mean square of the estimation error (nRMSE) $\left(\frac{1}{G}\frac{1}{K}\sum_{k=1}^{K}e_{k}^{2}\right)^{1/2}$, normalized mean absolute error (nMAE) $\frac{1}{G}\frac{1}{K}\sum_{k=1}^{K}|e_{k}|$, and normalized mean error (nME) $\frac{1}{G}\frac{1}{K}\sum_{k=1}^{K}e_{k}$, where $G$ is the total installed PV capacity (35.3 kWp) and $K$ is the samples number in the testing data set.

IV-A Data sets for training and testing

Time series are GHI, power flow at the PCC and PV generation measurements for 1 year from the real-life test facility described in the next paragraph. The first 2 are used for the training, while the last as the ground truth value to test the estimation performance. To preserve daily dynamics of the signals, time series are divided into daily sequences, then randomly shuffled to avoid to train and test the algorithms on different periods of the year. The series time resolutions are 10, 30, 60, 120, 300, 600 and 900 sec (downsampling by samples average) and are to assess algorithms performance with respect to sampling time. Selected resolutions include those normally implemented in existing metering systems, e.g., 900 s is the resolution of smart meters in Switzerland and Italy. Here, the intent is to verify whether such a sampling time is enough for the proposed application, or if performance would benefit from more densely sampled data. Each of the 7 datasets at different resolution is further split into 3 sub-sequences to perform a three-fold cross-validation, i.e., for each resolution, the algorithms are trained on the first fold and tested on the remaining 2; the process is repeated for all the 7 datasets, each time testing the algorithms on the part of the data which is not used for the training. In total, each algorithm is trained and tested 3 times for each resolution, for a total of 21 training and test runs. Measurements refer to days with a uniform mix of sky conditions: partly cloudy, clear sky and overcast. Algorithms performance is tested both when there are batteries performing PV self-consumption and when not, thus allowing to account for the case when the demand is correlated with PV generation.

IV-B Experimental Setup

Measurements are from a real-life setup in the region of Basel, Switzerland, with 4 private households, each equipped with a rooftop PV installation with different characteristics, as reported in Table I. PV converters operate in MPPT mode at unitary power factor. The total PV installed capacity is 35.3 kWp and the peak demand is 12 kWp. The households are also equipped with grid-connected battery systems with bidirectional power converter to implement PV energy self-consumption policies (actuated at 5 minutes resolution). Batteries specifications are summarized in Table I. Battery injections are monitored. They are removed from the power flow at the PCC by algebraic difference, such that, in the following analysis, it is possible to consider two cases: with and without battery action (self consumption).

PV and power flow measurements

The power flow at the PCC is the sum of the four households flows, measured synchronously at 10 s resolution. Similarly, the global PV production (used as the ground truth value to validate the estimation performance of the algorithms) is the sum of the single PV facilities power injections, measured at the converter level.

Global horizontal irradiance measurements

GHI measurements are from a pyranometer installed on the roof of the household ID1. The line distance of the remaining households from the GHI observation point is shown in the second to last column of Table I. All the measurements are synchronized and timestamped, and logged in a time series database.

V Results and Discussion

For a visual exemplification of the disaggregation process, the reader is referred to Section II-C and Fig. 2. In this section, we first assess the performance of the proposed methods individually. In V-E, we perform a joint performance assessment to compare the quality of the estimations of the different algorithms and support the assertion that they can be considered unsupervised. Key results are discussed and summarized in V-F. In V-G, the algorithms are tested in scenarios with a lower penetration of PV generation to verify if less prominent patterns of PV generation are detrimental to estimation performance. Finally in V-H, we discuss the computational performance.

V-A Method A

The estimation nRMSE as a function of the only parameter of Method A (i.e. sampling time) with and without PV self-consumption is shown in Fig. 3. With no self-consumption, the nRMSE stabilizes at around 2 kW for sampling times larger than 200 s; with self-consumption, performance is poorer and is best at around 150 s. In both cases, estimation performance when the input time series is densely sampled (large sampling time) is poor.

V-B Method B

On top of the series sampling time, Method B has a smoothing parameter $\lambda$ in (5) to weight the demand time variations $\widehat{L}_{k}-\widehat{L}_{k-1}$ in the cost function. Performance as a function of the two parameters is shown in Fig. 6. With self-consumption, best performance happens in the middle right region of the parameter space. This region is larger in the case without self-consumption. Performance degradation patterns do not have a well identifiable trend.

V-C Method C

The parameters of Method C are input time series sampling time and piecewise constant segment length $c$ (in number of samples). Their influence on the nRMSE is shown in Fig. 6. Estimation performance decreases when moving away from the axis origin, denoting that using densely sampled input time series and small $c$ values (the best performance is with 20 s resolution) are convenient. As mentioned in the formulation stage, this is to be expected because the choice of the two parameters affects the constant segment length of the demand piecewise constant model (i.e., the shorter it is, the better performance the persistence model has). Estimation performance is worse with self-consumption (a numerical quantification is given in the next paragraph). Contour lines of Fig. 6 denote that the performance degradation follows the same pattern when with and without self-consumption. Thus, even if estimation performance is different in the two cases, the optimal values of the parameters lay in the same parameters space area.

V-D Method D

Method D parameters are the lower and upper cut-off frequencies of the bandpass filter, and a tuning constant of the bisquare loss function $\rho(\cdot)$. The last was found not to impact substantially on the algorithm performance and is therefore excluded from the current analysis. The sensitivity of algorithm performance to upper and lower cut-off frequencies is shown in Fig. 6. Best performance happens in a well identifiable region in the lower left part of the parameter space, which however tend to shrink in the case with self-consumption.

V-E Joint Performance Comparison

The min, max, mean, and median statistics of the estimation nRMSE, nMAE and nME of the 4 algorithms are reported in Table II. For each algorithm, the reported statistics are calculated over all the combinations of the considered parameters values. They are to be interpreted in the following way:

• •

min: performance to be expected assuming to know a-priori the best performing set of parameters.

• •

max: performance to be expected when choosing the worst possible combination of parameters.

• •

mean: performance to be expected when choosing a random combination in the parameters space.

• •

median: to evaluate performance distribution skewness.

Table II shows that all the methods perform poorer under self-consumption regimes (e.g. mean nRMSE 4.6 to 5.6% and 5.3 to 7% for Method D and Method C, respectively), in other words when the demand includes a component anti-correlated with PV generation. In terms of mean and median statistics, Method D scores the best metrics, followed by C, B and A. In terms of min value, Method C outperforms the other, except for for the cases nMAE and nRMSE without self-consumption, where Method B is better, and nME with self-consumption, where Method D is the absolute best for all the metrics.

V-F Discussion

The previous results showed that the algorithms with the largest number of best scores is Method D, followed by C, B and A. If only sparsely sampled power flow observations are available (such as those from smart meters, typically at 15 minutes resolution), Method D should be selected because it keeps good performance even at low resolutions. If densely sampled observations are available, the performance of Method C and D are comparable. In this case, Method C has two advantages: i) parameters can be selected with an educated physical-based guess, ii) degradation patterns are the same for both with and without self-consumption.

When selecting the two best performing Methods (C and D), the mean nRMSE is 7% and 5.5% and nME is rather small, 0.5 and 0.8%. The latter metric is of importance because it denotes that estimators are almost unbiased, in other words, even if a single estimation in time is wrong, the estimated global PV production over a period is nearly correct.

It is worth noting that, in the proposed sensitivity analysis, PV measurements have been used to assess estimation performance. However, in practical applications when PV observations are not available, it is not possible to do so otherwise there would be no use for disaggregation algorithms. In case of Method C (min/max nRMSE in the range $4.0\div 13.6$%), parameters can be chosen with an educated guess in order to get closer to the best performance. As far as Method D is concerned, the min/max nRMSE gap is smaller ($4.2\div 9.6$%), and estimation performance is good over a wide area of the parameter space. In other words, although it is not possible to derive analytical criteria for performing an a-priori assignment of the parameters, empirical results showed that, in the proposed case, parameters can be chosen in a wide range of values without a sensible deterioration of the performance.

V-G Extension to cases with lower PV generation levels

In this section, the algorithms with the best performing parameters from the previous analysis are tested in scenarios with lower PV production capacity to verify how less prominent PV generation patterns are detrimental to their performance. In each scenario, synthetic time series of the power flow at the PCC are generated as the sum between the demand and a fraction of the original PV generation. Two additional scenarios are considered, with 50% and 25% of the original PV capacity (17.5 and 8.75 kW, respectively), which correspond to 146 and 76% PV penetration levels (i.e. installed PV capacity over the observed peak demand), respectively. Results are reported in Table III. Results show that lower PV penetration levels affects estimation performance of methods B and C (their performance worsens approximately by a factor 2 with a quarter of PV generation), whereas Method D is more robust and its performance is minimally affected.

V-H Computational aspects

In the disaggregation and estimation process, the only required real-time operation is the computation of (1), a cheap task which involves algebraic and trigonometric relationships. Computing $\boldsymbol{\alpha}^{o}$ is a training process without real-time requirements which can be executed off-line with historical data. The computational time is 183 s for Method A, 709 s for B, 103 s for C and 67 s for D.444Computational times refer to a workstation equipped with an Intel Xeon processor running at 2.70 GHz with Matlab on a virtualized operating system. Method C and Method A,B and D were executed on two different machines, machine 1 and 2. The computation time of Method C was adjusted by a factor $t_{2}/t_{1}$, where $t_{1}$ and $t_{2}$ is the computation time of a reference problem executed on machine 1 and 2, respectively.

VI Conclusions

The problem of disaggregating a sequence of active power flow measurements composed of unobserved PV generation and demand into the respective trajectories was considered. Four disaggregation algorithms were discussed. They attempt to explain similarities between the time series of the aggregated and estimated PV generation, three in the time domain and one in the frequency domain. Estimation algorithms leverage GHI measurements transposed onto a number of tilted planes with the objective of explaining PV production patterns from sites with potentially different configurations (a key feature in urban/suburban context where PV generation is mostly from rooftop PV facilities with tilt/azimuth configurations dictated by roof characteristics). The effect of the air temperature was modelled by preprocessing GHI values with a model-based approach. Algorithms require an offline optimization problem-based training phase with historical data. For three algorithms, the convexity of the underlying optimization problem, important to assure tractability, is verifiable a-priori by inspecting the input data. Reconstructing the PV power output requires computing on-line an algebraic relationship and is suitable for implementation with deterministic deadlines and low processing power, in real-time target devices. Algorithm performance was tested with data from a real-life setup, with PV generation from multiple sites with different configurations, different demand profiles, and battery systems for PV self-consumption. Results show that the best performing algorithms estimate PV generation with a root mean square and mean estimation errors in the ranges $3.4\div 8.8$% and $0.5\div 2.3$%, respectively, and that performance is minimally affected by the level of PV penetration in the prosumption mix. The practical utility of the proposed algorithms is envisaged in the context of power and energy management of distributed energy resources and data-driven PV generation forecasting in those situations where information from PV plants is not available due to issues such as privacy concerns or lack of adequate communication infrastructures.

Appendix A On the convexity of Method B and C

We discuss on the convexity of the problem (9)-(11). Let $\boldsymbol{P}=[P_{1},P_{2},\cdots,P_{K}]$, $\boldsymbol{I}^{\diagdown}_{j}=[I^{\diagdown}_{j,1},\cdots,I^{\diagdown}_{jK}],\ j=1,\dots,J$, $M\in\mathbf{R}^{K\times J}=\begin{pmatrix}\boldsymbol{I}^{\diagdown}_{1},\boldsymbol{I}^{\diagdown}_{2},\ldots,&\boldsymbol{I}^{\diagdown}_{J}\end{pmatrix}$ the matrix obtained by stacking horizontally the GNI columns. The estimated total PV production (1) is (matrix product) $\widehat{\boldsymbol{G}}=M\boldsymbol{\alpha}$, which replaced into (4) yields:

[TABLE]

where $\mathbb{1}$ is the $K\times K$ identity matrix, $S=\begin{pmatrix}\mathbb{1}_{K\times K},-M\end{pmatrix}\in\mathbf{R}^{K\times(K+J)}$ and $x=(\widehat{\boldsymbol{L}},\boldsymbol{\alpha})^{T}$. The least square cost (9) is:

[TABLE]

Minimizing (17) is the same as minimizing (minimization is invariant under sum with constants and scale factors):

[TABLE]

where (16) is used, $H=S^{T}S$ and $\boldsymbol{f}=S^{T}\boldsymbol{P}$. Eq. 19 is convex if $H$ is semidefinite positive. Since $H$ depends on input data, convexity cannot be enforced by construction, but it can be checked numerically. It was noted that adding a regularization term to the matrix quadratic matrix $H^{\prime}=H+\beta\cdot\mathbb{1}$ ($\beta=$1\text{\times}{10}^{-4}\text{,}$$) helps to achieve convexity while not impacting substantially on algorithms performance.