Model Estimation for Solar Generation Forecasting using Cloud Cover Data

Mo del Estima tion for Solar Generati on F orecasti ng using Clo ud Cov er Da ta Daniele Pepe , Gianni Bianchini , A ntonio Vicino Dip artimento di Inge gneria del l’Informazione e Scienze Matematiche, Universit` a di Siena, Via R oma 56, 53100 Siena, Italy Abstract This paper presen ts a parametric mo del appro a c h to addre s s the problem o f photov oltaic gener ation forecasting in a scenario where mea s uremen ts of meteorolo gical v a riables, i.e., s olar irradia nc e a nd temper ature, ar e not av ailable at the plant site. This s cenario is relev a nt to electric ity netw ork op eration, when a large num b er o f PV plan ts are deployed in the grid. The propo sed metho d makes use of r a w cloud cover da ta provided by a meteorolog ical ser vice co m bined with p ow er genera tion measurements, a nd is particularly suitable in PV plant in tegra tion on a large-scale bas is, due to low mo del co mplexit y and computationa l efficiency . An extensive v alida tion is p erformed using b oth simulated and rea l data. Keywor ds: Photov oltaic gener ation, Mo deling, Estimation, F oreca s ting, Cloud cover. 1. In tro duction A ma jor challenge in the integration of renewable energy source s into the grid is that p ow er generation is intermitten t, difficult to control, and strongly dep ending on the v ar iation of weather conditions. F or these reaso ns, forecasting of r enew able dis tributed g eneration has b ecome a funda- men tal need to gr id o p erato r s. In this resp ect, solar genera tion forecasts o n multiple time horizo ns are needed to satisfy grid co nstraints and demand. In particular , shor t-term fore c asts a r e req uired for the purp oses of p o wer pla n t op eration, grid balancing, rea l- time unit dispatching, automatic generation control, and energy trading. On the contrary , longer -term fore c asts ar e of interest to Distribution System Op erators (DSO) and T ransmission System Op erator s (TSO) for unit com- Email addr esses: pepe@di ism.unisi.it (Daniele P ep e), gian nibi@diism.u nisi.it (Gianni Bianchini), vicino@d iism.unisi.i t (Ant onio Vicino) Pr eprint submitte d to Solar Ener gy January 23, 2019 mitmen t, scheduling and for improving ba lance area control (see, e.g ., [1, 2, 3, 4] and refer ences therein). Concerning sola r p ow er generatio n, muc h attention has b een paid to the problem o f obtaining accurate day-ahead a nd hour -ahead for ecasts o f solar irradiance and/o r gener ated p ow er (see [5 , 6] for a co mprehensive overview on the s ub ject). Most contributions focus on solar irra diance forecasting [7]. Widely adopted approa c hes are based on Artificia l Neural Netw orks (ANNs) [8, 9, 10, 1 1] or Supp ort V ector Machines [12] with different t yp es of input data. Alternatively , clas sical linear time series foreca sting metho ds are used in [1 3, 14, 15], where the c onsidered time ser ie s is t ypically the global horizontal irr adiance (GHI) norma lized with a clear- sky mo del (see [1 6, 17] for a compr ehensiv e review). Glo bal radiation for ecasts are then fed a long with temp erature foreca sts to a simulation mo del of the plant [18] to compute the prediction of p ow er generation. In any case, co mputing re lia ble gener ation forecasts from pr e dicted meteor ological v ar iables hinges up on the av ailability of a reliable mo del of the pla nt, b e it physical or e s timated from da ta. Unfortunately , in man y practical scena rios, neither reliable plant mo dels, nor direct on-site mea- surements of solar irr adiance and other meteorolo gical v aria bles such as temp erature a re av ailable. This is the t y pical cas e of a a DSO dealing with hundreds o r thousands of distributed heterogeneo us and indep endently-oper a ted s o lar plants; in this case, the o nly av aila ble pla nt data is represe nted by generated pow er measuremen ts pr o vided by electronic meters. The con tribution of this pape r fo- cuses on the pr oblem of estimating reliable gener ation models in s uc h limited infor mation con texts. In [19, 2 0], tw o metho ds a re pro posed to estimate the parameters o f the PVUSA mo del [2 1] of a PV pla n t via a recur s iv e framework based only on measures o f genera ted p ow e r and temp erature forecasts. This solution is not alwa ys data- efficien t, s inc e g eneration data collec ted during cloudy days a re not explo ited. In order to obtain mor e ac curate re s ults, p o wer meas ur emen ts can b e combined with further data coming from a weather serv ice. Suc h data are typically av era ged over la rge geog raphic z o nes and therefore they ma y b e s c arcely informa tive for a sp e cific sp ot, yet they may pr o vide useful information when the goa l is to address the ag gregation of mult iple plants ov er a macro area . In this resp ect, clo ud cov er measurements der iv ed from satellite imag ing are an exceptio n, since they can b e made av ailable for s pecific lo cations with go o d spatia l and temp oral res o lution (up to 2 . 5 km × 2 . 5 km and 30 min, resp ectively) [22]. Ba sed on such data, suitable mo dels can b e estimated for irradia nce forec asting purp oses (see [23, 24, 2 5] for details). Irradia nce forecasts ca n 2 be used as inputs to a plan t g eneration mo del in o rder to provide energy productio n for ecasts. Of co urse, the latter step r equires that a reliable mo del of p o wer generatio n from ir radiance b e previously estima ted. In this pap er, we present a nov el a pproach for dir ect foreca sting of PV plant p ow er g eneration from cloudines s data. T o this purp ose, a clas s of para metric mo dels is intro duced which efficiently exploits the P VUSA mo del [2 1] and the notion o f Cloud Cover F actor (CCF) [26, 27] in a limited information scenar io. More sp ecifically , the mo del par ameters are estimated using generated p o wer data co m bined with additional informa tion provided b y a meteorolo gical service for the a rea where the pla nt is lo cated. Such da ta consist o f a time ser ies of raw cloud c over and temp erature re p orts. F or the prop osed mo dels, estima tion pr ocedur e s based on Recursive Least Squa res (RLS) a nd the Extended K alman Filter (E KF) [28], are devised. The pro perties and the p erformance of the prop osed metho d are demonstrated b oth in a simulated scenario and on ex perimental data from a plant curr e n tly in op eration. Pr eliminary results lea ding to this pap er were presented in [29]. The ma n uscript is or ganized a s follows. In Section 2 we int ro duce the mo delling to ols used; in Section 3 the prop osed mo dels a re pres en ted. Estimation pro cedures are presented in Sectio n 4 . Section 5 addre s ses the r elev ant for ecasting pro blems, while p erformance ev aluation is discussed in Section 6. Sim ulatio n results are illus tr ated in Section 7. Exp e rimen tal v a lidation re s ults a nd related disc us sion are rep orted in Sectio n 8. Finally , conclusions ar e drawn in Section 9. 2. Mo dels and metho ds 2.1. The PVUSA photovoltaic plant mo del A P V plant can b e efficiently mo delled using the PVUSA mo del [21], which expresses the instantaneous generated p ower a s a function of irradiance a nd air temper ature a ccording to the equation: P = µ 1 I + µ 2 I 2 + µ 3 I T , (1) where P , I , a nd T are the genera ted p o wer (kW), irradia nce (W / m 2 ), and air temp erature ( ◦ C), resp ectiv ely , and µ 1 , µ 2 , µ 3 , are the mo del para meters. Mo del (1) is linear and parsimonio us in terms of num ber of par ameters. Despite its simplicity , very go o d a ccuracy is obtained when this mo del is fit to real measured data. Mo reov er , due to their reliability , temper ature forecasts tak en from a meteoro logical ser vice can b e used efficien tly in place of actual measur emen ts in order to estimate the mo del parameters , as shown in [19]. This cannot 3 be do ne for irradia nc e , a s forecasts a r e muc h less reliable. W e find it co n venien t to expr e ss (1) also in the form P = µ 1 (1 + η 2 I + η 3 T ) I , (2) where µ 1 plays the role of the main power/irradiance ga in o f the plant, while η 2 = µ 2 /µ 1 and η 3 = µ 3 /µ 1 int ro duce corr ection terms. It is worth noticing that η 2 and η 3 in (2) ar e characterized by a na rrow ra nge of v aria bilit y a mong different PV technologies. Indeed, t ypical v alues o f these parameters a re g iv en by [21]: η 2 ∈ S 2 =  − 2 . 5 × 10 − 4 , − 1 . 9 × 10 − 5  , η 3 ∈ S 3 =  − 4 . 8 × 10 − 3 , − 1 . 7 × 10 − 3  . (3) 2.2. Cloud Cover In dex and Cloud Cover F actor The Clo ud Cover Index (CCI), here denoted b y N , refer s to the frac tion of the s ky obscur ed by clouds when o bserved from a particular lo cation. In meteor ology , N is mea s ured in okta [30], i.e., a n int eger ranging from 0 ok ta (clear sky) up to 8 o k ta (completely overcast). W eather s ervices often provide this informa tion in terms of p ercent age of cov ered s ky: from 0 . 0 (clear sky) to 1 . 0 (completely ov ercast) in steps o f 0 . 1 . A widely inv estigated issue is the relationship b et ween the CCI and the glo bal irr adiance o n a surface. In [27], a simple w ay of estimating solar ra dia tion in cloudy days is to m ultiply the clear-sky irradiance by the s o-called Cloud Cover F acto r (CCF), which is defined as a function of the CCI. Let I 0 be the clea r-sky irradia nce on a given surface. It turns o ut that the global irradia nce on the same sur face is given by I ( N ) = C ( N ) · I 0 , (4) where C ( N ) is the CCF, which can b e ex pressed in the form C ( N ) = 1 + µ 4 N + µ 5 N 2 . (5) In (5), the parameters µ 4 and µ 5 depe nd on the lo cal clima te. In particular , µ 5 is always negativ e. Depending on the sign of µ 4 , the curve C ( N ) may be non-monotonic (see Fig. 1). Indeed, in some lo cations C ( N ) may ris e ab ove 1 (i.e., the clear-s k y v alue) under slightly clo udy sk y (e.g , 0 < N < 0 . 3), since the diffuse ra diation b ecomes significant with middle-lo w cloudiness. Mo reov er , while the CCF is indepe ndent of solar elev ation, it is sub ject to seaso nal v aria tions. Alternativ e CCF mo dels ex ist in the liter ature [2 6]. 4 0 1 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 0.3 0.5 0.7 0.9 0 1 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 0.3 0.5 0.7 0.9 1.1 P S f r a g r e p la c e m e n t s N N C ( N ) C ( N ) Figure 1: Qualitative b eha vior of the CCF C ( N ) for µ 4 > 0 (left) and µ 4 < 0 (righ t). The CCF is usually estima ted for a g iven lo cation using on-site irradia nce mea suremen ts [3 1]. It is worth to remark, ho wev er, that the metho d pr esen ted in this pap er do es not rely on such measurements. 2.3. Cle ar-sky irr adianc e mo del In this pap er we need to explo it an estimate of the clea r-sky ir r adiance on a g iv en surface, a s opp osed to the GHI. As a first step, we mak e use o f the Helio don sim ulator mo del [3 2]. This mo del allows to compute the theoretical clear- s ky normal irradianc e (W / m 2 ) from the solar altitude h , i.e., the ang le over the ho rizon (rads), as: I cs,n =    A · 0 . 7 ( 1 sin h ) 0 . 678 if 0 < h < π / 2 0 otherwise, (6) where A = 135 3 W / m 2 denotes the appa ren t extra terrestrial irra dia nce. Given the theore tica l clear-sk y normal irradia nce I cs,n , the clea r-sky irra dia nce on an inclined surface I cs can b e derived from I cs,n and s ur face orientation with resp ect to the sun p osition. Denoting b y ζ the surfa c e azimuth a nd ψ the surfa c e tilt angle, o ne ha s that I cs = [sin( ψ ) cos( h ) cos( ζ − γ ) + cos( ψ ) sin( h )] I cs,n , (7) where γ is the solar azimuth. Clearly , I cs can b e co mputed for given v alues of ζ a nd ψ from la titude, longitude and time of day . F or ψ = 0, the irradia nce on an hor izon tal surface is o bta ined. In this study , we mak e the realistic assumption that the e x act o rien tation o f the PV panel surfaces of the c onsidered plant is not known a-prior i. Ho wev er, it is reasona ble to assume that the plant is efficient ly oriented for the sp ecific latitude acc ording to, e.g ., the guidelines given in [33]. 5 Therefore, we will use as the reference v a lue of the theoretical clea r-sky ir r adiance I 0 for a sp ecific plant, the v alue of (7) for ( ζ , ψ ) = ( ζ 0 , ψ 0 ), where ζ 0 and ψ 0 are ta ken fr om the ab ov e guidelines, i.e., I 0 = [sin( ψ 0 ) co s( h ) cos( ζ 0 − γ ) + cos( ψ 0 ) sin( h )] I cs,n . (8) The choice of the mo del in (6) is mainly due to its simplicity . More accurate mo dels ta king int o account atmo s pheric par ameters such a s aer osol o ptical depth, ozo ne, a nd water v ap or, could b e used up on av ailability o f such data in order to improv e the p erformance of the forec asting metho d prop osed in this work. 3. Plan t generation mo del The plant mo del (1) a nd the irr a diance mo del (4) − (5) introduce d in the pre v ious sectio n can be combined in order to obta in an expressio n of the genera ted p ow e r as a function of N , T , a nd I 0 , that is P = P ( I 0 , T , N ) = ( µ 1 + µ 2 I ( N ) + µ 3 T ) I ( N ) , (9) where I ( N ) =  1 + µ 4 N + µ 5 N 2  I 0 . (10) Therefore, the ov erall mo del is defined by the pa r ameter vector µ = h µ 1 µ 2 µ 3 µ 4 µ 5 i T , and (9) − (10) ca n b e rewr itten in a standard r e g ression form a s P = ϕ T ( I 0 , T , N ) θ ( µ ) , (11) where θ ( µ ) and ϕ ( I 0 , T , N ) are given by θ ( µ ) = h θ 1 ( µ ) θ 2 ( µ ) θ 3 ( µ ) θ 4 ( µ ) θ 5 ( µ ) θ 6 ( µ ) θ 7 ( µ ) θ 8 ( µ ) θ 9 ( µ ) θ 10 ( µ ) θ 11 ( µ ) i T = h µ 1 µ 1 µ 4 µ 1 µ 5 µ 2 2 µ 2 µ 4 µ 2 µ 2 4 + 2 µ 2 µ 5 2 µ 2 µ 4 µ 5 µ 2 µ 2 5 µ 3 µ 3 µ 4 µ 3 µ 5 i T , (12) ϕ ( I 0 , T , N ) = h I 0 I 0 N I 0 N 2 I 0 2 I 0 2 N I 0 2 N 2 I 0 2 N 3 I 0 2 N 4 T I 0 T I 0 N T I 0 N 2 i T . (13) 6 Notice that the mo del in (11)-(13) is not linear in the para meter vector µ . In order to estimate mo del para meters, tw o appr oaches will b e pursued in the following by comparing their resp ective per formances: (N) Kee p the parameter v ector size to a minimum and derive a nonlinear estimation algorithm based on the Extended Ka lman Filter (EKF); (L) Introduce an ov erpar ameterization by disreg arding the pa rameterization o f θ ( µ ) in µ a nd assuming θ ( · ) to b e a vector of independent parameters, thu s making the mo del linea r and amenable to standard least squares es tima tio n. The for mer approa c h has the a dv ant age of keeping the mo del par simonious and improving its ident ifiability features, while the latter is introduced with the purp ose of mak ing the estimatio n task simpler . 3.1. Nonline ar (N) mo dels The fir st idea prop osed here is to tackle the nonlinear estimation proble m explo iting the EKF. T o this pur pose, let us co nsider mo del (11) − (13) for the following tw o choices of the indep enden t parameter vector µ : 1. N5 mo del µ = h µ 1 µ 2 µ 3 µ 4 µ 5 i T , (14) θ ( µ ) = h µ 1 µ 1 µ 4 µ 1 µ 5 µ 2 2 µ 2 µ 4 µ 2 µ 2 4 + 2 µ 2 µ 5 2 µ 2 µ 4 µ 5 µ 2 µ 2 5 µ 3 µ 3 µ 4 µ 3 µ 5 i T , (15) which is the or ig inal nonlinea r mo del (11) − (13) characterized b y a minimal num b er of pa- rameters, and 2. N6 mo del µ = h µ 1 µ 2 µ 3 µ 4 µ 5 µ 6 i T , (16) θ ( µ ) = h µ 1 µ 1 µ 4 µ 1 µ 5 µ 2 2 µ 6 µ 4 µ 6 + 2 µ 2 µ 5 2 µ 5 µ 6 µ 2 µ 2 5 µ 3 µ 3 µ 4 µ 3 µ 5 i T , (17) which implies a slight overparameterization with resp ect to mo del N5, amounting to treating µ 2 , µ 4 and µ 2 µ 4 = µ 6 as indep endent para meter s. F or either choice, the mo del ta kes the fo rm (11) P = ϕ T ( I 0 , T , N ) θ ( µ ) . 7 3.2. Line ar overp ar ameterize d (L) mo del This mo del is derived b y trea ting θ ( µ ) in (12) as a vector θ of 11 independent par ameters. This wa y mo del (11) b ecomes P = ϕ T ( I 0 , T , N ) θ and is clear ly linear in the parameters . F or this mo del, a standa rd Recurs iv e Leas t Squa res (RLS) algorithm can b e us ed to es timate θ . 4. Mo del estim ation F or the models in tro duced in Sections 3.1 and 3.2, an o n-line up date of the estimate of the parameter vector µ (or θ ) is p erformed. T o the purp ose of illustrating the pro cedure, let us define the following qua ntities: • τ s : sa mpling time (mins); • k : discrete time index; • d : g eneric day; • τ d ( k ): time of day (TOD) corr esponding to index k (at the given long itude); • I 0 ( k ): theoretical clea r-sky sola r irradiance at time index k at the plant site, computed according to (6)-(8); • K d = { k d , . . . , k d } : set o f time indices p ertaining to light hour s in day d , i.e., the theor etical irradiance (6)-(8) is nonzero for τ d ( k ) corr esponding to k d ≤ k ≤ k d , and moreover k d +1 = k d + 1 , that is, only daylight time indices ar e consider ed; • ˆ µ ( k ) [ ˆ θ ( k )]: estima te o f the parameter vector at time k , after pro cessing the last data sa mple, being ˆ µ (0) [ ˆ θ (0 )] the initial g uess; • W ( k ): weather re p ort for time k relative to the the plant area provided by a meteorolog ical service. More specific a lly , W ( k ) = { N ( k ) , T ( k ) } , where N ( k ) and T ( k ) are the cloud c over and temp erature r epor ts , resp ectively 1 . 1 T ( k ) can be represente d by either a real-time estimate of the plan t site temperature or even some forecast 8 • P m ( k ): meas ured generated p ow er at time k ; • D ( k ) = { P m ( k ) , W ( k ) } = { P m ( k ) , N ( k ) , T ( k ) } : data av ailable at time k ; • φ ( k ) = ϕ ( I 0 ( k ) , T ( k ) , N ( k )): regression vector at time k (see (13)); The estima tio n pr o cedure is recursive. At ea c h time s tep k , a new p o wer-weather da ta s ample { P m ( k ) , W ( k ) } is acquired, the algo rithm is run and a curren t estimate ˆ µ ( k ) [ ˆ θ ( k )] of the parameter vector is computed acco rding to the criteria outlined b elow. • N5 and N 6 models . F or the purp ose of estimating µ we cons ider the standar d extended Kalman filter (EKF) formulation defined by the following state-space mo del [28]: µ ( k + 1) = µ ( k ) (18a) P ( k ) = φ T ( k ) θ ( µ ( k )) + w ( k ) , (18b) where µ ( k ) repr esen ts the state vector, P ( k ) is the output, and w ( k ) is a zero -mean Gaus- sian white no ise with given v ar iance r whic h describ es b o th the measur emen t nois e and the mo del uncer tain ty . According to the standard EK F theory , the following up date law for the parameter estima te ˆ µ ( k ) is obtained: ˆ µ ( k ) = ˆ µ ( k − 1) + K ( k )  P m ( k ) − φ T ( k ) θ ( ˆ µ ( k − 1))  , (19) where K ( k ) = R ( k − 1) H T ( k ) H ( k ) R ( k − 1) H T ( k ) + r , (20) and the estimate cov ariance matrix R ( k ) is up dated at each step accor ding to the r ule R ( k ) =  I − H T ( k ) K T ( k )  R ( k − 1) , (21) being H ( k ) = ∂ ∂ µ  φ T ( k ) θ ( µ )  µ = ˆ µ ( k − 1) . thereof, computed in adv ance. I ndeed, it has been s ho wn [19] that the reliabil it y of day-ahead temperature forecasts is comparable to that of actual measurements for model estimation purp oses. Obv iously , when on-si te measuremen ts of temperature T m ( k ) are a v ailable, then T ( k ) should be tak en equal to T m ( k ). In the exp erimen tal section of this paper, data from day -ahead forecasts are considered. 9 As in common pra ctice, the initial estimate cov aria nc e matrix R (0) is computed as R (0) = l (0) · I . The hig her the confidence on the initial guess ˆ µ (0), the low er l (0) can b e chosen. The initialization o f the para meter vector estima te ˆ µ (0) is discus s ed a t the end of this sec tio n. • L mo del . A re c ur siv e least-s quares (RLS) estima tio n s tep is p erformed, i.e., we compute L ( k ) = V ( k ) φ ( k ) , (22) where V ( k ) is the weigh t matrix a t time k , w hich in turn is r ecursively co mputed according to the standard RLS up date law: V ( k ) = V ( k − 1) − V ( k − 1) φ ( k ) φ T ( k ) V ( k − 1) 1 + φ T ( k ) V ( k − 1) φ ( k ) . (23) Finally the current parameter estimate is up dated a ccording to the law ˆ θ ( k ) = ˆ θ ( k − 1) + L ( k ) h P m ( k ) − φ T ( k ) ˆ θ ( k − 1) i . (24) As in the usual RLS practice, the initial weight matrix V (0) is s e t to V (0) = l (0) · I , whe r e I is the identit y ma trix and l (0) > 0 is chosen according to the co nfidence given to the initial parameter g ue s s ˆ θ (0) (higher l (0) meaning less co nfidence). As far a s the initialization of the parameter vector is co ncerned, a r eliable guess for the initia l v alues ˆ µ (0) in mo del N5 can b e computed fr om the following guidelines . • A go o d guess for the plant gain µ 1 is given b y ˆ µ 1 (0) = P nom / 1000 where P nom is the nominal plant p ow er in kW [19]; • the initial v alues ˆ µ 2 (0) a nd ˆ µ 3 (0) ca n b e chosen such that η 2 (0) = ˆ µ 2 (0) / ˆ µ 1 (0) and η 3 (0) = ˆ µ 3 (0) / ˆ µ 1 (0) are the central p oints o f the resp ective interv als in (3); • ˆ µ 4 (0) a nd ˆ µ 5 (0) may b e selected as the av erage v alues of the corr espo nding pa rameters for the climate of the mac r o-area where the plant is lo cated, whic h ar e usually av a ilable. F or example, for the Italia n regions considered in the experimental sec tion of this pap er such v alues are estimated in [34]. Given their r espective definitions, for the N6 mo del it can b e assumed that ˆ µ 6 (0) = ˆ µ 2 (0) ˆ µ 4 (0), 10 while for the L mo del a natur al choice (se e (12)) is g iv en by θ ( µ ) = h ˆ µ 1 (0) ˆ µ 1 (0) ˆ µ 4 (0) ˆ µ 1 (0) ˆ µ 5 (0) ˆ µ 2 (0) 2 ˆ µ 2 (0) ˆ µ 4 (0) ˆ µ 2 (0) ˆ µ 4 (0) 2 + 2 ˆ µ 2 (0) ˆ µ 5 (0) 2 ˆ µ 2 (0) ˆ µ 4 (0) ˆ µ 5 (0) ˆ µ 2 (0) ˆ µ 5 (0) 2 ˆ µ 3 (0) ˆ µ 3 (0) ˆ µ 4 (0) ˆ µ 3 (0) ˆ µ 5 (0) i T . 5. F orecasting In this section, the p erformance of mo dels N5, N6 and L , prop osed in Section 3 is ev a luated on the widely used Day-Ahead (DA) and Hour-Ahead (HA) forecas ts [35]. In order to suitably define them in o ur co n text, let d and k b e the generic day and time insta n t, resp ectively , in which a forecast is suppo sed to b e computed and submitted. F or a given time instant j ≥ k , let ˆ W ( j | k ) = { ˆ N ( j | k ) , ˆ T ( j | k ) } denote the weather forecast { cloud cov er, tempera ture } relative to time j av ailable at time k . W e denote by ˆ P ( j | k ; q ) the prediction of ge nerated p o wer for time instant j , computed a t time k using the para meter vector estimate av ailable at time q ≤ k according to the appropria te mo del, i.e., ˆ P ( j | k ; q ) =    ϕ T ( I 0 ( j ) , ˆ T ( j | k ) , ˆ N ( j | k )) · θ ( ˆ µ ( q )) (N5 , N6) ϕ T ( I 0 ( j ) , ˆ T ( j | k ) , ˆ N ( j | k )) · ˆ θ ( q ) (L) . (25) 5.1. Day-Ahe ad for e c ast The day-ahead forecast is us ua lly s ubmitted at 6 am on the day b efore ea c h op erating day , which b egins a t midnight on the day of submission, and cov ers all 24 hour s of that op erating day . Therefore, the day ahea d forecast is provided 19 to 4 2 hours prior to the op erating time. In this resp ect, we r ecall that the v ast ma jority of conv entional generation is scheduled in the D A market. F or a given day d , the DA forecast is submitted at time instant k d 0 corres p onding, e.g., to 6 am and consists o f the time series given by ˆ P DA ( d ) = n ˆ P ( j | k d 0 ; k d − 1 ) , j ∈ K d +1 o . (26) Such a forecast is computed on the bas is of the la st parameter es timate ˆ µ ( k d − 1 ) [ ˆ θ ( k d − 1 )] av ailable the day b efore (see Figur e 2 for the mea ning of k d 0 and the forecas t ho r izon). 11 P S f r a g r e p la c e m e n t s j K d +1 d + 1 d d − 1 ˆ P DA ( d ) 06:00 k d 0 (present) k d − 1 Figure 2: Da y-ahead forecasting. Note that K d +1 does not co ver a complete da y p erio d since k d +1 refers to the last sample pertaining to light hours in day d + 1. 5.2. Hour-Ahe ad for e c ast The hour-a head forec ast is usually submitted 105 min utes prio r to ea c h op erating hour and provides an a dvisory forecas t for the 7 hour s o f light (or the remaining ones, if les s) of the s a me day after the op erating hour. Therefore, such a forecast, a ssuming it is computed at the time index k just b efore submission, is g iv en by the time ser ie s ˆ P H A ( k ) = n ˆ P ( j | k ; k ) , j ∈ K H A ( k ) o (27) where K H A ( k ) deno tes the set of time indices p ertaining to the 7 -or-less-ho ur horizon star ting at the b eginning of the r elev ant op erating hour (see Figure 3 ). Note that ˆ µ ( k ) [ ˆ θ ( k )] represe nts the most rece nt parameter estimate av a ilable. 6. P erformance ev aluation In this section w e in tro duce the p erformance a ssessment indices that will be used to ev a lua te the fitness of mo dels to data a nd the efficac y of the pro posed metho ds in the for e casting pr oblems int ro duced in the previous sec tio n. It is worth to r emark that a single estimated mo del is used to provide b oth DA a nd HA foreca sts. 12 P S f r a g r e p la c e m e n t s j T O D K H A ( k ) d + 1 d d − 1 ˆ P H A ( k ) 0 6 : 0 0 k k d − 1 Figure 3: Hour- ahead forecasting 6.1. Err or me asur es F or the sake of simplicity , a g eneric definition o f the perfor mance indices tha t will b e used is given here. Details o n how such indices are computed for sp ecific problems (e.g., DA or HA forecasting) will b e pr o vided in the sequel. Let us consider a da ta set P P = n P m ( j ) , ˆ P ( j )  j ∈ K o (28) where K is a set of time indices of cardinality K . The set P is comp osed of pairs  P m ( j ) , ˆ P ( j )  with P m ( j ) > 0 , ˆ P ( j ) > 0, where ˆ P ( j ) represents the forecas ted p ow e r and P m ( j ) the corresp onding measured v alue. W e co nsider the fo llo wing s ta ndard er ror measures p ertaining to P : • Ro ot-Mean-Squar e Erro r (RMSE) RM S E = s 1 K X j ∈K  P m ( j ) − ˆ P ( j )  2 (29) • Mean-Bia s -Error (MBE) M B E = 1 K X j ∈K  P m ( j ) − ˆ P ( j )  (30) 13 • Mean-Absolute-Percentage-Err or (MAPE) M AP E = 1 K X j ∈K      P m ( j ) − ˆ P ( j ) P m ( j )      · 100 (31) • Determination c oefficient R 2 = 1 − P j ∈K  P m ( j ) − ˆ P ( j )  2 P j ∈K  P m ( j ) − P  2 (32) where P = 1 K P j ∈K P m ( j ) is the mean o f the measur ed power ov er the data set. Notice that here it is implicitly a ssumed that P j ∈K  P m ( j ) − ˆ P ( j )  2 ≤ P j ∈K  P m ( j ) − P  2 and in this case 1 − R 2 is the so-ca lled unex plained square er r or. • Normalized RMSE (NRMSE) N RM S E = v u u u t P j ∈K  P m ( j ) − ˆ P ( j )  2 P j ∈K  P m ( j ) − P  2 = p 1 − R 2 . (33) In a ddition to the ab o ve standa rd statistica l indice s , we also consider the following tw o er ror mea- sures, which are o f practica l interest for netw ork op eration, a s they are referred to the nominal plant p ow er P nom : • Normalized RMSE w.r .t. nominal p ow er ( RM S E N P ) RM S E N P = RM S E P nom (34) • Normalized MAPE w.r.t. nominal p ow er ( M AP E N P ) M AP E N P = 1 K X j ∈K      P m ( j ) − ˆ P ( j ) P nom      · 100 . (35) 6.2. One-day-ahe ad n aive pr e dictor As an additional ev aluation to ol, the per formance indices achiev ed using the prop osed approach will be compared to those o bta ined using the so- called One-Day-ahead Naive Predictor (ODNP), i.e., ˆ P OD N P ( j ) = P m d − 1 ( j ) , (36) where P m d − 1 ( j ) denotes the measur e of gener ated pow er re corded during the day befor e at the same time of day . 14 7. V alidation on simulated data Before assess ing the p erformance of the propose d appro ac h on mea sured data s e ts, in this section we present pre limina ry tests on simulated noisy da ta in order to highlight the pro perties of the prop osed approa c h. Nominal weather data (i.e., cloud cover N ( k ) and temper ature T ( k )) hav e been generated us ing a n empirical proc e dure which r eproduce s a realistic cloud cover and temper ature s e asonal trend in Italy . The co r resp o nding generated p ow er time series P m ( k ) is then computed a ccording to mo dels (1) and (4) using the v alues o f the physical para meters: µ 1 = 0 . 92 , µ 2 = − 1 . 237 × 10 − 4 , µ 3 = − 2 . 99 × 10 − 3 , µ 4 = − 0 . 3 , µ 5 = − 0 . 25 , ψ = 27 ° , ζ = 0 ° . (37) The pla n t nomina l p o wer P nom is a ssumed equal to 920 kW. Genera ted data ar e a ssumed to span the p erio d January 1s t (day d = 1) to Decem b er 31s t (day d = 365) and are sampled with a time step τ s = 15 min utes. The overall data se t is then g iv en by D = { D ( k ) = { P m ( k ) , T ( k ) , N ( k ) } , k ∈ K} , where the set of time indices K spans the time hor izon from the star ting day d to the final da y d . As previously sta ted, only the indices k corr esponding to hours of lig h t (i.e., for which the theo retical clear-sk y irra diance is nonzero) w ere consider ed. A mo del estimation s tep is therefore p erformed at each instant k us ing the data sample D ( k ). 7.1. Sc enarios Five different simulation scenar ios (i) − (v) hav e b een considered. Within each scenario, differen t noise and data pro cessing conditions identified by a Setup I D (SID) hav e b e e n simulated. The scenarios a re a rranged as follows: (i) Neither noise nor quantization is a dded. Raw data (SID 0) and data a veraged ov er 1 h (SID 1) (see T able 1) were consider ed. Av erag ing is intro duced in order to mimic the typical setting where hourly updates of weather and gener ation data ar e provided by the meteo rological service and the DSO, r e s pectively . (ii) A zero-mea n Gaussian white nois e with v ariance σ 2 N is added to N ( k ), N ( k ) is then quantized so that N ( k ) ∈ { 0 , 0 . 1 , 0 . 2 , . . . , 1 . 0 } , r esulting in a noise v a riance σ 2 N ,tot (SID 2 − 5 in T able 2). Data is finally av erage d ov er 1 h. 15 SID Description 0 No noi se added , raw data sampled at 15 mi n 1 No noi se added , data averaged ov e r 1 h T able 1: Setup descriptions for scenario (i). SID 2 3 4 5 3 σ N 0 . 0 0 . 1 0 . 5 1 . 0 σ N ,tot 0 . 289 0 . 2 91 0 . 33 3 0 . 441 T able 2: Setup description for scenario (ii). SID 6 7 8 9 10 11 12 3 σ N 0 . 0 0 . 0 0 . 0 0 . 0 0 . 0 0 . 0 0 . 3 3 σ P (kW) 10 50 1 00 0 0 0 50 3 σ T ( ◦ C) 0 0 0 1 3 5 3 T able 3: Setup description for scenarios (iii), (iv) and (v). (iii) A zero -mean Gaussian white noise with v ariance σ 2 P is added to P m ( k ), N ( k ) is quantized and data are av eraged ov er 1 h (SID 6 − 8, T a ble 3). (iv) A zero-mea n Gaussian white noise with v ariance σ 2 T is added to T ( k ), N ( k ) is qua ntized a nd data ar e averaged ov er 1 h (SID 9 − 11, T able 3). (v) A zer o-mean Ga ussian white nois e is a dded to all the re gression v ariables , N ( k ) is qua n tized and data are av eraged ov er 1 h (SID 12 , T a ble 3 ). Notice that here data av eraged ov er 1 h mea ns that for ea c h hour h there is only o ne sample computed a s mean of the data a cquired dur ing h . The per formance of the prop osed metho d on the scenario s intro duced above has b een ev alua ted with reference to the day-ahead (D A) forecast intro duce d in Section 5. F or c o mputing the forecasts, weather da ta from the da ta set D hav e been employ ed. In particular, at each time instan t k d 0 corres p onding to 6 am on each day , we ev aluate the D A fore c a st ˆ P DA ( d ) according to (2 5)-(26) with ˆ T ( j | k d 0 ) = T ( j ) , ˆ N ( j | k d 0 ) = N ( j ) , j ∈ K d . The ev aluation o f the p erformance indices star ts at day d 0 = 18, in or der to first ensure a rough degree of adapta tio n of the mo del par a meters. The following p erformance indices were consider e d: • RM S E d : the RMSE rela tiv e to the forecas t p ertaining to each day d , i.e., (29) e v aluated fo r ˆ P ( j ) = ˆ P ( j | k d 0 ; k d − 1 ), K = K d ; 16 • RMSE, MBE, MAPE , and determination co efficien t R 2 , ev a lua ted o ver the whole sim ulation, i.e., ov er the union of a ll DA foreca sts from day d 0 to day d . F or each SID, 10 s im ulations were p erformed, in which P m ( k ), T ( k ) a nd N ( k ) were obtained by using different no is e realiza tions. The average of p erformance indices over all the simulations were co nsidered. In each scenario, the initial v a lues ˆ µ 1 (0), ˆ µ 2 (0), ˆ µ 3 (0), ˆ µ 4 (0) and ˆ µ 5 (0) are chosen as 75% of the real v alues in (37), while ψ and ζ a re equal to their nominal v a lues. l (0) has b een initialized to 0 . 01 for b oth RLS algor ithms and EKF, r has b een initialized to 10 4 , which is a conserv ative choice in view of the v ar iance of the noise added to p ow e r measurements. 7.2. R esults Figure 4 shows one instance of the parameter estimate evolution for the N5 mo del for SID 0 and SID 1 . Using noise-free r a w data (SID 0) the para meters converge exactly to the real v alues (Figure 4a). In ca se of SID 1 the parameters c o n verge as well with sligh t deviations of ˆ µ 4 and ˆ µ 5 from the nominal v alues (Figure 4b). Note that in this c ase, the para meters take consider ably longer to converge. 0 200 100 300 1 0.6 0.8 0.7 0.9 0 20 40 10 30 0e00 −2e−04 2e−04 0 200 100 300 0 −0.005 0.005 0 200 400 100 300 −0.3 −0.25 0 200 400 100 300 −0.2 −0.26 −0.24 −0.22 −0.18 P S f r a g r e p la c e m e n t s ˆ µ 1 ˆ µ 2 ˆ µ 3 ˆ µ 4 ˆ µ 5 ψ ζ (a) Sc enario (i) - SID 0 , noise-fr ee data 0 2 000 4 000 1 0.6 0.8 1.2 1.4 0 2 000 4 000 0e00 −2e−04 2e−04 −1e−04 1e−04 0 2 000 4 000 0 −0.06 −0.04 −0.02 0 2 000 4 000 0 −0.4 −0.2 −0.3 −0.1 0 2 000 4 000 0 −0.6 −0.4 −0.2 P S f r a g r e p la c e m e n t s ˆ µ 1 ˆ µ 2 ˆ µ 3 ˆ µ 4 ˆ µ 5 ψ ζ (b) Sc enario (i) - SID 1 , noise-fr ee aver age d data Figure 4: N 5 mo del parameter estimation vs. time (iterations), scenario (i). Bl ac k dashe d lines are the nominal v alues of each parameter as rep orted in (3 7). Figures 5a, 5b and 5c show three examples of parameter es tima te evolution for the N5 model under the la rgest noise v ar iance on N , P a nd T , resp ectiv ely , Fig ure 5d shows the s ame plot when 17 noise is added to all the regress ion v ariables . Parameters ˆ µ 4 and ˆ µ 5 are a pparent ly q uite sensitive to the noise added to N . Indeed, under the conditions describe d b y the SID 5, for whic h σ N ,tot = 0 . 441, ˆ µ 4 is under estimated and ˆ µ 5 is overestimated (Figures 5a). All the o ther para meters always tend to conv erg e to their r eal v alues. 0 2 000 4 000 1 0.6 0.8 1.2 0 2 000 4 000 0e00 −2e−04 −3e−04 −1e−04 0 2 000 4 000 0 −0.02 0.02 0.04 0 2 000 4 000 0 −0.4 −0.2 −0.5 −0.3 −0.1 0 2 000 4 000 0 −0.6 −0.4 −0.2 P S f r a g r e p la c e m e n t s ˆ µ 1 ˆ µ 2 ˆ µ 3 ˆ µ 4 ˆ µ 5 ψ ζ (a) Sc enario (ii) - SID 5 , noise adde d only to N . 0 2 000 4 000 1 0.6 0.8 1.2 0 2 000 4 000 0e00 −2e−04 2e−04 −3e−04 −1e−04 1e−04 0 2 000 4 000 0 −0.04 −0.02 0 2 000 4 000 0 −0.2 −0.3 −0.1 0 2 000 4 000 0 −0.8 −0.6 −0.4 −0.2 P S f r a g r e p la c e m e n t s ˆ µ 1 ˆ µ 2 ˆ µ 3 ˆ µ 4 ˆ µ 5 ψ ζ (b) Sc enario (iii) - SID 8 , noise adde d only to P . 0 2 000 4 000 1 0.6 0.8 1.2 0 2 000 4 000 0e00 −2e−04 −1e−04 1e−04 0 2 000 4 000 0 −0.04 −0.02 −0.03 −0.01 0 2 000 4 000 0 −0.2 −0.3 −0.1 0.1 0 2 000 4 000 0 −0.8 −0.6 −0.4 −0.2 P S f r a g r e p la c e m e n t s ˆ µ 1 ˆ µ 2 ˆ µ 3 ˆ µ 4 ˆ µ 5 ψ ζ (c) Sc enario (iv) - SID 11 , noise adde d only to T . 0 2 000 4 000 1 0.6 0.8 1.2 0 2 000 4 000 0e00 −2e−04 −1e−04 1e−04 0 2 000 4 000 0 −0.04 −0.02 −0.03 −0.01 0 2 000 4 000 0 −0.4 −0.2 −0.3 −0.1 0.1 0 2 000 4 000 0 −0.8 −0.6 −0.4 −0.2 P S f r a g r e p la c e m e n t s ˆ µ 1 ˆ µ 2 ˆ µ 3 ˆ µ 4 ˆ µ 5 ψ ζ (d) Sc enario (v) - SID 12 , noise adde d to N , P and T . Figure 5: N5 model parameter estimation (blue line) vs. time (ite rations) and using da ta from scenario ( ii) to scenario (v). Black dashed li nes are the nominal v alues of eac h parameter as rep orted i n (37). Figure 6 depicts the evolution o f the daily RMSE ( RM S E d ) for incr easing SID over the whole simulation (365 da ys), r elativ e to L, N5, N6 mo dels and to the ODNP . Such er ror is compared with 18 the standard devia tion of simulated p o wer data. In vie w o f the fact that para meter convergence behaves differently for different SIDs and in o rder to hav e a consisten t comparison, the p erforma nce indices for a ll cases were computed star ting from day 15. Results show that the pr opos ed approach is more sensitive to noise added to N , than to noise added to P or T . The R M S E d computed with SID 5 (Figure 6a) is considerably higher than the RM S E d computed with SID 8 and 11 (Fig ures 6b and 6 c), which are really similar to each other. Cle a rly , under any of the tested nois e co nditions, the RM S E d achiev ed by a ll three models is by far low er than bo th the RM S E d obtained using the ODNP and the standar d devia tion of the simulated p o wer. Performance indices RMSE, MAPE, MB E , and R 2 ov er the who le simulation a re re p orted and compared for each SID a nd each mo del in Figure 7. R M S E a nd M AP E slo wly inc r ease a nd R 2 slowly decreases with σ N (SID 2 − 5) and σ P (SID 6 − 8), while they are almos t constant when σ T (SID 9 − 11 ) incr eases. Finally , note that M B E gr o ws when σ P is inc r eased and that it decrea ses bo th when σ T is inc r eased and w he n σ N grows. In Figure 8, the p o wer forecasts provided b y mo dels L, N5 and N6 during a simulation run with SID 12 a re qualitatively compar ed to mea sured data during three consecutive days. The picture refers to typical examples of pa rtially c loudy days. 19 0 200 400 100 300 20 40 60 80 120 140 160 180 220 240 260 280 320 340 360 380 0 200 100 300 50 150 250 P S f r a g r e p la c e m e n t s RM S E d (kW) time ( day) (a) SID 5 , noise adde d only to N . 0 200 400 100 300 20 40 60 80 120 140 160 180 220 240 260 280 320 340 360 380 0 200 100 50 150 250 P S f r a g r e p la c e m e n t s RM S E d (kW) time ( day) (b) SID 8 , noise adde d only to P . 0 200 400 100 300 20 40 60 80 120 140 160 180 220 240 260 280 320 340 360 380 0 200 100 300 50 150 250 P S f r a g r e p la c e m e n t s RM S E d (kW) time ( day) (c) SID 11 , noise adde d only to T . 0 200 400 100 300 20 40 60 80 120 140 160 180 220 240 260 280 320 340 360 380 0 200 100 50 150 250 P S f r a g r e p la c e m e n t s RM S E d (kW) time ( day) (d) SID 12 , noise adde d to N , P and T . Figure 6: Comparison b et ween the standard deviation of the simulated p o wer (dashed dot ted line), t he RM S E d computed using the ODNP ( dashed line), L model (blac k l ine), N5 mo del (blue line) and N6 mo del (red line). 20 0 1 2 3 4 5 6 7 8 9 10 11 12 0 20 40 60 0 1 2 3 4 5 6 7 8 9 10 11 12 0% 20% 40% 60% P S f r a g r e p la c e m e n t s RM S E (kW) M AP E M B E ( k W ) R 2 SID SID (a) RM S E and M AP E . 0 1 2 3 4 5 6 7 8 9 10 11 12 0 2 8 4 6 0 1 2 3 4 5 6 7 8 9 10 11 12 0.90 0.95 1.00 P S f r a g r e p la c e m e n t s R M S E ( k W ) M A P E M B E (kW) R 2 SID SID (b) M B E and R 2 . Figure 7: Performance indices vs. noise co mputed using simulat ed da ta. Red bars represen t the N5 model, blue bars the N6 model, and green bars th e L mo del. 62 64 63 61.4 61.6 61.8 62.2 62.4 62.6 62.8 63.2 63.4 63.6 63.8 0 200 400 600 800 100 300 500 700 0.0 0.5 1.0 0.0 0.5 1.0 P S f r a g r e p la c e m e n t s kW time ( day) Figure 8: Scenario (v) - SID 12, comparis on b et we en the simulat ed pow er (long dashed line), N5 mo del DA f orecast (red line), N6 model D A f orecast (blue li ne, whic h almost o verlaps with the red line) and L model DA f orecast (green line). S hort dashed line represen ts the CCI (scale to the r igh t). 21 8. V alidation on exp erimenta l data 8.1. Exp eriment set up The prop osed pr ocedur e has b een v a lidated using data fro m a photov oltaic plant of nominal power P nom = 920 kW p located in Sardinia (Italy). Av ailable data consist of t wo da tasets spanning F ebruar y 2nd (33 r d day of year) to May 1st (122nd day of year): • D 1 : a set o f data provided by the pr iv ate pr oducer running the plan t, which collects hourly samples o f averaged measure d p ow er and o ne day-ahead forecas ts of air temp e r ature; • D 2 : a set of weather rep orts including CCI, ev a luated by a meteoro logical station lo cated 20 km aw ay from the pla n t. This data s e t do es no t hav e a re g ular s ampling time. D 1 and D 2 hav e b een prepr oces sed in order to hav e a single dataset D o f hour ly data. If there a re more then one CCI v alues in D 2 for the same hour, an a veraged CCI is computed for that hour. On the contrary , if the CCI r eport in D 2 is missing for a given hour, then the measure d p ow er a nd the forecast of a ir temp erature in D 1 for that hour are dropp ed 2 . The da ta set employ ed is therefore given by the time ser ies D = { D ( k ) = { P m ( k ) , T ( k ) , N ( k ) } , k ∈ K} , where the set of time indices K spans, with a sampling time τ s = 1 h, fro m the starting day d = 33 to the final day d = 122. Also in this cas e, only the indices k corresp onding to hour s of light were considered. Mo del para meter up dates ar e c o mputed hourly using the data sample D ( k ). The initial v alues of the parameter s have b een chosen as follows 3 : ˆ µ 1 (0) = P nom / 1000 , ˆ µ 2 (0) = − 1 . 34 × 10 − 4 · ˆ µ 1 (0) , ˆ µ 3 (0) = − 3 . 25 × 10 − 3 · ˆ µ 1 (0) , ˆ µ 4 (0) = 0 . 784 , ˆ µ 5 (0) = − 1 . 344 , ψ = 27 ° , ζ = 12 ° , l (0) = 10 , r = 10 4 . 2 F or the sake of complete ness, and to further support the claim that one day-ahea d temper ature forecasts can b e used in our procedure in place of ac tual data, temp erature measuremen ts collect ed at the meteorological station ha v e been compared with the f orecasts included in D 1 . The RMSE and M BE computed on the ov erall data s et turned out to be equal to 1 . 9 ◦ C and 0 . 9 ◦ C , resp ectiv el y . 3 Concerning the ori en tation angles ψ and ζ , whic h were not av ailable, their v alues ha ve b een guessed b y visually comparing the generation curv e during a clear- sky da y with the clear-sky irradiance model. 22 8.2. Performanc e evaluation The p erformance of the prop osed metho d on the da ta set introduced a bov e has b een ev alua ted with re ference to the tw o standa r d fo recast t yp es int ro duced in Section 5, i.e., • DA foreca st . In this case, the for ecasts a nd the p erformance indices R M S E d , RMSE, MAPE, MBE and R 2 were computed using the sa me mo dality as in the previous section. • HA foreca st . The HA for ecasts ar e computed a ccording to (25) -(27) ba s ed on the da ta set D using ˆ T ( j | k ) = T ( j ) , ˆ N ( j | k ) = N ( j ) , j > k . In order to ev aluate the performa nce, a single 7 - hour HA forecast p er day has b een considered, and the e r ror measures RMSE, MAPE , MBE and R 2 hav e bee n computed over the the union of all such forecasts. The perfo rmance of the prop osed metho d has b een compared with that obtained using the following a pproaches. • ODNP (for DA forecasts only , clearly not suitable in the HA case). • Autoregr essiv e prediction models using only measur emen ts of generated pow er a nd only cloud cov er data, resp e ctiv ely . These metho ds r epresent a more realis tic benchmark with resp ect to ODNP a nd a re applica ble to b oth DA and HA. T o this end, the following tw o mo dels hav e bee n tested: – PVGM : a utoregressive (AR) mo del o f ge ne r ated p ow e r with regr ession horizon equal to 12 ho ur s. The para meter vector is g iven by a = [ a 1 . . . a 12 ] T and the regr essor at time j is e x pressed b y ψ ( j ) = [ P ( j − 1) . . . P ( j − 12)] T . Then the predictio n of genera ted pow er for time instant j , co mputed at time k using parameter s av ailable at time q ≤ k is: ˆ P ( j | k , q ) = ˆ a ( q ) T ψ ( j ) , where: ˆ a ( q ) is the estimate of a av ailable at time q , and for i = 1 . . . 12, the i -th term P ( j − i ) of ψ ( j ) is equal to P m ( j − i ) if j − i ≤ k , and equa l to ˆ P ( j − i | k , q ) otherw is e. – CCD: autor egressive with exogeno us input (ARX) mo del of gener a ted p ow er with re- gressio n hor izon equal to 2 hours. The inputs a re cloud cover data a nd the clear -sky 23 mo del (8) 4 . The parameter vector is given by b = [ b 1 . . . b 6 ] T and the regr essor at time j is expres s ed by ξ ( j ) =  I 0 ( j ) . . . I 0 ( j − 2 ) N ( j ) . . . N ( j − 2)  T . Then the prediction o f generated p ower ˆ P ( j | k , q ) is: ˆ P ( j | k , q ) = ˆ b ( q ) T ξ ( j ) , where: ˆ b ( q ) is the estimate of b av ailable at time q , a nd for i = 0 , 1 , 2, the i -th term N ( j − i ) o f ξ ( j ) is equa l to N m ( j − i ) if j − i ≤ k , and equa l to ˆ N ( j − i | k ) otherwis e. The para meters of the tw o models above have been estimated in the sa me recursive framew ork as the prop osed par a metric mo del by means of Linear Least Sq ua res. • ANN-based approa c h, using the MultiLayer Perceptro n (MLP) ar c hitecture. In this cas e we are lo oking for a fair compar ison in terms of mo del co mplexit y , therefor e tw o architectures with the s ame num be r of pa rameters as the N6 and L para metric mo dels, res p ectively , hav e bee n tested using the training alg orithms implemen ted in the MA TLAB Neur al Net work T o olb ox. Two ar c hitectures hav e b een selected: – MLP1: 1 hidden unit, 6 par ameters, – MLP2: 2 hidden units, 11 par ameters. Note that no s uc h comparison is poss ible with the N5 mo del since a 5 -parameter MLP cannot be synthesized. The netw o rks have b een trained in the same recursive framework as the prop osed par a metric mo del, but since a training step for the MLPs ca nnot be p erformed for each inco ming da ta sample D ( k ), the algo rithm ha s b een mo dified in o rder to build suitable learning sets fo r net work training . In pa rticular, a new learning se t is constructed e v ery 75 samples. 8.3. R esults and discussion In Figures 9 and 10 the ev olution of the par ameter es tima tes are rep orted for the N5 and N6 mo dels. All the para meters tend to conv erge for thes e as well a s for the L mo del, not depicted to save spa ce. The estimated parameter v alues and their standard devia tion at the end of the exp erimen t ar e shown b elow: 4 Cloud co ve r index alone cannot poss ibly pr o vi de en ough inform ation. 24 ( N 5) ˆ µ = h 1 . 03 − 2 . 54 × 10 − 4 − 2 . 12 × 10 − 4 1 . 29 × 10 − 1 − 8 . 42 × 10 − 1 i T , σ ˆ µ = h 1 . 72 × 10 − 2 1 . 58 × 10 − 5 7 . 30 × 10 − 4 6 . 02 × 10 − 2 8 . 00 × 10 − 2 i T , (38) ( N 6) ˆ µ = h 1 . 01 − 1 . 60 × 10 − 4 − 1 . 75 × 10 − 3 5 . 21 × 10 − 1 − 1 . 19 − 2 . 53 × 10 − 4 i T , σ ˆ µ = h 8 . 88 × 10 − 2 1 . 57 × 10 − 5 6 . 24 × 10 − 4 2 . 52 × 10 − 2 3 . 82 × 10 − 3 1 . 36 × 10 − 5 i T , (39) ( L ) ˆ µ = h 1 . 39 − 9 . 54 × 10 − 1 − 2 . 53 × 10 − 1 − 3 . 49 × 10 − 4 6 . 00 × 10 − 4 − 1 . 79 × 10 − 3 3 . 45 × 10 − 3 − 1 . 75 × 10 − 3 − 1 . 40 × 10 − 2 3 . 87 × 10 − 2 − 2 . 47 × 10 − 2 i T , σ ˆ µ = h 1 . 66 × 10 − 2 9 . 18 × 10 − 2 9 . 27 × 10 − 2 2 . 04 × 10 − 5 2 . 11 × 10 − 4 5 . 94 × 10 − 4 9 . 85 × 10 − 4 5 . 79 × 10 − 4 1 . 21 × 10 − 3 9 . 10 × 10 − 3 1 . 02 × 10 − 2 i T . (40) 0 200 400 600 800 0 1 0.5 1.5 0 200 400 600 800 0 −0.002 0.002 −0.001 0.001 0 200 400 600 800 0 0.02 0.04 0.01 0.03 0 200 400 600 800 0 1 0.5 0 200 400 600 800 0 −2 −1 1 P S f r a g r e p la c e m e n t s ˆ µ 1 ˆ µ 2 ˆ µ 3 ˆ µ 4 ˆ µ 5 ψ ζ Figure 9: N5 mo del parameter estimation. Paramete r estimates vs. time (iterations). Notice tha t all standard deviations are at least o ne o rder of magnitude lo wer than the r espective parameter estimates. The o nly ex ception concer ns the par ameter estimate ˆ µ 3 for the N5 mo del. This mo del, which is characterized by a minimal num b e r of pa rameters, a ppa ren tly do es no t fully capture the dep endency be tw een the genera ted p o wer and temp erature, thus gener ating p ersistent parameter drifts. On the contrary , the s ligh tly ov erpara meter ized N6 mo del shows m uch b etter conv ergence prope rties, as it is also apparent fr om the comparison of Figures 9 and 10. F ur thermore, the L mo del shows higher para meter standard deviatio ns, on av erage, with resp ect to N6. In order to analyze the impact of ov erpara metr ization on the consistency of the par a meter estimates w ith the physical mo del (2 ) , w e compa re the v alue s of ˆ µ 1 and the co rrection terms ˆ η 2 = ˆ µ 2 / ˆ µ 1 and ˆ η 3 = ˆ µ 3 / ˆ µ 1 obtained from (38), (39) and (40) . The estimates of ˆ µ 1 provided by 25 0 200 400 600 800 0 1 0.5 1.5 0 200 400 600 800 0 −0.002 0.002 −0.001 0.001 0 200 400 600 800 0 −0.1 −0.08 −0.06 −0.04 −0.02 0.02 0 200 400 600 800 1 0.2 0.4 0.6 0.8 0 200 400 600 800 −1 −1.4 −1.2 −1.5 −1.3 −1.1 0 200 400 600 800 0 −0.002 0.002 −0.001 0.001 P S f r a g r e p la c e m e n t s ˆ µ 1 ˆ µ 2 ˆ µ 3 ˆ µ 4 ˆ µ 5 ψ ζ ˆ µ 6 Figure 10: N6 mo del paramet er estimation, parameter estimates vs. t ime (iterations). mo dels N5 and N6 are slig h tly higher than the nominal p ow er / irradiance ga in P nom / 1000 = 0 . 92, while the v a lue obtained using mo del L is remar k ably higher . In T able 4, ˆ η 2 and ˆ η 3 are co mpared with their typical v alues rep orted in (3). Notice that N6 is the o nly mo del for which η 2 ∈ S 2 and η 3 ∈ S 3 , while for N5 ˆ η 3 is hig her than the upp er b ound of S 3 , which apparently means that µ 3 is ov erestimated. On the contrary , mo del L tends to underestimate µ 3 . ˆ η i S i N5 N6 L ˆ η 2  − 2 . 5 × 10 − 4 , − 1 . 9 × 10 − 5  − 2 . 47 × 10 − 4 − 1 . 58 × 10 − 4 − 2 . 51 × 10 − 4 ˆ η 3  − 4 . 8 × 10 − 3 , − 1 . 7 × 10 − 3  − 2 . 06 × 10 − 4 − 1 . 73 × 10 − 3 − 10 . 1 × 10 − 3 T able 4: Corr ection terms η 2 and η 3 computed using the estimated paramet er vect ors rep orted in (38), (39) and ( 40). Notice that in mo del N6 the ov erpar ametrization is given o nly b y the introductio n of parameter µ 6 , defined as the pro duct b et ween µ 2 and µ 4 . F rom (39), w e have that ˆ µ 2 = − 1 . 60 × 10 − 4 , ˆ µ 4 = − 1 . 19 , a nd ˆ µ 6 = − 2 . 53 × 10 − 4 . Therefore ˆ µ 2 ˆ µ 4 = − 9 . 12 × 10 − 4 , which is the same order of magnitude as ˆ µ 6 , thus suggesting that the estimates a re almo st consistent. A s imila r cons is tency is no t shown by mo del L, despite its go o d pe r formance on the DA and HA for ecasts, as shown in the seq uel. 26 In order to further ev alua te the c onsistency of mo dels N5 a nd N6, we ana lyze the CCF curve C ( N ) in (5 ) obta ined using the estimated v alues of µ 4 and µ 5 provided by such mo dels. In Figure 11 we rep ort the plo ts of the C ( N ) using the resp ective v alues of ˆ µ 4 and ˆ µ 5 . It turns out that bo th the estima ted CCFs are consistent with the seas o nal trend for the mediterr anean climate [3 4]. 0 0.2 0.4 0.6 0.8 1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 P S f r a g r e p la c e m e n t s C ( N ) N Figure 11: Comparison bet ween the CCF curv es obt ained using the estimated v alues of µ 4 and µ 5 prov ided b y models N 5 (blue li ne) and N6 (red line), x-marks deno te the maximum point s. As far as the forecasting per formance ev aluation is concerned, all e r ror meas ures were computed ov er the p erio d starting fr o m day d 0 = 57 in order to guara n tee, as in the simulated da ta case, at least a rough adaptation of the mo del pa rameters. T able 5 s ummarizes the per formance indices achiev ed by the prop osed parametric mo dels, com- pared with MLPs, ODNP , PVGM and CCD o n bo th DA and HA forec asts. Since the p erformance of the MLPs may b e quite sensitive to the initial co nditio ns, 30 sim ulations of MLP1 a nd MLP2 with the s ame data set hav e be e n per fo rmed with suitably genera ted random initial conditions 5 , a nd mean v alue/v ariance of all per formance indices ov er all sim ulations hav e b een r e ported. Clear ly , all parametric mo dels and MLPs p erform la rgely b etter than the ODNP and PVGM. On av erage , MLPs sho w er r ors 20% higher than parametric models. The CCD mo del y ie lds intermediate indices 5 MLP la yer weigh ts and biases are initialized according to the Nguy en-Widrow algo rithm. This algori thm c ho oses v alues in order to distribute the act ive region of each neuron in the lay er randomly but ev enly across the lay er’s input space. 27 betw een o ur mo dels and MLPs . The p erformances o f the thre e prop osed parametric mo dels are ov erall compar able. Notice that in netw ork op eration, p erformance indices qua lifying fo r ecast ac- curacy , e.g., RMSE and MAPE, are genera lly normalized with resp ect to the nominal plant p o wer. This normalizatio n is useful in view of how the DSO quantifies uncer tain ty in p o wer generatio n. Actually , when dealing with hundreds or thousa nds of distr ibuted units, either for maintenance activities or op erational r equiremen ts, the DSO takes into acc o un t how muc h each sp ecific unit is exp ected to pro duce with r espect to its nominal p o wer. Indeed, netw ork o pera tion requires the DSO to k no w at any time the overall p ow e r uncertaint y p ertaining to the whole distributed generatio n. Such uncertaint y can b e easily estimated from the nomina l capacity of each plant (known a-pr io ri) and the nor malized plant forecast erro r (computed o ff-line on historical data a nd usually with a smaller time r e s olution, typically , weekly or monthly). F or this reaso n, we consider the norma liz e d error s R M S E N P and M AP E N P . The results o btained for the la tter indices, whic h a r e around 8-10%, can b e consider ed acceptable fo r most problems in netw ork o pera tion. Figure 12 shows the daily RMSE ( R M S E d ) achiev ed by the propo s ed models in the D A foreca st, compared with ODNP , P V GM, and CCD. In Figures 13, s uch erro r is compared with those achiev ed by the MLPs in the DA forecast. Figure 14 depicts the comparis o n of the R M S E d achiev ed by the v ar ious models a nd benchmarks in the HA fo r ecast. In the a b ov e figures, all erro r s are also compared with the standar d deviatio n o f measured p ow er data. Both in T able 5 and Fig ure 14, data for the O DNP in the HA cas e are not rep orted a s the ODNP is unsuitable as a r eference mo del. The a pparent ly high normaliz e d forecast er rors (in particular the MAPE ) are mainly due to the fact that the CCI data refer to a weather station lo cated a b out 20 km awa y from the pla n t site. Another sour c e of for ecast disper s ion is the rough C CI data resolution. Bo th these tw o sources of uncertaint y deter io rate for ecasting accura cy as muc h as the clo udiness at the plant site changes ra pidly in time. It is worth remarking, how ever, that the pro posed mo dels use a very limited amount of informatio n which is easily a nd cheaply av a ilable, and that s uc h r a w data is informative enough for a reliable estimation o f the gener a tion mo del, as demonstra ted by the parameter conv ergence pr o perties highlighted in Figur e 1 0. F urthermore, it is expe c ted that the prop osed approach c a n b e fruitfully a dapted to large-s cale agg r egation of plants over macro areas, yielding significant p erformance improv ements with resp ect to the single plant ca se [3 6], with no additional da ta r equirements (except g eneration meas uremen ts). 28 Performance Indece s N5 N6 MLP1 Mean (V ariance) L MLP2 Mean (V ariance) ODNP PVGM CCD D A F or ecast RM S E (kW ) 109 11 0 133 (2 . 05) 110 127 (19 . 0) 227 283 121 M AP E 50 . 0% 5 0 . 7% 64 . 6% (5 . 17) 48 . 8% 66 . 9% (31 . 8) 87 . 4% 165 . 3 % 61 . 4% M B E (kW ) − 1 . 275 4 . 10 8 − 3 . 42 (1 . 38) 4 . 455 − 1 . 15 (37 . 7) 27 . 2 127 . 4 1 1 . 2 R 2 0 . 828 0 . 827 0 . 746 (2 . 98 × 10 − 5 ) 0 . 825 0 . 770 (2 . 71 × 10 − 4 ) 0 . 254 − 0 . 157 0 . 787 N RM S E 0 . 414 0 . 416 0 . 511 (2 . 94 × 10 − 5 ) 0 . 418 0 . 484 (2 . 28 × 10 − 4 ) 0 . 864 1 . 076 0 . 462 RM S E N P 0 . 119 0 . 119 0 . 146 (2 . 41 × 10 − 6 ) 0 . 120 0 . 139 (1 . 88 × 10 − 5 ) 0 . 248 0 . 308 0 . 132 M AP E N P 7 . 51% 7 . 7 8% 9 . 85% (1 . 70 × 10 − 2 ) 7 . 71% 9 . 28% (1 . 24 × 10 − 1 ) 15 . 28% 23 . 9 6% 9 . 49% HA F orecast RM S E (kW ) 103 10 4 129 (6 . 79) 105 125 (14 . 4) - 223 12 0 M AP E 45 . 1% 4 6 . 7% 66 . 9% (11 . 5) 42 . 7% 75 . 8% (183) - 1 77 . 3% 61 . 4% M B E (kW ) 11 . 74 11 . 97 − 2 . 03 (2 . 63) 12 . 53 − 4 . 01 (45 . 1) - 62 . 9 24 . 8 R 2 0 . 846 0 . 841 0 . 756 (9 . 73 × 10 − 5 ) 0 . 839 0 . 774 (1 . 89 × 10 − 4 ) - 0 . 274 0 . 788 N RM S E 0 . 393 0 . 399 0 . 498 (9 . 88 × 10 − 5 ) 0 . 401 0 . 477 (1 . 77 × 10 − 4 ) - 0 . 852 0 . 460 RM S E N P 0 . 112 0 . 114 0 . 142 (8 . 03 × 10 − 6 ) 0 . 115 0 . 136 (1 . 44 × 10 − 5 ) - 0 . 243 0 . 131 M AP E N P 6 . 88% 7 . 0 7% 9 . 05% (3 . 28 × 10 − 2 ) 7 . 16% 8 . 83% (5 . 28 × 10 − 2 ) - 1 8 . 66% 9 . 73% T able 5: Performance comparison of parametric mo dels, ODNP , PV GM, CCD and MLP computed starting from da y 57. Results are group ed by model complexity and purp ose in order to simplify the comparison: 5-parameter mo del (N5) i n red column, 6-parameter mo dels (N6 and MLP1) in blue columns, 11-parameter mo dels (L and MLP2) in green column, and benc hmark mo dels (ODNP , PVGM and CCD ) in gray columns. 29 100 60 80 120 50 70 90 110 130 55 65 75 85 95 105 115 125 0 200 400 100 60 80 120 50 70 90 110 130 55 65 75 85 95 105 115 125 0 200 400 P S f r a g r e p la c e m e n t s time ( day) time ( day) RM S E d (kW) RM S E d (kW) Figure 12: Comparison b et ween standard deviation of the measured p ow er (dashed dotted line) and the RM S E d computed on D A forecasts. T op: ODNP ( dashed line), PVGM (d asehd red line) and CCD (dashed blue line). Bottom: L mo del (green line), N5 mo del (red line) and N6 model (blue line). 100 60 80 120 50 70 90 110 130 55 65 75 85 95 105 115 125 0 200 100 300 50 150 250 350 P S f r a g r e p la c e m e n t s time ( day) RM S E d (kW) Figure 13: D A forecast. Comparison b et wee n standard deviation of the measured p ow er (dashed dotted li ne), the RM S E d computed using the N6 model (blue line), M LP1 (blue dashed l ine), L mo del (green line) and MLP2 (green dashed line). 100 60 80 120 50 70 90 110 130 55 65 75 85 95 105 115 125 0 200 400 100 60 80 120 50 70 90 110 130 55 65 75 85 95 105 115 125 0 200 400 P S f r a g r e p la c e m e n t s time ( day) time ( day) RM S E d (kW) RM S E d (kW) Figure 14: HA forecast. Comparison betw een the stan dard deviation of the measured p o wer (d ashed dotted line) and the RM S E d . T op: PVGM (dasehd red line) and CCD (dashed blue line). Bottom: computed using the N5 m odel (red line), N6 model (blue line), MLP1 (b lue dashed line), L model (green line) and MLP2 (green dashed l ine). 30 In Figure 1 5 the D A forecas ts provided by N6 mo del and MLP 1 during three different da ys and under three different weather conditions is co mpared with the measures of g enerated p ow er . 90 91 90.2 90.4 90.6 90.8 0 200 400 600 800 100 300 500 700 900 0.0 0.5 1.0 85 86 85.2 85.4 85.6 85.8 0 200 400 600 800 100 300 500 700 900 0.0 0.5 1.0 104 105 104.2 104.4 104.6 104.8 0 200 400 600 800 100 300 500 700 900 0.0 0.5 1.0 P S f r a g r e p la c e m e n t s kW kW kW time ( day) time ( day) time ( day) Figure 15: DA forecast. Comparis on b et ween the measured p o wer (black l ine), N6 mo del forecast (blue line), and MLP1 forecast (dashed bl ue line). Dashed dotted line is the m easured CCI. F rom right to left, an almost clear-s ky da y , an o vercast da y and a uniformly ov ercast day are de picted. The pr ocedur e has also b een tested using the last tw o weeks of data for v alidation purp oses only , that is , mo del parameter v alues were fro zen tw o weeks b efore the end of the data set, then offline D A and HA fo r ecasts for the who le tw o following weeks were gener ated a ccording to the mo dalities ab ov e, a nd the res p ective p erforma nc e indices were computed. It turned o ut that the per formances of the v arious mo dels in terms of all indices were a lmost ident ical to those obta ined with online par ameter adaptation (the difference in terms of all er r or measur es was b elow 1%). In this resp ect, it is worth r emarking that online parameter estimation has the adv antage of capturing seasona l pa rameter v aria tions (esp ecially a s far a s C ( N ) is concer ned). As far as the computational time is concerned, the MLP update algorithms take generally longer to execute with resp ect to the par ameter update pro cedure. More sp ecifically , a MLP training s tep using a K -sized sample data s e t re q uires long er computational time than K mo del par ameter up- dates using a single data sample. The ov er all computational times for all mode ls a nd MLP s are rep orted in T able 6 . All algo rithms were implemen ted in Scila b [37] version 5.5.2 and executed on a 2.4 Ghz Intel Xeon(R) v3 pro cessor running the Linu x op erating s ystem. Summarizing the 31 Performance Indice s N5 N6 MLP1 Mean (V ariance) L MLP2 Mean (V ariance) PVGM CCD Sim. Time (s ) 7 . 87 7 . 86 127 (285) 7 . 58 130 (433) 9 . 58 10 . 22 T able 6: Comput ational times. Results are group ed by complexity and purpose in order to simpli fy the comparison: 5 parameters (N5) in red column, 6 parameters (N6 and MLP1) in blue columns, 11 parameters (L and MLP 2) in green column, and benc hmark mo dels (ODNP , PVGM and CCD ) in gray columns. discussion ab ov e, we found that the prop osed parametr ic mo dels show co mparable foreca sting b e- havior and a g eneral per formance improv ement with r espect to MLPs, b oth from an accura cy and a computationa l vie wp oint. In par ticular, mo del N6 exhibits the b est compr o mise b etw e e n per - formance, complexity and co n vergence/consistency prop erties. This appar en tly shows that a slight ov erparameter ization with resp ect to a minimal mo del inspired by the theoretical mo dels of CCF and power gener a tion, turns out to b e b eneficial. On the co ntrary , excessive ov erpa rameterization impacts mo del ide ntifiabilit y and consistency . 9. Conclusions In this paper , we hav e prop osed an efficient pa rametric technique a imed at mo del estimation for dir ect for ecasting of PV p ower genera tion using cloud c over data. The appr oach is ba sed on r ecursive least squar es and Extended Ka lma n Filter. T o this purp ose, three mo dels have b een singled out: an 11-parameter linear model, and tw o nonlinear models, inv olving 5 and 6 parameters, resp ectiv ely . The pr ocedur e is esp ecially fit for the t ypical scena rio where the netw o rk op erator , due to the large n umber o f mana ged producers , has no access to on-site irradiance and tempera tur e measurements. The metho d exploits only the historical time s e r ies of generated p ow er, cloud cov er, and foreca st temper a ture, which ca n b e obtained from a meteor ological service . The proce dure inv olves mo dest memory and computational requirements. Its p erformance has been ev a luated on simulated data as well a s o n a single plant lo cated in Italy . The exp erimental results show a substantial impr ovemen t in terms of typical er ror measures with resp ect to alter nativ e approa c hes based on Neural Netw orks and standard linear autor egressive mo dels. All the prop osed mo dels show equiv alent fo r ecasting p erformance, while the 6-parameter nonlinear mo del gener ally outp erforms the other tw o as far a s data fitting prop erties ar e co ncerned. It is expec ted that the prop osed approach can b e fruitfully a da pted to the aggr e gation o f sev eral plan ts o ver macr o ar eas, with 32 significant per formance improv ements with resp ect to the single pla n t case. This issue is the sub ject of cur ren t inv estigation. References References [1] P . Denholm, R. M. Margolis, Ev aluating the limits of so lar pho to voltaics (pv) in traditional electric power systems, Energy Policy 35 (5) (2007 ) 2 852 – 28 61. doi:ht tp://dx. doi.org/10.1016/j.enpol.2006.10.014 . [2] P . Denholm, R. M. Margolis , E v aluating the limits of sola r photovoltaics (pv) in electric pow er systems utilizing energy stora g e and other ena bling technologies, E nergy Policy 35 (9) (2007 ) 4424 – 443 3. doi:h ttp://dx .doi.org/10.1016/j.enpol.2007.03.004 . [3] F. Albuyeh, Grid of the future, IE EE Po wer and Ener gy Mag azine 7 (2 ) (20 09) 5 2 –62. [4] C. Ba ldi, F. Cor ti, G. Di Lembo, F. Nebiaco lom b o, Monito ring a nd control o f a ctiv e distribu- tion grid, in: P ro c.CIRED W ork s hop 20 12, Lis bon, 201 2. [5] M. Diagne, M. David, P . Lauret, J. Boland, N. Schm utz, Review of solar irradiance forecasting metho ds and a prop osition fo r small-scale ins ula r grids, Renewable and Sustaina ble Energ y Reviews 27 (201 3) 65 –76. [6] C. F. Coimbra, J. Kleissl, R. Marquez, Overview of solar-for ecasting metho ds a nd a metric for accuracy ev alua tion, in: J. Kleissl (E d.), Solar E ner gy F o recasting and Reso urce Asses s men t, Elsevier Inc., W altham, Massa c husetts, USA, 2013 , Ch. 8 . [7] R. H. Inman, H. T. Pedro, C. F. Coim bra, Solar forecasting methods for renewable en- ergy in tegra tion, Progress in Energy and Combu stion Scienc e 39 (6) (2 0 13) 535 – 57 6. doi:ht tps://do i.org/10.1016/j.pecs.2013.06.002 . [8] G. Capizzi, C. Nap oli, F. Bonanno, Innov ative Second-Generation W av elets Construction With Recurrent Neura l Net works for Solar Radiation F orecasting, IE EE T ransactio ns o n Neural Net works and Lea r ning Systems 2 3 (11) (2012 ) 1 8 05–181 5 . 33 [9] J. W u, C. Chan, Pr ediction of hourly so lar radiation using a novel hybrid mo del, Sola r Energy 85 (2011 ) 8 08–817. [10] D. T ao , C. Shanxu, C. Changs ong, F orecasting p ow er output fo r grid-connected photov oltaic power s y stem without using solar radiation measurement, in: Pro c. 2nd IEE E International Sympo sium on Po wer Electronics for Distributed Gener ation Systems, 2 012, pp. 7 73–777. [11] C. Co rnaro, F. B ucci, M. Pierro, F. Del F rate, S. P ero naci, A. T arav at, So lar radia tion fore- cast using neural net works for the prediction of gr id connected pv pla nts energy pr oduction (dsp pro ject), in: P roc eedings of 28th Europ ean Photovoltaic Solar E ne r gy Conference and Exhibition, 20 13, pp. 39 92–399 9. [12] A. Ragnacc i, M. Pastorelli, P . V alig i, E . Ricci, Explo iting dimensionality r eduction techniques for photovoltaic power forecasting, in: Pro c. 2nd IEEE Ene r gyCon, Florence, Italy , 2012 . [13] G. Reik ard, Pre dicting s olar radia tion at high r esolutions: A comparison o f time series for ecasts, Solar Ener gy 8 3 (200 9) 34 2–349. [14] E. Lorenz, J. Hurk a, D. Heinemann, H. G. Beyer, Ir radiance F or ecasting fo r the Po wer P re- diction of Grid-Connected Photovoltaic Systems, IEEE Jo ur nal of selected topics in applied earth o bs erv ations and r emote sensing 2 (1) (200 9) 2– 10. [15] P . Bacher, H. Madsena, H. Nielsen, Online shor t-term solar p ow er forec asting, Sola r Energy 83 (2009 ) 1 772–178 3. [16] L. W ong, W. Chow, Sola r Radia tion Mo del, Applied Energy 6 9 (2001) 19 1–224. [17] C. A. Gueymard, D. Thevenard, Monthly average clear -sky bro a dband irr adiance database for worldwide s olar heat g ain and building c o oling lo a d ca lculations, Solar Energy 83 (2 0 09) 1998– 2018. [18] M. Patel, Wind a nd Solar P ow er Systems: Design, Analysis and Oper ation, T aylor and F r ancis, Bo ca Ra ton, FL, 2 006. [19] G. Bianchini, S. Paoletti, A. Vicino, F. Corti, F. Nebiaco lom b o, Mo del estimation of photo- voltaic power genera tion using par tia l information, in: Innov ative Smart Grid T e c hnologies Europ e (ISGT EUROPE), 2 013 4th IEE E/PES, IE EE, 2013, pp. 1– 5. 34 [20] G. Bianchini, S. Paoletti, A. Vicino, F. C o rti, F. Nebiacolo m b o , Estima tio n of a simple mo del of sola r power generation using partial information, in: 52 nd IEEE Conference on Decisio n and Control, 2013, pp. 99 6–1000. do i:10.110 9/CDC.2013.6760012 . [21] R. Do ws, E. Gough, P VUSA pr o curement, acce pta nce, and rating pra ctices for photovoltaic power plants, T ech. rep., Pacific Ga s and Electr ic Company , San Ramon, CA (19 95). [22] A. Hammer, D. Heinema nn, E. Lorenz , B. Lue ckehe, Short-ter m for e casting of solar radi- ation: a statis tica l appr o ac h using satellite data, Sola r Energ y 67 (1) (19 99) 139 – 150. doi:ht tp://dx. doi.org/10.1016/S0038- 092X(00)00038- 4 . [23] R. Perez, K . Mo ore , S. Wilcox, D. Renn ´ e, A. Zelenk a, F or ecasting s olar radiation - preliminar y ev aluation of an appr oach based upo n the national forecast database, So la r Ener gy 81 (6) (2007) 80 9 – 812. doi:htt p://dx.d oi.org/10.1016/j.solener.2006.09.009 . [24] R. Perez, S. K iv alov, J. Schlemmer, K. Hemk er, D. Renn´ e, T. E. Hoff, V alidatio n of short and medium term op erationa l s olar radiation forecas ts in the us , Solar Energ y 84 (1 2 ) (2 010) 2161– 2172. [25] B. O. Kang , K.-S. T am, New and improved metho ds to estimate day-ahead quantit y and quality o f sola r irr adiance, Applied Ener g y 137 (201 5) 240 – 249. doi:ht tp://dx. doi.org/10.1016/j.apenergy.2014.10.021 . [26] F. Ka sten, G. Czepla k , So la r and terr estrial radia tion dep enden t on the amount a nd type of cloud, Sola r energ y 24 (2) (19 80) 177 –189. [27] K. Kimura, D. Stephenson, Solar radiation on cloudy da ys, ASHRAE transa ctions 75 (1) (1969) 1–8. [28] F. Lewis, L. Xie, D. Popa, Optimal and robust estimation with an Intro duction to Sto c hastic Control Theor y , CRC P ress, 2007 . [29] D. Pepe, G. B ianc hini, A. Vicino, Mo del estimation of photov oltaic power gener ation using partial informatio n, in: IEE E Energy Co n 2016, IE EE, 2 016. [30] Wikip edia, O kta — Wikip edia, The Fr ee Encyclo pedia, [Online; accessed 10-Novem b er-201 5] (2015). URL ht tps://en .wikipedia.org/w/index.php?title=Okta$&$oldid=687869646 35 [31] H.-C. Kim, E. E . Hofmann, Ev aluation and deriv ation of cloud-cover algorithms for calculation of surface irra dia nce in sub-antarctic and a n tarctic environmen ts, Antarctic Science 1 7 (01) (2005) 13 5–150. [32] A. B. Meinel, M. P . Me ine l, Applied solar energy . An intro duction, Addison-W esley Publishing Co., Reading , MA, 1976. [33] A. Luque, S. Hege dus , Handb o ok of photovoltaic s c ience and engineering, J ohn Wiley & So ns , 2011. [34] A. Sp ena, G. D’Angiolini, C. Strati, First cor relations for solar radia tion o n clo udy days in italy , in: ASME-A TI-UIT 2 010 Conference o n Thermal a nd E nvironmental Issues in Energy Systems, Sor ren to, Italy , 20 10. [35] International Energy Age ncy, Photovoltaic a nd solar foreca sting: State o f the ar t, Photovoltaic Po wer Systems Pr ogramme Rep ort IEA P VPS T14 0 1, 2013 . [36] J. Wid ´ en, N. Carpman, V. Castellucci, D. Lingfors , J. O lauson, F. Remouit, M. Bergkv is t, M. Gra bbe, R. W ater s, V ariability assess men t and forecasting of renewables: A review for solar, wind, wa ve a nd tidal res o urces, Renewable and Sus tainable Energy Reviews 4 4 (201 5) 356 – 375. do i:http:/ /dx.doi.org/10.1016/j.rser.2014.12.019 . [37] Scilab. Op en sour ce softw are for numerical computation. URL ht tp://www .scilab.org 36