Test singlediode() implementations using a "gold" function and data set

[markcampanelli]

#409 Suggests Bishop [1988] as a truly explicit method to calculate i-from-v. Bishop's method is not subject to Newton convergence failure or to LambertW argument "blow up", handles ideal device parameters without extra logic, has no dependency on scipy, and naturally vectorizes via standard numpy functions' vectoization. The trade-off is that the explicit device terminal current computation occurs at the diode junction voltage (V_d := V + I R_s), not the device terminal voltage. (The device terminal voltage is "back-calculated" as V = V_d - I Rs.) Thus, additional computation work must be done to determine currents at specified terminal voltages, inc. I_sc @ V_sc=0, I_mp @ V_mp, and I_oc=0 @ V_oc.

The Bishop method is a very good candidate for a "gold" function implementation for testing the accuracy and speed of various implementations of i_from_v() that pvlib may wish to offer users or refactor to. By inverting the dataset, it could also be used to test v_from_i() implementations.

singlediode(photocurrent, saturation_current, resistance_series, resistance_shunt, nNsVth, ivcurve_pnts=None)

I think that after this bishop88 function is implemented (presumably by #409), we should create and commit a corresponding "gold" dataset that we can use for better testing, refactoring, and optimization of the singlediode() function. I think we should try to cover a large number (all?) of the devices in the CEC database over a combination of temperature and irradiance ranges (perhaps the IEC 61853-1 matrix?). This dataset should also be valuable to the greater PV modeling community who wish to test/contribute new computational algorithms.

Outstanding questions:
(1) What should the voltage range be? Specifically, the diode voltages will give some terminal voltages less than zero, and how much beyond V_oc should we try to extend?
(2) Can/should we also include I_sc, I_x, I_mp, I_xx, and V_oc points, and what about the slope at I_sc and slope at V_oc?
(3) Does python have an interval analysis package up to the task of verifying our computations, e.g., intPy?
(4) Is it feasible to run such a comprehensive test suite on every PR?

cc @mikofski @adriesse

[wholmgren]

I am a big +1 on testing pvlib against "gold datasets" with every PR. The implementation probably depends on the size and scope of a dataset. Adding numeric checks to the test suite to guard against unknown changes was a big effort, so I'd recommend starting small.

[adriesse]

I don't think singlediode() is the right hierarchy level for verifying algorithms and implementations. I would prefer to see functions like Voc(), Isc(), and especially MPP() in addition to I_from_V() and V_from_I(). For all of these there can be reference implementations and others. singlediode() just seems like a convenient way to call several of these in sequence.

[cwhanse]

I agree with Anton on the factorization.

…

Sent from my iPhone On Jan 31, 2018, at 1:12 PM, Anton Driesse <notifications@github.com<mailto:notifications@github.com>> wrote: I don't think singlediode() is the right hierarchy level for verifying algorithms and implementations. I would prefer to see functions like Voc(), Isc(), and especially MPP() in addition to I_from_V() and V_from_I(). For all of these there can be reference implementations and others. singlediode() just seems like a convenient way to call several of these in sequence, while most of the time they are non needed. - You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub<#411 (comment)>, or mute the thread<https://github.com/notifications/unsubscribe-auth/AFJNL38NAFtWN5CooIs98C8-OBSJyLMSks5tQMk0gaJpZM4RwzC5>.

[markcampanelli]

I’m getting started on this. I’m expecting a few PRs out of it as I start out small.

[markcampanelli]

I have decided to start with creating minimal "gold" I-V dataset for perhaps a half dozen or so model parameter sets.

Bear with me as this is my first real foray into reliable computing. Using pyinterval, below is a notebook snippet where I generate I-V data pairs using the bishop88 technique for a SunPower module, then I calculate the residuals using enclosing intervals for all the computed pairs and model parameters. This produces list of residual intervals, all of which contain zero and which have a width that gives one more/less confidence as to the validity of the solution pair. So, I think asserting zero in the interval and the interval endpoints not too far from zero could be the acceptance criterion for creating the gold dataset. Using the dataset would then entail running a particular choice of i_from_v() or v_from_i() and verifying that the answer is not too far from the gold values (no interval analysis used here, keeping this fairly straightforward). I am still wondering about the spacing and breadth of the I-V pairs.

Feedback here is welcomed.

v_oc_V = pvsystem.v_oc(rsh, rs, nnsvt, i_rs_A, i_ph_A)
v_d_V = np.linspace(0., v_oc_V, 11)
i_A = i_ph_A - i_rs_A*np.expm1(v_d_V/nnsvt) - g_p_S*v_d_V
v_V = v_d_V - r_s_Ohm*i_A
for sol_pair in zip(v_V, i_A):
    print("{:+20.16f} V, {:+20.16f} A".format(sol_pair[0], sol_pair[1]))
>>>
 -2.1696985296000002 V,  +5.8640500800000002 A
 +3.6145797904049686 V,  +5.8468400769024944 A
 +9.3988583628994498 V,  +5.8296293914008999 A
+15.1831387169629082 V,  +5.8124138908480232 A
+20.9674316417980293 V,  +5.7951644152365844 A
+26.7518132661460761 V,  +5.7776752112118412 A
+32.5368207553018678 V,  +5.7584944806796940 A
+38.3262443539114628 V,  +5.7273783191913070 A
+44.1468280728327542 V,  +5.6120456163199428 A
+50.1872779676257110 V,  +4.9024800056873206 A
+57.7791061885889121 V,  +0.0000000000000025 A

interval_res_A_list = [ interval(i_ph_A) - interval(i_rs_A)*\
    imath.expm1((interval(sol_pair[0]) + interval(sol_pair[1])*interval(rs))/interval(nnsvt)) - \
    interval(g_p_S)*(interval(sol_pair[0]) + interval(sol_pair[1])*interval(rs)) - interval(sol_pair[1]) \
    for sol_pair in zip(v_V, i_A) ]
interval_res_A_list
>>>
[interval([-0.0, 1.7763568394002505e-15]),
 interval([-8.881784197001252e-16, 8.881784197001252e-16]),
 interval([-0.0, 1.7763568394002505e-15]),
 interval([-1.7763568394002505e-15, 0.0]),
 interval([-8.881784197001252e-16, 8.881784197001252e-16]),
 interval([-8.881784197001252e-16, 8.881784197001252e-16]),
 interval([-0.0, 1.7763568394002505e-15]),
 interval([-8.881784197001252e-16, 8.881784197001252e-16]),
 interval([-8.881784197001252e-16, 1.7763568394002505e-15]),
 interval([-8.881784197001252e-16, 6.217248937900877e-15]),
 interval([-2.223221606811876e-14, 2.1316282072803006e-14])]

[mikofski]

Couldn't you also use the numerical precision test to generate exact values?

[markcampanelli]

Might be nothing more than a matter of approach, but my numerical background is more in terms of thinking in terms of intervals with exactly represented floating point endpoints containing a true value, rather than arbitrary precision “objects”.

[markcampanelli]

Looks like the built in decimal package is another alternative.

[mikofski]

No, IMO the decimal type is not good for science and engineering.

[mikofski]

Ah, my mistake, you mean for creating the gold data set, yes, decimal will be similar to sympy symbols, except that you need to specify the precision of decimal manually, whereas sympy will adjust that automatically because it auto differentiates

[markcampanelli]

Here's a gist of my beta I-V curve dataset creator: https://gist.github.com/thunderfish24/10bd9e57e0cb1414fc79601f53f16178

Note the tolerance interval on the residuals that I was able to achieve for the two modules considered (with ideal and degraded versions added) is [-1.e-12, 1.e-12]. I-V curves are engineered to fall "slightly" outside the first quadrant. (I have some unexpected import warnings to track down too.)

I really want to add a x-Si cell in the list of devices, then I will move on to saving the dataset in a lossless interchangeable format, as well as writing a generic function tester (possibly with timing info).

I think a second PR will add the "derived" parameters, such as Isc, Ix, R@Isc, Pmp, Ixx, R@Voc, Voc, and FF, but I will try to make the dataset format extensible so that adding things doesn't break existing consumers.

[mikofski]

@thunderfish24, this it's really great!.Thanks so much for your hard work.

[markcampanelli]

I have some first results comparing the existing i_from_v() and v_from_i() computations against the verified gold values computed using Bishop's explicit method. Preliminary results suggest that v precision can be surprising low for some parameter combo's. More details and code for others to run to follow.

[mikofski]

Voltage precision will be very difficult where di/dv is very large, eg approaching v_oc that's just a fact of life and why we don't bother trying hard in quadrant 4. Reverse bias is easier because at least you know it approaches an asymptote.

[cwhanse]

The v_from_i calculation is where the overflow issues occur in the lambertw scheme. There could certainly be a precision loss there. I also suspect that we can't expect the circle of (I, V from Bishop) then (recover V from I) to maintain precision. Suppose we know Vd precisely. The corresponding current will not be known as precisely because di/dv is very large as @mikofski points out. Precision in V will further declines when that current value I is fed back into v_from_i. I guess my point is that the v_from_i method has a precision, which I inferred from the precision of the lambertw calculation (verified to within 8 decimal places) but not the end-to-end algorithm.

[adriesse]

Perhaps we should define the allowable range for each of the parameters of the gold implementation. I think such ranges should extend beyond what the module database contains for commercial modules, for example, the resistances could both be allowed any non-negative value--perhaps even some negative ones. (I have a vague recollection than some fits produce negative resistances.)

[markcampanelli]

@adriesse I have seen my method try to fit negative series resistances for the SDM if I do not constrain to positive values. esp. in certain parts of the irradiance-temperature space. I am trying to see if a comparable double-diode model (DDM) suffers the same deficiency. This is why I think it can be tricky trying to implement a more physically realistic auxiliary equation for series resistance in the SDM, because the model itself may be so wrong that it confounds such efforts.

[mikofski]

I agree with both @adriesse and @thunderfish24 , and it's good practice. So far I haven't transformed any of the variable going into the Newton or bisection methods. Typically I use ustar = u**2 because that will always be positive, then u = sqrt(ustar) it's also positive

[markcampanelli]

K. I squashed some early bugs, and now I am 95% confident in this early result for checking i_from_v and v_from_i for a small area (2 cm x 2 cm) Si cell based on reference parameters in Bishop' paper.

'v_test_current_sum_max_abs_res': 0.009735069744962271 is the maximum absolute residual of the sum of currents at the diode node for a curve where 'G': 200.0, 'T': 50.0, over 25 points 21 points.

'v_test_max_abs_res': 0.27325988859491335 is the maximum absolute residual of the terminal voltage for the same curve where 'G': 200.0, 'T': 50.0 over 25 points. This is pretty big relative to the Voc ~0.46.

I will be diving into this anomaly soon, but I welcome independent verification using the parameters logged below. I really do hope to push some code up in a PR soon too. I think I'm close to a reasonable "starter" API for this that can be extended and implemented in unit tests.

UPDATE: I tried to make the output well-formatted, proper json (except for the out-of-spec "Infinity" entries, but I think the Python json parser can still accept those).

[{
    "name": "Bishop_ref_cell",
    "Technology": "Mono-c-Si",
    "N_s": 1,
    "i_from_v": {
        "current_sum": {
            "G": 1000.0,
            "T": 75.0,
            "i_l": 0.34609999999999996,
            "i_0": 3.0407536094957474e-05,
            "r_s": 26.6,
            "r_sh": Infinity,
            "nNsVth": 0.04500185335353259,
            "v_gold": [-1.5089440039918927, -1.1032699896018565, -0.7960826132220626, -0.5468527286276292, -0.3361162030859558, -0.1529423682914136, 0.009451536813227968, 0.1555783946319222, 0.2886000367572715, 0.4108218614392964, 0.5239790581185568, 0.6294119193683566, 0.7281782611349393, 0.8211282972540166, 0.9089561213900996, 0.992236079357277, 1.0714490760167272, 1.1470019953384938, 1.2192422969390602, 1.288469163776995, 1.3549421382270659],
            "i_gold": [0.0721330499225567, 0.05697325365344513, 0.045490738268200104, 0.036172832511647435, 0.028292879186989417, 0.021442719744583227, 0.015369040522771882, 0.00990327777849409, 0.00492731831036286, 0.00035503323230140493, -0.003878409236460545, -0.007823090311245784, -0.011518539033739894, -0.014996528720618463, -0.01828299987350146, -0.021399418816702875, -0.024363760414717495, -0.027191233566722706, -0.029894826532529506, -0.03248572343041983, -0.03497362691119921],
            "i_test": [0.07213304992255665, 0.05697325365344513, 0.045490738268200104, 0.03617283251164738, 0.028292879186989417, 0.021442719744583116, 0.015369040522771826, 0.009903277778494035, 0.004927318310362805, 0.0003550332323012939, -0.003878409236460656, -0.00782309031124584, -0.011518539033739894, -0.014996528720618518, -0.018282999873501515, -0.021399418816702875, -0.024363760414717495, -0.02719123356672276, -0.029894826532529506, -0.03248572343041983, -0.03497362691119926],
            "residual": [8.826273045769994e-15, 0.0, 0.0, 9.992007221626409e-15, 0.0, 2.0872192862952943e-14, 1.0658141036401503e-14, 1.0880185641326534e-14, 1.0935696792557792e-14, 2.2815083156046967e-14, 2.3148150063434514e-14, 1.1379786002407855e-14, 0.0, 1.2212453270876722e-14, 1.2378986724570495e-14, 0.0, 5.551115123125783e-16, 1.199040866595169e-14, -7.216449660063518e-16, -7.216449660063518e-16, 1.3600232051658168e-14],
            "max_abs_residual": 2.3148150063434514e-14
        },
        "terminal": {
            "G": 1100.0,
            "T": 75.0,
            "i_l": 0.38071,
            "i_0": 3.0407536094957474e-05,
            "r_s": 5.32,
            "r_sh": Infinity,
            "nNsVth": 0.04500185335353259,
            "v_gold": [-0.2862349706642596, -0.16501394774169376, -0.07096388773159051, 0.006751639105005369, 0.07343886021077428, 0.13212298794534888, 0.1847047299790537, 0.23246218834751264, 0.2762991532873451, 0.31688004472385467, 0.3547086582705502, 0.39017676033241677, 0.4235954625342937, 0.4552162479857097, 0.48524551139227873, 0.5138548863931274, 0.5411887520407054, 0.5673697999327987, 0.5925032368914325, 0.6166800078985031, 0.6399793026615911],
            "i_gold": [0.1300797193830317, 0.10800768510517172, 0.09084532768216752, 0.0766417656024439, 0.06443936983369158, 0.05369118682117602, 0.04405308373682687, 0.03529338396443271, 0.027248099552826754, 0.0197965724828173, 0.012847245652866268, 0.006328874281743335, 0.00018484733085000205, -0.005630620761174654, -0.011155125227200913, -0.016419943696101058, -0.0214513981060917, -0.02627186524759989, -0.03090053997367187, -0.03535402073282207, -0.03964676513218224],
            "i_test": [0.13007971938303164, 0.10800768510517172, 0.09084532768216752, 0.0766417656024439, 0.06443936983369164, 0.05369118682117596, 0.04405308373682687, 0.035293383964432656, 0.027248099552826643, 0.0197965724828173, 0.012847245652866268, 0.00632887428174328, 0.00018484733085000205, -0.0056306207611747094, -0.011155125227200913, -0.016419943696101058, -0.0214513981060917, -0.026271865247599835, -0.03090053997367187, -0.03535402073282201, -0.03964676513218229],
            "residual": [-5.551115123125783e-17, 0.0, 0.0, 0.0, 5.551115123125783e-17, -5.551115123125783e-17, 0.0, -5.551115123125783e-17, -1.1102230246251565e-16, 0.0, 0.0, -5.551115123125783e-17, 0.0, -5.551115123125783e-17, 0.0, 0.0, 0.0, 5.551115123125783e-17, 0.0, 5.551115123125783e-17, -5.551115123125783e-17],
            "max_abs_residual": 1.1102230246251565e-16
        }
    },
    "v_from_i": {
        "current_sum": {
            "G": 200.0,
            "T": 50.0,
            "i_l": 0.04661,
            "i_0": 1.072213666670876e-06,
            "r_s": 5.32,
            "r_sh": 3000.0,
            "nNsVth": 0.04177035447707613,
            "v_gold": [-0.008983183397948624, 0.028487135828642096, 0.057230737980813406, 0.08156525424438935, 0.10351080049220163, 0.12418597523443886, 0.14428116022719575, 0.16425640667690972, 0.18443643025668055, 0.20506096546862923, 0.2263135377900583, 0.2483389037053799, 0.2712541391388362, 0.2951559737626562, 0.3201258064136182, 0.3462332337899417, 0.3735385950441972, 0.4020948466829859, 0.43194897054437736, 0.463143049174246, 0.495715099714534],
            "i_gold": [0.04622058354841527, 0.04579656743487739, 0.04518025859754816, 0.04435787833493199, 0.04331857973726065, 0.04205326722983692, 0.040554036824137725, 0.038813859444244016, 0.0368263813869985, 0.03458578913575776, 0.03208671287232284, 0.029324154846327385, 0.026293434535373396, 0.02299014560587366, 0.019410121438134843, 0.01554940703406461, 0.011404235789082883, 0.006971010042442646, 0.0022462846113437216, -0.0027732472844729725, -0.008090766158497449],
            "v_test": [-0.008983183397944239, 0.028487135828646704, 0.05723073798081879, 0.08156525424439476, 0.10351080049220585, 0.12418597523444319, 0.1442811602272016, 0.1642564066769161, -0.0888234583382328, 0.20506096546863262, 0.22631353779006247, 0.24833890370537404, 0.271254139138847, 0.2951559737626468, 0.3201258064135999, 0.3462332337899312, 0.3735385950441952, 0.402094846682985, 0.43194897054436865, 0.4631430491742492, 0.4957150997145163],
            "residual": [-4.163336342344337e-17, -9.020562075079397e-17, -1.734723475976807e-16, -2.7755575615628914e-16, -3.191891195797325e-16, -4.579669976578771e-16, -8.326672684688674e-16, -1.1796119636642288e-15, 0.009735069744962271, -9.71445146547012e-16, -1.4363510381087963e-15, 2.4077961846558082e-15, -5.204170427930421e-15, 5.3013149425851225e-15, 1.1879386363489175e-14, 7.806255641895632e-15, 1.6861512186494565e-15, 8.465450562766819e-16, 9.273398021703017e-15, -3.761747857655706e-15, 2.306314861311165e-14],
            "max_abs_residual": 0.009735069744962271
        },
        "terminal": {
            "G": 200.0,
            "T": 50.0,
            "i_l": 0.04661,
            "i_0": 1.072213666670876e-06,
            "r_s": 5.32,
            "r_sh": 3000.0,
            "nNsVth": 0.04177035447707613,
            "v_gold": [-0.008983183397948624, 0.028487135828642096, 0.057230737980813406, 0.08156525424438935, 0.10351080049220163, 0.12418597523443886, 0.14428116022719575, 0.16425640667690972, 0.18443643025668055, 0.20506096546862923, 0.2263135377900583, 0.2483389037053799, 0.2712541391388362, 0.2951559737626562, 0.3201258064136182, 0.3462332337899417, 0.3735385950441972, 0.4020948466829859, 0.43194897054437736, 0.463143049174246, 0.495715099714534],
            "i_gold": [0.04622058354841527, 0.04579656743487739, 0.04518025859754816, 0.04435787833493199, 0.04331857973726065, 0.04205326722983692, 0.040554036824137725, 0.038813859444244016, 0.0368263813869985, 0.03458578913575776, 0.03208671287232284, 0.029324154846327385, 0.026293434535373396, 0.02299014560587366, 0.019410121438134843, 0.01554940703406461, 0.011404235789082883, 0.006971010042442646, 0.0022462846113437216, -0.0027732472844729725, -0.008090766158497449],
            "v_test": [-0.008983183397944239, 0.028487135828646704, 0.05723073798081879, 0.08156525424439476, 0.10351080049220585, 0.12418597523444319, 0.1442811602272016, 0.1642564066769161, -0.0888234583382328, 0.20506096546863262, 0.22631353779006247, 0.24833890370537404, 0.271254139138847, 0.2951559737626468, 0.3201258064135999, 0.3462332337899312, 0.3735385950441952, 0.402094846682985, 0.43194897054436865, 0.4631430491742492, 0.4957150997145163],
            "residual": [4.385380947269368e-15, 4.6074255521944e-15, 5.384581669432009e-15, 5.412337245047638e-15, 4.218847493575595e-15, 4.3298697960381105e-15, 5.856426454897701e-15, 6.38378239159465e-15, -0.27325988859491335, 3.3861802251067274e-15, 4.163336342344337e-15, -5.856426454897701e-15, 1.0769163338864018e-14, -9.381384558082573e-15, -1.8318679906315083e-14, -1.049160758270773e-14, -1.9984014443252818e-15, -8.881784197001252e-16, -8.715250743307479e-15, 3.219646771412954e-15, -1.765254609153999e-14],
            "max_abs_residual": 0.27325988859491335
        }
    }
}]

[cwhanse]

@thunderfish24 that's really odd behavior, at one voltage value the residual is 13 orders of magnitude greater than at any other value. Any chance you can fix the text wrapping in the posted code?

[markcampanelli]

Sorry I posted that hastily on my lunch break. I will try to fix soon.

[markcampanelli]

I updated the formatting. To be clear, these are the four I-V curves out of the (G,T) matrix that produced the worst maximum absolute residual across the I-V values in each curve for the specific function (i_from_v() or v_from_i()), for both the residual for the current sum at the diode node (ideally zero) and for the calculated terminal value's difference from the gold value.

I probably should rename i_test_abs_res to be i_test_terminal_res, etc. UPDATED names.

[markcampanelli]

I updated the json output above to have signed residuals that are consistent between the various computations.

[cwhanse]

@thunderfish24 I think there's something amiss in your calculation of v_test. In the 8th position you get a value of -0.0888234583382328. Using the existing pvlib.pvystem.v_from_i I get a value of 0.18443643025668521; all other values are exactly the same as you report. Whatever the mistake, it is causing the large residual.

[markcampanelli]

@cwhanse Thanks for checking. I cannot figure out why I'm computing these anomalous values. I have opened #426 so that others can run the same exact script and compare output, including comparison of python (3.6.4) and numpy (1.14.0) versions.

[cwhanse]

@pvlib/pvlib-maintainer expect a pull request in the near future from @reepoi that should close this issue. @markcampanelli