Add Huld PV model

[cwhanse]

- [ ] Closes #xxxx

• I am familiar with the contributing guidelines
• Tests added
• Updates entries in docs/sphinx/source/reference for API changes.
• Adds description and name entries in the appropriate "what's new" file in docs/sphinx/source/whatsnew for all changes. Includes link to the GitHub Issue with :issue:`num` or this Pull Request with :pull:`num`. Includes contributor name and/or GitHub username (link with :ghuser:`user`).
• New code is fully documented. Includes numpydoc compliant docstrings, examples, and comments where necessary.
• Pull request is nearly complete and ready for detailed review.
• Maintainer: Appropriate GitHub Labels (including remote-data) and Milestone are assigned to the Pull Request and linked Issue.

Adds Huld's model for DC power. This model is used in PVGIS.

[AdamRJensen]

I think it will be great addition having the Huld/PVGIS model available!

[adriesse]

FYI: This and other efficiency models are available in https://github.com/adriesse/pvpltools-python. (#1544)

[cwhanse]

Would pvsystem.pvgis be more widely recognized that huld?

[mikofski]

Would pvsystem.pvgis be more widely recognized that huld?

I’m not fan of pvgis as a subheading, module, subpackage, namespace.

also sometimes I use $25^{\circ}$

[adriesse]

Would pvsystem.pvgis be more widely recognized that huld?

I think so, in the same way pvsyst is more recognizable than... what's his name again?

[AdamRJensen]

Would pvsystem.pvgis be more widely recognized that huld?

I think so, in the same way pvsyst is more recognizable than... what's his name again?

I'm inclined to disagree as I wouldn't be surprised if pvgis one day switched to another PV performance model.

[mikofski]

Desoto is a woman

[adriesse]

Desoto is a woman

?? Mermoud is a man

[kandersolar]

@cwhanse I noticed that the current implementation doesn't exhibit the temperature dependence that I expect. Here's the temperature coefficient of power at STC:

pdc0 = 300
dT = 0.1
gamma = (huld(1000, 25 + dT/2, pdc0, cell_type='cSi') - huld(1000, 25 - dT/2, pdc0, cell_type='cSi')) / dT / pdc0
print(gamma)  # -1.56733e-05

I think this is because the k coefficient values from the PVGIS website are for a different form of Huld's equation which uses them differently from the 2011 paper. Compare:

• PVGIS website: $P(G, T) = G' \cdot A \cdot \text{eff}_\text{nom} \cdot (1 + k_1 \log G' + k_2 \log^2 G' + ...)$
• Huld et al. 2011: $P(G, T) = G' \cdot (P_\text{stc} + k_1 \log G' + k_2 \log^2 G' + ...)$

The two versions of $k_n$ differ by a factor of $P_\text{stc} = A \cdot \text{eff}_\text{nom}$. Right now we're using the website's coefficients with the paper's equation, thus the strange behavior. We need to choose which of the two forms to adopt here. Not sure which one is better.

[AdamRJensen]

It would be nice, as a minimum, to be able to replace the PVGIS website/API calculations.

[cwhanse]

The two versions of $k_n$ differ by a factor of $P_\text{stc} = A \cdot \text{eff}_\text{nom}$.

Good catch. I noticed the notation difference $k$ and $k'$ in the paper but overlooked what it meant.

Not sure which one is better.

I prefer Huld's paper, since the PVGIS technical document is harder to locate.

[cwhanse]

It would be nice, as a minimum, to be able to replace the PVGIS website/API calculations.

I agree, and my attempt fails. I don't know how to go about diagnosing where this code diverges from PVGIS.

from pvlib.iotools import pvgis
from pvarray import huld
from pvlib.temperature import faiman


res = pvgis.get_pvgis_hourly(45, 2, components=False, surface_tilt=35,
                             usehorizon=False, pvcalculation=True,
                             peakpower=1, pvtechchoice='crystSi')

data = res[0]

geff = data['poa_global']
ta = data['temp_air']
ws = data['wind_speed']
pdc = data['P']

tmod = faiman(geff, ta, ws, u0=26.9, u1=6.2)
expected = huld(geff, tmod, 1000, cell_type='csi')

print(pdc[:16])
print(expected[:16])

[kandersolar]

IAM and SMM, I'm guessing? Their docs say PVGIS uses Martin-Ruiz for IAM, so at least that is easily recreated. However, the spectral model (SPECMAGIC?) looks decidedly nontrivial to recreate, at least from a quick skim. https://doi.org/10.3390/rs4030622

[kandersolar]

