Photovoltaic Model

This page documents the functionality of the PV model available in SIMONA.

The initial parts of the model are presented in the paper Agent based approach to model photovoltaic feed-in in distribution network planning by Seack, Kays, and Rehtanz [1]. Since then several adaptions has been made that are documented as follows.

The PV Model is part of the SIMONA Simulation framework and represented by an agent.

Parameters

Attributes, Units and Remarks

Please refer to PowerSystemDataModel - PV Model for Attributes and Units used in this Model.

Implemented Behaviour

Sequence Diagram Behaviour PV Model

Calculations

The energy produced by a photovoltaic (PV) unit in a specific time step is based on the diffuse and direct radiance provided by the used weather data. The following steps are done to calculate (= estimate) the power feed in by the PV.

The model calculations are performed on the basis of radiance (power per area, symbol \(G\)). Some sources are primarily based on radiance [Myers.2017]. Others are developing models using radiation values (energy per area, symbols \(I\) or \(H\)), which are often times referencing a fixed time frame, e.g. of one hour. Conversions to radiance values can be easily made and can be used interchangeably ([Duffie.2013] p. 86, footnote).

Irradiance and irradiation describe power and energy arriving at a surface, while the terms radiance and radiation are used for general purposes ([Iqbal.1983] Ch. 2.6).

To calculate the overall feed in of the PV unit, the sum of the direct irradiance, diffuse irradiance and reflected irradiance is needed. In the following, the formulas to calculate each of these radiances are presented and explained. The sections end with the formula to calculate the corresponding power feed in.

The surface azimuth angle \(\alpha_e\) starts at negative values in the East and moves over 0° (South) towards positive values in the West (Source).

Declination Angle

The declination angle \(\delta\) (in radian!) is the day angle that represents the position of the earth in relation to the sun. To calculate this angle, we need to calculate the day angle \(J\). The day angle in radian is represented by:

\[ J = 2 \pi(\frac{n-1}{365}) \]

with
\(n\) = number of the day in the year (e.g. 1 January = 1, 20 February = 51)

Based on \(J\) the declination angle \(\delta\) (in radian!) can be calculated as follows:

\[\begin{split} \begin{aligned} \delta = 0.006918 - 0.399912 \cdot \cos(J) + 0.070257 \cdot \sin(J) \\ - 0.006758 \cdot \cos(2\cdot J) + 0.000907 \cdot \sin(2 \cdot J) \\ - 0.002697 \cdot \cos(3 \cdot J) + 0.00148 \cdot \sin(3 \cdot J) \end{aligned} \end{split}\]

References:

Hour Angle

The hour angle is a conceptual description of the rotation of the earth around its polar axis. According to Watter [Watter.2013], it is calculated (in radian!) as follows, with positive values in the morning:

\[ \omega = ((12 - T_\mathrm{sol}) \cdot 15) \cdot (\frac{\pi}{180}) \]

Since outside German literature the hour angle is defined as negative in the morning, arrives at 0° at noon (solar time) and ends with a positive value in the evening, we use the following adaption:

\[ \omega = ((T_\mathrm{sol} - 12) \cdot 15) \cdot (\frac{\pi}{180}) \]

with
\(T_\mathrm{sol}\) = local solar time
Note: The sun crosses the meridian of the observer at local noon (12:00 in local solar time), taking into account the eccentricity of earth’s orbit.

\[ T_\mathrm{sol} = T_\mathrm{mean} + \frac{T_\mathrm{ET}}{60} \]

with
\(T_\mathrm{mean}\) = local mean time, follows the average movement of the sun over a year.
\(T_\mathrm{ET}\) = equation of time (in Minutes, thus divided by 60), difference between \(T_\mathrm{mean}\) and \(T_\mathrm{sol}\) due to the eccentricity of the earth’s orbit.

\[ T_\mathrm{mean} = H_{\mathrm{UTC}} + \frac{M_{\mathrm{UTC}} + 4 \cdot (\lambda_{\mathrm{std}} + \lambda)}{60} \]

with
\(H_{\mathrm{UTC}}\) = Hour of the time basis (we choose UTC as time basis)
\(M_{\mathrm{UTC}}\) = Minute of the time basis (we choose UTC as time basis)
\(\lambda_{\mathrm{std}}\) = Standard meridian for the time basis (since we choose UTC as time basis, \(\lambda_{\mathrm{std}}\) is zero)
\(\lambda\) = longitude of the location (in Degrees) of the PV panel

