one penalty to rule them all

egret/model_library/transmission/tx_utils.py

@@ -239,6 +239,11 @@ def load_shed_limit(load, gens, gen_mins):
                             'flexible_ramp_up_price',
                             'flexible_ramp_down_price',
                             'supplemental_price',
+                            'regulation_penalty_price',
+                            'spinning_reserve_penalty_price',
+                            'non_spinning_reserve_penalty_price',
+                            'supplemental_reserve_penalty_price',
+                            'flexible_ramp_penalty_price',
                         ]
 
 ## TODO?: break apart by data that needed to be scaled down (capacity limits, power),
@@ -388,6 +393,9 @@ def load_shed_limit(load, gens, gen_mins):
                        ('system_attributes', None, None ) : [
                                                         'load_mismatch_cost',
                                                         'q_load_mismatch_cost',
+                                                        'transmission_flow_violation_cost',
+                                                        'contingency_flow_violation_cost',
+                                                        'interface_flow_violation_cost',
                                                         'reserve_shortfall_cost',
                                                      ] + \
                                                      ancillary_service_stack,

egret/model_library/unit_commitment/params.py

@@ -14,7 +14,7 @@
 from egret.model_library.transmission import tx_utils
 from egret.common.log import logger
 
-from .uc_utils import add_model_attr, uc_time_helper, SlackType
+from .uc_utils import add_model_attr, uc_time_helper, SlackType, make_penalty_rule, make_indexed_penalty_rule
 
 component_name = 'data_loader'
 
@@ -190,32 +190,52 @@ def load_params(model, model_data, slack_type):
     model.GenerationTimeInStage = Set(model.StageSet, within=model.TimePeriods,
                                         initialize={'Stage_1': list(), 'Stage_2': model.TimePeriods } )
 
-    ##########################################
-    # penalty costs for constraint violation #
-    ##########################################
-
-    ModeratelyBigPenalty = 1e3*system['baseMVA']
+    #################################################################################
+    # penalty costs for constraint violation
+    #
+    # While the user can specify these, by default we base all penalties
+    # off the "load_mismatch_cost", which always has the highest penalty
+    # value (default $1M/MWh). If the user sets "load_mismatch_cost"
+    # at $1000/MWh, the following penalties will be used:
+    #
+    # (defined here in params.py)
+    # "q_load_mismatch_cost"              : $500/MVh ("load_mismatch_cost"/2)
+    # "transmission_flow_violation_cost"  : $500/MWh ("load_mismatch_cost"/2)
+    # "contingency_flow_violation_cost"   : $500/MWh ("load_mismatch_cost"/2)
+    # "interface_flow_violation_cost"     : $300/MWh ("load_mismatch_cost"/(10/3))
+    # "reserve_shortfall_cost"            : $100/MWh ("load_mismatch_cost"/10)
+    #
+    # (defined in services.py)
+    # "regulation_penalty_price"          : $250/MWh ("load_mismatch_cost"/4)
+    # "spinning_reserve_penalty_price"    : $200/MWh ("load_mismatch_cost"/5)
+    # "non_spinning_reserve_penalty_price": $150/MWh ("load_mismatch_cost"/(20/3))
+    # "supplemental_reserve_penalty_price": $125/MWh ("load_mismatch_cost"/8)
+    # "flexible_ramp_penalty_price"       : $110/MWh ("load_mismatch_cost"/(100/11))
+    #
+    # Note these can be overridden by the user specifying the values themselves.
+    # Further, penalties on branch flows and interfaces can be set per-element.
+    ################################################################################
 
-    model.ReserveShortfallPenalty = Param(within=NonNegativeReals, default=ModeratelyBigPenalty, mutable=True, initialize=system.get('reserve_shortfall_cost', ModeratelyBigPenalty))
+    BigPenalty = 1e6*system['baseMVA']
 
-    BigPenalty = 1e4*system['baseMVA']
+    model.LoadMismatchPenalty = Param(within=NonNegativeReals, mutable=True, rule=lambda m : m.model_data.data['system'].get('load_mismatch_cost', BigPenalty))
+    model.LoadMismatchPenaltyReactive = Param(within=NonNegativeReals, mutable=True, rule=make_penalty_rule('q_load_mismatch_cost', 2.))
 
