- Email: [email protected]
Tilt to horizontal global solar irradiance conversion: An evaluation at high tilt angles and different orientations
Renewable Energy 113 (2017) 1529e1538
Contents lists available at ScienceDirect
Renewable Energy journal homepage: www.elsevier.com/locate/renene
Tilt to horizontal global solar irradiance conversion: An evaluation at high tilt angles and different orientations dric Bertrand Caroline Housmans*, Alessandro Ipe, Ce Royal Meteorological Institute of Belgium, Avenue Circulaire 3, B-1180, Brussels, Belgium
a r t i c l e i n f o
a b s t r a c t
Article history: Received 5 January 2017 Received in revised form 29 June 2017 Accepted 3 July 2017 Available online 5 July 2017
Many transposition models have been proposed in the literature to convert solar irradiance on the horizontal plane to that on a tilted plane. The inverse process, i.e., the conversion from tilted to horizontal is investigated here based on six months of in-plane global solar irradiance measurements recorded on the roof of the Royal Meteorological Institute of Belgium's radiation tower in Uccle (Latitude 50.79 N, Longitude 4.35 E). Up to three pyranometers mounted on inclined planes of different tilts and orientations were involved in the inverse transposition process. Our results indicate that (1) the tilt to horizontal irradiance conversion is improved by using a multi-pyranometers approach and (2) the improvement from using an isotropic model to anisotropic models in the inverse transposition problem is not significant. © 2017 Elsevier Ltd. All rights reserved.
Keywords: Solar radiation Transposition model Decomposition model
1. Introduction Knowledge of local solar irradiation is essential for many applications, including architectural design, solar energy systems, crop growth models, and evapotranspiration estimates (e.g., Wali et al. [1]; Trnka et al. [2]). The usual solar radiation parameters measured on ground are the global horizontal irradiance, Gh , the direct normal irradiance, Bn , the diffuse horizontal irradiance, Dh , and the sunshine duration. Traditionally, they are recorded by means of networks of meteorological stations. However, costs for installation and maintenance of such networks are very high and national networks comprise only few stations. For instance, the solar measurements network operated by the Royal Meteorological Institute of Belgium (RMI) has only 15 sites where various combinations of the solar irradiation components are recorded (e.g., e and Bertrand [3]). Consequently the availability of Journe observed solar radiation measurements has proven to be spatially and temporally inadequate for many applications. Over the last decades, satellite-based retrieval of solar radiation at ground level has proven to be valuable for delivering a global coverage of the global solar irradiance distribution at the Earth's et al. [7]). surface (e.g., Perez et al. [4]- [5]; Hammer et al. [6]; Renne The recent deployments of solar photovoltaic (PV) systems offer a
* Corresponding author. E-mail address: [email protected] (C. Housmans). http://dx.doi.org/10.1016/j.renene.2017.07.018 0960-1481/© 2017 Elsevier Ltd. All rights reserved.
potential opportunity of providing additional solar resource information as PV systems with global tilted irradiance, Gt , measurements become available in some locations. However, in most PV systems, modules are installed on a fixed plane to reduce installation and operation costs compared to tracking systems. To maximize array output, modules are installed at a tilt close to local latitude, or at some minimum tilt to ensure self-cleaning by rain. Gh estimates at the PV system location require therefore to convert the tilted solar irradiance recorded to that on a horizontal plane. If many transposition models have been proposed in the literature (see Yang [8] for a review) to convert solar irradiance on the horizontal plane to that on a tilted plane, the inverse process (i.e., converting from tilted to horizontal) is only poorly discussed in literature. The difficulty relies on the fact the procedure is analytically not invertible. As an example, single-sensor approach involving a numerical search method was proposed by Yang et al. [9] and evaluated on various combinations of decomposition (i.e., models that separate direct and diffuse solar components from the global one) and transposition models. They found their method sufficiently accurate for small zenith angles, but reported that the conversion error exponentially increases as zenith angle increases. More recently, Marion [10] presented an iterative method using a modified version of the DIRINT decomposition model (Perez et al. [11]) in combination with the transposition model of Perez et al. [5] to obtain horizontal irradiance from a single tilted sensor measurement. Performance of this method was found essentially the
1530
C. Housmans et al. / Renewable Energy 113 (2017) 1529e1538
same as for the DIRINT model used in the normal calculation mode for Gt measurements recorded at small tilt angle from the horizontal. At larger angles, deviations between modeled and measured direct and diffuse irradiance values increase as the tilt angle increases from the horizontal. Lastly, an alternative method based on the Olmo model (Olmo et al. [12]) that presents the particularly of being analytically invertible was proposed by Killinger et al. [13]. If the overall performance of the inverted Olmo model was found comparable with the other approaches, the results were slightly worse than those obtained by inverting the decomposition and transposition models in combination with an iterative solving process. The present study aims at (1) evaluate the tilt to horizontal irradiance conversion when using for model input Gt measurements recorded in angular configurations similar to those encountered in Belgian PV systems installations (e.g., with tilt angle as great as 50.79 ) and (2) determine if as suggested by Yang et al. [14] the use of a multi-pyranometers system could improve the conversion accuracy. The selection of the decomposition and transposition models used in our calculations is based on previous evaluation of popular decomposition and transposition models performance in Uccle by Demain et al. [15]- [16] and Bertrand et al. [17]. The paper is organized as follows: methods to perform the conversion from tilt to horizontal are presented in section 2. In situ measurements are briefly described in section 3. Performances of the different approaches are evaluated in section 4. Final remarks and conclusions are provided in section 5.
Transposition models have the general form:
(1)
where the tilted global solar irradiance, Gt , is expressed as the sum of the in-plane direct irradiance, Bt , in-plane diffuse irradiance, Dt , and the irradiance due to the ground reflection, Dg . The direct component, Bt , is obtained from:
Bt ¼ Bn cosqi ¼ Bh
cosqi ¼ Bh rb cosqz
Dt ¼ Dh Rd
(3)
Dg ¼ rGh Rr
(4)
where, Dh , is the diffuse horizontal irradiance, Gh , the global horizontal irradiance (i.e., Gh ¼ Dh þ Bh ), Rd , the diffuse transposition factor and r, the ground albedo. The transposition factor for ground reflection, Rr , can be modeled under the isotropic assumption (e.g., Gueymard [18]) as follows:
Riso r ¼
1 cosb 2
(2)
(5)
where, b, is the tilt angle of the inclined surface. See Fig. 1 for angles definition. 2.1. Single-sensor approach Considering the effective global horizontal transmittance, Kt , the direct normal transmittance, Kn , and the diffuse horizontal transmittance, Kd (i.e., Kd þ Kn ¼ Kt ):
Gh ¼ Kt Io cos qz Bn ¼ Kn Io Dh ¼ Kd Io cos qz
2. Conversion from tilt to horizontal
Gt ¼ Bt þ Dt þ Dg
with, Bn , the direct normal irradiance and, Bh , the direct irradiance on a horizontal surface, respectively. qi is the incidence angle and, qz , the solar zenith angle, respectively. Parameter rb ¼ cos qi =cos qz is a factor that accounts for the direction of the beam radiation. The diffuse component, Dt , and the irradiance due to the ground reflection, Dg , can be modeled as follows:
(6)
(where Io is the extraterrestrial normal incident irradiance), Eq. (1) can be rewritten as:
K K Gt ¼ Kt Io cos qi 1 d þ cos qz d Rd þ rRr Kt Kt
(7)
Eq. (7) indicates that the conversion of Gt to Gh requires the use of a decomposition model to estimate Kd from Kt . The decomposition model of Skartveit and Olseth [19] (hereafter referred to as OLS) was used in the present study. Formulation of this model is
Fig. 1. Angles to define the position of the sun and the orientation of a tilted plane: qz is the solar zenith angle, g, the solar elevation, qi , the incidence angle, b, the tilt angle of the inclined surface and the azimuth of the tilted plane is a.
C. Housmans et al. / Renewable Energy 113 (2017) 1529e1538
provided in Appendix A. It allows to estimate the direct normal irradiance, Bn , from the global solar irradiance, Gh , through an empirical equation depending on Kt . Performance of Eq. (7) has been evaluated for various geometric configurations using the isotropic transposition model proposed by Liu and Jordan [20] (hereafter referred to as LIU model) and the anisotropic models of Hay [21] (hereafter referred to as HAY model), Skartveit and Olseth [22] (hereafter referred to as SKA model) and Perez et al. [23] (hereafter referred to as PER model), respectively. Formulation of the selected transposition models is provided in Appendix B. Note that because the PER model coefficients are binned according to the sky clearness parameter, ε, (see Appendix (B.4)), the same strategy as in Yang et al. [14] was adopted here when running the PER model. Basi0 0 cally, Kt and Kd estimates are calculated for each ε bins using the associated model coefficients Fij . In each ε bins Gh 0 and Dh 0 esti0 0 mates are then retrieved from the corresponding Kt and Kd esti0 mates and used to determine the corresponding ε value. If the 0 0 assumed ε agrees with the calculated ε , the corresponding Gh and 0 Dh are selected as the true estimated Gh and Dh . Similarly, because the transposition factor defined by the SKA model depends on the Skartveit-Olseth's correction factor Z (see Eq. (24)), the conversion process is performed for the two Rd model formulations before selecting the appropriate solution. The global horizontal irradiance found has to satisfy the Z value definition assumed during the conversion. We solved Eq. (7) by a simple iteration procedure. Namely, Kt is adjusted using the bisection procedure until the Kt provides calculated Gt value that match the measured Gt (using as starting point Kt ¼ 0:5). Note that for the sake of simplicity, the foreground's albedo, r, was taken to 0.23 (Demain et al. [15]) in our calculations.
l¼
1531
Gt2 B1 Gt1 B2 Gt1 A2 Gt2 A1
(12)
Finally, substituting Eq. (12) into Eq. (10) leads to a first degree equation with one unknown: 0
A Dh þ B0 ¼ 0
(13)
where. 0
A ¼ B1 ðGt2 A1 Gt1 A2 Þ 0 B ¼ Gt1 ðGt1 A2 þ A1 B2 Þ Gt2 A1 B1 Gt1 Gt2 A2 Considering the anisotropic approach to the inverse transposition problem, Eq. (8) can be rewritten as:
Gt ¼ ABh Dh þ BBh þ CDh
(14)
where the formulation of the coefficients A, B, and C, depends on the considered transposition model. Following the same reasoning as for the isotropic approach, leads in the anisotropic case to solve a quadratic equation: 00
00
00
A D2h þ B Dh þ C ¼ 0
(15) 00
00
00
where the coefficients A , B and C depends on the considered transposition model. Given n tilted pyranometers (with different inclinations and/or orientations), Eq. (14) can be generalized as the matrix form (e.g., Yang et al. [14]):
(16)
2.2. Multi-sensors approach The idea that simultaneous readings of a multi-pyranometers system can be used to disangle the various components of solar radiation on inclined surfaces has been originally proposed to solve in remote locations the periodic adjustment required by normal incidence and shadow-band pyranometers to ensure that their readings remain accurate when long-term data acquisition is in progress (Faiman et al. [24]). By combining Eqs. (2)e(4) we rewrite Eq. (1) as:
Gt ¼ Bh rb þ Dh Rd þ rGh Rr
(8)
Assuming first Rd as isotropic, Eq. (8) becomes:
Gt ¼ Dh ðlA þ BÞ
(9)
where.
where L ¼ fAi g is a 2 n 2 third-order tensor, B ¼ fBi g is a n 2 matrix, ℭ is a n column vector, and x a column vector with 2 variables:
0 Ai 2<2n2 0 Ai 0 1 C1 B1 B C2 B2 C C2
(17)
Note that Eqs. (16) and (17) slightly differ from the formulations given in Yang et al. [14] due to the change in the column vector x, i.e., (Dh Bh ) instead of (Dh l). The least square (hereafter referred to as LS) solution to Eq. (16) is given by:
cos qi þ rRr cos qz 1þcos b þ rRr 2
A¼ B¼ l ¼ Bh =Dh
Supposing a two-pyranometers system with Gt1 and Gt2 the corresponding tilted global irradiance measurements, it follows from Eq. (9) that:
Gt1 ¼ Dh ðlA1 þ B1 Þ
(10)
Gt2 ¼ Dh ðlA2 þ B2 Þ
(11)
Dividing Eq. (10) by Eq. (11) to eliminate Dh yields, after rearranging:
(18) with jj:jj referring to the Euclidean norm. However, the LS is hard to solve and a standard technique to resolve Eq. (18) is to use a Newton type iteration method (e.g., Grosan and Abraham [25]). As an alternative, Eq. (16) can also be solved by minimizing the errors (this approach is hereafter denoted to as EM - Errors Minimization). In this case, the solution is to minimize
1532
C. Housmans et al. / Renewable Energy 113 (2017) 1529e1538
Table 1 Data availability percentage according to the sky condition for the considered period from July 15, 2015 to January 19, 2016. Sky condition
Kt
Data availability (%)
All sky Overcast Cloudy Partly cloudy Partly clear Clear
0.0e1.0 0.0e0.2 0.2e0.4 0.4e0.6 0.6e0.8 0.8e1.0
100% 32.07% 26.92% 18.01% 21.46% 1.54%
( min EðxÞ ¼
n X
) ε2i ðxÞ : x2<2
(19)
i¼1
where, εi ðxÞ ¼ ðAi Dh Bh þ Bi Bh þ Ci Dh Þ Gti , with i ¼ 1; …; n denoting the tilted pyranometer (n 2). Both ways to solve Eq. (16) have been considered and compared here using the Powell's quadratically convergent method (Powell [26]) to perform the minimization. It is a generic minimization method that allows to minimize a quadratic function of several variables without calculating derivatives. The key advantage of not requiring explicit solution of derivatives is the very fast execution time of the Powell method. In order to avoid the problem of linear dependence in Powell's algorithm, we adopted the modified Powell's method given in Acton [27] and implemented in Press et al. [28]. 3. In situ measurements The present analysis is based upon six months of global solar irradiance measurements performed on the roof of the radiometric tower of the Royal Meteorological Institute of Belgium located on the Brussels Uccle plateau (Latitude 50.79 N, Longitude 4.35 E at an elevation of 101 m above sea level). Belgian climate, like most of the Northern Europe (Peel et al. [29]) is maritime temperate with significant precipitation in all seasons. Cloud occurrence is high (see Table 1) and clouds are the dominant
source of small-scale variability in surface solar radiation. Fig. 2 compares the global solar irradiation (in kWh=m2 ) received in Uccle during the six months of the measurements campaign (black curve) with the corresponding climatological normals (i.e. the 1981e2010 average, in green). The highest and lowest values recorded since 1981 are also provided for illustration (in red and blue, respectively). Four data sets of global solar irradiance have been collected from July 15, 2015 to January 19, 2016. The first one was recorded by a Kipp & Zonen CM11 Secondary Standard pyranometer mounted on the horizontal plane (here after referred to as Plane 0). The second was also recorded by a Kipp & Zonen CM11 Secondary Standard pyranometer but mounted on a tilted plane of 50.79 (i.e., corresponding to the RMI radiometric tower's latitude) facing south (here after referred to as Plane 1). For the last two data sets, two additional Kipp & Zonen CMP-22 Secondary Standard pyranometers were installed on the tower with the same angular configurations than two neighboring residential PV installations, i.e., a tilted plane of 45 SW facing and a tilted plane of 50 facing E (here after referred to as Plane 2 and Plane 3, respectively). Regarding the accuracy of the measurements devices (see Fig. 3 for instruments illustration), the directional error (up to 80 with 1000 Wm2 beam) is less than 5 Wm2 for the CMP22 pyranometer and less than 10 Wm2 for the CM11 pyranometer, respectively. The pyranometer spectral selectivity (300e1500 nm) is smaller than 3% for the CMP22 instrument and smaller than 2% for the CM11 instrument, respectively. Both instruments have a thermopile detector and present a reduced thermal offset (i.e., about 2 Wm2 at 5 K/h temperature change) for the CM11 pyranometer and less than 1 Wm2 for the CMP22 pyranometer, respectively. Irradiance measurements were made with a 5-sec time step and then integrated to bring them to a 10-min time step. Note that because of the difference in the recording planes orientations and inclinations, the number of available measurements differs from one instrument to another (see Table 2). All 10-min data have undergone a series of automated quality control procedures similar to e and Bertrand [3] prior to be visually those described in Journe inspected by a human operator.
Fig. 2. Comparison between the global solar irradiation (in kWh=m2 ) received in Uccle during the six months (i.e., from July 15, 2015 to January 19, 2016) of the measurements campaign (black curve) with the corresponding climatological normals (i.e., the 1981e2010 average, in green). The highest and lowest values recorded since 1981 are also provided for illustration (in red and blue, respectively). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
C. Housmans et al. / Renewable Energy 113 (2017) 1529e1538
1533
Fig. 3. Panel A presents an illustration of the Kipp & Zonen CM11 Secondary Standard pyranometer instrument used to perform Plane 0 and Plane 1 global solar irradiance measurements. Similarly panel B shows the Kipp & Zonen CMP-22 Secondary Standard pyranometer instrument considered to perform Planes 2 and 3 global solar irradiance measurements. Note that in panel A, the instrument is illustrated in the mounting Plane 0 angular configuration on the roof of the RMI's radiometric tower while in panel B, the instrument is given in the mounting Plane 3 angular configuration.
Table 2 Angular configuration of the four recording planes considered in the study. Note that the surface azimuthal angle, a, is conventionally measured clockwise from the South. Also provided is the number of available data points for each of the different planes. Note that the total number of available simultaneous data points in all of the recording planes is 5816 (i.e., 100%). Pyranometer configuration
Azimuth Angle a [ ]
Tilt Angle b [ ]
Data Points Number
Plane Plane Plane Plane
/ 0.00 (S facing) 45.00 (SW facing) 90.00 (E facing)
0.00 50.79 45.00 50.00
10700 10578 9631 7269
0 1 2 3
4. Results The relative ability of the different approaches to predict horizontal irradiance from the tilted one was estimated by means of two statistical error indexes: Mean Bias Error (MBE) and Root Mean Square Error (RMSE).
MBE ¼
n 1X ðe Þ n i¼1 i
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u n u1 X e2 RMSE ¼ t n i¼1 i where ei ¼ ðGi;e Gi;o Þ is the residual value; Gi;e are the estimated values and Gi;o represent the observed measurements. A positive MBE (resp. a negative MBE) means that the model tends to overestimate (resp. underestimate) the observed measurements. To obtain dimensionless statistical indicators we expressed MBE and RMSE as fractions of mean solar global irradiance during the respective time interval,
MBE½% ¼
MBE M
RMSE½% ¼
RMSE M
P where M ¼ 1n ni¼1 ðGi;o Þ is the measurements mean. Note that for a proper estimation of the error statistics, only radiation data recorded with a solar zenithal angle, qz , smaller than 85 and an incidence angle on the plane used, qi , smaller than 90 were considered and tilted global solar irradiance records were further restricted to
non-zero values. Moreover when the retrieved global horizontal solar irradiance was negative or larger than the corresponding extraterrestrial irradiance the inverse transposition problem was considered unsuccessful. With these additional conditions the total number of data points available simultaneously in all of the recording planes reduces to 5816 (i.e. 100%). 4.1. Single-sensor approach Table 3 summarizes for each of the three tilted planes the performance of the inverted transposition models in the Gt to Gh conversion. To evaluate the transposition models on a same basis, only the 5117 data points conversions that have succeeded for all models have been considered in the MBE and RMSE calculations (e.g., 5117/5816 ¼ 87.98% of the total available data points, 5816 being the total number of measurements available simultaneously for the four planes). It is obvious from Table 3 that none of the considered transposition models performs best in terms of failure rate, RMSE and MBE over the three tilted planes. Moreover, the tilt angle and the surface's orientation have a major impact on the Gh estimation's reliability irrespective of the considered transposition model. The worst results in term of RMSE are for Plane 2 measurements conversion (i.e., rRMSE ranging from 28.7% to 30.8% vs 6.8%e12.2% reported for Plane 1 conversions). Furthermore, Plane 2 conversions underestimate Gh (i.e., rMBE of 14.1% to 11.7%) while a slight overestimation of 0.3%e5.3% is reported for Plane 1 and an overestimation of 2.3%e7.6% for Plane 3, respectively. Regarding the failure rate, the situation varies according to the considered transposition model. As an illustration, while almost all tilted irradiance measurements from Plane 1 are successfully converted when performing the inverse transposition with the HAY and SKA models (i.e., failure rate of 0.01%), the largest percentage of unsuccessfull conversions reported for the LIU and PER models is for the Plane 1 measurements conversion (i.e., a failure rate of 1.6% and 4%, respectively). For comparison, Marion (2015) reported a failure rate of 10% for its iterative method combining the DIRINT decomposition model with the PER transposition model. Yang et al. [9] reported rRMSE values of roughly 13% and 6% for the tilt to horizontal conversion of irradiance data collected in the tropical region of Singapore by sensors mounted on a 18.3 tilted NE facing plane and 6.1 SW facing plane, respectively. The overall performance of our Gt to Gh conversion appears relatively lower especially for Planes 2 and 3 measurements conversion. However, it is worth pointing out that (1) the tilt angles considered in Yang et al. [9] were significantly less than here and (2) Yang et al. [14] reported that the RMSE resulting from their conversion using a single tilted
1534
C. Housmans et al. / Renewable Energy 113 (2017) 1529e1538
resolution of Eq. (13) (i.e., a first degree equation with one unknown). Validation statistics were computed from the 1499 data points that succeeded the inverse transposition with all the considered transposition models and the two solving methods (e.g., 25.77% of the total available data points). It is clear from Table 4 that both solving methods produce very distinct results in terms of RMSE and MBE as well as in term of failure rate. Clearly, the EM method allows to reduce the occurrence of invalid solutions in the inverse transposition problem. As an example, the failure rate associated to Plane 1 - Plane 2 conversions carried out by running the LIU (HAY) model reduces from 8.19% (33.23%) with the quadratic formula to 2.03% (2.01%) with the EM method. However, it is worth pointing out that if the failure rate associated to the PER model is found to generally decrease with the EM method, the number of invalid solutions reported for this model is very large. Because the PER model depends on ε binning, it can appear that the retrieved ε value using the procedure described in section 2.1 does not agree with the expected bin, or that more than one retrieved ε value agree with the expected bin, increasing therefore the number of unsuccessful conversions associated to this model. Using the EM method instead of the LS procedure does not solve the ε binning formulation deficiencies; in contrary the occurrence of the problem is even enlarged with the EM method for the Plane 1 þ Plane 3 conversions (failure rate of about 37% with the EM method compared to 21% with the LS solution). By contrast, in depth-analysis of the conversion process indicated that for a number of data points, conversion performed with the HAY and SKA models using the EM method allows to produce a valid result when the discriminant or the roots of the quadratic equation are found negative or when two positive real solutions are found (as in this latter case, no objective criterion permitted us to select the appropriate solution and the conversion was considered unsuccesfull). Regarding the LIU isotropic transposition model, even if the inverse transposition reduces to a first degree equation with one unknown, use of the Powell's minimization algorithm tends to overdetermine the system allowing therefore to provide a solution (e.g., the best fitting solution) when no valid solution was found by the LS procedure. As a comparison, Yang et al. [14] reported a failure rate of about 2e3% for a two tilted pyranometers system solved by least squares method (LS) using the PER transposition model and
Table 3 Absolute (in W m2) and relative (in %) RMSE and MBE indexes for the single-sensor approach (one pyranometer). Results for the three pyranometer's mounting planes are provided for the four considered transposition models (see Appendix B) coupled with the OLS decomposition model (see Eq. (20)). Validation statistics are computed from 5117 data points (e.g., 87.98% of the total available data points). Transposition þ Decomposition Models
Pyranometer's mounting plane Failure rate (%) Plane 1 Plane 2 Plane 3 RMSE (W m2 ) Plane 1 Plane 2 Plane 3 rRMSE (%) Plane 1 Plane 2 Plane 3 MBE (W m2 ) Plane 1 Plane 2 Plane 3 rMBE (%) Plane 1 Plane 2 Plane 3
LIU þ OLS
HAY þ OLS
SKA þ OLS
PER þ OLS
1.61 0.27 0.32
0.01 0.18 2.43
0.01 0.18 2.43
4.03 2.81 3.08
33.88 85.83 59.14
19.03 85.76 65.29
19.20 83.61 63.54
25.44 80.16 61.70
12.15 30.78 21.21
6.82 30.75 23.41
6.89 29.98 22.78
9.12 28.74 22.13
14.85 37.45 6.46
0.90 39.18 10.55
4.15 35.51 14.71
5.63 32.72 21.27
5.32 13.43 2.32
0.32 14.05 3.78
1.49 12.73 5.28
2.02 11.73 7.63
sensor NE facing increased from 4% for a tilt angle of 10 e23% for a tilt angle of 40 . 4.2. Two tilted pyranometers Performance of the inverse transposition using two different tilted global irradiance measurements as input to the Gh computation are presented in Table 4 for the 1499 data points that have succeeded for all models as well as for both the EM method and the LS procedure. For each transposition model, Table 4 provides results obtained for both the quadratic formula (in brackets) and the EM method (i.e., Eq. (19)) to solve the quadratic equation. Note that for the LIU isotropic model, values in brackets are from the direct
Table 4 Absolute (in W m2) and relative (in %) RMSE and MBE indexes for a two pyranometers system. Results are provided for the four considered transposition models (see Appendix B). Validation statistics are computed from 1499 data points (e.g., 25.77% of the total available data points) for both the EM method and the LS procedure. Results of the LS procedure are given in brackets. Pyranometer's System
LIU
HAY
SKA
PER
2.03 (8.19) 0.74 (2.08) 0.67 (8.83)
2.01 (33.23) 0.74 (5.30) 0.67 (21.59)
3.34 (34.30) 3.90 (5.35) 0.68 (22.66)
30.88 (44.45) 36.97 (20.92) 26.44 (40.64)
2 3 3
13.89 (34.09) 11.04 (20.09) 8.16 (19.92)
13.89 (29.25) 11.04 (17.73) 8.16 (16.34)
12.51 (44.19) 10.08 (21.34) 6.44 (24.94)
13.20 (60.69) 9.49 (38.04) 7.09 (44.04)
2 3 3
12.14 (29.77) 9.65 (17.55) 7.13 (17.40)
12.14 (25.55) 9.65 (15.49) 7.13 (14.28)
10.92 (38.60) 8.80 (18.64) 5.62 (21.78)
11.53 (53.01) 8.29 (33.22) 6.19 (38.47)
2 3 3
7.80 (1.91) 5.90 (0.24) 5.35 (0.53)
7.80 (0.67) 5.90 (-0.45) 5.35 (-0.39)
3.11 (5.67) 0.06 (6.23) 0.05 (6.75)
4.47 (14.79) 2.09 (14.33) 1.04 (15.39)
2 3 3
6.81 (1.67) 5.15 (0.21) 4.67 (0.46)
6.81 (0.58) 5.15 (-0.40) 4.67 (-0.34)
2.72 (4.95) 0.05 (5.45) 0.04 (5.89)
3.90 (12.92) 1.83 (12.51) 0.91 (13.44)
Failure rate (%) Plane 1 þ Plane 2 Plane 1 þ Plane 3 Plane 2 þ Plane 3 RMSE (W m2 ) Plane 1 þ Plane Plane 1 þ Plane Plane 2 þ Plane rRMSE (%) Plane 1 þ Plane Plane 1 þ Plane Plane 2 þ Plane MBE (W m2 ) Plane 1 þ Plane Plane 1 þ Plane Plane 2 þ Plane rMBE (%) Plane 1 þ Plane Plane 1 þ Plane Plane 2 þ Plane
Transposition Model
C. Housmans et al. / Renewable Energy 113 (2017) 1529e1538
1535
Table 5 Absolute (in W m2) and relative (in %) RMSE and MBE indexes for a three pyranometers system. Results are provided for the four considered transposition models (see Appendix B). Validation statistics are computed from 4066 data points (e.g., 69.91% of the total available data points) using the EM method. Transposition Model
Failure rate (%)
RMSE (Wm2)
rRMSE (%)
MBE (W m2)
rMBE (%)
LIU HAY SKA PER
0.0 0.0 0.0 30.09
20.03 20.06 19.55 21.56
7.93 7.94 7.74 8.53
9.68 9.70 5.94 9.04
3.83 3.84 2.35 3.58
found that the isotropic approach had guaranteed solution with their tropical data. Globally, use of the EM method instead of the LS procedure to solve the quadratic equation decreases by a factor of two and even more in some angular configurations the RMSE in the Gh estimation. As an example, the rRMSE reported for the PER model in the Plane 1 þ Plane 2 conversion is reduced from 53% with the LS solution to 11.5% with the EM approach. For the bias, the benefit of using the EM approach instead of the LS procedure is less apparent in the sense that in one hand the magnitude of the bias is reduced for the more elaborate SKA and PER anisotropic transposition models while in the other hand the bias magnitude is enlarged for both the LIU isotropic model and the simple HAY anisotropic transposition model. To this respect, it is interesting to note that the LIU and HAY models behave similarly in terms of failure rate, RMSE and MBE when using the EM method. In overall, the accuracy of the tilt to horizontal conversion is improved when two tilted irradiance measurements are involved in the inverse transposition process. However, the conversion dependence to the angular configuration of the tilted input measurements does not totally vanish and they are still some data points for which the procedure fails to produce a valid estimation of the global horizontal solar irradiance. 4.3. A three tilted pyranometers system Table 5 summarizes the tilt to horizontal conversion errors using a three tilted pyranometers system (as obtained by the EM solving method) for the 4066 data points that have succeeded for all models. Clearly, involving three different tilted irradiance measurements in the conversion process always ensures solutions excepted for the PER model for which the conversion of about one third (i.e., 30.1%) of the data points are still unsuccessful. Globally, all models behave quite similarly in term of RMSE (overall rRMSE of 8%) and present a negative bias (i.e., rMBE ranging from 3.8% to 2.6%). Models performance as a function of the sky conditions and the solar elevation angle, g, are summarized in terms of RMSE and MBE
in Tables 6 and 7, respectively. Table 6 indicates that the models relative accuracy does not vary noticeably as the sky condition move from overcast to clear sky situations (i.e., a rRMSE variation of about 4.6%, 3.5% and 3% is reported for the PER, SKA and the LIU and HAY models, respectively). All models show a negative bias in clear sky conditions (i.e., 0.8 Kt 1.0) where the global radiation is mainly composed of direct radiation. As an illustration, an underestimation of about 43.7 W m2 (6.5%) is reported for the LIU, HAY and SKA models and of 21.1 W m2 (3.2%) for the PER model. In overcast condition (i.e., 0.0 Kt < 0.2) where the diffuse component largely domines, while the global radiation is underestimated by about 3.7 W m2 (4.9%) by the LIU and HAY models, it is in the opposite slightly overestimated by the SKA and PER models (i.e., MBE of 0.5 W m2 or 0.7% and 0.4 W m2 or 0.5%, respectively). If the SKA model is the best performing model in overcast and cloudy conditions (i.e., Kt < 0.4) with reported RMSE values of 4.5 W m2 (6%) and 9.6 W m2 (5.2%), respectively, its performances do not differ from those reported for the LIU and HAY models in partly cloudy, partly clear and clear sky situations (i.e., Kt 0.4) with RMSE values ranging from 30.2 W m2 (5.8%) in partly clear conditions to 58.2 W m2 (8.7%) in clear sky conditions. The PER model exhibits the lowest RMSE scores for the different sky conditions excepted in overcast situation where its performance (i.e., RMSE of 4.6 W m2 or 6.1%) is closely similar to that reported for the SKA model. Unsurprisingly, models accuracy shows a dependence on the solar elevation with an overall RMSE decreasing from 19.9 W m2 (or 14.1%) in the low solar elevation bin (i.e., g 30 ) to 13.5 W m2 (4.4%) in the high solar elevation bin (i.e., g > 60+). While the SKA model exhibits the best performance in term of RMSE irrespective of the solar elevation, Table 7 indicates that the difference in models accuracy within a given solar elevation bin is limited (i.e., 0.65%, 0.84% and 1.49% in the low, middle and high solar elevation bins, respectively). The magnitude of the models bias tends to reduce as the solar elevation angle increase. As an example, the underestimation of the global solar reported in Table 7 for the LIU and HAY models reduce from about 10.1 W m2 (7.2%) in low solar elevation conditions to 8.6 W m2 (2.8%) in high solar elevation conditions. The sign of the bias is even modified in the case of the
Table 6 Absolute (in W m2) and relative (in %) RMSE and MBE indexes for a three pyranometers system as a function of the sky condition defined in Table 1. Results are provided for the four considered transposition models (see Appendix B). The relative indexes are displayed in brackets. Validation statistics are computed from 4066 data points (e.g., 69.91% of the total available data points) using the EM method. Transposition model
RMSE LIU HAY SKA PER MBE LIU HAY SKA PER
Sky condition 0.0Kt < 0.2
0.2Kt < 0.4
0.4Kt < 0.6
0.6Kt < 0.8
0.8Kt 1.0
6.08 6.08 4.48 4.61
11.83 (6.39%) 11.83 (6.39%) 9.60 (5.19%) 12.82 (6.93%)
25.98 25.98 25.73 32.49
30.12 30.19 30.20 30.23
58.19 58.21 58.21 67.51
7.28 (3.94%) 7.28 (3.94%) 0.45 (0.24%) 3.94 (2.13%)
14.85 14.85 13.05 19.92
(8.09%) (8.09%) (5.97%) (6.14%)
3.69 (4.91%) 3.69 (4.91%) 0.51 (0.68%) 0.38 (0.51%)
(8.28%) (8.28%) (8.20%) (10.35%) (4.73%) (4.73%) (4.16%) (6.35%)
(5.76%) (5.77%) (5.77%) (5.78%)
16.42 16.48 16.47 21.85
(3.14%) (3.15%) (3.15%) (4.18%)
(8.70%) (8.70%) (8.70%) (10.09%)
43.66 43.67 43.67 21.12
(6.53%) (6.53%) (6.53%) (3.16%)
1536
C. Housmans et al. / Renewable Energy 113 (2017) 1529e1538
Table 7 Absolute (in W m2) and relative (in %) RMSE and MBE indexes for a three pyranometers system as a function of the solar elevation, g. Results are provided for the four considered transposition models (see Appendix B). The relative indexes are displayed in brackets. Validation statistics are computed from 4066 data points (e.g., 69.91% of the total available data points) using the EM method. Transposition model RMSE LIU HAY SKA PER MBE LIU HAY SKA PER
Solar elevation angle 5+ g < 30+
30+ g < 60+
60+ g 90+
19.70 19.77 19.55 20.47
20.41 20.40 19.59 22.75
(5.45%) (5.45%) (5.23%) (6.07%)
15.37 15.37 10.76 12.37
9.26 9.26 4.50 9.06
(2.47%) (2.47%) (1.20%) (2.42%)
8.62 (2.79%) 8.62 (2.79%) 3.42 (1.11%) 1.41 (0.46%)
(14.00%) (14.04%) (13.89%) (14.54%)
10.08 (7.16%) 10.11 (7.18%) 7.32 (5.20%) 9.07 (6.45%)
(4.98%) (4.98%) (3.49%) (4.01%)
SKA model for which a slight overestimation of the solar radiation of about 3.4 W m2 (1.1%) is reported at high solar elevation conditions while the model underestimates the global radiation by about 7.3 W m2 (5.2%) in low solar elevation conditions. Finally, comparing the performance between the isotropic and anisotropic approaches to the inverse transposition problem as reported in Tables 5e7 indicates that the improvement from using the LIU model to using anisotropic models is marginal in terms of RMSE and MBE and even inexistent in term of failure rate. 5. Conclusions Different approaches to convert global solar irradiance measured on a tilted plane to horizontal have been evaluated using solar radiation data measured on the roof of the radiometric tower of the Royal Meteorological Institute of Belgium in Uccle, Brussels. Up to three sensors have been considered in the horizontal solar irradiance reconstruction. When only a single tilted sensor is used, the conversion can be carried out with a decomposition model coupled with a transposition model to solve the inverse transposition problem. In this case, there is an additional error (additional to the inverse transposition problem) in the predicted horizontal irradiance. In addition, with only one tilted irradiance involved in the inverse modeling approach, the conversion's performance is very sensitive to the tilted pyranometer angular configuration (i.e., the tilt angle and the surface's orientation). None of the considered transposition models was found to outperform over the 3 pyranometers mounting plane configurations. When two (or more) tilted irradiance sensors are involved, only a transposition model is required and the inverse transposition problem can be solved more accurately. In practice, there are certain sun positions for which this procedure still fails to produce a valid estimation of the global horizontal irradiance. Consequently more instruments should be used so as to overdetermine the system. Here it has been shown that three tilted pyranometers set at different orientations is sufficient to guarantee solution for three of the four considered transposition models. Indeed, because of the non bijectivity of the anisotropic transposition model of Perez et al. [23] the conversion of a percentage of data points are unsuccessful. Sensitivity experiments conducted in the case of a two tilted pyranometers system indicate that solving the inverse transposition quadratic equation by minimizing the errors instead of using the least square procedure tends to significantly reduces the percentage of unsuccessful conversion irrespective of the considered transposition model. This is also true for the conversion performed with the isotropic model even if in this case the inverse
transposition problem reduces to a first degree equation with one unknown as the Powell's minimization algorithm overdetermines the system and provides the best fitting solution. Moreover the accuracy of the tilt to horizontal conversion is improved with the errors minimization method irrespective of the considered transposition model. Finally, comparing the performance between the isotropic and anisotropic approaches to the inverse transposition problem in angular configurations similar to those encountered in Belgian PV systems installations (e.g., with tilt angle as great as 50.79 ) indicates that the improvement from using the LIU isotropic model to using anisotropic models is not significant (e.g., a rRMSE of 7.9% and a rMBE of 9.7% are reported for the LIU model when considering a three tilted pyranometers system compared to a rRMSE of 7.7% and a rMBE of 5.9% for the best performing SKA anisotropic model) or even inexistent regarding the percentage of unsuccessful conversion (e.g., a failure rate of 0% is reported for the LIU model in the three tilted pyranometers system compared to a failure rate of 30.1% for the PER model). Acknowledgments This study was supported by the Belgian Science Policy Office (BELSPO) through the BRAIN-be research project:”Solar irradiation from the energy production of residential PV systems eSPIDER00 , contract No. BR/314/PI/SPIDER. Appendix A. Skartveit and Olseth (1987) decomposition model (OLS model) Skartveit and Olseth (1987) estimated the direct normal irradiance, Bn, from the global horizontal irradiance, Gh , and the solar elevation angle, g, for Bergen (Norway, 60.4oN) with the following equation based on hourly records of global and diffuse horizontal irradiances with mean solar elevation larger than 10 during 1965e1979:
Bn ¼
Gh ð1 JÞ sin g
(20)
where J is a function of the clearness index, Kt . The model was validated with data collected in 12 stations worldwide. The function J reads as:
8 1 for Kt < c1 > > i h > > 1=2 < 2 1 ð1 d1 Þ d2 c3 þ ð1 d2 Þ c3 for c1 Kt 1:09c2 J¼ > > > > : 1 1:09c 1 Y for Kt > 1:09c2 2 Kt (21) where: c1 ¼ 0:2 g c2 ¼ 0:87 0:56 0:06 c3 ¼ 0:5 1 þ sin p dc4 0:5 3 c 4 ¼ Kt c 1 0:06 g d1 ¼ 0:15 þ 0:43 e d2 ¼ 0:27 d3 ¼ c2 c1 1=2 Y ¼ 1 ð1 d1 Þ ðd2 cc3 þ ð1 d2 Þ cc23 Þ 4 0:5 cc3 ¼ 0:5 1 þ sin p cc d 3
cc4 ¼ 1:09 c2 c1
C. Housmans et al. / Renewable Energy 113 (2017) 1529e1538
Appendix B. Transposition models
1537
inclined surface. More specifically, the angular location of the circumsolar region is determined by the ratio a/b. They are calculated from the equations of solar geometry:
Each model develops the diffuse transposition factor (i.e., the ratio of diffuse radiation on a tilted surface to that of a horizontal, Rd ) according to specific assumptions.
a ¼ maxð0; cos qi Þ
B.1. LIU model (Liu and Jordan, 1962):
b ¼ maxðcos 85o ; cos qz Þ
The isotropic model (that assumes a uniform distribution of the diffuse radiation over the sky dome) was published in 1962:
The brightness coefficients F1 and F2 are derived from the so-called Perez coefficients:
RLIU d ¼
1 þ cosb 2
(22)
where b denotes the surface tilt angle with respect to the horizontal plane. B.2. HAY model (Hay, 1979): In the HAY model, diffuse radiation from the sky is composed of an isotropic component and a circumsolar one. Horizon brightening is not taken into account. An anisotropy index, FHay , is used to quantify a portion of the diffuse radiation treated as circumsolar with the remaining portion of diffuse radiation assumed to be isotropic, i.e.,
RHAY ¼ FHAY rb þ ð1 FHAY Þ d
1 þ cos b 2
(23)
where FHAY ¼ Bh =ðIo cosqz Þ is the Hay's sky-clarity factor and rb ¼ cosqi =cosqz the beam radiation conversion factor. qz is the solar zenith angle, qi is the incidence angle of the beam radiation on the tilted surface, Bh is the direct horizontal solar irradiance and Io is the extraterrestrial normal incident irradiance. The model reduces to the LIU model (Eq. (22)) for FHAY ¼ 0. B.3. SKA model (Skartveit and Olseth, 1986): Solar radiation measurements indicate that a significant part of sky diffuse radiation under overcast sky conditions comes from the sky region around the zenith. This effect vanishes when cloud cover disappears. Skartveit-Olseth (1986) modified the HAY model (Eq. (23)) in order to account for this effect,
RSKA ¼ FHAY rb þ Zcosb þ ð1 FHAY ZÞ d
1 þ cosb 2
(24)
where Z ¼ maxð0; 0:3 2FHAY Þ is the Skartveit-Olseth's correction factor. If FHAY 0:15, then Z ¼ 0 and the model reduces to the HAY model. B.4. PER model (Perez et al., 1987): Compared to the previously described transposition models, the model proposed by Perez et al. (1987) represents a more detailed analysis of the isotropic diffuse, circumsolar and horizon brightening radiation by using empirically derived coefficients. According to this model,
a 1 þ cosb þ F2 sinb RPER ¼ F1 þ ð1 F1 Þ d b 2
(25)
where F1 and F2 are sky brightness coefficients for the circumsolar region and the region above the horizon line, respectively. Note that if F1 ¼ F2 ¼ 0, it reduces to the LIU model (Eq. (22)). The coefficients a and b take into account the angle of incidence of the sun onto the
F1 ¼ F11 ðεÞ þ F12 ðεÞD þ F13 ðεÞqz F2 ¼ F21 ðεÞ þ F22 ðεÞD þ F23 ðεÞqz where the Perez coefficients Fij are function of the sky clearness parameter ε and the sky brightness parameter D. These factors are defined by:
ε¼
ðDh þ Bh Þ=Dh þ 1:041q3z
D¼m
1 þ 1:041q3z Dh Dh ¼ Io Io cosqz
where Dh is the diffuse horizontal solar irradiance, m is the optical air mass and qz is in radians. Many sets of Perez coefficient values have been determined by different studies. In this paper, we applied the set of coefficients from Perez et al. (1990).
References [1] M. Wali, F. Evrendilek, T. West, S. Watts, D. Pant, H. Gibbs, et al., Assessing terrestrial ecosystem sustainability: usefulness of regional carbon and nitrogen models, Nat. Resour. 35 (4) (1999) 20e33. [2] M. Trnka, Z. Zalud, J. Eitzinger, M. Dubrovsky, Global solar radiation in central european lowlands estimated by various empirical formulae, Agric. For. Meteorol. 131 (2005) 54e76. e, C. Bertrand, Quality control of solar radiation data within the rmib [3] M. Journe solar measurements network, Sol. Energy 85 (2011) 72e76. [4] R. Perez, R. Seals, R. Stewart, A. Zelenka, V. Estrada-Cajigal, Using satellitederived insolation data for the site/time specific simulation of solar energy systems, Sol. Energy 53 (1994) 491e495. [5] R. Perez, R. Seals, A. Zelenka, Comparing satellite remote sensing and ground network measurements for the production of site/time specific irradiance data, Sol. Energy 60 (1997) 89e96. [6] A. Hammer, D. Heinemann, C. Hoyer, R. Kuhlemann, E. Lorentz, R. Mueller, et al., Solar energy assessment using remote sensing technologies, Remote Sens. Environ. 86 (2003) 423e432. , R. George, S. Wilcox, T. Stoffel, D. Myers, D. Heimiller, Solar Resource [7] D. Renne Assessment. Technical Report Nrel/tp-581-42301, National Renewable Energy Laboratory, 2008, 1617 Cole Boulevard, Golden, Colorado:80401e3393. [8] D. Yang, Solar radiation on inclined surfaces: corrections and benchmarks, Sol. Energy 136 (2016) 288e302. [9] D. Yang, Z. Dong, A. Nobre, Y. Khoo, P. Jirutitijaroen, W. Walsh, Evaluation of transposition and decomposition models for converting global solar irradiance form tilted surface to horizontal in tropical regions, Sol. Energy 97 (2013) 369e387. [10] B. Marion, A model for deriving the direct normal and diffuse horizontal irradiance from the global tilted irradiance, Sol. Energy 122 (2015) 1037e1046. [11] R. Perez, P. Ineichen, E. Maxwell, R. Seals, A. Zelenka, Dynamic global-to-direct irradiance conversion models, ASHRAE Trans.-Res. Ser. 3578 (1992) 354e369. [12] F. Olmo, J. Vida, I. Foyo, Y. Castro-Diez, L. Alados-Arboledas, Prediction of global irradiance on inclined surfaces from horizontal global irradiance, Energy 24 (1999) 689e704. [13] S. Killinger, et al., Projection of power generation between differentlyoriented pv systems, Sol. Energy 136 (2016) 153e165. [14] D. Yang, Z. Ye, A. Nobre, H. Du, W. Walsh, L. Lim, et al., Bidirectional irradiance transposition based on the perez model, Sol. Energy 110 (2014) 768e780. e, C. Bertrand, Evaluation of different models to estimate [15] C. Demain, M. Journe the global solar radiation on inclined surfaces, Renew. Energy 50 (2013) 710e721. e, C. Bertrand, Corrigendum to ”evaluation of different [16] C. Demain, M. Journe
1538
[17]
[18] [19] [20] [21] [22] [23]
[24] [25] [26]
[27] [28]
C. Housmans et al. / Renewable Energy 113 (2017) 1529e1538 models to estimate the global solar radiation on inclined surfaces”[renewable energy 50 (2013) 710-721], Renew. Energy 101 (2017) 1401e1403. e, Evaluation of decomposition models C. Bertrand, G. Vanderveken, M. Journe of various complexity to estimate the direct solar irradiance over Belgium, Renew. Energy 74 (2015) 618e626. C.A. Gueymard, Direct and indirect uncertainties in the prediction of tilted irradiance for solar engineering applications, Sol. Energy 83 (2009) 432e444. A. Skartveit, J. Olseth, A model for the diffuse fraction of hourly global radiation, Sol. Energy 38 (1987) 271e274. B. Liu, R. Jordan, Daily insolation on surfaces tilted towards the equator, Trans. ASHRAE 67 (1962) 526e541. J. Hay, Study of Shortwave Radiation on Non-horizontal Surfaces. Rep No 7912, Atmospheric Environment Service, Downsview, Ontario, 1979. A. Skartveit, J. Olseth, Modelling slope irradiance at high latitudes, Sol. Energy 36 (1986) 333e344. R. Perez, R. Seals, P. Ineichen, R. Stewart, D. Menicucci, A new simplified version of the perez diffuse irradiance model for tilted surfaces, Sol. Energy 39 (1987) 221e231. D. Faiman, A. Zemel, A. Zangvil, A method for monitoring insolation in remote regions, Sol. Energy 38 (1987) 327e333. C. Grosan, A. Abraham, A new approach for solving nonlinear equations systems, IEEE Trans. Syst. Man. Cybern. Part A Syst. Humans 38 (2008) 698e714. M. Powell, An efficient method for finding the minimum of a function of several variables without calculating derivatives, Comput. J. 7 (2) (1964) 155e162. F. Acton, Numerical Methods that Work, 1990, corrected edition, Mathematical Association of America, Washington, 1970, pp. 464e467. W. Press, S. Teukolsky, W. Vetterling, B. Flannery, Numerical Recipes in C. The Art of Scientific Computing (Reprinted with corrections, 1997), second ed., Cambridge University Press, Cambrige, UK, 1992, p. 994.
€ppen-geiger [29] M. Peel, B. Finlayson, T. McMahon, Updated world map of the ko climate classification, Hydrol. Earth Syst. Sci. 11 (2007) 1633e1644.
Nomenclature Io : extraterrestrial normal incident irradiance ðW m2 Þ Gt : tilted global solar irradiance ðW m2 Þ Gh : global horizontal irradiance ðW m2 Þ Bt : in-plane direct irradiance ðW m2 Þ Bn : direct normal irradiance ðW m2 Þ Bh : direct irradiance on a horizontal surface ðW m2 Þ Dt : in-plane diffuse irradiance ðW m2 Þ Dh : diffuse horizontal irradiance ðW m2 Þ Dg : irradiance due to the ground reflection ðW m2 Þ rb : factor that accounts for the direction of beam radiation Rd : diffuse transposition factor Rr : transposition factor for ground reflection Kt : effective global horizontal transmittance (or clearness index) Kn : direct normal transmittance Kd : diffuse horizontal transmittance g: solar elevation angle (rad) a: azimuth of the tilted plane (rad) b: plane tilt angle (rad) qi : incidence angle (rad) qz : zenith of the Sun (rad) Р: ground albedo ε: sky clearness parameter Fij : Perez model coefficient Z: Skartveit-Olseth's correction factor
Contents lists available at ScienceDirect
Renewable Energy journal homepage: www.elsevier.com/locate/renene
Tilt to horizontal global solar irradiance conversion: An evaluation at high tilt angles and different orientations dric Bertrand Caroline Housmans*, Alessandro Ipe, Ce Royal Meteorological Institute of Belgium, Avenue Circulaire 3, B-1180, Brussels, Belgium
a r t i c l e i n f o
a b s t r a c t
Article history: Received 5 January 2017 Received in revised form 29 June 2017 Accepted 3 July 2017 Available online 5 July 2017
Many transposition models have been proposed in the literature to convert solar irradiance on the horizontal plane to that on a tilted plane. The inverse process, i.e., the conversion from tilted to horizontal is investigated here based on six months of in-plane global solar irradiance measurements recorded on the roof of the Royal Meteorological Institute of Belgium's radiation tower in Uccle (Latitude 50.79 N, Longitude 4.35 E). Up to three pyranometers mounted on inclined planes of different tilts and orientations were involved in the inverse transposition process. Our results indicate that (1) the tilt to horizontal irradiance conversion is improved by using a multi-pyranometers approach and (2) the improvement from using an isotropic model to anisotropic models in the inverse transposition problem is not significant. © 2017 Elsevier Ltd. All rights reserved.
Keywords: Solar radiation Transposition model Decomposition model
1. Introduction Knowledge of local solar irradiation is essential for many applications, including architectural design, solar energy systems, crop growth models, and evapotranspiration estimates (e.g., Wali et al. [1]; Trnka et al. [2]). The usual solar radiation parameters measured on ground are the global horizontal irradiance, Gh , the direct normal irradiance, Bn , the diffuse horizontal irradiance, Dh , and the sunshine duration. Traditionally, they are recorded by means of networks of meteorological stations. However, costs for installation and maintenance of such networks are very high and national networks comprise only few stations. For instance, the solar measurements network operated by the Royal Meteorological Institute of Belgium (RMI) has only 15 sites where various combinations of the solar irradiation components are recorded (e.g., e and Bertrand [3]). Consequently the availability of Journe observed solar radiation measurements has proven to be spatially and temporally inadequate for many applications. Over the last decades, satellite-based retrieval of solar radiation at ground level has proven to be valuable for delivering a global coverage of the global solar irradiance distribution at the Earth's et al. [7]). surface (e.g., Perez et al. [4]- [5]; Hammer et al. [6]; Renne The recent deployments of solar photovoltaic (PV) systems offer a
* Corresponding author. E-mail address: [email protected] (C. Housmans). http://dx.doi.org/10.1016/j.renene.2017.07.018 0960-1481/© 2017 Elsevier Ltd. All rights reserved.
potential opportunity of providing additional solar resource information as PV systems with global tilted irradiance, Gt , measurements become available in some locations. However, in most PV systems, modules are installed on a fixed plane to reduce installation and operation costs compared to tracking systems. To maximize array output, modules are installed at a tilt close to local latitude, or at some minimum tilt to ensure self-cleaning by rain. Gh estimates at the PV system location require therefore to convert the tilted solar irradiance recorded to that on a horizontal plane. If many transposition models have been proposed in the literature (see Yang [8] for a review) to convert solar irradiance on the horizontal plane to that on a tilted plane, the inverse process (i.e., converting from tilted to horizontal) is only poorly discussed in literature. The difficulty relies on the fact the procedure is analytically not invertible. As an example, single-sensor approach involving a numerical search method was proposed by Yang et al. [9] and evaluated on various combinations of decomposition (i.e., models that separate direct and diffuse solar components from the global one) and transposition models. They found their method sufficiently accurate for small zenith angles, but reported that the conversion error exponentially increases as zenith angle increases. More recently, Marion [10] presented an iterative method using a modified version of the DIRINT decomposition model (Perez et al. [11]) in combination with the transposition model of Perez et al. [5] to obtain horizontal irradiance from a single tilted sensor measurement. Performance of this method was found essentially the
1530
C. Housmans et al. / Renewable Energy 113 (2017) 1529e1538
same as for the DIRINT model used in the normal calculation mode for Gt measurements recorded at small tilt angle from the horizontal. At larger angles, deviations between modeled and measured direct and diffuse irradiance values increase as the tilt angle increases from the horizontal. Lastly, an alternative method based on the Olmo model (Olmo et al. [12]) that presents the particularly of being analytically invertible was proposed by Killinger et al. [13]. If the overall performance of the inverted Olmo model was found comparable with the other approaches, the results were slightly worse than those obtained by inverting the decomposition and transposition models in combination with an iterative solving process. The present study aims at (1) evaluate the tilt to horizontal irradiance conversion when using for model input Gt measurements recorded in angular configurations similar to those encountered in Belgian PV systems installations (e.g., with tilt angle as great as 50.79 ) and (2) determine if as suggested by Yang et al. [14] the use of a multi-pyranometers system could improve the conversion accuracy. The selection of the decomposition and transposition models used in our calculations is based on previous evaluation of popular decomposition and transposition models performance in Uccle by Demain et al. [15]- [16] and Bertrand et al. [17]. The paper is organized as follows: methods to perform the conversion from tilt to horizontal are presented in section 2. In situ measurements are briefly described in section 3. Performances of the different approaches are evaluated in section 4. Final remarks and conclusions are provided in section 5.
Transposition models have the general form:
(1)
where the tilted global solar irradiance, Gt , is expressed as the sum of the in-plane direct irradiance, Bt , in-plane diffuse irradiance, Dt , and the irradiance due to the ground reflection, Dg . The direct component, Bt , is obtained from:
Bt ¼ Bn cosqi ¼ Bh
cosqi ¼ Bh rb cosqz
Dt ¼ Dh Rd
(3)
Dg ¼ rGh Rr
(4)
where, Dh , is the diffuse horizontal irradiance, Gh , the global horizontal irradiance (i.e., Gh ¼ Dh þ Bh ), Rd , the diffuse transposition factor and r, the ground albedo. The transposition factor for ground reflection, Rr , can be modeled under the isotropic assumption (e.g., Gueymard [18]) as follows:
Riso r ¼
1 cosb 2
(2)
(5)
where, b, is the tilt angle of the inclined surface. See Fig. 1 for angles definition. 2.1. Single-sensor approach Considering the effective global horizontal transmittance, Kt , the direct normal transmittance, Kn , and the diffuse horizontal transmittance, Kd (i.e., Kd þ Kn ¼ Kt ):
Gh ¼ Kt Io cos qz Bn ¼ Kn Io Dh ¼ Kd Io cos qz
2. Conversion from tilt to horizontal
Gt ¼ Bt þ Dt þ Dg
with, Bn , the direct normal irradiance and, Bh , the direct irradiance on a horizontal surface, respectively. qi is the incidence angle and, qz , the solar zenith angle, respectively. Parameter rb ¼ cos qi =cos qz is a factor that accounts for the direction of the beam radiation. The diffuse component, Dt , and the irradiance due to the ground reflection, Dg , can be modeled as follows:
(6)
(where Io is the extraterrestrial normal incident irradiance), Eq. (1) can be rewritten as:
K K Gt ¼ Kt Io cos qi 1 d þ cos qz d Rd þ rRr Kt Kt
(7)
Eq. (7) indicates that the conversion of Gt to Gh requires the use of a decomposition model to estimate Kd from Kt . The decomposition model of Skartveit and Olseth [19] (hereafter referred to as OLS) was used in the present study. Formulation of this model is
Fig. 1. Angles to define the position of the sun and the orientation of a tilted plane: qz is the solar zenith angle, g, the solar elevation, qi , the incidence angle, b, the tilt angle of the inclined surface and the azimuth of the tilted plane is a.
C. Housmans et al. / Renewable Energy 113 (2017) 1529e1538
provided in Appendix A. It allows to estimate the direct normal irradiance, Bn , from the global solar irradiance, Gh , through an empirical equation depending on Kt . Performance of Eq. (7) has been evaluated for various geometric configurations using the isotropic transposition model proposed by Liu and Jordan [20] (hereafter referred to as LIU model) and the anisotropic models of Hay [21] (hereafter referred to as HAY model), Skartveit and Olseth [22] (hereafter referred to as SKA model) and Perez et al. [23] (hereafter referred to as PER model), respectively. Formulation of the selected transposition models is provided in Appendix B. Note that because the PER model coefficients are binned according to the sky clearness parameter, ε, (see Appendix (B.4)), the same strategy as in Yang et al. [14] was adopted here when running the PER model. Basi0 0 cally, Kt and Kd estimates are calculated for each ε bins using the associated model coefficients Fij . In each ε bins Gh 0 and Dh 0 esti0 0 mates are then retrieved from the corresponding Kt and Kd esti0 mates and used to determine the corresponding ε value. If the 0 0 assumed ε agrees with the calculated ε , the corresponding Gh and 0 Dh are selected as the true estimated Gh and Dh . Similarly, because the transposition factor defined by the SKA model depends on the Skartveit-Olseth's correction factor Z (see Eq. (24)), the conversion process is performed for the two Rd model formulations before selecting the appropriate solution. The global horizontal irradiance found has to satisfy the Z value definition assumed during the conversion. We solved Eq. (7) by a simple iteration procedure. Namely, Kt is adjusted using the bisection procedure until the Kt provides calculated Gt value that match the measured Gt (using as starting point Kt ¼ 0:5). Note that for the sake of simplicity, the foreground's albedo, r, was taken to 0.23 (Demain et al. [15]) in our calculations.
l¼
1531
Gt2 B1 Gt1 B2 Gt1 A2 Gt2 A1
(12)
Finally, substituting Eq. (12) into Eq. (10) leads to a first degree equation with one unknown: 0
A Dh þ B0 ¼ 0
(13)
where. 0
A ¼ B1 ðGt2 A1 Gt1 A2 Þ 0 B ¼ Gt1 ðGt1 A2 þ A1 B2 Þ Gt2 A1 B1 Gt1 Gt2 A2 Considering the anisotropic approach to the inverse transposition problem, Eq. (8) can be rewritten as:
Gt ¼ ABh Dh þ BBh þ CDh
(14)
where the formulation of the coefficients A, B, and C, depends on the considered transposition model. Following the same reasoning as for the isotropic approach, leads in the anisotropic case to solve a quadratic equation: 00
00
00
A D2h þ B Dh þ C ¼ 0
(15) 00
00
00
where the coefficients A , B and C depends on the considered transposition model. Given n tilted pyranometers (with different inclinations and/or orientations), Eq. (14) can be generalized as the matrix form (e.g., Yang et al. [14]):
(16)
2.2. Multi-sensors approach The idea that simultaneous readings of a multi-pyranometers system can be used to disangle the various components of solar radiation on inclined surfaces has been originally proposed to solve in remote locations the periodic adjustment required by normal incidence and shadow-band pyranometers to ensure that their readings remain accurate when long-term data acquisition is in progress (Faiman et al. [24]). By combining Eqs. (2)e(4) we rewrite Eq. (1) as:
Gt ¼ Bh rb þ Dh Rd þ rGh Rr
(8)
Assuming first Rd as isotropic, Eq. (8) becomes:
Gt ¼ Dh ðlA þ BÞ
(9)
where.
where L ¼ fAi g is a 2 n 2 third-order tensor, B ¼ fBi g is a n 2 matrix, ℭ is a n column vector, and x a column vector with 2 variables:
0 Ai 2<2n2 0 Ai 0 1 C1 B1 B C2 B2 C C2
(17)
Note that Eqs. (16) and (17) slightly differ from the formulations given in Yang et al. [14] due to the change in the column vector x, i.e., (Dh Bh ) instead of (Dh l). The least square (hereafter referred to as LS) solution to Eq. (16) is given by:
cos qi þ rRr cos qz 1þcos b þ rRr 2
A¼ B¼ l ¼ Bh =Dh
Supposing a two-pyranometers system with Gt1 and Gt2 the corresponding tilted global irradiance measurements, it follows from Eq. (9) that:
Gt1 ¼ Dh ðlA1 þ B1 Þ
(10)
Gt2 ¼ Dh ðlA2 þ B2 Þ
(11)
Dividing Eq. (10) by Eq. (11) to eliminate Dh yields, after rearranging:
(18) with jj:jj referring to the Euclidean norm. However, the LS is hard to solve and a standard technique to resolve Eq. (18) is to use a Newton type iteration method (e.g., Grosan and Abraham [25]). As an alternative, Eq. (16) can also be solved by minimizing the errors (this approach is hereafter denoted to as EM - Errors Minimization). In this case, the solution is to minimize
1532
C. Housmans et al. / Renewable Energy 113 (2017) 1529e1538
Table 1 Data availability percentage according to the sky condition for the considered period from July 15, 2015 to January 19, 2016. Sky condition
Kt
Data availability (%)
All sky Overcast Cloudy Partly cloudy Partly clear Clear
0.0e1.0 0.0e0.2 0.2e0.4 0.4e0.6 0.6e0.8 0.8e1.0
100% 32.07% 26.92% 18.01% 21.46% 1.54%
( min EðxÞ ¼
n X
) ε2i ðxÞ : x2<2
(19)
i¼1
where, εi ðxÞ ¼ ðAi Dh Bh þ Bi Bh þ Ci Dh Þ Gti , with i ¼ 1; …; n denoting the tilted pyranometer (n 2). Both ways to solve Eq. (16) have been considered and compared here using the Powell's quadratically convergent method (Powell [26]) to perform the minimization. It is a generic minimization method that allows to minimize a quadratic function of several variables without calculating derivatives. The key advantage of not requiring explicit solution of derivatives is the very fast execution time of the Powell method. In order to avoid the problem of linear dependence in Powell's algorithm, we adopted the modified Powell's method given in Acton [27] and implemented in Press et al. [28]. 3. In situ measurements The present analysis is based upon six months of global solar irradiance measurements performed on the roof of the radiometric tower of the Royal Meteorological Institute of Belgium located on the Brussels Uccle plateau (Latitude 50.79 N, Longitude 4.35 E at an elevation of 101 m above sea level). Belgian climate, like most of the Northern Europe (Peel et al. [29]) is maritime temperate with significant precipitation in all seasons. Cloud occurrence is high (see Table 1) and clouds are the dominant
source of small-scale variability in surface solar radiation. Fig. 2 compares the global solar irradiation (in kWh=m2 ) received in Uccle during the six months of the measurements campaign (black curve) with the corresponding climatological normals (i.e. the 1981e2010 average, in green). The highest and lowest values recorded since 1981 are also provided for illustration (in red and blue, respectively). Four data sets of global solar irradiance have been collected from July 15, 2015 to January 19, 2016. The first one was recorded by a Kipp & Zonen CM11 Secondary Standard pyranometer mounted on the horizontal plane (here after referred to as Plane 0). The second was also recorded by a Kipp & Zonen CM11 Secondary Standard pyranometer but mounted on a tilted plane of 50.79 (i.e., corresponding to the RMI radiometric tower's latitude) facing south (here after referred to as Plane 1). For the last two data sets, two additional Kipp & Zonen CMP-22 Secondary Standard pyranometers were installed on the tower with the same angular configurations than two neighboring residential PV installations, i.e., a tilted plane of 45 SW facing and a tilted plane of 50 facing E (here after referred to as Plane 2 and Plane 3, respectively). Regarding the accuracy of the measurements devices (see Fig. 3 for instruments illustration), the directional error (up to 80 with 1000 Wm2 beam) is less than 5 Wm2 for the CMP22 pyranometer and less than 10 Wm2 for the CM11 pyranometer, respectively. The pyranometer spectral selectivity (300e1500 nm) is smaller than 3% for the CMP22 instrument and smaller than 2% for the CM11 instrument, respectively. Both instruments have a thermopile detector and present a reduced thermal offset (i.e., about 2 Wm2 at 5 K/h temperature change) for the CM11 pyranometer and less than 1 Wm2 for the CMP22 pyranometer, respectively. Irradiance measurements were made with a 5-sec time step and then integrated to bring them to a 10-min time step. Note that because of the difference in the recording planes orientations and inclinations, the number of available measurements differs from one instrument to another (see Table 2). All 10-min data have undergone a series of automated quality control procedures similar to e and Bertrand [3] prior to be visually those described in Journe inspected by a human operator.
Fig. 2. Comparison between the global solar irradiation (in kWh=m2 ) received in Uccle during the six months (i.e., from July 15, 2015 to January 19, 2016) of the measurements campaign (black curve) with the corresponding climatological normals (i.e., the 1981e2010 average, in green). The highest and lowest values recorded since 1981 are also provided for illustration (in red and blue, respectively). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
C. Housmans et al. / Renewable Energy 113 (2017) 1529e1538
1533
Fig. 3. Panel A presents an illustration of the Kipp & Zonen CM11 Secondary Standard pyranometer instrument used to perform Plane 0 and Plane 1 global solar irradiance measurements. Similarly panel B shows the Kipp & Zonen CMP-22 Secondary Standard pyranometer instrument considered to perform Planes 2 and 3 global solar irradiance measurements. Note that in panel A, the instrument is illustrated in the mounting Plane 0 angular configuration on the roof of the RMI's radiometric tower while in panel B, the instrument is given in the mounting Plane 3 angular configuration.
Table 2 Angular configuration of the four recording planes considered in the study. Note that the surface azimuthal angle, a, is conventionally measured clockwise from the South. Also provided is the number of available data points for each of the different planes. Note that the total number of available simultaneous data points in all of the recording planes is 5816 (i.e., 100%). Pyranometer configuration
Azimuth Angle a [ ]
Tilt Angle b [ ]
Data Points Number
Plane Plane Plane Plane
/ 0.00 (S facing) 45.00 (SW facing) 90.00 (E facing)
0.00 50.79 45.00 50.00
10700 10578 9631 7269
0 1 2 3
4. Results The relative ability of the different approaches to predict horizontal irradiance from the tilted one was estimated by means of two statistical error indexes: Mean Bias Error (MBE) and Root Mean Square Error (RMSE).
MBE ¼
n 1X ðe Þ n i¼1 i
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u n u1 X e2 RMSE ¼ t n i¼1 i where ei ¼ ðGi;e Gi;o Þ is the residual value; Gi;e are the estimated values and Gi;o represent the observed measurements. A positive MBE (resp. a negative MBE) means that the model tends to overestimate (resp. underestimate) the observed measurements. To obtain dimensionless statistical indicators we expressed MBE and RMSE as fractions of mean solar global irradiance during the respective time interval,
MBE½% ¼
MBE M
RMSE½% ¼
RMSE M
P where M ¼ 1n ni¼1 ðGi;o Þ is the measurements mean. Note that for a proper estimation of the error statistics, only radiation data recorded with a solar zenithal angle, qz , smaller than 85 and an incidence angle on the plane used, qi , smaller than 90 were considered and tilted global solar irradiance records were further restricted to
non-zero values. Moreover when the retrieved global horizontal solar irradiance was negative or larger than the corresponding extraterrestrial irradiance the inverse transposition problem was considered unsuccessful. With these additional conditions the total number of data points available simultaneously in all of the recording planes reduces to 5816 (i.e. 100%). 4.1. Single-sensor approach Table 3 summarizes for each of the three tilted planes the performance of the inverted transposition models in the Gt to Gh conversion. To evaluate the transposition models on a same basis, only the 5117 data points conversions that have succeeded for all models have been considered in the MBE and RMSE calculations (e.g., 5117/5816 ¼ 87.98% of the total available data points, 5816 being the total number of measurements available simultaneously for the four planes). It is obvious from Table 3 that none of the considered transposition models performs best in terms of failure rate, RMSE and MBE over the three tilted planes. Moreover, the tilt angle and the surface's orientation have a major impact on the Gh estimation's reliability irrespective of the considered transposition model. The worst results in term of RMSE are for Plane 2 measurements conversion (i.e., rRMSE ranging from 28.7% to 30.8% vs 6.8%e12.2% reported for Plane 1 conversions). Furthermore, Plane 2 conversions underestimate Gh (i.e., rMBE of 14.1% to 11.7%) while a slight overestimation of 0.3%e5.3% is reported for Plane 1 and an overestimation of 2.3%e7.6% for Plane 3, respectively. Regarding the failure rate, the situation varies according to the considered transposition model. As an illustration, while almost all tilted irradiance measurements from Plane 1 are successfully converted when performing the inverse transposition with the HAY and SKA models (i.e., failure rate of 0.01%), the largest percentage of unsuccessfull conversions reported for the LIU and PER models is for the Plane 1 measurements conversion (i.e., a failure rate of 1.6% and 4%, respectively). For comparison, Marion (2015) reported a failure rate of 10% for its iterative method combining the DIRINT decomposition model with the PER transposition model. Yang et al. [9] reported rRMSE values of roughly 13% and 6% for the tilt to horizontal conversion of irradiance data collected in the tropical region of Singapore by sensors mounted on a 18.3 tilted NE facing plane and 6.1 SW facing plane, respectively. The overall performance of our Gt to Gh conversion appears relatively lower especially for Planes 2 and 3 measurements conversion. However, it is worth pointing out that (1) the tilt angles considered in Yang et al. [9] were significantly less than here and (2) Yang et al. [14] reported that the RMSE resulting from their conversion using a single tilted
1534
C. Housmans et al. / Renewable Energy 113 (2017) 1529e1538
resolution of Eq. (13) (i.e., a first degree equation with one unknown). Validation statistics were computed from the 1499 data points that succeeded the inverse transposition with all the considered transposition models and the two solving methods (e.g., 25.77% of the total available data points). It is clear from Table 4 that both solving methods produce very distinct results in terms of RMSE and MBE as well as in term of failure rate. Clearly, the EM method allows to reduce the occurrence of invalid solutions in the inverse transposition problem. As an example, the failure rate associated to Plane 1 - Plane 2 conversions carried out by running the LIU (HAY) model reduces from 8.19% (33.23%) with the quadratic formula to 2.03% (2.01%) with the EM method. However, it is worth pointing out that if the failure rate associated to the PER model is found to generally decrease with the EM method, the number of invalid solutions reported for this model is very large. Because the PER model depends on ε binning, it can appear that the retrieved ε value using the procedure described in section 2.1 does not agree with the expected bin, or that more than one retrieved ε value agree with the expected bin, increasing therefore the number of unsuccessful conversions associated to this model. Using the EM method instead of the LS procedure does not solve the ε binning formulation deficiencies; in contrary the occurrence of the problem is even enlarged with the EM method for the Plane 1 þ Plane 3 conversions (failure rate of about 37% with the EM method compared to 21% with the LS solution). By contrast, in depth-analysis of the conversion process indicated that for a number of data points, conversion performed with the HAY and SKA models using the EM method allows to produce a valid result when the discriminant or the roots of the quadratic equation are found negative or when two positive real solutions are found (as in this latter case, no objective criterion permitted us to select the appropriate solution and the conversion was considered unsuccesfull). Regarding the LIU isotropic transposition model, even if the inverse transposition reduces to a first degree equation with one unknown, use of the Powell's minimization algorithm tends to overdetermine the system allowing therefore to provide a solution (e.g., the best fitting solution) when no valid solution was found by the LS procedure. As a comparison, Yang et al. [14] reported a failure rate of about 2e3% for a two tilted pyranometers system solved by least squares method (LS) using the PER transposition model and
Table 3 Absolute (in W m2) and relative (in %) RMSE and MBE indexes for the single-sensor approach (one pyranometer). Results for the three pyranometer's mounting planes are provided for the four considered transposition models (see Appendix B) coupled with the OLS decomposition model (see Eq. (20)). Validation statistics are computed from 5117 data points (e.g., 87.98% of the total available data points). Transposition þ Decomposition Models
Pyranometer's mounting plane Failure rate (%) Plane 1 Plane 2 Plane 3 RMSE (W m2 ) Plane 1 Plane 2 Plane 3 rRMSE (%) Plane 1 Plane 2 Plane 3 MBE (W m2 ) Plane 1 Plane 2 Plane 3 rMBE (%) Plane 1 Plane 2 Plane 3
LIU þ OLS
HAY þ OLS
SKA þ OLS
PER þ OLS
1.61 0.27 0.32
0.01 0.18 2.43
0.01 0.18 2.43
4.03 2.81 3.08
33.88 85.83 59.14
19.03 85.76 65.29
19.20 83.61 63.54
25.44 80.16 61.70
12.15 30.78 21.21
6.82 30.75 23.41
6.89 29.98 22.78
9.12 28.74 22.13
14.85 37.45 6.46
0.90 39.18 10.55
4.15 35.51 14.71
5.63 32.72 21.27
5.32 13.43 2.32
0.32 14.05 3.78
1.49 12.73 5.28
2.02 11.73 7.63
sensor NE facing increased from 4% for a tilt angle of 10 e23% for a tilt angle of 40 . 4.2. Two tilted pyranometers Performance of the inverse transposition using two different tilted global irradiance measurements as input to the Gh computation are presented in Table 4 for the 1499 data points that have succeeded for all models as well as for both the EM method and the LS procedure. For each transposition model, Table 4 provides results obtained for both the quadratic formula (in brackets) and the EM method (i.e., Eq. (19)) to solve the quadratic equation. Note that for the LIU isotropic model, values in brackets are from the direct
Table 4 Absolute (in W m2) and relative (in %) RMSE and MBE indexes for a two pyranometers system. Results are provided for the four considered transposition models (see Appendix B). Validation statistics are computed from 1499 data points (e.g., 25.77% of the total available data points) for both the EM method and the LS procedure. Results of the LS procedure are given in brackets. Pyranometer's System
LIU
HAY
SKA
PER
2.03 (8.19) 0.74 (2.08) 0.67 (8.83)
2.01 (33.23) 0.74 (5.30) 0.67 (21.59)
3.34 (34.30) 3.90 (5.35) 0.68 (22.66)
30.88 (44.45) 36.97 (20.92) 26.44 (40.64)
2 3 3
13.89 (34.09) 11.04 (20.09) 8.16 (19.92)
13.89 (29.25) 11.04 (17.73) 8.16 (16.34)
12.51 (44.19) 10.08 (21.34) 6.44 (24.94)
13.20 (60.69) 9.49 (38.04) 7.09 (44.04)
2 3 3
12.14 (29.77) 9.65 (17.55) 7.13 (17.40)
12.14 (25.55) 9.65 (15.49) 7.13 (14.28)
10.92 (38.60) 8.80 (18.64) 5.62 (21.78)
11.53 (53.01) 8.29 (33.22) 6.19 (38.47)
2 3 3
7.80 (1.91) 5.90 (0.24) 5.35 (0.53)
7.80 (0.67) 5.90 (-0.45) 5.35 (-0.39)
3.11 (5.67) 0.06 (6.23) 0.05 (6.75)
4.47 (14.79) 2.09 (14.33) 1.04 (15.39)
2 3 3
6.81 (1.67) 5.15 (0.21) 4.67 (0.46)
6.81 (0.58) 5.15 (-0.40) 4.67 (-0.34)
2.72 (4.95) 0.05 (5.45) 0.04 (5.89)
3.90 (12.92) 1.83 (12.51) 0.91 (13.44)
Failure rate (%) Plane 1 þ Plane 2 Plane 1 þ Plane 3 Plane 2 þ Plane 3 RMSE (W m2 ) Plane 1 þ Plane Plane 1 þ Plane Plane 2 þ Plane rRMSE (%) Plane 1 þ Plane Plane 1 þ Plane Plane 2 þ Plane MBE (W m2 ) Plane 1 þ Plane Plane 1 þ Plane Plane 2 þ Plane rMBE (%) Plane 1 þ Plane Plane 1 þ Plane Plane 2 þ Plane
Transposition Model
C. Housmans et al. / Renewable Energy 113 (2017) 1529e1538
1535
Table 5 Absolute (in W m2) and relative (in %) RMSE and MBE indexes for a three pyranometers system. Results are provided for the four considered transposition models (see Appendix B). Validation statistics are computed from 4066 data points (e.g., 69.91% of the total available data points) using the EM method. Transposition Model
Failure rate (%)
RMSE (Wm2)
rRMSE (%)
MBE (W m2)
rMBE (%)
LIU HAY SKA PER
0.0 0.0 0.0 30.09
20.03 20.06 19.55 21.56
7.93 7.94 7.74 8.53
9.68 9.70 5.94 9.04
3.83 3.84 2.35 3.58
found that the isotropic approach had guaranteed solution with their tropical data. Globally, use of the EM method instead of the LS procedure to solve the quadratic equation decreases by a factor of two and even more in some angular configurations the RMSE in the Gh estimation. As an example, the rRMSE reported for the PER model in the Plane 1 þ Plane 2 conversion is reduced from 53% with the LS solution to 11.5% with the EM approach. For the bias, the benefit of using the EM approach instead of the LS procedure is less apparent in the sense that in one hand the magnitude of the bias is reduced for the more elaborate SKA and PER anisotropic transposition models while in the other hand the bias magnitude is enlarged for both the LIU isotropic model and the simple HAY anisotropic transposition model. To this respect, it is interesting to note that the LIU and HAY models behave similarly in terms of failure rate, RMSE and MBE when using the EM method. In overall, the accuracy of the tilt to horizontal conversion is improved when two tilted irradiance measurements are involved in the inverse transposition process. However, the conversion dependence to the angular configuration of the tilted input measurements does not totally vanish and they are still some data points for which the procedure fails to produce a valid estimation of the global horizontal solar irradiance. 4.3. A three tilted pyranometers system Table 5 summarizes the tilt to horizontal conversion errors using a three tilted pyranometers system (as obtained by the EM solving method) for the 4066 data points that have succeeded for all models. Clearly, involving three different tilted irradiance measurements in the conversion process always ensures solutions excepted for the PER model for which the conversion of about one third (i.e., 30.1%) of the data points are still unsuccessful. Globally, all models behave quite similarly in term of RMSE (overall rRMSE of 8%) and present a negative bias (i.e., rMBE ranging from 3.8% to 2.6%). Models performance as a function of the sky conditions and the solar elevation angle, g, are summarized in terms of RMSE and MBE
in Tables 6 and 7, respectively. Table 6 indicates that the models relative accuracy does not vary noticeably as the sky condition move from overcast to clear sky situations (i.e., a rRMSE variation of about 4.6%, 3.5% and 3% is reported for the PER, SKA and the LIU and HAY models, respectively). All models show a negative bias in clear sky conditions (i.e., 0.8 Kt 1.0) where the global radiation is mainly composed of direct radiation. As an illustration, an underestimation of about 43.7 W m2 (6.5%) is reported for the LIU, HAY and SKA models and of 21.1 W m2 (3.2%) for the PER model. In overcast condition (i.e., 0.0 Kt < 0.2) where the diffuse component largely domines, while the global radiation is underestimated by about 3.7 W m2 (4.9%) by the LIU and HAY models, it is in the opposite slightly overestimated by the SKA and PER models (i.e., MBE of 0.5 W m2 or 0.7% and 0.4 W m2 or 0.5%, respectively). If the SKA model is the best performing model in overcast and cloudy conditions (i.e., Kt < 0.4) with reported RMSE values of 4.5 W m2 (6%) and 9.6 W m2 (5.2%), respectively, its performances do not differ from those reported for the LIU and HAY models in partly cloudy, partly clear and clear sky situations (i.e., Kt 0.4) with RMSE values ranging from 30.2 W m2 (5.8%) in partly clear conditions to 58.2 W m2 (8.7%) in clear sky conditions. The PER model exhibits the lowest RMSE scores for the different sky conditions excepted in overcast situation where its performance (i.e., RMSE of 4.6 W m2 or 6.1%) is closely similar to that reported for the SKA model. Unsurprisingly, models accuracy shows a dependence on the solar elevation with an overall RMSE decreasing from 19.9 W m2 (or 14.1%) in the low solar elevation bin (i.e., g 30 ) to 13.5 W m2 (4.4%) in the high solar elevation bin (i.e., g > 60+). While the SKA model exhibits the best performance in term of RMSE irrespective of the solar elevation, Table 7 indicates that the difference in models accuracy within a given solar elevation bin is limited (i.e., 0.65%, 0.84% and 1.49% in the low, middle and high solar elevation bins, respectively). The magnitude of the models bias tends to reduce as the solar elevation angle increase. As an example, the underestimation of the global solar reported in Table 7 for the LIU and HAY models reduce from about 10.1 W m2 (7.2%) in low solar elevation conditions to 8.6 W m2 (2.8%) in high solar elevation conditions. The sign of the bias is even modified in the case of the
Table 6 Absolute (in W m2) and relative (in %) RMSE and MBE indexes for a three pyranometers system as a function of the sky condition defined in Table 1. Results are provided for the four considered transposition models (see Appendix B). The relative indexes are displayed in brackets. Validation statistics are computed from 4066 data points (e.g., 69.91% of the total available data points) using the EM method. Transposition model
RMSE LIU HAY SKA PER MBE LIU HAY SKA PER
Sky condition 0.0Kt < 0.2
0.2Kt < 0.4
0.4Kt < 0.6
0.6Kt < 0.8
0.8Kt 1.0
6.08 6.08 4.48 4.61
11.83 (6.39%) 11.83 (6.39%) 9.60 (5.19%) 12.82 (6.93%)
25.98 25.98 25.73 32.49
30.12 30.19 30.20 30.23
58.19 58.21 58.21 67.51
7.28 (3.94%) 7.28 (3.94%) 0.45 (0.24%) 3.94 (2.13%)
14.85 14.85 13.05 19.92
(8.09%) (8.09%) (5.97%) (6.14%)
3.69 (4.91%) 3.69 (4.91%) 0.51 (0.68%) 0.38 (0.51%)
(8.28%) (8.28%) (8.20%) (10.35%) (4.73%) (4.73%) (4.16%) (6.35%)
(5.76%) (5.77%) (5.77%) (5.78%)
16.42 16.48 16.47 21.85
(3.14%) (3.15%) (3.15%) (4.18%)
(8.70%) (8.70%) (8.70%) (10.09%)
43.66 43.67 43.67 21.12
(6.53%) (6.53%) (6.53%) (3.16%)
1536
C. Housmans et al. / Renewable Energy 113 (2017) 1529e1538
Table 7 Absolute (in W m2) and relative (in %) RMSE and MBE indexes for a three pyranometers system as a function of the solar elevation, g. Results are provided for the four considered transposition models (see Appendix B). The relative indexes are displayed in brackets. Validation statistics are computed from 4066 data points (e.g., 69.91% of the total available data points) using the EM method. Transposition model RMSE LIU HAY SKA PER MBE LIU HAY SKA PER
Solar elevation angle 5+ g < 30+
30+ g < 60+
60+ g 90+
19.70 19.77 19.55 20.47
20.41 20.40 19.59 22.75
(5.45%) (5.45%) (5.23%) (6.07%)
15.37 15.37 10.76 12.37
9.26 9.26 4.50 9.06
(2.47%) (2.47%) (1.20%) (2.42%)
8.62 (2.79%) 8.62 (2.79%) 3.42 (1.11%) 1.41 (0.46%)
(14.00%) (14.04%) (13.89%) (14.54%)
10.08 (7.16%) 10.11 (7.18%) 7.32 (5.20%) 9.07 (6.45%)
(4.98%) (4.98%) (3.49%) (4.01%)
SKA model for which a slight overestimation of the solar radiation of about 3.4 W m2 (1.1%) is reported at high solar elevation conditions while the model underestimates the global radiation by about 7.3 W m2 (5.2%) in low solar elevation conditions. Finally, comparing the performance between the isotropic and anisotropic approaches to the inverse transposition problem as reported in Tables 5e7 indicates that the improvement from using the LIU model to using anisotropic models is marginal in terms of RMSE and MBE and even inexistent in term of failure rate. 5. Conclusions Different approaches to convert global solar irradiance measured on a tilted plane to horizontal have been evaluated using solar radiation data measured on the roof of the radiometric tower of the Royal Meteorological Institute of Belgium in Uccle, Brussels. Up to three sensors have been considered in the horizontal solar irradiance reconstruction. When only a single tilted sensor is used, the conversion can be carried out with a decomposition model coupled with a transposition model to solve the inverse transposition problem. In this case, there is an additional error (additional to the inverse transposition problem) in the predicted horizontal irradiance. In addition, with only one tilted irradiance involved in the inverse modeling approach, the conversion's performance is very sensitive to the tilted pyranometer angular configuration (i.e., the tilt angle and the surface's orientation). None of the considered transposition models was found to outperform over the 3 pyranometers mounting plane configurations. When two (or more) tilted irradiance sensors are involved, only a transposition model is required and the inverse transposition problem can be solved more accurately. In practice, there are certain sun positions for which this procedure still fails to produce a valid estimation of the global horizontal irradiance. Consequently more instruments should be used so as to overdetermine the system. Here it has been shown that three tilted pyranometers set at different orientations is sufficient to guarantee solution for three of the four considered transposition models. Indeed, because of the non bijectivity of the anisotropic transposition model of Perez et al. [23] the conversion of a percentage of data points are unsuccessful. Sensitivity experiments conducted in the case of a two tilted pyranometers system indicate that solving the inverse transposition quadratic equation by minimizing the errors instead of using the least square procedure tends to significantly reduces the percentage of unsuccessful conversion irrespective of the considered transposition model. This is also true for the conversion performed with the isotropic model even if in this case the inverse
transposition problem reduces to a first degree equation with one unknown as the Powell's minimization algorithm overdetermines the system and provides the best fitting solution. Moreover the accuracy of the tilt to horizontal conversion is improved with the errors minimization method irrespective of the considered transposition model. Finally, comparing the performance between the isotropic and anisotropic approaches to the inverse transposition problem in angular configurations similar to those encountered in Belgian PV systems installations (e.g., with tilt angle as great as 50.79 ) indicates that the improvement from using the LIU isotropic model to using anisotropic models is not significant (e.g., a rRMSE of 7.9% and a rMBE of 9.7% are reported for the LIU model when considering a three tilted pyranometers system compared to a rRMSE of 7.7% and a rMBE of 5.9% for the best performing SKA anisotropic model) or even inexistent regarding the percentage of unsuccessful conversion (e.g., a failure rate of 0% is reported for the LIU model in the three tilted pyranometers system compared to a failure rate of 30.1% for the PER model). Acknowledgments This study was supported by the Belgian Science Policy Office (BELSPO) through the BRAIN-be research project:”Solar irradiation from the energy production of residential PV systems eSPIDER00 , contract No. BR/314/PI/SPIDER. Appendix A. Skartveit and Olseth (1987) decomposition model (OLS model) Skartveit and Olseth (1987) estimated the direct normal irradiance, Bn, from the global horizontal irradiance, Gh , and the solar elevation angle, g, for Bergen (Norway, 60.4oN) with the following equation based on hourly records of global and diffuse horizontal irradiances with mean solar elevation larger than 10 during 1965e1979:
Bn ¼
Gh ð1 JÞ sin g
(20)
where J is a function of the clearness index, Kt . The model was validated with data collected in 12 stations worldwide. The function J reads as:
8 1 for Kt < c1 > > i h > > 1=2 < 2 1 ð1 d1 Þ d2 c3 þ ð1 d2 Þ c3 for c1 Kt 1:09c2 J¼ > > > > : 1 1:09c 1 Y for Kt > 1:09c2 2 Kt (21) where: c1 ¼ 0:2 g c2 ¼ 0:87 0:56 0:06 c3 ¼ 0:5 1 þ sin p dc4 0:5 3 c 4 ¼ Kt c 1 0:06 g d1 ¼ 0:15 þ 0:43 e d2 ¼ 0:27 d3 ¼ c2 c1 1=2 Y ¼ 1 ð1 d1 Þ ðd2 cc3 þ ð1 d2 Þ cc23 Þ 4 0:5 cc3 ¼ 0:5 1 þ sin p cc d 3
cc4 ¼ 1:09 c2 c1
C. Housmans et al. / Renewable Energy 113 (2017) 1529e1538
Appendix B. Transposition models
1537
inclined surface. More specifically, the angular location of the circumsolar region is determined by the ratio a/b. They are calculated from the equations of solar geometry:
Each model develops the diffuse transposition factor (i.e., the ratio of diffuse radiation on a tilted surface to that of a horizontal, Rd ) according to specific assumptions.
a ¼ maxð0; cos qi Þ
B.1. LIU model (Liu and Jordan, 1962):
b ¼ maxðcos 85o ; cos qz Þ
The isotropic model (that assumes a uniform distribution of the diffuse radiation over the sky dome) was published in 1962:
The brightness coefficients F1 and F2 are derived from the so-called Perez coefficients:
RLIU d ¼
1 þ cosb 2
(22)
where b denotes the surface tilt angle with respect to the horizontal plane. B.2. HAY model (Hay, 1979): In the HAY model, diffuse radiation from the sky is composed of an isotropic component and a circumsolar one. Horizon brightening is not taken into account. An anisotropy index, FHay , is used to quantify a portion of the diffuse radiation treated as circumsolar with the remaining portion of diffuse radiation assumed to be isotropic, i.e.,
RHAY ¼ FHAY rb þ ð1 FHAY Þ d
1 þ cos b 2
(23)
where FHAY ¼ Bh =ðIo cosqz Þ is the Hay's sky-clarity factor and rb ¼ cosqi =cosqz the beam radiation conversion factor. qz is the solar zenith angle, qi is the incidence angle of the beam radiation on the tilted surface, Bh is the direct horizontal solar irradiance and Io is the extraterrestrial normal incident irradiance. The model reduces to the LIU model (Eq. (22)) for FHAY ¼ 0. B.3. SKA model (Skartveit and Olseth, 1986): Solar radiation measurements indicate that a significant part of sky diffuse radiation under overcast sky conditions comes from the sky region around the zenith. This effect vanishes when cloud cover disappears. Skartveit-Olseth (1986) modified the HAY model (Eq. (23)) in order to account for this effect,
RSKA ¼ FHAY rb þ Zcosb þ ð1 FHAY ZÞ d
1 þ cosb 2
(24)
where Z ¼ maxð0; 0:3 2FHAY Þ is the Skartveit-Olseth's correction factor. If FHAY 0:15, then Z ¼ 0 and the model reduces to the HAY model. B.4. PER model (Perez et al., 1987): Compared to the previously described transposition models, the model proposed by Perez et al. (1987) represents a more detailed analysis of the isotropic diffuse, circumsolar and horizon brightening radiation by using empirically derived coefficients. According to this model,
a 1 þ cosb þ F2 sinb RPER ¼ F1 þ ð1 F1 Þ d b 2
(25)
where F1 and F2 are sky brightness coefficients for the circumsolar region and the region above the horizon line, respectively. Note that if F1 ¼ F2 ¼ 0, it reduces to the LIU model (Eq. (22)). The coefficients a and b take into account the angle of incidence of the sun onto the
F1 ¼ F11 ðεÞ þ F12 ðεÞD þ F13 ðεÞqz F2 ¼ F21 ðεÞ þ F22 ðεÞD þ F23 ðεÞqz where the Perez coefficients Fij are function of the sky clearness parameter ε and the sky brightness parameter D. These factors are defined by:
ε¼
ðDh þ Bh Þ=Dh þ 1:041q3z
D¼m
1 þ 1:041q3z Dh Dh ¼ Io Io cosqz
where Dh is the diffuse horizontal solar irradiance, m is the optical air mass and qz is in radians. Many sets of Perez coefficient values have been determined by different studies. In this paper, we applied the set of coefficients from Perez et al. (1990).
References [1] M. Wali, F. Evrendilek, T. West, S. Watts, D. Pant, H. Gibbs, et al., Assessing terrestrial ecosystem sustainability: usefulness of regional carbon and nitrogen models, Nat. Resour. 35 (4) (1999) 20e33. [2] M. Trnka, Z. Zalud, J. Eitzinger, M. Dubrovsky, Global solar radiation in central european lowlands estimated by various empirical formulae, Agric. For. Meteorol. 131 (2005) 54e76. e, C. Bertrand, Quality control of solar radiation data within the rmib [3] M. Journe solar measurements network, Sol. Energy 85 (2011) 72e76. [4] R. Perez, R. Seals, R. Stewart, A. Zelenka, V. Estrada-Cajigal, Using satellitederived insolation data for the site/time specific simulation of solar energy systems, Sol. Energy 53 (1994) 491e495. [5] R. Perez, R. Seals, A. Zelenka, Comparing satellite remote sensing and ground network measurements for the production of site/time specific irradiance data, Sol. Energy 60 (1997) 89e96. [6] A. Hammer, D. Heinemann, C. Hoyer, R. Kuhlemann, E. Lorentz, R. Mueller, et al., Solar energy assessment using remote sensing technologies, Remote Sens. Environ. 86 (2003) 423e432. , R. George, S. Wilcox, T. Stoffel, D. Myers, D. Heimiller, Solar Resource [7] D. Renne Assessment. Technical Report Nrel/tp-581-42301, National Renewable Energy Laboratory, 2008, 1617 Cole Boulevard, Golden, Colorado:80401e3393. [8] D. Yang, Solar radiation on inclined surfaces: corrections and benchmarks, Sol. Energy 136 (2016) 288e302. [9] D. Yang, Z. Dong, A. Nobre, Y. Khoo, P. Jirutitijaroen, W. Walsh, Evaluation of transposition and decomposition models for converting global solar irradiance form tilted surface to horizontal in tropical regions, Sol. Energy 97 (2013) 369e387. [10] B. Marion, A model for deriving the direct normal and diffuse horizontal irradiance from the global tilted irradiance, Sol. Energy 122 (2015) 1037e1046. [11] R. Perez, P. Ineichen, E. Maxwell, R. Seals, A. Zelenka, Dynamic global-to-direct irradiance conversion models, ASHRAE Trans.-Res. Ser. 3578 (1992) 354e369. [12] F. Olmo, J. Vida, I. Foyo, Y. Castro-Diez, L. Alados-Arboledas, Prediction of global irradiance on inclined surfaces from horizontal global irradiance, Energy 24 (1999) 689e704. [13] S. Killinger, et al., Projection of power generation between differentlyoriented pv systems, Sol. Energy 136 (2016) 153e165. [14] D. Yang, Z. Ye, A. Nobre, H. Du, W. Walsh, L. Lim, et al., Bidirectional irradiance transposition based on the perez model, Sol. Energy 110 (2014) 768e780. e, C. Bertrand, Evaluation of different models to estimate [15] C. Demain, M. Journe the global solar radiation on inclined surfaces, Renew. Energy 50 (2013) 710e721. e, C. Bertrand, Corrigendum to ”evaluation of different [16] C. Demain, M. Journe
1538
[17]
[18] [19] [20] [21] [22] [23]
[24] [25] [26]
[27] [28]
C. Housmans et al. / Renewable Energy 113 (2017) 1529e1538 models to estimate the global solar radiation on inclined surfaces”[renewable energy 50 (2013) 710-721], Renew. Energy 101 (2017) 1401e1403. e, Evaluation of decomposition models C. Bertrand, G. Vanderveken, M. Journe of various complexity to estimate the direct solar irradiance over Belgium, Renew. Energy 74 (2015) 618e626. C.A. Gueymard, Direct and indirect uncertainties in the prediction of tilted irradiance for solar engineering applications, Sol. Energy 83 (2009) 432e444. A. Skartveit, J. Olseth, A model for the diffuse fraction of hourly global radiation, Sol. Energy 38 (1987) 271e274. B. Liu, R. Jordan, Daily insolation on surfaces tilted towards the equator, Trans. ASHRAE 67 (1962) 526e541. J. Hay, Study of Shortwave Radiation on Non-horizontal Surfaces. Rep No 7912, Atmospheric Environment Service, Downsview, Ontario, 1979. A. Skartveit, J. Olseth, Modelling slope irradiance at high latitudes, Sol. Energy 36 (1986) 333e344. R. Perez, R. Seals, P. Ineichen, R. Stewart, D. Menicucci, A new simplified version of the perez diffuse irradiance model for tilted surfaces, Sol. Energy 39 (1987) 221e231. D. Faiman, A. Zemel, A. Zangvil, A method for monitoring insolation in remote regions, Sol. Energy 38 (1987) 327e333. C. Grosan, A. Abraham, A new approach for solving nonlinear equations systems, IEEE Trans. Syst. Man. Cybern. Part A Syst. Humans 38 (2008) 698e714. M. Powell, An efficient method for finding the minimum of a function of several variables without calculating derivatives, Comput. J. 7 (2) (1964) 155e162. F. Acton, Numerical Methods that Work, 1990, corrected edition, Mathematical Association of America, Washington, 1970, pp. 464e467. W. Press, S. Teukolsky, W. Vetterling, B. Flannery, Numerical Recipes in C. The Art of Scientific Computing (Reprinted with corrections, 1997), second ed., Cambridge University Press, Cambrige, UK, 1992, p. 994.
€ppen-geiger [29] M. Peel, B. Finlayson, T. McMahon, Updated world map of the ko climate classification, Hydrol. Earth Syst. Sci. 11 (2007) 1633e1644.
Nomenclature Io : extraterrestrial normal incident irradiance ðW m2 Þ Gt : tilted global solar irradiance ðW m2 Þ Gh : global horizontal irradiance ðW m2 Þ Bt : in-plane direct irradiance ðW m2 Þ Bn : direct normal irradiance ðW m2 Þ Bh : direct irradiance on a horizontal surface ðW m2 Þ Dt : in-plane diffuse irradiance ðW m2 Þ Dh : diffuse horizontal irradiance ðW m2 Þ Dg : irradiance due to the ground reflection ðW m2 Þ rb : factor that accounts for the direction of beam radiation Rd : diffuse transposition factor Rr : transposition factor for ground reflection Kt : effective global horizontal transmittance (or clearness index) Kn : direct normal transmittance Kd : diffuse horizontal transmittance g: solar elevation angle (rad) a: azimuth of the tilted plane (rad) b: plane tilt angle (rad) qi : incidence angle (rad) qz : zenith of the Sun (rad) Р: ground albedo ε: sky clearness parameter Fij : Perez model coefficient Z: Skartveit-Olseth's correction factor