Note: The formula in Duffie [Duffie.2013] p. 11 is defining λ in degrees west being positive. However, since we are using the longitude directly, which is negative westwards, this changes the sign in the formula.

\[\begin{split} \begin{aligned} T_\mathrm{ET} = 0.0066 + 7.3525 \cdot \cos(J + 1.4992378274631293) \\ + 9.9359 \cdot \cos(2 \cdot J + 1.9006635554218247) \\+ 0.3387 \cdot \cos(3 \cdot J + 1.8360863730980346) \end{aligned} \end{split}\]

with
J = day angle (in radian!)

References:

Sunrise Angle

The hour angles at sunrise and sunset are very useful quantities to know. These two values have the same absolute value, however the sunset angle (\(\omega_{\mathrm{SS}}\)) is positive and the sunrise angle (\(\omega_{\mathrm{SR}}\)) is negative. Both can be calculated from:

\[ \omega_{\mathrm{SS}}=\cos^{-1}(-\tan (\phi) \cdot \tan (\delta)) \]
\[ \omega_{\mathrm{SR}}=-\omega_{\mathrm{SS}} \]

with
\(\delta\) = the declination angle
\(\phi\) = observer’s latitude

References:

Solar Altitude Angle

Represents the angle between the horizontal and the line to the sun, that is, the complement of the zenith angle.

\[ \sin(\alpha_s) = \sin (\phi) \cdot \sin (\delta) + \cos (\delta) \cdot \cos (\omega) \cdot \cos (\phi) \]

with
\(\delta\) = the declination angle
\(\phi\) = observer’s latitude
\(\omega\) = hour angle

References:

  • Mousavi Maleki, Hizam, and Gomes [Maleki.2017] p. 5

  • Duffie [Duffie.2013] p. 15 (formula 1.6.5) with \(\sin (\alpha_s) = \cos (\theta_z)\)

Zenith Angle

Represents the angle between the vertical and the line to the sun, that is, the angle of incidence of beam radiance on a horizontal surface.

\[ \theta_z = (\frac{\pi}{2}) - \alpha_s \]

with
\(\alpha_s\) = solar altitude angle

References: See Solar Altitude Angle

Incidence Angle

The angle of incidence is the angle between the Sun’s rays and the PV panel. It can be calculated as follows:

\[\begin{split} \begin{aligned} \theta_g = \arccos(\sin(\delta) \cdot \sin(\phi) \cdot \cos(\gamma_e) \\ - \sin(\delta) \cdot \cos(\phi) \cdot \sin(\gamma_e) \cdot \cos(\alpha_e) \\ + \cos(\delta) \cdot \cos(\phi) \cdot \cos(\gamma_e) \cdot \cos(\omega) \\ + \cos(\delta) \cdot \sin(\phi) \cdot \sin(\gamma_e) \cdot \cos(\alpha_e) \cdot \cos(\omega) \\ + \cos(\delta) \cdot \sin(\gamma_e) \cdot \sin(\alpha_e) \cdot \sin(\omega)) \end{aligned} \end{split}\]

with
\(\alpha_e\) = surface azimuth angle
\(\gamma_e\) = slope angle of the surface
\(\delta\) = the declination angle
\(\phi\) = observer’s latitude
\(\omega\) = hour angle

References:

Air Mass

Calculating the air mass ratio by dividing the radius of the earth with approx. effective height of the atmosphere (each in kilometer):

\[ \mathrm{airmassratio} = (\frac{6371 km}{9 km}) = 707.8\overline{8} \]
\[ \mathrm{airmass} = \sqrt{(707.8\overline{8} \cdot \cos({\theta_z})^2 +2 \cdot 707.8\overline{8} +1)} - 707.8\overline{8} \cdot \cos(\theta_z)) \]

References:

Extraterrestrial Radiance

The extraterrestrial radiance \(G_0\) is calculated by multiplying the eccentricity correction factor.

\[\begin{split} \begin{aligned} e = 1.00011 + 0.034221 \cdot \cos(J) + 0.001280 \cdot \sin(J) \\ + 0.000719 \cdot \cos(2 \cdot J) + 0.000077 \cdot \sin(2 \cdot J) \end{aligned} \end{split}\]

with the solar constant

\[ G_{\mathrm{SC}} = 1367 {\frac{W}{m^2}} \]

with
\(J\) = day angle

References:

