Model Estimation for Solar Generation Forecasting using Cloud Cover Data

TL;DR

This paper introduces a simple, efficient parametric model for forecasting solar power generation using cloud cover data, suitable for large-scale PV integration without needing local meteorological measurements.

Contribution

It proposes a novel model that leverages cloud cover data and power measurements, enabling large-scale PV forecasting without local meteorological sensors.

Findings

Model achieves accurate forecasts with low complexity.

Validated with both simulated and real data.

Suitable for large-scale PV plant integration.

Abstract

This paper presents a parametric model approach to address the problem of photovoltaic generation forecasting in a scenario where measurements of meteorological variables, i.e., solar irradiance and temperature, are not available at the plant site. This scenario is relevant to electricity network operation, when a large number of PV plants are deployed in the grid. The proposed method makes use of raw cloud cover data provided by a meteorological service combined with power generation measurements, and is particularly suitable in PV plant integration on a large-scale basis, due to low model complexity and computational efficiency. An extensive validation is performed using both simulated and real data.

Tables

Table 1. Table 1: Setup descriptions for scenario (i).

SID	Description
$0$	No noise added, raw data sampled at $15\min$
$1$	No noise added, data averaged over $1h$

Table 2. Table 3: Setup description for scenarios (iii), (iv) and (v).

SID	$6$	$7$	$8$	$9$	$10$	$11$	$12$
$3{\sigma }_N$	$0.0$	$0.0$	$0.0$	$0.0$	$0.0$	$0.0$	$0.3$
$3{\sigma }_P\left(\mathrm{kW}\right)$	$10$	$50$	$100$	$0$	$0$	$0$	$50$
$3{\sigma }_T\left(\mathrm{°C}\right)$	$0$	$0$	$0$	$1$	$3$	$5$	$3$

Table 3. Table 4: Correction terms η 2 subscript 𝜂 2 \eta_{2} and η 3 subscript 𝜂 3 \eta_{3} computed using the estimated parameter vectors reported in ( 42 ), ( 47 ) and ( 56 ).

${\hat{\eta }}_i$	${𝒮}_i$	N5	N6	L
${\hat{\eta }}_2$	$[-2.5\text{×}{10}^{-4},-1.9\text{×}{10}^{-5}]$	$-2.47\text{×}{10}^{-4}$	$-1.58\text{×}{10}^{-4}$	$-2.51\text{×}{10}^{-4}$
${\hat{\eta }}_3$	$[-4.8\text{×}{10}^{-3},-1.7\text{×}{10}^{-3}]$	$-2.06\text{×}{10}^{-4}$	$-1.73\text{×}{10}^{-3}$	$-10.1\text{×}{10}^{-3}$

Table 4. Table 5: Performance comparison of parametric models, ODNP, PVGM, CCD and MLP computed starting from day 57 57 57 . Results are grouped by model complexity and purpose in order to simplify the comparison: 5 5 5 -parameter model (N5) in red column, 6 6 6 -parameter models (N6 and MLP1) in blue columns, 11 11 11 -parameter models (L and MLP2) in green column, and benchmark models (ODNP, PVGM and CCD) in gray columns.

Performance Indeces		N5	N6	MLP1 Mean (Variance)	L	MLP2 Mean (Variance)	ODNP	PVGM	CCD
DA Forecast	$RMSE$ ($\mathrm{kW}$)	$109$	$110$	$133$ $\left(2.05\right)$	$110$	$127$ $\left(19.0\right)$	$227$	$283$	$121$
	$MAPE$	$50.0\%$	$50.7\%$	$64.6\%$ $\left(5.17\right)$	$48.8\%$	$66.9\%$ $\left(31.8\right)$	$87.4\%$	$165.3\%$	$61.4\%$
	$MBE$ ($\mathrm{kW}$)	$-1.275$	$4.108$	$-3.42$ $\left(1.38\right)$	$4.455$	$-1.15$ $\left(37.7\right)$	$27.2$	$127.4$	$11.2$
	$R^2$	$0.828$	$0.827$	$0.746$ $\left(2.98\text{×}{10}^{-5}\right)$	$0.825$	$0.770$ $\left(2.71\text{×}{10}^{-4}\right)$	$0.254$	$-0.157$	$0.787$
	$NRMSE$	$0.414$	$0.416$	$0.511$ $\left(2.94\text{×}{10}^{-5}\right)$	$0.418$	$0.484$ $\left(2.28\text{×}{10}^{-4}\right)$	$0.864$	$1.076$	$0.462$
	$RMS{E}_{NP}$	$0.119$	$0.119$	$0.146$ $\left(2.41\text{×}{10}^{-6}\right)$	$0.120$	$0.139$ $\left(1.88\text{×}{10}^{-5}\right)$	$0.248$	$0.308$	$0.132$
	$MAP{E}_{NP}$	$7.51\%$	$7.78\%$	$9.85\%$ $\left(1.70\text{×}{10}^{-2}\right)$	$7.71\%$	$9.28\%$ $\left(1.24\text{×}{10}^{-1}\right)$	$15.28\%$	$23.96\%$	$9.49\%$
HA Forecast	$RMSE$ ($\mathrm{kW}$)	$103$	$104$	$129$ $\left(6.79\right)$	$105$	$125$ $\left(14.4\right)$	-	$223$	$120$
	$MAPE$	$45.1\%$	$46.7\%$	$66.9\%$ $\left(11.5\right)$	$42.7\%$	$75.8\%$ $\left(183\right)$	-	$177.3\%$	$61.4\%$
	$MBE$ ($\mathrm{kW}$)	$11.74$	$11.97$	$-2.03$ $\left(2.63\right)$	$12.53$	$-4.01$ $\left(45.1\right)$	-	$62.9$	$24.8$
	$R^2$	$0.846$	$0.841$	$0.756$ $\left(9.73\text{×}{10}^{-5}\right)$	$0.839$	$0.774$ $\left(1.89\text{×}{10}^{-4}\right)$	-	$0.274$	$0.788$
	$NRMSE$	$0.393$	$0.399$	$0.498$ $\left(9.88\text{×}{10}^{-5}\right)$	$0.401$	$0.477$ $\left(1.77\text{×}{10}^{-4}\right)$	-	$0.852$	$0.460$
	$RMS{E}_{NP}$	$0.112$	$0.114$	$0.142$ $\left(8.03\text{×}{10}^{-6}\right)$	$0.115$	$0.136$ $\left(1.44\text{×}{10}^{-5}\right)$	-	$0.243$	$0.131$
	$MAP{E}_{NP}$	$6.88\%$	$7.07\%$	$9.05\%$ $\left(3.28\text{×}{10}^{-2}\right)$	$7.16\%$	$8.83\%$ $\left(5.28\text{×}{10}^{-2}\right)$	-	$18.66\%$	$9.73\%$

Table 5. Table 6: Computational times. Results are grouped by complexity and purpose in order to simplify the comparison: 5 5 5 parameters (N5) in red column, 6 6 6 parameters (N6 and MLP1) in blue columns, 11 11 11 parameters (L and MLP2) in green column, and benchmark models (ODNP, PVGM and CCD) in gray columns.

