fixes #50 - generates pvcell coeffs from iec 61853

.gitignore

@@ -26,6 +26,8 @@ _*/
 testPV/
 benchmark_*/
 .coverage
+.cache/
+__pycache__/
 
 # virtualenv
 venv/

.travis.yml

@@ -5,9 +5,6 @@ before_install:
     - sudo apt-get -qq update
     - sudo apt-get install -y gfortran
     - sudo apt-get install -y libatlas-dev liblapack-dev libblas-dev
-before_script:
-    - pip uninstall -y dulwich
-    - pip install -e git+https://github.com/mikofski/dulwich.git@versioneer#egg=dulwich
 script:
     - pytest
     - python setup.py bdist_wheel

pvmismatch/contrib/__init__.py

@@ -0,0 +1,4 @@
+"""
+Contributions to pvmismatch that are considered useful enough to be distributed
+but are not necessarily considered part of pvmismatch core.
+"""

pvmismatch/contrib/gen_coeffs/__init__.py

@@ -0,0 +1,210 @@
+"""
+Methods to generate diode coefficients.
+"""
+
+from pvlib.pvsystem import sapm
+import numpy as np
+from scipy import optimize
+from pvmismatch.contrib.gen_coeffs import diode, two_diode
+from pvmismatch.pvmismatch_lib.pvcell import ISAT1_T0, ISAT2, RS, RSH
+
+# IEC 61853 test matrix
+TC_C = [15.0, 25.0, 50.0, 75.0]
+IRR_W_M2 = [100.0, 200.0, 400.0, 600.0, 800.0, 1000.0, 1100.0]
+TEST_MAT = np.meshgrid(TC_C, IRR_W_M2)
+
+def gen_iec_61853_from_sapm(pvmodule):
+    """
+    Generate an IEC 61853 test from Sandia Array Performance Model (sapm).
+
+    :param pvmodule: PV module to be tested
+    :type pvmodule: dict
+
+    Module is a dictionary according to :def:`pvlib.pvsystem.sapm`.
+    """
+    tc, irr = TEST_MAT
+    return sapm(irr / 1000.0, tc, pvmodule)
+
+
+def gen_two_diode(isc, voc, imp, vmp, nseries, nparallel,
+                  tc, x0 = None, *args, **kwargs):
+    """
+    Generate two-diode model parameters for ``pvcell`` given 
+    """
+    isc_cell = isc / nparallel
+    voc_cell = voc / nseries
+    imp_cell = imp / nparallel
+    vmp_cell = vmp / nseries
+    if x0 is None:
+        isat1 = ISAT1_T0  # [A]
+        isat2 = ISAT2
+        rs = RS  # [ohms]
+        rsh = RSH  # [ohms]
+    else:
+        isat1 = x0[0]
+        isat2 = x0[1]
+        rs = x0[2]
+        rsh = x0[3]
+    x = np.array([np.log(isat1), np.log(isat2), np.sqrt(rs), np.sqrt(rsh)])
+    sol = optimize.root(
+        fun=residual_two_diode, x0=x,
+        args=(isc_cell, voc_cell, imp_cell, vmp_cell, tc),
+        jac=True,
+        *args, **kwargs
+    )
+    if sol.success:
+        isat1 = np.exp(sol.x[0])
+        isat2 = np.exp(sol.x[1])
+        rs = sol.x[2] ** 2.0
+        rsh = sol.x[3] ** 2.0
+    return (isat1, isat2, rs, rsh), sol
+
+
+def gen_sapm(iec_61853):
+    i_sc = iec_61853['i_sc']
+    # calculate Isc0 and alpha_Isc
+    # given Isc = Ee * Isc0 * (1 + alpha_Isc * (Tc - T0))
+    # as Ee * Isc0 + Ee * Isc0 * alpha_Isc * (Tc - T0) = Isc
+    # so Ax = B
+    # where x0 = Isc0 and x1 = Isc0 * alpha_Isc
+    # and A = [Ee, Ee * (Tc - T0)] and B = Isc
+    tc, irr = TEST_MAT
+    ee = (irr / 1000.0).flatten()
+    delta_t = (tc - 25.0).flatten()
+    a = np.array([ee, ee * delta_t])
+    b = i_sc.flatten()
+    x, res, rank, s = np.linalg.lstsq(a.T, b.T)
+    isc0, alpha_isc = x[0], x[1] / x[0]
+    return isc0, alpha_isc
+
+
+
+def residual_two_diode(x, isc, voc, imp, vmp, tc):
+    """
+    Objective function to solve 2-diode model.
+    :param x: parameters isat1, isat2, rs and rsh
+    :param isc: short circuit current [A] at tc [C]
+    :param voc: open circuit voltage [V] at tc [C]
+    :param imp: max power current [A] at tc [C]
+    :param vmp: max power voltage [V] at tc [C]
+    :param tc: cell temperature [C]
+    :return: norm of the residuals its sensitivity
+    """
+    # Constants
+    q = diode.QE  # [C/electron] elementary electric charge
+    # (n.b. 1 Coulomb = 1 A * s)
+    kb = diode.KB  # [J/K/molecule] Boltzmann's constant
+    tck = tc + 273.15  # [K] reference temperature
+    # Governing Equation
+    vt = kb * tck / q  # [V] thermal voltage
+    # Rescale Variables
+    isat1_t0 = np.exp(x[0])
+    isat2 = np.exp(x[1])
+    rs = x[2] ** 2.0
+    rsh = x[3] ** 2.0
+    # first diode saturation current
+    isat1 = diode.isat_t(tc, isat1_t0)
+    # Short Circuit
+    vd_isc, _ = diode.fvd(vc=0.0, ic=isc, rs=rs)
+    id1_isc, _ = diode.fid(isat=isat1, vd=vd_isc, m=1.0, vt=vt)
+    id2_isc, _ = diode.fid(isat=isat2, vd=vd_isc, m=2.0, vt=vt)
+    ish_isc, _ = diode.fish(vd=vd_isc, rsh=rsh)
+    # Photo-generated Current
+    iph = isc + id1_isc + id2_isc + ish_isc  # [A]
+    # Open Circuit
+    vd_voc, jvd_voc = diode.fvd(vc=voc, ic=0.0, rs=rs)
+    id1_voc, jid1_voc = diode.fid(isat=isat1, vd=vd_voc, m=1.0, vt=vt)
+    id2_voc, jid2_voc = diode.fid(isat=isat2, vd=vd_voc, m=2.0, vt=vt)
+    ish_voc, jish_voc = diode.fish(vd=vd_voc, rsh=rsh)
+    # Max Power Point
+    vd_mpp, jvd_mpp = diode.fvd(vc=vmp, ic=imp, rs=rs)
+    id1_mpp, jid1_mpp = diode.fid(isat=isat1, vd=vd_mpp, m=1.0, vt=vt)
+    id2_mpp, jid2_mpp = diode.fid(isat=isat2, vd=vd_mpp, m=2.0, vt=vt)
+    ish_mpp, jish_mpp = diode.fish(vd=vd_mpp, rsh=rsh)
+    # Slope at Max Power Point
+    dpdv, jdpdv = two_diode.fdpdv(
+        isat1=isat1, isat2=isat2, rs=rs, rsh=rsh, ic=imp, vc=vmp, vt=vt
+    )
+    # Shunt Resistance
+    frsh, jrsh = two_diode.fjrsh(
+        isat1=isat1, isat2=isat2, rs=rs, rsh=rsh, vt=vt, isc=isc
+    )
+    # Residual
+    # should be (M, ) array with M residual equations (constraints)
+    f2 = np.stack([
+        (iph - id1_voc - id2_voc - ish_voc).T,  # Open Circuit
+        (iph - id1_mpp - id2_mpp - ish_mpp - imp).T,  # Max Power Point
+        dpdv.T,  # Slope at Max Power Point
+        frsh.T  # Shunt Resistance
+    ], axis=0).flatten()
+    # Jacobian
+    # should be (M, N) array with M residuals and N variables
+    # [[df1/dx1, df1/dx2, ...], [df2/dx1, df2/dx2, ...]]
+    jvoc = np.stack((
+        -jid1_voc[0],  # d/disat1
+        -jid2_voc[0],  # d/disat2
+        -jvd_voc[2] * (jid1_voc[1] + jid2_voc[1] + jish_voc[0]),  # d/drs
+        -jish_voc[1]  # d/drsh
+    ), axis=0).T.reshape(-1, 4)
+    jmpp = np.stack((
+        -jid1_mpp[0],  # d/disat1
+        -jid2_mpp[0],  # d/disat2
+        -jvd_mpp[2] * (jid1_mpp[1] + jid2_mpp[1] + jish_mpp[0]),  # d/drs
+        -jish_mpp[1]  # d.drsh
+    ), axis=0).T.reshape(-1, 4) 
+    # Scaling Factors
+    scale_fx = np.array([np.exp(x[0]), np.exp(x[1]), 2 * x[2], 2 * x[3]])
+    # scales each column by the corresponding element
+    j2 = np.concatenate(
+        (jvoc, jmpp, jdpdv[:4].T.reshape(-1, 4), jrsh[:4].T.reshape(-1, 4)),
+        axis=0
+    ) * scale_fx
+    return f2, j2
+
+
+PVMODULES = {
+    "SunPower_SPR_E20_435": {
+        "Vintage": "2013-01-01",
+        "Area": 2.16,
+        "Material": "c-Si",
+        "Cells_in_Series": 128,
+        "Parallel_Strings": 1,
+        "Isco": 6.4293,
+        "Voco": 86.626,
+        "Impo": 6.0102,
+        "Vmpo": 72.3771,
+        "Aisc": 0.00037,
+        "Aimp": 7.02e-05,
+        "C0": 1.0115,
+        "C1": -0.0115,
+        "C2": 0.218474,
+        "C3": -7.224183,
+        "Bvoco": -0.248,
+        "Mbvoc": 0.0,
+        "Bvmpo": -0.2584,
+        "Mbvmp": 0.0,
+        "N": 1.011,
+        "A0": 0.957,
+        "A1": 0.0402,
+        "A2": -0.008515,
+        "A3": 0.0007141,
+        "A4": -2.132e-05,
+        "B0": 1.0002,
+        "B1": -0.000213,
+        "B2": 3.63416e-05,
+        "B3": -2.175e-06,
+        "B4": 5.2796e-08,
+        "B5": -4.4351e-10,
+        "DTC": 3.0,
+        "FD": 1.0,
+        "A": -3.46,
+        "B": -0.07599,
+        "IXO": 6.1717,
+        "IXXO": 4.3997,
+        "C4": 0.9891,
+        "C5": 0.0109,
+        "C6": 1.0869,
+        "C7": -0.0869,
+        "Notes": "Source: Estimated"
+    }
+}