Beam Irradiance on Sloped Surface

For our use case, \(\omega_2\) is normally set to the hour angle one hour after \(\omega_1\). Within one hour distance to sunrise/sunset, we adjust \(\omega_1\) and \(\omega_2\) accordingly:

\[\begin{split} \begin{aligned} (\omega_1, \omega_2) = \begin{cases} (\omega_{\mathrm{SR}}, \omega_{\mathrm{SR}} + \Delta\omega), & \text{for} (\omega_{\mathrm{SR}}-\frac{\Delta \omega}{2}) < \omega < \omega_{\mathrm{SR}} \\ (\omega, \omega+ \Delta\omega), & \text{for } \omega_{\mathrm{SR}} \le \omega \le (\omega_{\mathrm{SS}}- \Delta\omega) \\ (\omega_{\mathrm{SS}}-\Delta\omega,\omega_{\mathrm{SS}}), & \text{for }(\omega_{\mathrm{SR}}-\Delta\omega) < \omega < (\omega_{\mathrm{SS}}-\frac{\Delta\omega}{2}) \end{cases} \end{aligned} \end{split}\]

Additionally, the condition \(\theta_g < 90°\) must be met (the sun must not be behind the surface).

with
\(\omega\) = hour angle
\(\omega_{\mathrm{SS}}\) = hour angle \(\omega\) at sunset
\(\omega_{\mathrm{SR}}\) = hour angle \(\omega\) at sunrise
\(\Delta\omega\) = \(15^\circ \cdot (\frac {\pi}{180^\circ})\) (one hour worth of \(\omega\))

From here on, formulas from given reference below are used:

\[\begin{split} \begin{aligned} a = (\sin(\delta) \cdot \sin(\phi) \cdot \cos(\gamma_e) - \sin(\delta) \cdot \cos(\phi) \cdot \sin(\gamma_e) \cdot \cos(\alpha_e)) \cdot (\omega_2 - \omega_1) \\ + (\cos(\delta) \cdot \cos(\phi) \cdot \cos(\gamma_e) + \cos(\delta) \cdot \sin(\phi) \cdot \sin(\gamma_e) \cdot \cos(\alpha_e)) \cdot (\sin(\omega_2) \\ - \sin(\omega_1)) - (\cos(\delta) \cdot \sin(\gamma_e) \cdot \sin(\alpha_e)) \cdot (\cos(\omega_2) - \cos(\omega_1)) \end{aligned} \end{split}\]
\[ b = (\cos(\phi) \cdot \cos(\delta)) \cdot (\sin(\omega_2) - \sin(\omega_1)) + (\sin(\phi) \cdot \sin(\delta)) \cdot (\omega_2 - \omega_1) \]
\[ G_{\mathrm{beam,S}} = G_{\mathrm{beam,H}} \cdot \frac{a}{b} \]

with
\(\delta\) = the declination angle
\(\phi\) = observer’s latitude
\(\gamma_e\) = slope angle of the surface
\(\omega_1\) = hour angle \(\omega\)
\(\omega_2\) = hour angle \(\omega\) + 1 hour
\(\alpha_e\) = surface azimuth angle
\(G_{\mathrm{beam,H}}\) = beam irradiance (horizontal surface)

Please note:

  1. \(\frac{1}{180}\pi\) is omitted from these formulas, as we are already working with data in radians.

  2. In contrast to the primary source [Duffie.2013], radiance values instead of radiation values are used here, as described above.

Reference:

Diffuse Irradiance on Sloped Surface

The diffuse irradiance is computed using the Perez model, which divides the irradiance into three parts. First, there is an intensified irradiance from the direct vicinity of the sun. Furthermore, there is Rayleigh scattering, backscatter (which lead to increased in intensity on the horizon) and isotropic irradiance considered.

A cloud index is defined by

\[ \epsilon = \frac{\frac{G_{\mathrm{dif,H}} + G_{\mathrm{beam,N}}}{G_{\mathrm{dif,H}}} + 5.535 \cdot 10^{-6} \cdot (\theta_z)^3}{1 + 5.535 \cdot 10^{-6} \cdot (\theta_z)^3} \]

with angle \(\theta_z\) values in degrees ([Duffie.2013] p. 94) and \(G_{\mathrm{beam,N}} = \frac{G_{\mathrm{beam,H}}}{\cos (\theta_z)}\) ([Duffie.2013] p. 95).

Calculating a brightness index

