Auxiliary equation for PVsyst shunt resistance may be incorrect

[markcampanelli]

Describe the bug
The implementation of PVsyst's auxiliary equation for the shunt resistance as a function of irradiance at

https://github.com/pvlib/pvlib-python/blob/v0.8.0/pvlib/pvsystem.py#L1326-L1331

does not match the description given by equation (5) in this reference

https://www.pvsyst.com/wp-content/uploads/2019/01/thf_article_valence_2010.pdf

It is unclear from the code what pvlib's version is trying to compensate for and why it deviates from at least one of the references given and PVsyst's online documentation at

https://www.pvsyst.com/help/index.html?pvmodule_rshexp.htm

My guess is that it might be trying to make shunt resistance at the reference irradiance precisely equal to Rshunt. See below.

The PVsyst documentation/notation is confusing because Rsh(Gref) (also denoted as Rshunt) is not precisely the shunt resistance at reference conditions. Rather, it appears that it is a reasonable approximation of the shunt resistance at the reference irradiance for larger values of the exponent parameter, e.g., for the common PVsyst default value of 5.5. However, this approximation can be quite poor, as evidenced by the second example below taken from an actual PAN file:

import numpy

def R_sh_at_G(*, G, R_sh_ref, R_sh_zero, R_sh_exp, G_ref=1000.):
    return R_sh_ref + (R_sh_zero - R_sh_ref) * numpy.exp(-R_sh_exp * (G / G_ref))

G = 1000.

R_sh_ref, R_sh_zero, R_sh_exp = 9000, 10000, 5.5

print(f"R_sh_at_G(G={G}, R_sh_ref={R_sh_ref}, R_sh_zero={R_sh_zero}, R_sh_exp={R_sh_exp}) = {R_sh_at_G(G=G, R_sh_ref=R_sh_ref, R_sh_zero=R_sh_zero, R_sh_exp=R_sh_exp)} vs. Rshunt={R_sh_ref}")

R_sh_ref, R_sh_zero, R_sh_exp = 900, 960000, 1.10

print(f"R_sh_at_G(G={G}, R_sh_ref={R_sh_ref}, R_sh_zero={R_sh_zero}, R_sh_exp={R_sh_exp}) = {R_sh_at_G(G=G, R_sh_ref=R_sh_ref, R_sh_zero=R_sh_zero, R_sh_exp=R_sh_exp)} vs. Rshunt={R_sh_ref}")

which gives

R_sh_at_G(G=1000.0, R_sh_ref=9000, R_sh_zero=10000, R_sh_exp=5.5) = 9004.086771438464 vs. Rshunt=9000
R_sh_at_G(G=1000.0, R_sh_ref=900, R_sh_zero=960000, R_sh_exp=1.1) = 320156.6563748281 vs. Rshunt=900

Furthermore, for these parameters, pvlib.pvsystem.calcparams_pvsyst() calculates 9000.000000000002 and 319556.2403501564, respectively

This discrepancy is particularly concerning for algorithms that make the assumption that Rshunt in the PAN file is actually the value at the reference irradiance, such as (cc @frivollier)

https://github.com/frivollier/pvsyst_tools/blob/master/pvsyst/module.py#L275-L278

To Reproduce
See code snippets above.

Expected behavior
Implement the auxiliary equation for shunt resistance given in the PVsyst reference and online documentation, or better justify and document the difference.

Screenshots
N/A

Versions:

• pvlib.__version__: 0.8.0
• pandas.__version__: 1.1.3
• python: 3.8.5

Additional context
N/A

[cwhanse]

The Rshunt equation coded in pvlib comes from the first reference

[1] K. Sauer, T. Roessler, C. W. Hansen, Modeling the Irradiance and
       Temperature Dependence of Photovoltaic Modules in PVsyst,
       IEEE Journal of Photovoltaics v5(1), January 2015

Ken (first author) confirmed its correctness with Pvsyst staff before we published that paper.

[markcampanelli]

Thanks for clarifying @cwhanse. I would say that the online PVsyst documentation is misleading.

Would you/anyone be opposed to factoring out the R_sh calculation into a separately accessible function?

[markcampanelli]

As an aside, I wonder if a simple 2-parameter model such as

def R_sh_at_G(*, G, R_sh_zero, R_sh_ref, G_ref=1000.):
    return R_sh_zero * (R_sh_ref / R_sh_zero)**(G/G_ref)

would have sufficed?

[cwhanse]

I'm not opposed on principle to a R_sh helper function but what advantage does that have? No other module model uses that particular expression.

[markcampanelli]

It would help to compute the actual shunt resistance at STC in situations like this: frivollier/pvsyst_tools#2.

It appears that the current implementation does not always make the shunt resistance at STC equal/close to Rshunt.

[cwhanse]

I see the point. There's a helper function here so pvlib has the same calculation twice. I'm open to ideas to improve this.

Under what conditions is Rsh@G=1000 not equal to Rsh_ref? The algebra works out that way. Maybe its a matter of roundoff error?

[markcampanelli]

@cwhanse Sorry this detail got lost in my above examples, because I was focused on the discrepancy between the PVsyst documentation and the the actual implementation duplicated in pvlib. My understanding is that the example below shows a considerable discrepancy between the Rshunt=900 Ohms taken from the PAN file and the reference value of 319556.2403501564 Ohms computed by pvlib.pvsystem.calcparams_pvsyst:

pvsyst_params = {'alpha_sc': 0.00039522258414766555, 'gamma_ref': 1.017, 'mu_gamma': -0.0003, 'I_L_ref': 18.402482849369548, 'I_o_ref': 5.1690496892812056e-11, 'R_sh_ref': 900.0, 'R_sh_0': 960000.0, 'R_s': 0.15, 'cells_in_series': 60, 'R_sh_exp': 1.1, 'EgRef': 1.12, 'irrad_ref': 1000, 'temp_ref': 25}
pvlib.pvsystem.calcparams_pvsyst(effective_irradiance=1000., temp_cell=25., **pvsyst_params)

which for me produced

(18.402482849369548,
 5.1690496892812056e-11,
 0.15,
 319556.2403501564,
 1.5677606661864085)

Unfortunately, I don't have any insight into how the PAN values for the shunt resistance were generated.

[cwhanse]

I get the same results, unsurprisingly. What's odd is the magnitude of the Rsh_0 in that PAN file, and the departure from the default of 5.5 for the exponent. If the fitting tools was misused...

I'll write Ken and Bruno and see what I can learn about the correctness of that expression. The intend of calcparams_pvsyst is to implement what Pvsyst does (well, did as of v6).

[markcampanelli]

I can verify that the values are correct from the PAN file. I'm checking with my source to see if I can share the full PAN file and a supporting PDF file. Hopefully I can, then I could perhaps ask them to track down the origin of the file, but first it would be nice to have all the cards on the table.

  RShunt=900
  Rp_0=960000
  Rp_Exp=1.10

UPDATE: This is all the information that I will be able to provide.

[adriesse]

Good discussion. There are data and PAN files from CFV labs here if it useful:

(https://pvpmc.sandia.gov/pv-research/pv-lifetime-project/pv-lifetime-modules/)

[markcampanelli]

@cwhanse I notice here that the "fitting" function and the "forward model" function each implement their own copy of the shunt resistance calculation. So, we could DRY up the code by using (and also publicly exposing) a helper function as the single source of truth. Thoughts?

Also, were you able to verify that this implementation remains valid in PVsyst v7?

[markcampanelli]

Good discussion. There are data and PAN files from CFV labs here if it useful:

(https://pvpmc.sandia.gov/pv-research/pv-lifetime-project/pv-lifetime-modules/)

Thanks @adriesse. All of these PAN files have the PVsyst-recommended value of 5.5 for R_sh_exp. However, it is still a free-parameter in PVsyst itself, so I think accommodation of possible variation is indeed proper in pvlib-python and pvsyst_tools.

[cwhanse]

Helper functions would be an improvement, but let's make sure they are traceable to a specific module performance model (and possibly to a version of that model).

[adriesse]

After rereading this long-forgotten thread, I wonder whether @cwhanse got any clarification or confirmation from Bruno and/or Ken...

@markcampanelli do you find that this exponent gives you a useful degree of freedom for fitting?

[campanelli-sunpower]

@adriesse Good question. I am not quite there yet with PVfit. (Thanks to @cwhanse's detective work, the actual implementation of the PVsyst model appears to be nailed down better in pvlib than in PVsyst's own help documentation!) Nevertheless, when I cross that bridge I'm very inclined to take PVsyst's advice here: "The RshExp exponential coefficient is always -5.5, according to all our measurements of Rshunt at sun. This should not be modified." ref

I was never able to track down the origin of non-standard value (Rp_Exp=1.10) that triggered this issue.

My attention is refocused here because I am working on a translation from the PVsyst model (via PAN file) to SAPM (for SunPower's PVAPI simulation tool). This is why you'll see me posting using two different Github users :).

[cwhanse]

After rereading this long-forgotten thread, I wonder whether @cwhanse got any clarification or confirmation from Bruno and/or Ken...

I didn't get any further clarification about v6. I haven't asked about v7.

[markcampanelli]

@cwhanse I would like to know if the maintainers of pvlib would be open to moving the https://github.com/frivollier/pvsyst_tools into pvlib. I don't think @frivollier has much time now to maintain/improve, and there is duplicated functionality.

[cwhanse]

@markcampanelli I'd say yes. It would be necessary to bring some of the PAN and OND files for test purposes, and to explicitly state the Pvsyst versions for which this read capability works. The code should go in a new module in pvlib-python/iotools.

[adriesse]

Looks a like substantial amount of work, but of course it will be useful to many people!

[AdamRJensen]

What functions from pvsyst_tools should be moved? It seems pan_to_module_param is the main function.

[markcampanelli]

@AdamRJensen There is the associated ond_to_inverter_param, but it does not appear to do any "model parameter reconstruction" like the module code does (e.g., reverse saturation current). In terms of upgrades, I have noticed a few entries in the PAN files I'm working with are missing in the Python dictionary. We might also think about adding some sort of formal file/model versioning scheme. I have reached out to @frivolier for his blessing to make the move.

[kandersolar]

Under what conditions is Rsh@G=1000 not equal to Rsh_ref? The algebra works out that way. Maybe its a matter of roundoff error?

My understanding is that the example below shows a considerable discrepancy between the Rshunt=900 Ohms taken from the PAN file and the reference value of 319556.2403501564 Ohms computed by pvlib.pvsystem.calcparams_pvsyst

The algebra works out that way iff the Rsh_base = np.maximum(0.0, Rsh_tmp) line does not clip negative Rsh_tmp values to zero. The large Rsh0 and small Rshexp values in @markcampanelli's PAN file indeed end up triggering that clip, resulting in the variance between Rsh(Gstc) and Rsh_ref.

Since Rsh_base is what the exponential curve converges towards as irradiance tends to infinity, it must be greater than zero for the modeled shunt resistance to be positive at all irradiances. The given example parameters are a good example: without the clip, the enormous Rsh0 and steep Rshexp result in negative Rsh from 1001.7 W/m2 onwards.