---
title: "Power Predictive Model Comparison"
author: "Tyler Burleyson"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{model_comparison}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

```{r include = FALSE}
library(PVplr)
library(knitr)
library(broom)
library(purrr)
library(tidyr)

# using the package's example data
var_list <- PVplr::plr_build_var_list(time_var = "timestamp", 
                               power_var = "power", 
                               irrad_var = "g_poa", 
                               temp_var = "mod_temp", 
                               wind_var = NA)

test_dfc <- plr_cleaning(test_df, var_list)
```
The Solar Durability and Lifetime Extension (SDLE) center currently uses several power predictive models in analysis. Power data cannot be easily compared between systems due to variations in conditions, so these models allow meaningful comparisons by predicting power as a function of weather. Each model uses unique formulas for this function, with the end goal of producing predicted power values that are independent of weather.

These methods are discussed in-depth at an academic level in the SDLE paper, "Performance Loss Rate Consistency and Uncertainty across Multiple Methods and Filtering Criteria". This vignette will focus on practical differences between the models that you should know in order to make the most use of them. Before beginning, please notice the following result from the study:

> "The PLR with the lowest uncertainty is not necessarily the most accurate PLR. It has been observed that PLR results on the same system can exist outside of each others error bars based on filtering and predictive model. While one method may give consistent results it is important to compare between methods to ensure that the result is an accurate representation on the actual performance of the system."

All of the methods share several qualities. Options exist to group data for regression by day, week, or month. Data is grouped by these time periods and linear models are created for each group. Predicted values (either generated by the method or supplied by the user) are given to these models. If you are interested in correcting power at a particular set of standard conditions, use this option. Finally, fitted values and standard deviations are attached to the output in the .fitted and .se.fit columns (in all cases other than the 6k model; see its entry below). For more in-depth analysis of uncertainty in PLR calculations, see the bootstrap functions.

## XbX Method ##
The simplest model is the plr_xbx_model. It is data driven, meaning it does not take in outside information for the creation of its linear regression models. This offers a great deal of flexibility - XbX works well with either POA or GHI irradiance, and either module or ambient temperature readings. It also lends itself to change point analysis.

### Equation ###
$P_{pred.} = \beta_{0} + \beta_{1} G + \beta_{2} T + \epsilon$

This model uses a simplistic, linear equation relating power to irradiance (G) and temperature (T) with constant coefficients and an error term ($\epsilon$). It is this simplicity that enables its use for change point PLR, as well as the options to use POA or GHI irradiance and module or air temperature.

### Implementation ###
```{r}
test_xbx_wbw_res <- plr_xbx_model(test_dfc, var_list, by = "week", data_cutoff = 30, predict_data = NULL)

knitr::kable(test_xbx_wbw_res[1:5, ], caption = "XbX Model: Week-by-Week Implementation")
```

Due to its data-driven nature, this method will typically show the strongest seasonality of those offered in this package. This is likely due to extrapolation done when predicting representative conditions affects seasonality in the data. Higher irradiance filters will reduce this seasonality, but also reduce the amount of data being used.

## XbX + UTC Method ##

Modelling temperature can be challenging due to low variation on a daily basis and high variation across seasons. Models like the XbX are thus prone to extrapolation, as mentioned above. It is possible to introduce a Universal Temperature Correction (UTC) to solve this issue by converting measured temperature to representative temperature, resulting in a corrected power. Then, this corrected power can be used to model based off of irradiance, with temperature as a variable removed.

#### Determining UTC ####

Data is subset to a narrow irradiance range so that the effect of temperature on generated power can be studied in isolation. A linear regression is performed between temperature and power, giving the change in power per degree celsius. This is then converted to percent power loss per degree celcius ($\gamma_{T}$).

### Equation ###
$P_{cor} = \frac{P_{obs}}{1 + \gamma_{T} (T_{obs} - T_{rep}) (\frac{G_{obs}}{G_{rep}})}$

$P_{cor} = \beta_{0} + \beta_{1} G + \epsilon$

Note the similarity in the second equation to that in the XbX method; 

### Implementation ###
```{r}
test_xbxutc_wbw_res <- plr_xbx_utc_model(test_dfc, var_list, by = "week", data_cutoff = 30, predict_data = NULL, ref_irrad = 900, irrad_range = 10)

knitr::kable(test_xbxutc_wbw_res[1:5, ], caption = "XbX + UTC: Week-by-Week Implementation")
```


## PVUSA Method ##

This method is well known in the industry. It's physics-based, meaning it's formula is derived from Ohm's law, $P = I * V$. The model assumes current to be a function of plane-of-array (POA) irradiance, and voltage to be a function of irradiance and module temperature, predicted from the ambient temperature and wind speed (when available). The equation, then, becomes the following:

### Equation ###
$P = G_{POA} (\beta_{0} + \beta_{1} G + \beta_{2} T_{amb} + \beta_{3} W)$

It is worth noting that the method requires POA irradiance and ambient temperature specifically - it lacks that flexibility found in the XbX method.

### Implementation ###
```{r}
test_pvusa_wbw_res <- plr_pvusa_model(test_dfc, var_list, by = "week", data_cutoff = 30, predict_data = NULL)

knitr::kable(test_pvusa_wbw_res[1:5, ], caption = "PVUSA: Week-by-Week Implementation")
```

See [this paper](https://www.nrel.gov/docs/fy09osti/45376.pdf) as well for an in-depth look at uncertainties in the model.

## 6k Method ##

The most complicated of the SDLE models, the 6k model makes use of nameplate power and six constants for modeling. Similarly to the PVUSA model, it is an implementation of Ohm's law, $P = I * V$. Here, irradiance as a fraction of standard irradiance represents current. Voltage is modeled by a complex statement involving nameplate power, the fraction of irradiance over standard irradiance, and the difference between module temperature and reference temperature. 

### Equation ###

$G^{\prime} = G_{POA} / G_{STC}$

$T^{\prime} = T_{mod} / T_{STC}$

$P = G^{\prime} * (P_{NP} + k_1 \log(G^{\prime}) + k_2 \log(G^{\prime})^2 + k_3 T^{\prime} + k_4 T^{\prime} log(G^{\prime}) + k_5 T^{\prime} (G^{\prime})^2 + k_6 (T^{\prime})^2)$


The $G^{\prime}$ and $T^{\prime}$ terms are defined for the purposes of temperature and irradiance correction. 

### Implementation ###
```{r}
test_6k_wbw_res <- plr_6k_model(test_dfc, var_list, nameplate_power = 230, by = "week", data_cutoff = 30, predict_data = NULL)

knitr::kable(test_6k_wbw_res[1:5, ], caption = "6k: Week-by-Week Implementation")
```

Another important distinction between the 6k and other models is its linear regression scheme. While the others use standard least-squares linear regression through the base R 'lm()' function, the 6k implements the Levenberge-Marquardt method through minpack.lm (specifically, the 'nlsLM()' function). This causes several disparities. Firstly, the returned dataframe lacks .se.fit and sigma columns. Secondly, the model attempts to fit the data by adjusting its six constants, so the fit improves dramatically with a greater quantity of data.

The 6k model has the greatest uncertainty of any of the models as well as the most extreme outliers, and so it is recommended that data given to it is filtered as little as possible. In particular, an SDLE publication found that the lowest deviation was achieved when irradiance is filtered only below 0 watts per meter squared.

One final note on the 6k model: since it uses nameplate power in its fitting, if it is given standard conditions (STC) for predict_data, it will simply return the nameplate power value at all points in time. Be sure to either let predicted data be calculated by the model, or supply values specific to the actual location of the system.

#### Potential Errors ####

Note that the model can fail to fit the data properly, most often with shorter time scales such as day-by-day modeling. In this case, it returns a rather esoteric error message. Usually there are many warning messages as well: here they are suppressed for readability.s
```{r, error = TRUE, warning = FALSE}
# This code should create an error
test6k_dbd_res <- plr_6k_model(test_dfc, var_list, nameplate_power = 230, by = "day", data_cutoff = 30, predict_data = NULL)
```

Sources
------------

Alan J. Curran, et al. Performance Loss Rate Consistency and Uncertainty Across Multiple Methods and Filtering Criteria. IEEE PVSC 46, Chicago, IL, USA.

A. J. Curran, R. Zhang, Y. Hu, R. Haddadian, J. L. Braid, T. J. Peshek,
and R. H. French, “Determining the Power Rate of Change of 353 Solar
Plant Inverters Using a Month-by-Month Analysis and Common Data
Science Applications to Power Time Series,” CWRU, Aug. 2017.

D. L. King, J. A. Kratochvil, and W. E. Boyson, “Field
experience with a new performance characterization procedure for
photovoltaic arrays,” Sandia National Laboratories, Tech. Rep., 1997. 
Available: http://www.osti.gov/energycitations/servlets/purl/
629484-omiEsH/webviewable/

T. Huld, G. Friesen, A. Skoczek, R. P. Kenny, T. Sample, M. Field,
and E. D. Dunlop, “A power-rating model for crystalline silicon pv
modules,” Solar Energy Materials and Solar Cells, vol. 95, no. 12, pp.
3359–3369, 2011.