Merge pull request #1213 from ie3-institute/sp/#1212-PvModel-irradiance

Basing `PvModel` on irradiance, not irradiation

CHANGELOG.md

@@ -140,6 +140,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Replaced Java Durations with Scala Durations [#1068](https://github.com/ie3-institute/simona/issues/1068)
 - Typo and format of `ThermalGrid` and `ThermalHouse` ScalaDocs [#1196](https://github.com/ie3-institute/simona/issues/1196)
 - Refactor `EmRuntimeConfig` [#1181](https://github.com/ie3-institute/simona/issues/1181)
+- Based `PvModel` calculations on irradiance (power per area) instead of irradiation (energy per area) [#1212](https://github.com/ie3-institute/simona/issues/1212)
 
 ### Fixed
 - Fix rendering of references in documentation [#505](https://github.com/ie3-institute/simona/issues/505)

docs/readthedocs/models/pv_model.md

@@ -4,7 +4,8 @@
 
 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](https://ieeexplore.ieee.org/abstract/document/7038345). Since then several adaptions has been made that are documented as follows.
+The initial parts of the model are presented in the paper _Agent based approach to model photovoltaic feed-in in distribution network planning_ by {cite:cts}`Seack.2014`.
+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.
 
@@ -20,21 +21,29 @@ Please refer to {doc}`PowerSystemDataModel - PV Model <psdm:models/input/partici
 
 ![Sequence Diagram Behaviour PV Model](http://www.plantuml.com/plantuml/proxy?cache=no&src=https://raw.githubusercontent.com/ie3-institute/simona/dev/docs/uml/main/models/pv_model/BehaviourPvModel.puml)
 
-## Output visualization
-
 ## Calculations
 
-The energy produced by a photovoltaic (pv) unit in a specific time step is based on the diffuse and direct radiation provided by the used weather data. The following steps are done to calculate (= estimate) the power feed by the pv.
+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 {cite:p}`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 ({cite:p}`Duffie.2013` p. 86, footnote).
 
-To calculate the overall feed in of the pv unit, the sum of the direct radiation, diffuse radiation and reflected radiation is needed. In the following, the formulas to calculate each of these radiations are presented and explained. The sections end with the formula to calculate the corresponding power feed in.
+_Irradiance_ and _irradiation_ describe power and energy arriving at a surface, while the terms _radiance_ and _radiation_ are used for general purposes ({cite:p}`Iqbal.1983` Ch. 2.6).
 
-**Caution:** all angles are given in radian!
+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)](https://www.photovoltaik.org/wissen/azimutwinkel)
+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](https://www.photovoltaik.org/wissen/azimutwinkel)).
 
 ### 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:
+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})
@@ -59,7 +68,9 @@ $$
 
 ### Hour Angle
 
-The hour angle is a conceptual description of the rotation of the earth around its polar axis. It starts with a negative value in the morning, arrives at 0° at noon (solar time) and ends with a positive value in the evening. The hour angle (in radian!) is calculated as follows
+The hour angle is a conceptual description of the rotation of the earth around its polar axis.
+It starts with a negative value in the morning, arrives at 0° at noon (solar time) and ends with a positive value in the evening.
+The hour angle (in radian!) is calculated as follows
 
 $$
 \omega = ((12 - ST) \cdot 15) \cdot (\frac{\pi}{180})
@@ -99,7 +110,8 @@ $$
 *with*\
 **J** = day angle (in radian!)
 
-**Note:** The used formulas are based on *\"DIN 5034-2: Tageslicht in Innenräumen, Grundlagen.\"* and therefore valid especially for Germany and Europe. For international calculations a more general formulation that can be found in [Maleki, S.A., Hizam, H., & Gomes, C. (2017). Estimation of Hourly, Daily and Monthly Global Solar Radiation on Inclined Surfaces: Models Re-Visited.](https://res.mdpi.com/d_attachment/energies/energies-10-00134/article_deploy/energies-10-00134-v2.pdf) might be used.
+**Note:** The used formulas are based on *\"DIN 5034-2: Tageslicht in Innenräumen, Grundlagen.\"* and therefore valid especially for Germany and Europe.
+For international calculations a more general formulation that can be found in {cite:p}`Maleki.2017` might be used.
 
 **References:**
 
@@ -110,7 +122,9 @@ $$
 
 ### 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_{SS}$) is positive and the sunrise angle ($\omega_{SR}$) is negative. Both can be calculated from:
+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_{SS}$) is positive and the sunrise angle ($\omega_{SR}$) is negative.
+Both can be calculated from:
 
 $$
 \omega_{SS}=\cos^{-1}(-\tan (\phi) \cdot \tan (\delta))
@@ -151,7 +165,7 @@ $$
 
 ### Zenith Angle
 
-Represents the angle between the vertical and the line to the sun, that is, the angle of incidence of beam radiation on a horizontal surface.
+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}
@@ -165,7 +179,8 @@ 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:
+The angle of incidence is the angle between the Sun\'s rays and the PV panel.
+It can be calculated as follows:
 
 $$
 \begin{eqnarray*}
