Optimal surface_azimuth is not 180° using pvlib??

[Casey McGonigle]

Hi all and thank you for taking the time to read this and help me out,

I'm building a model using the open-source pvlib software (and CEC Modules) to estimate yearly photovoltaic energy output. I have encountered some inconsistencies in the model and would appreciate any troubleshooting the community can offer.

My main problem is this: the model tells me that the ideal Northern hemisphere surface_azimuth for energy production (ie. surface_azimuth with the highest energy output) is around 76° (just North of due-East) while the worst surface_azimuth for energy production is around 270° (due West). However, I understand that the ideal surface_azimuth in the Northern hemisphere should be 180° (due South) with the worst surface_azimuth at 0° (due North).

I've attached a photo that illustrates how energy output changes with surface_azimuth 

Can anyone help me rectify this issue or correct my understanding?

I've copied my code here as well in order to help troubleshooting:

```

import os

import pandas as pd

import numpy as np

import os

import os.path

import matplotlib.pyplot as plt

import pvlib

from geopy.exc import GeocoderTimedOut 

from geopy.geocoders import Nominatim 

from IPython.display import Image

## GET THE LATITUDE & LONGITUDE OF A GIVEN CITY

geolocator = Nominatim(user_agent="____'s app") 

geo = geolocator.geocode("Berkeley") 

## CHECK THAT CITY IS CORRECT (by Country, State, etc.)

print(geo.address)

# CHECK THE LAT, LON order

print(geo.latitude, geo.longitude)

## SELECT THE YEAR & TIMES YOU'D LIKE TO MODEL OFF

YEAR = 2019

STARTDATE = '%d-01-01T00:00:00' % YEAR

ENDDATE = '%d-12-31T23:59:59' % YEAR

TIMES = pd.date_range(start=STARTDATE, end=ENDDATE, freq='H')

## ACCESS THE NREL API TO EXTRACT WEATHER DATA

NREL_API_KEY = os.getenv('NREL_API_KEY', 'DEMO_KEY')

## FILL IN THE BLANK WITH YOUR EMAIL ADRRESS

EMAIL = os.getenv('EMAIL', '_______.com')

##NEED TO COMMENT OUT THIS LINE BELOW -- if you call it too many times within an hour, it will break your code

header, data = pvlib.iotools.get_psm3(LATITUDE, LONGITUDE, NREL_API_KEY, EMAIL)

## SELECT THE PVLIB PANEL & INTERVTER YOU'D LIKE TO USE

## CAN ALSO CHOOSE FROM SANDIA LABS' DATASET OF PANELS & INVERTERS (check out the function)

## WE CHOSE THE CECMods because they highlighted the ones that were BIPV

INVERTERS = pvlib.pvsystem.retrieve_sam('CECInverter')

INVERTER_10K = INVERTERS['SMA_America__SB10000TL_US__240V_']

CECMODS = pvlib.pvsystem.retrieve_sam('CECMod')

## SELECT THE PANEL YOU'D LIKE TO USE (NOTE: THE PEVAFERSA MODEL IS A BIPV PANEL)

CECMOD_MONO = CECMODS['Pevafersa_America_IP_235_GG']

## CREATING AN ARRAY TO ITERATE THROUGH IN ORDER TO TEST DIFFERENT SURFACE_AZIMUTHS

heading_array = np.arange(0, 361, 2)

heading_array

heading_DF = pd.DataFrame(heading_array).rename(columns = {0: "Heading"})

heading_DF.head()

# geo IS AN OBJECT (the given city) CREATED ABOVE

LATITUDE, LONGITUDE = geo.latitude, geo.longitude

# data IS AN OBJECT (the weather patterns) CREATED ABOVE

# TIMES IS ALSO CREATED ABOVE, AND REPRESENTS TIME

data.index = TIMES

dni = data.DNI.values

ghi = data.GHI.values

dhi = data.DHI.values

surface_albedo = data['Surface Albedo'].values

temp_air = data.Temperature.values

dni_extra = pvlib.irradiance.get_extra_radiation(TIMES).values

# GET SOLAR POSITION

sp = pvlib.solarposition.get_solarposition(TIMES, LATITUDE, LONGITUDE)

solar_zenith = sp.apparent_zenith.values

solar_azimuth = sp.azimuth.values

# CREATING THE ARRY TO STORE THE DAILY ENERGY OUTPUT BY SOLAR AZIMUTH

e_by_az = []

# IDEAL surface_tilt ANGLE IN NORTHERN HEMISPHERE IS ~25

surface_tilt = 25

# ITERATING THROUGH DIFFERENT SURFACE_AZIMUTH VALUES

for heading in heading_DF["Heading"]:

    surface_azimuth = heading

    poa_sky_diffuse = pvlib.irradiance.get_sky_diffuse(

        surface_tilt, surface_azimuth, solar_zenith, solar_azimuth,

        dni, ghi, dhi, dni_extra=dni_extra, model='haydavies')

    # calculate the angle of incidence using the surface_azimuth and (hardcoded) surface_tilt

    aoi = pvlib.irradiance.aoi(

        surface_tilt, surface_azimuth, solar_zenith, solar_azimuth)

    # https://pvlib-python.readthedocs.io/en/stable/generated/pvlib.irradiance.aoi.html

    # https://pvlib-python.readthedocs.io/en/stable/generated/pvlib.pvsystem.PVSystem.html

    poa_ground_diffuse = pvlib.irradiance.get_ground_diffuse(

        surface_tilt, ghi, albedo=surface_albedo)

    poa = pvlib.irradiance.poa_components(

        aoi, dni, poa_sky_diffuse, poa_ground_diffuse)

    poa_direct = poa['poa_direct']

    poa_diffuse = poa['poa_diffuse']

    poa_global = poa['poa_global']

    iam = pvlib.iam.ashrae(aoi)

    effective_irradiance = poa_direct*iam + poa_diffuse

    temp_cell = pvlib.temperature.pvsyst_cell(poa_global, temp_air)

    # THIS IS THE MAGIC

    cecparams = pvlib.pvsystem.calcparams_cec(

        effective_irradiance, temp_cell,

        CECMOD_MONO.alpha_sc, CECMOD_MONO.a_ref,

        CECMOD_MONO.I_L_ref, CECMOD_MONO.I_o_ref,

        CECMOD_MONO.R_sh_ref, CECMOD_MONO.R_s, CECMOD_MONO.Adjust)

    # mpp is the list of energy output by hour for the whole year using a single panel

    mpp = pvlib.pvsystem.max_power_point(*cecparams, method='newton')

    mpp = pd.DataFrame(mpp, index=TIMES)

    first48 = mpp[:48]

    Edaily = mpp.p_mp.resample('D').sum()

    # Edaily is the list of energy output by day for the whole year using a single panel

    Eyearly = sum(Edaily)

    e_by_az.append(Eyearly)

## LINKING THE Heading (ie. surface_azimuth) AND THE Eyearly (ie. yearly energy output) IN A DF

heading_DF["Eyearly"] = e_by_az

heading_DF.head()

## VISUALIZE ENERGY OUTPUT BY SURFACE_AZIMUTH

plt.plot(heading_DF["Heading"], heading_DF["Eyearly"])

plt.xlabel("Surface_Azimuth Angle")

plt.ylabel("Yearly Energy Output with tilt @ " + str(surface_tilt))

plt.title("Yearly Energy Output by Solar_Azimuth Angle using surface_tilt = " + str(surface_tilt) + " in Berkeley, CA");

# FIND SURFACE_AZIMUTH THAT YIELDS THE MAX ENERGY OUTPUT

heading_DF[heading_DF["Eyearly"] == max(heading_DF["Eyearly"])]

# FIND SURFACE_AZIMUTH THAT YIELDS THE MIN ENERGY OUTPUT

heading_DF[heading_DF["Eyearly"] == min(heading_DF["Eyearly"])]

```

