add marion's adjustment to pvwatts_dc

pvlib/pvsystem.py

@@ -2878,14 +2878,34 @@ def scale_voltage_current_power(data, voltage=1, current=1):
 
 @renamed_kwarg_warning(
     "0.13.0", "g_poa_effective", "effective_irradiance")
-def pvwatts_dc(effective_irradiance, temp_cell, pdc0, gamma_pdc, temp_ref=25.):
+def pvwatts_dc(effective_irradiance, temp_cell, pdc0, gamma_pdc, temp_ref=25.,
+               k=None, cap_adjustment=False):
     r"""
-    Implements NREL's PVWatts DC power model. The PVWatts DC model [1]_ is:
+    Implements NREL's PVWatts (Version 5) DC power model. The PVWatts Version
+    5 DC model [1]_ is:
 
     .. math::
 
         P_{dc} = \frac{G_{poa eff}}{1000} P_{dc0} ( 1 + \gamma_{pdc} (T_{cell} - T_{ref}))
 
+    This model has also been referred to as the power temperature coefficient
+    model.
+
+    This function accepts an optional irradiance adjustment factor, `k`, based
+    on [2]_. This applies a piece-wise adjustment to power based on irradiance,
+    where `k` is the reduction in actual power at 200 W/m^2 relative to power
+    calculated at 200 W/m^2 as 0.2*`pdc0`. For example, a 500 W module that
+    produces 95 W at 200 W/m^2 (a 5% relative reduction in efficiency) would
+    have a value of `k` = 0.01.
+
+    .. math::
+
+        k=\frac{0.2P_{dc0}-P_{200}}{P_{dc0}}
+
+    This adjustment increases relative efficiency for irradiance above 1000
+    W/m^2, which may not be desired. An optional input, `capped_adjustment`,
+    modifies the adjustment from [2]_ to only apply below 1000 W/m^2.
+
     Note that ``pdc0`` is also used as a symbol in
     :py:func:`pvlib.inverter.pvwatts`. ``pdc0`` in this function refers to the DC
     power of the modules at reference conditions. ``pdc0`` in
@@ -2909,6 +2929,10 @@ def pvwatts_dc(effective_irradiance, temp_cell, pdc0, gamma_pdc, temp_ref=25.):
     temp_ref: numeric, default 25.0
         Cell reference temperature. PVWatts defines it to be 25 C and
         is included here for flexibility. [C]
+    k: numeric, optional
+        Irradiance correction factor, defined in [2]_. [unitless]
+    cap_adjustment: Boolean, default False
+        If True, apply the optional adjustment at and below 1000 W/m^2.
 
     Returns
     -------
@@ -2920,11 +2944,34 @@ def pvwatts_dc(effective_irradiance, temp_cell, pdc0, gamma_pdc, temp_ref=25.):
     .. [1] A. P. Dobos, "PVWatts Version 5 Manual"
            http://pvwatts.nrel.gov/downloads/pvwattsv5.pdf
            (2014).
+    .. [2] B. Marion, "Comparison of Predictive Models for
+           Photovoltaic Module Performance,"
+           https://doi.org/10.1109/PVSC.2008.4922586,
+           https://docs.nrel.gov/docs/fy08osti/42511.pdf
+           (2008).
     """  # noqa: E501
 
     pdc = (effective_irradiance * 0.001 * pdc0 *
            (1 + gamma_pdc * (temp_cell - temp_ref)))
 
+    # apply Marion's correction if k is anything but zero
+    if k is not None:
+        err_1 = (k * (1 - (1 - effective_irradiance / 200)**4) /
+                 (effective_irradiance / 1000))
+        err_2 = (k * (1000 - effective_irradiance) / (1000 - 200))
+
+        pdc_marion = np.where(effective_irradiance <= 200,
+                              pdc * (1 - err_1),
+                              pdc * (1 - err_2))
+
+        # "cap" Marion's correction at 1000 W/m^2
+        if cap_adjustment is True:
+            pdc_marion = np.where(effective_irradiance >= 1000,
+                                  pdc,
+                                  pdc_marion)
+
+        pdc = pdc_marion
+
     return pdc