Implement dynamic faiman model and function to fit this to measurements

docs/examples/operating-temperature/README.rst

@@ -0,0 +1,3 @@
+Operating Temperature
+---------------------
+

docs/examples/operating-temperature/plot_dynamic_faiman.py

@@ -0,0 +1,154 @@
+"""
+A simple way to incorporate thermal inertia
+===========================================
+
+What can be simpler than a moving average?
+
+"""
+
+# %%
+#
+# Applying a moving average filter to the PV array operating conditions
+# is a simple technique to compensate for the thermal inertia of the module,
+# which delays and dampens temperature fluctuations.
+# It is useful for simulating at small time steps, but even more useful for
+# fitting models to field data as demonstrated in [1]_.
+# The functions :py:func:`pvlib.temperature.faiman_dyn` and
+# :py:func:`pvlib.temperature.fit_faiman_dyn` incorporate this moving average
+# technique.
+#
+# This example reads a csv file containing one-minute average monitoring data.
+# The function :py:func:`pvlib.temperature.fit_faiman_dyn` is used to determine
+# the model parameters and the function :py:func:`pvlib.temperature.faiman_dyn`
+# is used to demonstrate how well it worked.
+#
+# Contributed by Anton Driesse, PV Performance Labs, October 2023.
+#
+# References
+# ----------
+# .. [1] Driesse, A. (2022) "Module operating temperature model parameter
+#    determination." Sandia National Laboratories, Albuquerque NM.
+#    :doi:`10.5281/zenodo.10003736`
+#
+
+import os
+import pandas as pd
+import matplotlib.pyplot as plt
+
+import pvlib
+from pvlib.temperature import faiman, faiman_dyn, fit_faiman_dyn
+from pvlib.temperature import GenericLinearModel
+
+# %%
+#
+# Read a CSV file containing one week of weather data and module temperature
+#
+
+PVLIB_DIR = pvlib.__path__[0]
+DATA_FILE = os.path.join(PVLIB_DIR, 'data', 'tmod_sample_data_subset.csv')
+
+df = pd.read_csv(DATA_FILE, index_col=0, parse_dates=True)
+
+print(df.head())
+
+# %%
+#
+# Estimate the dynamic Faiman model parameters using data where Gpoa > 10
+# to eliminate night time data.
+#
+# The fitting procedure takes a simple approach to finding the optimal
+# thermal_inertia: it tries a range of values and chooses the one that
+# produces the lowest RMSE.
+#
+
+dff = df.copy()
+dff = dff.where(df.poa_global > 10)
+
+params, details = fit_faiman_dyn(dff.temp_pv, dff.poa_global,
+                                 dff.temp_air, dff.wind_speed,
+                                 thermal_inertia=(0, 15, 0.5),
+                                 full_output=True)
+
+for k, v in params.items():
+    print('%-15s = %5.2f' % (k, v))
+
+# %%
+#
+# With the full_output option you can obtain the results for all the values
+# of thermal_inertia that were evaluated.  The optimal point is clearly visible
+# below, but a minute shorter or longer actually doesn't make much difference
+# in the RMSE.
+#
+
+plt.figure()
+plt.plot(details.thermal_inertia, details.rmse, '.-')
+plt.grid(alpha=0.5)
+plt.xlabel('Thermal inertia [minutes]')
+plt.ylabel('RMSE [°C]')
+plt.title('Optimal values: u0=%.2f, u1=%.2f, thermal_inertia=%.1f'
+          % tuple(params.values()))
+plt.show()
+
+# %%
+#
+# Now calculate the modeled operating temperature of the PV modules.  The
+# u0 and u1 values found for the dynamic model can be used with the
+# regular Faiman model too.
+#
+
+df['temp_pv_faiman'] = faiman(df.poa_global, df.temp_air, df.wind_speed,
+                              u0=params['u0'], u1=params['u1'])
+df['temp_pv_faiman_dyn'] = faiman_dyn(df.poa_global, df.temp_air,
+                                      df.wind_speed, **params)
+
+DAY = slice('2020-03-20 7:00', '2020-03-20 19:00')
+# sphinx_gallery_thumbnail_number = 2
+plt.figure()
+plt.plot(df.temp_pv[DAY])
+plt.plot(df.temp_pv_faiman[DAY], alpha=0.5, zorder=0)
+plt.plot(df.temp_pv_faiman_dyn[DAY])
+plt.legend(['measured', 'faiman', 'faiman_dyn'])
+plt.grid(alpha=0.5)
+plt.xlabel('2020-03-20')
+plt.ylabel('PV temperature [°C]')
+plt.show()
+
+# %%
+
+plt.figure()
+l1 = plt.plot(df['temp_pv'], df['temp_pv_faiman'], '.', color='C1')
+l2 = plt.plot(df['temp_pv'], df['temp_pv_faiman_dyn'], '.', color='C2')
+plt.legend(['faiman', 'faiman_dyn'])
+l1[0].set_alpha(0.5)
+l2[0].set_alpha(0.25)
+plt.grid(alpha=0.5)
+plt.xlabel('Measured temperature [°C]')
+plt.ylabel('Modeled temperature [°C]')
+plt.show()
+
+# %%
+#
+# Both graphs above demonstrate that substantial improvement in modeled
+# operating temperature is obtained by this simple technique.  Perhaps more
+# important than this, however, is the fact that parameter values can be
+# extracted from field data with minimal or no filtering.
+#
+# Finally, translate the Faiman model parameters to other model parameters
+# using :py:func:`pvlib.temperature.GenericLinearModel()`
+
+glm = GenericLinearModel(module_efficiency=0.19, absorptance=0.88)
+glm = glm.use_faiman(u0=params['u0'], u1=params['u1'])
+
+# %% PVsyst paramaters
+
+params_pvsyst = glm.to_pvsyst()
+
+for k, v in params_pvsyst.items():
+    print('%-17s = %5.2f' % (k, v))
+
+# %% SAPM parameters
+
+params_sapm = glm.to_sapm()
+
+for k, v in params_sapm.items():
+    print('%-1s = %5.2f' % (k, v))