pvmismatch/contrib/gen_coeffs/diode.py

@@ -0,0 +1,119 @@
+"""
+Common diode model equations.
+"""
+
+import numpy as np
+from scipy.constants import elementary_charge as QE, Boltzmann as KB
+
+T0 = 25.0  # [C]
+E0 = 1000.0  # [W/m^2]
+EG = 1.1  # [eV]]
+
+
+def fid(isat, vd, m, vt):
+    """
+    Diode current, I_d, and its derivatives w.r.t. I_sat, V_d, m and V_t.
+
+    :param isat: diode saturation current [A]
+    :type isat: float
+    :param vd: diode voltage [V]
+    :type vd: float
+    :param m: diode ideality factor
+    :type m: float
+    :param vt: thermal voltage [V]
+    :type vt: float
+    :return: diode current [A] and its derivatives
+    :rtype: float
+
+    Diode current is given by Shockley's equation ...
+
+    .. math::
+        I_d = I_{sat} \\left(\\exp\\left(\\frac{V_d}{m V_t} \\right) - 1\\right)
+
+    ... where I_d is the diode current in amps, I_sat is the saturation current
+    in amps, V_d is the diode voltage in volts, m is the diode ideality factor
+    and V_t is the thermal voltage in volts ...
+
+    .. math::
+        V_t = \\frac{k_B T}{q_e}
+
+    ... in which k_B is the Boltzmann constant, T is the ambient temperature in
+    Kelvin and q_e is teh elementary charge in coulombs per electron.
+
+    https://en.wikipedia.org/wiki/Shockley_diode_equation
+
+    """
+    vact = m * vt  # activation voltage [V]
+    growth = np.exp(vd / vact)  # growth term
+    expfact = (growth - 1.0)  # exponential factor
+    isat_growth = isat * growth  # combination parameter
+    vd_isat_growth = -vd * isat_growth  # combination parameter
+    # diode current
+    id_ = np.atleast_1d(isat * expfact)
+    # derivatives
+    d_isat = np.atleast_1d(expfact)  # df w.r.t. isat
+    d_vd = np.atleast_1d(isat_growth / vact)  # df w.r.t. vd
+    d_m = np.atleast_1d(vd_isat_growth / (m ** 2.0 * vt))  # df w.r.t. m
+    d_vt = np.atleast_1d(vd_isat_growth / (m * vt ** 2.0))  # df w.r.t. vt
+    jac = np.array([d_isat, d_vd, d_m, d_vt])  # jacobian
+    return id_, jac
+
+
+def fish(vd, rsh):
+    """
+    Shunt current, I_sh, and its derivatives w.r.t. V_d and R_sh.
+
+    :param vd: diode voltage [V]
+    :param rsh: shunt resistance [Ohms]
+    :return: shunt current [A]
+    """
+    ish = np.atleast_1d(vd / rsh)
+    shaper = np.ones(ish.shape)  # make sure scalars are the right shape
+    d_vd = np.atleast_1d(1.0 / rsh) * shaper
+    d_rsh = np.atleast_1d(-ish * d_vd)
+    jac = np.array([d_vd, d_rsh])
+    return ish, jac
+
+
+def fvd(vc, ic, rs):
+    """
+    Diode voltage, V_d, and its derivatives w.r.t. V_c, I_c, R_s.
+
+    :param vc: cell voltage [V]
+    :param ic: cell current [A]
+    :param rs: series resistance [Ohms]
+    :return: diode voltage [V]
+    """
+    vd = np.atleast_1d(vc + rs * ic)
+    shaper = np.ones(vd.shape)  # make sure scalars are the right shape
+    jac = np.array([np.ones(vd.shape),
+                    np.atleast_1d(rs) * shaper,  # d/dIc
+                    np.atleast_1d(ic) * shaper])  # d/dRs
+    return vd, jac
+
+
+def isat_t(tc, isat0):
+    tck = tc + 273.15
+    t0k = T0 + 273.15
+    tstar = tck ** 3.0 / t0k ** 3.0
+    inv_delta_t = 1.0 / t0k - 1.0 / tck  # [1/K]
+    exp_tstar = np.exp(EG * QE / KB * inv_delta_t)
+    isat_t = isat0 * tstar * exp_tstar
+    return isat_t
+
+
+def isc_t(tc, isc0, alpha_isc):
+    delta_t = tc - T0
+    isc_t = isc0 * (1.0 + alpha_isc * delta_t)
+    return isc_t
+
+
+def aph(tc, isc0, alpha_isc, isat1, isat2, vt, rs, rsh):
+    isc = isc_t(tc, isc0, alpha_isc)
+    vd_sc, _ = fvd(0, isc, rs)
+    isat1_t = isat_t(tc, isat0=isat1)
+    id1_sc, _ = fid(isat1_t, vd_sc, 1.0, vt)
+    id2_sc, _ = fid(isat2, vd_sc, 2.0, vt)
+    ish_sc, _ = fish(vd_sc, rsh)
+    aph = 1.0 + (id1_sc + id2_sc + ish_sc) / isc
+    return aph