Thank you for any help you can offer,

Casey

[Kevin Anderson]

Hi Casey,

When summary results don't seem right, it's often a good idea to examine the underlying timeseries data.  The first plot below shows data from your script and shows a couple days of power, irradiance, and calculated solar position.  In particular, notice the mismatch between irradiance and solar zenith -- irradiance peaks around noon as expected, but solar zenith is lowest around 19:00 or so, which is incorrect -- zenith should be lowest around solar noon as well.  This is because your TIMES variable is not timezone-aware and so the solar position calculation is assuming UTC.  The second plot below shows a couple days where TIMES is tz-localized with tz='Etc/GMT+8' to match the PSM3 data and things look much better.

Also, you mentioned the "ideal surface_azimuth in the Northern hemisphere should be 180°" -- this is mostly true, but local weather patterns might shift the optimal azimuth off of 180 a bit.  For example a place might tend to have cloudy mornings and sunny afternoons, in which case it would be better to face somewhat west of south.

One other thing, when using hourly weather data, it is usually better to calculate solar position for the middle of the interval so that it's more representative of the actual solar position across each hour.  See for example the 30-minute adjustments here: https://pvlib-python.readthedocs.io/en/stable/auto_examples/plot_transposition_gain.html

Cheers,

Kevin

--
You received this message because you are subscribed to the Google Groups "pvlib-python" group.
To unsubscribe from this group and stop receiving emails from it, send an email to pvlib-python...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/pvlib-python/d251dc74-2b36-497e-b3a6-1e9d533d1fe4n%40googlegroups.com.

[Casey McGonigle]

Kevin, 

Thanks for getting back to me so quickly. I'll make sure to localize the timezone with tz='Etc/GMT+8'.

However, this doesn't address the main issue I'm running into. I understand and accept that the optimal surface_azimuth may vary from due South by a couple of degrees, but this code claims the optimal surface azimuth is more than 90° away from due South. How do I reconcile that and what could be causing this large discrepancy between what we expect the optimal surface azimuth to be and what our code is outputting?

Thanks,

Casey

[kevina...@gmail.com]

Hi Casey,

I'm confused -- do you still get ~76 as the optimal azimuth after correcting the timezone localization?  I think the timezone issue actually is the main problem here.  Because of that inconsistency, there is a big timing mismatch between the solar position and your weather data, meaning the transposed irradiance is incorrect and the modeled power is therefore also incorrect.  Take a look at the first plot I posted above -- power output clearly doesn't match irradiance, so the power and anything that depends on it is bogus. As long as that modeling issue is present, any model results are untrustworthy.  After correcting the timezone issue, the power timeseries is much more reasonable (see second plot above) and I get a much more reasonable optimization curve:

> How do I reconcile that and what could be causing this large discrepancy between what we expect the optimal surface azimuth to be and what our code is outputting?

Put another way -- the code you originally posted is modeling a hypothetical system where solar position and irradiance are timeshifted with respect to each other.  That's a big departure from reality, and so real-life expectations don't apply to your model.

Kevin

[Casey McGonigle]

Hi Kevin,

I just reran my code with your timezone correction (apologies for responding before rerunning) and you're right. That change made the maximum energy output occur near (within 15° of) due South (ie. 180°). 

That was a huge step for me. I really appreciate your continued help.

-Casey