@@ -195,7 +210,7 @@ $$
 
 ### Air Mass
 
-Calculating the air mass ratio by dividing the radius of the earth with approx. effective height of the atmosphere (each in kilometer)
+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}
@@ -211,9 +226,9 @@ $$
 * {cite:cts}`WikiAirMass`
 
 
-### Extraterrestrial Radiation
+### Extraterrestrial Radiance
 
-The extraterrestrial radiation $I_0$ is calculated by multiplying the eccentricity correction factor
+The extraterrestrial radiance $G_0$ is calculated by multiplying the eccentricity correction factor.
 
 $$
 \begin{eqnarray*}
@@ -238,9 +253,10 @@ $$
 * {cite:cts}`Duffie.2013` (the fifth ed. seems to have a typo in formula (1.4.1b): factor $0.000719$ is missing a zero)
 
 
-### Beam Radiation on Sloped Surface
+### 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:
+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{eqnarray*}
@@ -276,41 +292,45 @@ b = (\cos(\phi) \cdot \cos(\delta)) \cdot (\sin(\omega_{2}) - \sin(\omega_{1}))
 $$
 
 $$
-E_{\mathrm{beam},S} = E_{\mathrm{beam},H} \cdot \frac{a}{b}
+G_{\mathrm{beam},S} = G_{\mathrm{beam},H} \cdot \frac{a}{b}
 $$
 
-**Please note:** $\frac{1}{180}\pi$ is omitted from these formulas, as we are already working with data in *radians*.
-
 *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\
-**$E_{\mathrm{beam},H}$** = beam radiation (horizontal surface)
+**$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 {cite:p}`Duffie.2013`, radiance values instead of radiation values are used here, as described above.
 
 **Reference:**
 
 * {cite:cts}`Duffie.2013` p. 88
 
 
-### Diffuse Radiation on Sloped Surface
+### Diffuse Irradiance on Sloped Surface
 
-The diffuse radiation is computed using the Perez model, which divides the radiation in three parts. First, there is an intensified radiation from the direct vicinity of the sun. Furthermore, there is Rayleigh scattering, backscatter (which lead to increased in intensity on the horizon) and isotropic radiation considered.
+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{E_{\mathrm{dif},H} + E_{\mathrm{beam},N}}{E_{\mathrm{dif},H}} + 5.535 \cdot 10^{-6} \cdot \theta_{z}^3}{1 + 5.535 \cdot 10^{-6} \cdot \theta_{z}^3}
+\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** ({cite:cts}`Duffie.2013` p. 94) and $E_{\mathrm{beam},N} = \frac{E_{\mathrm{beam},H}}{\cos (\theta_z)}$ ({cite:cts}`Duffie.2013` p. 95).
+with angle $\theta_z$ values in **degrees** ({cite:p}`Duffie.2013` p. 94) and $G_{\mathrm{beam},N} = \frac{G_{\mathrm{beam},H}}{\cos (\theta_z)}$ ({cite:p}`Duffie.2013` p. 95).
 
 Calculating a brightness index
 
 $$
-\Delta = m \cdot \frac{E_{\mathrm{dif},H}}{I_{0}}
+\Delta = m \cdot \frac{G_{\mathrm{dif},H}}{G_{0}}
 $$
 
 **Perez Fij coefficients (Myers 2017):**
@@ -385,10 +405,10 @@ $$
 b = max(0.087, \sin(\alpha_{s}))
 $$
 
-the diffuse radiation can be calculated:
+the diffuse irradiance can be calculated:
 
 $$
-E_{\mathrm{dif},S} = E_{\mathrm{dif},H} \cdot (\frac{1}{2} \cdot (1 +
+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}))
 $$
