Maximum power point tracking for mismatched devices (single-diode model)

pvlib/pvsystem.py

@@ -12,6 +12,7 @@
 from urllib.request import urlopen
 import numpy as np
 from scipy import constants
+import scipy.optimize
 import pandas as pd
 from dataclasses import dataclass
 from abc import ABC, abstractmethod
@@ -3038,3 +3039,91 @@ def combine_loss_factors(index, *losses, fill_method='ffill'):
         combined_factor *= (1 - loss)
 
     return 1 - combined_factor
+
+
+def _negative_total_power(current, *args):
+    """
+    Compute negative of total power generated by devices in series at
+    specified current. Designed for use by scipy.optimize.minimize.
+    """
+    return -np.sum(current * v_from_i(current, *args))
+
+
+def max_power_point_mismatched(
+    photocurrent,
+    saturation_current,
+    resistance_series,
+    resistance_shunt,
+    nNsVth,
+    *,
+    i_mp_ic=None,
+):
+    """
+    FIXME Replace this with proper docstring.
+
+    Compute maximum power info for (possibly) mismatched set of devices in
+    series. When using this serially on time-series data, passing i_mp_ic from
+    previous step may speed up computation. Algorithm falls back to automated
+    computation of i_mp_ic if solution fails with provided i_mp_ic. The value
+    of i_mp_ic used is returned along with i_mp (same value for all devices),
+    v_mp, p_mp, i_mp_string, v_mp_string, and p_mp_string.
+    """
+    if i_mp_ic is None:
+        retry_ic = False
+
+        i_mp_ic = max_power_point(
+            np.mean(photocurrent),
+            np.mean(saturation_current),
+            np.mean(resistance_series),
+            np.mean(resistance_shunt),
+            np.mean(nNsVth),
+        )["i_mp"]
+    else:
+        retry_ic = True
+
+    args = (
+        photocurrent,
+        saturation_current,
+        resistance_series,
+        resistance_shunt,
+        nNsVth,
+    )
+
+    sol = scipy.optimize.minimize(
+        _negative_total_power, i_mp_ic, args=args, jac='3-point'
+    )
+
+    if sol.success:
+        i_mp = sol.x[0]
+        v_mp = v_from_i(
+            i_mp,
+            photocurrent,
+            saturation_current,
+            resistance_series,
+            resistance_shunt,
+            nNsVth,
+        )
+
+        return {
+            "i_mp_ic": i_mp_ic,
+            "i_mp": i_mp,
+            "v_mp": v_mp,
+            "p_mp": i_mp * v_mp,
+            "i_mp_string": i_mp,
+            "v_mp_string": np.sum(v_mp),
+            "p_mp_string": -sol.fun,
+        }
+
+    if retry_ic:
+        # Try solution one more time using automated inital condition.
+        # Caller can detect this occurance by seeing change in i_mp_ic in
+        # return value.
+        return max_power_point_mismatched(
+            photocurrent,
+            saturation_current,
+            resistance_series,
+            resistance_shunt,
+            nNsVth,
+        )
+
+    raise RuntimeError(f"unsuccessful solution: {sol}")

pvlib/tests/test_pvsystem.py

@@ -4,6 +4,8 @@
 import numpy as np
 from numpy import nan, array
 import pandas as pd
+import scipy.constants
+import scipy.optimize
 
 import pytest
 from .conftest import assert_series_equal, assert_frame_equal
@@ -2509,3 +2511,141 @@ def test_Array_temperature_missing_parameters(model, keys):
         array.temperature_model_parameters = params
         with pytest.raises(KeyError, match=match):
             array.get_cell_temperature(irrads, temps, winds, model)