-    model.LoadMismatchPenalty = Param(within=NonNegativeReals, mutable=True, initialize=system.get('load_mismatch_cost', BigPenalty))
-    model.LoadMismatchPenaltyReactive = Param(within=NonNegativeReals, mutable=True, initialize=system.get('q_load_mismatch_cost', BigPenalty/2.))
+    model.ReserveShortfallPenalty = Param(within=NonNegativeReals, mutable=True, rule=make_penalty_rule('reserve_shortfall_cost', 10.))
 
     model.Contingencies = Set(initialize=contingencies.keys())
 
     # leaving this unindexed for now for simpility
     model.SystemContingencyLimitPenalty = Param(within=NonNegativeReals,
-                                          initialize=system.get('contingency_flow_violation_cost', BigPenalty/2.),
+                                          rule=make_penalty_rule('contingency_flow_violation_cost', 2.),
                                           mutable=True)
 
     model.SystemTransmissionLimitPenalty = Param(within=NonNegativeReals,
-                                           initialize=system.get('transmission_flow_violation_cost', BigPenalty/2.),
+                                           rule=make_penalty_rule('transmission_flow_violation_cost', 2.),
                                            mutable=True)
 
     model.SystemInterfaceLimitPenalty = Param(within=NonNegativeReals,
-                                        initialize=system.get('interface_flow_violation_cost', BigPenalty/4.),
+                                        rule=make_penalty_rule('interface_flow_violation_cost', (10/3.)), #3.333
                                         mutable=True)
 
     ##############################################
@@ -258,7 +278,6 @@ def _warn_neg_impedence(m, v, l):
     model.HVDCLineOutOfService = Param(model.HVDCLines, model.TimePeriods, within=Boolean, default=False,
                                        initialize=TimeMapper(dc_branch_attrs.get('planned_outage', dict())))
 
-    _branch_penalties = {}
     _branches_with_slack = []
     for bn, branch in branches.items():
         if 'violation_penalty' in branch:
@@ -269,7 +288,6 @@ def _warn_neg_impedence(m, v, l):
                 if slack_type == SlackType.NONE:
                     logger.warning("Ignoring slacks on individual transmission constraints because SlackType.NONE was specified")
                     break
-                _branch_penalties[bn] = val
                 _branches_with_slack.append(bn)
                 if val <= 0:
                     logger.warning("Branch {} has a non-positive penalty {}, this will cause its limits to be ignored!".format(bn,val))
@@ -280,9 +298,8 @@ def _warn_neg_impedence(m, v, l):
 
     model.BranchLimitPenalty = Param(model.BranchesWithSlack,
                                      within=NonNegativeReals,
-                                     default=value(model.SystemTransmissionLimitPenalty),
-                                     mutable=True,
-                                     initialize=_branch_penalties)
+                                     rule=make_indexed_penalty_rule('branch', model.SystemTransmissionLimitPenalty),
+                                     mutable=True)
 
     ## Interfaces
     model.Interfaces = Set(initialize=interface_attrs['names'])
@@ -311,13 +328,11 @@ def get_interface_line_pairs(m):
 
     model.InterfaceLineOrientation = Param(model.InterfaceLinePairs, initialize=_interface_line_orientation_dict, within=set([-1,0,1]))
 
-    _interface_penalties = {}
     _interfaces_with_slack = []
     for i_n, interface in interfaces.items():
         if 'violation_penalty' in interface:
             val = interface['violation_penalty']
             if val is not None:
-                _interface_penalties[i_n] = val
                 _interfaces_with_slack.append(i_n)
                 if val <= 0:
                     logger.warning("Interface {} has a non-positive penalty {}, this will cause its limits to be ignored!".format(i_n,val))
@@ -328,9 +343,8 @@ def get_interface_line_pairs(m):
 
     model.InterfaceLimitPenalty = Param(model.InterfacesWithSlack,
                                         within=NonNegativeReals,
-                                        default=value(model.SystemInterfaceLimitPenalty),
                                         mutable=True,
-                                        initialize=_interface_penalties)
+                                        rule=make_indexed_penalty_rule('interface', model.SystemInterfaceLimitPenalty))
 
     ##########################################################
     # string indentifiers for the set of thermal generators. #