Performance Indices	N5	N6	MLP1 Mean (Variance)	L	MLP2 Mean (Variance)	PVGM	CCD
Sim. Time ($s$)	$7.87$	$7.86$	$127$ $\left(285\right)$	$7.58$	$130$ $\left(433\right)$	$9.58$	$10.22$

Equations

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

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

$$\begin{array}[]{l}\eta_{2}\in{\cal S}_{2}=\left[$-2.5\text{\times}{10}^{-4}$,$-1.9\text{\times}{10}^{-5}$\right],\\ \eta_{3}\in{\cal S}_{3}=\left[$-4.8\text{\times}{10}^{-3}$,$-1.7\text{\times}{10}^{-3}$\right].\end{array}$$

$$I(N)=C(N)\cdot I^0,$$

$$C(N)=1+{\mu }_{4}N+{\mu }_5{N}^2.$$

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

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

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

$$P=P(I^0,T,N)=({\mu }_1+{\mu }_{2}I(N)+{\mu }_{3}T)I(N),$$

$$I(N)=\left(1+{\mu }_{4}N+{\mu }_5{N}^2\right)I^0.$$

$$\mu ={\left[\begin{array}{ccccc}{\mu }_1& {\mu }_2& {\mu }_3& {\mu }_4& {\mu }_5\end{array}\right]}^T,$$

$$P={\phi }^T(I^0,T,N)\theta (\mu ),$$

$$\begin{array}[]{rcl}\theta(\mu)&=&{\left[\begin{array}[]{ccccc}\theta_{1}(\mu)&\theta_{2}(\mu)&\theta_{3}(\mu)&\theta_{4}(\mu)&\theta_{5}(\mu)\par\end{array}\right.}\\ &&{\left.\begin{array}[]{cccccc}\theta_{6}(\mu)&\theta_{7}(\mu)&\theta_{8}(\mu)&\theta_{9}(\mu)&\theta_{10}(\mu)&\theta_{11}(\mu)\end{array}\right]^{T}}\\ \nobreak\leavevmode\hfil&=&{\left[\begin{array}[]{ccccc}\mu_{1}&\mu_{1}\mu_{4}&\mu_{1}\mu_{5}&\mu_{2}&2\mu_{2}\mu_{4}\par\end{array}\right.}\\ &&{\left.\begin{array}[]{cccccc}\mu_{2}\mu_{4}^{2}+2\mu_{2}\mu_{5}&2\mu_{2}\mu_{4}\mu_{5}&\mu_{2}\mu_{5}^{2}&\mu_{3}&\mu_{3}\mu_{4}&\mu_{3}\mu_{5}\end{array}\right]^{T},\par\par}\end{array}$$

$$\begin{array}[]{l}\varphi(I^{0},T,N)={\left[\begin{array}[]{ccccc}I^{0}&I^{0}N&I^{0}N^{2}&{I^{0}}^{2}&{I^{0}}^{2}N\end{array}\right.}\\ {\left.\begin{array}[]{cccccc}{I^{0}}^{2}N^{2}&{I^{0}}^{2}N^{3}&{I^{0}}^{2}N^{4}&TI^{0}&TI^{0}N&TI^{0}N^{2}\end{array}\right]^{T}}.\end{array}$$

$$\mu ={\left[\begin{array}{ccccc}{\mu }_1& {\mu }_2& {\mu }_3& {\mu }_4& {\mu }_5\end{array}\right]}^T,$$

$$\begin{array}[]{rcl}\theta(\mu)&=&{\left[\begin{array}[]{ccccc}\mu_{1}&\mu_{1}\mu_{4}&\mu_{1}\mu_{5}&\mu_{2}&2\mu_{2}\mu_{4}\par\end{array}\right.}\\ &&{\left.\begin{array}[]{cccccc}\mu_{2}\mu_{4}^{2}+2\mu_{2}\mu_{5}&2\mu_{2}\mu_{4}\mu_{5}&\mu_{2}\mu_{5}^{2}&\mu_{3}&\mu_{3}\mu_{4}&\mu_{3}\mu_{5}\end{array}\right]^{T},\par\par}\end{array}$$

$$\mu ={\left[\begin{array}{cccccc}{\mu }_1& {\mu }_2& {\mu }_3& {\mu }_4& {\mu }_5& {\mu }_6\end{array}\right]}^T,$$

$$\begin{array}[]{rcl}\theta(\mu)&=&{\left[\begin{array}[]{ccccc}\mu_{1}&\mu_{1}\mu_{4}&\mu_{1}\mu_{5}&\mu_{2}&2\mu_{6}\par\end{array}\right.}\\ &&{\left.\begin{array}[]{cccccc}\mu_{4}\mu_{6}+2\mu_{2}\mu_{5}&2\mu_{5}\mu_{6}&\mu_{2}\mu_{5}^{2}&\mu_{3}&\mu_{3}\mu_{4}&\mu_{3}\mu_{5}\end{array}\right]^{T},\par\par}\end{array}$$

$$P={\phi }^T(I^0,T,N)\theta (\mu ).$$

$$P={\phi }^T(I^0,T,N)\theta$$

$$\mu (k+1)$$

$$P(k)$$

$$\hat{\mu }(k)=\hat{\mu }(k-1)+K(k)\left[P^m(k)-{\varphi }^T(k)\theta (\hat{\mu }(k-1))\right],$$

$$\text{vspace}\ast -2mmK(k)=\frac{R(k-1)H^T(k)}{H(k)R(k-1)H^T(k)+r},$$

$$R(k)=\left(I-H^T(k)K^T(k)\right)R(k-1),$$

$$H(k)=\frac{\partial }{\partial \mu }{\left[{\varphi }^T(k)\theta (\mu )\right]}_{\mu =\hat{\mu }(k-1)}.$$

$$L(k)=V(k)\varphi (k),$$

$$V(k)=V(k-1)-\frac{V(k-1)\varphi (k){\varphi }^T(k)V(k-1)}{1+{\varphi }^T(k)V(k-1)\varphi (k)}.$$

$$\hat{\theta }(k)=\hat{\theta }(k-1)+L(k)\left[P^m(k)-{\varphi }^T(k)\hat{\theta }(k-1)\right].$$

$$\begin{array}[]{rcl}\theta(\mu)&=&{\left[\begin{array}[]{ccccc}\hat{\mu}_{1}(0)&\hat{\mu}_{1}(0)\hat{\mu}_{4}(0)&\hat{\mu}_{1}(0)\hat{\mu}_{5}(0)&\hat{\mu}_{2}(0)&2\hat{\mu}_{2}(0)\hat{\mu}_{4}(0)\par\end{array}\right.}\\ &&{\left.\begin{array}[]{ccc}\hat{\mu}_{2}(0)\hat{\mu}_{4}(0)^{2}+2\hat{\mu}_{2}(0)\hat{\mu}_{5}(0)&2\hat{\mu}_{2}(0)\hat{\mu}_{4}(0)\hat{\mu}_{5}(0)&\hat{\mu}_{2}(0)\hat{\mu}_{5}(0)^{2}\\ \end{array}\right.}\\ &&{\left.\begin{array}[]{ccc}\hat{\mu}_{3}(0)&\hat{\mu}_{3}(0)\hat{\mu}_{4}(0)&\hat{\mu}_{3}(0)\hat{\mu}_{5}(0)\end{array}\right]^{T}.}\end{array}$$

