add fit_huld

pvlib · cwhanse · Feb 23, 2024 · Feb 23, 2024 · Feb 23, 2024 · Feb 23, 2024
commit 55c95305c63b26b6476f35b1551581332765dd5b
diff --git a/pvlib/pvarray.py b/pvlib/pvarray.py
@@ -11,6 +11,7 @@
 import numpy as np
 from scipy.optimize import curve_fit
 from scipy.special import exp10
+import statsmodels.api as sm
 
 
 def pvefficiency_adr(effective_irradiance, temp_cell,
@@ -348,3 +349,33 @@ 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]).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]).T
+    return ee, tc
+
+
+def fit_huld(effective_irradiance, temp_mod, pdc):
+    gprime = effective_irradiance / 1000
+    tprime = temp_mod - 25
+    # accomodate gprime<=0
+    with np.errstate(divide='ignore'):
+        logGprime = np.log(gprime, out=np.zeros_like(gprime),
+                           where=np.array(gprime > 0))
+        Y = np.divide(pdc, gprime, out=np.zeros_like(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)
+
+    rlm_model = sm.RLM(Y, X)
+    rlm_result = rlm_model.fit()
+    pdc0 = rlm_result.params[0]
+    k = rlm_result.params[1:]
+    return pdc0, k
diff --git a/pvlib/tests/test_pvarray.py b/pvlib/tests/test_pvarray.py
@@ -69,3 +69,18 @@ def test_huld():
     with pytest.raises(ValueError,
                        match='Either k or cell_type must be specified'):
         res = pvarray.huld(1000, 25, 100)
+
+
+def test_fit_huld():
+    # 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)
+        expected = np.array([pdc0,] + [v for v in k0], dtype=float)
+        modeled = np.hstack((m_pdc0, k))
+        assert_allclose(expected, modeled, rtol=1e-8)