+
+
+@pytest.mark.parametrize(
+    'inputs',
+    [
+        # FIXME Need many more argument combinations in additional test cases.
+        {
+            "photocurrent": 6.2,
+            "saturation_current": 1.0e-8,
+            "n": 1.1,
+            "resistance_series": 0.0001,
+            "resistance_shunt": 5000.0,
+            "Ns": 60,
+            "T": 25.0,
+        },
+        {
+            "photocurrent": 6.2,
+            "saturation_current": 1.0e-8,
+            "n": 1.1,
+            "resistance_series": 0.0001,
+            "resistance_shunt": 5000.0,
+            "Ns": 60,
+            "T": 25.0,
+            "i_mp_ic": 5.8,
+        },
+        {
+            "photocurrent": np.array([5.8, 6.2]),
+            "saturation_current": 1.0e-8,
+            "n": 1.1,
+            "resistance_series": 0.0001,
+            "resistance_shunt": 5000.0,
+            "Ns": 60,
+            "T": 25.0,
+        },
+        {
+            "photocurrent": np.array([5.8, 6.2]),
+            "saturation_current": 1.0e-8,
+            "n": 1.1,
+            "resistance_series": 0.0001,
+            "resistance_shunt": 5000.0,
+            "Ns": 60,
+            "T": 25.0,
+            "i_mp_ic": None,
+        },
+        {
+            "photocurrent": np.array([5.8, 6.2]),
+            "saturation_current": 1.0e-8,
+            "n": 1.1,
+            "resistance_series": 0.0001,
+            "resistance_shunt": 5000.0,
+            "Ns": 60,
+            "T": 25.0,
+            "i_mp_ic": -1.0e14,
+        },
+        {
+            "photocurrent": np.array([5.8, 6.2]),
+            "saturation_current": 1.0e-8,
+            "n": 1.1,
+            "resistance_series": 0.0001,
+            "resistance_shunt": 5000.0,
+            "Ns": 60,
+            "T": 25.0,
+            "i_mp_ic": 5.6,
+        },
+    ]
+)
+def test_max_power_point_mismatched(inputs):
+    """Test max power point computation for mismatched devices in series."""
+
+    photocurrent = inputs["photocurrent"]
+    saturation_current = inputs["saturation_current"]
+    resistance_series = inputs["resistance_series"]
+    resistance_shunt = inputs["resistance_shunt"]
+    q_C = scipy.constants.value("elementary charge")
+    k_B_J_per_K = scipy.constants.value("Boltzmann constant")
+    T_K = scipy.constants.convert_temperature(inputs["T"], "Celsius", "Kelvin")
+    nNsVth = inputs["n"] * inputs["Ns"] * k_B_J_per_K * T_K / q_C
+
+    if "i_mp_ic" in inputs:
+        result = pvsystem.max_power_point_mismatched(
+            photocurrent,
+            saturation_current,
+            resistance_series,
+            resistance_shunt,
+            nNsVth,
+            i_mp_ic=inputs["i_mp_ic"],
+        )
+    else:
+        result = pvsystem.max_power_point_mismatched(
+            photocurrent,
+            saturation_current,
+            resistance_series,
+            resistance_shunt,
+            nNsVth,
+        )
+
+    # FIXME Replace this with test assertions.
+    print(result)
+
+
+def test_max_power_point_mismatched_unsuccessful_solver(monkeypatch):
+    """
+    Test mismatched max power point computation where solver is unsuccessful.
+    """
+    photocurrent = 6.2
+    saturation_current = 1.0e-8
+    resistance_series = 0.0001
+    resistance_shunt = 5000.0
+    n = 1.1
+    Ns = 60
+    T = 25.0
+    T_K = scipy.constants.convert_temperature(T, "Celsius", "Kelvin")
+    k_B_J_per_K = scipy.constants.value("Boltzmann constant")
+    q_C = scipy.constants.value("elementary charge")
+    nNsVth = n * Ns * k_B_J_per_K * T_K / q_C
+
+    def minimize_monkeypatched(*_, **__):
+        """Return an unsuccessful solution from solver."""
+        return scipy.optimize.OptimizeResult(success=False)
+
+    # Monkey patch solver to return unsuccessfully.
+    monkeypatch.setattr(
+        scipy.optimize,
+        "minimize",
+        minimize_monkeypatched,
+    )
+
+    with pytest.raises(RuntimeError) as e_info:
+        pvsystem.max_power_point_mismatched(
+            photocurrent,
+            saturation_current,
+            resistance_series,
+            resistance_shunt,
+            nNsVth,
+        )
+
+    # FIXME Replace this with test assertions.
+    print(e_info)