ENH: porting pvsyst and single diode param est

from PVLIB_MATLAB by Matt Guttenberg Aug-2016, see pvlib#229 and pvlib#227

Added ported functions and tests to pvlib

Added ported versions of calc_theta_phi_exact, lambertw, singlediode,
update_io_known_n, update_rsh_fixed_pt and associated tests to pvl-python.
Also added a __init__ file ot the tests folder to describe the file folder
and allow the tests to be called if necessary. PVsyst_parameter_estimation
and singlediode are not complete versions

Changed the tests to properly call the appropriate functions

Finished the singlediode script

Added Test script for singlediode

Added tests for singlediode functions and updated scripts

Updated the lambertw and singlediode scripts to fix some small errors that
were found while testing. Added 38 tests for the singlediode script
functions althought more tests will be added.

Finished singlediode and corresponding tests

Finished adding tests for singlediode that tests all of the
functionalities of the script, including properly raising errors. Updated
the singlediode script to fix any errors found while testing. Added a
statement to the lambertw function to make the calculations slightly
faster in the event a -inf case.

Added a new file for the Schumaker_QSpline script

Finished Schumaker_QSpline script

Fixed Schumaker_QSpline script based on test

Adjusted all index references to be integer values to avoid the warnings
that were displayed. Fixed all logic errors. Fixed the sorting mechanism
for the final array.

Added test file for Schumaker_QSpline script

Added a Schumaker_QSpline test

Added a Schumaker_QSpline test (more to be added) and fixed the
singlediode_tests script to adhere to PEP8 standards

Added a couple more tests for Schumaker_QSpline

Added the est_single_diode_parameter port and fixed typos

Added the est_single_diode_parameter ported script to pvlib. Fixed a typo
in the Schumaker_QSpline script that caused an error under a very
particular circumstance. Added brackets in the Schumaker_QSpline test
script to get rid of syntax errors.

Added tests for the est_single_diode_param script

Delete __init__.py

Changed one test case for lambertw

Fixed some ambiguities in the imports

Updated the imports in many of the regular and test scripts to fix some
ambiguity issues that might have been occurring during testing

Started adding code to the PVsyst_parameter_estimation script

Added Classes to singlediode and PVsyst_Parameter_estimation

Added a class to move on with the code for the PVsyst_paramter_estimation
script and added a class to the singlediode script to wrap up the answers
in a more succinct format. Adjusted the singlediode tests to reflect this
change.

Ported over some of the code for PVsyst_parameter_estimation

Ported over some of the sections from PVsyst_parameter_estimation. Added
numdiff function to the PVsyst_parameter_estimation script. Since the
script uses matplotlib, added the requirement to the setup.

Finished PVsyst_parameter_estimation and adjusted outputs of functions

Finished porting PVsyst_parameter_estimation except for the section with
the robust_fit. Changed the outputs of this function and singlediode from
a class to a ordereddict and changed the inputs for
PVsyst_parameter_estimation from classes to ordereddicts as well.

Changed the setup to include scipy as a required toolbox

Added robust fit algorithm to PVsyst_parameter_estimation

Added the robust fit algorithm where it was needed to determine desired
information. Since this alrogithm used the statsmodels toolbox, the setup
was changed to reflect this addition.

Removed lambertw since scipy already had a working copy. Script Changes

Removed lambertw and associated test script since scipy.special already
had a working lambertw that has been tested and verified. Changed
associated scripts to reflect this change. Some bug fixes for
Pvsyst_paramter_estimation

Removed some coded scripts and updated PVsyst_parameter_estimation

Removed i_from_v, v_from_i and singlediode and used the already coded and
tested versions within pvlib. Updated PVsyst_parameter_estimation based on
some bug fixes.

More bug fixes

Fixed additional typos and other errors that have prevented the code from
completing. Now that the code runs to completion, additional checks will
be run to determine whether the converged values are appropriate.

Updated calc_theta_phi_exact to handle an edge case

modified calc_theta_phi_exact to handle the edge case where nnsvth was
equal to 0.

Brought back lambertw script for use in update_rsh_fixed_pt

Reintroduced lambertw script for use in update_rsh_fixed_pt since the
lambertw script was having problems for cases where rsh was extremely high
in calc_theta_phi_exact

A couple mroe PVsyst_parameter_estimation bug fixes

Fixed import errors and removed unused lines

Fixed an import error in update_io_known_n and removed v_from_i from the
tests since the code is using the code that already exists and has been
tested in pvsystem. Removed an unused import in update_rsh_fixed_pt

Fixed typos / errors

