Pull Request for PVsyst_parameter_estimation

[mattguttenberg]

This is the pull request for #227. Since the issue was opened, I have finished the singlediode port with accompanying tests and the Schumaker_QSpline port. I still need to write some tests to make sure that Schumaker_QSpline is fully functional; so far I have only tested it with out nose by using the python command console. Once the tests for Schumaker_QSpline are finished, I'll begin porting est_single_diode_param and write some tests for that code and then finish up with the PVsyst_parameter_estimation port with some accompanying tests. With that said, I wanted to start this pull request since this is such a large merge and wanted to start the feedback / discussion process.

[wholmgren]

Thanks @mattguttenberg! This looks great. I really appreciate that you wrote tests, too. Let's first address a few bigger and more general things. After that, I'll go through the code more carefully and make some line by line comments.

style: please install flake8 and pep8-naming and run them against all of your new code. This way I don't have to comment on every little style issue. Also run them against some of the existing modules so that you can get a feel for what to ignore. The big things that I see right now are line length and capitalization.

module structure: pvlib python prefers to consolidate similar functions into modules. We'll also want to make some (all?) of the helper functions into private methods. This is so that they don't become automatically documented by sphinx, and so that people don't rely on them for their own code.

lambertw: Is there a reason that you implemented lambertw instead of using the scipy lambertw? See e.g. pvsystem.v_from_i for an example of its usage in pvlib.

schumaker_qspline: I'm not familiar with this algorithm, but I'm wondering if scipy and/or statsmodels already has it (maybe under a different name or with the right parameters).

plotting: I'd like to see a better separation between plotting code and the rest of the code. Move the plotting stuff into a separate functions. You should also import matplotlib within those functions rather than at the top of the module. This way we don't need to make matplotlib a required dependency.

tests: put the calculations inside the test functions. numpy now recommends using assert_allclose instead of assert_array_almost_equal

documentation: Compare your documentation syntax to the numpydoc standard and/or the rest of pvlib.

[cwhanse]

Ø schumaker_qspline: I'm not familiar with this algorithm, but I'm wondering if scipy and/or statsmodels already has it (maybe under a different name or with the right parameters).
I very much doubt it. The algorithm produces a quadratic spline that preserves both monotonicity and concavity. I didn’t find any code when I looked 2 years ago. There’s a nascent R package now which might give us some ideas about making my matlab code more general and perhaps become a submission to a python package.

Cliff

[mattguttenberg]

@wholmgren Thank you for all of the great feedback. I did notice that lambertw was already added to scipy and had made use of it in most of the scripts I had wrote. When I was running tests for calc_theta_phi_exact, however, I found that in some of the test cases I written, an edge case was observed were lambertw from scipy returned infinity while the lambertw script produced a very large value but not infinity which lead to an assertion error. Thus I had chosen to keep the function in on this port instead of 100% switching to the scipy version of the function to try and minimize discrepancy between the MATLAB version and ported version I had wrote.

[wholmgren]

Can you use pvsystem.v_from_i instead?

[mattguttenberg]

I thought I removed the v_from_i script I wrote and changed to using the pvsystem function in one of my latest commits.

[mattguttenberg]

After double checking with my latest commit, I had already removed the v_from_i script I had written after I realized that pvsystem already had it and changed the scripts to use the one that already existed

[wholmgren]

It seemed to me that you could avoid the whole lambertw issue if you relied on pvsystem.v_from_i (which already handles overflow internally). Maybe I didn't read carefully enough, though.

[mattguttenberg]

Oh I see. The error occurred in calc_theta_phi_exact. v_from_i could be used to calculated phi from the looks of it since they share the same argw formula but the formula for theta would still need the script as its argw formula is different from that of v_from_i and if I remember correctly, it was the theta term that was throwing the error mentioned above.

[wholmgren]

Could we add similar code to pvsystem.i_from_v? If so, we might want to move the log solver into its own function that lives in, say, tools.py. Could be something like...

# in tools.py
def lambertw_maybe_log(argW, logargW):
    # scipy's result
    lambertwterm = lambertw(argW).real

    # 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.
    w = logargW
    for i in range(0, 3):
        w = w * (1 - np.log(w) + logargW) / (1 + w)
    lambertwterm_log = w

    lambertwterm = np.where(np.isfinite(lambertwterm), lambertwterm,
                            lambertwterm_log)

    return lambertwterm

def v_from_i(...):
    # ...
    argW = I0 * Rsh / nNsVth * np.exp(Rsh * (-I + IL + I0) / nNsVth)

    # Calculate using log(argW) in case argW is really big
    logargW = (np.log(I0) + np.log(Rsh) - np.log(nNsVth) +
               Rsh * (-I + IL + I0) / nNsVth)

    lambertwterm = tools.lambertw_maybe_log(argW, logargW)
    # ...

