Dividing by cos(solar_zenith) in haydavies and reindl functions

[josricha]

A few issues with the pvlib package seem to be caused by dividing by cos(solar_zenith) since solar zenith can vary from 0-180 degrees. As the solar_zenith approaches 90 degrees, these results will approach infinity. The problem becomes very apparent when the value of solar_zenith is between 89 and 90 degrees which is rare unless working with minute interval (or smaller) because the denominator then becomes a large scaling factor. I have attached a quick sample code to illustrate the issue.
cos(solar_zenith)_errors.txt

A few questions and comments to note:

1. This seems to have been addressed in a similar context in the DISC Model and Zenith Threshold issue but appear again in the haydavies algorithm implementation on line 729 of irradiance.py when calculating the projection_ratio and could potentially appear elsewhere in the package.

2. The MATLAB version of the haydavies algorithm, pvl_haydavies1980.m, uses the formula "RB = max(COSTT, 0) ./ max(cosd(SunZen), 0.01745))" where 0.01745 corresponds to a zenith angle of 89 degrees to avoid the issue.

3. Sandia National Laboratories documentation seems to suggest that dividing by cos(solar_zenith) is not part of the Hay and Davies diffuse model.

Is there a documented reason for the differences between these formulas or is this a mistake in implementation?

[wholmgren]

This seems to have been addressed in a similar context in the DISC Model and Zenith Threshold

A related issue was discussed and a potential solution was proposed, but it was not merged.

could potentially appear elsewhere in the package.

I see that a similar issue could come up in reindl. I don't see any other problems in the irradiance.py module. I'm going to go ahead and change this issue to be specific to the transposition functions. Please post separate issues if you find them in the rest of the library.

I'll post a few haydavies and reindl calculations that use only slightly unphysical input data:

# surface_tilt, surface_azimuth, dhi, dni, dni_et, solar_zenith, solar_azimuth
In [148]: pvlib.irradiance.haydavies(30, 180, 50, 0, 1360, 90, 270)
Out[148]: 46.65063509461097

In [149]: pvlib.irradiance.haydavies(30, 180, 50, 1, 1360, 90, 270)
Out[149]: 46.66655467923936

In [150]: pvlib.irradiance.haydavies(30, 180, 50, 1, 1360, 90, 230)
Out[150]: 192969086526870.06

In [151]: pvlib.irradiance.haydavies(30, 180, 50, 0, 1360, 90, 230)
Out[151]: 46.65063509461097

In [152]: pvlib.irradiance.haydavies(30, 180, 50, 1, 1360, 89.9, 230)
Out[152]: 53.41820537086939

In [153]: pvlib.irradiance.haydavies(30, 180, 50, 1, 1360, 89.99, 230)
Out[153]: 114.34857082698177

In [154]: pvlib.irradiance.haydavies(30, 180, 50, 1, 1360, 89.999, 230)
Out[154]: 723.6521641382168

In [158]: pvlib.irradiance.reindl(30, 180, 50, 1, 50, 1360, 89.999, 230)
Out[158]: 723.6526416460301

So the functions behave ok for some problematic inputs but not for others. I think it's in the scope of pvlib to be able to gracefully handle these inputs.

The MATLAB version of the haydavies algorithm, pvl_haydavies1980.m, uses the formula "RB = max(COSTT, 0) ./ max(cosd(SunZen), 0.01745))" where 0.01745 corresponds to a zenith angle of 89 degrees to avoid the issue.

We could do this here for consistency. It seems somewhat arbitrary but pragmatic to me. I looked at the git history and I saw that the original port did contain this code, but I removed it as one of many changes to clean up and simplify the library. My bad.

Sandia National Laboratories documentation seems to suggest that dividing by cos(solar_zenith) is not part of the Hay and Davies diffuse model.

The pvlib implementation roughly follows the Loutzenhiser 2007 description of the model. From a numerical standpoint, the PVPMC website's description would be a better implementation. There may be some legacy here that I'm not aware of, such as the MATLAB function could have at one point returned the total POA irradiance, rather than just the sky diffuse component.

Is there a documented reason for the differences between these formulas or is this a mistake in implementation?

See Loutzenhiser 2007. The implementation is correct, it's just not numerically the best way to do it.

[adriesse]

Sandia National Laboratories documentation seems to suggest that dividing by cos(solar_zenith) is not part of the Hay and Davies diffuse model.

The pvlib implementation roughly follows the Loutzenhiser 2007 description of the model. From a numerical standpoint, the PVPMC website's description would be a better implementation. There may be some legacy here that I'm not aware of, such as the MATLAB function could have at one point returned the total POA irradiance, rather than just the sky diffuse component.

I think the Sandia National Laboratories documentation is wrong on this one. Just imagine 100% circumsolar diffuse: tilted diffuse could never be greater than horizontal diffuse with this formula, even with the array pointed at the sun.

[adriesse]

I see that a similar issue could come up in reindl. I don't see any other problems in the irradiance.py module. I'm going to go ahead and change this issue to be specific to the transposition functions. Please post separate issues if you find them in the rest of the library.

Essentially the same problem occurs in the function dni (Determine DNI from GHI and DHI). Practical solutions for that particular case were discussed at length in #250.

[adriesse]

So the functions behave ok for some problematic inputs but not for others. I think it's in the scope of pvlib to be able to gracefully handle these inputs.

I agree with this. Maybe focus on the ratio AI/cos(zenith) which should become 0/0 at the horizon.

[cwhanse]

I think the Sandia National Laboratories documentation is wrong on this one. Just imagine 100% circumsolar diffuse: tilted diffuse could never be greater than horizontal diffuse with this formula, even with the array pointed at the sun.

@adriesse you are correct, thanks for catching that. We wrote cos(AOI) where we should have written Rb. And we should give the equation for Rb = cos(AOI) / cos(zenith).

Loutzenhiser doesn't give an explicit expression for Rb. I could dig up the original Hay and Davies paper if we need to confirm, but others have written out the expression Rb = cos(AOI) / cos(zenith)

https://www.researchgate.net/publication/275276907_Comparative_study_of_isotropic_and_anisotropic_sky_models_to_estimate_solar_radiation_incident_on_tilted_surface_A_case_study_for_Bhopal_India

for example.

[cwhanse]

The handling of zenith near 90 in PVLib for Matlab was our addition to Hay and Davies. An alternative value for Hay and Davies at zenith = 90 would be the isotropic sky diffuse component, reasoning that at zenith = 90 there is not direct or circumsolar irradiance that reaches a tilted surface near ground level. We could linearly interpolate from 89 (or 88) to 90.

I am confident that the Hay and Davies model did not anticipate being applied at zenith = 90 and that we won't find any guidance in the original paper.

[wholmgren]

The Hay-Davies model is equal to the isotropic model when DNI is 0. So, I think the matlab approach of ensuring that Rb doesn't explode while DNI goes to 0 at the horizon is a kind of implicit interpolation. Maybe that's not the right way to put it, but hopefully it gets the point across.

There is one more difference with the matlab implementation. In pvl_haydavies1989.m:

RB = max(COSTT,0)./max(cosd(SunZen),0.01745);

Our python haydavies does not restrict the numerator to positive values:

        cos_tt = aoi_projection(surface_tilt, surface_azimuth,
                                solar_zenith, solar_azimuth)
        cos_solar_zenith = tools.cosd(solar_zenith)
        Rb = cos_tt / cos_solar_zenith

The same discrepancy is also present in the Reindl implementations. So, that's another problem with the python implementations.