Fixed an error in PVsyst_parameter_estimation in the numdiff function
where the summations were not occuring properly. Fixed a typo in pvsystem
where a parameter was accidentally changed.

Final bug fixes for PVsyst_parameter_estimation

Fixed a typo in one of the plots

Style Update

Modified all written scripts to adhere to the flake8 and pep8-naming
guidelines. Updated tests as suggested by Will

Documentation Update

Updated script documentation for each of the written scripts to fit better
with the pvlib documentation

Added Documents from Cliff

Author: Matt Guttenberg

.gitignore

@@ -86,3 +86,4 @@ coverage.xml
 #Ignore some notebooks
 *.ipynb
 !docs/tutorials/*.ipynb
+*.idea

pvlib/PVsyst_parameter_estimation.py

@@ -0,0 +1,238 @@
+import numpy as np
+
+
+def schumaker_qspline(x, y):
+    """
+    Schumaker_QSpline fits a quadratic spline which preserves monotonicity and
+    convexity in the data.
+
+    Syntax
+    ------
+    outa, outxk, outy, kflag = schumaker_qspline(x, y)
+
+    Description
+    -----------
+    Calculates coefficients for C1 quadratic spline interpolating data X, Y
+    where length(x) = N and length(y) = N, which preserves monotonicity and
+    convexity in the data.
+
+    Parameters
+    ----------
+    x, y: numpy arrays of length N containing (x, y) points between which the
+          spline will interpolate.
+
+    Returns
+    -------
+    outa: a Nx3 matrix of coefficients where the ith row defines the quadratic
+          interpolant between xk_i to xk_(i+1), i.e., y = A[i, 0] *
+          (x - xk[i]] ** 2 + A[i, 1] * (x - xk[i]) + A[i, 2]
+    outxk: an ordered vector of knots, i.e., values xk_i where the spline
+           changes coefficients. All values in x are used as knots. However
+           the algorithm may insert additional knots between data points in x
+           where changes in convexity are indicated by the (numerical)
+           derivative. Consequently output outxk has length >= length(x).
+    outy: y values corresponding to the knots in outxk. Contains the original
+          data points, y, and also y-values estimated from the spline at the
+          inserted knots.
+    kflag: a vector of length(outxk) of logicals, which are set to true for
+           elements of outxk that are knots inserted by the algorithm.
+
+    References
+    ----------
+    [1] PVLib MATLAB
+    [2] L. L. Schumaker, "On Shape Preserving Quadratic Spline Interpolation",
+        SIAM Journal on Numerical Analysis 20(4), August 1983, pp 854 - 864
+    [3] M. H. Lam, "Monotone and Convex Quadratic Spline Interpolation",
+        Virginia Journal of Science 41(1), Spring 1990
+    """
+
+    # A small number used to decide when a slope is equivalent to zero
+    eps = 1e-6
+
+    # Make sure vectors are 1D arrays
+    if x.ndim != 1.:
+        x = x.flatten([range(x.size)])
+    if y.ndim != 1.:
+        y = y.flatten([range(y.size)])
+
+    n = len(x)
+
+    # compute various values used by the algorithm: differences, length of line
+    # segments between data points, and ratios of differences.
+    delx = np.diff(x)  # delx[i] = x[i + 1] - x[i]
+    dely = np.diff(y)
+
+    delta = dely / delx
+
+    # Calculate first derivative at each x value per [3]
+
+    s = np.zeros(x.shape)
+
+    left = np.append(0., delta)
+    right = np.append(delta, 0.)
+
+    pdelta = left * right
+
+    u = pdelta > 0.
+
+    # [3], Eq. 9 for interior points
+    # fix tuning parameters in [2], Eq 9 at chi = .5 and eta = .5
+    s[u] = pdelta[u] / (.5 * left[u] + .5 * right[u])
+
+    # [3], Eq. 7 for left endpoint
+    if delta[0] * (2. * delta[0] - s[1]) > 0.:
+        s[0] = 2. * delta[0] - s[1]
+
+    # [3], Eq. 8 for right endpoint
+    if delta[n - 2] * (2. * delta[n - 2] - s[n - 2]) > 0.:
+        s[n - 1] = 2. * delta[n - 2] - s[n - 2]
+
+    # determine knots. Start with initial pointsx
+    # [2], Algorithm 4.1 first 'if' condition of step 5 defines intervals
+    # which won't get internal knots
+    tests = s[0.:(n - 1)] + s[1:n]
+    u = np.abs(tests - 2. * delta[0:(n - 1)]) <= eps
+    # u = true for an interval which will not get an internal knot
+
+    k = n + sum(~u)  # total number of knots = original data + inserted knots
+
+    # set up output arrays
+    # knot locations, first n - 1 and very last (n + k) are original data
+    xk = np.zeros(k)
+    yk = np.zeros(k)  # function values at knot locations
+    # logicals that will indicate where additional knots are inserted
+    flag = np.zeros(k, dtype=bool)
+    a = np.zeros((k, 3.))
+
+    # structures needed to compute coefficients, have to be maintained in
+    # association with each knot
+
+    tmpx = x[0:(n - 1)]
+    tmpy = y[0:(n - 1)]
+    tmpx2 = x[1:n]
+    tmps = s[0.:(n - 1)]
+    tmps2 = s[1:n]
+    diffs = np.diff(s)
+
+    # structure to contain information associated with each knot, used to
+    # calculate coefficients
+    uu = np.zeros((k, 6.))
+
+    uu[0:(n - 1), :] = np.array([tmpx, tmpx2, tmpy, tmps, tmps2, delta]).T
+
+    # [2], Algorithm 4.1 subpart 1 of Step 5
+    # original x values that are left points of intervals without internal
+    # knots
+    xk[u] = tmpx[u]
+    yk[u] = tmpy[u]
+    # constant term for each polynomial for intervals without knots
+    a[u, 2] = tmpy[u]
+    a[u, 1] = s[u]
+    a[u, 0] = .5 * diffs[u] / delx[u]  # leading coefficients
+
+    # [2], Algorithm 4.1 subpart 2 of Step 5
+    # original x values that are left points of intervals with internal knots
+    xk[~u] = tmpx[~u]
+    yk[~u] = tmpy[~u]
+
+    aa = s[0:(n - 1)] - delta[0:(n - 1)]
+    b = s[1:n] - delta[0:(n - 1)]
+
+    sbar = np.zeros(k)
+    eta = np.zeros(k)
+    # will contain mapping from the left points of intervals containing an
+    # added knot to each inverval's internal knot value
+    xi = np.zeros(k)
+
+    t0 = aa * b >= 0
+    # first 'else' in Algorithm 4.1 Step 5
+    v = np.logical_and(~u, t0[0:len(u)])
+    q = np.sum(v)  # number of this type of knot to add
+
+    if q > 0.:
+        xk[(n - 1):(n + q - 1)] = .5 * (tmpx[v] + tmpx2[v])  # knot location
+        uu[(n - 1):(n + q - 1), :] = np.array([tmpx[v], tmpx2[v], tmpy[v],
+                                               tmps[v], tmps2[v], delta[v]]).T
+        xi[v] = xk[(n - 1):(n + q - 1)]
+
+    t1 = np.abs(aa) > np.abs(b)
+    w = np.logical_and(~u, ~v)  # second 'else' in Algorithm 4.1 Step 5
+    w = np.logical_and(w, t1)
+    r = np.sum(w)
+
+    if r > 0.:
+        xk[(n + q - 1):(n + q + r - 1)] = tmpx2[w] + aa[w] * delx[w] / diffs[w]
+        uu[(n + q - 1):(n + q + r - 1), :] = np.array([tmpx[w], tmpx2[w],
+                                                       tmpy[w], tmps[w],
+                                                       tmps2[w], delta[w]]).T
+        xi[w] = xk[(n + q - 1):(n + q + r - 1)]
+
+    z = np.logical_and(~u, ~v)  # last 'else' in Algorithm 4.1 Step 5
+    z = np.logical_and(z, ~w)
+    ss = np.sum(z)
+
+    if ss > 0.:
+        xk[(n + q + r - 1):(n + q + r + ss - 1)] = tmpx[z] + b[z] * delx[z] / \
+                                                             diffs[z]
+        uu[(n + q + r - 1):(n + q + r + ss - 1), :] = \
+            np.array([tmpx[z], tmpx2[z], tmpy[z], tmps[z], tmps2[z],
+                      delta[z]]).T
+        xi[z] = xk[(n + q + r - 1):(n + q + r + ss - 1)]
+
+    # define polynomial coefficients for intervals with added knots
+    ff = ~u
+    sbar[ff] = (2 * uu[ff, 5] - uu[ff, 4]) + \
+               (uu[ff, 4] - uu[ff, 3]) * (xi[ff] - uu[ff, 0]) / (uu[ff, 1] -
+                                                                 uu[ff, 0])
+    eta[ff] = (sbar[ff] - uu[ff, 3]) / (xi[ff] - uu[ff, 0])
+
+    sbar[(n - 1):(n + q + r + ss - 1)] = \
+        (2 * uu[(n - 1):(n + q + r + ss - 1), 5] -
+         uu[(n - 1):(n + q + r + ss - 1), 4]) + \
+        (uu[(n - 1):(n + q + r + ss - 1), 4] -
+         uu[(n - 1):(n + q + r + ss - 1), 3]) * \
+        (xk[(n - 1):(n + q + r + ss - 1)] -
+         uu[(n - 1):(n + q + r + ss - 1), 0]) / \
+        (uu[(n - 1):(n + q + r + ss - 1), 1] -
+         uu[(n - 1):(n + q + r + ss - 1), 0])
+    eta[(n - 1):(n + q + r + ss - 1)] = \
+        (sbar[(n - 1):(n + q + r + ss - 1)] -
+         uu[(n - 1):(n + q + r + ss - 1), 3]) / \
+        (xk[(n - 1):(n + q + r + ss - 1)] -
+         uu[(n - 1):(n + q + r + ss - 1), 0])
+
+    # constant term for polynomial for intervals with internal knots
+    a[~u, 2] = uu[~u, 2]
+    a[~u, 1] = uu[~u, 3]
+    a[~u, 0] = .5 * eta[~u]  # leading coefficient
+
+    a[(n - 1):(n + q + r + ss - 1), 2] = \
+        uu[(n - 1):(n + q + r + ss - 1), 2] + \
+        uu[(n - 1):(n + q + r + ss - 1), 3] * \
+        (xk[(n - 1):(n + q + r + ss - 1)] -
+         uu[(n - 1):(n + q + r + ss - 1), 0]) + \
+        .5 * eta[(n - 1):(n + q + r + ss - 1)] * \
+        (xk[(n - 1):(n + q + r + ss - 1)] -
+         uu[(n - 1):(n + q + r + ss - 1), 0]) ** 2.
+    a[(n - 1):(n + q + r + ss - 1), 1] = sbar[(n - 1):(n + q + r + ss - 1)]
+    a[(n - 1):(n + q + r + ss - 1), 0] = \
+        .5 * (uu[(n - 1):(n + q + r + ss - 1), 4] -
+              sbar[(n - 1):(n + q + r + ss - 1)]) / \
+        (uu[(n - 1):(n + q + r + ss - 1), 1] -
+         uu[(n - 1):(n + q + r + ss - 1), 0])
+
+    yk[(n - 1):(n + q + r + ss - 1)] = a[(n - 1):(n + q + r + ss - 1), 2]
+
+    xk[n + q + r + ss - 1] = x[n - 1]
+    yk[n + q + r + ss - 1] = y[n - 1]
+    flag[(n - 1):(n + q + r + ss - 1)] = True  # these are all inserted knots
+
+    tmp = np.vstack((xk, a.T, yk, flag)).T
+    # sort output in terms of increasing x (original plus added knots)
+    tmp2 = tmp[tmp[:, 0].argsort(kind='mergesort')]
+    outxk = tmp2[:, 0]
+    outn = len(outxk)
+    outa = tmp2[0:(outn - 1), 1:4]
+    outy = tmp2[:, 4]
+    kflag = tmp2[:, 5]
+    return outa, outxk, outy, kflag

pvlib/Schumaker_QSpline.py

@@ -0,0 +1,141 @@
+import numpy as np
+from pvlib.lambertw import lambertw
+
+
+def calc_theta_phi_exact(imp, il, vmp, io, nnsvth, rs, rsh):
+    """
+    CALC_THETA_PHI_EXACT computes Lambert W values appearing in the analytic
+    solutions to the single diode equation for the max power point.
+
+    Syntax
+    ------
+    theta, phi = calc_theta_phi_exact(imp, il, vmp, io, nnsvth, rs, rsh)
+
+    Description
+    -----------
+    calc_theta_phi_exact calculates values for the Lambert W function which
+    are used in the analytic solutions for the single diode equation at the
+    maximum power point. For V=V(I),
+    phi = W(Io*Rsh/n*Vth * exp((IL + Io - Imp)*Rsh/n*Vth)). For I=I(V),
+    theta = W(Rs*Io/n*Vth *
+    Rsh/ (Rsh+Rs) * exp(Rsh/ (Rsh+Rs)*((Rs(IL+Io) + V)/n*Vth))
+
+    Parameters
+    ----------
+    imp: a numpy array of length N of values for Imp (A)
+    il: a numpy array of length N of values for the light current IL (A)
+    vmp: a numpy array of length N of values for Vmp (V)
+    io: a numpy array of length N of values for Io (A)
+    nnsvth: a numpy array of length N of values for the diode factor x
+            thermal voltage for the module, equal to Ns
+            (number of cells in series) x Vth
+            (thermal voltage per cell).
+    rs: a numpy array of length N of values for the series resistance (ohm)
+    rsh: a numpy array of length N of values for the shunt resistance (ohm)
+
+    Returns
+    -------
+    theta: a numpy array of values for the Lamber W function for solving
+           I = I(V)
+    phi: a numpy array of values for the Lambert W function for solving
+         V = V(I)
+
+    References
+    ----------
+    [1] PVLib MATLAB
+    [2] C. Hansen, Parameter Estimation for Single Diode Models of Photovoltaic
+        Modules, Sandia National Laboratories Report SAND2015-XXXX
+    [3] A. Jain, A. Kapoor, "Exact analytical solutions of the parameters of
+        real solar cells using Lambert W-function", Solar Energy Materials and
+        Solar Cells, 81 (2004) 269-277.
+    """
+
+    # Argument for Lambert W function involved in V = V(I) [2] Eq. 12; [3]
+    # Eq. 3
+    if any(nnsvth == 0.):
+        argw = []
+        for i in nnsvth:
+            if i == 0.:
+                argw.append(float("Inf"))
+            else:
+                argw.append(rsh * io / i * np.exp(rsh * (il + io - imp) / i))
+        argw = np.array(argw)
+    else:
+        argw = rsh * io / nnsvth * np.exp(rsh * (il + io - imp) / nnsvth)
+    u = argw > 0
+    w = np.zeros(len(u))
+    w[~u] = float("Nan")
+    if any(argw[u] == float("Inf")):
+        tmp = []
+        for i in argw[u]:
+            if i == float("Inf"):
+                tmp.append(float("Nan"))
+            else:
+                tmp.append(lambertw(i).real)
+        tmp = np.array(tmp)
+    else:
+        tmp = lambertw(argw[u]).real
+    ff = np.isnan(tmp)
+
+    # NaN where argw overflows. Switch to log space to evaluate
+    if any(ff):
+        logargw = np.log(rsh[u]) + np.log(io[u]) - np.log(nnsvth[u]) + \
+                  rsh[u] * (il[u] + io[u] - imp[u]) / nnsvth[u]
+        # Three iterations of Newton-Raphson method to solve w+log(w)=logargW.
+        # The initial guess is w=logargW. Where direct evaluation (above)
+        # results in NaN from overflow, 3 iterations of Newton's method gives
+        # approximately 8 digits of precision.
+        x = logargw
+        for i in range(3):
+            x *= ((1. - np.log(x) + logargw) / (1. + x))
+        tmp[ff] = x[ff]
+    w[u] = tmp
+    phi = np.transpose(w)
+
+    # Argument for Lambert W function involved in I = I(V) [2] Eq. 11; [3]
+    # E1. 2
+    if any(nnsvth == 0.):
+        argw = []
+        for i in nnsvth:
+            if i == 0.:
+                argw.append(float("Inf"))
+            else:
+                argw.append(rsh / (rsh + rs) * rs * io / i *
+                            np.exp(rsh / (rsh + rs) * (rs * (il + io) + vmp) /
+                                   i))
+        argw = np.array(argw)
+    else:
+        argw = rsh / (rsh + rs) * rs * io / nnsvth * \
+               np.exp(rsh / (rsh + rs) * (rs * (il + io) + vmp) / nnsvth)
+    u = argw > 0
+    w = np.zeros(len(u))
+    w[~u] = float("Nan")
+    if any(argw[u] == float("Inf")):
+        tmp = []
+        for i in argw[u]:
+            if i == float("Inf"):
+                tmp.append(float("Nan"))
+            else:
+                tmp.append(lambertw(i).real)
+        tmp = np.array(tmp)
+    else:
+        tmp = lambertw(argw[u]).real
+    ff = np.isnan(tmp)
+
+    # NaN where argw overflows. Switch to log space to evaluate
+    if any(ff):
+        logargw = np.log(rsh[u]) / (rsh[u] + rs[u]) + np.log(rs[u]) + \
+                  np.log(io[u]) - np.log(nnsvth[u]) + \
+                  (rsh[u] / (rsh[u] + rs[u])) * \
+                  (rs[u] * (il[u] + io[u]) + vmp[u]) / nnsvth[u]
+        # Three iterations of Newton-Raphson method to solve w+log(w)=logargW.
+        # The initial guess is w=logargW. Where direct evaluation (above)
+        # results in NaN from overflow, 3 iterations of Newton's method gives
+        # approximately 8 digits of precision.
+        x = logargw
+        for i in range(3):
+            x *= ((1. - np.log(x) + logargw) / (1. + x))
+        tmp[ff] = x[ff]
+    w[u] = tmp
+    theta = np.transpose(w)
+    return theta, phi