Estimation of Photovoltaic Generation Forecasting Models using Limited Information

TL;DR

This paper introduces a low-complexity algorithm for estimating photovoltaic generation models using limited data, specifically only power measurements, a clear-sky irradiance model, and temperature forecasts, validated on real data.

Contribution

A novel estimation method for PV models that operates without on-site meteorological measurements, relying solely on power data and forecasted variables.

Findings

Effective estimation with limited data

Validated on real photovoltaic plant data

Low computational complexity

Abstract

This work deals with the problem of estimating a photovoltaic generation forecasting model in scenarios where measurements of meteorological variables (i.e. solar irradiance and temperature) at the plant site are not available. A novel algorithm for the estimation of the parameters of the well-known PVUSA model of a photovoltaic plant is proposed. Such a method is characterized by a low computational complexity, and efficiently exploits only power generation measurements, a theoretical clear-sky irradiance model, and temperature forecasts provided by a meteorological service. An extensive experimental validation of the proposed method on real data is also presented.

Figures

Tables

Table 1. Table 1: Quality indices of irradiance and air temperature forecasts for data set D2.

Table 2. Table 2: Initial parameters and panel orientation angles.

Table 3. Table 3: Algorithm parameters for validation on measured data (D1)

Table 4. Table 4: Performance comparison of CSD, SRLS and ODNP computed starting from day 28 28 28 (D1).

Performance Indices		CSD	SRLS	ODNP
DA Forecast	$RMSE$ ($\mathrm{kW}$)	$31.0$	$23.1$	$143.2$
	$MAPE$	$31\%$	$26\%$	$109\%$
	$MBE$ ($\mathrm{kW}$)	$-7.01$	$-6.73$	$3.00$
	$R^2$	$0.98$	$0.99$	$0.65$
	$NRMSE$	$0.13$	$0.10$	$0.59$
	$RMS{E}_{NP}$	$0.032$	$0.024$	$0.15$
	$MAP{E}_{NP}$	$2.2\%$	$1.5\%$	$8.4\%$

Table 5. Table 5: Data sets used in the evaluation of the effect of β 0 subscript 𝛽 0 \beta_{0} . POD denotes the percentage of scaled clear-sky data introduced.

Data set ID	D1	D11	D12
POD	$0\%$	$5\%$	$14\%$

Table 6. Table 6: Mean values and standard deviations of the parameters estimated using D 1 2 subscript 1 2 1_{2} and three different values of β 0 subscript 𝛽 0 \beta_{0} (computed from day 28 28 28 ).

Table 7. Table 7: Algorithm parameters

Table 8. Table 8: Performance comparison of CSD, SRLS and ODNP computed starting from day 28 28 28 (D2).

Performance Indices		CSD	SRLS	ODNP
DA Forecast	$RMSE$ ($\mathrm{kW}$)	$117.9$	$118.5$	$193.3$
	$MAPE$	$58.8\%$	$55.2\%$	$85.6\%$
	$MBE$ ($\mathrm{kW}$)	$-7.69$	$35.6$	$-5.6$
	$R^2$	$0.799$	$0.797$	$0.458$
	$NRMSE$	$0.448$	$0.451$	$0.736$
	$RMS{E}_{NP}$	$0.128$	$0.129$	$0.201$
	$MAP{E}_{NP}$	$8.3\%$	$9.8\%$	$12.4\%$
HA Forecast	$RMSE$ ($\mathrm{kW}$)	$138.2$	$136.2$	-
	$MAPE$	$52.1\%$	$46.0\%$	-
	$MBE$ ($\mathrm{kW}$)	$-25.8$	$33.0$	-
	$R^2$	$0.655$	$0.665$	-
	$NRMSE$	$0.588$	$0.579$	-
	$RMS{E}_{NP}$	$0.150$	$0.148$	-
	$MAP{E}_{NP}$	$10.0\%$	$11.9\%$	-

Equations

$$P={\mu }_{1}I+{\mu }_2{I}^2+{\mu }_{3}IT,$$

$$P={\mu }_1\cdot \alpha (I,T)\cdot I,$$

$$\alpha (I,T)=1+{\eta }_{2}I+{\eta }_{3}T,$$

$${\eta }_2={\mu }_2/{\mu }_1,{\eta }_3={\mu }_3/{\mu }_1.$$

$$\begin{split}\eta_{2}&\in\left[\underline{\eta}_{2},\overline{\eta}_{2}\right]=\left[$-2.5\text{\times}{10}^{-4}$,$-1.9\text{\times}{10}^{-5}$\right],\\ \eta_{3}&\in\left[\underline{\eta}_{3},\overline{\eta}_{3}\right]=\left[$-4.8\text{\times}{10}^{-3}$,$-1.7\text{\times}{10}^{-3}$\right].\end{split}$$

$$P={\varphi }^{\prime }(I,T)\text{leavevmode}\mu ,$$

$$\varphi (I,T)=[I\text{leavevmode}\text{leavevmode}{I}^2\text{leavevmode}\text{leavevmode}IT{]}^{\prime },$$

$$\mu =[{\mu }_1\text{leavevmode}\text{leavevmode}{\mu }_2\text{leavevmode}\text{leavevmode}{\mu }_3{]}^{\prime }.$$

$$\hat{P}={\varphi }^{\prime }(\hat{I},\hat{T})\cdot \hat{\mu }.$$

$${\hat{P}}^{cs}={\varphi }^{\prime }(I^{cs},\hat{T})\cdot \hat{\mu }.$$