def i_from_v(...):
    # ...
    argW = (Rs*I0*Rsh *
            np.exp(Rsh*(Rs*(IL+I0)+V) / (nNsVth*(Rs+Rsh))) /
            (nNsVth*(Rs + Rsh)))

    logargW = np.log(...) + np.log(...) + ...

    lambertwterm = tools.lambertw_maybe_log(argW, logargW)
    # ...

[cwhanse]

Here's the IV curve data we provide with PVLib matlab, in csv format. Also, python code to read the data into a dict with the structure currently used by the parameter estimation code.
PVsyst_demo.txt

import numpy as np

specs = {key: None for key in ['ns', 'aisc','bvoc','type']}
keylist  = ['isc', 'imp', 'vmp', 'voc', 'poa', 'tc', 'ee']
ivcurves = {key: None for key in keylist}

with open('PVsyst_demo.txt', 'r') as f:

    Ns, aIsc, bVoc, descr = f.readline().split(',')

    specs['ns'] = int(float(Ns))
    specs['aisc'] = float(aIsc)
    specs['bvoc'] = float(bVoc)
    specs['type'] = descr

    strN, strM = f.readline().split(',')
    N = int(strN)
    M = int(strM)

    isc = np.empty(N)
    imp = np.empty_like(isc)
    vmp = np.empty_like(isc)
    voc = np.empty_like(isc)
    ee = np.empty_like(isc)
    poa = np.empty_like(isc)
    tc = np.empty_like(isc)
    v = np.empty((N, M))    
    i = np.empty_like(v)    

    for k in range(0,N):
        isc[k], imp[k], vmp[k], voc[k], poa[k], tc[k], ee[k] = f.readline().split(',')
        tmp = [str.strip(',') for str in f.readline().split()]
        while len(tmp)<M:
            tmp.append('NaN')
        v[k,] = [float(x) for x in tmp]
        tmp = [str.strip(',') for str in f.readline().split()]
        while len(tmp)<M:
            tmp.append('NaN')
        i[k,] = [float(x) for x in tmp]

ivcurves['v'] = v
ivcurves['i'] = i

[cwhanse]

Here's the solution for pvsyst from Matlab, and python code to read the csv file into an OrderedDict.

PVsyst_demo_model.txt

import numpy as np
from collections import OrderedDict


specs = {key: None for key in ['ns', 'aisc','bvoc','type']}
paramlist  = ['IL_ref', 'Io_ref', 'eG', 'Rsh_ref', 'Rsh0', 'Rshexp', 'Rs_ref', 'gamma_ref', 'mugamma']
varlist = ['Iph', 'Io', 'Rsh', 'Rs', 'u']
pvsyst = OrderedDict(key = paramlist + varlist)

with open('PVsyst_demo_model.txt', 'r') as f:

    Ns, aIsc, bVoc, descr = f.readline().split(',')

    specs['ns'] = int(float(Ns))
    specs['aisc'] = float(aIsc)
    specs['bvoc'] = float(bVoc)
    specs['type'] = descr

    IL_ref, Io_ref, eG, Rsh_ref, Rsh0, Rshexp, Rs_ref, gamma_ref, mugamma = f.readline().split(',')
    for key in paramlist:
        pvsyst[key] = float(eval(key))

    strN = f.readline()
    N = int(strN)

    Iph = np.empty(N)
    Io = np.empty_like(Iph)
    Rsh = np.empty_like(Iph)
    Rs = np.empty_like(Iph)
    u = np.empty_like(Iph)

    for k in range(0,N):
        Iph[k], Io[k], Rsh[k], Rs[k], u[k] = f.readline().split(',')

for var in varlist:
    pvsyst[var] = eval(var)

[wholmgren]

@mattguttenberg @cwhanse any plans for finishing this PR? Maybe someone else wants to pick it up? I don't personally have a need for this feature or the time to do it.

[mikofski]

I will try to pick this up and fix it soon, maybe by the end of the year?

[mikofski]

continuing this in a new branch, also PV_LIB MATLAB works fine in Octave on Linux.

[mikofski]

@wholmgren and @cwhanse do you want to close this in favor of #708

• I've squashed all of these commits into c00d419 and rebased them on v0.4 then rebased that on master, to preserve the history
• then isolated them in a new temporary folder until we decide what to do with them
• then got them working in name only, there are some tests, all tests pass, but I think the testing is still insufficient