Add Prilliman et al transience model to pvlib.temperature

.github/PULL_REQUEST_TEMPLATE.md

@@ -3,7 +3,7 @@
  - [ ] Closes #xxxx
  - [ ] I am familiar with the [contributing guidelines](https://pvlib-python.readthedocs.io/en/latest/contributing.html)
  - [ ] Tests added
- - [ ] Updates entries to [`docs/sphinx/source/api.rst`](https://github.com/pvlib/pvlib-python/blob/master/docs/sphinx/source/api.rst) for API changes.
+ - [ ] Updates entries in [`docs/sphinx/source/reference`](https://github.com/pvlib/pvlib-python/blob/master/docs/sphinx/source/reference) for API changes.
  - [ ] Adds description and name entries in the appropriate "what's new" file in [`docs/sphinx/source/whatsnew`](https://github.com/pvlib/pvlib-python/tree/master/docs/sphinx/source/whatsnew) for all changes. Includes link to the GitHub Issue with `` :issue:`num` `` or this Pull Request with `` :pull:`num` ``. Includes contributor name and/or GitHub username (link with `` :ghuser:`user` ``).
  - [ ] New code is fully documented. Includes [numpydoc](https://numpydoc.readthedocs.io/en/latest/format.html) compliant docstrings, examples, and comments where necessary.
  - [ ] Pull request is nearly complete and ready for detailed review.

docs/sphinx/source/reference/pv_modeling.rst

@@ -43,6 +43,7 @@ PV temperature models
    temperature.fuentes
    temperature.ross
    temperature.noct_sam
+   temperature.prilliman
    pvsystem.PVSystem.get_cell_temperature
 
 Temperature Model Parameters

docs/sphinx/source/whatsnew/v0.9.1.rst

@@ -11,6 +11,8 @@ Deprecations
 
 Enhancements
 ~~~~~~~~~~~~
+* Added :py:func:`pvlib.temperature.prilliman` for modeling cell temperature
+  at short time steps (:issue:`1081`, :pull:`1391`)
 
 Bug fixes
 ~~~~~~~~~

pvlib/clearsky.py

@@ -679,23 +679,6 @@ def _to_centered_series(vals, idx, samples_per_window):
     return pd.Series(index=idx, data=vals).shift(shift)
 
 
-def _get_sample_intervals(times, win_length):
-    """ Calculates time interval and samples per window for Reno-style clear
-    sky detection functions
-    """
-    deltas = np.diff(times.values) / np.timedelta64(1, '60s')
-
-    # determine if we can proceed
-    if times.inferred_freq and len(np.unique(deltas)) == 1:
-        sample_interval = times[1] - times[0]
-        sample_interval = sample_interval.seconds / 60  # in minutes
-        samples_per_window = int(win_length / sample_interval)
-        return sample_interval, samples_per_window
-    else:
-        raise NotImplementedError('algorithm does not yet support unequal '
-                                  'times. consider resampling your data.')
-
-
 def _clear_sample_index(clear_windows, samples_per_window, align, H):
     """
     Returns indices of clear samples in clear windows
@@ -849,8 +832,8 @@ def detect_clearsky(measured, clearsky, times=None, window_length=10,
     else:
         clear = clearsky
 
-    sample_interval, samples_per_window = _get_sample_intervals(times,
-                                                                window_length)
+    sample_interval, samples_per_window = \
+        tools._get_sample_intervals(times, window_length)
 
     # generate matrix of integers for creating windows with indexing
     H = hankel(np.arange(samples_per_window),

pvlib/temperature.py

@@ -7,6 +7,10 @@
 import pandas as pd
 from pvlib.tools import sind
 from pvlib._deprecation import warn_deprecated
+from pvlib.tools import _get_sample_intervals
+import scipy
+import warnings
+
 
 TEMPERATURE_MODEL_PARAMETERS = {
     'sapm': {
@@ -821,3 +825,155 @@ def noct_sam(poa_global, temp_air, wind_speed, noct, module_efficiency,
     heat_loss = 1 - module_efficiency / tau_alpha
     wind_loss = 9.5 / (5.7 + 3.8 * wind_adj)
     return temp_air + cell_temp_init * heat_loss * wind_loss
+
+
+def prilliman(temp_cell, wind_speed, unit_mass=11.1, coefficients=None):
+    """
+    Smooth short-term cell temperature transients using the Prilliman model.
+
+    The Prilliman et al. model [1]_ applies a weighted moving average to
+    the output of a steady-state cell temperature model to account for
+    a module's thermal inertia by smoothing the cell temperature's
+    response to changing weather conditions.
+
+    .. warning::
+        This implementation requires the time series inputs to be regularly
+        sampled in time with frequency less than 20 minutes.  Data with
+        irregular time steps should be resampled prior to using this function.
+
+    Parameters
+    ----------
+    temp_cell : pandas.Series with DatetimeIndex
+        Cell temperature modeled with steady-state assumptions. [C]
+
+    wind_speed : pandas.Series
+        Wind speed, adjusted to correspond to array height [m/s]
+
+    unit_mass : float, default 11.1
+        Total mass of module divided by its one-sided surface area [kg/m^2]
+
+    coefficients : 4-element list-like, optional
+        Values for coefficients a_0 through a_3, see Eq. 9 of [1]_
+
+    Returns
+    -------
+    temp_cell : pandas.Series
+        Smoothed version of the input cell temperature. Input temperature
+        with sampling interval >= 20 minutes is returned unchanged. [C]
+
+    Notes
+    -----
+    This smoothing model was developed and validated using the SAPM
+    cell temperature model for the steady-state input.
+
+    Smoothing is done using the 20 minute window behind each temperature
+    value. At the beginning of the series where a full 20 minute window is not
+    possible, partial windows are used instead.
+
+    Output ``temp_cell[k]`` is NaN when input ``wind_speed[k]`` is NaN, or
+    when no non-NaN data are in the input temperature for the 20 minute window
+    preceding index ``k``.
+
+    References
+    ----------
+    .. [1] M. Prilliman, J. S. Stein, D. Riley and G. Tamizhmani,
+       "Transient Weighted Moving-Average Model of Photovoltaic Module
+       Back-Surface Temperature," IEEE Journal of Photovoltaics, 2020.
+       :doi:`10.1109/JPHOTOV.2020.2992351`
+    """
+
+    # `sample_interval` in minutes:
+    sample_interval, samples_per_window = \
+        _get_sample_intervals(times=temp_cell.index, win_length=20)
+
+    if sample_interval >= 20:
+        warnings.warn("temperature.prilliman only applies smoothing when "
+                      "the sampling interval is shorter than 20 minutes "
+                      f"(input sampling interval: {sample_interval} minutes);"
+                      " returning input temperature series unchanged")
+        # too coarsely sampled for smoothing to be relevant
+        return temp_cell
+
+    # handle cases where the time series is shorter than 20 minutes total
+    samples_per_window = min(samples_per_window, len(temp_cell))
+
+    # prefix with NaNs so that the rolling window is "full",
+    # even for the first actual value:
+    prefix = np.full(samples_per_window, np.nan)
+    temp_cell_prefixed = np.append(prefix, temp_cell.values)
+
+    # generate matrix of integers for creating windows with indexing
+    H = scipy.linalg.hankel(np.arange(samples_per_window),
+                            np.arange(samples_per_window - 1,
+                                      len(temp_cell_prefixed) - 1))
+    # each row of `subsets` is the values in one window
+    subsets = temp_cell_prefixed[H].T
+
+    # `subsets` now looks like this (for 5-minute data, so 4 samples/window)
+    # where "1." is a stand-in for the actual temperature values
+    # [[nan, nan, nan, nan],
+    #  [nan, nan, nan,  1.],
+    #  [nan, nan,  1.,  1.],
+    #  [nan,  1.,  1.,  1.],
+    #  [ 1.,  1.,  1.,  1.],
+    #  [ 1.,  1.,  1.,  1.],
+    #  [ 1.,  1.,  1.,  1.],
+    #  ...
+
+    # calculate weights for the values in each window
+    if coefficients is not None:
+        a = coefficients
+    else:
+        # values from [1], Table II
+        a = [0.0046, 0.00046, -0.00023, -1.6e-5]
+
+    wind_speed = wind_speed.values
+    p = a[0] + a[1]*wind_speed + a[2]*unit_mass + a[3]*wind_speed*unit_mass
+    # calculate the time lag for each sample in the window, paying attention
+    # to units (seconds for `timedeltas`, minutes for `sample_interval`)
+    timedeltas = np.arange(samples_per_window, 0, -1) * sample_interval * 60
+    weights = np.exp(-p[:, np.newaxis] * timedeltas)
+
+    # Set weights corresponding to the prefix values to zero; otherwise the
+    # denominator of the weighted average below would be wrong.
+    # Weights corresponding to (non-prefix) NaN values must be zero too
+    # for the same reason.
+
+    # Right now `weights` is something like this
+    # (using 5-minute inputs, so 4 samples per window -> 4 values per row):
+    # [[0.0611, 0.1229, 0.2472, 0.4972],
+    #  [0.0611, 0.1229, 0.2472, 0.4972],
+    #  [0.0611, 0.1229, 0.2472, 0.4972],
+    #  [0.0611, 0.1229, 0.2472, 0.4972],
+    #  [0.0611, 0.1229, 0.2472, 0.4972],
+    #  [0.0611, 0.1229, 0.2472, 0.4972],
+    #  [0.0611, 0.1229, 0.2472, 0.4972],
+    #  ...
+
+    # After the next line, the NaNs in `subsets` will be zeros in `weights`,
+    # like this (with more zeros for any NaNs in the input temperature):
+
+    # [[0.    , 0.    , 0.    , 0.    ],
+    #  [0.    , 0.    , 0.    , 0.4972],
+    #  [0.    , 0.    , 0.2472, 0.4972],
+    #  [0.    , 0.1229, 0.2472, 0.4972],
+    #  [0.0611, 0.1229, 0.2472, 0.4972],
+    #  [0.0611, 0.1229, 0.2472, 0.4972],
+    #  [0.0611, 0.1229, 0.2472, 0.4972],
+    #  ...
+
+    weights[np.isnan(subsets)] = 0
+
+    # change the first row of weights from zero to nan -- this is a
+    # trick to prevent div by zero warning when dividing by summed weights
+    weights[0, :] = np.nan
+
+    # finally, take the weighted average of each window:
+    # use np.nansum for numerator to ignore nans in input temperature, but
+    # np.sum for denominator to propagate nans in input wind speed.
+    numerator = np.nansum(subsets * weights, axis=1)
+    denominator = np.sum(weights, axis=1)
+    smoothed = numerator / denominator
+    smoothed[0] = temp_cell.values[0]
+    smoothed = pd.Series(smoothed, index=temp_cell.index)
+    return smoothed

pvlib/tests/test_temperature.py

@@ -8,6 +8,8 @@
 from pvlib import temperature, tools
 from pvlib._deprecation import pvlibDeprecationWarning
 
+import re
+
 
 @pytest.fixture
 def sapm_default():
@@ -293,3 +295,71 @@ def test_noct_sam_options():
 def test_noct_sam_errors():
     with pytest.raises(ValueError):
         temperature.noct_sam(1000., 25., 1., 34., 0.2, array_height=3)
+
+
+def test_prilliman():
+    # test against values calculated using pvl_MAmodel_2, see pvlib #1081
+    times = pd.date_range('2019-01-01', freq='5min', periods=8)
+    cell_temperature = pd.Series([0, 1, 3, 6, 10, 15, 21, 27], index=times)
+    wind_speed = pd.Series([0, 1, 2, 3, 2, 1, 2, 3])
+
+    # default coeffs
+    expected = pd.Series([0, 0, 0.7047457, 2.21176412, 4.45584299, 7.63635512,
+                          12.26808265, 18.00305776], index=times)
+    actual = temperature.prilliman(cell_temperature, wind_speed, unit_mass=10)
+    assert_series_equal(expected, actual)
+
+    # custom coeffs
+    coefficients = [0.0046, 4.5537e-4, -2.2586e-4, -1.5661e-5]
+    expected = pd.Series([0, 0, 0.70716941, 2.2199537, 4.47537694, 7.6676931,
+                          12.30423167, 18.04215198], index=times)
+    actual = temperature.prilliman(cell_temperature, wind_speed, unit_mass=10,
+                                   coefficients=coefficients)
+    assert_series_equal(expected, actual)
+
+    # even very short inputs < 20 minutes total still work
+    times = pd.date_range('2019-01-01', freq='1min', periods=8)
+    cell_temperature = pd.Series([0, 1, 3, 6, 10, 15, 21, 27], index=times)
+    wind_speed = pd.Series([0, 1, 2, 3, 2, 1, 2, 3])
+    expected = pd.Series([0, 0, 0.53557976, 1.49270094, 2.85940173,
+                          4.63914366, 7.09641845, 10.24899272], index=times)
+    actual = temperature.prilliman(cell_temperature, wind_speed, unit_mass=12)
+    assert_series_equal(expected, actual)
+
+
+def test_prilliman_coarse():
+    # if the input series time step is >= 20 min, input is returned unchanged,
+    # and a warning is emitted
+    times = pd.date_range('2019-01-01', freq='30min', periods=3)
+    cell_temperature = pd.Series([0, 1, 3], index=times)
+    wind_speed = pd.Series([0, 1, 2])
+    msg = re.escape("temperature.prilliman only applies smoothing when the "
+                    "sampling interval is shorter than 20 minutes (input "
+                    "sampling interval: 30.0 minutes); returning "
+                    "input temperature series unchanged")
+    with pytest.warns(UserWarning, match=msg):
+        actual = temperature.prilliman(cell_temperature, wind_speed)
+    assert_series_equal(cell_temperature, actual)
+
+
+def test_prilliman_nans():
+    # nans in inputs are handled appropriately; nans in input tcell
+    # are ignored but nans in wind speed cause nan in output
+    times = pd.date_range('2019-01-01', freq='1min', periods=8)
+    cell_temperature = pd.Series([0, 1, 3, 6, 10, np.nan, 21, 27], index=times)
+    wind_speed = pd.Series([0, 1, 2, 3, 2, 1, np.nan, 3])
+    actual = temperature.prilliman(cell_temperature, wind_speed)
+    expected = pd.Series([True, True, True, True, True, True, False, True],
+                         index=times)
+    assert_series_equal(actual.notnull(), expected)
+
+    # check that nan temperatures do not mess up the weighted average;
+    # the original implementation did not set weight=0 for nan values,
+    # so the numerator of the weighted average ignored nans but the
+    # denominator (total weight) still included the weight for the nan.
+    cell_temperature = pd.Series([1, 1, 1, 1, 1, np.nan, 1, 1], index=times)
+    wind_speed = pd.Series(1, index=times)
+    actual = temperature.prilliman(cell_temperature, wind_speed)
+    # original implementation would return some values < 1 here
+    expected = pd.Series(1., index=times)
+    assert_series_equal(actual, expected)

pvlib/tools.py

@@ -344,3 +344,20 @@ def _golden_sect_DataFrame(params, VL, VH, func):
             raise Exception("EXCEPTION:iterations exceeded maximum (50)")
 
     return func(df, 'V1'), df['V1']
+
+
+def _get_sample_intervals(times, win_length):
+    """ Calculates time interval and samples per window for Reno-style clear
+    sky detection functions
+    """
+    deltas = np.diff(times.values) / np.timedelta64(1, '60s')
+
+    # determine if we can proceed
+    if times.inferred_freq and len(np.unique(deltas)) == 1:
+        sample_interval = times[1] - times[0]
+        sample_interval = sample_interval.seconds / 60  # in minutes
+        samples_per_window = int(win_length / sample_interval)
+        return sample_interval, samples_per_window
+    else:
+        raise NotImplementedError('algorithm does not yet support unequal '
+                                  'times. consider resampling your data.')