diff --git a/docs/sphinx/source/reference/pv_modeling/system_models.rst b/docs/sphinx/source/reference/pv_modeling/system_models.rst index 11137bead2..459a1083f9 100644 --- a/docs/sphinx/source/reference/pv_modeling/system_models.rst +++ b/docs/sphinx/source/reference/pv_modeling/system_models.rst @@ -56,3 +56,4 @@ PVGIS model :toctree: ../generated/ pvarray.huld + pvarray.fit_huld diff --git a/docs/sphinx/source/whatsnew/v0.11.0.rst b/docs/sphinx/source/whatsnew/v0.11.0.rst index 37cbad3bf9..8d95672dfe 100644 --- a/docs/sphinx/source/whatsnew/v0.11.0.rst +++ b/docs/sphinx/source/whatsnew/v0.11.0.rst @@ -9,6 +9,7 @@ Breaking changes ~~~~~~~~~~~~~~~~ * The deprecated ``pvlib.modelchain.basic_chain`` has now been removed. (:pull:`1862`) * Remove the `poa_horizontal_ratio` function and all of its references. (:issue:`1697`, :pull:`2021`) +* Parameter ``temp_mod`` in :py:func:`pvlib.pvarray.huld` is renamed to ``temp_module``. (:pull:`1979`) * The `leap_day` parameter in :py:func:`~pvlib.iotools.get_psm3` now defaults to True instead of False. (:issue:`1481`, :pull:`1991`) @@ -23,7 +24,7 @@ Enhancements shade perpendicular to ``axis_azimuth``. The function is applicable to both fixed-tilt and one-axis tracking systems. (:issue:`1689`, :pull:`1725`, :pull:`1962`) - +* Added a method to fit the Huld PV model. (:pull:`1979`) Bug fixes ~~~~~~~~~ diff --git a/pvlib/pvarray.py b/pvlib/pvarray.py index ab7530f3e5..e79b96f8a7 100644 --- a/pvlib/pvarray.py +++ b/pvlib/pvarray.py @@ -242,7 +242,7 @@ def _infer_k_huld(cell_type, pdc0): return k -def huld(effective_irradiance, temp_mod, pdc0, k=None, cell_type=None): +def huld(effective_irradiance, temp_module, pdc0, k=None, cell_type=None): r""" Power (DC) using the Huld model. @@ -263,13 +263,13 @@ def huld(effective_irradiance, temp_mod, pdc0, k=None, cell_type=None): ---------- effective_irradiance : numeric The irradiance that is converted to photocurrent. [:math:`W/m^2`] - temp_mod: numeric + temp_module: numeric Module back-surface temperature. [C] pdc0: numeric Power of the modules at reference conditions 1000 :math:`W/m^2` and :math:`25^{\circ}C`. [W] k : tuple, optional - Empirical coefficients used in the power model. Length 6. If ``k`` is + Empirical coefficients used in the Huld model. Length 6. If ``k`` is not provided, ``cell_type`` must be specified. cell_type : str, optional If provided, must be one of ``'cSi'``, ``'CIS'``, or ``'CdTe'``. @@ -332,6 +332,10 @@ def huld(effective_irradiance, temp_mod, pdc0, k=None, cell_type=None): E. Dunlop. A power-rating model for crystalline silicon PV modules. Solar Energy Materials and Solar Cells 95, (2011), pp. 3359-3369. :doi:`10.1016/j.solmat.2011.07.026`. + + See Also + -------- + pvlib.pvarray.fit_huld """ if k is None: if cell_type is not None: @@ -340,7 +344,7 @@ def huld(effective_irradiance, temp_mod, pdc0, k=None, cell_type=None): raise ValueError('Either k or cell_type must be specified') gprime = effective_irradiance / 1000 - tprime = temp_mod - 25 + tprime = temp_module - 25 # accomodate gprime<=0 with np.errstate(divide='ignore'): logGprime = np.log(gprime, out=np.zeros_like(gprime), @@ -350,3 +354,99 @@ def huld(effective_irradiance, temp_mod, pdc0, k=None, cell_type=None): k[2] * tprime + k[3] * tprime * logGprime + k[4] * tprime * logGprime**2 + k[5] * tprime**2) return pdc + + +def _build_iec61853(): + ee = np.array([100, 100, 200, 200, 400, 400, 400, 600, 600, 600, 600, + 800, 800, 800, 800, 1000, 1000, 1000, 1000, 1100, 1100, + 1100], dtype=np.float64).T + tc = np.array([15, 25, 15, 25, 15, 25, 50, 15, 25, 50, 75, + 15, 25, 50, 75, 15, 25, 50, 75, 25, 50, 75], + dtype=np.float64).T + return ee, tc + + +def fit_huld(effective_irradiance, temp_module, pdc, method='ols'): + r''' + Fit the Huld model to the input data. + + Parameters + ---------- + effective_irradiance : array-like + The irradiance that is converted to photocurrent. [:math:`W/m^2`] + temp_module: array-like + Module back-surface temperature. [C] + pdc: array-like + DC power at ``effective_irradiance`` and ``temp_module``. [W] + method : string, default 'ols' + The regression method. Must be either ``'ols'`` or ``'robust'``. + + Returns + ------- + pdc0: numeric + Power of the modules at reference conditions 1000 :math:`W/m^2` + and :math:`25^{\circ}C`. [W] + k : tuple + Empirical coefficients used in the Huld model. Length 6. + + Notes + ----- + Requires ``statsmodels``. + + Input data ``effective_irradiance``, ``temp_module`` and ``pdc`` must be + equal length. Data are aligned as columns and any row with a NaN in any + column is dropped. Remaining data must be at least length 7. + + Huld [1]_ describes fitting by least-squares but doesn't detail the + procedure, although it is reasonable that ordinary least-squares regression + was used. Ordinary least-squares regression can be unduly influenced by + data points with large deviations from the mean of the dependent variable + ``pdc``. For example, these data may occur when short-duration shading is + present. If filtering is applied to exclude such data, a suitable model + can be obtained with ordinary least-squares regression. As an alternative, + this function also provides a 'robust' regression option (an iteratively + reweighted least-squares method) that is more tolerant of outliers in the + dependent variable. + + References + ---------- + .. [1] T. Huld, G. Friesen, A. Skoczek, R. Kenny, T. Sample, M. Field, + E. Dunlop. A power-rating model for crystalline silicon PV modules. + Solar Energy Materials and Solar Cells 95, (2011), pp. 3359-3369. + :doi:`10.1016/j.solmat.2011.07.026`. + + See Also + -------- + pvlib.pvarray.huld + ''' + + try: + import statsmodels.api as sm + except ImportError: + raise ImportError( + 'fit_huld requires statsmodels') + + gprime = effective_irradiance / 1000 + tprime = temp_module - 25 + # accomodate gprime<=0 + logGprime = np.log(gprime, out=np.nan*gprime, + where=np.array(gprime > 0)) + with np.errstate(divide='ignore'): + Y = np.divide(pdc, gprime, out=np.nan*gprime, + where=np.array(gprime > 0)) + + X = np.stack((logGprime, logGprime**2, tprime, tprime*logGprime, + tprime*logGprime**2, tprime**2), axis=0).T + X = sm.add_constant(X) + + if method == 'ols': + model = sm.OLS(Y, X, missing='drop') + elif method == 'robust': + model = sm.RLM(Y, X, missing='drop') + else: + raise ValueError("method must be ols or robust") + + result = model.fit() + pdc0 = result.params[0] + k = result.params[1:] + return pdc0, k diff --git a/pvlib/tests/test_pvarray.py b/pvlib/tests/test_pvarray.py index 693ef78b2a..309c492c53 100644 --- a/pvlib/tests/test_pvarray.py +++ b/pvlib/tests/test_pvarray.py @@ -3,7 +3,7 @@ from numpy.testing import assert_allclose from .conftest import assert_series_equal import pytest - +from pvlib.tests.conftest import requires_statsmodels from pvlib import pvarray @@ -69,3 +69,40 @@ def test_huld(): with pytest.raises(ValueError, match='Either k or cell_type must be specified'): res = pvarray.huld(1000, 25, 100) + + +@pytest.mark.parametrize('method', ['ols', 'robust']) +@requires_statsmodels +def test_fit_huld(method): + # test is to recover the parameters in _infer_huld_k for each cell type + # IEC61853 conditions to make data for fitting + ee, tc = pvarray._build_iec61853() + techs = ['csi', 'cis', 'cdte'] + pdc0 = 250 + for tech in techs: + k0 = pvarray._infer_k_huld(tech, pdc0) + pdc = pvarray.huld(ee, tc, pdc0, cell_type=tech) + m_pdc0, k = pvarray.fit_huld(ee, tc, pdc, method=method) + expected = np.array([pdc0, ] + [v for v in k0], dtype=float) + modeled = np.hstack((m_pdc0, k)) + assert_allclose(expected, modeled, rtol=1e-8) + # once more to check that NaNs are handled + ee[7] = np.nan + tc[9] = np.nan + k0 = pvarray._infer_k_huld('csi', pdc0) + pdc = pvarray.huld(ee, tc, pdc0, cell_type='csi') + pdc[11] = np.nan + m_pdc0, k = pvarray.fit_huld(ee, tc, pdc, method='ols') + expected = np.array([pdc0, ] + [v for v in k0], dtype=float) + modeled = np.hstack((m_pdc0, k)) + assert_allclose(expected, modeled, rtol=1e-8) + + +@requires_statsmodels +def test_fit_huld_method_error(): + ee, tc = pvarray._build_iec61853() + pdc0 = 250 + pdc = pvarray.huld(ee, tc, pdc0, cell_type='csi') + method = 'brute_force' + with pytest.raises(ValueError, match="method must be ols or robust"): + m_pdc0, k = pvarray.fit_huld(ee, tc, pdc, method=method)