Add example 5bus case, contingency screening capability, and (#128)

* initial pass at 5bus example * allowing default forecastables * relaxing line limits * adding runner script * correcting issue with REAL_TIME period resolution * Revert "relaxing line limits" This reverts commit b7d9eb3. * updating network and adding generator * ensure reserve factor is honored even if not pre-existing * first pass at enabling contingency monitoring * considering more flexible thermal units * adding contingency reporting * updating price thresholds * switching penalty logic to use grid-parity-exchange/Egret#259 * NFC: removing finished TODO * fixing names * adding python simulator script * more reasonable mipgap * removing violation, for now * Fix spelling, update variable name Co-authored-by: Darryl Melander <[email protected]>
grid-parity-exchange · Jan 12, 2022 · 8fef0fa · 8fef0fa
1 parent 05492ba
commit 8fef0fa
Show file tree

Hide file tree

Showing 22 changed files with 228,714 additions and 48 deletions.
diff --git a/examples/5bus/DAY_AHEAD_load.csv b/examples/5bus/DAY_AHEAD_load.csv
diff --git a/examples/5bus/DAY_AHEAD_renewables.csv b/examples/5bus/DAY_AHEAD_renewables.csv
diff --git a/examples/5bus/REAL_TIME_load.csv b/examples/5bus/REAL_TIME_load.csv
diff --git a/examples/5bus/REAL_TIME_renewables.csv b/examples/5bus/REAL_TIME_renewables.csv
diff --git a/examples/5bus/branch.csv b/examples/5bus/branch.csv
@@ -0,0 +1,8 @@
+UID,From Bus,To Bus,R,X,B,Cont Rating,LTE Rating,STE Rating,Tr Ratio
+"branch_2_3",2,3,0.00108,0.0108,0.01852,80.0,92.0,104.0,0
+"branch_1_2",1,2,0.00281,0.0281,0.00712,66.6,76.6,86.6,0
+"branch_1_4",1,4,0.00304,0.0304,0.00658,66.6,76.6,86.6,0
+"branch_4_10",4,10,0.00297,0.0297,0.00674,66.6,76.6,86.6,0
+"branch_1_10",1,10,0.00064,0.0064,0.03126,66.6,76.6,86.6,0
+"branch_3_4_0",3,4,0.00297,0.0297,0.00337,66.6,76.6,86.6,1.05
+"branch_3_4_1",3,4,0.003274425,0.03274425,0.003056689,28.4,32.7,36.9,0.952380952
diff --git a/examples/5bus/bus.csv b/examples/5bus/bus.csv
@@ -0,0 +1,6 @@
+Bus ID,Bus Name,BaseKV,Bus Type,MW Load,MVAR Load,V Mag,V Angle,Area,Zone
+1,bus1,230,PV,0,0,1,0.048935018,1,1
+4,bus4,230,Ref,28.5714286,9.3907143,1,0,1,1
+10,bus5,230,PV,0,0,1,0.06266308,1,1
+2,bus2,230,PQ,21.4285714,7.0435714,1.04407,-0.012822061,2,1
+3,bus3,230,PV,21.4285714,7.0435714,1,-0.009768957,2,1
diff --git a/examples/5bus/gen.csv b/examples/5bus/gen.csv
@@ -0,0 +1,9 @@
+GEN UID,Bus ID,Unit Type,Fuel,MW Inj,MVAR Inj,V Setpoint p.u.,PMax MW,PMin MW,QMax MVAR,QMin MVAR,Min Down Time Hr,Min Up Time Hr,Ramp Rate MW/Min,Start Time Cold Hr,Start Time Warm Hr,Start Time Hot Hr,Start Heat Cold MBTU,Start Heat Warm MBTU,Start Heat Hot MBTU,Non Fuel Start Cost $,Fuel Price $/MMBTU,Output_pct_0,Output_pct_1,Output_pct_2,Output_pct_3,Output_pct_4,HR_avg_0,HR_incr_1,HR_incr_2,HR_incr_3,HR_incr_4
+3_CT,3,CT,G,0,0,1,20,8,20,-20,1,1,3,,,1,,,51.75,0,0.75,0.4,0.6,0.8,1,,135722.5,97862.5,98072.5,107135,
+10_STEAM,10,STEAM,C,0,0,1,76,30,76,-76,4,8,2,12,10,4,11383.41,10488.35,7355.42,0,1,0.394736842,0.596052632,0.798684211,1,,30166.66667,14402.61434,17182.46753,18283.66017,
+4_CC,4,CC,G,0,0,1,100,10,55,-55,3,3,3.7,,,3,,,5665.23,0,0.5,0.4,0.6,0.8,1,,51019.54545,26818.18182,29550.90909,30308.18182,
+4_STEAM,4,STEAM,O,0,0,1,12,5,12,-12,2,4,1,12,4,2,703.76,455.37,393.28,0,1.5,0.416666667,0.608333333,0.808333333,1,,179457.9999,124504.3483,125050,133643.478,
+10_PV,10,PV,S,0,0,1,25.9,0,0,0,,,,,,,,,,0,,,,,,,,,,,
+2_RTPV,2,RTPV,S,0,0,1,9.3,0,0,0,,,,,,,,,,0,,,,,,,,,,,
+1_HYDRO,1,HYDRO,H,0,0,1,50,0,16,-10,,,,,,,,,,0,,,,,,,,,,,
+4_WIND,4,WIND,W,0,0,1,120,0,0,0,,,,,,,,,,0,,,,,,,,,,,
diff --git a/examples/5bus/initial_status.csv b/examples/5bus/initial_status.csv
@@ -0,0 +1,3 @@
+3_CT,10_STEAM,4_CC,4_STEAM
+24,-24,24,-24,
+14,0,40,0
diff --git a/examples/5bus/reserves.csv b/examples/5bus/reserves.csv
@@ -0,0 +1 @@
+Reserve Product,Requirement (MW)
diff --git a/examples/5bus/simulation_objects.csv b/examples/5bus/simulation_objects.csv
@@ -0,0 +1,7 @@
+Simulation_Parameters,Description,DAY_AHEAD,REAL_TIME
+Periods_per_Step,the number of discrete periods represented in each simulation step,24,1
+Period_Resolution,seconds per period,3600,300
+Date_From,simulation beginning period,01/01/2020 0:00,12/31/2020 0:00
+Date_To,simulation ending period (must account for lookahed data availability),01/01/2020 0:00,12/31/2020 0:00
+Look_Ahead_Periods_per_Step,the number of look ahead periods included in each optimization step,24,2
+Look_Ahead_Resolution,look-ahead period resolution,3600,300
diff --git a/examples/5bus/timeseries_pointers.csv b/examples/5bus/timeseries_pointers.csv
@@ -0,0 +1,17 @@
+Simulation,Category,Object,Parameter,Data File
+DAY_AHEAD,Generator,1_HYDRO,PMin MW,DAY_AHEAD_renewables.csv
+DAY_AHEAD,Generator,2_RTPV,PMin MW,DAY_AHEAD_renewables.csv
+DAY_AHEAD,Generator,1_HYDRO,PMax MW,DAY_AHEAD_renewables.csv
+DAY_AHEAD,Generator,2_RTPV,PMax MW,DAY_AHEAD_renewables.csv
+DAY_AHEAD,Generator,10_PV,PMax MW,DAY_AHEAD_renewables.csv
+DAY_AHEAD,Generator,4_WIND,PMax MW,DAY_AHEAD_renewables.csv
+DAY_AHEAD,Area,1,MW Load,DAY_AHEAD_load.csv
+DAY_AHEAD,Area,2,MW Load,DAY_AHEAD_load.csv
+REAL_TIME,Generator,1_HYDRO,PMin MW,REAL_TIME_renewables.csv
+REAL_TIME,Generator,2_RTPV,PMin MW,REAL_TIME_renewables.csv
+REAL_TIME,Generator,1_HYDRO,PMax MW,REAL_TIME_renewables.csv
+REAL_TIME,Generator,2_RTPV,PMax MW,REAL_TIME_renewables.csv
+REAL_TIME,Generator,10_PV,PMax MW,REAL_TIME_renewables.csv
+REAL_TIME,Generator,4_WIND,PMax MW,REAL_TIME_renewables.csv
+REAL_TIME,Area,1,MW Load,REAL_TIME_load.csv
+REAL_TIME,Area,2,MW Load,REAL_TIME_load.csv
diff --git a/examples/simulate_5bus.py b/examples/simulate_5bus.py
@@ -0,0 +1,28 @@
+from prescient.simulator import Prescient
+
+# set some options
+prescient_options = {
+        "data_path":"./5bus/",
+        "input_format":"rts-gmlc",
+        "simulate_out_of_sample":True,
+        "run_sced_with_persistent_forecast_errors":True,
+        "output_directory":"5bus_output",
+        "start_date":"07-10-2020",
+        "num_days":7,
+        "sced_horizon":1,
+        "ruc_mipgap":0.01,
+        "reserve_factor":0.1,
+        "deterministic_ruc_solver":"cbc",
+        "deterministic_ruc_solver_options":{"feas":"off", "DivingF":"on",},
+        "sced_solver":"cbc",
+        "sced_frequency_minutes":5,
+        "ruc_horizon":36,
+        "compute_market_settlements":True,
+        "monitor_all_contingencies":True,
+        "output_solver_logs":False,
+        "price_threshold":1000,
+        "contingency_price_threshold":100,
+        "reserve_price_threshold":5,
+        }
+# run the simulator
+Prescient().simulate(**prescient_options)
diff --git a/examples/simulate_5bus.txt b/examples/simulate_5bus.txt
@@ -0,0 +1,22 @@
+command/exec simulator.py
+--data-directory=5bus
+--input-format=rts-gmlc
+--simulate-out-of-sample 
+--run-sced-with-persistent-forecast-errors 
+--output-directory=5bus_output
+--start-date=07-10-2020
+--num-days=7
+--sced-horizon=1
+--ruc-mipgap=0.01
+--reserve-factor=0.1
+--deterministic-ruc-solver=cbc
+--deterministic-ruc-solver-options="feas=off DivingG=on"
+--sced-solver=cbc
+--sced-frequency-minutes=5
+--ruc-horizon=36
+--compute-market-settlements
+--monitor-all-contingencies
+#--output-solver-logs
+--price-threshold=1000
+--contingency-price-threshold=100
+--reserve-price-threshold=5
diff --git a/prescient/engine/data_extractors.py b/prescient/engine/data_extractors.py
@@ -173,6 +173,20 @@ def get_max_power_available(self, sced: OperationsModel, g: G) -> float:
     def get_flow_level(self, sced: OperationsModel, line: L) -> float:
         pass
 
+    @abstractmethod
+    def get_flow_violation_level(self, sced: OperationsModel, line: L) -> float:
+        pass
+
+    @abstractmethod
+    def get_all_contingency_flow_levels(self, sced: OperationsModel) -> Dict[Tuple[L,L], float]:
+        """Get the flows for the monitored contingencies."""
+        pass
+
+    @abstractmethod
+    def get_all_contingency_flow_violation_levels(self, sced: OperationsModel) -> Dict[Tuple[L,L], float]:
+        """Get the flow violations for the monitored contingencies."""
+        pass
+
     @abstractmethod
     def get_storage_input_dispatch_level(self, sced: OperationsModel, storage: S) -> float:
         pass
@@ -425,6 +439,10 @@ def get_all_flow_levels(self, sced: OperationsModel) -> Dict[L, float]:
         return {l: self.get_flow_level(sced, l)
                 for l in self.get_transmission_lines(sced)}
 
+    def get_all_flow_violation_levels(self, sced: OperationsModel) -> Dict[L, float]:
+        return {l: self.get_flow_violation_level(sced, l)
+                for l in self.get_transmission_lines(sced)}
+
     def get_all_bus_demands(self, sced: OperationsModel) -> Dict[B, float]:
         return {b: self.get_bus_demand(sced, b)
                 for b in self.get_buses(sced)}

diff --git a/prescient/engine/egret/data_extractors.py b/prescient/engine/egret/data_extractors.py
@@ -10,7 +10,7 @@
 from __future__ import annotations
 from typing import TYPE_CHECKING
 if TYPE_CHECKING:
-    from typing import Iterable
+    from typing import Iterable, Dict, Tuple
     from prescient.engine.abstract_types import *
 
 import numpy as np
@@ -139,7 +139,7 @@ def get_thermal_headroom(self, sced: OperationsModel, g: G) -> float:
         return val
 
     def get_min_downtime(self, sced: OperationsModel, g: G) -> float:
-        return sced.data['elements']['generator'][g]['min_up_time']
+        return sced.data['elements']['generator'][g]['min_down_time']
 
     def get_scaled_startup_ramp_limit(self, sced: OperationsModel, g: G) -> float:
         return sced.data['elements']['generator'][g]['startup_capacity']
@@ -187,6 +187,30 @@ def get_generator_cost(self, sced: OperationsModel, g: G) -> float:
     def get_flow_level(self, sced: OperationsModel, line: L) -> float:
         return sced.data['elements']['branch'][line]['pf']['values'][0]
 
+    def get_flow_violation_level(self, sced: OperationsModel, line: L) -> float:
+        pf_violation = sced.data['elements']['branch'][line].get('pf_violation', 0.)
+        if pf_violation != 0.:
+            pf_violation = pf_violation['values'][0]
+        return pf_violation
+
+    def get_all_contingency_flow_levels(self, sced: OperationsModel) -> Dict[Tuple[L,L], float]:
+        contingency_dict = {}
+        for c_dict in sced.data['elements'].get('contingency', {}).values():
+            line_out = c_dict['branch_contingency']
+            monitored_branches = c_dict.get('monitored_branches',{'values':[{}]})
+            for bn, b_dict in monitored_branches['values'][0].items():
+                contingency_dict[line_out, bn] = b_dict['pf']
+        return contingency_dict
+
+    def get_all_contingency_flow_violation_levels(self, sced: OperationsModel) -> Dict[Tuple[L,L], float]:
+        contingency_viol = {}
+        for c_dict in sced.data['elements'].get('contingency', {}).values():
+            line_out = c_dict['branch_contingency']
+            monitored_branches = c_dict.get('monitored_branches',{'values':[{}]})
+            for bn, b_dict in monitored_branches['values'][0].items():
+                contingency_viol[line_out, bn] = b_dict.get('pf_violation', 0.)
+        return contingency_viol
+
     def get_bus_mismatch(self, sced: OperationsModel, bus: B) -> float:
         return self.get_load_mismatch(sced, bus)
 

diff --git a/prescient/engine/egret/egret_plugin.py b/prescient/engine/egret/egret_plugin.py
@@ -21,8 +21,9 @@
 from egret.data.model_data import ModelData
 from egret.parsers.prescient_dat_parser import get_uc_model, create_model_data_dict_params
 from egret.models.unit_commitment import _time_series_dict, _preallocated_list, _solve_unit_commitment, \
-                                         _save_uc_results, create_tight_unit_commitment_model, \
-                                         _get_uc_model
+                                        _save_uc_results, create_tight_unit_commitment_model, \
+                                        _get_uc_model
+from egret.model_library.transmission.tx_calc import construct_connection_graph, get_N_minus_1_branches
 
 from prescient.util import DEFAULT_MAX_LABEL_LENGTH
 from prescient.util.math_utils import round_small_values
@@ -146,6 +147,7 @@ def create_sced_instance(data_provider:DataProvider,
                 sced_data[t] = forecast[t] * forecast_error_ratio
 
     _ensure_reserve_factor_honored(options, sced_md, range(sced_horizon))
+    _ensure_contingencies_monitored(options, sced_md)
 
     # Set generator commitments & future state
     for g, g_dict in sced_md.elements(element_type='generator', generator_type='thermal'):
@@ -368,8 +370,10 @@ def create_deterministic_ruc(options,
     # Create a new model
     md = data_provider.get_initial_forecast_model(options, ruc_horizon, 60)
 
+    initial_ruc = current_state is None or current_state.timestep_count == 0
+
     # Populate the T0 data
-    if current_state is None or current_state.timestep_count == 0:
+    if initial_ruc:
         data_provider.populate_initial_state_data(options, md)
     else:
         _copy_initial_state_into_model(options, current_state, md)
@@ -404,6 +408,8 @@ def create_deterministic_ruc(options,
     # Ensure the reserve requirement is satisfied
     _ensure_reserve_factor_honored(options, md, range(ruc_horizon))
 
+    _ensure_contingencies_monitored(options, md, initial_ruc)
+
     return md
 
 
@@ -483,14 +489,12 @@ def solve_deterministic_day_ahead_pricing_problem(solver, ruc_results, options,
     reserve_requirement = ('reserve_requirement' in pricing_instance.data['system'])
 
     system = pricing_instance.data['system']
-    # In case of demand shortfall, the price skyrockets, so we threshold the value.
-    if ('load_mismatch_cost' not in system) or (system['load_mismatch_cost'] > options.price_threshold):
-        system['load_mismatch_cost'] = options.price_threshold
 
-    # In case of reserve shortfall, the price skyrockets, so we threshold the value.
-    if reserve_requirement:
-        if ('reserve_shortfall_cost' not in system) or (system['reserve_shortfall_cost'] > options.reserve_price_threshold):
-            system['reserve_shortfall_cost'] = options.reserve_price_threshold
+    # In case of shortfall, the price skyrockets, so we threshold the value.
+    for system_key, threshold_value in get_attrs_to_price_option(options):
+        if threshold_value is not None and ((system_key not in system) or
+                (system[system_key] > threshold_value)):
+            system[system_key] = threshold_value
 
     ptdf_manager.mark_active(pricing_instance)
     pyo_model = create_pricing_model(pricing_instance, relaxed=True,
@@ -675,6 +679,38 @@ def _ensure_reserve_factor_honored(options:Options, md:EgretModel, time_periods:
             if reserve_reqs[t] < min_reserve:
                 reserve_reqs[t] = min_reserve
 
+def _ensure_contingencies_monitored(options:Options, md:EgretModel, initial_ruc:bool = False) -> None:
+    ''' Add contingency screening, if that option is enabled '''
+    if initial_ruc:
+        _ensure_contingencies_monitored.contingency_dicts = {}
+
+    for bn, b in md.elements('branch'): 
+        if not b.get('in_service', True):
+            raise RuntimeError(f"Remove branches from service by setting the `planned_outage` attribute. "
+                    f"Branch {bn} has `in_service`:False")
+    for bn, b in md.elements('bus'): 
+        if not b.get('in_service', True):
+            raise RuntimeError(f"Buses cannot be removed from service in Prescient")
+
+    if options.monitor_all_contingencies:
+        key = []
+        for bn, b in md.elements('branch'):
+            if 'planned_outage' in b:
+                if isinstance(b['planned_outage'], dict):
+                    if any(b['planned_outage']['values']):
+                        key.append(b)
+                elif b['planned_outage']:
+                    key.append(b)
+        key = tuple(key)
+        if key not in _ensure_contingencies_monitored.contingency_dicts:
+            mapping_bus_to_idx = { k : i for i,k in enumerate(md.data['elements']['bus'].keys())}
+            graph = construct_connection_graph(md.data['elements']['branch'], mapping_bus_to_idx)
+            contingency_list = get_N_minus_1_branches(graph, md.data['elements']['branch'], mapping_bus_to_idx)
+            contingency_dict = { cn : {'branch_contingency':cn} for cn in contingency_list} 
+            _ensure_contingencies_monitored.contingency_dicts[key] = contingency_dict
+
+        md.data['elements']['contingency'] = _ensure_contingencies_monitored.contingency_dicts[key]
+
 def _copy_initial_state_into_model(options:Options, 
                                    current_state:SimulationState, 
                                    md:EgretModel):
@@ -683,3 +719,21 @@ def _copy_initial_state_into_model(options:Options,
         g_dict['initial_p_output']  = current_state.get_initial_power_generated(g)
     for s,s_dict in md.elements('storage'):
         s_dict['initial_state_of_charge'] = current_state.get_initial_state_of_charge(s)
+
+def get_attrs_to_price_option(options:Options):
+    '''
+    Create a map from internal attributes to various price thresholds
+    for the LMP SCED
+    '''
+    return {
+            'load_mismatch_cost' : options.price_threshold,
+            'contingency_flow_violation_cost' : options.contingency_price_threshold,
+            'transmission_flow_violation_cost' : options.transmission_price_threshold,
+            'interface_flow_violation_cost' : options.interface_price_threshold,
+            'reserve_shortfall_cost' : options.reserve_price_threshold,
+            'regulation_penalty_price' : options.regulation_price_threshold,
+            'spinning_reserve_penalty_price' : options.spinning_reserve_price_threshold,
+            'non_spinning_reserve_penalty_price' : options.non_spinning_reserve_price_threshold,
+            'supplemental_reserve_penalty_price' : options.supplemental_reserve_price_threshold,
+            'flexible_ramp_penalty_price' : options.flex_ramp_price_threshold,
+            }.items()
diff --git a/prescient/engine/egret/engine.py b/prescient/engine/egret/engine.py
@@ -33,7 +33,7 @@
 from egret.common.lazy_ptdf_utils import uc_instance_binary_relaxer
 
 from prescient.engine.modeling_engine import SlackType as EngineSlackType, NetworkType as EngineNetworkType
-from egret.model_library.unit_commitment.uc_utils import SlackType as EgretSlackType
+from egret.model_library.unit_commitment.uc_utils import SlackType as EgretSlackType, reset_unit_commitment_penalties
 
 
 _network_type_to_egret_network_constraints = {
@@ -259,17 +259,11 @@ def create_and_solve_lmp(self,
         lmp_sced_instance = sced_instance.clone()
 
         # In case of demand shortfall, the price skyrockets, so we threshold the value.
-        if 'load_mismatch_cost' not in lmp_sced_instance.data['system'] or \
-                lmp_sced_instance.data['system']['load_mismatch_cost'] > \
-                    options.price_threshold:
-            lmp_sced_instance.data['system']['load_mismatch_cost'] = options.price_threshold
-
-        # In case of reserve shortfall, the price skyrockets, so we threshold the value.
-        if 'reserve_shortfall_cost' not in lmp_sced_instance.data['system'] or \
-                lmp_sced_instance.data['system']['reserve_shortfall_cost'] > \
-                    options.reserve_price_threshold:
-            lmp_sced_instance.data['system']['reserve_shortfall_cost'] = \
-                    options.reserve_price_threshold
+        system = lmp_sced_instance.data['system']
+        for system_key, threshold_value in self._p.get_attrs_to_price_option(options):
+            if threshold_value is not None and ((system_key not in system) or
+                    (system[system_key] > threshold_value)):
+                system[system_key] = threshold_value
 
         if self._last_sced_pyo_model is None:
             self._ptdf_manager.mark_active(lmp_sced_instance)
@@ -314,23 +308,10 @@ def _transform_for_lmp(self, pyo_model, pyo_solver, lmp_sced_instance):
         uc_instance_binary_relaxer(pyo_model, pyo_solver)
 
         ## reset the penalites
-        system = lmp_sced_instance.data['system']
-
-        update_obj = False
-
-        new_load_penalty = system['baseMVA'] * system['load_mismatch_cost']
-        if not math.isclose(new_load_penalty, pyo_model.LoadMismatchPenalty.value):
-            pyo_model.LoadMismatchPenalty.value = new_load_penalty
-            update_obj = True
-
-        new_reserve_penalty =  system['baseMVA'] * system['reserve_shortfall_cost']
-        if not math.isclose(new_reserve_penalty, pyo_model.ReserveShortfallPenalty.value):
-            pyo_model.ReserveShortfallPenalty.value = new_reserve_penalty
-            update_obj = True
-
         pyo_model.model_data = lmp_sced_instance
+        reset_unit_commitment_penalties(pyo_model)
 
-        if update_obj and isinstance(pyo_solver, PersistentSolver):
+        if isinstance(pyo_solver, PersistentSolver):
             pyo_solver.set_objective(pyo_model.TotalCostObjective)
 
     def _print_sced_info(self,
@@ -456,8 +437,8 @@ def _get_solver_list(name):
         if not pe.SolverFactory(self._sced_solver).available():
             raise RuntimeError(f"Solver {self._sced_solver} is not available to Pyomo")
 
-
-    def _print_persistence_warning(self, solver):
+    @staticmethod
+    def _print_persistence_warning(solver):
         print(f"WARNING: Solver {solver} supports persistence, which "
                "improves the performance of Prescient. Consider installing the "
               f"python bindings for {solver}.")
@@ -487,6 +468,7 @@ def __init__(self, options: Options):
             self.create_simulation_actuals = egret_plugin.create_simulation_actuals
             self.solve_deterministic_day_ahead_pricing_problem = egret_plugin.solve_deterministic_day_ahead_pricing_problem
             self._zero_out_costs = egret_plugin._zero_out_costs
+            self.get_attrs_to_price_option = egret_plugin.get_attrs_to_price_option
 
             if options.simulator_plugin != None:
                 try: