# Economic Dispatch

_**[Power Systems Optimization](https://github.com/east-winds/power-systems-optimization)**_

_by Michael R. Davidson and Jesse D. Jenkins (last updated: August 31, 2021)_

This notebook will introduce a key operational model&mdash;economic dispatch (ED)&mdash;which minimizes the short-run production costs of meeting electricity demand from a given set of generators subject to various technical constraints.

We build up the model in several stages. We begin with a single-time period and simple generator constraints just related to allowable output ranges. This allows us to compare the model with the intuition provided by the merit order supply curve (also known as the "dispatch stack"). 

We then add multiple time periods by considering the economic dispatch over an entire day. Due to the simple generator constraints, this simple model does not introduce any coupling across time.

More realism is added by considering complex engineering constraints that introduce time coupling&mdash;namely, ramp rates, which are an important limitation of many generators.

Finally, we examine the famous "duck curve" by adding more solar onto the system and explore the role that ramp rates play in the solution.

## Introduction to economic dispatch

Every day, system operators need to decide how to meet demand from a large variety of generators with different costs and engineering requirements. Prior to large-scale computational capabilities, this was accomplished by "priority lists" or similar heuristics to determine the ordering of plants. However, with modern optimization algorithms, operators can find the global minimum of production costs&mdash;and even incorporate coupling of engineering processes across time periods.

**Economic dispatch** (ED) is the problem of minimizing short-run costs of production in order to meet a given demand, considering relevant engineering constraints of the generators. It *does not* fully capture the physics of electricity flows, which will further constrain the feasible space and which we will come back to in later notebooks on optimal power flow, nor does it account for decisions and constraints related to turning on or "committing" large thermal generators, which we'll discuss in our next notebook on unit commitment. Typically, *it does* incorporate some type of network representation, which we have ignored here by considering only a single bus. More on the networks in the notebook on optimal power flow.

## Single-time period, simple generator constraints

We will first examine the case where we are optimizing dispatch for a single snapshot in time, with only very simple constraints on the generators.
    
$$
\begin{align}
\min \ & \sum_{g \in G} VarCost_g \times GEN_g & \\
\text{s.t.} & \\
 & \sum_{g} GEN_g = Demand & \\
 & GEN_g \leq Pmax_g & \forall \quad g \in G \\
 & GEN_g \geq Pmin_g & \forall \quad g \in G \\
\end{align}
$$

The **decision variable** in the above problem is:

- $GEN_{g}$, the generation (in MW) produced by each generator, $g$

The **parameters** are:

- $Pmin_g$, the minimum operating bounds for the generator (based on engineering or natural resource constraints)
- $Pmax_g$, the maximum operating bounds for the generator (based on engineering or natural resource constraints)
- $Demand$, the demand (in MW)

(For this simple problem, we will let $Pmin_g=0$ and revisit this later in the unit commitment notebook.)

In addition, we have:

$$
VarCost_g = VarOM_g + HeatRate_g \times FuelCost_g
$$

Note that, in contrast with the basic capacity expansion problem, we are not concerned with fixed costs. Why not?

The simplest response is that these costs have already been incurred and regardless of how much a generator produces, its fixed costs will not change. These costs are thus "[sunk costs](https://en.wikipedia.org/wiki/Sunk_cost)" and are constant in the objective function. Our optimal decision variables would not change by adjusting this constant up or down.

Therefore, we can safely ignore fixed costs for the purposes of finding optimal dispatch. (We will still have to consider them to calculate producer profits, however.)

[A famous electricity regulator](https://en.wikipedia.org/wiki/Ignacio_J._P%C3%A9rez_Arriaga), noting that it is frequently misunderstood and ignored even by experienced market participants, has called the application of the "[sunk cost fallacy](https://en.wikipedia.org/wiki/Sunk_cost#Fallacy_effect)" to electricity systems "[Grandma's Inheritance Theorem](http://link.springer.com/content/pdf/bbm%3A978-1-4471-5034-3%2F1.pdf)." In this hypothetical scenario, if you were to inherit a diamond ring or perhaps a financial contract-for-differences that pays out relative to cleared market prices (you have a very sophisticated grandma in this scenario!), you should not be tempted to change your production strategy.

Now, let's implement ED.

### 1. Load packages

We will try out another open source solver GLPK in this tutorial.


```julia
# New packages introduced in this tutorial
import Pkg; Pkg.add("PrettyTables"); Pkg.add("VegaLite"); Pkg.add("GLPK"); Pkg.add("FileIO")
```

    [32m[1m    Updating[22m[39m registry at `~/.julia/registries/General`
    [32m[1m    Updating[22m[39m git-repo `https://github.com/JuliaRegistries/General.git`
    [32m[1m   Resolving[22m[39m package versions...
    [32m[1m  No Changes[22m[39m to `~/.julia/environments/v1.6/Project.toml`
    [32m[1m  No Changes[22m[39m to `~/.julia/environments/v1.6/Manifest.toml`
    [32m[1m   Resolving[22m[39m package versions...
    [32m[1m  No Changes[22m[39m to `~/.julia/environments/v1.6/Project.toml`
    [32m[1m  No Changes[22m[39m to `~/.julia/environments/v1.6/Manifest.toml`
    [32m[1m   Resolving[22m[39m package versions...
    [32m[1m  No Changes[22m[39m to `~/.julia/environments/v1.6/Project.toml`
    [32m[1m  No Changes[22m[39m to `~/.julia/environments/v1.6/Manifest.toml`
    [32m[1m   Resolving[22m[39m package versions...
    [32m[1m  No Changes[22m[39m to `~/.julia/environments/v1.6/Project.toml`
    [32m[1m  No Changes[22m[39m to `~/.julia/environments/v1.6/Manifest.toml`



```julia
using JuMP, GLPK
using Plots; plotly();
using VegaLite  # to make some nice plots
using DataFrames, CSV, PrettyTables
using FileIO
ENV["COLUMNS"]=120; # Set so all columns of DataFrames and Matrices are displayed
```

    ┌ Info: For saving to png with the Plotly backend PlotlyBase has to be installed.
    └ @ Plots /Users/michd/.julia/packages/Plots/4UbNP/src/backends.jl:435


### 2. Load and format data

We will use data based on San Diego Gas and Electric (SDG&E, via the [PowerGenome](https://github.com/gschivley/PowerGenome) data platform) plus a few neighboring generators, consisting of:

- 25 generators (including some clustering of smaller generators, and excluding behind-the-meter solar)
- estimated hourly demand for 2020 (net load at the transmission substation level after subtracting 600MW of behind-the-meter solar from original demand)
- variable generation capacity factors
- estimated natural gas fuel costs


```julia
datadir = joinpath("ed_data") 
# Note: joinpath is a good way to create path reference that is agnostic
# to what file system you are using (e.g. whether directories are denoted 
# with a forward or backwards slash).
gen_info = CSV.read(joinpath(datadir,"Generators_data.csv"), DataFrame);
fuels = CSV.read(joinpath(datadir,"Fuels_data.csv"), DataFrame);
loads = CSV.read(joinpath(datadir,"Demand.csv"), DataFrame);
gen_variable = CSV.read(joinpath(datadir,"Generators_variability.csv"), DataFrame);

# Rename all columns to lowercase (by convention)
for f in [gen_info, fuels, loads, gen_variable]
    rename!(f,lowercase.(names(f)))
end
```

**Construct generator dataframe**


```julia
# Keep only the columns relevant to our ED model 
# (We'll come back to other columns in our unit commitment notebooks)
select!(gen_info, 1:26, :stor) 
gen_df = outerjoin(gen_info,  fuels, on = :fuel) # load in fuel costs and add to data frame
rename!(gen_df, :cost_per_mmbtu => :fuel_cost)   # rename column for fuel cost
gen_df[ismissing.(gen_df[:,:fuel_cost]), :fuel_cost] .= 0

# create "is_variable" column to indicate if this is a variable generation source (e.g. wind, solar)
# Note: Julia's strict vectorization syntax requires us to create a new variable using the '!' indexing operator
# and the '.=' broadcasting assignment. See: https://julia.guide/broadcasting
gen_df[!, :is_variable] .= false
gen_df[in(["onshore_wind_turbine","small_hydroelectric","solar_photovoltaic"]).(gen_df.resource),
    :is_variable] .= true;

# create full name of generator (including geographic location and cluster number)
#  for use with variable generation dataframe
gen_df.gen_full = lowercase.(gen_df.region .* "_" .* gen_df.resource .* "_" .* string.(gen_df.cluster) .* ".0");

# remove generators with no capacity (e.g. new build options that we'd use if this was capacity expansion problem)
gen_df = gen_df[gen_df.existing_cap_mw .> 0,:];

```

### Long and wide data format example


```julia
df = gen_df[1:5,1:5]
```




<div class="data-frame"><p>5 rows × 5 columns</p><table class="data-frame"><thead><tr><th></th><th>r_id</th><th>resource</th><th>region</th><th>existing_cap_mw</th><th>num_units</th></tr><tr><th></th><th title="Union{Missing, Int64}">Int64?</th><th title="Union{Missing, String}">String?</th><th title="Union{Missing, String}">String?</th><th title="Union{Missing, Float64}">Float64?</th><th title="Union{Missing, Int64}">Int64?</th></tr></thead><tbody><tr><th>1</th><td>1</td><td>biomass</td><td>WEC_SDGE</td><td>21.4</td><td>20</td></tr><tr><th>2</th><td>2</td><td>hydroelectric_pumped_storage</td><td>WEC_SDGE</td><td>42.0</td><td>2</td></tr><tr><th>3</th><td>3</td><td>natural_gas_fired_combined_cycle</td><td>WEC_SDGE</td><td>607.1</td><td>1</td></tr><tr><th>4</th><td>4</td><td>natural_gas_fired_combined_cycle</td><td>WEC_SDGE</td><td>30.0</td><td>1</td></tr><tr><th>5</th><td>5</td><td>natural_gas_fired_combined_cycle</td><td>WEC_SDGE</td><td>49.4</td><td>1</td></tr></tbody></table></div>




```julia
df = copy(gen_df)
select!(df, 1:5)
df = df[1:5,:]
```




<div class="data-frame"><p>5 rows × 5 columns</p><table class="data-frame"><thead><tr><th></th><th>r_id</th><th>resource</th><th>region</th><th>existing_cap_mw</th><th>num_units</th></tr><tr><th></th><th title="Union{Missing, Int64}">Int64?</th><th title="Union{Missing, String}">String?</th><th title="Union{Missing, String}">String?</th><th title="Union{Missing, Float64}">Float64?</th><th title="Union{Missing, Int64}">Int64?</th></tr></thead><tbody><tr><th>1</th><td>1</td><td>biomass</td><td>WEC_SDGE</td><td>21.4</td><td>20</td></tr><tr><th>2</th><td>2</td><td>hydroelectric_pumped_storage</td><td>WEC_SDGE</td><td>42.0</td><td>2</td></tr><tr><th>3</th><td>3</td><td>natural_gas_fired_combined_cycle</td><td>WEC_SDGE</td><td>607.1</td><td>1</td></tr><tr><th>4</th><td>4</td><td>natural_gas_fired_combined_cycle</td><td>WEC_SDGE</td><td>30.0</td><td>1</td></tr><tr><th>5</th><td>5</td><td>natural_gas_fired_combined_cycle</td><td>WEC_SDGE</td><td>49.4</td><td>1</td></tr></tbody></table></div>




```julia
stack(df, 
    [:existing_cap_mw, :num_units], 
    variable_name=:var,
    value_name=:val)
```




<div class="data-frame"><p>10 rows × 5 columns</p><table class="data-frame"><thead><tr><th></th><th>r_id</th><th>resource</th><th>region</th><th>var</th><th>val</th></tr><tr><th></th><th title="Union{Missing, Int64}">Int64?</th><th title="Union{Missing, String}">String?</th><th title="Union{Missing, String}">String?</th><th title="String">String</th><th title="Union{Missing, Float64}">Float64?</th></tr></thead><tbody><tr><th>1</th><td>1</td><td>biomass</td><td>WEC_SDGE</td><td>existing_cap_mw</td><td>21.4</td></tr><tr><th>2</th><td>2</td><td>hydroelectric_pumped_storage</td><td>WEC_SDGE</td><td>existing_cap_mw</td><td>42.0</td></tr><tr><th>3</th><td>3</td><td>natural_gas_fired_combined_cycle</td><td>WEC_SDGE</td><td>existing_cap_mw</td><td>607.1</td></tr><tr><th>4</th><td>4</td><td>natural_gas_fired_combined_cycle</td><td>WEC_SDGE</td><td>existing_cap_mw</td><td>30.0</td></tr><tr><th>5</th><td>5</td><td>natural_gas_fired_combined_cycle</td><td>WEC_SDGE</td><td>existing_cap_mw</td><td>49.4</td></tr><tr><th>6</th><td>1</td><td>biomass</td><td>WEC_SDGE</td><td>num_units</td><td>20.0</td></tr><tr><th>7</th><td>2</td><td>hydroelectric_pumped_storage</td><td>WEC_SDGE</td><td>num_units</td><td>2.0</td></tr><tr><th>8</th><td>3</td><td>natural_gas_fired_combined_cycle</td><td>WEC_SDGE</td><td>num_units</td><td>1.0</td></tr><tr><th>9</th><td>4</td><td>natural_gas_fired_combined_cycle</td><td>WEC_SDGE</td><td>num_units</td><td>1.0</td></tr><tr><th>10</th><td>5</td><td>natural_gas_fired_combined_cycle</td><td>WEC_SDGE</td><td>num_units</td><td>1.0</td></tr></tbody></table></div>



**Modify load and variable generation dataframes**

1. Convert from GMT to GMT-8

The load and variable generation files are in GMT. We want to convert to local California time. Here, we will ignore daylight savings and simply subtract 8 hours. (A more proper solution would involve the use of the [TimeZones](https://github.com/JuliaTime/TimeZones.jl/) package.)

To accomplish this, we 'wrap' the last 8 hours of the data series (i.e., the last 8 hours of the year) around to the front of the series.


```julia
gen_variable.hour = mod.(gen_variable.hour .- 9, 8760) .+ 1 
sort!(gen_variable, :hour)
loads.hour = mod.(loads.hour .- 9, 8760) .+ 1
sort!(loads, :hour);
```

2. Convert from "wide" to "long" format:

The file is in "wide" format, which has separate columns for values from different types of generators. This has more human readability, but often we want to have it in "long" format for computation.

"Long" format refers to a dataframe with a separate row entry for every value. The rest of the columns are "indexes" that uniquely identify the row. In Julia, we use the DataFrame package's [`stack`](https://juliadata.github.io/DataFrames.jl/stable/man/reshaping_and_pivoting/) function to accomplish this. (This is similar to the [`melt`](https://www.rdocumentation.org/packages/reshape2/versions/1.4.4/topics/melt) function for those familiar with DataFrames in R)


```julia
describe(gen_variable)
```




<div class="data-frame"><p>31 rows × 7 columns</p><table class="data-frame"><thead><tr><th></th><th>variable</th><th>mean</th><th>min</th><th>median</th><th>max</th><th>nmissing</th><th>eltype</th></tr><tr><th></th><th title="Symbol">Symbol</th><th title="Float64">Float64</th><th title="Real">Real</th><th title="Float64">Float64</th><th title="Real">Real</th><th title="Int64">Int64</th><th title="DataType">DataType</th></tr></thead><tbody><tr><th>1</th><td>hour</td><td>4380.5</td><td>1</td><td>4380.5</td><td>8760</td><td>0</td><td>Int64</td></tr><tr><th>2</th><td>wec_sdge_biomass_1.0</td><td>1.0</td><td>1</td><td>1.0</td><td>1</td><td>0</td><td>Int64</td></tr><tr><th>3</th><td>wec_sdge_hydroelectric_pumped_storage_1.0</td><td>1.0</td><td>1</td><td>1.0</td><td>1</td><td>0</td><td>Int64</td></tr><tr><th>4</th><td>wec_sdge_natural_gas_fired_combined_cycle_1.0</td><td>1.0</td><td>1</td><td>1.0</td><td>1</td><td>0</td><td>Int64</td></tr><tr><th>5</th><td>wec_sdge_natural_gas_fired_combined_cycle_2.0</td><td>1.0</td><td>1</td><td>1.0</td><td>1</td><td>0</td><td>Int64</td></tr><tr><th>6</th><td>wec_sdge_natural_gas_fired_combined_cycle_3.0</td><td>1.0</td><td>1</td><td>1.0</td><td>1</td><td>0</td><td>Int64</td></tr><tr><th>7</th><td>wec_sdge_natural_gas_fired_combined_cycle_4.0</td><td>1.0</td><td>1</td><td>1.0</td><td>1</td><td>0</td><td>Int64</td></tr><tr><th>8</th><td>wec_sdge_natural_gas_fired_combustion_turbine_1.0</td><td>1.0</td><td>1</td><td>1.0</td><td>1</td><td>0</td><td>Int64</td></tr><tr><th>9</th><td>wec_sdge_natural_gas_fired_combustion_turbine_2.0</td><td>1.0</td><td>1</td><td>1.0</td><td>1</td><td>0</td><td>Int64</td></tr><tr><th>10</th><td>wec_sdge_natural_gas_fired_combustion_turbine_3.0</td><td>1.0</td><td>1</td><td>1.0</td><td>1</td><td>0</td><td>Int64</td></tr><tr><th>11</th><td>wec_sdge_natural_gas_fired_combustion_turbine_4.0</td><td>1.0</td><td>1</td><td>1.0</td><td>1</td><td>0</td><td>Int64</td></tr><tr><th>12</th><td>wec_sdge_natural_gas_fired_combustion_turbine_5.0</td><td>1.0</td><td>1</td><td>1.0</td><td>1</td><td>0</td><td>Int64</td></tr><tr><th>13</th><td>wec_sdge_natural_gas_fired_combustion_turbine_6.0</td><td>1.0</td><td>1</td><td>1.0</td><td>1</td><td>0</td><td>Int64</td></tr><tr><th>14</th><td>wec_sdge_natural_gas_fired_combustion_turbine_7.0</td><td>1.0</td><td>1</td><td>1.0</td><td>1</td><td>0</td><td>Int64</td></tr><tr><th>15</th><td>wec_sdge_natural_gas_fired_combustion_turbine_8.0</td><td>1.0</td><td>1</td><td>1.0</td><td>1</td><td>0</td><td>Int64</td></tr><tr><th>16</th><td>wec_sdge_natural_gas_fired_combustion_turbine_9.0</td><td>1.0</td><td>1</td><td>1.0</td><td>1</td><td>0</td><td>Int64</td></tr><tr><th>17</th><td>wec_sdge_natural_gas_fired_combustion_turbine_10.0</td><td>1.0</td><td>1</td><td>1.0</td><td>1</td><td>0</td><td>Int64</td></tr><tr><th>18</th><td>wec_sdge_natural_gas_fired_combustion_turbine_11.0</td><td>1.0</td><td>1</td><td>1.0</td><td>1</td><td>0</td><td>Int64</td></tr><tr><th>19</th><td>wec_sdge_natural_gas_fired_combustion_turbine_12.0</td><td>1.0</td><td>1</td><td>1.0</td><td>1</td><td>0</td><td>Int64</td></tr><tr><th>20</th><td>wec_sdge_onshore_wind_turbine_1.0</td><td>0.253519</td><td>0.0</td><td>0.19465</td><td>0.9089</td><td>0</td><td>Float64</td></tr><tr><th>21</th><td>wec_sdge_small_hydroelectric_1.0</td><td>0.41567</td><td>0.2167</td><td>0.38475</td><td>0.6337</td><td>0</td><td>Float64</td></tr><tr><th>22</th><td>wec_sdge_solar_photovoltaic_1.0</td><td>0.247288</td><td>0.0</td><td>0.0</td><td>0.7959</td><td>0</td><td>Float64</td></tr><tr><th>23</th><td>wec_sdge_naturalgas_ccccsavgcf_0.0</td><td>1.0</td><td>1</td><td>1.0</td><td>1</td><td>0</td><td>Int64</td></tr><tr><th>24</th><td>wec_sdge_naturalgas_ccavgcf_0.0</td><td>1.0</td><td>1</td><td>1.0</td><td>1</td><td>0</td><td>Int64</td></tr><tr><th>25</th><td>wec_sdge_naturalgas_ctavgcf_0.0</td><td>1.0</td><td>1</td><td>1.0</td><td>1</td><td>0</td><td>Int64</td></tr><tr><th>26</th><td>wec_sdge_battery_0.0</td><td>1.0</td><td>1</td><td>1.0</td><td>1</td><td>0</td><td>Int64</td></tr><tr><th>27</th><td>wec_sdge_naturalgas_ccs100_0.0</td><td>1.0</td><td>1</td><td>1.0</td><td>1</td><td>0</td><td>Int64</td></tr><tr><th>28</th><td>wec_sdge_landbasedwind_ltrg1_1.0</td><td>0.341167</td><td>0.0</td><td>0.249565</td><td>0.93217</td><td>0</td><td>Float64</td></tr><tr><th>29</th><td>wec_sdge_landbasedwind_ltrg1_2.0</td><td>0.257347</td><td>0.0</td><td>0.194376</td><td>0.912162</td><td>0</td><td>Float64</td></tr><tr><th>30</th><td>wec_sdge_utilitypv_losangeles_1.0</td><td>0.346479</td><td>0.0</td><td>0.0</td><td>1.0</td><td>0</td><td>Float64</td></tr><tr><th>&vellip;</th><td>&vellip;</td><td>&vellip;</td><td>&vellip;</td><td>&vellip;</td><td>&vellip;</td><td>&vellip;</td><td>&vellip;</td></tr></tbody></table></div>




```julia
gen_variable_long = stack(gen_variable, 
                        Not(:hour), 
                        variable_name=:gen_full,
                        value_name=:cf);
# Now we have a "long" dataframe; 
# let's look at the first 6 entries of a wind resource for example
first(gen_variable_long[gen_variable_long.gen_full.=="wec_sdge_onshore_wind_turbine_1.0",:],6)
```




<div class="data-frame"><p>6 rows × 3 columns</p><table class="data-frame"><thead><tr><th></th><th>hour</th><th>gen_full</th><th>cf</th></tr><tr><th></th><th title="Int64">Int64</th><th title="String">String</th><th title="Float64">Float64</th></tr></thead><tbody><tr><th>1</th><td>1</td><td>wec_sdge_onshore_wind_turbine_1.0</td><td>0.1694</td></tr><tr><th>2</th><td>2</td><td>wec_sdge_onshore_wind_turbine_1.0</td><td>0.1658</td></tr><tr><th>3</th><td>3</td><td>wec_sdge_onshore_wind_turbine_1.0</td><td>0.1502</td></tr><tr><th>4</th><td>4</td><td>wec_sdge_onshore_wind_turbine_1.0</td><td>0.1455</td></tr><tr><th>5</th><td>5</td><td>wec_sdge_onshore_wind_turbine_1.0</td><td>0.0932</td></tr><tr><th>6</th><td>6</td><td>wec_sdge_onshore_wind_turbine_1.0</td><td>0.0717</td></tr></tbody></table></div>



**Extract single time from the data**
(we will get to the full set later)

- demand
- variable generation

Note: variable generation profiles are often encoded in terms of hourly maximum capacity factor, expressed as a per unit value or percentage of installed capacity. Hourly capacity factor is then later multiplied by installed capacity to yield a maximum hourly generation for each resource.


```julia
hr = 2416  # pick 4pm on a spring day
loads_single = loads[loads[:,:hour] .== hr, Not(:hour)];
var_cf_single = gen_variable_long[
    gen_variable_long.hour .== hr, 
    Not(:hour)]
```




<div class="data-frame"><p>30 rows × 2 columns</p><table class="data-frame"><thead><tr><th></th><th>gen_full</th><th>cf</th></tr><tr><th></th><th title="String">String</th><th title="Float64">Float64</th></tr></thead><tbody><tr><th>1</th><td>wec_sdge_biomass_1.0</td><td>1.0</td></tr><tr><th>2</th><td>wec_sdge_hydroelectric_pumped_storage_1.0</td><td>1.0</td></tr><tr><th>3</th><td>wec_sdge_natural_gas_fired_combined_cycle_1.0</td><td>1.0</td></tr><tr><th>4</th><td>wec_sdge_natural_gas_fired_combined_cycle_2.0</td><td>1.0</td></tr><tr><th>5</th><td>wec_sdge_natural_gas_fired_combined_cycle_3.0</td><td>1.0</td></tr><tr><th>6</th><td>wec_sdge_natural_gas_fired_combined_cycle_4.0</td><td>1.0</td></tr><tr><th>7</th><td>wec_sdge_natural_gas_fired_combustion_turbine_1.0</td><td>1.0</td></tr><tr><th>8</th><td>wec_sdge_natural_gas_fired_combustion_turbine_2.0</td><td>1.0</td></tr><tr><th>9</th><td>wec_sdge_natural_gas_fired_combustion_turbine_3.0</td><td>1.0</td></tr><tr><th>10</th><td>wec_sdge_natural_gas_fired_combustion_turbine_4.0</td><td>1.0</td></tr><tr><th>11</th><td>wec_sdge_natural_gas_fired_combustion_turbine_5.0</td><td>1.0</td></tr><tr><th>12</th><td>wec_sdge_natural_gas_fired_combustion_turbine_6.0</td><td>1.0</td></tr><tr><th>13</th><td>wec_sdge_natural_gas_fired_combustion_turbine_7.0</td><td>1.0</td></tr><tr><th>14</th><td>wec_sdge_natural_gas_fired_combustion_turbine_8.0</td><td>1.0</td></tr><tr><th>15</th><td>wec_sdge_natural_gas_fired_combustion_turbine_9.0</td><td>1.0</td></tr><tr><th>16</th><td>wec_sdge_natural_gas_fired_combustion_turbine_10.0</td><td>1.0</td></tr><tr><th>17</th><td>wec_sdge_natural_gas_fired_combustion_turbine_11.0</td><td>1.0</td></tr><tr><th>18</th><td>wec_sdge_natural_gas_fired_combustion_turbine_12.0</td><td>1.0</td></tr><tr><th>19</th><td>wec_sdge_onshore_wind_turbine_1.0</td><td>0.3006</td></tr><tr><th>20</th><td>wec_sdge_small_hydroelectric_1.0</td><td>0.4351</td></tr><tr><th>21</th><td>wec_sdge_solar_photovoltaic_1.0</td><td>0.5695</td></tr><tr><th>22</th><td>wec_sdge_naturalgas_ccccsavgcf_0.0</td><td>1.0</td></tr><tr><th>23</th><td>wec_sdge_naturalgas_ccavgcf_0.0</td><td>1.0</td></tr><tr><th>24</th><td>wec_sdge_naturalgas_ctavgcf_0.0</td><td>1.0</td></tr><tr><th>25</th><td>wec_sdge_battery_0.0</td><td>1.0</td></tr><tr><th>26</th><td>wec_sdge_naturalgas_ccs100_0.0</td><td>1.0</td></tr><tr><th>27</th><td>wec_sdge_landbasedwind_ltrg1_1.0</td><td>0.482802</td></tr><tr><th>28</th><td>wec_sdge_landbasedwind_ltrg1_2.0</td><td>0.302244</td></tr><tr><th>29</th><td>wec_sdge_utilitypv_losangeles_1.0</td><td>0.690915</td></tr><tr><th>30</th><td>wec_sdge_utilitypv_losangeles_2.0</td><td>0.743176</td></tr></tbody></table></div>



### 3. Create solver function

Here, we introduce a functionalized version of the optimization problem. We pass the problem parameters to the function, and it constructs and solves the model, returning the result. This will be useful when we wish to call the same basic model multiple times, e.g. for conducting sensitivities on parameters (among other things).


```julia
#=
Function to solve economic dispatch problem (single-time period, single-zone)
Inputs:
    gen_df -- dataframe with generator info
    loads  -- dataframe with load info
    gen_variable -- capacity factors of variable generators (in "long" format)
Note: it is always a good idea to include a comment blog describing your
function's inputs clearly!
=#
function economic_dispatch_single(gen_df, loads, gen_variable)
    ED = Model(GLPK.Optimizer) # You could use Clp as well, with Clp.Optimizer
    
    # Define sets based on data
      # A set of all variable generators
    G_var = gen_df[gen_df[!,:is_variable] .== 1,:r_id] 
      # A set of all non-variable generators
    G_nonvar = gen_df[gen_df[!,:is_variable] .== 0,:r_id]
      # Set of all generators
    G = gen_df.r_id
    # Extract some parameters given the input data
      # Generator capacity factor time series for variable generators
    gen_var_cf = innerjoin(gen_variable, 
                    gen_df[gen_df.is_variable .== 1 , 
                        [:r_id, :gen_full, :existing_cap_mw]], 
                    on = :gen_full)
        
    # Decision variables   
    @variables(ED, begin
        GEN[G]  >= 0     # generation
        # Note: we assume Pmin = 0 for all resources for simplicty here
    end)
                
    # Objective function
    @objective(ED, Min, 
        sum( (gen_df[i,:heat_rate_mmbtu_per_mwh] * gen_df[i,:fuel_cost] +
            gen_df[i,:var_om_cost_per_mwh]) * GEN[i] 
                        for i in G_nonvar) + 
        sum(gen_df[i,:var_om_cost_per_mwh] * GEN[i] 
                        for i in G_var)
    )

    # Demand constraint
    @constraint(ED, cDemand, 
        sum(GEN[i] for i in G) == loads[1,:demand])

    # Capacity constraint (non-variable generation)
    for i in G_nonvar
        @constraint(ED, GEN[i] <= gen_df[i,:existing_cap_mw])
    end

    # Variable generation capacity constraint
    for i in 1:nrow(gen_var_cf)
        @constraint(ED, GEN[gen_var_cf[i,:r_id] ] <= 
                        gen_var_cf[i,:cf] *
                        gen_var_cf[i,:existing_cap_mw])
    end

    # Solve statement (! indicates runs in place)
    optimize!(ED)

    # Dataframe of optimal decision variables
    solution = DataFrame(
        r_id = gen_df.r_id,
        resource = gen_df.resource,
        gen = value.(GEN).data
        )

    # Return the solution and objective as named tuple
    return (
        solution = solution, 
        cost = objective_value(ED),
    )
end
```




    economic_dispatch_single (generic function with 1 method)



### 4. Solve and print data


```julia
solution = economic_dispatch_single(gen_df, loads_single, var_cf_single);
solution.solution
```




<div class="data-frame"><p>25 rows × 3 columns</p><table class="data-frame"><thead><tr><th></th><th>r_id</th><th>resource</th><th>gen</th></tr><tr><th></th><th title="Union{Missing, Int64}">Int64?</th><th title="Union{Missing, String}">String?</th><th title="Float64">Float64</th></tr></thead><tbody><tr><th>1</th><td>1</td><td>biomass</td><td>21.4</td></tr><tr><th>2</th><td>2</td><td>hydroelectric_pumped_storage</td><td>42.0</td></tr><tr><th>3</th><td>3</td><td>natural_gas_fired_combined_cycle</td><td>607.1</td></tr><tr><th>4</th><td>4</td><td>natural_gas_fired_combined_cycle</td><td>30.0</td></tr><tr><th>5</th><td>5</td><td>natural_gas_fired_combined_cycle</td><td>49.4</td></tr><tr><th>6</th><td>6</td><td>natural_gas_fired_combined_cycle</td><td>570.0</td></tr><tr><th>7</th><td>7</td><td>natural_gas_fired_combined_cycle</td><td>336.0</td></tr><tr><th>8</th><td>8</td><td>natural_gas_fired_combined_cycle</td><td>336.0</td></tr><tr><th>9</th><td>9</td><td>natural_gas_fired_combustion_turbine</td><td>0.0</td></tr><tr><th>10</th><td>10</td><td>natural_gas_fired_combustion_turbine</td><td>49.9</td></tr><tr><th>11</th><td>11</td><td>natural_gas_fired_combustion_turbine</td><td>0.0</td></tr><tr><th>12</th><td>12</td><td>natural_gas_fired_combustion_turbine</td><td>0.0</td></tr><tr><th>13</th><td>13</td><td>natural_gas_fired_combustion_turbine</td><td>86.8545</td></tr><tr><th>14</th><td>14</td><td>natural_gas_fired_combustion_turbine</td><td>0.0</td></tr><tr><th>15</th><td>15</td><td>natural_gas_fired_combustion_turbine</td><td>0.0</td></tr><tr><th>16</th><td>16</td><td>natural_gas_fired_combustion_turbine</td><td>0.0</td></tr><tr><th>17</th><td>17</td><td>natural_gas_fired_combustion_turbine</td><td>0.0</td></tr><tr><th>18</th><td>18</td><td>natural_gas_fired_combustion_turbine</td><td>0.0</td></tr><tr><th>19</th><td>19</td><td>natural_gas_fired_combustion_turbine</td><td>0.0</td></tr><tr><th>20</th><td>20</td><td>natural_gas_fired_combustion_turbine</td><td>0.0</td></tr><tr><th>21</th><td>21</td><td>natural_gas_fired_combustion_turbine</td><td>0.0</td></tr><tr><th>22</th><td>22</td><td>natural_gas_fired_combustion_turbine</td><td>0.0</td></tr><tr><th>23</th><td>23</td><td>onshore_wind_turbine</td><td>65.4106</td></tr><tr><th>24</th><td>24</td><td>small_hydroelectric</td><td>3.18493</td></tr><tr><th>25</th><td>25</td><td>solar_photovoltaic</td><td>284.75</td></tr></tbody></table></div>



We will now plot the **supply curve** (also known as the **dispatch stack**) for this single period to illustrate the concept of **merit order**. First, let's create a dataframe with the variable cost as a column, and then sort on it from lowest to highest.


```julia
supply_curve = leftjoin(gen_df,
                var_cf_single, 
                on = :gen_full)
supply_curve[!, :varcost] .= 0.0
supply_curve[!, :cap] .= 0.0

# Store varcost of non-variable generators:
I = supply_curve[!,:is_variable] .== 0  # `I` contains indexes to the non-variable generators
supply_curve[I,:varcost] .= 
    supply_curve[I,:heat_rate_mmbtu_per_mwh] .* 
        supply_curve[I,:fuel_cost] .+
    supply_curve[I,:var_om_cost_per_mwh]

# Calculate available capacity for each generator
supply_curve[I,:cap] = supply_curve[I,:existing_cap_mw]

# Store varcost of variable generators (in this case, 0)
I = (supply_curve[:,:is_variable] .== 1) # `I` contains indexes to the variable generators
supply_curve[I,:varcost] = supply_curve[I,:var_om_cost_per_mwh]

# Calculate available capacity for each generator
# (adjusted for variable generation)
supply_curve[I,:cap] = supply_curve[I,:existing_cap_mw] .* supply_curve[I,:cf]

sort!(supply_curve, :varcost);
```

Now, let's plot this as a supply curve where:

- width = capacity of the generator (i.e., $Pmax$)
- height = variable cost of the generator ($VarCost$)


```julia
rectangle(w, h, x, y) = Shape(x .+ [0,w,w,0], y .+ [0,0,h,h])

p = plot(nrow(supply_curve))

marg_gen = 0
marg_price = 0
x = 0
for i in 1:nrow(supply_curve)
    plot!(p, rectangle(supply_curve[i,:cap],
            supply_curve[i,:varcost],x,0), opacity=.5, 
            label = "")
    if (x < loads_single[1,:demand]) && (x + supply_curve[i,:cap] > loads_single[1,:demand])
        marg_gen = i
        marg_price = supply_curve[i,:varcost]
    end
    x = x + supply_curve[i,:cap]
end 
vline!([loads_single[1,:demand]],linecolor="black",linewidth=5,
label = "demand")
hline!([marg_price],linecolor="blue",linewidth=5,
label = "SRMC")
title!("Supply and demand curves")
xlabel!("Capacity")
ylabel!("Marginal cost")
p
```




    <div id="588be611-416d-495f-b008-f9a578204c1e" style="width:600px;height:400px;"></div>
    




In the above, everything to the left of demand is dispatched at full capacity, except for the last generator (the "marginal generator") where demand stops. Everything to the right is not dispatched (=0). We stored the marginal generator in `marg_gen`:


```julia
supply_curve[marg_gen,:]
```




<p>DataFrameRow (34 columns)</p><div class="data-frame"><table class="data-frame"><thead><tr><th></th><th>r_id</th><th>resource</th><th>region</th><th>existing_cap_mw</th><th>num_units</th><th>cap_size</th><th>var_om_cost_per_mwh</th><th>var_om_cost_per_mwh_in</th></tr><tr><th></th><th title="Union{Missing, Int64}">Int64?</th><th title="Union{Missing, String}">String?</th><th title="Union{Missing, String}">String?</th><th title="Union{Missing, Float64}">Float64?</th><th title="Union{Missing, Int64}">Int64?</th><th title="Union{Missing, Float64}">Float64?</th><th title="Union{Missing, Float64}">Float64?</th><th title="Union{Missing, Int64}">Int64?</th></tr></thead><tbody><tr><th>13</th><td>13</td><td>natural_gas_fired_combustion_turbine</td><td>WEC_SDGE</td><td>96.6</td><td>3</td><td>32.2</td><td>10.8</td><td>0</td></tr></tbody></table></div>



Which is generating $GEN_{marg\_gen}=$


```julia
solution.solution[solution.solution.r_id .== supply_curve[marg_gen,:r_id],
    :gen][1]
```




    86.85450800000007



**Combine by generation type**

We are often interested in the total generation *by type*, not individual generators. We use `combine` and `groupby` methods to sum over resource:


```julia
sol_gen = combine(groupby(solution.solution, :resource), 
            :gen => sum)
# We use the PrettyTable package here to round numbers when displaying
pretty_table(sol_gen, formatters = ft_printf("%5.3f")) 
```

    ┌──────────────────────────────────────┬──────────┐
    │[1m                             resource [0m│[1m  gen_sum [0m│
    │[90m                              String? [0m│[90m  Float64 [0m│
    ├──────────────────────────────────────┼──────────┤
    │                              biomass │   21.400 │
    │         hydroelectric_pumped_storage │   42.000 │
    │     natural_gas_fired_combined_cycle │ 1928.500 │
    │ natural_gas_fired_combustion_turbine │  136.755 │
    │                 onshore_wind_turbine │   65.411 │
    │                  small_hydroelectric │    3.185 │
    │                   solar_photovoltaic │  284.750 │
    └──────────────────────────────────────┴──────────┘


Finally, let's add back the 600MW of behind-the-meter solar to get a full picture of all the electricity on the grid.


```julia
# we use the copy function to create a copy of the dataframe (as opposed to working with the original)
sol_gen_btm = copy(sol_gen)

# We'll use the same variability profile for BTM solar as the utility-scale solar here (in reality, they would probably differ)
btm = DataFrame(resource = ["solar_photovoltaic_btm"], gen_sum = var_cf_single[var_cf_single.gen_full .== "wec_sdge_solar_photovoltaic_1.0",:cf] * 600)
append!(sol_gen_btm, btm)
pretty_table(sol_gen_btm, formatters = ft_printf("%5.3f")) # We use the PrettyTable package to round numbers when displaying
```

    ┌──────────────────────────────────────┬──────────┐
    │[1m                             resource [0m│[1m  gen_sum [0m│
    │[90m                              String? [0m│[90m  Float64 [0m│
    ├──────────────────────────────────────┼──────────┤
    │                              biomass │   21.400 │
    │         hydroelectric_pumped_storage │   42.000 │
    │     natural_gas_fired_combined_cycle │ 1928.500 │
    │ natural_gas_fired_combustion_turbine │  136.755 │
    │                 onshore_wind_turbine │   65.411 │
    │                  small_hydroelectric │    3.185 │
    │                   solar_photovoltaic │  284.750 │
    │               solar_photovoltaic_btm │  341.700 │
    └──────────────────────────────────────┴──────────┘


Hence, the majority of energy is being met by combined cycle natural gas turbines (CCGT), second by solar PV, and third by behind-the-meter solar.

## Multiple-time period, simple generator constraints

We neglected the time dimension above and now we will add it back in. The key changes to the model arise from adding an additional index for time:

$$
\begin{align}
\min \ & \sum_{g \in G, t \in T} VarCost_g \times GEN_{g,t} & \\
\text{s.t.} & \\
 & \sum_{g} GEN_{g,t} = Demand_t & \forall \quad t \in T \\
 & GEN_{g,t} \leq Pmax_{g,t} & \forall \quad g \in G , t \in T \\
 & GEN_{g,t} \geq Pmin_{g,t} & \forall \quad g \in G , t \in T
\end{align}
$$


Note that for conventional generators, we often have $Pmax_{g,t} = Pmax_{g}$ and $Pmin_{g,t} = Pmin_{g}$ constant over time. (We will consider unit commitment later, which will alter this logic.)

For variable renewable generators, it is natural that $Pmax_{g,t}$ varies with time (i.e., based on solar irradiation or wind speeds). $Pmin_{g,t}$ for wind and solar resources is usually 0, but for hydropower resources, minimum streamflow constraints can produce a time-variant parameter.

### 3. Create solver function

(We reuse steps 1 and 2 above to load packages and data.)


```julia
#=
Function to solve economic dispatch problem (multi-time period, single-zone)
Inputs:
    gen_df -- dataframe with generator info
    loads  -- load by time
    gen_variable -- capacity factors of variable generators (in "long" format)
=#
function economic_dispatch_multi(gen_df, loads, gen_variable)
    ED = Model(GLPK.Optimizer)

    # Define sets based on data
    G_var = gen_df[gen_df[!,:is_variable] .== 1,:r_id] 
    G_nonvar = gen_df[gen_df[!,:is_variable] .== 0,:r_id]
    G = gen_df.r_id
    T = loads.hour

    # Generator capacity factor time series for variable generators
    gen_var_cf = innerjoin(gen_variable, 
                    gen_df[gen_df.is_variable .== 1 , 
                        [:r_id, :gen_full, :existing_cap_mw]], 
                    on = :gen_full)
        
    # Decision variables   
    @variables(ED, begin
        GEN[G, T]  >= 0     # generation
        # Note: we assume Pmin = 0 for all resources for simplicty here
    end)
                
    # Objective function
    @objective(ED, Min, 
        sum( (gen_df[gen_df.r_id .== i,:heat_rate_mmbtu_per_mwh][1] * gen_df[gen_df.r_id .== i,:fuel_cost][1] +
            gen_df[gen_df.r_id .== i,:var_om_cost_per_mwh][1]) * GEN[i,t] 
                        for i in G_nonvar for t in T) + 
        sum(gen_df[gen_df.r_id .== i,:var_om_cost_per_mwh][1] * GEN[i,t] 
                        for i in G_var for t in T)
    )
    
    # Demand constraint
    @constraint(ED, cDemand[t in T], 
        sum(GEN[i,t] for i in G) == loads[loads.hour .== t,:demand][1])

    # Capacity constraints (non-variable generation)
    @constraint(ED, Cap_nonvar[i in G_nonvar, t in T], 
        GEN[i,t] <= gen_df[gen_df.r_id .== i,:existing_cap_mw][1])

    # Variable generation capacity constraints
    @constraint(ED, Cap_var[i in 1:nrow(gen_var_cf)], 
            GEN[gen_var_cf[i,:r_id], gen_var_cf[i,:hour] ] <= 
                        gen_var_cf[i,:cf] *
                        gen_var_cf[i,:existing_cap_mw])
    
    # Solve statement (! indicates runs in place)
    optimize!(ED)

    # Dataframe of optimal decision variables
    # This requires a few extra steps because we are indexing decision variables over two sets: generator and time
    solution = DataFrame(value.(GEN).data, :auto)  # :auto argument asks DataFrames to pick the variable names automatically
    ax1 = value.(GEN).axes[1]
    ax2 = value.(GEN).axes[2]
    cols = names(solution)
    insertcols!(solution, 1, :r_id => ax1)
    solution = stack(solution, Not(:r_id), variable_name=:hour)
    solution.hour = foldl(replace, [cols[i] => ax2[i] for i in 1:length(ax2)], init=solution.hour)
    rename!(solution, :value => :gen)
    solution.hour = convert.(Int64,solution.hour)
    
    # Return the solution and objective as named tuple
    return (
        solution = solution, 
        cost = objective_value(ED),
        status = termination_status(ED)
    )

end
```




    economic_dispatch_multi (generic function with 1 method)



### 4. Solve a day's economic dispatch

We will subset our year-long dataset to a single spring day and solve.


```julia
n=100
T_period = (n*24+1):((n+1)*24)

loads_multi = loads[in.(loads.hour,Ref(T_period)),:]
gen_variable_multi = gen_variable_long[in.(gen_variable_long.hour,Ref(T_period)),:]

gen_df_sens = copy(gen_df)
gen_df_sens[gen_df_sens.resource .== "solar_photovoltaic",
    :existing_cap_mw] .= 3500

solution = economic_dispatch_multi(gen_df_sens, loads_multi, gen_variable_multi);
```


```julia
# **Add back in BTM solar and plot results**

sol_gen = innerjoin(solution.solution, 
                    gen_df[!, [:r_id, :resource]], 
                    on = :r_id)
sol_gen = combine(groupby(sol_gen, [:resource, :hour]), 
            :gen => sum)
sol_gen_btm = copy(sol_gen)
sol_gen_btm[sol_gen_btm.resource .== "solar_photovoltaic", :resource] .= "_solar_photovoltaic"
sol_gen_btm[sol_gen_btm.resource .== "onshore_wind_turbine", :resource] .= "_onshore_wind_turbine"
sol_gen_btm[sol_gen_btm.resource .== "small_hydroelectric", :resource] .= "_small_hydroelectric"

btm = DataFrame(resource = repeat(["_solar_photovoltaic_btm"]; outer=length(T_period)), 
    hour = T_period,
    gen_sum = gen_variable_multi[gen_variable_multi.gen_full .== "wec_sdge_solar_photovoltaic_1.0",:cf] * 600)
append!(sol_gen_btm, btm)

# Reorder resources
sort!(sol_gen_btm, [:hour, :resource])
CSV.write("results/sol_gen_btm.csv",sol_gen_btm)

sol_gen_sum = combine(groupby(sol_gen_btm, [:resource]), 
            :gen_sum => sum)
```




<div class="data-frame"><p>8 rows × 2 columns</p><table class="data-frame"><thead><tr><th></th><th>resource</th><th>gen_sum_sum</th></tr><tr><th></th><th title="Union{Missing, String}">String?</th><th title="Float64">Float64</th></tr></thead><tbody><tr><th>1</th><td>_onshore_wind_turbine</td><td>2314.42</td></tr><tr><th>2</th><td>_small_hydroelectric</td><td>76.4384</td></tr><tr><th>3</th><td>_solar_photovoltaic</td><td>20546.0</td></tr><tr><th>4</th><td>_solar_photovoltaic_btm</td><td>3522.18</td></tr><tr><th>5</th><td>biomass</td><td>513.6</td></tr><tr><th>6</th><td>hydroelectric_pumped_storage</td><td>1008.0</td></tr><tr><th>7</th><td>natural_gas_fired_combined_cycle</td><td>29749.6</td></tr><tr><th>8</th><td>natural_gas_fired_combustion_turbine</td><td>1224.89</td></tr></tbody></table></div>



**Plot charts**

We use the [VegaLite](https://www.queryverse.org/VegaLite.jl/stable/) library to make some fancier plots.

We also introduce the pipe operator `|>`, a helpful syntax that "pipes" the data on the left into the first argument of the function on the right. (Some more examples are [here](https://syl1.gitbook.io/julia-language-a-concise-tutorial/useful-packages/pipe).)


```julia
sol_gen_btm |>
@vlplot(:line, 
    x=:hour, y=:gen_sum, 
    color={"resource:n", scale={scheme="category20"}})

```


```julia
g = sol_gen_btm |>
@vlplot(:area, 
    x=:hour, y={:gen_sum, stack=:zero}, 
    color={"resource:n", scale={scheme="category20"}},
    width=600,height=600)
g |> FileIO.save("results/ED_multitime.png")
```

For reference, here is the load shape (after subtracting BTM solar) we are trying solve:


```julia
loads_multi |>
@vlplot(:line, 
    x=:hour, y=:demand)
```

## Multiple-time period, complex generator constraints with time coupling

The above problem did not include more than one time index per constraint. Hence, we could actually solve this each time period separately without regard to what is happening before or after. (These are known as "separable" problems.)

We now introduce time coupling into the model, by considering ramp rates, which limit how much generators can change output from one period to the next (an engineering constraint).

The new model:

$$
\begin{align}
\min \ & \sum_{g \in G, t \in T} VarCost_g \times GEN_{g,t} & \\
\text{s.t.} & \\
 & \sum_{g} GEN_{g,t} = Demand_t & \forall \quad t \in T \\
 & GEN_{g,t} \leq Pmax_{g,t} & \forall \quad g \in G , t \in T \\
 & GEN_{g,t} \geq Pmin_{g,t} & \forall \quad g \in G , t \in T \\
 & GEN_{g,t+1} - GEN_{g,t} \leq RampUp_{g} & \forall \quad g \in G , t = 1..T-1 \\
  & GEN_{g,t} - GEN_{g,t+1} \leq RampDn_{g} & \forall \quad g \in G , t = 1..T-1
\end{align}
$$



### 3. Create solver function

(We reuse steps 1 and 2 above to load packages and data.)

First, we define a function to convert two-dimensional outputs from JuMP into easy-to-read dataframes.


```julia
#=
Function to convert JuMP outputs (technically, AxisArrays) with two-indexes to a dataframe
Inputs:
    var -- JuMP AxisArray (e.g., value.(GEN))
Reference: https://jump.dev/JuMP.jl/v0.19/containers/
=#
function value_to_df_2dim(var)
    solution = DataFrame(var.data, :auto)
    ax1 = var.axes[1]
    ax2 = var.axes[2]
    cols = names(solution)
    insertcols!(solution, 1, :r_id => ax1)
    solution = stack(solution, Not(:r_id), variable_name=:hour)
    solution.hour = foldl(replace, [cols[i] => ax2[i] for i in 1:length(ax2)], init=solution.hour)
    rename!(solution, :value => :gen)
    solution.hour = convert.(Int64,solution.hour)
    return solution
end
```




    value_to_df_2dim (generic function with 1 method)




```julia
#=
Function to solve economic dispatch problem (multi-time period, single-zone)
Inputs:
    gen_df -- dataframe with generator info
    loads  -- load by time
    gen_variable -- capacity factors of variable generators (in "long" format)
=#
function economic_dispatch_multi_time(gen_df, loads, gen_variable)
    ED = Model(GLPK.Optimizer)

    # Define sets based on data
    G_var = gen_df[gen_df[!,:is_variable] .== 1,:r_id] 
    G_nonvar = gen_df[gen_df[!,:is_variable] .== 0,:r_id]
    G = gen_df.r_id
    T = loads.hour
    T_red = loads.hour[1:end-1]  # time periods used for ramp constraints
    
    # Generator capacity factor time series for variable generators
    gen_var_cf = innerjoin(gen_variable, 
                    gen_df[gen_df.is_variable .== 1 , 
                        [:r_id, :gen_full, :existing_cap_mw]], 
                    on = :gen_full)
        
    # Decision variables   
    @variables(ED, begin
        GEN[G, T]  >= 0     # generation
        # Note: we assume Pmin = 0 for all resources for simplicty here
    end)
                
    # Objective function
    @objective(ED, Min, 
        sum( (gen_df[gen_df.r_id .== i,:heat_rate_mmbtu_per_mwh][1] * gen_df[gen_df.r_id .== i,:fuel_cost][1] +
            gen_df[gen_df.r_id .== i,:var_om_cost_per_mwh][1]) * GEN[i,t] 
                        for i in G_nonvar for t in T) + 
        sum(gen_df[gen_df.r_id .== i,:var_om_cost_per_mwh][1] * GEN[i,t] 
                        for i in G_var for t in T)
    )
    
    # Demand constraint
    @constraint(ED, cDemand[t in T], 
        sum(GEN[i,t] for i in G) == loads[loads.hour .== t,:demand][1])

    # Capacity constraints (non-variable generation)
    @constraint(ED, Cap_nonvar[i in G_nonvar, t in T], 
        GEN[i,t] <= gen_df[gen_df.r_id .== i,:existing_cap_mw][1])

    # Variable generation capacity constraints
    @constraint(ED, Cap_var[i in 1:nrow(gen_var_cf)], 
            GEN[gen_var_cf[i,:r_id], gen_var_cf[i,:hour] ] <= 
                        gen_var_cf[i,:cf] *
                        gen_var_cf[i,:existing_cap_mw])
    
    # Ramp up constraints
    @constraint(ED, RampUp[i in G, t in T_red], 
        GEN[i,t+1] - GEN[i,t] <= gen_df[gen_df.r_id .== i,:existing_cap_mw][1] * 
                                 gen_df[gen_df.r_id .== i,:ramp_up_percentage][1] )

    # Ramp down constraints
    @constraint(ED, RampDn[i in G, t in T_red], 
        GEN[i,t] - GEN[i,t+1] <= gen_df[gen_df.r_id .== i,:existing_cap_mw][1] * 
                                 gen_df[gen_df.r_id .== i,:ramp_dn_percentage][1] )

    
    # Solve statement (! indicates runs in place)
    optimize!(ED)

    # Dataframe of optimal decision variables
    solution = value_to_df_2dim(value.(GEN))
    
    # Dataframe of ramp up constraint marginal values 
    # (aka "dual values" of the ramp constraint)
    # (to be used in solar sensitivity case)
    marginals = value_to_df_2dim(dual.(RampUp))
    
    # Return the solution, objective, and marginals as named tuple
    return (
        solution = solution, 
        cost = objective_value(ED),
        marginals = marginals
    )

end
```




    economic_dispatch_multi_time (generic function with 1 method)



### 4. Solve a day's economic dispatch


```julia
solution = economic_dispatch_multi_time(gen_df, loads_multi, gen_variable_multi);

sol_gen = innerjoin(solution.solution, 
                    gen_df[!, [:r_id, :resource]], 
                    on = :r_id)
sol_gen = combine(groupby(sol_gen, [:resource, :hour]), 
            :gen => sum)
sol_gen_btm = sol_gen
sol_gen_btm[sol_gen_btm.resource .== "solar_photovoltaic", :resource] .= "_solar_photovoltaic"
sol_gen_btm[sol_gen_btm.resource .== "onshore_wind_turbine", :resource] .= "_onshore_wind_turbine"
sol_gen_btm[sol_gen_btm.resource .== "small_hydroelectric", :resource] .= "_small_hydroelectric"

btm = DataFrame(resource = repeat(["_solar_photovoltaic_btm"]; outer=length(T_period)), 
    hour = T_period,
    gen_sum = gen_variable_multi[gen_variable_multi.gen_full .== "wec_sdge_solar_photovoltaic_1.0",:cf] * 600)
append!(sol_gen_btm, btm)

# Reorder resources
sort!(sol_gen_btm, [:hour, :resource])

sol_gen_btm |>
@vlplot(:area, 
    x=:hour, y={:gen_sum, stack=:zero}, 
    color={"resource:n", scale={scheme="category20"}},
    width=600,height=600) |> 
    FileIO.save("results/ED_moderate.png")
```

### 5. Run sensitivity on solar capacity

Let's now examine what happens if we increase utility solar capacity from the current 500MW up to 3,500MW.


```julia
# Increase solar photovoltaic capacity
gen_df_sens = copy(gen_df)
gen_df_sens[gen_df_sens.resource .== "solar_photovoltaic",
    :existing_cap_mw] .= 3500

solution = economic_dispatch_multi_time(gen_df_sens, loads_multi, gen_variable_multi);

sol_gen = innerjoin(solution.solution, 
                    gen_df[!, [:r_id, :resource]], 
                    on = :r_id)

sol_gen = combine(groupby(sol_gen, [:resource, :hour]), 
            :gen => sum)
sol_gen_btm = copy(sol_gen)
sol_gen_btm[sol_gen_btm.resource .== "solar_photovoltaic", :resource] .= "_solar_photovoltaic"
sol_gen_btm[sol_gen_btm.resource .== "onshore_wind_turbine", :resource] .= "_onshore_wind_turbine"
sol_gen_btm[sol_gen_btm.resource .== "small_hydroelectric", :resource] .= "_small_hydroelectric"

btm = DataFrame(resource = repeat(["_solar_photovoltaic_btm"]; outer=length(T_period)), 
    hour = T_period,
    gen_sum = gen_variable_multi[gen_variable_multi.gen_full .== "wec_sdge_solar_photovoltaic_1.0",:cf] * 600)
append!(sol_gen_btm, btm)

# Reorder resources
sort!(sol_gen_btm, [:hour, :resource])

# Plot and save the resulting generation figure
p = sol_gen_btm |>
@vlplot(:area, 
    x=:hour, y={:gen_sum, stack=:zero}, 
    color={"resource:n", scale={scheme="category20"}},
    width=600,height=600)
save("results/sol_gen_btm.png",p)
p
```

The large reduction in net load during the day and the loss of solar (orange) in the late afternoon timed with the rise in the evening peak is sometimes called the **"duck curve"** because of their passing resemblance to a duck's belly and head. Squint at it and then examine this (https://www.energy.gov/eere/articles/confronting-duck-curve-how-address-over-generation-solar-energy):



You may notice that the model still uses combustion turbines (light pink above) during hours when there is sufficient capacity of combined cycle plants, which are more efficient. This is because of ramp limits, which we will examine now.

All constraints (in a linear program) have a **"marginal value"** (also known as "shadow price" or "dual value") which is how much the objective would change if the right-hand-side (RHS) of the constraint were increased by one unit. Let's examine the marginals of the ramp up constraint for the combined cycle plants using the values from the new equation we added to `economic_dispatch_multi_time` (line 72).


```julia
marg_gen = innerjoin(solution.marginals, 
                    gen_df[!, [:r_id, :resource]], 
                    on = :r_id)
marg_gen[marg_gen.resource .== "natural_gas_fired_combined_cycle",:] |>
@vlplot(:line, 
    x=:hour, y={:gen, title="Marginal"}, 
    color={"r_id:n", scale={scheme="category20"}})

```

Thus, during the ramp hours, the marginal value becomes negative&mdash;in this case, reaching ~ -9 for the combined cycle plant with `r_id=5`. This is the interpretation:

**If the RHS of the RampUp constraint for `r_id=5` combined cycle plant were increased by 1 unit, the problem objective would reduce by about 9.** 

The precise value of the largest magnitude ramp marginal is:


```julia
marg_gen[argmin(marg_gen.gen),:]
```




<p>DataFrameRow (4 columns)</p><div class="data-frame"><table class="data-frame"><thead><tr><th></th><th>r_id</th><th>hour</th><th>gen</th><th>resource</th></tr><tr><th></th><th title="Union{Missing, Int64}">Int64?</th><th title="Int64">Int64</th><th title="Float64">Float64</th><th title="Union{Missing, String}">String?</th></tr></thead><tbody><tr><th>430</th><td>5</td><td>2418</td><td>-8.9672</td><td>natural_gas_fired_combined_cycle</td></tr></tbody></table></div>



**Why is this the case?** Due to the ramp constraint, we have to substitute some generation from `r_id=5` with a more expensive unit, in this case the more flexible `r_id=9` combustion turbine. 

We can, in fact, recover this marginal value by re-examining the variable costs of these generators at the top of the notebook:


```julia
varcost_comp = supply_curve[in.(supply_curve.r_id,Ref([5,9])),[:r_id,:varcost]]
varcost_comp
```




<div class="data-frame"><p>2 rows × 2 columns</p><table class="data-frame"><thead><tr><th></th><th>r_id</th><th>varcost</th></tr><tr><th></th><th title="Union{Missing, Int64}">Int64?</th><th title="Float64">Float64</th></tr></thead><tbody><tr><th>1</th><td>5</td><td>29.3061</td></tr><tr><th>2</th><td>9</td><td>38.2733</td></tr></tbody></table></div>



If our guess is true, the difference between these should be our marginal value, and voila:


```julia
varcost_comp[2,2] - varcost_comp[1,2]
```




    8.967199999999998



Also notice that the marginals are zero for most of the day. Recalling from Notebook 1 our discussion on "binding" constraints, this is because the RampUp constraint is **not binding** during these hours. Alternatively, changing the RHS of the RampUp constraint does not affect the problem objective.

10fa612bf113d41cb448951b022e6e5edd82a222
ipynb
Jupyter Notebook
Notebooks/04-Economic-Dispatch.ipynb
tyson-dao/power-systems-optimization
e35fb71be414fd115588a85e3a94ea3261af549a
CC-BY-4.0 MIT
2020-09-12T00:21:20.000Z
2022-02-28T12:57:28.000Z
Notebooks/04-Economic-Dispatch.ipynb
tyson-dao/power-systems-optimization
e35fb71be414fd115588a85e3a94ea3261af549a
CC-BY-4.0 MIT
2020-09-10T13:50:07.000Z
2021-09-21T20:53:49.000Z
Notebooks/04-Economic-Dispatch.ipynb
tyson-dao/power-systems-optimization
e35fb71be414fd115588a85e3a94ea3261af549a
CC-BY-4.0 MIT
2020-09-02T20:41:27.000Z
2022-03-25T14:56:39.000Z
Qwen/Qwen-72B
1. YES
2. YES