$$\hat{P}(j|k;q)=\left\{\begin{array}[]{ll}\varphi^{T}(I^{0}(j),\hat{T}(j|k),\hat{N}(j|k))\cdot\theta(\hat{\mu}(q))&{\rm(N5,\leavevmode\nobreak\ N6)}\\ \varphi^{T}(I^{0}(j),\hat{T}(j|k),\hat{N}(j|k))\cdot\hat{\theta}(q)&{\rm(L)}.\end{array}\right.$$

Full text

††journal: Solar Energy

1 Introduction

A major challenge in the integration of renewable energy sources into the grid is that power generation is intermittent, difficult to control, and strongly depending on the variation of weather conditions. For these reasons, forecasting of renewable distributed generation has become a fundamental need to grid operators. In this respect, solar generation forecasts on multiple time horizons are needed to satisfy grid constraints and demand. In particular, short-term forecasts are required for the purposes of power plant operation, grid balancing, real-time unit dispatching, automatic generation control, and energy trading. On the contrary, longer-term forecasts are of interest to Distribution System Operators (DSO) and Transmission System Operators (TSO) for unit commitment, scheduling and for improving balance area control (see, e.g., [1, 2, 3, 4] and references therein).

Concerning solar power generation, much attention has been paid to the problem of obtaining accurate day-ahead and hour-ahead forecasts of solar irradiance and/or generated power (see [5, 6] for a comprehensive overview on the subject). Most contributions focus on solar irradiance forecasting [7]. Widely adopted approaches are based on Artificial Neural Networks (ANNs) [8, 9, 10, 11] or Support Vector Machines [12] with different types of input data. Alternatively, classical linear time series forecasting methods are used in [13, 14, 15], where the considered time series is typically the global horizontal irradiance (GHI) normalized with a clear-sky model (see [16, 17] for a comprehensive review). Global radiation forecasts are then fed along with temperature forecasts to a simulation model of the plant [18] to compute the prediction of power generation. In any case, computing reliable generation forecasts from predicted meteorological variables hinges upon the availability of a reliable model of the plant, be it physical or estimated from data.

Unfortunately, in many practical scenarios, neither reliable plant models, nor direct on-site measurements of solar irradiance and other meteorological variables such as temperature are available. This is the typical case of a a DSO dealing with hundreds or thousands of distributed heterogeneous and independently-operated solar plants; in this case, the only available plant data is represented by generated power measurements provided by electronic meters. The contribution of this paper focuses on the problem of estimating reliable generation models in such limited information contexts. In [19, 20], two methods are proposed to estimate the parameters of the PVUSA model [21] of a PV plant via a recursive framework based only on measures of generated power and temperature forecasts. This solution is not always data-efficient, since generation data collected during cloudy days are not exploited.

In order to obtain more accurate results, power measurements can be combined with further data coming from a weather service. Such data are typically averaged over large geographic zones and therefore they may be scarcely informative for a specific spot, yet they may provide useful information when the goal is to address the aggregation of multiple plants over a macro area. In this respect, cloud cover measurements derived from satellite imaging are an exception, since they can be made available for specific locations with good spatial and temporal resolution (up to $2.5\leavevmode\nobreak\$\mathrm{km}$\times 2.5\leavevmode\nobreak\$\mathrm{km} and $30\leavevmode\nobreak\ $\mathrm{min}, respectively) [22]. Based on such data, suitable models can be estimated for irradiance forecasting purposes (see [23, 24, 25] for details). Irradiance forecasts can be used as inputs to a plant generation model in order to provide energy production forecasts. Of course, the latter step requires that a reliable model of power generation from irradiance be previously estimated.

In this paper, we present a novel approach for direct forecasting of PV plant power generation from cloudiness data. To this purpose, a class of parametric models is introduced which efficiently exploits the PVUSA model [21] and the notion of Cloud Cover Factor (CCF) [26, 27] in a limited information scenario. More specifically, the model parameters are estimated using generated power data combined with additional information provided by a meteorological service for the area where the plant is located. Such data consist of a time series of raw cloud cover and temperature reports. For the proposed models, estimation procedures based on Recursive Least Squares (RLS) and the Extended Kalman Filter (EKF) [28], are devised. The properties and the performance of the proposed method are demonstrated both in a simulated scenario and on experimental data from a plant currently in operation. Preliminary results leading to this paper were presented in [29].

The manuscript is organized as follows. In Section 2 we introduce the modelling tools used; in Section 3 the proposed models are presented. Estimation procedures are presented in Section 4. Section 5 addresses the relevant forecasting problems, while performance evaluation is discussed in Section 6. Simulation results are illustrated in Section 7. Experimental validation results and related discussion are reported in Section 8. Finally, conclusions are drawn in Section 9.

2 Models and methods

2.1 The PVUSA photovoltaic plant model

A PV plant can be efficiently modelled using the PVUSA model [21], 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.

Model (1) is linear and parsimonious in terms of number of parameters. Despite its simplicity, very good accuracy is obtained when this model is fit to real measured data. Moreover, due to their reliability, temperature forecasts taken from a meteorological service can be used efficiently in place of actual measurements in order to estimate the model parameters, as shown in [19]. This cannot be done for irradiance, as forecasts are much less reliable.

We find it convenient to express (1) also in the form

[TABLE]

where $\mu_{1}$ plays the role of the main power/irradiance gain of the plant, while $\eta_{2}=\mu_{2}/\mu_{1}$ and $\eta_{3}=\mu_{3}/\mu_{1}$ introduce correction terms. It is worth noticing that $\eta_{2}$ and $\eta_{3}$ in (2) are characterized by a narrow range of variability among different PV technologies. Indeed, typical values of these parameters are given by [21]:

[TABLE]

2.2 Cloud Cover Index and Cloud Cover Factor

The Cloud Cover Index (CCI), here denoted by $N$, refers to the fraction of the sky obscured by clouds when observed from a particular location. In meteorology, $N$ is measured in okta [30], i.e., an integer ranging from $0\text{\,}\mathrm{o}\mathrm{k}\mathrm{t}\mathrm{a}$ (clear sky) up to $8\text{\,}\mathrm{o}\mathrm{k}\mathrm{t}\mathrm{a}$ (completely overcast). Weather services often provide this information in terms of percentage of covered sky: from $0.0$ (clear sky) to $1.0$ (completely overcast) in steps of $0.1$.

A widely investigated issue is the relationship between the CCI and the global irradiance on a surface. In [27], a simple way of estimating solar radiation in cloudy days is to multiply the clear-sky irradiance by the so-called Cloud Cover Factor (CCF), which is defined as a function of the CCI. Let $I^{0}$ be the clear-sky irradiance on a given surface. It turns out that the global irradiance on the same surface is given by

[TABLE]

where $C(N)$ is the CCF, which can be expressed in the form

[TABLE]

In (5), the parameters $\mu_{4}$ and $\mu_{5}$ depend on the local climate. In particular, $\mu_{5}$ is always negative. Depending on the sign of $\mu_{4}$, the curve $C(N)$ may be non-monotonic (see Fig. 1).

Indeed, in some locations $C(N)$ may rise above 1 (i.e., the clear-sky value) under slightly cloudy sky (e.g, $0<N<0.3$), since the diffuse radiation becomes significant with middle-low cloudiness. Moreover, while the CCF is independent of solar elevation, it is subject to seasonal variations. Alternative CCF models exist in the literature [26].

The CCF is usually estimated for a given location using on-site irradiance measurements [31]. It is worth to remark, however, that the method presented in this paper does not rely on such measurements.

2.3 Clear-sky irradiance model

In this paper we need to exploit an estimate of the clear-sky irradiance on a given surface, as opposed to the GHI. As a first step, we make use of the Heliodon simulator model [32]. This model allows to compute the theoretical 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 surface $I^{cs}$ can be derived from $I^{cs,n}$ and surface orientation 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 an horizontal surface is obtained.

In this study, we make the realistic assumption that the exact orientation of the PV panel surfaces of the considered plant is not known a-priori. However, it is reasonable to assume that the plant is efficiently oriented for the specific latitude according to, e.g., the guidelines given in [33]. Therefore, we will use as the reference value of the theoretical clear-sky irradiance $I^{0}$ for a specific plant, the value of (7) for $(\zeta,\psi)=(\zeta_{0},\psi_{0})$, where $\zeta_{0}$ and $\psi_{0}$ are taken from the above guidelines, i.e.,

[TABLE]

The choice of the model in (6) is mainly due to its simplicity. More accurate models taking into account atmospheric parameters such as aerosol optical depth, ozone, and water vapor, could be used upon availability of such data in order to improve the performance of the forecasting method proposed in this work.

3 Plant generation model

The plant model (1) and the irradiance model (4)$-$(5) introduced in the previous section can be combined in order to obtain an expression of the generated power as a function of $N$, $T$, and $I^{0}$, that is

[TABLE]

where

[TABLE]

Therefore, the overall model is defined by the parameter vector

[TABLE]

and (9)$-$(10) can be rewritten in a standard regression form as

[TABLE]

where $\theta(\mu)$ and $\varphi(I^{0},T,N)$ are given by

[TABLE]

[TABLE]

Notice that the model in (11)-(13) is not linear in the parameter vector $\mu$. In order to estimate model parameters, two approaches will be pursued in the following by comparing their respective performances:

• (N)

Keep the parameter vector size to a minimum and derive a nonlinear estimation algorithm based on the Extended Kalman Filter (EKF);

• (L)

Introduce an overparameterization by disregarding the parameterization of $\theta(\mu)$ in $\mu$ and assuming $\theta(\cdot)$ to be a vector of independent parameters, thus making the model linear and amenable to standard least squares estimation.

The former approach has the advantage of keeping the model parsimonious and improving its identifiability features, while the latter is introduced with the purpose of making the estimation task simpler.

3.1 Nonlinear (N) models

The first idea proposed here is to tackle the nonlinear estimation problem exploiting the EKF. To this purpose, let us consider model (11)$-$(13) for the following two choices of the independent parameter vector $\mu$:

N5 model

[TABLE]

[TABLE]

which is the original nonlinear model (11)$-$(13) characterized by a minimal number of parameters, and 2. 2.

N6 model

[TABLE]

[TABLE]

which implies a slight overparameterization with respect to model N5, amounting to treating $\mu_{2}$, $\mu_{4}$ and $\mu_{2}\mu_{4}=\mu_{6}$ as independent parameters.

For either choice, the model takes the form (11)

[TABLE]

3.2 Linear overparameterized (L) model

This model is derived by treating $\theta(\mu)$ in (12) as a vector $\theta$ of 11 independent parameters. This way model (11) becomes

[TABLE]

and is clearly linear in the parameters. For this model, a standard Recursive Least Squares (RLS) algorithm can be used to estimate $\theta$.

4 Model estimation

For the models introduced in Sections 3.1 and 3.2, an on-line update of the estimate of the parameter vector $\mu$ (or $\theta$) is performed. To the purpose of illustrating the procedure, let us define the following quantities:

$\tau_{s}$: sampling time (mins);

• 2.

$k$: discrete time index;

• 3.

$d$: generic day;

• 4.

$\tau_{d}(k)$: time of day (TOD) corresponding to index $k$ (at the given longitude);

• 5.

$I^{0}(k)$: theoretical clear-sky solar irradiance at time index $k$ at the plant site, computed according to (6)-(8);

• 6.

${\cal K}_{d}=\{\underline{k}_{d},\dots,\overline{k}_{d}\}$: set of time indices pertaining to light hours in day $d$, i.e., the theoretical irradiance (6)-(8) is nonzero for $\tau_{d}(k)$ corresponding to $\underline{k}_{d}\leq k\leq\overline{k}_{d}$, and moreover $\underline{k}_{d+1}=\overline{k}_{d}+1$, that is, only daylight time indices are considered;

• 7.

$\hat{\mu}(k)$ [$\hat{\theta}(k)$]: estimate of the parameter vector at time $k$, after processing the last data sample, being $\hat{\mu}(0)$ [$\hat{\theta}(0)$] the initial guess;

• 8.

${W}(k)$: weather report for time $k$ relative to the the plant area provided by a meteorological service. More specifically, $W(k)=\{N(k),T(k)\}$, where $N(k)$ and $T(k)$ are the cloud cover and temperature reports, respectively111$T(k)$ can be represented by either a real-time estimate of the plant site temperature or even some forecast thereof, computed in advance. Indeed, it has been shown [19] that the reliability of day-ahead temperature forecasts is comparable to that of actual measurements for model estimation purposes. Obviously, when on-site measurements of temperature $T^{m}(k)$ are available, then $T(k)$ should be taken equal to $T^{m}(k)$. In the experimental section of this paper, data from day-ahead forecasts are considered..

• 9.

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

• 10.

$D(k)=\{P^{m}(k),W(k)\}=\{P^{m}(k),N(k),T(k)\}$: data available at time $k$;

• 11.

$\phi(k)=\varphi(I^{0}(k),T(k),N(k)$): regression vector at time $k$ (see (13));

The estimation procedure is recursive. At each time step $k$, a new power-weather data sample $\{P^{m}(k),W(k)\}$ is acquired, the algorithm is run and a current estimate $\hat{\mu}(k)$ [$\hat{\theta}(k)$] of the parameter vector is computed according to the criteria outlined below.

N5 and N6 models. For the purpose of estimating $\mu$ we consider the standard extended Kalman filter (EKF) formulation defined by the following state-space model [28]:

[TABLE]

where $\mu(k)$ represents the state vector, $P(k)$ is the output, and $w(k)$ is a zero-mean Gaussian white noise with given variance $r$ which describes both the measurement noise and the model uncertainty. According to the standard EKF theory, the following update law for the parameter estimate $\hat{\mu}(k)$ is obtained:

[TABLE]

where

[TABLE]

and the estimate covariance matrix $R(k)$ is updated at each step according to the rule

[TABLE]

being

[TABLE]

As in common practice, the initial estimate covariance matrix $R(0)$ is computed as $R(0)=l(0)\cdot\cal I$. The higher the confidence on the initial guess $\hat{\mu}(0)$, the lower $l(0)$ can be chosen. The initialization of the parameter vector estimate $\hat{\mu}(0)$ is discussed at the end of this section.

• 2.

L model. A recursive least-squares (RLS) estimation step is performed, i.e., we compute

[TABLE]

where $V(k)$ is the weight matrix at time $k$, which in turn is recursively computed according to the standard RLS update law:

[TABLE]

Finally the current parameter estimate is updated according to the law

[TABLE]

As in the usual RLS practice, the initial weight matrix $V(0)$ is set to $V(0)=l(0)\cdot\leavevmode\nobreak\ \cal I$, where $\cal I$ is the identity matrix and $l(0)>0$ is chosen according to the confidence given to the initial parameter guess $\hat{\theta}(0)$ (higher $l(0)$ meaning less confidence).

As far as the initialization of the parameter vector is concerned, a reliable guess for the initial values $\hat{\mu}(0)$ in model N5 can be computed from the following guidelines.

A good guess for the plant gain $\mu_{1}$ is given by $\hat{\mu}_{1}(0)=P_{nom}/1000$ where $P_{nom}$ is the nominal plant power in $\mathrm{kW}$ [19];

• 2.

the initial values $\hat{\mu}_{2}(0)$ and $\hat{\mu}_{3}(0)$ can be chosen such that $\eta_{2}(0)=\hat{\mu}_{2}(0)/\hat{\mu}_{1}(0)$ and $\eta_{3}(0)=\hat{\mu}_{3}(0)/\hat{\mu}_{1}(0)$ are the central points of the respective intervals in (3);

• 3.

$\hat{\mu}_{4}(0)$ and $\hat{\mu}_{5}(0)$ may be selected as the average values of the corresponding parameters for the climate of the macro-area where the plant is located, which are usually available. For example, for the Italian regions considered in the experimental section of this paper such values are estimated in [34].

Given their respective definitions, for the N6 model it can be assumed that $\hat{\mu}_{6}(0)=\hat{\mu}_{2}(0)\hat{\mu}_{4}(0)$, while for the L model a natural choice (see (12)) is given by

[TABLE]

5 Forecasting

In this section, the performance of models N5, N6 and L, proposed in Section 3 is evaluated on the widely used Day-Ahead (DA) and Hour-Ahead (HA) forecasts [35]. In order to suitably define them in our context, let $d$ and $k$ be the generic day and time instant, respectively, in which a forecast is supposed to be computed and submitted. For a given time instant $j\geq k$, let $\hat{W}(j|k)=\{\hat{N}(j|k),\hat{T}(j|k)\}$ denote the weather forecast {cloud cover, temperature} relative to time $j$ available at time $k$. We denote by $\hat{P}(j|k;q)$ the prediction of generated power for time instant $j$, computed at time $k$ using the parameter vector estimate available at time $q\leq k$ according to the appropriate model, i.e.,

[TABLE]

5.1 Day-Ahead forecast

The day-ahead 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. Therefore, the day ahead forecast is provided 19 to 42 hours prior to the operating time. In this respect, we recall that the vast majority of conventional generation is scheduled in the DA market. For a given day $d$, the DA forecast is submitted at time instant $k_{0}^{d}$ corresponding, e.g., to $6$ am and consists of the time series given by

[TABLE]

Such a forecast is computed on the basis of the last parameter estimate $\hat{\mu}(\overline{k}_{d-1})$ [$\hat{\theta}(\overline{k}_{d-1})$] available the day before (see Figure 2 for the meaning of $k_{0}^{d}$ and the forecast horizon).

5.2 Hour-Ahead forecast

The hour-ahead 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. Therefore, such a forecast, assuming it is computed at the time index $k$ just before submission, is given by the time series

[TABLE]

where ${\cal K}_{HA}(k)$ denotes the set of time indices pertaining to the 7-or-less-hour horizon starting at the beginning of the relevant operating hour (see Figure 3). Note that $\hat{\mu}(k)$ [$\hat{\theta}(k)$] represents the most recent parameter estimate available.

6 Performance evaluation

In this section we introduce the performance assessment indices that will be used to evaluate the fitness of models to data and the efficacy of the proposed methods in the forecasting problems introduced in the previous section. It is worth to remark that a single estimated model is used to provide both DA and HA forecasts.

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 for specific problems (e.g., DA or HA forecasting) will be provided in the sequel.

Let us consider a data set $\cal P$

[TABLE]

where $\cal K$ is a set of time indices of cardinality $K$. The set ${\cal P}$ is composed of pairs $\left(P^{m}(j),\,\hat{P}(j)\right)$ with $P^{m}(j)>0,\leavevmode\nobreak\ \hat{P}(j)>0$, where $\hat{P}(j)$ represents the forecasted power and $P^{m}(j)$ the corresponding measured value. We consider the following standard error measures pertaining to ${\cal P}$:

Root-Mean-Square Error (RMSE)

[TABLE]

• 2.

Mean-Bias-Error (MBE)

[TABLE]

• 3.

Mean-Absolute-Percentage-Error (MAPE)

[TABLE]

• 4.

Determination coefficient

[TABLE]

where $\overline{P}=\frac{1}{K}\sum_{j\in{\cal K}}P^{m}(j)$ is the mean of the measured power over the data set. Notice that here it is implicitly assumed that \sum_{j\in{\cal K}}\left(P^{m}(j)-\hat{P}(j)\right)^{2}\leq\sum_{j\in{\cal K}}\Bigl{(}P^{m}(j)-\overline{P}\Bigr{)}^{2} and in this case $1-R^{2}$ is the so-called unexplained square error.

• 5.

Normalized RMSE (NRMSE)

[TABLE]

In addition to the above standard statistical indices, we also consider the following two error measures, which are of practical interest for network operation, as they are referred to the nominal plant power $P_{nom}$:

Normalized RMSE w.r.t. nominal power ($RMSE_{NP}$)

[TABLE]

• 2.

Normalized MAPE w.r.t. nominal power ($MAPE_{NP}$)

[TABLE]

6.2 One-day-ahead naive predictor

As an additional evaluation tool, the performance indices achieved using the proposed approach will be compared to those obtained using the so-called One-Day-ahead Naive Predictor (ODNP), 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.

7 Validation on simulated data

Before assessing the performance of the proposed approach on measured data sets, in this section we present preliminary tests on simulated noisy data in order to highlight the properties of the proposed approach. Nominal weather data (i.e., cloud cover $N(k)$ and temperature $T(k)$) have been generated using an empirical procedure which reproduces a realistic cloud cover and temperature seasonal trend in Italy. The corresponding generated power time series $P^{m}(k)$ is then computed according to models (1) and (4) using the values of the physical parameters:

[TABLE]

The plant nominal power $P_{nom}$ is assumed equal to $920\text{\,}\mathrm{kW}$. Generated data are assumed to span the period January 1st (day $\underline{d}=1$) to December 31st (day $\overline{d}=365$) and are sampled with a time step $\tau_{s}=15$ minutes. The overall data set is then given by

[TABLE]

where the set of time indices $\cal K$ spans the time horizon from the starting day $\underline{d}$ to the final day $\overline{d}$. As previously stated, only the indices $k$ corresponding to hours of light (i.e., for which the theoretical clear-sky irradiance is nonzero) were considered. A model estimation step is therefore performed at each instant $k$ using the data sample $D(k)$.

7.1 Scenarios

Five different simulation scenarios (i)$-$(v) have been considered. Within each scenario, different noise and data processing conditions identified by a Setup ID (SID) have been simulated. The scenarios are arranged as follows:

1. (i)

Neither noise nor quantization is added. Raw data (SID [math]) and data averaged over $1\text{\,}\mathrm{h}$ (SID $1$) (see Table 2) were considered. Averaging is introduced in order to mimic the typical setting where hourly updates of weather and generation data are provided by the meteorological service and the DSO, respectively.

2. (ii)

A zero-mean Gaussian white noise with variance ${\sigma^{2}_{N}}$ is added to $N(k)$, $N(k)$ is then quantized so that $N(k)\in\left\{0,\,0.1,\,0.2,\,\ldots,\,1.0\right\}$, resulting in a noise variance $\sigma^{2}_{N,tot}$ (SID $2-5$ in Table 2). Data is finally averaged over $1\text{\,}\mathrm{h}$. 3. (iii)

A zero-mean Gaussian white noise with variance $\sigma_{P}^{2}$ is added to $P^{m}(k)$, $N(k)$ is quantized and data are averaged over $1\text{\,}\mathrm{h}$ (SID $6-8$, Table 3). 4. (iv)

A zero-mean Gaussian white noise with variance $\sigma_{T}^{2}$ is added to $T(k)$, $N(k)$ is quantized and data are averaged over $1\text{\,}\mathrm{h}$ (SID $9-11$, Table 3). 5. (v)

A zero-mean Gaussian white noise is added to all the regression variables, $N(k)$ is quantized and data are averaged over $1\text{\,}\mathrm{h}$ (SID $12$, Table 3).

Notice that here data averaged over $1\text{\,}\mathrm{h}$ means that for each hour $h$ there is only one sample computed as mean of the data acquired during $h$.

The performance of the proposed method on the scenarios introduced above has been evaluated with reference to the day-ahead (DA) forecast introduced in Section 5. For computing the forecasts, weather data from the data set $\cal D$ have been employed. In particular, at each time instant $k_{0}^{d}$ corresponding to 6 am on each day, we evaluate the DA forecast $\hat{P}^{DA}(d)$ according to (25)-(26) with

[TABLE]

The evaluation of the performance indices starts at day $d_{0}=18$, in order to first ensure a rough degree of adaptation of the model parameters. The following performance indices were considered:

$RMSE_{d}$: the RMSE relative to the forecast pertaining to each day $d$, i.e., (29) evaluated for $\hat{P}(j)=\hat{P}(j|k_{0}^{d};\overline{k}_{d-1})$, ${\cal K}={\cal K}_{d}$;

• 2.

RMSE, MBE, MAPE, and determination coefficient $R^{2}$, evaluated over the whole simulation, i.e., over the union of all DA forecasts from day $d_{0}$ to day $\overline{d}$.

For each SID, $10$ simulations were performed, in which $P^{m}(k)$, $T(k)$ and $N(k)$ were obtained by using different noise realizations. The average of performance indices over all the simulations were considered.

In each scenario, the initial values $\hat{\mu}_{1}(0)$, $\hat{\mu}_{2}(0)$, $\hat{\mu}_{3}(0)$, $\hat{\mu}_{4}(0)$ and $\hat{\mu}_{5}(0)$ are chosen as $75\%$ of the real values in (37), while $\psi$ and $\zeta$ are equal to their nominal values. $l(0)$ has been initialized to $0.01$ for both RLS algorithms and EKF, $r$ has been initialized to $10^{4}$, which is a conservative choice in view of the variance of the noise added to power measurements.

7.2 Results

Figure 4 shows one instance of the parameter estimate evolution for the N5 model for SID [math] and SID $1$. Using noise-free raw data (SID 0) the parameters converge exactly to the real values (Figure 4a). In case of SID $1$ the parameters converge as well with slight deviations of $\hat{\mu}_{4}$ and $\hat{\mu}_{5}$ from the nominal values (Figure 4b). Note that in this case, the parameters take considerably longer to converge.

Figures 5a, 5b and 5c show three examples of parameter estimate evolution for the N5 model under the largest noise variance on $N$, $P$ and $T$, respectively, Figure 5d shows the same plot when noise is added to all the regression variables. Parameters $\hat{\mu}_{4}$ and $\hat{\mu}_{5}$ are apparently quite sensitive to the noise added to $N$. Indeed, under the conditions described by the SID $5$, for which $\sigma_{N,tot}=0.441$, $\hat{\mu}_{4}$ is underestimated and $\hat{\mu}_{5}$ is overestimated (Figures 5a). All the other parameters always tend to converge to their real values.

Figure 6 depicts the evolution of the daily RMSE ($RMSE_{d}$) for increasing SID over the whole simulation (365 days), relative to L, N5, N6 models and to the ODNP. Such error is compared with the standard deviation of simulated power data. In view of the fact that parameter convergence behaves differently for different SIDs and in order to have a consistent comparison, the performance indices for all cases were computed starting from day $15$. Results show that the proposed approach is more sensitive to noise added to $N$, than to noise added to $P$ or $T$. The $RMSE_{d}$ computed with SID $5$ (Figure 6a) is considerably higher than the $RMSE_{d}$ computed with SID $8$ and $11$ (Figures 6b and 6c), which are really similar to each other. Clearly, under any of the tested noise conditions, the $RMSE_{d}$ achieved by all three models is by far lower than both the $RMSE_{d}$ obtained using the ODNP and the standard deviation of the simulated power.

Performance indices RMSE, MAPE, MBE, and $R^{2}$ over the whole simulation are reported and compared for each SID and each model in Figure 7. $RMSE$ and $MAPE$ slowly increase and $R^{2}$ slowly decreases with $\sigma_{N}$ (SID $2-5$) and $\sigma_{P}$ (SID $6-8$), while they are almost constant when $\sigma_{T}$ (SID $9-11$) increases. Finally, note that $MBE$ grows when $\sigma_{P}$ is increased and that it decreases both when $\sigma_{T}$ is increased and when $\sigma_{N}$ grows.

In Figure 8, the power forecasts provided by models L, N5 and N6 during a simulation run with SID $12$ are qualitatively compared to measured data during three consecutive days. The picture refers to typical examples of partially cloudy days.

8 Validation on experimental data

8.1 Experiment set up

The proposed procedure has been validated using data from a photovoltaic plant of nominal power $P_{nom}=$920\text{,}\mathrm{kW}$$p located in Sardinia (Italy). Available data consist of two datasets spanning February 2nd (33rd day of year) to May 1st (122nd day of year):

${\cal D}_{1}$: a set of data provided by the private producer running the plant, which collects hourly samples of averaged measured power and one day-ahead forecasts of air temperature;

• 2.

${\cal D}_{2}$: a set of weather reports including CCI, evaluated by a meteorological station located $20\text{\,}\mathrm{km}$ away from the plant. This data set does not have a regular sampling time.

${\cal D}_{1}$ and ${\cal D}_{2}$ have been preprocessed in order to have a single dataset ${\cal D}$ of hourly data. If there are more then one CCI values in ${\cal D}_{2}$ for the same hour, an averaged CCI is computed for that hour. On the contrary, if the CCI report in ${\cal D}_{2}$ is missing for a given hour, then the measured power and the forecast of air temperature in ${\cal D}_{1}$ for that hour are dropped222For the sake of completeness, and to further support the claim that one day-ahead temperature forecasts can be used in our procedure in place of actual data, temperature measurements collected at the meteorological station have been compared with the forecasts included in ${\cal D}_{1}$. The RMSE and MBE computed on the overall data set turned out to be equal to $1.9^{\circ}C$ and $0.9^{\circ}C$, respectively.. The data set employed is therefore given by the time series

[TABLE]

where the set of time indices $\cal K$ spans, with a sampling time $\tau_{s}=$1\text{,}\mathrm{h}$$, from the starting day $\underline{d}=33$ to the final day $\overline{d}=122$. Also in this case, only the indices $k$ corresponding to hours of light were considered. Model parameter updates are computed hourly using the data sample $D(k)$.

The initial values of the parameters have been chosen as follows 333Concerning the orientation angles $\psi$ and $\zeta$, which were not available, their values have been guessed by visually comparing the generation curve during a clear-sky day with the clear-sky irradiance model.:

[TABLE]

8.2 Performance evaluation

The performance of the proposed method on the data set introduced above has been evaluated with reference to the two standard forecast types introduced in Section 5, i.e.,

DA forecast. In this case, the forecasts and the performance indices $RMSE_{d}$, RMSE, MAPE, MBE and $R^{2}$ were computed using the same modality as in the previous section.

• 2.

HA forecast. The HA forecasts are computed according to (25)-(27) based on the data set $\cal D$ using

[TABLE]

In order to evaluate the performance, a single 7-hour HA forecast per day has been considered, and the error measures RMSE, MAPE, MBE and $R^{2}$ have been computed over the the union of all such forecasts.

The performance of the proposed method has been compared with that obtained using the following approaches.

ODNP (for DA forecasts only, clearly not suitable in the HA case).

• 2.

Autoregressive prediction models using only measurements of generated power and only cloud cover data, respectively. These methods represent a more realistic benchmark with respect to ODNP and are applicable to both DA and HA. To this end, the following two models have been tested:

• (a)

PVGM: autoregressive (AR) model of generated power with regression horizon equal to 12 hours. The parameter vector is given by $a=\left[a_{1}\ldots a_{12}\right]^{T}$ and the regressor at time $j$ is expressed by $\psi(j)=\left[P(j-1)\ldots P(j-12)\right]^{T}$. Then the prediction of generated power for time instant $j$, computed at time $k$ using parameters available at time $q\leq k$ is:

[TABLE]

where: $\hat{a}(q)$ is the estimate of $a$ available at time $q$, and for $i=1\,\ldots\,12$, the $i$-th term $P(j-i)$ of $\psi(j)$ is equal to $P^{m}(j-i)$ if $j-i\leq k$, and equal to $\hat{P}(j-i|k,q)$ otherwise.

• (b)

CCD: autoregressive with exogenous input (ARX) model of generated power with regression horizon equal to 2 hours. The inputs are cloud cover data and the clear-sky model (8)444Cloud cover index alone cannot possibly provide enough information.. The parameter vector is given by $b=\left[b_{1}\ldots b_{6}\right]^{T}$ and the regressor at time $j$ is expressed by $\xi(j)=\left[I^{0}(j)\ldots I^{0}(j-2)\leavevmode\nobreak\ \leavevmode\nobreak\ N(j)\ldots N(j-2)\right]^{T}$. Then the prediction of generated power $\hat{P}(j|k,q)$ is:

[TABLE]

where: $\hat{b}(q)$ is the estimate of $b$ available at time $q$, and for $i=0,\,1,\,2$, the $i$-th term $N(j-i)$ of $\xi(j)$ is equal to $N^{m}(j-i)$ if $j-i\leq k$, and equal to $\hat{N}(j-i|k)$ otherwise.

The parameters of the two models above have been estimated in the same recursive framework as the proposed parametric model by means of Linear Least Squares.

• 3.

ANN-based approach, using the MultiLayer Perceptron (MLP) architecture. In this case we are looking for a fair comparison in terms of model complexity, therefore two architectures with the same number of parameters as the N6 and L parametric models, respectively, have been tested using the training algorithms implemented in the MATLAB Neural Network Toolbox. Two architectures have been selected:

• (a)

MLP1: $1$ hidden unit, $6$ parameters,

• (b)

MLP2: $2$ hidden units, $11$ parameters.

Note that no such comparison is possible with the N5 model since a 5-parameter MLP cannot be synthesized. The networks have been trained in the same recursive framework as the proposed parametric model, but since a training step for the MLPs cannot be performed for each incoming data sample $D(k)$, the algorithm has been modified in order to build suitable learning sets for network training. In particular, a new learning set is constructed every $75$ samples.

8.3 Results and discussion

In Figures 9 and 10 the evolution of the parameter estimates are reported for the N5 and N6 models. All the parameters tend to converge for these as well as for the L model, not depicted to save space. The estimated parameter values and their standard deviation at the end of the experiment are shown below:

[TABLE]

[TABLE]

[TABLE]

Notice that all standard deviations are at least one order of magnitude lower than the respective parameter estimates. The only exception concerns the parameter estimate $\hat{\mu}_{3}$ for the N5 model. This model, which is characterized by a minimal number of parameters, apparently does not fully capture the dependency between the generated power and temperature, thus generating persistent parameter drifts. On the contrary, the slightly overparameterized N6 model shows much better convergence properties, as it is also apparent from the comparison of Figures 9 and 10. Furthermore, the L model shows higher parameter standard deviations, on average, with respect to N6.

In order to analyze the impact of overparametrization on the consistency of the parameter estimates with the physical model (2), we compare the values of $\hat{\mu}_{1}$ and the correction terms $\hat{\eta}_{2}=\hat{\mu}_{2}/\hat{\mu}_{1}$ and $\hat{\eta}_{3}=\hat{\mu}_{3}/\hat{\mu}_{1}$ obtained from (42), (47) and (56). The estimates of $\hat{\mu}_{1}$ provided by models N5 and N6 are slightly higher than the nominal power/irradiance gain $P_{nom}/1000=0.92$, while the value obtained using model L is remarkably higher. In Table 4, $\hat{\eta}_{2}$ and $\hat{\eta}_{3}$ are compared with their typical values reported in (3). Notice that N6 is the only model for which $\eta_{2}\in{\cal S}_{2}$ and $\eta_{3}\in{\cal S}_{3}$, while for N5 $\hat{\eta}_{3}$ is higher than the upper bound of ${\cal S}_{3}$, which apparently means that $\mu_{3}$ is overestimated. On the contrary, model L tends to underestimate $\mu_{3}$.

Notice that in model N6 the overparametrization is given only by the introduction of parameter $\mu_{6}$, defined as the product between $\mu_{2}$ and $\mu_{4}$. From (47), we have that $\hat{\mu}_{2}=$-1.60\text{\times}{10}^{-4}, $\hat{\mu}_{4}=-1.19$, and $\hat{\mu}_{6}=$-2.53\text{\times}{10}^{-4}. Therefore $\hat{\mu}_{2}\hat{\mu}_{4}=$-9.12\text{\times}{10}^{-4}$$, which is the same order of magnitude as $\hat{\mu}_{6}$, thus suggesting that the estimates are almost consistent. A similar consistency is not shown by model L, despite its good performance on the DA and HA forecasts, as shown in the sequel.

In order to further evaluate the consistency of models N5 and N6, we analyze the CCF curve $C(N)$ in (5) obtained using the estimated values of $\mu_{4}$ and $\mu_{5}$ provided by such models. In Figure 11 we report the plots of the $C(N)$ using the respective values of $\hat{\mu}_{4}$ and $\hat{\mu}_{5}$. It turns out that both the estimated CCFs are consistent with the seasonal trend for the mediterranean climate [34].

As far as the forecasting performance evaluation is concerned, all error measures were computed over the period starting from day $d_{0}=57$ in order to guarantee, as in the simulated data case, at least a rough adaptation of the model parameters.

Table 5 summarizes the performance indices achieved by the proposed parametric models, compared with MLPs, ODNP, PVGM and CCD on both DA and HA forecasts. Since the performance of the MLPs may be quite sensitive to the initial conditions, 30 simulations of MLP1 and MLP2 with the same data set have been performed with suitably generated random initial conditions555MLP layer weights and biases are initialized according to the Nguyen-Widrow algorithm. This algorithm chooses values in order to distribute the active region of each neuron in the layer randomly but evenly across the layer’s input space., and mean value/variance of all performance indices over all simulations have been reported. Clearly, all parametric models and MLPs perform largely better than the ODNP and PVGM. On average, MLPs show errors 20% higher than parametric models. The CCD model yields intermediate indices between our models and MLPs. The performances of the three proposed parametric models are overall comparable. Notice that in network operation, performance indices qualifying forecast accuracy, e.g., RMSE and MAPE, are generally normalized with respect to the nominal plant power. This normalization is useful in view of how the DSO quantifies uncertainty in power generation. Actually, when dealing with hundreds or thousands of distributed units, either for maintenance activities or operational requirements, the DSO takes into account how much each specific unit is expected to produce with respect to its nominal power. Indeed, network operation requires the DSO to know at any time the overall power uncertainty pertaining to the whole distributed generation. Such uncertainty can be easily estimated from the nominal capacity of each plant (known a-priori) and the normalized plant forecast error (computed off-line on historical data and usually with a smaller time resolution, typically, weekly or monthly). For this reason, we consider the normalized errors $RMSE_{NP}$ and $MAPE_{NP}$. The results obtained for the latter indices, which are around 8-10%, can be considered acceptable for most problems in network operation.

Figure 12 shows the daily RMSE ($RMSE_{d}$) achieved by the proposed models in the DA forecast, compared with ODNP, PVGM, and CCD. In Figures 13, such error is compared with those achieved by the MLPs in the DA forecast. Figure 14 depicts the comparison of the $RMSE_{d}$ achieved by the various models and benchmarks in the HA forecast. In the above figures, all errors are also compared with the standard deviation of measured power data.

Both in Table 5 and Figure 14, data for the ODNP in the HA case are not reported as the ODNP is unsuitable as a reference model.

The apparently high normalized forecast errors (in particular the MAPE) are mainly due to the fact that the CCI data refer to a weather station located about $20\text{\,}\mathrm{km}$ away from the plant site. Another source of forecast dispersion is the rough CCI data resolution. Both these two sources of uncertainty deteriorate forecasting accuracy as much as the cloudiness at the plant site changes rapidly in time. It is worth remarking, however, that the proposed models use a very limited amount of information which is easily and cheaply available, and that such raw data is informative enough for a reliable estimation of the generation model, as demonstrated by the parameter convergence properties highlighted in Figure 10. Furthermore, it is expected that the proposed approach can be fruitfully adapted to large-scale aggregation of plants over macro areas, yielding significant performance improvements with respect to the single plant case [36], with no additional data requirements (except generation measurements).

In Figure 15 the DA forecasts provided by N6 model and MLP1 during three different days and under three different weather conditions is compared with the measures of generated power.

The procedure has also been tested using the last two weeks of data for validation purposes only, that is, model parameter values were frozen two weeks before the end of the data set, then offline DA and HA forecasts for the whole two following weeks were generated according to the modalities above, and the respective performance indices were computed. It turned out that the performances of the various models in terms of all indices were almost identical to those obtained with online parameter adaptation (the difference in terms of all error measures was below 1%). In this respect, it is worth remarking that online parameter estimation has the advantage of capturing seasonal parameter variations (especially as far as $C(N)$ is concerned).

As far as the computational time is concerned, the MLP update algorithms take generally longer to execute with respect to the parameter update procedure. More specifically, a MLP training step using a $K$-sized sample data set requires longer computational time than $K$ model parameter updates using a single data sample. The overall computational times for all models and MLPs are reported in Table 6. All algorithms were implemented in Scilab [37] version 5.5.2 and executed on a 2.4 Ghz Intel Xeon(R) v3 processor running the Linux operating system.

Summarizing the discussion above, we found that the proposed parametric models show comparable forecasting behavior and a general performance improvement with respect to MLPs, both from an accuracy and a computational viewpoint. In particular, model N6 exhibits the best compromise between performance, complexity and convergence/consistency properties. This apparently shows that a slight overparameterization with respect to a minimal model inspired by the theoretical models of CCF and power generation, turns out to be beneficial. On the contrary, excessive overparameterization impacts model identifiability and consistency.

9 Conclusions

In this paper, we have proposed an efficient parametric technique aimed at model estimation for direct forecasting of PV power generation using cloud cover data. The approach is based on recursive least squares and Extended Kalman Filter. To this purpose, three models have been singled out: an 11-parameter linear model, and two nonlinear models, involving 5 and 6 parameters, respectively. The procedure is especially fit for the typical scenario where the network operator, due to the large number of managed producers, has no access to on-site irradiance and temperature measurements. The method exploits only the historical time series of generated power, cloud cover, and forecast temperature, which can be obtained from a meteorological service. The procedure involves modest memory and computational requirements. Its performance has been evaluated on simulated data as well as on a single plant located in Italy. The experimental results show a substantial improvement in terms of typical error measures with respect to alternative approaches based on Neural Networks and standard linear autoregressive models. All the proposed models show equivalent forecasting performance, while the 6-parameter nonlinear model generally outperforms the other two as far as data fitting properties are concerned. It is expected that the proposed approach can be fruitfully adapted to the aggregation of several plants over macro areas, with significant performance improvements with respect to the single plant case. This issue is the subject of current investigation.