Incorporate a model of demand flexibility

README.md

@@ -1,7 +1,7 @@
 # REISE.jl
 Renewable Energy Integration Simulation Engine.
 
-This repository contains, in the **src** folder, the Julia scripts to run the power-flow study in the U.S. electric grid. The simulation engine relies on [Gurobi] as the optimization solver.
+This repository contains, in the `src` folder, the Julia scripts to run the power-flow study in the U.S. electric grid. The simulation engine relies on the use of an external optimization solver; any solver [compatible with JuMP] can be used, though code that makes efficient use of the capabilities of [Gurobi] is included in this repository.
 
 ## Table of Contents
 1. [Dependencies](#dependencies)
@@ -43,8 +43,8 @@ If Gurobi is to be used as the solver, this will need to be installed as well:
 - [Gurobi Installation Guide]
 
 The package itself has two components that require installation:
-- [`Julia` package](#julia-package-installation) to run the simulation
-- optional [`python` scripts](#python-requirements-installation) for some additional pre- and post-processing
+- [Julia package](#julia-package-installation) to run the simulation
+- optional [Python scripts](#python-requirements-installation) for some additional pre- and post-processing
 
 Instead of installing locally, this package can also be used with the included [Docker](#Docker) image. 
 
@@ -216,7 +216,7 @@ pkg> activate /PATH/TO/REISE.jl
 
 (REISE) pkg> instantiate
 ```
-In order for the below `python` scripts to use this environment, set a
+In order for the below Python scripts to use this environment, set a
 `JULIA_PROJECT` environment variable to the path to `/PATH/TO/REISE.jl`.
 
 For more information about the different installation options, please see the
@@ -233,15 +233,15 @@ using REISE
 ### Python Packages
 
 #### Python Requirements Installation
-The dependencies of the python scripts contained in `pyreisejl` are not
+The dependencies of the Python scripts contained in `pyreisejl` are not
 automatically installed. See `requirements.txt` for details. These requirements
 can be installed using pip:
 ```bash
 pip install -r requirements.txt
 ```
 
 #### Python Installation Verification
-To verify that the included python scripts can successfully run `REISE`, open a python interpreter and run the following commands. They should return with no errors.
+To verify that the included Python scripts can successfully run `REISE`, open a Python interpreter and run the following commands. They should return with no errors.
 ```python
 from julia.api import Julia
 Julia(compiled_modules=False)
@@ -290,9 +290,9 @@ Optional arguments for `REISE.run_scenario` can still be passed as desired.
 
 ## Usage (Python)
 
-The python scripts included in `pyreisejl` perform some additional input validation for the Julia engine before running the simulation and extract data from  the resulting `.mat` files to `.pkl` files.
+The Python scripts included in `pyreisejl` perform some additional input validation for the Julia engine before running the simulation and extract data from  the resulting `.mat` files to `.pkl` files.
 
-There are two main python scripts included in `pyreisejl`:
+There are two main Python scripts included in `pyreisejl`:
 - `pyreisejl/utility/call.py`
 - `pyreisejl/utility/extract_data.py`
 
@@ -317,7 +317,7 @@ pyreisejl/utility/extract_data.py -s '2016-01-01' -e '2016-01-07' -x '/PATH/TO/O
 
 **Note** To see the available options for the `call.py` or `extract_data.py` script, use the `-h, --help` flag when calling the script.
 
-To run the `REISE.jl` simulation from python, using Gurobi as the solver, run `call.py` with the following required options:
+To run the `REISE.jl` simulation from Python, using Gurobi as the solver, run `call.py` with the following required options:
 ```
   -s, --start-date START_DATE
                         The start date for the simulation in format
@@ -341,7 +341,7 @@ Note that the start and end dates need to match dates contained in the input
 profiles (demand, hydro, solar, wind).
 
 
-This python script will validate some of the inputs and translate them into the
+This Python script will validate some of the inputs and translate them into the
 required Julia inputs listed below. By default, the Julia engine creates
 `result_*.mat` files in an `output` folder created in the given input directory.
 To save the matlab files to a different directory, there is an optional flag to
@@ -382,11 +382,13 @@ matlab files generated by the Julia engine:
 * CONGL.pkl (congestion, lower flow limit)
 * AVERAGED_CONG.pkl (time averaged congestion)
 
-If the grid used in the simulation contains DC lines and/or energy storage devices, the following files will also be extracted as necessary:
+If the grid used in the simulation contains DC lines, energy storage devices, or flexible demand resources, the following files will also be extracted as necessary:
 
 * PF_DCLINE.pkl (power flow on DC lines)
 * STORAGE_PG.pkl (power generated by storage units)
 * STORAGE_E.pkl (energy state of charge)
+* LOAD_SHIFT_DN.pkl (demand that is curtailed)
+* LOAD_SHIFT_UP.pkl (demand that is added)
 
 If one or more intervals of the simulation were found to be infeasible without shedding load, the following file will also be extracted:
 * LOAD_SHED.pkl (load shed profile for each load bus)
@@ -451,7 +453,7 @@ automatic extraction flag and to `extract_data.py`:
 
 ### Compatibility with PowerSimData
 
-Within the python code in this repo, there is some code to maintain
+Within the Python code in this repo, there is some code to maintain
 compatibility with the `PowerSimData` framework.
 
 Both `call.py` and `extract_data.py` can be called using a positional
@@ -462,7 +464,7 @@ Manager and Execute List Manager. This option is not intended for manual
 simulation runs.
 
 Note also the different naming convention for various directories by
-`PowerSimData` as compared to the options for the python scripts within
+`PowerSimData` as compared to the options for the Python scripts within
 this repository.
 
 ## Docker
@@ -484,7 +486,7 @@ docker run -it -v /LOCAL/PATH/TO/GUROBI.LIC:/usr/share/gurobi_license -v /LOCAL/
 ```
 
 The following command will start a bash shell session within the container,
-using the `python` commands described above.
+using the Python commands described above.
 
 ```bash
 python pyreisejl/utility/call.py -s '2016-01-01' -e '2016-01-07' -int 24 -i '/usr/share/data'
@@ -556,6 +558,12 @@ at time ![t](https://render.githubusercontent.com/render/math?math=t).
 - ![v_{l,t}](https://render.githubusercontent.com/render/math?math=v_%7Bl%2Ct%7D):
 Branch limit violation for branch ![l](https://render.githubusercontent.com/render/math?math=l)
 at time ![t](https://render.githubusercontent.com/render/math?math=t).
+- ![\delta_{b,t}^{\text{down}}](https://render.githubusercontent.com/render/math?math=%5Cdelta_%7Bb%2Ct%7D%5E%7B%5Ctext%7Bdown%7D%7D):
+Amount of flexible demand curtailed at bus ![b](https://render.githubusercontent.com/render/math?math=b)
+at time ![t](https://render.githubusercontent.com/render/math?math=t).
+- ![\delta_{b,t}^{\text{up}}](https://render.githubusercontent.com/render/math?math=%5Cdelta_%7Bb%2Ct%7D%5E%7B%5Ctext%7Bup%7D%7D):
+Amount of flexible demand added at bus ![b](https://render.githubusercontent.com/render/math?math=b)
+at time ![t](https://render.githubusercontent.com/render/math?math=t).
 - ![{\theta}_{b,t}](https://render.githubusercontent.com/render/math?math=%7B%5Ctheta%7D_%7Bb%2Ct%7D):
 Voltage angle of bus ![b](https://render.githubusercontent.com/render/math?math=b)
 at time ![t](https://render.githubusercontent.com/render/math?math=t).
@@ -619,6 +627,14 @@ Power available from time-varying generator (hydro, wind, solar) ![i](https://re
 at time ![t](https://render.githubusercontent.com/render/math?math=t).
 - ![x_{l}](https://render.githubusercontent.com/render/math?math=x_%7Bl%7D):
 Impedance of branch ![l](https://render.githubusercontent.com/render/math?math=l).
+- ![\underline{\delta}_{b, t}](https://render.githubusercontent.com/render/math?math=%5Cunderline%7B%5Cdelta%7D_%7Bb%2Ct%7D):
+Demand flexibility curtailments available (in MW) at bus ![b](https://render.githubusercontent.com/render/math?math=b)
+at time ![t](https://render.githubusercontent.com/render/math?math=t).
+- ![\overline{\delta}_{b, t}](https://render.githubusercontent.com/render/math?math=%5Coverline%7B%5Cdelta%7D_%7Bb%2Ct%7D):
+Demand flexibility additions available (in MW) at bus ![b](https://render.githubusercontent.com/render/math?math=b)
+at time ![t](https://render.githubusercontent.com/render/math?math=t).
+- ![\Delta^{\text{balance}}](https://render.githubusercontent.com/render/math?math=\Delta^{\text{balance}}):
+The length of the rolling load balance window (in hours), used to account for the duration that flexible demand is deviating from the base demand.
 - ![\eta_{b}^{\text{chg}}](https://render.githubusercontent.com/render/math?math=%5Ceta_%7Bb%7D%5E%7B%5Ctext%7Bchg%7D%7D):
 Charging efficiency of storage device at bus ![b](https://render.githubusercontent.com/render/math?math=b).
 - ![\eta_{b}^{\text{dis}}](https://render.githubusercontent.com/render/math?math=%5Ceta_%7Bb%7D%5E%7B%5Ctext%7Bdis%7D%7D):
@@ -630,12 +646,12 @@ Discharging efficiency of storage device at bus ![b](https://render.githubuserco
 All equations apply over all entries in the indexed sets unless otherwise listed.
 
 - ![0 \le g_{i,s,t} \le g_{i,s,t}^{\text{max}}](https://render.githubusercontent.com/render/math?math=0%5Cle%20g_%7Bi%2Cs%2Ct%7D%5Cle%20g%5E%7B%5Ctext%7Bmax%7D%7D_%7Bi%2Cs%2Ct%7D):
-generator segment power is non-negative and less than the segment width.
-- ![0 \le s_{b,t} \le a^{\text{shed}} \cdot d_{b,t}](https://render.githubusercontent.com/render/math?math=0%20%5Cle%20s_%7Bb%2Ct%7D%20%5Cle%20a%5E%7B%5Ctext%7Bshed%7D%7D%20%5Ccdot%20d_%7Bb%2Ct%7D):
-load shed is non-negative and less than the demand at that bus, if load shedding is enabled.
+Generator segment power is non-negative and less than the segment width.
+- ![0 \le s_{b,t} \le a^{\text{shed}} \cdot ( d_{b,t} + \delta_{b,t}^{\text{up}} - \delta_{b,t}^{\text{down}} )](https://render.githubusercontent.com/render/math?math=0%20%5Cle%20s_%7Bb%2Ct%7D%20%5Cle%20a%5E%7B%5Ctext%7Bshed%7D%7D%20%5Ccdot%20%28%20d_%7Bb%2Ct%7D%20%2B%20%5Cdelta_%7Bb%2Ct%7D%5E%7B%5Ctext%7Bup%7D%7D%20-%20%5Cdelta_%7Bb%2Ct%7D%5E%7B%5Ctext%7Bdown%7D%7D%20%29):
+Load shed is non-negative and less than the demand at that bus (including the impact of demand flexibility), if load shedding is enabled.
 If not, load shed is fixed to 0.
 - ![0 \le v_{b,t} \le a^{\text{viol}} \cdot M](https://render.githubusercontent.com/render/math?math=0%20%5Cle%20v_%7Bb%2Ct%7D%20%5Cle%20a%5E%7B%5Ctext%7Bshed%7D%7D%20%5Ccdot%20M):
-transmission violations are non-negative, if they are enabled
+Transmission violations are non-negative, if they are enabled
 (![M](https://render.githubusercontent.com/render/math?math=M) is a sufficiently large constant that there is no effective upper limit when ![a^{\text{shed}} = 1](https://render.githubusercontent.com/render/math?math=a%5E%7B%5Ctext%7Bshed%7D%7D%20%3D%201)).
 If not, they are fixed to zero.
 - ![0 \le J_{b,t}^{\text{chg}} \le J_{b}^{\text{max}}](https://render.githubusercontent.com/render/math?math=0%20%5Cle%20J_%7Bb%2Ct%7D%5E%7B%5Ctext%7Bchg%7D%7D%20%5Cle%20J_%7Bb%7D%5E%7B%5Ctext%7Bmax%7D%7D):
@@ -648,8 +664,7 @@ Storage state-of-charge is non-negative and limited by the maximum state of char
 Hydro generator power is fixed to the profiles.
 - ![0 \le g_{i,t} \le w_{i,t} \forall i \in I^{\text{S}} \cup I^{\text{W}}](https://render.githubusercontent.com/render/math?math=0%20%5Cle%20g_%7Bi%2Cs%7D%20%5Cle%20w_%7Bi%2Ct%7D%20%5Cforall%20i%20%5Cin%20I%5E%7B%5Ctext%7BS%7D%7D%20%5Ccup%20I%5E%7B%5Ctext%7BW%7D%7D):
 Solar and wind generator power is non-negative and not greater than the availability profiles.
-
-- ![\sum_{i \in I} m_{i,b}^{\text{unit}} g_{i,t} + \sum_{l \in L} m_{l,b}^{\text{line}} f_{l,t} + J_{b,t}^{\text{dis}} + s_{b,t} = d_{b,t} + J_{b,t}^{\text{chg}}](https://render.githubusercontent.com/render/math?math=%5Csum_%7Bi%20%5Cin%20I%7D%20m_%7Bi%2Cb%7D%5E%7B%5Ctext%7Bunit%7D%7D%20g_%7Bi%2Ct%7D%20%2B%20%5Csum_%7Bl%20%5Cin%20L%7D%20m_%7Bl%2Cb%7D%5E%7B%5Ctext%7Bline%7D%7D%20f_%7Bl%2Ct%7D%20%2B%20J_%7Bb%2Ct%7D%5E%7B%5Ctext%7Bdis%7D%7D%20%2B%20s_%7Bb%2Ct%7D%20%3D%20d_%7Bb%2Ct%7D%20%2B%20J_%7Bb%2Ct%7D%5E%7B%5Ctext%7Bchg%7D%7D):
+- ![\sum_{i \in I} m_{i,b}^{\text{unit}} g_{i,t} + \sum_{l \in L} m_{l,b}^{\text{line}} f_{l,t} + J_{b,t}^{\text{dis}} + s_{b,t} + \delta_{b,t}^{\text{down}} = d_{b,t} + J_{b,t}^{\text{chg}} + \delta_{b,t}^{\text{up}}](https://render.githubusercontent.com/render/math?math=%5Csum_%7Bi%20%5Cin%20I%7D%20m_%7Bi%2Cb%7D%5E%7B%5Ctext%7Bunit%7D%7D%20g_%7Bi%2Ct%7D%20%2B%20%5Csum_%7Bl%20%5Cin%20L%7D%20m_%7Bl%2Cb%7D%5E%7B%5Ctext%7Bline%7D%7D%20f_%7Bl%2Ct%7D%20%2B%20J_%7Bb%2Ct%7D%5E%7B%5Ctext%7Bdis%7D%7D%20%2B%20s_%7Bb%2Ct%7D%20%2B%20%5Cdelta_%7Bb%2Ct%7D%5E%7B%5Ctext%7Bdown%7D%7D%20%3D%20d_%7Bb%2Ct%7D%20%2B%20J_%7Bb%2Ct%7D%5E%7B%5Ctext%7Bchg%7D%7D%20%2B%20%5Cdelta_%7Bb%2Ct%7D%5E%7B%5Ctext%7Bup%7D%7D%20):
 Power balance at each bus ![b](https://render.githubusercontent.com/render/math?math=b)
 at time ![t](https://render.githubusercontent.com/render/math?math=t).
 - ![g_{i,t} = g_{i}^{\text{min}} + \sum_{s \in S} g_{i,s,t}](https://render.githubusercontent.com/render/math?math=g_%7Bi%2Ct%7D%20%3D%20g_%7Bi%7D%5E%7B%5Ctext%7Bmin%7D%7D%20%2B%20%5Csum_%7Bs%20%5Cin%20S%7D%20g_%7Bi%2Cs%2Ct%7D):
@@ -664,6 +679,14 @@ Ramp-down constraint.
 Power flow over each branch is limited by the branch power limit, and can only exceed this value by using the 'violation' variable (if enabled), which is penalized in the objective function.
 - ![f_{l,t} = \frac{1}{x_{l}} \sum_{b \in B} m_{l,b}^{\text{line}} \theta_{b,t}](https://render.githubusercontent.com/render/math?math=f_%7Bl%2Ct%7D%20%3D%20%5Cfrac%7B1%7D%7Bx_%7Bl%7D%7D%20%5Csum_%7Bb%20%5Cin%20B%7D%20m_%7Bl%2Cb%7D%5E%7B%5Ctext%7Bline%7D%7D%20%5Ctheta_%7Bb%2Ct%7D):
 Power flow over each branch is proportional to the admittance and the angle difference.
+- ![0 \le \delta_{b,t}^{\text{down}} \le \underline{\delta_{b,t}}](https://render.githubusercontent.com/render/math?math=0%20%5Cle%20%5Cdelta_%7Bb%2Ct%7D%5E%7B%5Ctext%7Bdown%7D%7D%20%5Cle%20%5Cunderline%7B%5Cdelta%7D_%7Bb%2Ct%7D):
+Bound on the amount of demand that flexible demand resources can curtail.
+- ![0 \le \delta_{b,t}^{\text{up}} \le \overline{\delta_{b,t}}](https://render.githubusercontent.com/render/math?math=0%20%5Cle%20%5Cdelta_%7Bb%2Ct%7D%5E%7B%5Ctext%7Bup%7D%7D%20%5Cle%20%5Coverline%7B%5Cdelta%7D_%7Bb%2Ct%7D):
+Bound on the amount of demand that flexible demand resources can add.
+- ![\sum_{t = k}^{k + \Delta^{\text{balance}}} \delta_{b,t}^{\text{up}} - \delta_{b,t}^{\text{down}} \ge 0, \quad \forall b \in B, \quad k = 1, ..., |T| - \Delta^{\text{balance}}](https://render.githubusercontent.com/render/math?math=\sum_{t%20=%20k}^{k%20%2B%20\Delta^{\text{balance}}}%20\delta_{b,t}^{\text{up}}%20-%20\delta_{b,t}^{\text{down}}%20\ge%200,%20\quad%20\forall%20b%20\in%20B,%20\quad%20k%20=%201,%20...,%20|T|%20-%20\Delta^{\text{balance}}):
+Rolling load balance for flexible demand resources; used to restrict the time that flexible demand resources can deviate from the base demand.
+- ![\sum_{t \in T} \delta_{b,t}^{\text{up}} - \delta_{b,t}^{\text{down}} \ge 0, \quad \forall b \in B](https://render.githubusercontent.com/render/math?math=\sum_{t%20\in%20T}%20\delta_{b,t}^{\text{up}}%20-%20\delta_{b,t}^{\text{down}}%20\ge%200,%20\quad%20\forall%20b%20\in%20B):
+Interval load balance for flexible demand resources.
 
 ### Objective function
 

pyreisejl/utility/extract_data.py

@@ -113,6 +113,12 @@ def extract_data(results):
             extraction_vars |= {"load_shed"}
         except KeyError:
             pass
+        try:
+            temps["load_shift_up"] = output_mpc["flexible_demand"]["load_shift_up"].T
+            temps["load_shift_dn"] = output_mpc["flexible_demand"]["load_shift_dn"].T
+            extraction_vars |= {"load_shift_up", "load_shift_dn"}
+        except KeyError:
+            pass
 
         # Extract which number result currently being processed
         i = result_num(filename)
@@ -233,6 +239,8 @@ def _get_outputs_from_converted(matfile):
         "pf": case.mpc.branchid,
         "lmp": case.mpc.bus[:, 0].astype(np.int64),
         "load_shed": case.mpc.bus[:, 0].astype(np.int64),
+        "load_shift_up": case.mpc.bus[:, 0].astype(np.int64),
+        "load_shift_dn": case.mpc.bus[:, 0].astype(np.int64),
         "congu": case.mpc.branchid,
         "congl": case.mpc.branchid,
     }

src/REISE.jl

@@ -9,9 +9,9 @@ import Requires
 import SparseArrays: sparse, SparseMatrixCSC
 
 
-include("types.jl")         # Defines Case, Results, Storage,
+include("types.jl")         # Defines Case, Results, Storage, DemandFlexibility,
                             #     VariablesOfInterest
-include("read.jl")          # Defines read_case, read_storage
+include("read.jl")          # Defines read_case, read_storage, read_demand_flexibility
 include("prepare.jl")       # Defines reise_data_mods
 include("model.jl")         # Defines _build_model (used in interval_loop)
 include("loop.jl")          # Defines interval_loop
@@ -61,12 +61,14 @@ function run_scenario(;
     stderr_filepath = joinpath(outputfolder, "stderr.err")
     case = read_case(inputfolder)
     storage = read_storage(inputfolder)
+    demand_flexibility = read_demand_flexibility(inputfolder)
     println("All scenario files loaded!")
     case = reise_data_mods(case, num_segments=num_segments)
     save_input_mat(case, storage, inputfolder, outputfolder)
     model_kwargs = Dict(
         "case" => case,
         "storage" => storage,
+        "demand_flexibility" => demand_flexibility,
         "interval_length" => interval,
         )
     # If a number of threads is specified, add to solver settings dict

src/loop.jl

@@ -40,6 +40,7 @@ function interval_loop(factory_like, model_kwargs::Dict,
     # Constant parameters
     case = model_kwargs["case"]
     storage = model_kwargs["storage"]
+    demand_flexibility = model_kwargs["demand_flexibility"]
     sets = _make_sets(case, storage)
     storage_enabled = (sets.num_storage > 0)
     unused_load_shed_intervals_turnoff = 14
@@ -52,6 +53,11 @@ function interval_loop(factory_like, model_kwargs::Dict,
         interval_start = start_index + (i - 1) * interval
         interval_end = interval_start + interval - 1
         model_kwargs["start_index"] = interval_start
+        if demand_flexibility.enabled
+            bus_flex_amt = _make_bus_demand_flexibility_amount(
+                case, demand_flexibility, interval_start, interval_end
+            )
+        end
         if i == 1
             # Build a model with no initial ramp constraint
             if storage_enabled
@@ -86,8 +92,9 @@ function interval_loop(factory_like, model_kwargs::Dict,
             if (("load_shed_enabled" in keys(model_kwargs))
                 && (model_kwargs["load_shed_enabled"] == true))
                 for t in 1:interval, i in 1:length(sets.load_bus_idx)
-                    JuMP.set_upper_bound(voi.load_shed[i, t],
-                                         bus_demand[sets.load_bus_idx[i], t])
+                    JuMP.set_normalized_rhs(
+                        voi.load_shed_ub[i, t], bus_demand[sets.load_bus_idx[i], t]
+                    )
                 end
             end
             for t in 1:interval, g in 1:sets.num_hydro
@@ -115,6 +122,16 @@ function interval_loop(factory_like, model_kwargs::Dict,
                     JuMP.set_normalized_rhs(voi.initial_soc[s], storage_e0[s])
                 end
             end
+            if demand_flexibility.enabled
+                for t in 1:interval, i in 1:length(sets.load_bus_idx)
+                    JuMP.set_upper_bound(
+                        voi.load_shift_up[i, t], bus_flex_amt[sets.load_bus_idx[i], t]
+                    )
+                    JuMP.set_upper_bound(
+                        voi.load_shift_dn[i, t], bus_flex_amt[sets.load_bus_idx[i], t]
+                    )
+                end
+            end
         end
 
         while true