\[ \Delta = m \cdot \frac{G_{\mathrm{dif,H}}}{G_0} \]

Perez Fij coefficients (Myers 2017):

\(\epsilon\) is sorted into one of eight bins according to its value:

\(\epsilon\) low

\(\epsilon\) high

Bin number\(x\)

1

1.065

1

1.065

1.230

2

1.230

1.500

3

1.500

1.950

4

1.950

2.800

5

2.800

4.500

6

4.500

6.200

7

6.200

\(\infty\)

8

In order to calculate indexes representing the horizon brightness and the brightness in the vicinity of the sun, the following factors are computed.

\[ F_{11}(x) = -0.0161 \cdot x^3 + 0.1840 \cdot x^2 - 0.3806 \cdot x + 0.2324 \]
\[ F_{12}(x) = 0.0134 \cdot x^4 - 0.1938 \cdot x^3 + 0.8410 \cdot x^2 - 1.4018 \cdot x + 1.3579 \]
\[ F_{13}(x) = 0.0032 \cdot x^3 - 0.028 \cdot x^2 - 0.0056 \cdot x - 0.0385 \]
\[ F_{21}(x) = -0.0048 \cdot x^3 + 0.0536 \cdot x^2 - 0.1049 \cdot x + 0.0034 \]
\[ F_{22}(x) = 0.0012 \cdot x^3 - 0.0067 \cdot x^2 + 0.0091 \cdot x - 0.0269 \]
\[ F_{23}(x) = 0.0052 \cdot x^3 - 0.0971 \cdot x^2 + 0.2856 \cdot x - 0.1389 \]

Horizon brightness index:

\[ F_1 = F_{11}(x) + F_{12}(x) \cdot \Delta + F_{13}(x) \cdot \theta_z \]

Sun ambient brightness index:

\[ F_2 = F_{21}(x) + F_{22}(x) \cdot \Delta + F_{23}(x) \cdot \theta_z \]

Using the factors

\[ a = max(0, \cos(\theta_g)) \]

and

\[ b = max(0.087, \sin(\alpha_s)) \]

the diffuse irradiance can be calculated:

\[ G_{\mathrm{dif,S}} = G_{\mathrm{dif,H}} \cdot (\frac{1}{2} \cdot (1 + cos(\gamma_e)) \cdot (1- F_1) + \frac{a}{b} \cdot F_1 + F_2 \cdot \sin(\gamma_e)) \]

with
\(\theta_z\) = zenith angle
\(\theta_g\) = angle of incidence
\(\alpha_s\) = solar altitude angle
\(\gamma_e\) = slope angle of the surface
\(G_0\) = extraterrestrial radiance
\(m\) = air mass
\(G_{\mathrm{beam,H}}\) = beam irradiance (horizontal surface)
\(G_{\mathrm{beam,N}}\) = beam irradiance (normal incidence, thus irradiance on a plane normal to the direction of the beam)
\(G_{\mathrm{dif,H}}\) = diffuse irradiance (horizontal surface)

Please note: In contrast to the primary source [Duffie.2013], radiance values instead of radiation values are used here, as described above.

References:

Reflected Irradiance on Sloped Surface

\[ G_{\mathrm{ref},S} = G_{\mathrm{Ges,H}} \cdot \frac{\rho}{2} \cdot (1- \cos(\gamma_e)) \]

with
\(G_{\mathrm{Ges,H}}\) = total horizontal irradiance (\(G_{\mathrm{beam,H}} + G_{\mathrm{dif,H}})\)
\(\gamma_e\) = slope angle of the surface
\(\rho\) = albedo

Reference:

  • Mousavi Maleki, Hizam, and Gomes [Maleki.2017] p. 19

Output

Received energy is calculated as the sum of all three types of irradiance.

\[ G_{\mathrm{total}} = G_{\mathrm{beam,S}} + G_{\mathrm{dif,S}} + G_{\mathrm{ref,S}} \]

with
\(G_{\mathrm{beam,S}}\) = Beam irradiance
\(G_{\mathrm{dif,S}}\) = Diffuse irradiance
\(G_{\mathrm{ref,S}}\) = Reflected irradiance

A generator correction factor (depending on month surface slope \(\gamma_e\)) and a temperature correction factor (depending on month) multiplied on top.

It is checked whether proposed output exceeds maximum (\(p_{\mathrm{max}}\)), in which case a warning is logged. If output falls below activation threshold, it is set to 0.