Stochastic simulation model of hourly total solar radiation

  • Published on

  • View

  • Download


  • Solar Energy Vol. 37, No. 2, pp. 119-126, 1986 0038-092X/86 $3.00 + .00 Printed in the U.S.A. 1986 Pergamon Journals Ltd.


    A. BALOUKTSIS and PH. TSALIDES Department of Electrical Engineering, School of Engineering, Democritus University of Thrace,

    67100 Xanthi, Greece

    (Received 17 April 1985; revision received 24 February 1986)

    Abstract--A stochastic simulation model of hourly global solar radiation is presented in this paper. It is developed by introducing the concept of "time dependent frequency distribution" (TDFD) of hourly insolation values. In this model the two most critical aspects of time series simulation, i.e., the reproduced time series values which have the appropriate time dependent frequency distribution for the parameter being simulated and the correlation between successive values, are taken into ac- count. The elimination of the TDFD of the data and the transformation of the data distribution to a Gaussian distribution (required for the stationary time series analysis) were carried out using a mapping technique. The autocorrelation function of the transformed data showed that the produced time series is stationary. Then, an antimapping coefficient matrix is developed, which provides a simple yet an effective simulation device. The described model has been applied in Athens (Greece) where hourly insolation data covering a period of two years are used. The theoretical results obtained using this simulation model, regarding both the TDFD and the correlation, are in agreement with the measured data.

    1. I N T R O D U C T I O N

    In the design of Photovoltaics (PV) systems it is desirable to simulate their performance on an hour to hour basis. Therefore, there is a need for an ap- propriate stochastic simulation model of hourly global solar radiation in order to generate realistic operating conditions.

    Previous studies carried out by Baker et al.[1] and by Bennett[2] on weekly and daily solar radia- tion models have shown that the distribution is neg- atively skewed and occasionally bimodal. It has also been shown by Engels et al.[3] that the prob- ability distribution of hourly irradiation also pos- sesses these characteristics. Therefore, because of the non-Gaussian nature of the distributions, the works of Brinkworth[4] and Sfeir[5] are not the most appropriate in the design of solar conversion systems. Moreover, the seasonal variation of ir- radiation (TDFD) and the correlation between suc- cessive hourly irradiation values are very important in the optimum design of PV systems. The Markov transition matrix developed by Mustacchi et al.[6] does not take into account this seasonal variation. Also, the mathematical model developed by Exell[7] does not take into account the correlation, and therefore, both the models may not be the most suitable in the optimum design of PV systems.

    To incorporate the seasonal variation of solar irradiation as well as the correlation between suc- cessive hourly irradiation values in the optimum de- sign of PV systems, a stochastic simulation model of hourly global solar radiation is developed in this paper. In this model the two most critical aspects of time series simulation, (a) the reproduced time series values which have the appropriate TDFD for

    the parameter being simulated and (b) the corre- lation between successive hourly values, are taken into account.

    2. D E S C R I P T I O N OF T H E MODEL

    2.1 Annual and daily periodicity removal It is well known that hourly data represent a pe-

    riodicity during the day as well as a periodicity throughout the year for each hour of the day. The Fourier coefficientss of the annual periodicity cor- responding to hth hour are given[8] by minimising the quantity


    2 d = l

    [Ia h - {ao h + a h cos(2~r(d - 1)/365)

    + b{' sin(2~r(d - 1)/365)}] 2 (1)

    where a h, al h, bl h are the Fourier coefficients cor- responding to hth hour of the day, and they are given by


    a0 h = 2/365 ~ Ia h, d = l


    a h = 2/365 ~ I h cos(2"rr(d- 1)/365), d = l


    bl h = 2/365 ~ I h sin(2~r(d - 1)/365). d = l


    Ia h is the measured hourly global solar radiation dur- ing the hth hour of day d of the year (in kJ/m2), and h = 1, 2, . . . , k (1 corresponding to sunrise hour



    and k corresponding to sunset hour. k is taken to be equal to the max. number of hours of a day from sunrise to sunset, assumed to be constant through- out the year).

    Once the above Fourier coefficients have been calculated, both the periodicities are then removed by simply subtracting the quantity (Ia h) from luh, i.e.,

    Ha h = I,~ - (Ia h) (3)

    where (Iuh) is given by

    (I~) = ah/2 + a h cos(Z~r(d - 1)/365)

    + b~ h sin(2~r(d - 1)/365). (4)

    2.2 The m a p p i n g a n d a n t i m a p p i n g t echn ique The time series of the transformed data Hu h , after

    the removal of the two periodicities, is not yet sta- tionary, because the time dependence of the fre- quency distribution of solar radiation still exists. Therefore, there is a need for an appropriate trans- formation which leads to a Gaussian distribution (constant throughout the year) required for the sta- t ionary time series analysis. This transformation was carried out assuming that the frequency dis- tribution of solar radiation is almost constant for the period of one month (the variation being very small, as it has also been verified in the example), and thus it was considered that twelve monthly fre- quency distributions describe very satisfactorily the above mentioned TDFD. In the process of the transformation a mapping technique has been used, as follows: The range of values of the transformed data HUh for the ith month was divided into l inter- vals and then the cumulative frequency distribu- tion, CFD(x~,), was obtained, where x~ corresponds to the midpoint of the nth interval of the ith month. Subsequently, the discrete values y~, of a new sta- tistical variable uuh were determined[9]:

    Fv(y I - q (Y~,) q) = CFD(xl,) (5)

    where i = 1, 2 . . . . . 12 and

    F v ( y ~ ) - F v ( - q ) = CFD(x~,) (6) F v ( q ) - F u ( - q)

    where Ftz(y~) is the standard normal distribution (with zero mean and unit variance). From eqn (6), it follows that

    F v ( y ~ ) = F u ( - q ) + CFD(x~)

    x [ F v ( q ) - F u ( - q ) ]

    1 - F u ( y i ) = 1 - F v ( - q ) - CFD(x~,) (7)

    x [ F v ( q ) - F v ( - q ) ]

    1 - F u ( y ~ ) = F v ( q ) - CFD(x~,)

    [ 2 F v ( q ) - l].

    and PH. TSALIDES

    The discrete points y~, are given[10] by:

    go + g~t + g2t 2 . = - + e t p ) , (8) VI, l 1 + g3 t + g4/2 + gs t 3

    where e (p ) is the approximation error, less than 4.5 x 10-4; t = [ln(I/p2)]v2; and

    go = 2.515517 g3 = 1.432788

    gj = 0.802853 g4 = 0.189269

    g2 = 0.010328 g5 = 0.001308.

    Ifp~ = F u ( q ) - CFD(xl,) [2Fu(q) - 1], then

    P = Pn, if Pn ~< 0.5, and

    p = 1 - p~, y~, = -y~,, ifP~ > 0.5.

    The points y~ thus obtained have a one-to-one correspondence with i xn. The pairs (x~, y~,) were ap- proximated by a polynomial MAX-MIN method[ l l ] . The transition from Ha h to u h is achieved with the mapping polynomial, and that from uuh to HUh with the antimapping polynomial, re- spectively.

    The limits - q, q in eqn (5) were selected in order to achieve the greatest accurancy in the polynomial approximation of the discrete points of the mapping and antimapping technique.

    After that, the mapping and antimapping coef- ficient matrices (12, j) were developed (where j is the coefficient number of the polynomials), to pro- vide the transition from h h Ha to Ud and vice versa.

    The new variable u h follows a Gaussian distri- bution (constant throughout the year) and its au- tocorrelation function, ac f (uh ) , is shown in Fig. 1. Therefore, the stochastic process of the variable u h is a weakly stationary Gaussian random process and, hence, the method (see Section 2.3) of cal- culating its spectra as well as the method (see Sec- tion 2.4) of simulation described below can be used.

    2.3 Calcu la t ion o f spec t ra The spectra was calculated from the transformed

    data u h using the truncated autocovariance function method[12] with the Tukey-Hanning window and an appropriate truncation point M (M has been taken to be equal to 50 in the example described in Section 3). The spectra was determined from the autocorrelation coefficients rather than by direct Fourier analysis of uuh. These spectra estimates were smoothed by a Direct Search least-squares method[13]. The calculated spectra, S(to), (of the example) and the corresponding least-squares fit- ting curve are shown in Fig. 2.

    2.4 S i m u l a t i o n A comprehensive method of simulating a sta-

    t ionary random process, such as uuh, with a spectral density function S(co) specified in the region I co I <

  • Stochastic simulation model of hourly total solar radiation 121

    (J t~ olll, I l l l l l l l l l l , , , ] l J l , , . . . . . . . . . . . . . . . . . .




    Fig. I. Autocorrelation function of the normally distributed variable u~L



    i 0

    + calculated estimates o fitting curve

    o +

    +o o

    ~ " 9 0 9 r O e O , e e e e 6 ~ ~ 6 ~ ~ 8 n/2 n


    Fig. 2. Spectra of u~ ( + ) and the corresponding least-squares fitting curve (O).

    ~r, is described by Shinozuka[14]. Defining Ato = ~r/N, the random process can be simulated by a se- ries of cosine waves at almost evenly spaced fre- quencies and with amplitudes weighted according to the spectral energy at the wave frequency, as follows:


    u h = (2) 1/2 ~ [2S(toi)Ato] 1/2 cos(to~t + qoi) (10) i = l


    toi = [i - (1/2)]Ato ~oi: independent random phase angles uniformly

    distributed between 0 and 2~r to" = toi + ~to

    and Bto is a small random frequency introduced to avoid the periodicity of the simulated process, uni- formly distributed between -Ato ' /2 and Aa~'/2, with Ato' < A~.

    The simulated process is Gaussian by virtue of the central limit theorem. In the present work, N has been taken to be equal to 20 and Ato' equal to Ato/20.

    Once u h has been simulated, the quantity H h is then calculated using the antimapping coefficient matrix. Finally, the hourly values of global irradia- tion are obtained by simply adding the quantity (I h) to H h, i.e.,

    I h = H h + (I h) (11)

    The process of simulation is shown in the flow- chart of Fig. 3.

    3. EXAMPLE

    The model described in the previous sections has been applied in Athens (Greece), where hourly global insolation data covering the period of two years from 1980 to 1981 are used[15, 16].

    Figure 4 shows the quantity (lh). It is obvious


    h u d

    simulation (Shinozuka's method)

    h u d

    -[ d=1 to 365 [ I

    selection of the month (i=I .... 12)

    -] h=1 to K(i) I I

    random numbers


    j=ll3 " coefficients cij (i:I , . .12, j:1 .... m)

    Fig. 3. F|ow-chart of the simulation process (k(i): hours/day for each month).

  • Z

    Stochastic simulation model of hourly total solar radiation



    Fig. 4. A 2-D plot representing the annual and the daily periodicity of hourly global insolation data.

    3 . 5

    C='O 0

    - 3 . 5 - 2 7 2 0

    / .!

    // / / / /'

    / 4"

    #f ~ ~ ~ ~--" I'/ | ! ,J

    0 1 0 4 6 h

    Hd ( k J / m 2 )

    Fig. 5. Correspondence between H h (after the removal of periodicities) and the normally distributed variable u h in Athens (Greece).

    , A P R I L

    - - - J U L Y

    . . . . . D E C E M B E R





    o o

    A A o

    o A


    m 0 0 A A A

    A O o


    0 6 O 0 0 0 0 I J l

    0 1046 - 2 7 2 o

    Fig. 6. Measured (O) and simulated (&) frequency distributions covering a period of one year in Athens (Greece).



    : E " m

    ..s cr o


    0 - 2 5 1 0

    0 0

    A A

    A A 0

    A o

    , ,s68g 0

    H : ( k j / m 2 )

    o: DATA




    Fig. 7. Measured (O) and simulated (A) frequency distributions for April in Athens (Greece).

    c ._o r.~



    .35 o, DATA

    o A, S I M U L A T I O N o


    A A



    A A


    0 a a A 6 # ~ 6 6 6 6 6 6 6 - 2 7 2 0 0 418

    Hdh( kJ/m 2)

    Fig. 8. Measured (O) and simulated (/x) frequency distributions for July in Athens (Greece).

  • Stochastic simulation model of hourly total solar radiation 125

    35 o,DATA

    A , S I M U L A T I O N


    f o ~ A O

    o ,~ o L. A A 6 A A ,, AoA 8oO 6

    A A O O O

    0 t I 6 ,

    - 1 1 7 2 0 711 H : ( k J l m 2 )

    Fig. 9. Measured () and simulated (A) frequency distributions for December in Athens (Greece).

    that this quantity represents a 2-D digital low-pass filter.

    The variable q in eqn (5) has been taken equal to 3.5 throughout the year (greater values must be selected for the summer months and smaller for the winter months), a value which is expected to be the same for any other location with similar distribution forms of solar irradiation. Figure 5 shows plots of the correspondence between u h and Hhe in Athens for three different months, April, July and Decem- ber.

    The following anti-mapping matrix, C(i, j ) , (i = 12, j = 7) provides a simple yet an effective simu- lation device:

    C(i, 31 =

    - 113.45 454.58 -58.544 201.62 305.27 - 100.31

    -33.059 643.94 - 118.89 149.81 404.33 -204.55 147.43 171.36 - 195.3

    -29.042 -41.973 -39.211 - 172.03 -29.879 296.95 -46.827 36.365 123.16

    56.996 220.7 -213.84 230.66 363.02 -273.97 75.074 297.99 -80 .639

    -65.532 426.25 - 116.84

    shown in Fig. 6. Figures 7, 8 and 9 show three fre- quency distributions of the sequence HI) covering the period of the months April, July and December, respectively, (which are considered to be the most representative months of the year) and the corre- sponding frequency distributions after the simula- tion. It is obvious that the agreement between the actual data and the theoretical results is very good in all cases. The results for the other months of the year were found to be similar.

    Furthermore, it is obvious from Figs. 7, 8 and 9 that the described model is fully cdmpatible with the concept of the "t ime dependent frequency dis- t r ibution" (TDFD) introduced in this paper.

    -23.183 5.0635 0.6277 -(I.1251 35.528 5.8586 -2.9293 -0.0251

    -8.1183 9.29 -0 .5859 -0 .1988 52.141 17.115 -3 .8499 -0 .4938 119.56 8.3276 -7 .8672 0.0749 145.54 -24.271 -8 .9134 1.8413 126.42 - 112.57 -6 .8629 6.6955 54.82 - 38.918 -2 .5108 2.1342 69.131 18.747 -4.8961 -(I.544 52.853 33.603 -4 .3102 -1.3391 11.55 4.1429 -0 .8788 0.0083

    - 13.935 17.785 0.036 -0 .7784

    The autocorrelation function of the normally dis- tributed variable uha, acf(uha), is shown in Fig. 1, where the stationarity of the time series is apparent. Its power spectra and the corresponding least- squares fitting curve are shown in Fig. 2.

    Hourly global radiation sequences covering a pe- riod of one year have been obtained using the simu- lation model shown in Fig. 3. The actual frequency distribution of the sequence H,~ and that corre- sponding to the same sequence after the simulation (by means of the model described in this paper) are


    A simple yet an effective stochastic simulation model of hourly global solar radiation is presented in this paper

    The introduction of the " t ime dependent fre- quency distribution" concept leads to a reliable simulator

    It is also remarkable that the spectra of the nor- mally distributed variable u,~ is expected to be al- most the same for any other location with similar climatological conditions.


    The t ime series thus obtained exhibits all the sto- chastic features (correlation, seasonal variation) of the solar irradiation of a given location.

    Acknowledgements--We wish to thank Prof. A. Thanai- lakis for the stimulating discussions and his helpful com- ments. We also wish to thank Professors J. Antonopoulos, E. Galoussis and D. Tsanakas for their interest in the preparation of this work.


    1. D. G. Baker and J. C. Klink, Solar radiation reception, probabilities, and areal distribution in the north-cen- tral region. North Central Reg. Res. Publ., p. 225. Agricultural Experiment Station, University of Min- nesota (1975).

    2. I. Bennett, Frequency of daily insolation in Anglo North America during June and December. Solar En- ergy 11(1), 41-55 (1967).

    3. J. Dann Engels, S. M. Pollock and J. A. Clark, Sto- chastic modelling of hourly terrestrial solar radiation. Prelim. Rep., Office of Energy Research, University of Michigan, Ann Arbor, MI (Jan. 1980).

    4. B. J. Brinkworth, Autocorrelation and stochastic modelling of insolation sequences. Solar Energy 19, 343-347 (1977).

    5. A. A. Sfeir, A stochastic model for predicting solar

    system performance. Solar Energy 25(2), 149-154 (1980).

    6. C. Mustacchi, V. Cena and M. Rocchi, Stochastic simulation of hourly global radiation sequences. Solar Energy 23, 47-51 (1979).

    7. R. H. B. Exell, A mathematical model for solar ra- diation in South-East Asia (Thailand). Solar Energy 26, 161-168 (1981).

    8. P. Bloomfield, Fourier Analysis of Time Series. John Wiley, New York (1976).

    9. A. Papoulis, Probability, Random Variables and Pro- cesses. McGraw-Hill, New York (1965).

    10. M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions. Dover Publications, New York (1972).

    11. F. Scheid, Theory and Problems of Numerical An- alysis. Schaum's Outline Series. McGraw-Hill, New York (1968).

    12. C. Chatfield, The Analysis of Time Series: Theory and Practice. Chapman and Hall Ltd., London (1975).

    13. D. M. Himmelblau, Applied Non-Linear Program- ming. McGraw-Hill, New York (1972).

    14. M. Shinozuka and C. M. Jan, Digital simulation of random processes and its applications. J. of Sound and Vibration 25(1), 111-128 (1972).

    15. National Observatory of Athens Meteorological In- stitute, Climatological Bulletin (1980).

    16. National Observatory of Athens Meteorological In- stitute, Climatological Bulletin (1981).


View more >