$$I^{cs,n}=\{\begin{array}{ll}A\cdot 0.7^{{\left(\frac{1}{\sin h}\right)}^{0.678}}& \text{if }0<h<\pi /2\\ 0& \text{otherwise,}\end{array}$$

$$I^{cs}=\left[\sin (\psi )\cos (h)\cos (\zeta -\gamma )+\cos (\psi )\sin (h)\right]I^{cs,n},$$

$$\underline{\alpha }(I,T)\le \alpha (I,T)\le \overline{\alpha }(I,T),$$

$$\underline{\alpha }(I,T)$$

$$\overline{\alpha }(I,T)$$

$${\mu }_1\cdot I\cdot \underline{\alpha }(I,T)\le P\le {\mu }_1\cdot I\cdot \overline{\alpha }(I,T).$$

$$\begin{array}{rl}\Delta P(j)& =P(j)-P(j-1)\\ & ={\mu }_1\left[I(j-1)\Delta \alpha (j)+\Delta I(j)\alpha (I(j),T(j))\right],\end{array}$$

$$\Delta I(j)=I(j)-I(j-1),$$

$$\Delta \alpha (j)=\alpha (I(j),T(j))-\alpha (I(j-1),T(j-1)).$$

$$\underline{\Delta \alpha }(j)\le \Delta \alpha (j)\le \overline{\Delta \alpha }(j),$$

$$\left[\begin{array}[]{c}\underline{\Delta\alpha}(j)\\ \overline{\Delta\alpha}(j)\end{array}\right]=Q(j)\left[\begin{array}[]{c}\Delta I(j)\\ \Delta T(j)\end{array}\right]$$

$$\begin{array}[]{c|c|c}\hline\cr\hline\cr Q(j)&\Delta T(j)\geq 0&\Delta T(j)<0\\ \hline\cr\Delta I(j)\geq 0&\begin{bmatrix}\underline{\eta}_{2}&\underline{\eta}_{3}\\ \overline{\eta}_{2}&\overline{\eta}_{3}\end{bmatrix}&\begin{bmatrix}\underline{\eta}_{2}&\overline{\eta}_{3}\\ \overline{\eta}_{2}&\underline{\eta}_{3}\end{bmatrix}\\ \hline\cr\Delta I(j)<0&\begin{bmatrix}\overline{\eta}_{2}&\overline{\eta}_{3}\\ \underline{\eta}_{2}&\underline{\eta}_{3}\end{bmatrix}&\begin{bmatrix}\overline{\eta}_{2}&\overline{\eta}_{3}\\ \underline{\eta}_{2}&\overline{\eta}_{3}\end{bmatrix}\\ \hline\cr\hline\cr\end{array}$$

$${\mu }_1{\underline{\delta }}_P(j)$$

$$\left[\begin{array}[]{c}\underline{\delta}_{P}(j)\\ \overline{\delta}_{P}(j)\end{array}\right]=R(j)\left[\begin{array}[]{c}I(j-1)\\ \Delta I(j)\end{array}\right]$$

$$\begin{array}[]{c|c}\hline\cr\hline\cr R(j)&\nobreak\leavevmode\hfil\\ \hline\cr\Delta I(j)\geq 0&\begin{bmatrix}\underline{\Delta\alpha}(j)&\underline{\alpha}\left(I(j),T(j)\right)\\ \overline{\Delta\alpha}(j)&\overline{\alpha}\left(I(j),T(j)\right)\end{bmatrix}\\ \hline\cr\Delta I(j)<0&\begin{bmatrix}\underline{\Delta\alpha}(j)&\overline{\alpha}\left(I(j),T(j)\right)\\ \overline{\Delta\alpha}(j)&\underline{\alpha}\left(I(j),T(j)\right)\end{bmatrix}\\ \hline\cr\hline\cr\end{array}$$

$${\mathcal{P}}_{\mathcal{J}}=\left\{\{P^m(j),T(j),P^{cs}(j)\right\},j\in \mathcal{J}\},$$

$$P^{cs}(j)$$

$$\begin{array}{r}{\mu }_1\cdot I^{cs}(j)\cdot \underline{\alpha }\left(I^{cs}(j),T(j)\right)\le P^{cs}(j)\\ \le {\mu }_1\cdot I^{cs}(j)\cdot \overline{\alpha }\left(I^{cs}(j),T(j)\right).\end{array}$$

$$j_{max}$$

$$I_{max}^{cs}$$

$$P_{max}^{cs}$$

$$\begin{array}{r}\frac{I^{cs}(j)\cdot \underline{\alpha }\left(I^{cs}(j),T(j)\right)}{I_{max}^{cs}\cdot \alpha (I_{max}^{cs},T(j_{max}))}\le \frac{P^{cs}(j)}{P_{max}^{cs}}\\ \le \frac{I^{cs}(j)\cdot \overline{\alpha }\left(I^{cs}(j),T(j)\right)}{I_{max}^{cs}\cdot \alpha \left(I_{max}^{cs},T(j_{max})\right)},\end{array}$$

Full text

††journal: Automatica

1 Introduction

The electrical grid can be no longer considered a unidirectional means of distributing energy from conventional plants to the final users, but a Smart Grid, where strong interaction between producers and users takes place [1]. A major challenge in the integration of renewable energy sources into the grid [2] is that power generation is intermittent, difficult to control, and strongly dependent on the variation of weather conditions. For this reason, forecasting of renewable distributed generation has become a fundamental requirement in order to reliably manage conventional power plant operation, grid balancing, real-time unit dispatching [3], demand constraints [4], and energy market requirements. In this respect, renewable generation forecasts on different time horizons are of special interest to various players that operate in the active grid, in particular to Distribution System Operators (DSO) and Transmission System Operators (TSO) (see [1, 5, 6] and references therein).

Concerning photovoltaic (PV) power generation, researchers have devoted much attention to the problem of obtaining accurate generation forecasts over different time horizons, e.g., day-ahead and hour-ahead [7, 8]. Most contributions, however, focus on the problem of solar irradiance prediction [9, 10, 11, 12]. To tackle this problem, several approaches based on Artificial Neural Networks (ANNs) [13, 14, 15] or Support Vector Machines [16] can be found in the literature. Alternatively, classical linear time series prediction methods are used in [17, 18, 19], where the considered time series is typically the global horizontal irradiance (GHI) [20, 21]. GHI forecasts are typically used along with temperature forecasts in a simulation model of the PV plant [22] in order to calculate generated power predictions. In all cases, computing reliable forecasts from predicted meteorological variables hinges upon the availability of an accurate model of the plant, be it physical or estimated from data.

Unfortunately, in many common scenarios, neither a plant model, nor direct on-site measurements of solar irradiance and other meteorological variables (e.g., temperature) are available. This is always the case with a DSO dealing with hundreds or thousands of heterogeneous, independently owned and operated PV plants; in this case, the only available data consists of generated power measurements provided by electronic meters, and of irradiance and temperature forecasts provided by a meteorological service. The problem of forecasting power generation in this case is addressed in [23] by means of a neural network and in [24, 25] using a parametric model. In these approaches, however, further information on the cloud cover index at the plant site is assumed to be available. In [26, 27], a heuristic method for the estimation of the parameters of well-known PVUSA model [28] based on theoretical clear-sky irradiance is presented, while in [29], a recursive procedure based on the clear-sky criteria proposed in [30] is devised. However, the former approach does not allow for capturing possible parameter variations or seasonal drifts, and moreover both approaches require trial-and-error in order to manually tune a number of algorithm parameters whose values may vary significantly according to the climate zone.

1.1 Paper contribution

In this paper, a novel approach to the problem of estimating the parameters of the PVUSA model in the partial information case is presented. The only historical data used by the method consist of generated power, and temperature (but not irradiance) forecasts. Parameter estimates are obtained by carefully exploiting the information contained in portions of the generated power data which turn out to be meaningful if combined with theoretical clear-sky irradiance over the same period. More specifically, we introduce three tests to be performed on generated power data in order to detect portions of such data that were generated under clear-sky conditions. The information contained in such portions is then exploited in a recursive parameter estimation algorithm in combination with theoretical clear-sky irradiance provided by a suitable model. The method proposed in this paper improves over [26, 27, 29], since it is able to adapt to parameter variations and requires the tuning of a single threshold coefficient whose physical role can be interpreted in terms of the cloud cover factor (CCF) [31].

The paper is structured as follows: in Section 2 the modeling tools are introduced; in Section 3 the proposed clear-sky detection tests are developed; the model estimation procedure is presented in Section 4. In Section 5, the relevant forecasting problems are recalled, and performance evaluation criteria are discussed in Section 6. Experimental validation results are reported in Section 7, and conclusions are drawn in Section 8.

2 Models and methods

2.1 The PVUSA photovoltaic plant model

A PV plant can be efficiently modelled using the PVUSA model [32], which expresses the instantaneous generated power as a function of irradiance and air temperature according to the equation:

[TABLE]

where $P$, $I$, and $T$ are the generated power ($\mathrm{kW}$), irradiance ($\mathrm{W}\mathrm{/}\mathrm{m}^{2}$), and air temperature ($\mathrm{\SIUnitSymbolCelsius}$), respectively, and $\mu_{1}$, $\mu_{2}$, $\mu_{3}$, are the model parameters. It is important to notice that model (1) is linear in the parameters. For the purpose of this work, it is useful to express (1) in the form

[TABLE]

where

[TABLE]

being

[TABLE]

From (2), it is apparent that $\mu_{1}$ represents the main power/irradiance gain of the plant, while $\alpha(I,T)$ in (3) can be seen as correction term. In this respect, it is worth noticing that the ratios $\eta_{2}$ and $\eta_{3}$ in (4) are characterized by well-established variability ranges among different PV technologies (see [32]). Such ranges are given by:

[TABLE]

This property will be exploited in the proposed estimation procedure. It is also appropriate to represent (1) also in the standard regressive form

[TABLE]

where the regressor is given by

[TABLE]

and the parameter vector is

[TABLE]

The PVUSA model can be fruitfully exploited for the purpose of computing forecasts of generated power on the basis of predicted meteorological variables. Indeed, once a correct estimate $\hat{\mu}$ of the parameter vector is available, a reliable power generation forecast $\hat{P}$ can be obtained by substituting predicted irradiance $\hat{I}$ and temperature $\hat{T}$, provided by a meteorological service, into the model equation (6), i.e.,

[TABLE]

Similarly, generation forecasts under clear-sky conditions can be obtained by using the theoretical irradiance $I^{cs}$ at the plant location, as provided by a suitable model, and temperature forecasts, i.e.,

[TABLE]

Notice that clear-sky generation forecasts provide an upper bound on the power that can be generated by a plant, and as such they can be used by the DSO, for instance, when scheduling the maintenance of the portion of the grid where the plant is located. Despite its simplicity, very good forecasting accuracy is obtained from the PVUSA model when the parameter vector $\mu$ is estimated using measured irradiance and temperature data via, e.g., standard least squares fitting (see, e.g., [26]).

A DSO that manages a high number of independent generation facilities may not have access to time series of irradiance and temperature measured on the premises of each plant, while power generation data are always available through meters. In order to estimate model parameters, replacing the measured values of $I$ and $T$ with forecasts $\hat{I}$ and $\hat{T}$ provided by a meteorological service is not a viable solution, due to the fact that forecasting errors on the irradiance are in general too large. On the contrary, temperature forecasts are quite reliable and can be used in place of actual measurements (see [26, 27] for details).

2.2 Clear-sky irradiance model

In this paper, a theoretical estimate of the global clear-sky irradiance on a given surface is required. To this aim, although several different models are present in the literature [33], the Heliodon simulator model [34] is used. This model is characterized by a high degree of simplicity and allows to compute the theoretical global clear-sky normal irradiance ($\mathrm{W}\mathrm{/}\mathrm{m}^{2}$) from the solar altitude $h$, i.e., the angle over the horizon (rads), as:

[TABLE]

where $A=1353$ $\mathrm{W}\mathrm{/}\mathrm{m}^{2}$ denotes the apparent extraterrestrial irradiance. Given the theoretical clear-sky normal irradiance $I^{cs,n}$, the clear-sky irradiance on an inclined panel surface $I^{cs}$ can be derived from $I^{cs,n}$ and the orientation of the surface with respect to the sun position. Denoting by $\zeta$ the surface azimuth and $\psi$ the surface tilt angle, one has that

[TABLE]

where $\gamma$ is the solar azimuth. Clearly, $I^{cs}$ can be computed for given values of $\zeta$ and $\psi$ from latitude, longitude and time of day. For $\psi=0$, the irradiance on a horizontal surface is obtained.

In this study, it is assumed that the exact orientation of the PV panel surfaces of the considered plant is not known a-priori. However, it is reasonable to suppose that the plant is efficiently oriented for the specific latitude according to, e.g., the guidelines given in [35]. Therefore, the value of (12) with $(\zeta,\psi)$ taken from the above guidelines will be used as a reference for the theoretical clear-sky irradiance $I^{cs}$ for a given plant.

3 Clear-sky data detection

In this paper, the following key idea is exploited for the purpose of estimating the parameters of the PVUSA model (1) of a PV plant without resorting to on-site irradiance measurements. Given a time series composed of generated power measurements and temperature forecasts (or measurements, if available), suitable tests can be performed on the data in order to detect portions of the power curve which have been generated under a clear-sky condition; this allows for fitting the parameters of the PVUSA model to such data by using theoretical clear-sky irradiance (e.g., via the model (11),(12)) in the regressor of (7) in place of the actual measured irradiance. This section deals with the derivation of such tests, which will be referred to as CS tests in the sequel.

In view of (5), suitable bounds can be derived on $\alpha(I,T)$ and $P$ in the PVUSA model (2)-(3). Indeed, from (3) and (5), it is easily checked that

[TABLE]

where

[TABLE]

Moreover, for realistic values of $I$ and $T$, it always holds that $\underline{\alpha}(I,T)>0$ and $\overline{\alpha}(I,T)<1$. From (13) and (2), the following bound on $P$ is obtained:

[TABLE]

Let us now consider a time series $\{P(j)$, $I(j)$, $T(j)\}$ of the variables in (1), where $j$ represents a discrete time index. The increment of $P(j)$ can be expressed as

[TABLE]

where

[TABLE]

Let $\Delta T(j)=T(j)-T(j-1)$. Taking into account (14)-(15), it is easily checked that the following bounds on $\Delta\alpha(j)$ hold:

[TABLE]

where

[TABLE]

and the matrix $Q(j)$ depends on the signs of $\Delta I(j)$ and $\Delta T(j)$ according to the following table:

[TABLE]

In view of (17), this allows to derive the following bounds on $\Delta P(j)$:

[TABLE]

where

[TABLE]

and the matrix $R(j)$, depending on the sign of $\Delta I(j)$, is given by

[TABLE]

The bounds (16) and (20) allow to devise the sought CS tests. Let us consider a time interval $\mathcal{J}$, and the following associated time series

[TABLE]

where, for each time instant $j$, $P^{m}(j)$ represents the measured plant power reported by meters, $T(j)$ is a temperature forecast (or measurement), and $P^{cs}(j)$ is the clear-sky generated power predicted by a PVUSA model characterized by given values of the parameters $\mu_{1},\mu_{2}$, and $\mu_{3}$, i.e.,

[TABLE]

where the clear-sky irradiance $I^{cs}(j)$ is computed, e.g., via (11). Clearly, by (16),

[TABLE]

Now let

[TABLE]

The quantities $I^{cs}_{max}$, $P^{cs}_{max}$, and $j_{max}$ define, respectively, the maximum clear-sky irradiance, the maximum predicted clear-sky generated power, and the time index for which this maximum value occurs within the given time window $\mathcal{J}$.

Normalizing (24) with respect to $P^{cs}_{max}$ yields

[TABLE]

and hence the following bounds on the ratio $\frac{P^{cs}(j)}{P^{cs}_{max}}$ hold:

[TABLE]

where

[TABLE]

It is important to observe that (29)-(30) define bounds on the clear-sky power time series which do not depend on the model parameters. Condition (29) can be exploited in order to classify a time window $\mathcal{J}$ of measured power data points as generated under clear-sky conditions. Indeed, given the time series $\{P^{m}(j)$,$T(j)$, $j\in{\mathcal{J}}\}$, the following test is introduced:

CS Test 1**.**

[TABLE]

The satisfaction of CS test 1 is in general not sufficient to classify power data within ${\mathcal{J}}$ as having been generated under a clear-sky condition. Specifically, if the sky is partially cloudy during the time interval ${\mathcal{J}}$, the measured power may heavily oscillate, but could remain quite close to the clear-sky power at the maximum [30], thus satisfying (31). To overcome this issue, a further condition on the normalized increment of the power time series is derived. Let $\underline{\delta}_{P}^{cs}(j)$ and $\overline{\delta}_{P}^{cs}(j)$ be defined by (21) evaluated for $I(j)=I^{cs}(j)$ and $\Delta I(j)=\Delta I^{cs}(j)=I^{cs}(j)-I^{cs}(j-1)$. The increment of $P^{cs}(j)$ is given by

[TABLE]

and satisfies

[TABLE]

by (20). Normalizing  (32) with respect to $P^{cs}_{max}$, the following bounds on the normalized increment $\frac{\Delta P^{cs}(j)}{P^{cs}_{max}}$ are obtained:

[TABLE]

i.e.,

[TABLE]

where

[TABLE]

Note that the bounds (34)-(35), similarly to (29)-(30), do not depend on the model parameters. Condition (34) provides the second criterion for classifying a time window ${\mathcal{J}}$ of measured power data points as clear-sky. The following test is introduced:

CS Test 2**.**

[TABLE]

where $\Delta P^{m}(j)$ is the increment of the measured power, i.e., $\Delta P^{m}(j)=P^{m}(j)-P^{m}(j-1)$.

CS tests 1 and 2 detect deviations in the shape of the normalized power curve from the clear-sky condition caused by cloudiness in different scenarios. However, due to normalization, such conditions may turn out to be fulfilled on a given time window ${\mathcal{J}}$ when the corresponding data are generated under perfectly uniform cloudiness, i.e., when the actual irradiance satisfies

[TABLE]

where $0<\beta<1$ is a constant that represents a uniform cloud cover factor (see [36]) in the time window $\mathcal{J}$. If the data collected within such a time window are used to perform a model parameter adaptation step in a recursive estimation procedure, the algorithm may tend to underestimate the power/irradiance gain of the plant at such step. This fact may be detrimental when a long series of data collected under uniform cloudiness is processed. To mitigate this effect, a further test is introduced. Suppose that a current estimate $\hat{\mu}$ of the model parameters is available. Accordingly, a current estimate of the generated power under clear-sky conditions is given by

[TABLE]

Let $\hat{P}^{cs}_{max}$ be the peak value of $\hat{P}^{cs}(j)$ in $\mathcal{J}$, i.e.,

[TABLE]

Provided that CS tests 1 and 2 are passed by the data in time window $\mathcal{J}$, the following further condition is introduced, which involves a comparison of the maximum currently predicted clear-sky power $\hat{P}^{cs}_{max}$ with the corresponding generated power $P^{m}(j_{max})$ as follows:

CS Test 3**.**

[TABLE]

where $0<\epsilon<1$ is a parameter chosen by the designer, typically a number slightly higher than 0, whose role and choice is discussed later on. CS test 3 has the specific role of detecting, under the condition that CS tests 1 and 2 are satisfied, whether the peak value of measured power within the considered time window lies above a given fraction of the clear-sky power currently estimated by the model. Condition (38) can be satisfied in the following cases:

$\frac{P^{m}(j_{max})}{\hat{P}^{cs}_{max}}\geq 1$: in this case the model is currently underestimating clear-sky power;

• 2.

$1-\epsilon\leq\frac{P^{m}(j_{max})}{\hat{P}^{cs}_{max}}<1$ and the current model is overestimating the generated power by a small amount;

• 3.

$1-\epsilon\leq\frac{P^{m}(j_{max})}{\hat{P}^{cs}_{max}}<1$ and uniform cloudiness is present within the given time window, so that generation is marginally lower than the clear-sky power currently predicted by the model.

With the exception of the latter case, the simultaneous satisfaction of CS tests 1,2, and 3 requires that the model parameters be adapted in order to fit the measured power series with the predicted one within $\mathcal{J}$.

Remark 1**.**

The parameter $\epsilon$ plays a key role in detecting whether the clear-sky curve provided by the model matches or underestimates power data satisfying CS tests 1 and 2, which are related to the shape of the normalized power curve. Setting this value very close to zero allows for good adaptation when the model is underestimating the clear-sky power (for this reason it is advisable to choose an underestimate of $\mu_{1}$ as the initial guess in the estimation procedure, as detailed in the next section). Higher values, on the other hand, allow for adjusting the model when it is overestimating; the latter case is very important for capturing possible slow parameter drifts as well as seasonal variations in the accuracy of the theoretical clear-sky model. However, increasing $\epsilon$ may cause adaptation to long series of data generated under uniform cloudiness. To further clarify this aspect, let us assume that the true plant is described by a PVUSA model (2) characterized by $\mu_{1}=\mu_{1}^{0}$ and $\alpha(I,T)=\alpha^{0}(I,T)$, and that uniform cloudiness is present within $\mathcal{J}$ so that (37) holds for some $0<\beta<1$. It follows that

[TABLE]

Therefore condition (38) becomes:

[TABLE]

For given $\epsilon$, a rough estimate of the values of the uniform cloud cover factor $\beta$ for which CS test 3 is satisfied can be obtained by observing that $\frac{\alpha^{0}(I^{cs}_{max}\beta,T(j_{max}))}{\hat{\alpha}(I^{cs}_{max},T(j_{max}))}\approx 1$ (especially for $\beta$ close to 1) and that a rough approximation of the main power/irradiance gain $\mu_{1}^{0}$ is given by the ratio $P_{nom}/1000$, where $P_{nom}$ denotes the nominal plant power in kW. In view of (39), CS test 3 is passed when $\beta$ approximately satisfies

[TABLE]

The relationship (40) provides an interpretation of the parameter $\epsilon$ and represents a possible guideline for tuning such parameter on the basis of the minimum value of the cloud cover factor for which the designer allows data generated under uniform cloudiness to be considered for parameter adaptation. However, the effect of the choice of $\epsilon$ will be extensively discussed on the basis of real data in the experimental section of this paper. **

4 Model estimation

According to the observations in the previous sections, we now introduce the proposed PVUSA plant model estimation method, which yields an on-line update of the parameter vector estimate $\hat{\mu}$ by relying only on the information contained on a time series composed by measured power $P^{m}$ and forecast (or measured) temperature $T$. The model estimation procedure is recursive, and combines CS tests $1-3$ with a standard Recursive Least-Squares (RLS) algorithm using a dynamical time window.

The following definitions are instrumental for building up the procedure:

$k$: present time index;

• 2.

$d$: present day;

• 3.

${\mathcal{I}}_{d}=[\underline{k}_{d},\overline{k}_{d}]$: time interval corresponding to light hours in day $d$, i.e., $I^{cs}(k)>0$ for all $\underline{k}_{d}\leq k\leq\overline{k}_{d}$;

• 4.

${\mathcal{J}}_{k,l}$: set of time indices corresponding to a time window of given length $l$ ending at $k$, i.e., ${\mathcal{J}}_{k,l}=\{k-l+1,\dots,k\}$;

• 5.

$\hat{\mu}(k)$: estimate of the parameter vector at time $k$, being $\hat{\mu}(0)$ the initial guess;

• 6.

$I^{cs}(j)$: theoretical clear-sky solar irradiance at time step $j$, computed according to (11),(12), or a different model for the plant site;

• 7.

${T}(j)$: temperature forecast (or measurement, if available) at time $j$ at the plant site, provided by a meteorological service;

• 8.

$P^{m}(j)$: measured generated power at time $j$;

• 9.

$D(j)=\{P^{m}(j),T(j),I^{cs}(j)\}$: data sample at time $j$;

• 10.

${\mathcal{D}}({\mathcal{J}})=\left\{D(j),\leavevmode\nobreak\ j\in{\mathcal{J}}\right\}$: data set pertaining to time window ${\mathcal{J}}$;

• 11.

$\epsilon$: fixed threshold value $(0<\epsilon<1$);

• 12.

$l_{min}$: minimum time window length.

The estimation algorithm is constructed as follows (see Figure 1). The procedure is reset on each day $d$ at time $k=\underline{k}_{d}$. The current parameter estimate $\hat{\mu}(\underline{k}_{d})$ is initialized with the last estimate obtained on day $d-1$. An initial data set ${\mathcal{D}}({\mathcal{J}}_{k,l_{min}})$ is constructed at time $k=\underline{k}_{d}+l_{min}-1$ corresponding to an initial time window ${\mathcal{J}}_{k,l_{min}}$ of length $l_{min}$. If ${\mathcal{D}}({\mathcal{J}}_{k,l_{min}})$ does not pass CS tests $1-3$, then the procedure is reset at time $k=\underline{k}_{d}+1$. Otherwise (i.e., if ${\mathcal{D}}({\mathcal{J}}_{k,l_{min}})$ is recognized as generated under clear-sky), a new data sample $D(k)$ is acquired at each following step $k$ and added to the current data set ${\mathcal{D}}({\mathcal{J}}_{k,l})$, incrementing the length of the time window ${\mathcal{J}}_{k,l}$ by one. Then, CS tests $1-3$ are performed on ${\mathcal{D}}({\mathcal{J}}_{k,l})$. If tests are passed, then further data samples are added to the data set until one of the tests fails (or the end of the day is reached) at some time $k^{\prime}$. When this occurs, the data set ${\mathcal{D}}({\mathcal{J}}_{k^{\prime}-1,l-1})$ is deemed to be generated under clear-sky conditions and an RLS adaptation step is performed using such data in order to obtain an updated parameter estimate $\hat{\mu}(k^{\prime})$. Then, the algorithm is reset at time $k=k^{\prime}$ and repeated. A detailed description of the procedure is reported in Algorithm 1.

Concerning the selection of the initial parameter guess $\hat{\mu}(0)$, the following observations are in order.

As previously stated, a good guess for the main power/irradiance gain $\mu_{1}$ is represented by $\hat{\mu}_{1}(0)=P_{nom}/1000$, where $P_{nom}$ denotes the nominal plant power [24, 27]. As pointed out in Remark 1, it is appropriate to start with an underestimate of this value, e.g., $75\%$, to ensure faster parameter adaptation.

• 2.

As for the initial values $\hat{\mu}_{2}(0)$ and $\hat{\mu}_{3}(0)$, it is convenient to choose them so that $\mu_{2}(0)/\mu_{1}(0)$ and $\mu_{3}(0)/\mu_{1}(0)$ are equal to the central values of the intervals ${\mathcal{S}}_{2}$ and ${\mathcal{S}}_{3}$ in (5), respectively [24].

5 Forecasting

In this section, we briefly describe how the estimated PVUSA model can be used in order to provide the generation forecasts used in the experimental part of this work. Let $k$ be a generic time instant in which a forecast is supposed to be computed and submitted, e.g., to the DSO. For a given time instant $j\geq k$, let $\hat{W}(j|k)=\{\hat{I}(j|k),\hat{T}(j|k)\}$ denote the weather forecast (irradiance and temperature) relative to time $j$ available at time $k$, where the irradiance forecast is projected on the panel surface using a-priori information on the plant orientation, if available, or a guess thereof taken from guidelines such as those in [35]. The prediction of generated power for time instant $j$, computed at time $k$ using the parameter vector estimate $\hat{\mu}(q)$ available at time $q\leq k$, is given by

[TABLE]

In the following section, the forecasting performance of the PVUSA model estimated using the procedure detailed in Section 4 will be evaluated on the widely used Day-Ahead (DA) and Hour-Ahead (HA) forecasts [37]. The DA forecast is usually submitted at 6 am on the day before each operating day, which begins at midnight on the day of submission, and covers all 24 hours of that operating day. The HA forecast is usually submitted 105 minutes prior to each operating hour and provides an advisory forecast for the 7 hours of light (or the remaining ones, if less) of the same day after the operating hour. The time series representing the DA and HA forecasts can be constructed from the pointwise forecast (41) by letting $j$ and $q$ vary in suitable sets. The details are omitted here for the sake of brevity and the reader is referred to Section 5 of [24].

6 Performance evaluation

In this section we introduce the performance assessment indices that will be used to evaluate the efficacy of the proposed method in the forecasting problems sketched in the previous section.

6.1 Error measures

For the sake of simplicity, a generic definition of the performance indices that will be used is given here. Details on how such indices are computed using a predictor such as (41) in the specific contexts of DA or HA forecasting are provided in [24]. Let $\hat{P}(j)$ represent the forecasted power and $P^{m}(j)$ the corresponding measured value. The following standard error measures are considered:

[TABLE]

where ${\mathcal{K}}=\{1,\dots,K\}$ denotes the time span of the data set and $\bar{P}$ is the sample mean of the measured power. The last two indices, i.e., $RMSE_{NP}$ and $MAPE_{NP}$, are normalized with respect to the nominal plant power $P_{nom}$ and are of practical interest for network operation. In particular, values lower than $10\%$ are considered acceptable for network operation [38, 8].

6.2 Benchmarks

As an additional evaluation tool, the performance indices achieved using the proposed approach will be compared to those obtained using:

ODNP: the One-Day-ahead Naive Predictor, i.e.,

[TABLE]

where $P^{m}_{d-1}(j)$ denotes the measure of generated power recorded during the day before at the same time of day,

• 2.

SRLS: a PVUSA model estimated via a standard RLS algorithm in the complete information case, i.e., using actual measurements of generated power, irradiance and temperature.

7 Experimental results

In this section an extensive validation of the proposed procedure is presented. Two experiments have been run to evaluate the performance of the algorithm. In the first one, both model estimation and validation have been conducted using measured data (power and temperature for estimation, irradiance and temperature for forecasting) in order to assess the performance of the estimation procedure net of errors due to inaccuracies of weather forecasts. In the second, meteorological predictions have been used both for model parameter fitting and generation forecasting. The latter scenario corresponds to a typical DSO use case.

7.1 Experiment set up

For the two experiments performed, the following data sets have been used, respectively:

• D1:

data from a PV plant P1 with nominal power $P_{nom}=$960\text{,}\mathrm{kW}$$p located in the campus of the University of Salento, in Monteroni di Lecce, Puglia, Italy (see [39] for details). Data, ranging from March 5th, 2012 to December 31st, 2013, consist of hourly samples of averaged measured power, air temperature and normal irradiance (the latter used only for comparison in the SRLS benchmark) ;

• D2:

data from a PV plant P2 with nominal power $P_{nom}=$920\text{,}\mathrm{kW}$$p located in Sardinia. Data, ranging from February 2nd, 2012, to May 1st, 2012, consist of hourly samples of averaged measured power, one day-ahead forecasts of air temperature and one day-ahead forecasts of normal irradiance. Information about the quality of such forecasts is reported in Table 1.

Therefore the data sets used for model estimation in the two cases are given by:

[TABLE]

where the sets of time indices $\mathcal{K}_{1}$ and $\mathcal{K}_{2}$ span the entire periods reported above for D1 and D2, respectively, with a sampling time $\tau_{s}=$1\text{,}\mathrm{h}$$, and $I^{cs}(k)$ is generated using (11),(12). Clearly, only time indices $k$ corresponding to hours of light were considered.

The initial parameter vector has been chosen according to the criteria in Section 4, i.e., $\hat{\mu}_{1}(0)=0.75\leavevmode\nobreak\ P_{nom}/1000$, $\hat{\mu}_{2}(0)=$-1.34\text{\times}{10}^{-4}$\cdot\hat{\mu}_{1}(0)$, and $\hat{\mu}_{3}(0)=$-3.25\text{\times}{10}^{-3}$\cdot\hat{\mu}_{1}(0)$. Concerning the panel orientation angles $(\zeta,\psi)$ used in (12), they have been chosen using a-priori knowledge: measurements of the panel angles and location, for P1 and P2, respectively. In particular, plant P1 is actually composed of two arrays with different orientations; for this plant an equivalent orientation has been estimated by averaging the respective angles, considering the nominal powers as weights. The parameters just described are summarized in Table 2.

Concerning the choice of $\epsilon$, it is worth recalling (see Remark 1) that in order for CS test 3 to reject uniformly cloudy data with a CCF $\beta\leq\beta_{0}$, $\epsilon$ can be chosen approximately as

[TABLE]

where $\hat{\mu}_{1}$ represents the currently available estimate of $\mu_{1}$. Therefore, we find it convenient to fix the CCF bound $\beta_{0}$ and adjust $\epsilon$ dynamically via (43) as soon as a new estimate $\hat{\mu}_{1}$ is computed. In this respect, we observe that the range of variability of the CCF depends on the climate of the macro-area where the plant is located, which is usually available. For the Italian case, typical values of the CCF range from $0.5$ to $1$ [40]. In the experiments of this section, we choose $\beta_{0}=0.9$. However, higher/smaller values of $\beta_{0}$ within the typical variability range make the CS detection algorithm more/less selective. Therefore, an evaluation of this effect is also in order.

7.2 Validation on measured data (D1)

The proposed method (denoted as CSD) has been evaluated with reference to day-ahead (DA) forecasts [24] by taking actual measurements of meteorological variables as the respective forecasts. The performance is compared with that of both the ODNP and the SRLS. Initialization data are summarized in Table 3.

The time evolution of the parameters estimated using CSD and SRLS algorithms are shown in Figure 2. Since the two algorithms use different data, namely theoretical irradiance for CSD and measured irradiance for SRLS, it is not surprising that parameters tend to slightly different values.

As far as the forecasting performance is concerned, all error measures on DA predictions were computed over the period starting from day $28$, in order to guarantee at least a rough adaptation of the model parameters. In Table 4 the performance indices achieved by the proposed CSD approach are compared with SRLS and ODNP. Errors computed on CSD and SRLS are comparable and clearly show better performance with respect to the ODNP. In Figure 3, the time evolution of the daily RMSE ($RMSE_{d}$) is, shown .

A visual representation of the algorithm behavior with special attention to clear-sky detection is shown in Figure 4. In those graphs, sequences of red markers denote time windows in which the measured power is detected as being generated under a clear-sky condition. The adaptation of model parameters is apparent from a comparison of the measured and predicted power in successive clear-sky periods. With reference to Figure 4, in day $8$ the first clear-sky window is detected: note that $\hat{P}^{cs}$ is much lower then $P^{m}$. During day $9$ the second clear-sky window is detected, in this case the model overestimates the actual generated power. On day $10$ the model fit has largely improved. The remaining plots show other three, non consecutive days: days $32$ and $419$ are completely clear-sky; day $91$ is a partially clear-sky day, in which about a half of the samples is rejected by the algorithm.

Finally, in Figure 5, DA forecasts provided by CSD and SRLS during three different days and under three different weather conditions are compared with the measures of generated power.

7.3 Influence of $\beta_{0}$

To show the influence of the choice of $\beta_{0}$ on the behavior of the algorithm, the following experiment on data set D1 has been performed. The presence of a certain amount of power data generated under uniform cloudiness has been simulated by scaling down a given fraction of the power data collected during days which appeared to be clear-sky by inspection. The power curve related to each of such days has been scaled by a factor ranging from $0.5$ to $0.9$. Three different datasets have been generated, each containing a different fraction of scaled data, as described in Table 5. For each data set, the model estimation experiment has been repeated several times by varying $\beta_{0}$ from $0.40$ to $0.95$ with steps of $0.05$.

For the sake of fairness, original data from D1 have been used in all cases to compute forecasting errors.

Figure 6 shows the percentage of generated power measurements detected as clear-sky by the algorithm in the different data sets for varying $\beta_{0}$. Figure 7 depicts the corresponding value of the $MAPE_{NP}$ on DA forecasts. When $\beta_{0}$ increases, the CSD algorithm becomes more selective. This fact is reflected in the $MAPE_{NP}$, which is lower in general for higher $\beta_{0}$. For given $\beta_{0}$, the error increases with the percentage of uniformly cloudy days. This phenomenon becomes less apparent as $\beta_{0}$ increases.

In Figure 8, the evolution of the parameter estimates performed on D12 for varying $\beta_{0}$ is reported. Parameter estimates tend to become almost stationary in all cases. Mean values and variances of $\hat{\mu}$ are reported in Table 6. Notice that $\hat{\mu}_{1}$ shows lower sensitivity to $\beta_{0}$ compared to $\hat{\mu}_{2}$ and $\hat{\mu}_{3}$.

Figure 9 depicts measured power and predicted clear-sky power during three uniformly cloudy days belonging to D$1_{2}$ for $\beta_{0}=0.9$. All data satisfy CS tests 1 and 2. CS test 3 is satisfied only for the last day, in which power data was obtained using a scaling factor greater or equal to 0.9.

7.4 Validation on predicted data (D2)

In this section a typical DSO scenario is reproduced, in which it is assumed that measurements of weather variables are not available at the plant site. Therefore, measured power and temperature forecasts are used to estimate the plant parameters, while irradiance and temperature forecasts are used to cast predictions of generated power. The algorithm configuration parameters are reported in Table 7.

In this scenario, the performance of the proposed method has been evaluated with reference to both day-ahead (DA) and hour-ahead (HA) forecasts, and compared with the performance achieved by SRLS and ODNP. Forecasting error measures are reported in Table 8 and Figure 10. While ODNP still has the worst performance indices, CSD performs even better then SRLS. However, it should be observed that forecasting errors in this case are to a large extent due to the quality of weather reports (see Table 1).

Three examples of DA forecast computed using CSD approach and SRLS during different weather conditions are shown in Figure 11.

7.5 Further remarks

With reference to Tables 4 and 8, it is important to observe that the normalized errors ($MAPE_{NP})$ computed on DA forecasts are below $10\%$, which demonstrates viability for network operation. Furthermore, the performance indices achieved by CSD are very close to those obtained by SRLS, i.e., via a PVUSA model estimated using measured irradiance.

Concerning the role of the tunable parameter $\beta_{0}$, results in Section 7.3 show that the estimate of the main power/irradiance gain $\mu_{1}$ is quite robust with respect to $\beta_{0}$, and moreover the parameter estimates $\hat{\mu}$ tend to converge regardless of the value of $\beta_{0}$. Even for small $\beta_{0}$, i.e., when the algorithm is not selective and CS test 3 is satisfied even for heavy uniform cloudiness, CSD is able to provide reasonably accurate forecasts. Increasing the values of $\beta_{0}$, the algorithm tends to reject more and more data measured under a uniformly cloudy sky, resulting in an improvement of the forecast quality.

The proposed algorithm has been implemented in Scilab [41]. Each iteration took on average less than one second on an i7 2.6 Ghz processor, thus demonstrating that the approach carries an extremely low computational burden.

8 Conclusions

In this paper, an efficient technique for estimating a forecasting model of photovoltaic power generation from limited information has been proposed. The approach is based on a set of tests performed on power data combined with a recursive estimation framework. It only exploits the time series of generated power and forecasts of temperature, the latter obtained from a meteorological service. The procedure especially fits the typical scenario where the network operator has no access to on-site measurements of irradiance and temperature, due to the large number of plants connected to the grid.

The algorithm has been extensively validated on two plants located in Italy, both on measured data and on forecasts of weather variables. The latter case reproduces a typical DSO scenario. Experiments worked out show very good forecasting performance, with limited computational burden.

Ongoing work addresses the aggregation of several plants covering large geographic areas. Due to a better quality of weather forecasts in this case, a significant accuracy improvement is expected. The integration of PV power generation forecasting in smart buildings and in microgrids will also be considered.