docs/sphinx/source/reference/pv_modeling/temperature.rst

@@ -12,6 +12,8 @@ PV temperature models
    temperature.sapm_cell_from_module
    temperature.pvsyst_cell
    temperature.faiman
+   temperature.faiman_dyn
+   temperature.fit_faiman_dyn
    temperature.faiman_rad
    temperature.fuentes
    temperature.ross

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

@@ -12,6 +12,10 @@ Enhancements
 * :py:func:`pvlib.bifacial.infinite_sheds.get_irradiance` and
   :py:func:`pvlib.bifacial.infinite_sheds.get_irradiance_poa` now include
   shaded fraction in returned variables. (:pull:`1871`)
+* Add a simple dynamic version of the Faiman temperature model
+  :py:func:`pvlib.temperature.faiman_dyn`
+  and a function to obtain parameters from field data
+  :py:func:`pvlib.temperature.faiman_dyn`(:issue:`1877`, :pull:`1878`)
 
 Bug fixes
 ~~~~~~~~~

pvlib/temperature.py

@@ -6,8 +6,9 @@
 import numpy as np
 import pandas as pd
 from pvlib.tools import sind
-from pvlib._deprecation import warn_deprecated
+# from pvlib._deprecation import warn_deprecated
 from pvlib.tools import _get_sample_intervals
+from scipy.optimize import curve_fit
 import scipy
 import scipy.constants
 import warnings
@@ -459,6 +460,248 @@ def faiman(poa_global, temp_air, wind_speed=1.0, u0=25.0, u1=6.84):
     return temp_air + temp_difference
 
 
