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: November 6, 2023)_\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 (no [unit commitment](05-Unit-Commitment.ipynb));\n",
    "- Linear network flow constraints between three aggregate model regions (no [DC OPF](06-Optimal-Power-Flow.ipynb)); and\n",
    "- Representative time periods (10 sample days selected via a clustering method 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, and [available open source on Github here](https://github.com/GenXProject/GenX/) \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",
    "- GridPath available at https://www.gridpath.io/, with [documentation](https://gridpath.readthedocs.io/en/latest/) and [available open source at Github here](https://github.com/blue-marble/gridpath)\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": [],
   "source": [
    "# Ensure these packages are installed; if not, use the Pkg package and Pkg.add() function to install\n",
    "using JuMP\n",
    "using HiGHS\n",
    "using Plots\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": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div><div style = \"float: left;\"><span>2×10 DataFrame</span></div><div style = \"clear: both;\"></div></div><div class = \"data-frame\" style = \"overflow-x: scroll;\"><table class = \"data-frame\" style = \"margin-bottom: 6px;\"><thead><tr class = \"header\"><th class = \"rowNumber\" style = \"font-weight: bold; text-align: right;\">Row</th><th style = \"text-align: left;\">Network_lines</th><th style = \"text-align: left;\">z1</th><th style = \"text-align: left;\">z2</th><th style = \"text-align: left;\">z3</th><th style = \"text-align: left;\">Line_Max_Flow_MW</th><th style = \"text-align: left;\">Line_Min_Flow_MW</th><th style = \"text-align: left;\">Line_Loss_Percentage</th><th style = \"text-align: left;\">Line_Max_Reinforcement_MW</th><th style = \"text-align: left;\">Line_Reinforcement_Cost_per_MW_yr</th><th style = \"text-align: left;\">Line_Fixed_Cost_per_MW_yr</th></tr><tr class = \"subheader headerLastRow\"><th class = \"rowNumber\" style = \"font-weight: bold; text-align: right;\"></th><th title = \"Union{Missing, Float64}\" style = \"text-align: left;\">Float64?</th><th title = \"Union{Missing, Float64}\" style = \"text-align: left;\">Float64?</th><th title = \"Union{Missing, Float64}\" style = \"text-align: left;\">Float64?</th><th title = \"Union{Missing, Float64}\" style = \"text-align: left;\">Float64?</th><th title = \"Union{Missing, Float64}\" style = \"text-align: left;\">Float64?</th><th title = \"Union{Missing, Float64}\" style = \"text-align: left;\">Float64?</th><th title = \"Union{Missing, Float64}\" style = \"text-align: left;\">Float64?</th><th title = \"Union{Missing, Float64}\" style = \"text-align: left;\">Float64?</th><th title = \"Union{Missing, Float64}\" style = \"text-align: left;\">Float64?</th><th title = \"Float64\" style = \"text-align: left;\">Float64</th></tr></thead><tbody><tr><td class = \"rowNumber\" style = \"font-weight: bold; text-align: right;\">1</td><td style = \"text-align: right;\">1.0</td><td style = \"text-align: right;\">1.0</td><td style = \"text-align: right;\">0.0</td><td style = \"text-align: right;\">-1.0</td><td style = \"text-align: right;\">3702.0</td><td style = \"text-align: right;\">-3702.0</td><td style = \"text-align: right;\">0.0196</td><td style = \"text-align: right;\">3702.0</td><td style = \"text-align: right;\">21032.0</td><td style = \"text-align: right;\">1051.6</td></tr><tr><td class = \"rowNumber\" style = \"font-weight: bold; text-align: right;\">2</td><td style = \"text-align: right;\">2.0</td><td style = \"text-align: right;\">0.0</td><td style = \"text-align: right;\">1.0</td><td style = \"text-align: right;\">-1.0</td><td style = \"text-align: right;\">5529.0</td><td style = \"text-align: right;\">-10555.0</td><td style = \"text-align: right;\">0.0256</td><td style = \"text-align: right;\">5529.0</td><td style = \"text-align: right;\">27595.0</td><td style = \"text-align: right;\">1379.75</td></tr></tbody></table></div>"
      ],
      "text/latex": [
       "\\begin{tabular}{r|ccccccc}\n",
       "\t& Network\\_lines & z1 & z2 & z3 & Line\\_Max\\_Flow\\_MW & Line\\_Min\\_Flow\\_MW & \\\\\n",
       "\t\\hline\n",
       "\t& Float64? & Float64? & Float64? & Float64? & Float64? & Float64? & \\\\\n",
       "\t\\hline\n",
       "\t1 & 1.0 & 1.0 & 0.0 & -1.0 & 3702.0 & -3702.0 & $\\dots$ \\\\\n",
       "\t2 & 2.0 & 0.0 & 1.0 & -1.0 & 5529.0 & -10555.0 & $\\dots$ \\\\\n",
       "\\end{tabular}\n"
      ],
      "text/plain": [
       "\u001b[1m2×10 DataFrame\u001b[0m\n",
       "\u001b[1m Row \u001b[0m│\u001b[1m Network_lines \u001b[0m\u001b[1m z1       \u001b[0m\u001b[1m z2       \u001b[0m\u001b[1m z3       \u001b[0m\u001b[1m Line_Max_Flow_MW \u001b[0m\u001b[1m Line_Min\u001b[0m ⋯\n",
       "     │\u001b[90m Float64?      \u001b[0m\u001b[90m Float64? \u001b[0m\u001b[90m Float64? \u001b[0m\u001b[90m Float64? \u001b[0m\u001b[90m Float64?         \u001b[0m\u001b[90m Float64?\u001b[0m ⋯\n",
       "─────┼──────────────────────────────────────────────────────────────────────────\n",
       "   1 │           1.0       1.0       0.0      -1.0            3702.0           ⋯\n",
       "   2 │           2.0       0.0       1.0      -1.0            5529.0\n",
       "\u001b[36m                                                               5 columns omitted\u001b[0m"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lines"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Vector{String3}:\n",
       " \"z1\"\n",
       " \"z2\"\n",
       " \"z3\""
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "zones"
   ]
  },
  {
   "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",
    "# SETS\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 transmission lines (or paths)"
   ]
  },
  {
   "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) also returns the intersection of two vectors A and B and could be used here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## DEFINE MODEL"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# LP model using HiGHS solver\n",
    "Expansion_Model =  Model(HiGHS.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",
    "#NOTE that since we are modeling representative time periods, \n",
    "# all operations-related costs have to be weighted by the hourly weight \n",
    "# (indicating how many hours in the full year the modeled hour represents)\n",
    "# in order to ensure operations-related costs estimate annual costs and\n",
    "# are appropriately equivalent to annualized investment/fixed costs.\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, weighted by hourly sample weight to ensure non-served energy costs estimate annual 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": [
      "Running HiGHS 1.5.3 [date: 1970-01-01, git hash: 45a127b78]\n",
      "Copyright (c) 2023 HiGHS under MIT licence terms\n",
      "Presolving model\n",
      "37442 rows, 14141 cols, 110872 nonzeros\n",
      "37438 rows, 14141 cols, 110864 nonzeros\n",
      "Presolve : Reductions: rows 37438(-5339); columns 14141(-1815); elements 110864(-12532)\n",
      "Solving the presolved LP\n",
      "Using EKK dual simplex solver - serial\n",
      "  Iteration        Objective     Infeasibilities num(sum)\n",
      "          0     0.0000000000e+00 Pr: 482(3.17783e+06); Du: 0(1.92253e-10) 0s\n",
      "      13904     1.4196651126e+10 Pr: 0(0); Du: 0(2.29937e-11) 1s\n",
      "Solving the original LP from the solution after postsolve\n",
      "Model   status      : Optimal\n",
      "Simplex   iterations: 13904\n",
      "Objective value     :  1.4196651126e+10\n",
      "HiGHS run time      :          1.23\n",
      "  1.612698 seconds (592.20 k allocations: 47.414 MiB, 1.37% gc time, 14.35% compilation time: 78% of which was recompilation)\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": 25,
   "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.9.3",
   "language": "julia",
   "name": "julia-1.9"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.9.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}