NPV_and_LCOE_and_IRR.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Implements the functions for the net present value and levelized cost of energy\n",
    "\n",
    "- NPV and NPV2 differ in how wholesale price is assumed to fluctuate. \n",
    "- In NPV it is an uncertain average number over the project lifetime. \n",
    "- In NPV2 it varies within given window during every year of the project. In NPV2 one could replace this randomness with the average wholesale price without large error.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.optimize as opt\n",
    "#from scipy.optimize import brentq  #root finding"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def NPV(construction_period,n,r,cf,It0,OMkwh,OMfixed,wholesale_price):\n",
    "    prod=cf*1*365*24\n",
    "    constructiont=int(round(construction_period))\n",
    "    It=np.zeros(n+constructiont)\n",
    "    #It[0]=It0\n",
    "    Mt=np.ones(n+constructiont)*(OMkwh*prod+OMfixed)\n",
    "    Et=np.ones(n+constructiont)*prod  #kwH\n",
    "    if constructiont>0:  #if construction time zero, capex is the overnight cost         \n",
    "        for x in range(constructiont):\n",
    "            It[x]=It0/constructiont  #spread the investment evenly over the construction period\n",
    "            Mt[x]=0   #No costs or production during construction\n",
    "            Et[x]=0\n",
    "    else:\n",
    "        It[0]=It0    \n",
    "    discounted_benefits=0.0\n",
    "    discounted_costs=0.0\n",
    "    for x in range(n+constructiont):\n",
    "        t=x\n",
    "        benefit=Et[x]*(wholesale_price/1000)/np.power((1+r),t)\n",
    "        cost=(It[x]+Mt[x])/np.power((1+r),t)\n",
    "        discounted_benefits=discounted_benefits+benefit\n",
    "        discounted_costs=discounted_costs+cost\n",
    "    return discounted_benefits-discounted_costs\n",
    "\n",
    "#Assume that wholesale price simply fluctuates in time in the given range...i.e. permanently low wholesale price\n",
    "def NPV2(construction_period,n,r,cf,It0,OMkwh,OMfixed,wholesalemin,wholesalemax):\n",
    "    prod=cf*1*365*24\n",
    "    constructiont=int(round(construction_period))\n",
    "    It=np.zeros(n+constructiont)\n",
    "    #It[0]=It0\n",
    "    Mt=np.ones(n+constructiont)*(OMkwh*prod+OMfixed)\n",
    "    Et=np.ones(n+constructiont)*prod  #kwH\n",
    "    if constructiont>0.5:\n",
    "        for x in range(constructiont):\n",
    "            It[x]=It0/constructiont  #spread the investment evenly over the construction period\n",
    "            Mt[x]=0   #No costs or production during construction\n",
    "            Et[x]=0\n",
    "    else:\n",
    "        It[0]=It0\n",
    "    discounted_benefits=0.0\n",
    "    discounted_costs=0.0\n",
    "    for x in range(n+constructiont):\n",
    "        t=x\n",
    "        wholesale_price=np.random.uniform(low=wholesalemin,high=wholesalemax)\n",
    "        #wholesale_price=0.5*(wholesalemin+wholesalemax)\n",
    "        benefit=Et[x]*(wholesale_price/1000)/np.power((1+r),t)\n",
    "        cost=(It[x]+Mt[x])/np.power((1+r),t)\n",
    "        discounted_benefits=discounted_benefits+benefit\n",
    "        discounted_costs=discounted_costs+cost\n",
    "    return discounted_benefits-discounted_costs\n",
    "\n",
    "#Assume that wholesale price constant\n",
    "def NPV3(construction_period,n,r,cf,It0,OMkwh,OMfixed,wholesale):\n",
    "    prod=cf*1*365*24\n",
    "    constructiont=int(round(construction_period))\n",
    "    It=np.zeros(n+constructiont)\n",
    "    #It[0]=It0\n",
    "    Mt=np.ones(n+constructiont)*(OMkwh*prod+OMfixed)\n",
    "    Et=np.ones(n+constructiont)*prod  #kwH\n",
    "    if constructiont>0.5:\n",
    "        for x in range(constructiont):\n",
    "            It[x]=It0/constructiont  #spread the investment evenly over the construction period\n",
    "            Mt[x]=0   #No costs or production during construction\n",
    "            Et[x]=0\n",
    "    else:\n",
    "        It[0]=It0\n",
    "    discounted_benefits=0.0\n",
    "    discounted_costs=0.0\n",
    "    for x in range(n+constructiont):\n",
    "        t=x\n",
    "        wholesale_price=wholesale\n",
    "        #wholesale_price=0.5*(wholesalemin+wholesalemax)\n",
    "        benefit=Et[x]*(wholesale_price/1000)/np.power((1+r),t)\n",
    "        cost=(It[x]+Mt[x])/np.power((1+r),t)\n",
    "        discounted_benefits=discounted_benefits+benefit\n",
    "        discounted_costs=discounted_costs+cost\n",
    "    return discounted_benefits-discounted_costs\n",
    "\n",
    "\n",
    "def LCOE(construction_period,n,r,cf,It0,OMkwh,OMfixed):\n",
    "    #Returns LCOE of the energy source in currency/kWh\n",
    "    #Distribute investment evenly over the construction time...production starts after that\n",
    "    prod=cf*1*365*24\n",
    "    constructiont=int(round(construction_period))\n",
    "    It=np.zeros(n+constructiont)\n",
    "    #It[0]=It0\n",
    "    Mt=np.ones(n+constructiont)*(OMkwh*prod+OMfixed)\n",
    "    Et=np.ones(n+constructiont)*prod  #kwH\n",
    "    if constructiont>0.5:\n",
    "        for x in range(constructiont):\n",
    "            It[x]=It0/constructiont  #spread the investment evenly over the construction period\n",
    "            Mt[x]=0   #No costs or production during construction\n",
    "            Et[x]=0\n",
    "    else:\n",
    "        It[0]=It0\n",
    "        \n",
    "    discounted_costs=0.0\n",
    "    discounted_benefits=0.0\n",
    "    for x in range(n+constructiont):\n",
    "        t=x\n",
    "        benefit=Et[x]/np.power((1+r),t)\n",
    "        cost=(It[x]+Mt[x])/np.power((1+r),t)        \n",
    "        discounted_benefits=discounted_benefits+benefit\n",
    "        discounted_costs=discounted_costs+cost\n",
    "    LCOEval=discounted_costs/discounted_benefits\n",
    "    return LCOEval\n",
    "\n",
    "#def NPV3(construction_period,n,r,cf,It0,OMkwh,OMfixed,wholesale):\n",
    "#Find the discount rate at which NPV is zero\n",
    "def IRR(construction_period,n,r,cf,It0,OMkwh,OMfixed,wholesale):\n",
    "    #Root of NPV3\n",
    "    a=-0.99\n",
    "    b=0.5\n",
    "    lower=NPV3(construction_period,n,a,cf,It0,OMkwh,OMfixed,wholesale)\n",
    "    upper=NPV3(construction_period,n,b,cf,It0,OMkwh,OMfixed,wholesale)\n",
    "    if lower*upper<0.0:\n",
    "        r=opt.brentq(lambda discount: NPV3(construction_period,n,discount,cf,It0,OMkwh,OMfixed,wholesale),a,b)\n",
    "        #r=brentq(NPV3, a, b, args=(), xtol=2e-6,rtol=1.0e-6, maxiter=50,full_output=False, disp=True\n",
    "    else:\n",
    "        r=-2.0\n",
    "    return r"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}