Notebooks/07-Complex-Capacity-Expansion.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Complex Capacity Expansion\n",
    "\n",
    "_**[Power Systems Optimization](https://github.com/east-winds/power-systems-optimization)**_\n",
    "\n",
    "_by Jesse D. Jenkins and Michael R. Davison (last updated: August 31, 2021)_\n",
    "\n",
    "This notebook, the final in our series, combines elements of several prior models, including [Basic Capacity Expansion](03-Basic-Capacity-Expansion.ipynb), [Economic Dispatch](04-Economic-Dispatch.ipynb), and a [Network Flow](06-Optimal-Power-Flow.ipynb) model, to demonstrate a more complex capacity expansion planning model that includes chronologically sequential hourly economic dispatch decisions with time coupling constraints (ramp limits and energy storage) and transport flow constraints to represent power transmission limits between multiple geospatial regions.\n",
    "\n",
    "Note that detailed capacity expansion planning models can easily become large-scale constrainted optimization problems that push the limits of computational tractability. As such, these models commonly employ a range of abstraction methods along several dimensions, including temporal resolution, operational detail, and network/geospatial detail."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"img/dimensionality.png\" style=\"width: 450px; height: auto\" align=\"left\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Image source: [Jenkins & Sepulveda (2017)](https://energy.mit.edu/publication/enhanced-decision-support-changing-electricity-landscape/)*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, to keep this tutorial model relatively digestible and to keep the solution time rapid for use in an interactive notebook like this, the model presented here includes:\n",
    "\n",
    "- Economic dispatch with chronological ramp & storage constraints;\n",
    "- Linear network flow constraints between three aggregate model regions; and\n",
    "- Representative time periods (10 sample days selected via a clustering adapted from Mallapragada et al. (2018), \"[Impact of model resolution on scenario outcomes for electricity sector system expansion](https://doi.org/10.1016/j.energy.2018.08.015\n",
    ")\" *Energy* 163)\n",
    "\n",
    "Additionally all capacity investment are continuous, rather than discrete, keeping this model a linear program (LP). \n",
    "\n",
    "In the `complex_expansion_data/` path, we provide several different sets of inputs using different sample time periods, for your use and experimentation, including 10 days (used as default below), 4 weeks, 8 weeks, 16 weeks, and 52 weeks (full 8760 hours). Alter the `inputs_path` parameter below to select a different time series if desired.\n",
    "\n",
    "All input data for this model are from the open source [PowerGenome](https://github.com/gschivley/PowerGenome) data platform.\n",
    "\n",
    "This 'core' model can be extended in several directions to produce a more complicated/realistic capacity expansion planning model, incorporating discrete generator and/or transmission expansion, unit commitment constraints for thermal generators, other resources (e.g. hydropower, flexible demand, and much more), DC optimal power flow constraints, and many other elements. \n",
    "\n",
    "For examples of fully-developed complex capacity expansion models, see:\n",
    "- GenX, described in Jenkins & Sepulveda (2017), \"[Enhanced Decision Support for a Changing Electricity Landscape: The GenX Configurable Electricity Resource Capacity Expansion Model](https://energy.mit.edu/publication/enhanced-decision-support-changing-electricity-landscape/),\" MIT Energy Initiative Working Paper 2017-10\n",
    "- PyPSA, described in Brown et al. (2018), \"[PyPSA: Python for Power System Analysis](https://doi.org/10.5334/jors.188),\" *Journal of Open Research Software* and available open source at https://pypsa.org/\n",
    "- SWITCH 2.0, described in Johnston et al. (2019), \"[Switch 2.0: A modern platform for planning high-renewable power systems](https://doi.org/10.1016/j.softx.2019.100251),\" *SoftwareX* 10 and available open source at http://switch-model.org/\n",
    "\n",
    "We will dispense with a written mathematical formulation for this model, and proceed directly to the implementation in JuMP below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## LOAD PACKAGES"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "┌ Info: For saving to png with the Plotly backend PlotlyBase has to be installed.\n",
      "└ @ Plots /Users/michd/.julia/packages/Plots/4UbNP/src/backends.jl:435\n"
     ]
    }
   ],
   "source": [
    "# Uncomment this line if you need to install these packages:\n",
    "# import Pkg; Pkg.add(\"Cbc\"); Pkg.add(\"Statistics\")\n",
    "using JuMP, Clp, Cbc, GLPK\n",
    "using Plots; plotly()\n",
    "using DataFrames, CSV"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## LOAD INPUTS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read input data for a case with 10 sample days of data\n",
    "inputs_path = \"complex_expansion_data/10_days/\"\n",
    "  # Generators (and storage) data:\n",
    "generators = DataFrame(CSV.File(joinpath(inputs_path, \"Generators_data.csv\")))\n",
    "  # Many of the columns in the input data will be unused (this is input format for the GenX model)\n",
    "  # Select the ones we want for this model\n",
    "generators = select(generators, :R_ID, :Resource, :zone, :THERM, :DISP, :NDISP, :STOR, :HYDRO, :RPS, :CES,\n",
    "                    :Commit, :Existing_Cap_MW, :Existing_Cap_MWh, :Cap_size, :New_Build, :Max_Cap_MW,\n",
    "                    :Inv_cost_per_MWyr, :Fixed_OM_cost_per_MWyr, :Inv_cost_per_MWhyr, :Fixed_OM_cost_per_MWhyr,\n",
    "                    :Var_OM_cost_per_MWh, :Start_cost_per_MW, :Start_fuel_MMBTU_per_MW, :Heat_rate_MMBTU_per_MWh, :Fuel,\n",
    "                    :Min_power, :Ramp_Up_percentage, :Ramp_Dn_percentage, :Up_time, :Down_time,\n",
    "                    :Eff_up, :Eff_down);\n",
    "  # Set of all generators\n",
    "G = generators.R_ID;\n",
    "# Uncomment this line to explore the data if you wish:\n",
    "# show(generators, allrows=true, allcols=true)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "  # Read demand input data and record parameters\n",
    "demand_inputs = DataFrame(CSV.File(joinpath(inputs_path, \"Load_data.csv\")))\n",
    "# Value of lost load (cost of involuntary non-served energy)\n",
    "VOLL = demand_inputs.Voll[1]\n",
    "  # Set of price responsive demand (non-served energy) segments\n",
    "S = convert(Array{Int64}, collect(skipmissing(demand_inputs.Demand_segment))) \n",
    "#NOTE:  collect(skipmising(input)) is needed here in several spots because the demand inputs are not 'square' (different column lengths)\n",
    "\n",
    "  # Data frame for price responsive demand segments (nse)\n",
    "  # NSE_Cost = opportunity cost per MWh of demand curtailment\n",
    "  # NSE_Max = maximum % of demand that can be curtailed in each hour\n",
    "  # Note that nse segment 1 = involuntary non-served energy (load shedding) at $9000/MWh\n",
    "  # and segment 2 = one segment of voluntary price responsive demand at $600/MWh (up to 7.5% of demand)\n",
    "nse = DataFrame(Segment=S, \n",
    "                NSE_Cost = VOLL.*collect(skipmissing(demand_inputs.Cost_of_demand_curtailment_perMW)),\n",
    "                NSE_Max = collect(skipmissing(demand_inputs.Max_demand_curtailment)))\n",
    "\n",
    "  # Set of sequential hours per sub-period\n",
    "hours_per_period = convert(Int64, demand_inputs.Hours_per_period[1])\n",
    "  # Set of time sample sub-periods (e.g. sample days or weeks)\n",
    "P = convert(Array{Int64}, 1:demand_inputs.Subperiods[1])\n",
    "  # Sub period cluster weights = number of hours represented by each sample period\n",
    "W = convert(Array{Int64}, collect(skipmissing(demand_inputs.Sub_Weights)))\n",
    "  # Set of all time steps\n",
    "T = convert(Array{Int64}, demand_inputs.Time_index)\n",
    "  # Create vector of sample weights, representing how many hours in the year\n",
    "  # each hour in each sample period represents\n",
    "sample_weight = zeros(Float64, size(T,1))\n",
    "t=1\n",
    "for p in P\n",
    "    for h in 1:hours_per_period\n",
    "        sample_weight[t] = W[p]/hours_per_period\n",
    "        t=t+1\n",
    "    end\n",
    "end\n",
    "\n",
    "  # Set of zones \n",
    "Z = convert(Array{Int64}, 1:3)\n",
    "# Notes on zones: \n",
    "# Zone 1 is the Texas Panhandle, home to good wind resource but no local demand (not part of ERCOT)\n",
    "# Zone 2 is eastern half of ERCOT, home to majority of Texas population and major cities like Houston, Dallas-Forth Worth, Austin, and San Antonio\n",
    "# Zone 3 is western half of ERCOT, less populated, but great wind and solar resources\n",
    "\n",
    "  # Load/demand time series by zone (TxZ array)\n",
    "demand = select(demand_inputs, :Load_MW_z1, :Load_MW_z2, :Load_MW_z3);\n",
    "# Uncomment this line to explore the data if you wish:\n",
    "# show(demand, allrows=true, allcols=true)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "  # Read generator capacity factors by hour (used for variable renewables)\n",
    "  # There is one column here for each resource (row) in the generators DataFrame\n",
    "variability = DataFrame(CSV.File(joinpath(inputs_path, \"Generators_variability.csv\")))\n",
    "  # Drop the first column with row indexes, as these are unecessary\n",
    "variability = variability[:,2:ncol(variability)];\n",
    "# Uncomment this line to explore the data if you wish:\n",
    "# show(variability, allrows=true, allcols=true)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "  # Read fuels data\n",
    "fuels = DataFrame(CSV.File(joinpath(inputs_path, \"Fuels_data.csv\")));\n",
    "# Uncomment this line to explore the data if you wish:\n",
    "# show(fuels, allrows=true, allcols=true);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "  # Read network data\n",
    "network = DataFrame(CSV.File(joinpath(inputs_path, \"Network.csv\")));\n",
    "  #Again, there is a lot of entries in here we will not use (formatted for GenX inputs), so let's select what we want\n",
    "  # Array of network zones (z1, z2, z3)\n",
    "zones = collect(skipmissing(network.Network_zones))\n",
    "  # Network map showing lines connecting zones\n",
    "lines = select(network[1:2,:], \n",
    "    :Network_lines, :z1, :z2, :z3, \n",
    "    :Line_Max_Flow_MW, :Line_Min_Flow_MW, :Line_Loss_Percentage, \n",
    "    :Line_Max_Reinforcement_MW, :Line_Reinforcement_Cost_per_MW_yr)\n",
    "  # Add fixed O&M costs for lines = 1/20 of reinforcement cost\n",
    "lines.Line_Fixed_Cost_per_MW_yr = lines.Line_Reinforcement_Cost_per_MW_yr./20\n",
    "  # Set of all lines\n",
    "L = convert(Array{Int64}, lines.Network_lines);\n",
    "# Uncomment this line to explore the data if you wish:\n",
    "# show(lines, allrows=true, allcols=true)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate generator (and storage) total variable costs, start-up costs, \n",
    "# and associated CO2 per MWh and per start\n",
    "generators.Var_Cost = zeros(Float64, size(G,1))\n",
    "generators.CO2_Rate = zeros(Float64, size(G,1))\n",
    "generators.Start_Cost = zeros(Float64, size(G,1))\n",
    "generators.CO2_Per_Start = zeros(Float64, size(G,1))\n",
    "for g in G\n",
    "    # Variable cost ($/MWh) = variable O&M ($/MWh) + fuel cost ($/MMBtu) * heat rate (MMBtu/MWh)\n",
    "    generators.Var_Cost[g] = generators.Var_OM_cost_per_MWh[g] +\n",
    "        fuels[fuels.Fuel.==generators.Fuel[g],:Cost_per_MMBtu][1]*generators.Heat_rate_MMBTU_per_MWh[g]\n",
    "    # CO2 emissions rate (tCO2/MWh) = fuel CO2 content (tCO2/MMBtu) * heat rate (MMBtu/MWh)\n",
    "    generators.CO2_Rate[g] = fuels[fuels.Fuel.==generators.Fuel[g],:CO2_content_tons_per_MMBtu][1]*generators.Heat_rate_MMBTU_per_MWh[g]\n",
    "    # Start-up cost ($/start/MW) = start up O&M cost ($/start/MW) + fuel cost ($/MMBtu) * start up fuel use (MMBtu/start/MW) \n",
    "    generators.Start_Cost[g] = generators.Start_cost_per_MW[g] +\n",
    "        fuels[fuels.Fuel.==generators.Fuel[g],:Cost_per_MMBtu][1]*generators.Start_fuel_MMBTU_per_MW[g]\n",
    "    # Start-up CO2 emissions (tCO2/start/MW) = fuel CO2 content (tCO2/MMBtu) * start up fuel use (MMBtu/start/MW) \n",
    "    generators.CO2_Per_Start[g] = fuels[fuels.Fuel.==generators.Fuel[g],:CO2_content_tons_per_MMBtu][1]*generators.Start_fuel_MMBTU_per_MW[g]\n",
    "end\n",
    "# Note: after this, we don't need the fuels Data Frame again..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Drop hydropower and biomass plants from generators set for simplicity \n",
    "# (these are a small share of total ERCOT capacity, ~500 MW\n",
    "G = intersect(generators.R_ID[.!(generators.HYDRO.==1)],G)\n",
    "G = intersect(generators.R_ID[.!(generators.NDISP.==1)],G);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SETS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Organize all of our data in one place...\n",
    "# This code block is unecessary, but after all of the input steps above\n",
    "# writing it all here helps us see all of the sets and parameter DataFrames.\n",
    "\n",
    "# SUBSETS\n",
    "  # By naming convention, all sets are single capital letters\n",
    "G # Set of all generators\n",
    "S # Set of all non-served energy (price responsive demand) segments\n",
    "P  # Set of time sample sub-periods (e.g. sample days or weeks)\n",
    "W  # Sub period cluster weights = number of periods (days/weeks) represented by each sample period\n",
    "T  # Set of all time steps \n",
    "Z  # Set of zones (by number)\n",
    "L;  # Set of all lines"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## PARAMETERS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PARAMETER DataFrames\n",
    "  # By naming convention, all parameter data frames are lowercase\n",
    "nse        # non-served energy parameters (by s in S)\n",
    "generators # generation (and storage) parameters (by g in G)\n",
    "demand     # demand parameters (by t in T)\n",
    "zones      # network zones (by z in Z)\n",
    "lines;      # transmission lines (by l in L)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SUBSETS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "#SUBSETS\n",
    "  # By naming convention, all subsets are UPPERCASE\n",
    "\n",
    "  # Subset of G of all thermal resources subject to unit commitment constraints\n",
    "UC = intersect(generators.R_ID[generators.Commit.==1], G)\n",
    "  # Subset of G NOT subject to unit commitment constraints\n",
    "ED = intersect(generators.R_ID[.!(generators.Commit.==1)], G)\n",
    "  # Subset of G of all storage resources\n",
    "STOR = intersect(generators.R_ID[generators.STOR.>=1], G)\n",
    "  # Subset of G of all variable renewable resources\n",
    "VRE = intersect(generators.R_ID[generators.DISP.==1], G)\n",
    "  # Subset of all new build resources\n",
    "NEW = intersect(generators.R_ID[generators.New_Build.==1], G)\n",
    "  # Subset of all existing resources\n",
    "OLD = intersect(generators.R_ID[.!(generators.New_Build.==1)], G)\n",
    "  # Subset of all RPS qualifying resources\n",
    "RPS = intersect(generators.R_ID[generators.RPS.==1], G);\n",
    "\n",
    "# Notes: findall(x->x in A, B) returns the intersection of two vectors A and B"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## DEFINE MODEL"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# LP model using Clp solver\n",
    "Expansion_Model =  Model(Clp.Optimizer);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## DECISION VARIABLES"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# DECISION VARIABLES\n",
    "  # By naming convention, all decision variables start with v and then are in UPPER_SNAKE_CASE\n",
    "\n",
    "# Capacity decision variables\n",
    "@variables(Expansion_Model, begin\n",
    "        vCAP[g in G]            >= 0     # power capacity (MW)\n",
    "        vRET_CAP[g in OLD]      >= 0     # retirement of power capacity (MW)\n",
    "        vNEW_CAP[g in NEW]      >= 0     # new build power capacity (MW)\n",
    "        \n",
    "        vE_CAP[g in STOR]       >= 0     # storage energy capacity (MWh)\n",
    "        vRET_E_CAP[g in intersect(STOR, OLD)]   >= 0     # retirement of storage energy capacity (MWh)\n",
    "        vNEW_E_CAP[g in intersect(STOR, NEW)]   >= 0     # new build storage energy capacity (MWh)\n",
    "        \n",
    "        vT_CAP[l in L]          >= 0     # transmission capacity (MW)\n",
    "        vRET_T_CAP[l in L]      >= 0     # retirement of transmission capacity (MW)\n",
    "        vNEW_T_CAP[l in L]      >= 0     # new build transmission capacity (MW)\n",
    "end)\n",
    "\n",
    "# Set upper bounds on capacity for renewable resources \n",
    "# (which are limited in each resource 'cluster')\n",
    "for g in NEW[generators[NEW,:Max_Cap_MW].>0]\n",
    "    set_upper_bound(vNEW_CAP[g], generators.Max_Cap_MW[g])\n",
    "end\n",
    "\n",
    "# Set upper bounds on transmission capacity expansion\n",
    "for l in L\n",
    "    set_upper_bound(vNEW_T_CAP[l], lines.Line_Max_Reinforcement_MW[l])\n",
    "end\n",
    "\n",
    "# Operational decision variables\n",
    "@variables(Expansion_Model, begin\n",
    "        vGEN[T,G]       >= 0  # Power generation (MW)\n",
    "        vCHARGE[T,STOR] >= 0  # Power charging (MW)\n",
    "        vSOC[T,STOR]    >= 0  # Energy storage state of charge (MWh)\n",
    "        vNSE[T,S,Z]     >= 0  # Non-served energy/demand curtailment (MW)\n",
    "        vFLOW[T,L]      # Transmission line flow (MW); \n",
    "          # note line flow is positive if flowing\n",
    "          # from source node (indicated by 1 in zone column for that line) \n",
    "          # to sink node (indicated by -1 in zone column for that line); \n",
    "          # flow is negative if flowing from sink to source.\n",
    "end)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## CONSTRAINTS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# CONSTRAINTS\n",
    "  # By naming convention, all constraints start with c and then are TitleCase\n",
    "\n",
    "# (1) Supply-demand balance constraint for all time steps and zones\n",
    "@constraint(Expansion_Model, cDemandBalance[t in T, z in Z], \n",
    "        sum(vGEN[t,g] for g in intersect(generators[generators.zone.==z,:R_ID],G)) +\n",
    "        sum(vNSE[t,s,z] for s in S) - \n",
    "        sum(vCHARGE[t,g] for g in intersect(generators[generators.zone.==z,:R_ID],STOR)) -\n",
    "        demand[t,z] - \n",
    "        sum(lines[l,Symbol(string(\"z\",z))] * vFLOW[t,l] for l in L) == 0\n",
    ");\n",
    "# Notes: \n",
    "# 1. intersect(generators[generators.zone.==z,:R_ID],G) is the subset of all \n",
    "# generators/storage located at zone z in Z.\n",
    "# 2. sum(lines[l,Symbol(string(\"z\",z))].*FLOW[l,t], l in L) is the net sum of \n",
    "# all flows out of zone z (net exports) \n",
    "# 3. We use Symbol(string(\"z\",z)) to convert the numerical reference to z in Z\n",
    "# to a Symbol in set {:z1, :z2, :z3} as this is the reference to the columns\n",
    "# in the lines data for zone z indicating which whether z is a source or sink\n",
    "# for each line l in L."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (2-6) Capacitated constraints:\n",
    "@constraints(Expansion_Model, begin\n",
    "# (2) Max power constraints for all time steps and all generators/storage\n",
    "    cMaxPower[t in T, g in G], vGEN[t,g] <= variability[t,g]*vCAP[g]\n",
    "# (3) Max charge constraints for all time steps and all storage resources\n",
    "    cMaxCharge[t in T, g in STOR], vCHARGE[t,g] <= vCAP[g]\n",
    "# (4) Max state of charge constraints for all time steps and all storage resources\n",
    "    cMaxSOC[t in T, g in STOR], vSOC[t,g] <= vE_CAP[g]\n",
    "# (5) Max non-served energy constraints for all time steps and all segments and all zones\n",
    "    cMaxNSE[t in T, s in S, z in Z], vNSE[t,s,z] <= nse.NSE_Max[s]*demand[t,z]\n",
    "# (6a) Max flow constraints for all time steps and all lines\n",
    "    cMaxFlow[t in T, l in L], vFLOW[t,l] <= vT_CAP[l]\n",
    "# (6b) Min flow constraints for all time steps and all lines\n",
    "    cMinFlow[t in T, l in L], vFLOW[t,l] >= -vT_CAP[l]\n",
    "end)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (7-9) Total capacity constraints:\n",
    "@constraints(Expansion_Model, begin\n",
    "# (7a) Total capacity for existing units\n",
    "    cCapOld[g in OLD], vCAP[g] == generators.Existing_Cap_MW[g] - vRET_CAP[g]\n",
    "# (7b) Total capacity for new units\n",
    "    cCapNew[g in NEW], vCAP[g] == vNEW_CAP[g]\n",
    "        \n",
    "# (8a) Total energy storage capacity for existing units\n",
    "    cCapEnergyOld[g in intersect(STOR, OLD)], \n",
    "        vE_CAP[g] == generators.Existing_Cap_MWh[g] - vRET_E_CAP[g]\n",
    "# (8b) Total energy storage capacity for existing units\n",
    "    cCapEnergyNew[g in intersect(STOR, NEW)], \n",
    "        vE_CAP[g] == vNEW_E_CAP[g]\n",
    "        \n",
    "# (9) Total transmission capacity\n",
    "    cTransCap[l in L], vT_CAP[l] == lines.Line_Max_Flow_MW[l] - vRET_T_CAP[l] + vNEW_T_CAP[l]\n",
    "end)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Because we are using time domain reduction via sample periods (days or weeks),\n",
    "# we must be careful with time coupling constraints at the start and end of each\n",
    "# sample period. \n",
    "\n",
    " # First we record a subset of time steps that begin a sub period \n",
    " # (these will be subject to 'wrapping' constraints that link the start/end of each period)\n",
    "STARTS = 1:hours_per_period:maximum(T)        \n",
    " # Then we record all time periods that do not begin a sub period \n",
    "# (these will be subject to normal time couping constraints, looking back one period)\n",
    "INTERIORS = setdiff(T,STARTS)\n",
    "\n",
    "# (10-12) Time coupling constraints\n",
    "@constraints(Expansion_Model, begin\n",
    "    # (10a) Ramp up constraints, normal\n",
    "    cRampUp[t in INTERIORS, g in G], \n",
    "        vGEN[t,g] - vGEN[t-1,g] <= generators.Ramp_Up_percentage[g]*vCAP[g]\n",
    "    # (10b) Ramp up constraints, sub-period wrapping\n",
    "    cRampUpWrap[t in STARTS, g in G], \n",
    "        vGEN[t,g] - vGEN[t+hours_per_period-1,g] <= generators.Ramp_Up_percentage[g]*vCAP[g]    \n",
    "    \n",
    "    # (11a) Ramp down, normal\n",
    "    cRampDown[t in INTERIORS, g in G], \n",
    "        vGEN[t-1,g] - vGEN[t,g] <= generators.Ramp_Dn_percentage[g]*vCAP[g] \n",
    "    # (11b) Ramp down, sub-period wrapping\n",
    "    cRampDownWrap[t in STARTS, g in G], \n",
    "        vGEN[t+hours_per_period-1,g] - vGEN[t,g] <= generators.Ramp_Dn_percentage[g]*vCAP[g]     \n",
    "   \n",
    "    # (12a) Storage state of charge, normal\n",
    "    cSOC[t in INTERIORS, g in STOR], \n",
    "        vSOC[t,g] == vSOC[t-1,g] + generators.Eff_up[g]*vCHARGE[t,g] - vGEN[t,g]/generators.Eff_down[g]\n",
    "    # (12a) Storage state of charge, wrapping\n",
    "    cSOCWrap[t in STARTS, g in STOR], \n",
    "        vSOC[t,g] == vSOC[t+hours_per_period-1,g] + generators.Eff_up[g]*vCHARGE[t,g] - vGEN[t,g]/generators.Eff_down[g]\n",
    "end)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## OBJECTIVE FUNCTION"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The objective function is to minimize the sum of fixed costs associated with\n",
    "# capacity decisions and variable costs associated with operational decisions\n",
    "\n",
    "# Create expressions for each sub-component of the total cost (for later retrieval)\n",
    "@expression(Expansion_Model, eFixedCostsGeneration,\n",
    "     # Fixed costs for total capacity \n",
    "    sum(generators.Fixed_OM_cost_per_MWyr[g]*vCAP[g] for g in G) +\n",
    "     # Investment cost for new capacity\n",
    "    sum(generators.Inv_cost_per_MWyr[g]*vNEW_CAP[g] for g in NEW)\n",
    ")\n",
    "@expression(Expansion_Model, eFixedCostsStorage,\n",
    "     # Fixed costs for total storage energy capacity \n",
    "    sum(generators.Fixed_OM_cost_per_MWhyr[g]*vE_CAP[g] for g in STOR) + \n",
    "     # Investment costs for new storage energy capacity\n",
    "    sum(generators.Inv_cost_per_MWhyr[g]*vNEW_E_CAP[g] for g in intersect(STOR, NEW))\n",
    ")\n",
    "@expression(Expansion_Model, eFixedCostsTransmission,\n",
    "     # Investment and fixed O&M costs for transmission lines\n",
    "    sum(lines.Line_Fixed_Cost_per_MW_yr[l]*vT_CAP[l] +\n",
    "        lines.Line_Reinforcement_Cost_per_MW_yr[l]*vNEW_T_CAP[l] for l in L)\n",
    ")\n",
    "@expression(Expansion_Model, eVariableCosts,\n",
    "     # Variable costs for generation, weighted by hourly sample weight\n",
    "    sum(sample_weight[t]*generators.Var_Cost[g]*vGEN[t,g] for t in T, g in G)\n",
    ")\n",
    "@expression(Expansion_Model, eNSECosts,\n",
    "     # Non-served energy costs\n",
    "    sum(sample_weight[t]*nse.NSE_Cost[s]*vNSE[t,s,z] for t in T, s in S, z in Z)\n",
    ")\n",
    "  \n",
    "@objective(Expansion_Model, Min,\n",
    "    eFixedCostsGeneration + eFixedCostsStorage + eFixedCostsTransmission +\n",
    "    eVariableCosts + eNSECosts\n",
    ");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## RUN MODEL"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " 10.876692 seconds (20.38 M allocations: 1.212 GiB, 4.36% gc time, 0.87% compilation time)\n",
      "Coin0506I Presolve 37438 (-5339) rows, 14139 (-1817) columns and 110862 (-12534) elements\n",
      "Clp0006I 0  Obj 0 Primal inf 23505530 (1200)\n",
      "Clp0006I 412  Obj 4973414.7 Primal inf 24018742 (2126)\n",
      "Clp0006I 824  Obj 3.20547e+08 Primal inf 26141575 (2938)\n",
      "Clp0006I 1236  Obj 1.0280627e+09 Primal inf 26650748 (3649)\n",
      "Clp0006I 1648  Obj 2.1599333e+09 Primal inf 25467993 (4241)\n",
      "Clp0006I 2060  Obj 3.4673683e+09 Primal inf 25216069 (4699)\n",
      "Clp0006I 2472  Obj 3.8453983e+09 Primal inf 25791420 (5183)\n",
      "Clp0006I 2884  Obj 5.0810923e+09 Primal inf 28256500 (5549)\n",
      "Clp0006I 3296  Obj 6.4672279e+09 Primal inf 27641049 (5808)\n",
      "Clp0006I 3708  Obj 6.7660465e+09 Primal inf 19201947 (6165)\n",
      "Clp0006I 4120  Obj 7.3238473e+09 Primal inf 26933070 (6439)\n",
      "Clp0006I 4532  Obj 7.7672718e+09 Primal inf 43669479 (7517)\n",
      "Clp0006I 4944  Obj 8.2188563e+09 Primal inf 22268741 (6771)\n",
      "Clp0006I 5356  Obj 8.4668867e+09 Primal inf 17126438 (6882)\n",
      "Clp0006I 5768  Obj 8.8812114e+09 Primal inf 38830060 (7821)\n",
      "Clp0006I 6180  Obj 9.2049772e+09 Primal inf 37928485 (8392)\n",
      "Clp0006I 6592  Obj 9.5690571e+09 Primal inf 8778011.6 (6316)\n",
      "Clp0006I 7004  Obj 9.9196412e+09 Primal inf 10735438 (6539)\n",
      "Clp0006I 7416  Obj 1.0267453e+10 Primal inf 8157148.2 (6339)\n",
      "Clp0006I 7828  Obj 1.0394518e+10 Primal inf 25850554 (7674)\n",
      "Clp0006I 8240  Obj 1.0609354e+10 Primal inf 10619101 (7419)\n",
      "Clp0006I 8652  Obj 1.0746493e+10 Primal inf 14210861 (6603)\n",
      "Clp0006I 9064  Obj 1.0971659e+10 Primal inf 1.1940173e+08 (8475)\n",
      "Clp0006I 9476  Obj 1.1235954e+10 Primal inf 11379723 (6478)\n",
      "Clp0006I 9888  Obj 1.145209e+10 Primal inf 12039141 (6973)\n",
      "Clp0006I 10300  Obj 1.1750363e+10 Primal inf 8.8581977e+08 (7632)\n",
      "Clp0006I 10712  Obj 1.1993498e+10 Primal inf 16939044 (6911)\n",
      "Clp0006I 11124  Obj 1.2207603e+10 Primal inf 5989987.1 (6179)\n",
      "Clp0006I 11536  Obj 1.2479164e+10 Primal inf 1.6482241e+08 (7290)\n",
      "Clp0006I 11948  Obj 1.2627664e+10 Primal inf 2.5085544e+08 (8804)\n",
      "Clp0006I 12360  Obj 1.2850092e+10 Primal inf 33697694 (7932)\n",
      "Clp0006I 12772  Obj 1.3000754e+10 Primal inf 1.7648947e+08 (8197)\n",
      "Clp0006I 13184  Obj 1.3210868e+10 Primal inf 6086538 (5900)\n",
      "Clp0006I 13596  Obj 1.3488819e+10 Primal inf 3487703.9 (5654)\n",
      "Clp0006I 14008  Obj 1.3692361e+10 Primal inf 8120269.9 (6291)\n",
      "Clp0006I 14420  Obj 1.3820945e+10 Primal inf 24564816 (4923)\n",
      "Clp0006I 14832  Obj 1.3949836e+10 Primal inf 8130172.2 (5003)\n",
      "Clp0006I 15244  Obj 1.4105273e+10 Primal inf 8548259 (4512)\n",
      "Clp0006I 15656  Obj 1.4159479e+10 Primal inf 983694 (2788)\n",
      "Clp0006I 16068  Obj 1.4193138e+10 Primal inf 651020.38 (1705)\n",
      "Clp0006I 16480  Obj 1.4196599e+10 Primal inf 34.861779 (425)\n",
      "Clp0006I 16892  Obj 1.4196631e+10 Primal inf 11.160862 (151)\n",
      "Clp0006I 17043  Obj 1.4196651e+10\n",
      "Clp0000I Optimal - objective value 1.4196651e+10\n",
      "Coin0511I After Postsolve, objective 1.4196651e+10, infeasibilities - dual 0 (0), primal 0 (0)\n",
      "Clp0032I Optimal objective 1.419665113e+10 - 17043 iterations time 1.972, Presolve 0.05\n"
     ]
    }
   ],
   "source": [
    "@time optimize!(Expansion_Model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## EXTRACT RESULTS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Record generation capacity and energy results\n",
    "generation = zeros(size(G,1))\n",
    "for i in 1:size(G,1)\n",
    "    # Note that total annual generation is sumproduct of sample period weights and hourly sample period generation \n",
    "    generation[i] = sum(sample_weight.*value.(vGEN)[:,G[i]].data) \n",
    "end\n",
    "\n",
    "# Total annual demand is sumproduct of sample period weights and hourly sample period demands\n",
    "total_demand = sum(sum.(eachcol(sample_weight.*demand)))\n",
    "# Maximum aggregate demand is the maximum of the sum of total concurrent demand in each hour\n",
    "peak_demand = maximum(sum(eachcol(demand)))\n",
    "MWh_share = generation./total_demand.*100\n",
    "cap_share = value.(vCAP).data./peak_demand.*100\n",
    "generator_results = DataFrame(\n",
    "    ID = G, \n",
    "    Resource = generators.Resource[G],\n",
    "    Zone = generators.zone[G],\n",
    "    Total_MW = value.(vCAP).data,\n",
    "    Start_MW = generators.Existing_Cap_MW[G],\n",
    "    Change_in_MW = value.(vCAP).data.-generators.Existing_Cap_MW[G],\n",
    "    Percent_MW = cap_share,\n",
    "    GWh = generation/1000,\n",
    "    Percent_GWh = MWh_share\n",
    ")\n",
    "\n",
    "# Record energy storage energy capacity results (MWh)\n",
    "storage_results = DataFrame(\n",
    "    ID = STOR, \n",
    "    Zone = generators.zone[STOR],\n",
    "    Resource = generators.Resource[STOR],\n",
    "    Total_Storage_MWh = value.(vE_CAP).data,\n",
    "    Start_Storage_MWh = generators.Existing_Cap_MWh[STOR],\n",
    "    Change_in_Storage_MWh = value.(vE_CAP).data.-generators.Existing_Cap_MWh[STOR],\n",
    ")\n",
    "\n",
    "\n",
    "# Record transmission capacity results\n",
    "transmission_results = DataFrame(\n",
    "    Line = L, \n",
    "    Total_Transfer_Capacity = value.(vT_CAP).data,\n",
    "    Start_Transfer_Capacity = lines.Line_Max_Flow_MW,\n",
    "    Change_in_Transfer_Capacity = value.(vT_CAP).data.-lines.Line_Max_Flow_MW,\n",
    ")\n",
    "\n",
    "\n",
    "## Record non-served energy results by segment and zone\n",
    "num_segments = maximum(S)\n",
    "num_zones = maximum(Z)\n",
    "nse_results = DataFrame(\n",
    "    Segment = zeros(num_segments*num_zones),\n",
    "    Zone = zeros(num_segments*num_zones),\n",
    "    NSE_Price = zeros(num_segments*num_zones),\n",
    "    Max_NSE_MW = zeros(num_segments*num_zones),\n",
    "    Total_NSE_MWh = zeros(num_segments*num_zones),\n",
    "    NSE_Percent_of_Demand = zeros(num_segments*num_zones)\n",
    ")\n",
    "i=1\n",
    "for s in S\n",
    "    for z in Z\n",
    "        nse_results.Segment[i]=s\n",
    "        nse_results.Zone[i]=z\n",
    "        nse_results.NSE_Price[i]=nse.NSE_Cost[s]\n",
    "        nse_results.Max_NSE_MW[i]=maximum(value.(vNSE)[:,s,z].data)\n",
    "        nse_results.Total_NSE_MWh[i]=sum(sample_weight.*value.(vNSE)[:,s,z].data)\n",
    "        nse_results.NSE_Percent_of_Demand[i]=sum(sample_weight.*value.(vNSE)[:,s,z].data)/total_demand*100\n",
    "        i=i+1\n",
    "    end\n",
    "end\n",
    "\n",
    "# Record costs by component (in million dollars)\n",
    " # Note: because each expression evaluates to a single value, \n",
    " # value.(JuMPObject) returns a numerical value, not a DenseAxisArray;\n",
    " # We thus do not need to use the .data extension here to extract numeric values\n",
    "cost_results = DataFrame(\n",
    "    Total_Costs = objective_value(Expansion_Model)/10^6,\n",
    "    Fixed_Costs_Generation = value.(eFixedCostsGeneration)/10^6,\n",
    "    Fixed_Costs_Storage = value.(eFixedCostsStorage)/10^6,\n",
    "    Fixed_Costs_Transmission = value.(eFixedCostsTransmission)/10^6,\n",
    "    Variable_Costs = value.(eVariableCosts)/10^6,\n",
    "    NSE_Costs = value.(eNSECosts)/10^6\n",
    ");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Write results to file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Output path [set path to desired output directory here]\n",
    "outpath = \"results\" #/YOUR/PATH/HERE\n",
    "\n",
    "# If output directory does not exist, create it\n",
    "if !(isdir(outpath))\n",
    "    mkdir(outpath)\n",
    "end\n",
    "\n",
    "CSV.write(joinpath(outpath, \"generator_results.csv\"), generator_results)\n",
    "CSV.write(joinpath(outpath, \"storage_results.csv\"), storage_results)\n",
    "CSV.write(joinpath(outpath, \"transmission_results.csv\"), transmission_results)\n",
    "CSV.write(joinpath(outpath, \"nse_results.csv\"), nse_results)\n",
    "CSV.write(joinpath(outpath, \"cost_results.csv\"), cost_results);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Julia 1.6.2",
   "language": "julia",
   "name": "julia-1.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}