Model Chain example PVLIB - don't trust in 1 axis tracking AC output

Question

I'm trying to use PVLIB to estimate output power for a PV System installed in the west of my country.

As an example I've got 2 days of hourly GHI, 2m Temperature and 10m wind speed from MERRA2 reanalysis.

I want to estimate how much power a fixed PV System or 1 axis tracking system would generate using the forementioned dataset, and ModelChain function from PVLIB. I first estimate DNI and DHI from GHI data using DISC model to obtain DNI and then DHI is the difference between GHI and DNI*cos(Z)

a) First behaviour I am not completely sure if it is Ok. Here is the plot of GHI, DNI , DHI, T2m and Wind Speed. It seems that DNI is shifted with its maximum occurring 1 hour before GHI maximum.

Weather Figure

After preparing irradiance data I calculated AC using Model Chain, specifying the fixed PV System and 1 axis single tracking system. The thing is that I don't trust in the AC output for a 1-single axis system. I expected a plateau shape of AC output and i found a kind of weird behaviour.

Here is the otuput values of power generation i expected to see:

Expectation

And here is the estimated output by PVLIB

Reality

I hope someone can help me to find the error on my proccedure.

Here is the code:

# =============================================================================
# Example of using MERRA2 data and PVLIB
# =============================================================================
import numpy as np
import pandas as pd
import pandas as pd
import matplotlib.pyplot as plt
import pvlib
from pvlib.pvsystem import PVSystem
from pvlib.location import Location
from pvlib.modelchain import ModelChain

# =============================================================================
# 1) Create small data set extracted from MERRA
# =============================================================================
GHI             = np.array([0,0,0,0,0,0,0,0,0,10.8,148.8,361,583,791.5,998.5,1105.5,1146.5,1118.5,1023.5,
                            860.2,650.2,377.1,165.1,16,0,0,0,0,0,0,0,0,0,11.3,166.2,395.8,624.5,827,986,
                            1065.5,1079,1025.5,941.5,777,581.5,378.9,156.2,20.6,0,0,0,0])
temp_air        = np.array([21.5,20.5,19.7,19.6,18.8,17.9,17.1,16.5,16.2,16.2,17,21.3,24.7,26.9,28.8,30.5,
                            31.6,32.4,33,33.3,32.9,32,30.6,28.7,25.4,23.9,22.6,21.2,20.3,19.9,19.5,19.1,18.4,
                            17.7,18.3,23,25.1,27.3,29.5,31.2,32.1,32.6,32.6,32.5,31.8,30.7,29.6,28.1,24.6,22.9,
                            22.3,23.2])
wind_speed      = np.array([3.1,2.7,2.5,2.6,2.8,3,3,3,2.8,2.5,2.1,1,2.2,3.7,4.8,5.6,6.1,6.4,6.5,6.6,6.3,5.8,5.3,
                            3.7,3.9,4,3.6,3.4,3.4,3,2.6,2.3,2.1,2,2.2,2.7,3.2,4.3,5.1,5.6,5.7,5.8,5.8,5.7,5.4,4.8,
                            4.4,3.1,2.7,2.3,1.1,0.6])
local_timestamp = pd.DatetimeIndex(start='1979-12-31 21:00', end='1980-01-03 00:00', freq='1h',tz='America/Argentina/Buenos_Aires')
d               = {'ghi':GHI,'temp_air':temp_air,'wind_speed':wind_speed}
data            = pd.DataFrame(data=d)
data.index      = local_timestamp

lat             = -31.983   
lon             = -68.530
location        = Location(latitude  = lat, 
                           longitude = lon,
                           tz        = 'America/Argentina/Buenos_Aires',
                            altitude = 601)

# =============================================================================
# 2) SOLAR POSITION AND ATMOSPHERIC MODELING
# =============================================================================
solpos          = pvlib.solarposition.get_solarposition(time      = local_timestamp, 
                                                        latitude  = lat,
                                                        longitude = lon,
                                                        altitude  = 601)

# DNI and DHI calculation from GHI data
DNI             = pvlib.irradiance.disc(ghi             = data.ghi, 
                                        solar_zenith    = solpos.zenith, 
                                        datetime_or_doy = local_timestamp)
DHI             = data.ghi - DNI.dni*np.cos(np.radians(solpos.zenith.values))

d               = {'ghi': data.ghi,'dni': DNI.dni,'dhi': DHI,'temp_air':data.temp_air,'wind_speed':data.wind_speed }
weather         = pd.DataFrame(data=d)
plt.plot(weather)

# =============================================================================
# 3) SYSTEM SPECIFICATIONS
# =============================================================================
# load some module and inverter specifications
sandia_modules    = pvlib.pvsystem.retrieve_sam('SandiaMod')
cec_inverters     = pvlib.pvsystem.retrieve_sam('cecinverter')
sandia_module     = sandia_modules['Canadian_Solar_CS5P_220M___2009_']
cec_inverter      = cec_inverters['Power_Electronics__FS2400CU15__645V__645V__CEC_2018_']

# Fixed system with tilt=abs(lat)-10
f_system          = PVSystem(      surface_tilt                                = abs(lat)-10, 
                                   surface_azimuth                             = 0,
                                   module                                      = sandia_module,
                                   inverter                                    = cec_inverter,
                                   module_parameters                           = sandia_module,
                                   inverter_parameters                         = cec_inverter,
                                   albedo                                      = 0.20,
                                   modules_per_string                          = 100,
                                   strings_per_inverter                        = 100)
# 1 axis tracking system
t_system          = pvlib.tracking.SingleAxisTracker(axis_tilt                 = 0, #abs(-33.5)-10
                                                     axis_azimuth              = 0,
                                                     max_angle                 = 52,
                                                     backtrack                 = True,
                                                     module                    = sandia_module,
                                                     inverter                  = cec_inverter,                                                       
                                                     module_parameters         = sandia_module,
                                                     inverter_parameters       = cec_inverter,
                                                     name                      = 'tracking',
                                                     gcr                       = .3,
                                                     modules_per_string        = 100,
                                                     strings_per_inverter      = 100)

# =============================================================================
# 4) MODEL CHAIN USING ALL THE SPECIFICATIONS for a fixed and 1 axis tracking systems
# =============================================================================
mc_f       = ModelChain(f_system, location)
mc_t       = ModelChain(t_system, location)


# Next, we run a model with some simple weather data.
mc_f.run_model(times=weather.index, weather=weather)
mc_t.run_model(times=weather.index, weather=weather)

# =============================================================================
# 5) Get only AC output form a fixed and 1 axis tracking systems and assign 
# 0 values to each NaN
# =============================================================================
d              = {'fixed':mc_f.ac,'tracking':mc_t.ac}
AC             = pd.DataFrame(data=d)
i              = np.isnan(AC.tracking)
AC.tracking[i] = 0
i              = np.isnan(AC.fixed)
AC.fixed[i]    = 0

plt.plot(AC)

I hope anyone could help me with the intepretation of the results and debugging of the code.

Thanks a lot!

Answer 1

I suspect your issue is due to the way the hourly data is treated. Be sure that you're consistent with the interval labeling (beginning/end) and treatment of instantaneous vs. average data. One likely cause is using hourly average GHI data to derive DNI data. pvlib.solarposition.get_solarposition returns the solar position at the instants in time that are passed to it. So you're mixing up hourly average GHI values with instantaneous solar position values when you use pvlib.irradiance.disc to calculate DNI and when you calculate DHI. Shifting your time index by 30 minutes will reduce, but not eliminate, the error. Another approach is to resample the input data to be of 1-5 minute resolution.