@@ -398,11 +418,13 @@ $$
 **$\theta_{g}$** = angle of incidence\
 **$\alpha_{s}$** = solar altitude angle\
 **$\gamma_{e}$** = slope angle of the surface\
-**$I_{0}$** = Extraterrestrial Radiation\
+**$G_{0}$** = extraterrestrial radiance\
 **$m$** = air mass\
-**$E_{\mathrm{beam},H}$** = beam radiation (horizontal surface)\
-**$E_{\mathrm{beam},N}$** = beam radiation (normal incidence, thus radiation on a plane normal to the direction of the beam)\
-**$E_{\mathrm{dif},H}$** = diffuse radiation (horizontal surface)
+**$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 {cite:p}`Duffie.2013`, radiance values instead of radiation values are used here, as described above.
 
 **References:**
 
@@ -412,15 +434,15 @@ $$
 * {cite:cts}`Duffie.2013` p. 95f
 
 
-### Reflected Radiation on Sloped Surface
+### Reflected Irradiance on Sloped Surface
 
 $$
-E_{\mathrm{ref},S} = E_{\mathrm{Ges},H} \cdot \frac{\rho}{2} \cdot (1-
+G_{\mathrm{ref},S} = G_{\mathrm{Ges},H} \cdot \frac{\rho}{2} \cdot (1-
 \cos(\gamma_{e}))
 $$
 
 *with*\
-**$E_{\mathrm{Ges},H}$** = total horizontal radiation ($E_{\mathrm{beam},H} + E_{\mathrm{dif},H})$\
+**$G_{\mathrm{Ges},H}$** = total horizontal irradiance ($G_{\mathrm{beam},H} + G_{\mathrm{dif},H})$\
 **$\gamma_e$** = slope angle of the surface\
 **$\rho$** = albedo
 
@@ -431,17 +453,18 @@ $$
 
 ### Output
 
-Received energy is calculated as the sum of all three types of irradiation.
+Received energy is calculated as the sum of all three types of irradiance.
 
 $$
-E_{\mathrm{total}} = E_{\mathrm{beam},S} + E_{\mathrm{dif},S} + E_{\mathrm{ref},S}
+G_{\mathrm{total}} = G_{\mathrm{beam},S} + G_{\mathrm{dif},S} + G_{\mathrm{ref},S}
 $$
 
 *with*\
-**$E_{\mathrm{beam},S}$** = Beam radiation\
-**$E_{\mathrm{dif},S}$** = Diffuse radiation\
-**$E_{\mathrm{ref},S}$** = Reflected radiation
+**$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_{max}$), in which case a warning is logged. If output falls below activation threshold, it is set to 0.
+It is checked whether proposed output exceeds maximum ($p_{max}$), in which case a warning is logged.
+If output falls below activation threshold, it is set to 0.

src/main/scala/edu/ie3/simona/agent/participant/pv/PvAgentFundamentals.scala

@@ -50,7 +50,6 @@ import edu.ie3.simona.ontology.messages.flex.FlexibilityMessage.{
   FlexResponse,
 }
 import edu.ie3.simona.ontology.messages.services.WeatherMessage.WeatherData
-import edu.ie3.simona.service.weather.WeatherService.FALLBACK_WEATHER_STEM_DISTANCE
 import edu.ie3.simona.util.TickUtil.TickLong
 import edu.ie3.util.quantities.PowerSystemUnits.PU
 import edu.ie3.util.quantities.QuantityUtils.RichQuantityDouble
@@ -203,17 +202,6 @@ protected trait PvAgentFundamentals
       baseStateData.startDate
     val dateTime = tick.toDateTime
 
-    val tickInterval =
-      baseStateData.receivedSecondaryDataStore
-        .lastKnownTick(tick - 1) match {
-        case Some(dataTick) =>
-          tick - dataTick
-        case _ =>
-          /* At the first tick, we are not able to determine the tick interval from last tick
-           * (since there is none). Then we use a fallback pv stem distance. */
-          FALLBACK_WEATHER_STEM_DISTANCE
-      }
-
     // take the last weather data, not necessarily the one for the current tick:
     // we might receive flex control messages for irregular ticks
     val (_, secondaryData) = baseStateData.receivedSecondaryDataStore
@@ -240,7 +228,6 @@ protected trait PvAgentFundamentals
 
     PvRelevantData(
       dateTime,
-      tickInterval,
       weatherData.diffIrr,
       weatherData.dirIrr,
     )