+def faiman_dyn(poa_global, temp_air, wind_speed=1.0, u0=25.0, u1=6.84,
+               thermal_inertia=0):
+    r'''
+    Calculate cell or module temperature using the Faiman model augmented
+    with a thermal inertia term.
+
+    The Faiman model is an empirical heat loss factor model [1]_ and is
+    adopted in the IEC 61853 standards [2]_ and [3]_.  The method of
+    accounting for thermal inertia was proposed and demonstrated
+    by Driesse in [4]_.
+
+    The model can be used to represent cell or module temperature.
+
+    Parameters
+    ----------
+    poa_global : pandas.Series with DatetimeIndex
+        Total incident irradiance [W/m^2].
+
+    temp_air : numeric
+        Ambient dry bulb temperature [C].
+
+    wind_speed : numeric, default 1.0
+        Wind speed in m/s measured at the same height for which the wind loss
+        factor was determined.  The default value 1.0 m/s is the wind
+        speed at module height used to determine NOCT. [m/s]
+
+    u0 : numeric, default 25.0
+        Combined heat loss factor coefficient. The default value is one
+        determined by Faiman for 7 silicon modules
+        in the Negev desert on an open rack at 30.9° tilt.
+        :math:`\left[\frac{\text{W}/{\text{m}^2}}{\text{C}}\right]`
+
+    u1 : numeric, default 6.84
+        Combined heat loss factor influenced by wind. The default value is one
+        determined by Faiman for 7 silicon modules
+        in the Negev desert on an open rack at 30.9° tilt.
+        :math:`\left[ \frac{\text{W}/\text{m}^2}{\text{C}\ \left( \text{m/s} \right)} \right]`
+
+    thermal_inertia : scalar numeric, default 0.0
+        Represents the delay for a change in module temperature
+        in response to changes in operating conditions.
+        Typical values are 5-10 minutes. [minutes]
+
+    Returns
+    -------
+    numeric, values in degrees Celsius
+
+    Notes
+    -----
+    Thermal inertia is simulated using a simple moving average on the
+    operating conditions.
+
+
+    References
+    ----------
+    .. [1] Faiman, D. (2008). "Assessing the outdoor operating temperature of
+       photovoltaic modules." Progress in Photovoltaics 16(4): 307-315.
+       :doi:`10.1002/pip.813`
+
+    .. [2] "IEC 61853-2 Photovoltaic (PV) module performance testing and energy
+       rating - Part 2: Spectral responsivity, incidence angle and module
+       operating temperature measurements". IEC, Geneva, 2018.
+
+    .. [3] "IEC 61853-3 Photovoltaic (PV) module performance testing and energy
+       rating - Part 3: Energy rating of PV modules". IEC, Geneva, 2018.
+
+    .. [4] Driesse, A. (2022) "Module operating temperature model parameter
+       determination." Sandia National Laboratories, Albuquerque NM.
+       :doi:`10.5281/zenodo.10003736`
+
+    See also
+    --------
+    pvlib.temperature.fit_faiman_dyn
+    pvlib.temperature.faiman
+
+    '''   # noQA: E501
+
+    if not isinstance(poa_global, pd.Series):
+        raise ValueError('poa_global must be a pandas Series')
+
+    if not isinstance(poa_global.index, pd.DatetimeIndex):
+        raise ValueError('poa_global must have a Datetime index')
+
+    # assume the first time increment represent the series time step
+    # later gaps or large steps will not cause the function to fail
+    # only some minor glitches in the module temperature may result
+    time = poa_global.index
+    timestep = time[1] - time[0]
+    timestep_minutes = timestep.seconds / 60
+
+    window = round(2 * thermal_inertia / timestep_minutes)
+
+    if window > 1:
+        roll_options = dict(window=window, min_periods=1, center=False)
+
+        poa_global = pd.Series(poa_global).rolling(**roll_options).mean()
+        temp_air = pd.Series(temp_air).rolling(**roll_options).mean()
+        wind_speed = pd.Series(wind_speed).rolling(**roll_options).mean()
+
+    total_loss_factor = u0 + u1 * wind_speed
+    heat_input = poa_global
+    temp_difference = heat_input / total_loss_factor
+    return temp_air + temp_difference
+
+
+def fit_faiman_dyn(temp_pv, poa_global, temp_air, wind_speed,
+                   thermal_inertia=(0.0, 15.0, 1.0),
+                   full_output=False, **kwargs):
+    '''
+    Determine the optimal parameters for the faiman_dyn temperature model
+    from a set of measurements using a method proposed and demonstrated
+    by Driesse in [1]_.
+
+    Parameters
+    ----------
+    temp_pv : numeric
+        Cell or module temperature [C].
+
+    poa_global : pandas.Series with DatetimeIndex
+        Total incident irradiance [W/m^2].
+
+    temp_air : numeric
+        Ambient dry bulb temperature [C].
+
+    wind_speed : numeric, default 1.0
+        Wind speed in m/s measured at the same height for which the wind loss
+        factor was determined.  The default value 1.0 m/s is the wind
+        speed at module height used to determine NOCT. [m/s]
+
+    thermal_inertia : tuple of numeric, default (0.0, 15.0, 1.0)
+        Represents the delay for a change in module temperature
+        in response to changes in operating conditions.
+        This tuple specifies the range of values to be evaluated
+        in the fitting process: (minimum value, maximum value, increment).
+        [minutes]
+
+    full_output : boolean, default False
+        non-zero to return optional output
+
+    kwargs :
+        Additional arguments passed to scipy curvefit.
+
+    Returns
+    -------
+    p: dict
+        A dictionary that has the three keys u0, 01, thermal_inertia
+        and contains the optimal parameter values.
+
+    results : pandas DataFrame
+        Optional table of model parameters for each value of
+        thermal_inertia that was evaluated.
+
+    Notes
+    -----
+    Thermal inertia is simulated using a simple moving average of the
+    operating conditions. The fitting process is repeated for
+    the range of values specified by the thermal_inertia parameter.
+    From the collected results, the set of parameters producing the
+    smallest RMSE is identified as optimal.
+
+    References
+    ----------
+    .. [1] Driesse, A. (2022) "Module operating temperature model parameter
+       determination." Sandia National Laboratories, Albuquerque NM.
+       :doi:`10.5281/zenodo.10003736`
+
+    See also
+    --------
+    pvlib.temperature.faiman_dyn
+    pvlib.temperature.faiman
+    pvlib.temperature.GenericLinearModel
+    '''
+
+    if not isinstance(poa_global, pd.Series):
+        raise ValueError('poa_global must be a pandas Series')
+
+    if not isinstance(poa_global.index, pd.DatetimeIndex):
+        raise ValueError('poa_global must have a Datetime index')
+
+    # assume the first time increment represent the series time step
+    # later gaps or large steps will not cause the function to fail
+    # only some minor glitches in the module temperature may result
+    time = poa_global.index
+    timestep = time[1] - time[0]
+    timestep_minutes = timestep.seconds / 60
+
+    # prepare for fitting
+    def fitfun(xdata, *params):
+        return faiman(*xdata, *params)
+
+    fit_options = dict(p0=[20, 5])
+    fit_options.update(kwargs)
+
+    min_window = max(1, round(2 * thermal_inertia[0] / timestep_minutes))
+    max_window = max(1, round(2 * thermal_inertia[1] / timestep_minutes))
+    window_inc = max(1, round(2 * thermal_inertia[2] / timestep_minutes))
+
+    results = []
+
+    # repeat fits using a range of averaging windows
+    for window in range(min_window, max_window + 1, window_inc):
+
+        roll_options = dict(window=window, center=False)
+
+        poa_global_mean = pd.Series(poa_global).rolling(**roll_options).mean()
+        temp_air_mean = pd.Series(temp_air).rolling(**roll_options).mean()
+        wind_speed_mean = pd.Series(wind_speed).rolling(**roll_options).mean()
+
+        valid = True
+        for v in temp_pv, poa_global_mean, temp_air_mean, wind_speed_mean:
+            valid &= np.isfinite(v)
+
+        xdata = [poa_global_mean[valid],
+                 temp_air_mean[valid],
+                 wind_speed_mean[valid]]
+        ydata = temp_pv[valid]
+
+        popt, pcov = curve_fit(fitfun, xdata, ydata, **fit_options)
+
+        # calculate RMSE
+        ymodel = fitfun(xdata, *popt)
+        rmse = np.sqrt(np.nanmean(np.square(ymodel - ydata)))
+
+        inertia = (window / 2) * timestep_minutes
+        results.append((window, *popt, inertia, rmse))
+
+    # organize the results
+    param_names = ['u0', 'u1', 'thermal_inertia']
+    result_columns = ['window', *param_names, 'rmse']
+
+    results = pd.DataFrame(results, columns=result_columns)
+
+    # choose the solution with the smallest RMSE
+    best_window = results.rmse.idxmin()
+    best_params = results.loc[best_window, param_names].to_dict()
+
+    if full_output:
+        return best_params, results
+    else:
+        return best_params
+
+
 def faiman_rad(poa_global, temp_air, wind_speed=1.0, ir_down=None,
                u0=25.0, u1=6.84, sky_view=1.0, emissivity=0.88):
     r'''

pvlib/tests/test_temperature.py

@@ -127,6 +127,72 @@ def test_faiman_ndarray():
     assert_allclose(expected, result, atol=0.001)
 
 
+def test_faiman_dyn_good():
+    time = pd.date_range('3/14/15 9:00', periods=13, freq='5T')
+    df = pd.DataFrame(index=time)
+
+    df['poa_global'] = 0
+    df['poa_global'][3:] = 1000
+    df['temp_air'] = 25
+    df['temp_air'][8:] = -15
+    df['wind_speed'] = 1
+
+    # without thermal inertia
+    temp_pv = temperature.faiman_dyn(df.poa_global, df.temp_air, df.wind_speed,
+                                     u0=20, u1=5, thermal_inertia=0)
+    expected = [25, 25, 25, 65, 65, 65, 65, 65, 25, 25, 25, 25, 25]
+    assert_allclose(expected, temp_pv, atol=0.001)
+
+    # with thermal inertia
+    temp_pv = temperature.faiman_dyn(df.poa_global, df.temp_air, df.wind_speed,
+                                     u0=20, u1=5, thermal_inertia=10)
+    expected = [25, 25, 25, 35, 45, 55, 65, 65, 55, 45, 35, 25, 25]
+    assert_allclose(expected, temp_pv, atol=0.001)
+
+
+def test_faiman_dyn_bad():
+    with pytest.raises(ValueError, match='pandas Series'):
+        temperature.faiman_dyn(1000, 25)
+
+    with pytest.raises(ValueError, match='Datetime index'):
+        temperature.faiman_dyn(pd.Series(1000), 25)
+
+
+def test_fit_faiman_dyn_bad():
+    with pytest.raises(ValueError, match='pandas Series'):
+        temperature.fit_faiman_dyn(50, 1000, 25, 1)
+
+    with pytest.raises(ValueError, match='Datetime index'):
+        temperature.fit_faiman_dyn(50, pd.Series(1000), 25, 1)
+
+
+def test_fit_faiman_dyn_good():
+
+    time = pd.date_range('3/14/15 9:00', periods=30, freq='2T')
+    df = pd.DataFrame(index=time)
+
+    df['poa_global'] = np.random.uniform(0, 1000, len(time))
+    df['temp_air'] = np.random.uniform(10, 30, len(time))
+    df['wind_speed'] = np.random.uniform(0, 5, len(time))
+
+    expected = dict(u0=26.5, u1=3.5, thermal_inertia=8)
+    temp_pv = temperature.faiman_dyn(df.poa_global, df.temp_air, df.wind_speed,
+                                     **expected)
+
+    params = temperature.fit_faiman_dyn(temp_pv, df.poa_global,
+                                        df.temp_air, df.wind_speed,
+                                        thermal_inertia=(0, 15, 1))
+    for k, v in expected.items():
+        assert_allclose(params[k], v, atol=1e-4)
+
+    result = temperature.fit_faiman_dyn(temp_pv, df.poa_global,
+                                        df.temp_air, df.wind_speed,
+                                        thermal_inertia=(0, 15, 1),
+                                        full_output=True)
+    params, detail = result
+    assert isinstance(detail, pd.DataFrame)
+
+
 def test_faiman_rad_no_ir():
     expected = temperature.faiman(900, 20, 5)
     result = temperature.faiman_rad(900, 20, 5)