Indeed including IAM gets a lot closer:

Code:

Click to expand!

from pvlib.iotools import pvgis
from pvarray import huld
from pvlib.temperature import faiman
from pvlib.iam import martin_ruiz
import matplotlib.pyplot as plt


# use surface_tilt = 0 so that zenith = aoi.  Means we don't have to compute solar position for IAM
df, _, meta = pvgis.get_pvgis_hourly(45, 2, components=True, surface_tilt=0,
                                     usehorizon=False, pvcalculation=True,
                                     peakpower=1, pvtechchoice='crystSi')

# %%

fig, axes = plt.subplots(1, 2, sharex=True, sharey=True, figsize=(8, 4))

for (iam, label), ax in zip([(1, "no IAM"),
                             (martin_ruiz(90 - df['solar_elevation']), "with Martin-Ruiz")], axes):
    smm = 1
    geff = smm * (iam * df['poa_direct'] + df['poa_sky_diffuse'] + df['poa_ground_diffuse'])
    ta = df['temp_air']
    ws = df['wind_speed']
    pdc = df['P']
    
    tmod = faiman(geff, ta, ws, u0=26.9, u1=6.2)
    expected = huld(geff, tmod, 1000, cell_type='csi')
    
    ax.scatter(expected, pdc, s=1)
    ax.axline((0, 0), slope=1, c='k', ls=':')
    ax.set_title(label)

axes[0].set_xlim(0, 1000)
axes[0].set_ylim(0, 1000)
axes[0].set_ylabel('PVGIS $P_{dc}$')
for ax in axes:
    ax.set_xlabel('pvlib $P_{dc}$')

fig.tight_layout()

[kandersolar]

We should consider documenting (in the huld docstring) which convention the input k is supposed to follow. The way it is now, there is little to alert a user that manually passing the k values from the PVGIS document is incorrect.

[AdamRJensen]

@AdamRJensen I don't think we can reproduce the PVGIS API's results, since the spectral adjustment is not exposed nor documented in detail (it appears that PVGIS has an internal look-up table of sorts for spectral adjustment factors). The scope of pvarray.huld is only the conversion of effective_irradiance and cell_temp to pdc. Perhaps a full replication of PVGIS could be realized later, with something like ModelChain.with_pvwatts. Agree?

Agreed! I imagine that some future users will ask why they don't match. So I suggest adding an admonition stating that the PVGIS API also applies a spectral correction.

[adriesse]

A thought for consideration: As Huld notes in his paper, the model can also be expressed as an efficiency model (eq.2), which has Pstc factored out. A function pvefficiency_huld() would not only be a nice complement to the existing functions, but could also use the published parameter values as-is.

[adriesse]

Suggested change

	def huld(effective_irradiance, temp_mod, pdc0, k=None, cell_type=None):
	def pvpower_huld(effective_irradiance, temp_mod, pdc0, k=None, cell_type=None):

[cwhanse]

I'm not in favor of this renaming. Efficiency can be easily calculated, and this function can be wired into ModelChain in the same way that other DC models are.

[adriesse]

Wasn't there was some discussion about avoiding function names consisting only of an author's name?

[adriesse]

Most model functions have separate arguments for each parameter. Are there pros and cons?

[cwhanse]

My personal preference is to use tuples for polynomial coefficients. It would be tedious to list k1, k2, etc. and define each as the "nth of 6 coefficients..."

[kandersolar]

One consideration: ModelChain's ability to infer models based on supplied parameters relies on the models having distinct sets of parameter names. k1, k2 etc seems less likely to run into future conflict than a single k does. But neither seems particularly likely to me.

I also favor a single k. The individual values are not very meaningful on their own (IMHO), so a single bundled parameter seems appropriate to me.

[adriesse]

Does the model use cell or module temperature?

[cwhanse]

Module. Good catch.

[adriesse]

Technically logG tends to -inf as G tends to 0.

[kandersolar]

I notice that small but nonzero irradiance results in negative power values:

I'm not immediately seeing any provision for handling this in the reference. We could constrain the output to be non-negative. Either way I think it deserves a docstring note.

[adriesse]

Suggested change

	:math:`G' < 20 W/m^2`, :math:`P_{dc}` may be negative.
	:math:`G' < 20 W/m^2`, :math:`P_{dc}` may be negative.
	This is caused by the :math:`\log(G')` term, which tends toward
	negative infinity as irradiance approaches zero.

I can provide a reference that discusses this, if it is of interest.

[kandersolar]

Discussion has died down, so I assume no further changes are needed. Thanks @cwhanse and reviewers!