project-valuation.Rmd

Energy Storage Technology Valuation
========================================================

An R version of the method described in [Energy Storage Technology Valuation Primer](http://www.epri.com/abstracts/Pages/ProductAbstract.aspx?ProductId=000000000001008810) from EPRI.

## Deterministic Model

NPV is calculated in the `npv()` function in `func.r`

Source [here](http://tolstoy.newcastle.edu.au/R/help/05/12/16765.html)

   >t = c(-8000, 100, 100, 100, 2000, 3000, 4000, 5000)

   >r = 0.1

   >npv(r, t)

   >[1] 301.1624


```{r}
require(plyr)
require(pastecs)
require(ggplot2)
require(reshape)
require(FME)
require(xtable)
source('~/Code/energy_projects/func.r')
options(scipen=5)


percentify <- function(x, r=1, keep.zeros=FALSE){
         ifelse (as.numeric(x) != 0 | keep.zeros==TRUE,  paste( format(round ( as.numeric(x) * 100, r) , nsmall=r) , '%', sep='') 
          	,
        		'')       		
}


inputs <-data.frame(
  default=c(
  80,20,5,5,10,250,100000,2
  )
)
row.names(inputs) <- c(
    'on.peak.rate',
    'off.peak.rate',
    'demand.charge',
    'capacity',
    'load.shift.duration',
    'load.shifting.periods',
    'cost.per.volt.sag',
    'volt.sag.count'
    )
inputs$var.name <- row.names(inputs)

# From http://www.tybecenergy.com/pricehistory/pjm_settle.php
on.peak.rate <- 44.40 # $
off.peak.rate <- 31.13 # $


b.total <- calc.benefits(inputs)

sprintf("$%.f", b.total)
```

Now exercise the model for sensitivity analysis

```{r figure.3-1}

percent(exercise(inputs, 'capacity',1.10) / b.total - 1)
percent(exercise(inputs, 'on.peak.rate',1.10) / b.total - 1)

df.delta10 <- ddply(inputs, .(var.name), function(x){ 
    c( delta10=(exercise(inputs, x$var.name, 1.10 )-b.total)/b.total
       )
     }
  )

df.delta10 <- transform(df.delta10, 
                          var.name = reorder(var.name, delta10))

ggplot(df.delta10, aes(var.name, delta10)) + geom_bar(stat="identity") + coord_flip() + scale_y_continuous(labels=percent, limits=c(-0.05,0.1))

```

Now examine the change in benefits over a range of different inputs

### Two-Dimensional Data Table, 5MW Energy Storage Example, Total Annual Benefits


```{r table.3-1}
on.peak.rate.alt <- c(50,60,70,80,90,100)
load.shift.duration.alt <- c(8,9,10,11) 

alt.vars <- expand.grid("on.peak.rate.alt"=on.peak.rate.alt, "load.shift.duration.alt"=load.shift.duration.alt)

for(i in 1:nrow(alt.vars)){
  alt.inputs <- inputs
  alt.inputs['on.peak.rate','default']  <- alt.vars$on.peak.rate.alt[i]
  alt.inputs['load.shift.duration','default']  <- alt.vars$load.shift.duration.alt[i]
  alt.vars$benefit[i] <- calc.benefits(alt.inputs) 
  }
alt.vars
cast(alt.vars, on.peak.rate.alt~load.shift.duration.alt)
```


### Two-Dimensional Data Table, 10MW Energy Storage Example, Total Annual Benefits



```{r table.3-2}
alt.inputs <- inputs
alt.inputs['capacity','default'] <- 10
on.peak.rate.alt <- c(50,60,70,80,90,100)
load.shift.duration.alt <- c(8,9,10,11) 

alt.vars <- expand.grid("on.peak.rate.alt"=on.peak.rate.alt, "load.shift.duration.alt"=load.shift.duration.alt)

for(i in 1:nrow(alt.vars)){
  alt.inputs['on.peak.rate','default']  <- alt.vars$on.peak.rate.alt[i]
  alt.inputs['load.shift.duration','default']  <- alt.vars$load.shift.duration.alt[i]
  alt.vars$benefit[i] <- calc.benefits(alt.inputs) 
  }
alt.vars
cast(alt.vars, on.peak.rate.alt~load.shift.duration.alt)
```

## Estimating Uncertainty

### Example 2 -- Calculating Total Uncertainty for Energy Storage Annual Savings

#### Table 3-3

#### Nominal Values for Energy Storage Example

```{r table.3.3,results='asis'}
print(xtable(inputs), type='html') 
```


```{r uncertainty}
inputs$uncertainty.pct <- c(0.2,0.2,0.2,0.03,0.1,0.1,0.2,.5)
inputs$uncertainty.abs <- inputs$default * inputs$uncertainty.pct
inputs$low <- inputs$default - inputs$uncertainty.abs
inputs$high <- inputs$default + inputs$uncertainty.abs

```

#### Table 3-4

#### Uncertainty Estimates for Energy Storage Example
```{r table.3.4, results='asis'} 
print(xtable(inputs), type='html') 
``` 


#### Eq. 3-3

Vary each input by 1% and measure the change in output

```{r kline.mclintock}

alt.inputs <- inputs[,c('default','var.name','uncertainty.abs')]



varied <- ddply(alt.inputs, .(var.name, uncertainty.abs), function(x) c(value=x$default, sensitivity=calc.sensitivity(alt.inputs, x$var.name, 1.01) ))
varied$km <- (varied$uncertainty.abs * varied$sensitivity)^2

varied

varied <- ddply(varied, .(),  transform, contribution=km/sum(km) )
varied$contribution.pct <- percentify(varied$contribution)

```

#### Table 3-5

#### Uncertainty Analysis, Energy Storage Example

```{r table.3.5, results='asis'}
print(xtable(varied), type='html') 

cat('Total Uncertainty: ' , sprintf("$%.0f", sum(varied$km)^.5) )

```

#### Figure 3-3

#### Uncertainty Analysis, Energy Storage System
```{r figure.3.3}
ggplot(varied, aes(var.name, contribution)) + geom_bar(stat="identity") + coord_flip() + scale_y_continuous(labels=percent, limits=c(-0.05,1))

```

Run the scenario again, with the assumption of 10% uncertainty for On-Peak Energy rate.

#### Table 3-6

#### Uncertainty Analysis Energy Storage Example, Reducing the On-Peak Energy Rate Uncertainty Assumption from +/- 20% to +/- 10%

```{r table.3.6, results='asis'}
inputs.b <- inputs
inputs.b['on.peak.rate','uncertainty.pct'] <- 0.1
inputs.b$uncertainty.abs <- inputs.b$default * inputs.b$uncertainty.pct

varied.b <- ddply(inputs.b, .(var.name, uncertainty.abs), function(x) c(value=x$default, sensitivity=calc.sensitivity(inputs.b, x$var.name, 1.01) ))

varied.b$km <- (varied.b$uncertainty.abs * varied.b$sensitivity)^2

varied.b <- ddply(varied.b, .(),  transform, contribution=km/sum(km) )
varied.b$contribution.pct <- percentify(varied.b$contribution)
print(xtable(varied.b), type='html') 

cat('Total Uncertainty: ' , sprintf("$%.0f", sum(varied.b$km)^.5) )

```

## 4

## Monte Carlo Modeling

Assign standard deviations to the inputs
```{r monte.carlo}
inputs$sd <- (inputs$default * inputs$uncertainty.pct)/2

```

#### Figure 4-1

#### Assumed Probability Distribution, On-Peak Energy Rate

```{r fig.4.1}
x=seq(56,104,length=200)
y=dnorm(x,mean=80,sd=8)
plot(x,y,type="l",lwd=2,col="red")
```
#### Figure 4-2

#### Assumed Probability Distribution, System Size

```{r fig.4.2}
x=seq(4.775,5.225,length=200)
y=dnorm(x,mean=5,sd=0.075)
plot(x,y,type="l",lwd=2,col="red")
```

#### Table 4-1

#### Monte Carlo Model Inputs, Energy Storage Example

```{r table.4.1, results='asis'}
print(xtable(inputs[,c('var.name','default','uncertainty.abs')] ), type='html')
```

Run multiple iterations across the input distributions

#### Figure 4-3
#### Monte Carlo Simulation, Frequency Distribution for Total Annual Benefits
Corresponds to Figures 4-3 and 4-4 in the source document

```{r mc.inputs}

mc.results <- mc.benefits(inputs, 5000)

t.test(mc.results, conf.level=0.95)

hist(mc.results, freq=FALSE, breaks=50)
data.frame(stat.desc(mc.results))

# Uncertainty:
cat(
    '$', format(stat.desc(mc.results)['std.dev'] * 2 , nsmall = 2, big.mark=',')
    ,sep='')


varied.mc <- ddply(inputs, .(var.name, uncertainty.abs), function(x) c(value=x$default, sensitivity=calc.sensitivity.mc(inputs, x$var.name, 1.01, 5000) ))

varied.mc$km <- (varied.mc$uncertainty.abs * varied.mc$sensitivity)^2
varied.mc <- ddply(varied.mc, .(),  transform, contribution=km/sum(km) )
varied.mc$contribution.pct <- percentify(varied.mc$contribution)

varied.mc[,c('var.name','contribution.pct')]

```