Adding user-configurable slack type

egret/common/lazy_ptdf_utils.py

@@ -693,7 +693,7 @@ def _add_contingency_violations(lazy_violations, flows, mb, md, solver, ptdf_opt
         contingencies_monitored.append((cn, i_b))
         if new_slacks:
             m = model
-            obj_coef = pyo.value(m.TimePeriodLengthHours*m.ContingencyLimitPenalty)
+            obj_coef = pyo.value(m.TimePeriodLengthHours*m.SystemContingencyLimitPenalty)
 
             if persistent_solver:
                 m_model = m.model()

egret/model_library/transmission/branch.py

@@ -1090,7 +1090,7 @@ def declare_ineq_p_contingency_branch_thermal_bounds(model, index_set,
                 pos_slack = m.pfc_slack_pos[contingency_name, branch_name]
                 uc_model = slack_cost_expr.parent_block()
                 slack_cost_expr.expr += (uc_model.TimePeriodLengthHours
-                                         * uc_model.ContingencyLimitPenalty
+                                         * uc_model.SystemContingencyLimitPenalty
                                          * (neg_slack + pos_slack) )
                 assert len(m.pfc_slack_pos) == len(m.pfc_slack_neg)
             else:

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
+from .uc_utils import add_model_attr, uc_time_helper, SlackType
 
 component_name = 'data_loader'
 
@@ -70,7 +70,7 @@ def initial_time_periods_offline_rule(m, g):
     model.InitialTimePeriodsOffLine = Param(model.ThermalGenerators, within=NonNegativeIntegers, initialize=initial_time_periods_offline_rule, mutable=True)
 
 @add_model_attr(component_name)
-def load_params(model, model_data):
+def load_params(model, model_data, slack_type):
 
     '''
     This loads unit commitment params from a GridModel object
@@ -189,6 +189,34 @@ def load_params(model, model_data):
                                         initialize={'Stage_1':model.TimePeriods, 'Stage_2': list() } ) 
     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']
+
+    model.ReserveShortfallPenalty = Param(within=NonNegativeReals, default=ModeratelyBigPenalty, mutable=True, initialize=system.get('reserve_shortfall_cost', ModeratelyBigPenalty))
+
+    BigPenalty = 1e4*system['baseMVA']
+
+    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.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.),
+                                          mutable=True)
+
+    model.SystemTransmissionLimitPenalty = Param(within=NonNegativeReals,
+                                           initialize=system.get('transmission_flow_violation_cost', BigPenalty/2.),
+                                           mutable=True)
+
+    model.SystemInterfaceLimitPenalty = Param(within=NonNegativeReals,
+                                        initialize=system.get('interface_flow_violation_cost', BigPenalty/4.),
+                                        mutable=True)
 
     ##############################################
     # Network definition (S)
@@ -230,18 +258,26 @@ 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 = dict()
-    _md_violation_penalties = branch_attrs.get('violation_penalty')
-    if _md_violation_penalties is not None:
-        for i, val in _md_violation_penalties.items():
+    _branch_penalties = {}
+    _branches_with_slack = []
+    for bn, branch in branches.items():
+        if 'violation_penalty' in branch:
+            val = branch['violation_penalty']
             if val is not None:
-                _branch_penalties[i] = val
+                _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(i,val))
+                    logger.warning("Branch {} has a non-positive penalty {}, this will cause its limits to be ignored!".format(bn,val))
+        elif slack_type == SlackType.TRANSMISSION_LIMITS:
+            _branches_with_slack.append(bn)
 
-    model.BranchesWithSlack = Set(within=model.TransmissionLines, initialize=_branch_penalties.keys())
+    model.BranchesWithSlack = Set(within=model.TransmissionLines, initialize=_branches_with_slack)
 
-    model.BranchLimitPenalty = Param(model.BranchesWithSlack, within=NonNegativeReals, initialize=_branch_penalties)
+    model.BranchLimitPenalty = Param(model.BranchesWithSlack,
+                                     within=NonNegativeReals,
+                                     default=value(model.SystemTransmissionLimitPenalty),
+                                     mutable=True,
+                                     initialize=_branch_penalties)
 
     ## Interfaces
     model.Interfaces = Set(initialize=interface_attrs['names'])
@@ -270,18 +306,26 @@ def get_interface_line_pairs(m):
 
     model.InterfaceLineOrientation = Param(model.InterfaceLinePairs, initialize=_interface_line_orientation_dict, within=set([-1,0,1]))
 
-    _interface_penalties = dict()
-    _md_violation_penalties = interface_attrs.get('violation_penalty')
-    if _md_violation_penalties is not None:
-        for i, val in _md_violation_penalties.items():
+    _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] = val
+                _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,val))
+                    logger.warning("Interface {} has a non-positive penalty {}, this will cause its limits to be ignored!".format(i_n,val))
+        elif slack_type == SlackType.TRANSMISSION_LIMITS:
+            _interfaces_with_slack.append(bn)
 
-    model.InterfacesWithSlack = Set(within=model.Interfaces, initialize=_interface_penalties.keys())
+    model.InterfacesWithSlack = Set(within=model.Interfaces, initialize=_interfaces_with_slack)
 
-    model.InterfaceLimitPenalty = Param(model.InterfacesWithSlack, within=NonNegativeReals, initialize=_interface_penalties)
+    model.InterfaceLimitPenalty = Param(model.InterfacesWithSlack,
+                                        within=NonNegativeReals,
+                                        default=value(model.SystemInterfaceLimitPenalty),
+                                        mutable=True,
+                                        initialize=_interface_penalties)
 
     ##########################################################
     # string indentifiers for the set of thermal generators. #
@@ -1257,23 +1301,6 @@ def power_generation_piecewise_points_rule(m, g):
 
     ## END PRODUCTION COST CALCULATIONS
 
-    #########################################
-    # penalty costs for constraint violation #
-    #########################################
-
-    ModeratelyBigPenalty = 1e3*system['baseMVA']
-
-    model.ReserveShortfallPenalty = Param(within=NonNegativeReals, default=ModeratelyBigPenalty, mutable=True, initialize=system.get('reserve_shortfall_cost', ModeratelyBigPenalty))
-
-    BigPenalty = 1e4*system['baseMVA']
-
-    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.Contingencies = Set(initialize=contingencies.keys())
-
-    # leaving this unindexed for now for simpility
-    model.ContingencyLimitPenalty = Param(within=NonNegativeReals, initialize=system.get('contingency_flow_violation_cost', BigPenalty/2.), mutable=True)
 
     #
     # STORAGE parameters

egret/model_library/unit_commitment/power_balance.py

@@ -12,7 +12,7 @@
 from pyomo.core.expr.numeric_expr import LinearExpression
 import math
 
-from .uc_utils import add_model_attr
+from .uc_utils import add_model_attr, SlackType
 from .power_vars import _add_reactive_power_vars
 from .generation_limits import _add_reactive_limits
 
@@ -591,39 +591,44 @@ def qg_expr_rule(block,b):
     return qg_expr_rule
 
 ## Defines generic interface for egret tramsmission models
-def _add_egret_power_flow(model, network_model_builder, reactive_power=False, slacks=True):
+def _add_egret_power_flow(model, network_model_builder,
+                          reactive_power=False,
+                          slack_type=SlackType.TRANSMISSION_LIMITS):
 
     ## save flag for objective
     model.reactive_power = reactive_power
 
-    system_load_mismatch = (network_model_builder in \
+    if (slack_type==SlackType.BUS_BALANCE) and \
+            (network_model_builder in \
                             [_copperplate_approx_network_model, \
                              _copperplate_relax_network_model, \
-                            ]
-                           )
+                            ]):
+        # only one slack as there's only a
+        # single power-balance constraint
+        slack_type = SlackType.TRANSMISSION_LIMITS
 
-    if slacks:
-        if system_load_mismatch:
-            _add_system_load_mismatch(model)
-        else:
-            _add_load_mismatch(model)
+    if slack_type == SlackType.BUS_BALANCE:
+        _add_load_mismatch(model)
+    elif slack_type == SlackType.TRANSMISSION_LIMITS:
+        _add_system_load_mismatch(model)
+    elif slack_type == SlackType.NONE:
+        _add_blank_load_mismatch(model)
     else:
-        if system_load_mismatch:
-            _add_blank_system_load_mismatch(model)
-        else:
-            _add_blank_load_mismatch(model)
+        raise ValueError(f"Unrecognized slack_type: {slack_type}")
 
     _add_hvdc(model)
 
     if reactive_power:
-        if system_load_mismatch:
-            raise Exception("Need to implement system mismatch for reactive power")
         _add_reactive_power_vars(model)
         _add_reactive_limits(model)
-        if slacks:
+        if slack_type == SlackType.BUS_BALANCE:
             _add_q_load_mismatch(model)
-        else:
+        elif slack_type == SlackType.TRANSMISSION_LIMITS:
+            raise NotImplementedError("Need to implement transmission slacks for reactive power")
+        elif slack_type == SlackType.NONE:
             _add_blank_q_load_mistmatch(model)
+        else:
+            raise ValueError(f"Unrecognized slack_type: {slack_type}")
 
     # for transmission network
     model.TransmissionBlock = Block(model.TimePeriods, concrete=True)
@@ -653,40 +658,44 @@ def _add_egret_power_flow(model, network_model_builder, reactive_power=False, sl
                                             'non_dispatchable_vars': None,
                                             'storage_service': None,
                                             })
-def copperplate_power_flow(model, slacks=True):
-    _add_egret_power_flow(model, _copperplate_approx_network_model, reactive_power=False, slacks=slacks)
+def copperplate_power_flow(model, slack_type=SlackType.TRANSMISSION_LIMITS):
+    _add_egret_power_flow(model, _copperplate_approx_network_model, reactive_power=False,
+                            slack_type=slack_type)
 
 @add_model_attr(component_name, requires = {'data_loader': None,
                                             'power_vars': None,
                                             'non_dispatchable_vars': None,
                                             'storage_service': None,
                                             })
-def copperplate_relaxed_power_flow(model, slacks=True):
-    _add_egret_power_flow(model, _copperplate_relax_network_model, reactive_power=False, slacks=slacks)
+def copperplate_relaxed_power_flow(model, slack_type=SlackType.TRANSMISSION_LIMITS):
+    _add_egret_power_flow(model, _copperplate_relax_network_model, reactive_power=False,
+                            slack_type=slack_type)
 
 @add_model_attr(component_name, requires = {'data_loader': None,
                                             'power_vars': None,
                                             'non_dispatchable_vars': None,
                                             'storage_service': None,
                                             })
-def ptdf_power_flow(model, slacks=True):
-    _add_egret_power_flow(model, _ptdf_dcopf_network_model, reactive_power=False, slacks=slacks)
+def ptdf_power_flow(model, slack_type=SlackType.TRANSMISSION_LIMITS):
+    _add_egret_power_flow(model, _ptdf_dcopf_network_model, reactive_power=False,
+                            slack_type=slack_type)
 
 @add_model_attr(component_name, requires = {'data_loader': None,
                                             'power_vars': None,
                                             'non_dispatchable_vars': None,
                                             'storage_service': None,
                                             })
-def btheta_power_flow(model, slacks=True):
-    _add_egret_power_flow(model, _btheta_dcopf_network_model, reactive_power=False, slacks=slacks)
+def btheta_power_flow(model, slack_type=SlackType.TRANSMISSION_LIMITS):
+    _add_egret_power_flow(model, _btheta_dcopf_network_model, reactive_power=False,
+                            slack_type=slack_type)
 
 
 @add_model_attr(component_name, requires = {'data_loader': None,
                                             'power_vars': None,
                                             'non_dispatchable_vars': None,
                                             'storage_service': None,
                                             })
-def power_balance_constraints(model, slacks=True):
+def power_balance_constraints(model, slack_type=SlackType.TRANSMISSION_LIMITS):
     '''
     adds the demand and network constraints to the model
     '''
@@ -767,12 +776,16 @@ def interface_violation_cost_rule(m,t):
     # for contingency violation costs at a time step
     # This model does not work with contingencies, but we need the Expression for the objective
     model.ContingencyViolationCost = Expression(model.TimePeriods, rule=lambda m,t:0.)
-    
-    if slacks:
+
+    if slack_type == SlackType.BUS_BALANCE:
         _add_load_mismatch(model)
-    else:
+    elif slack_type == SlackType.TRANSMISSION_LIMITS:
+        _add_system_load_mismatch(model)
+    elif slack_type == SlackType.NONE:
         _add_blank_load_mismatch(model)
-
+    else:
+        raise ValueError(f"Unrecognized slack_type: {slack_type}")
+
     # Power balance at each node (S)
     def power_balance(m, b, t):
         # bus b, time t (S)

egret/model_library/unit_commitment/uc_model_generator.py

@@ -18,6 +18,7 @@
         uptime_downtime, startup_costs, \
         services, power_balance, reserve_requirement, \
         objective, fuel_supply, fuel_consumption, security_constraints
+from .uc_utils import SlackType
 from egret.model_library.transmission.tx_utils import scale_ModelData_to_pu
 from collections import namedtuple
 import pyomo.environ as pe
@@ -37,7 +38,7 @@
                              ]
                             )
 
-def generate_model( model_data, uc_formulation, relax_binaries=False, ptdf_options=None, PTDF_matrix_dict=None ):
+def generate_model( model_data, uc_formulation, relax_binaries=False, ptdf_options=None, PTDF_matrix_dict=None, slack_type=SlackType.BUS_BALANCE):
     """
     returns a UC uc_formulation as an abstract model with the 
     components specified in a UCFormulation, with the option
@@ -57,6 +58,10 @@ def generate_model( model_data, uc_formulation, relax_binaries=False, ptdf_optio
         Dictionary of egret.data.ptdf_utils.PTDFMatrix objects for use in model construction.
         WARNING: Nearly zero checking is done on the consistency of this object with the
                  model_data. Use with extreme caution!
+    slack_type : SlackType, optional
+        Types of slacks to use in the unit commitment model. By default,
+        a global load/generation mismatch slack is placed at the reference
+        bus and all transmission limits are enforced with soft constraints.
 
     Returns
     -------
@@ -66,7 +71,7 @@ def generate_model( model_data, uc_formulation, relax_binaries=False, ptdf_optio
 
     md = model_data.clone_in_service()
     scale_ModelData_to_pu(md, inplace=True)
-    return _generate_model( md, *_get_formulation_from_UCFormulation( uc_formulation ), relax_binaries , ptdf_options, PTDF_matrix_dict )
+    return _generate_model( md, *_get_formulation_from_UCFormulation( uc_formulation ), relax_binaries , ptdf_options, PTDF_matrix_dict, slack_type )
 
 def _generate_model( model_data,
                     _status_vars,
@@ -84,6 +89,7 @@ def _generate_model( model_data,
                     _relax_binaries = False,
                     _ptdf_options = None,
                     _PTDF_matrix_dict = None,
+                    _slack_type = SlackType.TRANSMISSION_LIMITS,
                     ):
 
     model = pe.ConcreteModel()
@@ -114,7 +120,7 @@ def _generate_model( model_data,
     ## to relax binaries
     model.relax_binaries = _relax_binaries
 
-    params.load_params(model, model_data)
+    params.load_params(model, model_data, _slack_type)
     getattr(status_vars, _status_vars)(model)
     getattr(power_vars, _power_vars)(model)
     getattr(reserve_vars, _reserve_vars)(model)
@@ -126,7 +132,7 @@ def _generate_model( model_data,
     getattr(startup_costs, _startup_costs)(model)
     services.storage_services(model)
     services.ancillary_services(model)
-    getattr(power_balance, _power_balance)(model)
+    getattr(power_balance, _power_balance)(model, _slack_type)
     getattr(reserve_requirement, _reserve_requirement)(model)
 
     if 'fuel_supply' in model_data.data['elements'] and bool(model_data.data['elements']['fuel_supply']):

egret/model_library/unit_commitment/uc_utils.py

@@ -14,6 +14,7 @@
 """
 
 ## 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.core.expr.numeric_expr import LinearExpression
@@ -23,6 +24,16 @@
 import logging
 logger = logging.getLogger('egret.model_library.unit_commitment.uc_utils')
 
+class SlackType(Enum):
+    '''
+    BUS_BALANCE: Slacks at every bus balance constraint
+    TRANSMISSION_LIMITS: Slacks at the reference bus and every transmission limit
+    NONE: Slacks nowhere (model may be infeasible)
+    '''
+    BUS_BALANCE = 1
+    TRANSMISSION_LIMITS = 2
+    NONE = 3
+
 def add_model_attr(attr, requires = {}):
     def actual_decorator(func):
         @wraps(func)