So, alot is made of the functionality of PVSYST as a standard tool to model PV systems. No arguments here. That said, it does not have a programming interface. Thankfully there are some great alternatives out there developed through Sandia National Laboratories, particularly the PV-LIB Toolbox (https://pvpmc.sandia.gov/applications/pv_lib-toolbox/), which contains some really excellent functions for modeling solar position, irradiance, and PV systems. The tool was initially created in MATLAB and subsequently migrated to Python (pvlib: http://pvlib-python.readthedocs.io/en/latest/).
Goals and Objectives
Model Plane of Array Irradiance for a fictitious fixed tilt and tracking PV system scenarios using pvlib scripted through R via the reticulate package. I realize it would probably be a bit smoother to use Python directly, but to be honest I love R and also wanted to test out a hybrid language workflow using a combination of R and Python available through reticulate package https://rstudio.github.io/reticulate/. Lastly, I wanted to see how much solar irradiance I could expect from an array mounted on my roof at home.
Let’s walk through the steps and see how this all works.
Load R Libraries
First, let’s load the r libraries and define a timestamp convention.
library(tidyverse) # data analysis and visualization library(reticulate) # loading python libraries in R library(dygraphs) # timeseries visualization library(timetk) # tidy timeseries tools library(leaflet) # mapping library(knitr) # tables library(viridis) # color sheme library(ggridges) # distribution plot Sys.setenv(TZ='UTC')
load NSRDB data
To do this, we’ll be querying data from the National Renewable Energy Laboratory (NREL) Physical Solar Model v3 https://nsrdb.nrel.gov/current-version via their developer API (https://nsrdb.nrel.gov/api-instructions). I’m going to use my own home in MA as an example.
leaflet(data = data.frame(lat=lat,lon=lon,name='home')) %>% addProviderTiles('Esri.WorldImagery') %>% addMarkers(lng = ~lon,lat=~lat) %>% addPopups(lng=~lon,lat=~lat,popup=~name)
query the data
You’ll need your own API key to run this, which you can get for free through NREL’s website https://nsrdb.nrel.gov/download-instructions. For this example I’m just going to grab a random year of data.
# information on dataset: https://nsrdb.nrel.gov/current-version attributes = c('air_temperature,clearsky_dhi,clearsky_dni,clearsky_ghi,cloud_type,dew_point,dhi,dni,fill_flag,ghi,relative_humidity,solar_zenith_angle,surface_albedo,surface_pressure,total_precipitable_water,wind_direction,wind_speed') leap_year = 'false' utc='true' year<-2010 interval='60' local_tzone<-'America/New_York' url = paste('http://developer.nrel.gov/api/solar/nsrdb_psm3_download.csv?wkt=POINT(',lon,'%20',lat,')&names=',year,'&leap_day=false&interval=',interval,'&utc=',utc, '&full_name=',name,'&email=',email,'&affiliation=',affiliation,'&mailing_list=',mailing_list,'&reason=',reason,'&api_key=',api_key,'&attributes=',attributes,sep='') #download data as a dataframe df<-read_csv(url,skip=2,col_types = cols()) df$timestamp_utc=as.POSIXct(paste(df$Year,df$Month,df$Day,df$Hour,df$Minute,'00',sep='-'),format='%Y-%m-%d-%H-%M-%S',tz='UTC') # tidy names names(df)<-gsub(pattern = ' ',replacement = '',x = names(df))
here’s what the data looks like.
## # A tibble: 6 x 23 ## Year Month Day Hour Minute Temperature ClearskyDHI ClearskyDNI ## <int> <int> <int> <int> <int> <int> <int> <int> ## 1 2010 1 1 0 30 -3 0 0 ## 2 2010 1 1 1 30 -3 0 0 ## 3 2010 1 1 2 30 -3 0 0 ## 4 2010 1 1 3 30 -3 0 0 ## 5 2010 1 1 4 30 -3 0 0 ## 6 2010 1 1 5 30 -3 0 0 ## # ... with 15 more variables: ClearskyGHI <int>, CloudType <int>, ## # DewPoint <int>, DHI <int>, DNI <int>, FillFlag <int>, GHI <int>, ## # RelativeHumidity <dbl>, SolarZenithAngle <dbl>, SurfaceAlbedo <dbl>, ## # Pressure <dbl>, PrecipitableWater <dbl>, WindDirection <dbl>, ## # WindSpeed <dbl>, timestamp_utc <dttm>
Load Python Libraries
Next, I’ll load the python libraries that I’ll need for this example using reticulate, including pvlib. Note, the “convert=F” option for PVLIB in which I’d like to keep the output from each pvlib function in it’s native python format rather than auto-converting it to its R equivalent.
pvlib<-import("pvlib",convert=F) np<-import("numpy") pd<-import("pandas")
PV system setup
I’ve got pretty decent sun exposure on my roof (170 degrees azimuth, so almost due south) and my roof tilt is almost in lock stop with my latitude at 45 degrees from horizon. So let’s model a fixed tilt system on my roof. Also, since this is somewhat of a tutorial illustrating functionality, for kicks let’s model a tracking system that I’ll put on my front lawn. I’m sure the neighbors would LOVE it.
location = pvlib$location$Location(latitude = lat,longitude = lon,tz = 'UTC',altitude = 100,name = 'home') pvsystem_fixed = pvlib$pvsystem$PVSystem(surface_tilt=45,surface_azimuth=170,albedo=.2,name='home fixed tilt') pvsystem_tracking=pvlib$tracking$SingleAxisTracker(axis_tilt = 0,axis_azimuth = 170,max_angle = 52,backtrack = TRUE,gcr = .5) python_timestamp<-pd$DatetimeIndex(start='2010-01-01 00:30', end='2010-12-31 23:30', freq='1h',tz=location$tz)
Solar Position and Atmospheric Modeling
To convert the three components of irradiance provided from NSRDB (GHI,DHI and DNI) to the plane of array (POA) of my rooftop, I first need to know where the sun is in the sky at any given point in the year, some default assumptions about the amount of extraterrestrial radiation and the atmosphere it is travelling through and ultimately the angle of incidence (AOI) on the collector plane of my rooftop solar modules.
python_timestamp<-pd$DatetimeIndex(start=as.character(min(df$timestamp_utc)),end=as.character(max(df$timestamp_utc)), freq='1h',tz=location$tz) solpos = location$get_solarposition(python_timestamp) dni_extra = pvlib$irradiance$extraradiation(python_timestamp) airmass = pvlib$atmosphere$relativeairmass(solpos$apparent_zenith) aoi = pvlib$irradiance$aoi(pvsystem_fixed$surface_tilt,pvsystem_fixed$surface_azimuth,solpos$apparent_zenith, solpos$azimuth) # convert output back to R-equivalent and append to data frame. df<-cbind(df,py_to_r(solpos),dni_extra=py_to_r(dni_extra$values),airmass=py_to_r(airmass),aoi=py_to_r(aoi))
Plane of Array Transposition
Calculating POA from GHI, DNI and DHI is known as transposition and standard transposition model is the ‘perez’ model which thankfully pvlib has a handy function for. R. Perez, P. Ineichen, R. Seals, and J. Michalsky, “Modeling daylight availability and irradiance components from direct and global irradiance,” Sol. Energy, vol. 44, pp. 271-289, 1990.
Fixed Tilt Example
We previously calculated the AOI on the module surface for all periods in our dataset as well as other relevant parameters, so we have everything we need to model POA for the fixed tilt scenario.
total_irrad = pvlib$irradiance$total_irrad(surface_tilt=pvsystem_fixed$surface_tilt, surface_azimuth=pvsystem_fixed$surface_azimuth, apparent_zenith= solpos$apparent_zenith, albedo=pvsystem_fixed$albedo, azimuth= solpos$azimuth, dni= pd$Series(df$DNI,python_timestamp), # pvlib requires a data series with an index ghi= pd$Series(df$GHI,python_timestamp), dhi= pd$Series(df$DHI,python_timestamp), dni_extra= dni_extra, model='perez', solar_zenith= solpos$zenith, airmass= airmass) df<-cbind(df,py_to_r(total_irrad)) df<-replace_na(df,list(poa_global=0,poa_direct=0,poa_diffuse=0,poa_sky_diffuse=0,poa_ground_diffuse=0)) %>% dplyr::rename(poa_global_fixed=poa_global)
Single Axis Tracking Example
The single axis tracking example is a bit more complicated than the fixed tilt as we’ll first need to model how the tracker follows the sun such that we can know what the AOI will be on the trackers at any given time. But once we have that, the
tracker_data = pvlib$tracking$singleaxis(solpos$apparent_zenith, solpos$azimuth, axis_tilt=pvsystem_tracking$axis_tilt, axis_azimuth=pvsystem_tracking$axis_azimuth, max_angle=52, backtrack=TRUE, gcr=.5) perez_diffuse = pvlib$irradiance$perez(tracker_data$surface_tilt, tracker_data$surface_azimuth, pd$Series(df$DHI,python_timestamp), pd$Series(df$DNI, python_timestamp), dni_extra, solpos$apparent_zenith, solpos$azimuth, airmass=airmass) ground_irrad = pvlib$irradiance$grounddiffuse(tracker_data$surface_tilt, pd$Series(df$GHI,python_timestamp), albedo=.2) cosd<-pvlib$tools$cosd(tracker_data$aoi) poa_global_tracking = py_to_r(cosd)*df$DNI + py_to_r(perez_diffuse) + py_to_r(ground_irrad) df<-cbind(df,poa_global_tracking)%>% replace_na(list(poa_global_tracking=0))
Results and Comparisons
# convert to local time df<-mutate(df,timestamp_local=timestamp_utc) attributes(df$timestamp_local)$tzone<-local_tzone df %>% select(timestamp_local,poa_global_fixed,poa_global_tracking,GHI,DNI,DHI) %>% tk_xts(date_var=timestamp_local,silent = T) %>% dygraph() %>% dyRoller() %>% dyRangeSelector() %>% dyOptions(useDataTimezone = TRUE)
Example Clear Day
You can see on a clear day that the POA irradiance at my latitude in the northern hemisphere is typically a bit higher than it’s GHI equivalent. You can also see a slight skew with higher POA values in the morning for the fixed tilt scenario than for the afternoon and vice versa for the tracker scenario owed to thenon-ideal azimuthal orientation of my home to due south (170 degrees, instead of 180).
df %>% filter(timestamp_local>=as.POSIXct('2010-09-21',tz=local_tzone) & timestamp_local<as.POSIXct('2010-09-22',tz=local_tzone)) %>% select(timestamp_local,poa_global_fixed,poa_global_tracking,GHI,DNI,DHI) %>% gather(type,`W/m^2`,poa_global_fixed:DHI) %>% ggplot(aes(timestamp_local,`W/m^2`,color=type))+ geom_line(lwd=1.2)+ theme_minimal()+ scale_color_viridis(discrete = T)+ labs(title='Fixed Tilt and Tracking Sample Clear Day')
aggregates<-df %>% group_by(Month) %>% summarise(GHI=sum(GHI)*.001, POA_fixed=sum(poa_global_fixed)*.001, POA_tracking=sum(poa_global_tracking)*.001 ) kable(aggregates,digits = 1,align = 'c',caption='Monthly kWh/m^2')
kable(colSums(select(aggregates,-Month)),digits=1,align='c',caption = 'Annual kWh/m^2')
Regardless of my PV system choice (fixed tilt vs. tracking) I would gain extra insolation over the course of the year as compared to just placing the panels on the ground (GHI). Looking at a distribution of irradiance by month and mounting type, it becomes clear that the tracker mount would pay the most dividends over the summer months (higher irradiance values) and the fixed tilt mounting would win out in the winter months with the crossover point occuring during spring/fall.
df %>% select(timestamp_local,poa_global_fixed,poa_global_tracking,GHI) %>% gather(type,`W/m^2`,poa_global_fixed:GHI) %>% mutate(month=factor(format(timestamp_local,'%m'))) %>% filter(`W/m^2`>0) %>% ggplot(aes(x=`W/m^2`,y=month,fill=type))+ geom_density_ridges(alpha=.25,color='black')+ theme_ridges(font_size = 12)+ scale_fill_viridis(discrete = T)+ theme(legend.position = 'bottom')+ labs(title='Irradiance Distribution by Month',caption='Filtered for W/m^2 > 0')
OK, so while ultimately mounting my fictitious PV system on a tracker rack would give me ~6% more POA insolation over the course of the year, the juice is probably not worth the squeeze and I should stick with just mounting the panels on my roof like everyone else. In the next iteration I’ll try and walk through modeling energy capture of my fictitous fixed tilt and tracking PV systems in pvlib.