egret/model_library/unit_commitment/services.py

@@ -11,7 +11,7 @@
 from pyomo.environ import *
 import math
 
-from .uc_utils import add_model_attr, uc_time_helper 
+from .uc_utils import add_model_attr, uc_time_helper, make_penalty_rule
 from .status_vars import _is_relaxed
 
 @add_model_attr('storage_service', requires = {'data_loader': None,
@@ -272,30 +272,50 @@ def _check_for_requirement( requirement ):
         raise Exception('Exception adding ancillary_services! ancillary_services requires one of: garver_3bin_vars, garver_2bin_vars, garver_3bin_relaxed_stop_vars, ALS_state_transition_vars, to be used for the status_vars.')
 
     ## set some penalties by default based on the other model penalties
-    default_reg_pen = value(model.LoadMismatchPenalty+model.ReserveShortfallPenalty)/2.
     ## set these penalties in relation to each other, from higher quality service to lower
+    #################################################################################
+    # penalty costs for constraint violation
+    #
+    # While the user can specify these, by default we base all penalties
+    # off the "load_mismatch_cost", which always has the highest penalty
+    # value (default $1M/MWh). If the user sets "load_mismatch_cost"
+    # at $1000/MWh, the following penalties will be used:
+    #
+    # (defined in params.py)
+    # "q_load_mismatch_cost"              : $500/MVh ("load_mismatch_cost"/2)
+    # "transmission_flow_violation_cost"  : $500/MWh ("load_mismatch_cost"/2)
+    # "contingency_flow_violation_cost"   : $500/MWh ("load_mismatch_cost"/2)
+    # "interface_flow_violation_cost"     : $300/MWh ("load_mismatch_cost"/(10/3))
+    # "reserve_shortfall_cost"            : $100/MWh ("load_mismatch_cost"/10)
+    #
+    # (defined here in services.py)
+    # "regulation_penalty_price"          : $250/MWh ("load_mismatch_cost"/4)
+    # "spinning_reserve_penalty_price"    : $200/MWh ("load_mismatch_cost"/5)
+    # "non_spinning_reserve_penalty_price": $150/MWh ("load_mismatch_cost"/(20/3))
+    # "supplemental_reserve_penalty_price": $125/MWh ("load_mismatch_cost"/8)
+    # "flexible_ramp_penalty_price"       : $110/MWh ("load_mismatch_cost"/(100/11))
+    #
+    # Note these can be overridden by the user specifying the values themselves.
+    # Further, penalties on branch flows and interfaces can be set per-element.
+    ################################################################################
     model.RegulationPenalty = Param(within=NonNegativeReals,
-            initialize=system.get('regulation_penalty_price', default_reg_pen),
+            rule=make_penalty_rule('regulation_penalty_price', 4.),
             mutable=True)
 
-    default_spin_pen = value(model.RegulationPenalty+model.ReserveShortfallPenalty)/2.
     model.SpinningReservePenalty = Param(within=NonNegativeReals, 
-            initialize=system.get('spinning_reserve_penalty_price', default_spin_pen),
+            rule=make_penalty_rule('spinning_reserve_penalty_price', 5.),
             mutable=True)
 
-    default_nspin_pen = value(model.SpinningReservePenalty+model.ReserveShortfallPenalty)/2.
     model.NonSpinningReservePenalty = Param(within=NonNegativeReals,
-            initialize=system.get('non_spinning_reserve_penalty_price', default_nspin_pen),
+            rule=make_penalty_rule('non_spinning_reserve_penalty_price', (20/3.)), #6.667
             mutable=True)
 
-    default_supp_pen = value(model.NonSpinningReservePenalty+model.ReserveShortfallPenalty)/2.
     model.SupplementalReservePenalty = Param(within=NonNegativeReals,
-            initialize=system.get('supplemental_reserve_penalty_price', default_supp_pen),
+            rule=make_penalty_rule('supplemental_reserve_penalty_price', 8.),
             mutable=True)
 
-    default_flex_pen = value(model.NonSpinningReservePenalty+model.SpinningReservePenalty)/2.
     model.FlexRampPenalty = Param(within=NonNegativeReals,
-            initialize=system.get('flexible_ramp_penalty_price', default_flex_pen),
+            rule=make_penalty_rule('flexible_ramp_penalty_price', (100/11.)), #9.09
             mutable=True)
 
     thermal_gen_attrs = md.attributes(element_type='generator', generator_type='thermal')

egret/model_library/unit_commitment/uc_utils.py

@@ -16,14 +16,17 @@
 ## some useful functions and function decorators for building these dynamic models
 from enum import Enum
 from functools import wraps
-from pyomo.environ import Var, quicksum
+from pyomo.environ import Param, Var, quicksum, value
 from pyomo.core.expr.numeric_expr import LinearExpression
+from pyomo.core.base.initializer import ScalarCallInitializer, IndexedCallInitializer
 
 import warnings
 
 import logging
 logger = logging.getLogger('egret.model_library.unit_commitment.uc_utils')
 
+from egret.model_library.transmission.tx_utils import scale_ModelData_to_pu, unscale_ModelData_to_pu
+
 class SlackType(Enum):
     '''
     BUS_BALANCE: Slacks at every bus balance constraint
@@ -124,3 +127,30 @@ def get_linear_expr(*args):
         if not is_var(arg):
             return linear_summation
     return _linear_expression
+
+# Helpers for making penalty factors "commonly" mutable.
+# E.g., change LoadMismatchPenalty and the rest adjust
+# automatically if not directly specified
+def make_penalty_rule(penalty_key, divisor):
+    def penalty_rule(m):
+        return m.model_data.data['system'].get(penalty_key, value(m.LoadMismatchPenalty/divisor))
+    return penalty_rule
+
+def make_indexed_penalty_rule(element_key, base_penalty):
+    def penalty_rule(m, idx):
+        return m.model_data.data['elements'][element_key][idx].get('violation_penalty', base_penalty._rule(m, None))
+    return penalty_rule
+
+def _reconstruct_pyomo_component(component):
+    component.clear()
+    component._constructed = False
+    component.construct()
+
+def reset_unit_commitment_penalties(m):
+    scale_ModelData_to_pu(m.model_data, inplace=True)
+    _reconstruct_pyomo_component(m.LoadMismatchPenalty)
+    for param in m.component_objects(Param):
+        if param.mutable and isinstance(param._rule, (ScalarCallInitializer, IndexedCallInitializer)) \
+                and (param._rule._fcn.__name__ == 'penalty_rule'):
+            _reconstruct_pyomo_component(param)
+    unscale_ModelData_to_pu(m.model_data, inplace=True)

egret/models/tests/uc_test_instances/tiny_uc_12_results.json

@@ -1 +1 @@
-{"elements": {"load": {}, "bus": {"Arne": {"id": "113", "base_kv": 230.0, "matpower_bustype": "ref", "vm": 1.03943, "va": {"data_type": "time_series", "values": [-0.0, -0.0, -0.0, -0.0]}, "v_min": 0.95, "v_max": 1.05, "area": "1", "zone": "14.0", "p_balance_violation": {"data_type": "time_series", "values": [-22.0, -22.0, -0.0, -0.0]}, "pl": {"data_type": "time_series", "values": [0.0, 0.0, 0.0, 0.0]}}}, "generator": {"gen": {"bus": "Arne", "in_service": true, "mbase": 100.0, "pg": {"data_type": "time_series", "values": [22.0, 22.0, 0.0, 0.0]}, "qg": 10.99, "p_min": 22.0, "p_max": 55.00000000000001, "q_min": -15.0, "q_max": 19.0, "ramp_q": 3.7000000000000006, "fuel": "NG", "unit_type": "CT", "area": "3", "zone": "32.0", "generator_type": "thermal", "p_fuel": {"data_type": "fuel_curve", "values": [[22.0, 338.69], [33.0, 434.31], [44.0, 542.93], [55.00000000000001, 652.26]]}, "startup_fuel": [[2.2, 14.574]], "non_fuel_startup_cost": 0.0, "shutdown_cost": 0.0, "agc_capable": true, "p_min_agc": 22.0, "p_max_agc": 55.00000000000001, "ramp_agc": 3.7000000000000006, "ramp_up_60min": 222.00000000000003, "ramp_down_60min": 222.00000000000003, "fuel_cost": 3.88722, "startup_capacity": 22.0, "shutdown_capacity": 22.0, "min_up_time": 2.2, "min_down_time": 2.2, "initial_status": 1, "initial_p_output": 0.0, "initial_q_output": 0.0, "future_status": 0.0, "startup_curve": [], "shutdown_curve": [], "commitment": {"data_type": "time_series", "values": [1, 1, 0, 0]}, "commitment_cost": {"data_type": "time_series", "values": [1316.5625418, 1316.5625418, 0.0, 0.0]}, "production_cost": {"data_type": "time_series", "values": [0.0, 0.0, 0.0, 0.0]}, "headroom": {"data_type": "time_series", "values": [33.00000000000001, 0.0, 0.0, 0.0]}}}, "branch": {}, "interface": {}, "storage": {}, "dc_branch": {}, "zone": {}, "area": {}}, "system": {"name": "RTS-GMLC", "baseMVA": 100.0, "reference_bus": "Arne", "reference_bus_angle": 0, "time_period_length_minutes": 60, "time_keys": ["1", "2", "3", "4"], "total_cost": 442633.1250836}}
+{"elements": {"load": {}, "bus": {"Arne": {"id": "113", "base_kv": 230.0, "matpower_bustype": "ref", "vm": 1.03943, "va": {"data_type": "time_series", "values": [-0.0, -0.0, -0.0, -0.0]}, "v_min": 0.95, "v_max": 1.05, "area": "1", "zone": "14.0", "p_balance_violation": {"data_type": "time_series", "values": [-22.0, -22.0, -0.0, -0.0]}, "pl": {"data_type": "time_series", "values": [0.0, 0.0, 0.0, 0.0]}}}, "generator": {"gen": {"bus": "Arne", "in_service": true, "mbase": 100.0, "pg": {"data_type": "time_series", "values": [22.0, 22.0, 0.0, 0.0]}, "qg": 10.99, "p_min": 22.0, "p_max": 55.00000000000001, "q_min": -15.0, "q_max": 19.0, "ramp_q": 3.7000000000000006, "fuel": "NG", "unit_type": "CT", "area": "3", "zone": "32.0", "generator_type": "thermal", "p_fuel": {"data_type": "fuel_curve", "values": [[22.0, 338.69], [33.0, 434.31], [44.0, 542.93], [55.00000000000001, 652.26]]}, "startup_fuel": [[2.2, 14.574]], "non_fuel_startup_cost": 0.0, "shutdown_cost": 0.0, "agc_capable": true, "p_min_agc": 22.0, "p_max_agc": 55.00000000000001, "ramp_agc": 3.7000000000000006, "ramp_up_60min": 222.00000000000003, "ramp_down_60min": 222.00000000000003, "fuel_cost": 3.88722, "startup_capacity": 22.0, "shutdown_capacity": 22.0, "min_up_time": 2.2, "min_down_time": 2.2, "initial_status": 1, "initial_p_output": 0.0, "initial_q_output": 0.0, "future_status": 0.0, "startup_curve": [], "shutdown_curve": [], "commitment": {"data_type": "time_series", "values": [1, 1, 0, 0]}, "commitment_cost": {"data_type": "time_series", "values": [1316.5625418, 1316.5625418, 0.0, 0.0]}, "production_cost": {"data_type": "time_series", "values": [0.0, 0.0, 0.0, 0.0]}, "headroom": {"data_type": "time_series", "values": [33.00000000000001, 0.0, 0.0, 0.0]}}}, "branch": {}, "interface": {}, "storage": {}, "dc_branch": {}, "zone": {}, "area": {}}, "system": {"name": "RTS-GMLC", "baseMVA": 100.0, "reference_bus": "Arne", "reference_bus_angle": 0, "time_period_length_minutes": 60, "time_keys": ["1", "2", "3", "4"], "total_cost": 44002633.1250836}}