Merge pull request #238 from bknueven/priced_load

Priced load and renewables

egret/common/lazy_ptdf_utils.py

@@ -196,6 +196,11 @@ def _add_violations( self, name, other_name, viol_array, viol_indices):
             # to be a parallel constraint.
             # If we haven't added any constraints yet
             # any(()) is False, so this won't fire
+            # TODO: since this object gets re-created each iteration,
+            #       and the violations_store just has violations we're
+            #       adding *this* iteration, it's possible to add parallel
+            #       lines, which may not be necessary in may cases (e.g.,
+            #       when the line is binding but not over the limit)
             close_to_existing = any( math.isclose( val, existing ) for existing in self.violations_store.values() )
             if close_to_existing:
                 viol_indices.pop(idx)
@@ -676,7 +681,6 @@ def _add_interface_violations(lazy_violations, flows, mb, md, solver, ptdf_optio
                                                         obj_coef*mb.pfi_slack_neg[i_n] )
         if persistent_solver:
             solver.add_constraint(constr[i_n])
-        print(f"adding constraint for interface {i_n}")
 
 def _add_contingency_violations(lazy_violations, flows, mb, md, solver, ptdf_options,
                                 PTDF, model, baseMVA, persistent_solver, rel_ptdf_tol, abs_ptdf_tol,
@@ -710,51 +714,38 @@ def _add_contingency_violations(lazy_violations, flows, mb, md, solver, ptdf_opt
         if persistent_solver:
             solver.add_constraint(constr[cn,bn])
 
-## helper for generating pf
-def _iter_over_initial_set(branches, branches_in_service, PTDF):
+
+def add_initial_monitored_constraints(mb, md, branches_in_service, ptdf_options, PTDF, time=None): 
+
+    viol_not_lazy = set()
     for bn in branches_in_service:
-        branch = branches[bn]
+        branch = md.data['elements']['branch'][bn] 
         if 'lazy' in branch and not branch['lazy']:
             if bn in PTDF.branchname_to_index_masked_map:
-                i = PTDF.branchname_to_index_masked_map[bn]
-                yield i, bn
+                viol_not_lazy.add(PTDF.branchname_to_index_masked_map[bn])
             else:
                 logger.warning("Branch {0} has flag 'lazy' set to False but is excluded from monitored set based on kV limits".format(bn))
 
-### initial monitored set adder
-def add_initial_monitored_branches(mb, branches, branches_in_service, ptdf_options, PTDF):
-    ## static information between runs
-    rel_ptdf_tol = ptdf_options['rel_ptdf_tol']
-    abs_ptdf_tol = ptdf_options['abs_ptdf_tol']
+    int_viol_not_lazy = set()
+    if 'interface' in md.data['elements']:
+        for i_n, interface in md.data['elements']['interface'].items():
+            if 'lazy' in interface and not interface['lazy']:
+                int_viol_not_lazy.add(PTDF.interfacename_to_index_map[i_n])
 
-    constr = mb.ineq_pf_branch_thermal_bounds
-    viol_in_mb = mb._idx_monitored
-    for i, bn in _iter_over_initial_set(branches, branches_in_service, PTDF):
-        thermal_limit = PTDF.branch_limits_array_masked[i]
-        mb.pf[bn] = libbranch.get_power_flow_expr_ptdf_approx(mb, bn, PTDF, abs_ptdf_tol=abs_ptdf_tol, rel_ptdf_tol=rel_ptdf_tol)
-        constr[bn], _ = _generate_branch_thermal_bounds(mb, bn, thermal_limit)
-        viol_in_mb.append(i)
+    # not easy to support in the current
+    # set-up, as 'lazy' would need to be
+    # set on a branch, branch basis
+    cont_viol_not_lazy = set()
 
-def _iter_over_initial_set_interfaces(interfaces, PTDF):
-    for i_n, interface in interfaces.items():
-        if 'lazy' in interface and not interface['lazy']:
-            i = PTDF.interfacename_to_index_map[i_n]
-            yield i, i_n
+    #blank flows
+    flows = _CalculatedFlows()
+    lazy_violations = _LazyViolations(branch_lazy_violations=viol_not_lazy,
+                                interface_lazy_violations=int_viol_not_lazy,
+                                contingency_lazy_violations=cont_viol_not_lazy)
 
-### initial monitored set adder
-def add_initial_monitored_interfaces(mb, interfaces, ptdf_options, PTDF):
-    ## static information between runs
-    rel_ptdf_tol = ptdf_options['rel_ptdf_tol']
-    abs_ptdf_tol = ptdf_options['abs_ptdf_tol']
+    add_violations(lazy_violations, flows, mb, md, None,
+            ptdf_options, PTDF, time=time, prepend_str="[Initial Set] ", obj_multi=None)
 
-    constr = mb.ineq_pf_interface_bounds
-    int_viol_in_mb = mb._interfaces_monitored
-    for i, i_n in _iter_over_initial_set_interfaces(interfaces, PTDF):
-        minimum_limit = PTDF.interface_min_limits[i]
-        maximum_limit = PTDF.interface_max_limits[i]
-        mb.pfi[i_n] = libbranch.get_power_flow_interface_expr_ptdf(mb, i_n, PTDF, abs_ptdf_tol=abs_ptdf_tol, rel_ptdf_tol=rel_ptdf_tol)
-        constr[i_n], _ = _generate_interface_bounds(mb, i_n, minimum_limit, maximum_limit)
-        int_viol_in_mb.append(i)
 
 def copy_active_to_next_time(m, b_next, PTDF_next, slacks, slacks_I, slacks_C):
     active_slack_tol = m._ptdf_options['active_flow_tol']

egret/data/model_data.py

@@ -221,12 +221,16 @@ def elements(self, element_type, **kwargs):
 
         # additional attributes have been specified
         for name, elem in self.data['elements'][element_type].items():
-            include = True
             for k, v in kwargs.items():
-                if k not in elem or elem[k] != v:
-                    include = False
+                if k not in elem:
                     break
-            if include:
+                if type(v) is tuple:
+                    if elem[k] not in v:
+                        break
+                else:
+                    if elem[k] != v:
+                        break
+            else: # no break
                 yield name, elem
 
     def attributes(self, element_type, **kwargs):

egret/data/tests/test_model_data.py

@@ -34,7 +34,8 @@
                 {
                     'B1': {'bus_type': 'PV'},
                     'B2': {'bus_type': 'PQ'},
-                    'B3': {'bus_type': 'PQ'}
+                    'B3': {'bus_type': 'PQ'},
+                    'B4': {'bus_type': 'ref'},
                 },
             'branch':
                 {
@@ -135,17 +136,23 @@ def test_elements():
             assert e["generator_type"] == 'solar'
 
     buses = dict(md.elements(element_type='bus'))
-    assert len(buses) == 3
+    assert len(buses) == 4
     assert 'B1' in buses
     assert 'B2' in buses
     assert 'B3' in buses
+    assert 'B4' in buses
     assert buses['B1']['bus_type'] == 'PV'
 
     buses = dict(md.elements(element_type='bus', bus_type='PQ'))
     assert len(buses) == 2
     assert 'B2' in buses
     assert 'B3' in buses
 
+    buses = dict(md.elements(element_type='bus', bus_type=('PV','ref')))
+    assert len(buses) == 2
+    assert 'B1' in buses
+    assert 'B4' in buses
+
 def test_attributes():
     md = ModelData(dict(testdata))
 

egret/model_library/transmission/branch.py

@@ -972,7 +972,7 @@ def declare_ineq_p_branch_thermal_lbub(model, index_set,
                                             )
                 else:
                     expr = m.pf[branch_name] - pos_slack
-                m.ineq_pf_branch_thermal_lb[branch_name] = \
+                m.ineq_pf_branch_thermal_ub[branch_name] = \
                     (None, expr, p_thermal_limits[branch_name])
             else:
                 m.ineq_pf_branch_thermal_ub[branch_name] = \

egret/model_library/transmission/tx_utils.py

@@ -312,6 +312,8 @@ def load_shed_limit(load, gens, gen_mins):
                                                       'q_load',
                                                       'p_load_shed',
                                                       'q_load_shed',
+                                                      'p_price',
+                                                      'q_price',
                                                      ],
                        ('element_type','branch', None) : [
                                                        'rating_long_term',

egret/model_library/unit_commitment/non_dispatchable_vars.py

@@ -23,4 +23,9 @@ def nd_bounds_rule(m,n,t):
         return (m.MinNondispatchablePower[n,t], m.MaxNondispatchablePower[n,t])
     model.NondispatchablePowerUsed = Var(model.AllNondispatchableGenerators, model.TimePeriods, within=Reals, bounds=nd_bounds_rule)
 
+    def _nd_power_cost_rule(m,n,t):
+        return m.TimePeriodLengthHours*m.NondispatchableMarginalCost[n,t]*m.NondispatchablePowerUsed[n,t]
+    model.NondispatchableProductionCost = Expression(model.AllNondispatchableGenerators, model.TimePeriods,
+                                                    rule=_nd_power_cost_rule)
+
     return

egret/model_library/unit_commitment/objective.py

@@ -130,7 +130,11 @@ def generation_stage_cost_expression_rule(m, st):
               sum(m.BranchViolationCost[t] for t in m.GenerationTimeInStage[st]) + \
               sum(m.InterfaceViolationCost[t] for t in m.GenerationTimeInStage[st]) + \
               sum(m.ContingencyViolationCost[t] for t in m.GenerationTimeInStage[st]) + \
-              sum(m.StorageCost[s,t] for s in m.Storage for t in m.GenerationTimeInStage[st])
+              sum(m.StorageCost[s,t] for s in m.Storage for t in m.GenerationTimeInStage[st]) + \
+              sum(m.PriceResponsiveLoadCost[l,t] for l in m.PriceResponsiveLoad
+                      for t in m.GenerationTimeInStage[st]) + \
+              sum(m.NondispatchableProductionCost[n,t] for n in m.AllNondispatchableGenerators
+                      for t in m.GenerationTimeInStage[st])
         if m.reactive_power:
             cc += sum(m.LoadMismatchCostReactive[t] for t in m.GenerationTimeInStage[st])
         if m.security_constraints:

egret/model_library/unit_commitment/params.py

@@ -99,7 +99,7 @@ def load_params(model, model_data):
 
     loads = dict(md.elements(element_type='load'))
     thermal_gens = dict(md.elements(element_type='generator', generator_type='thermal'))
-    renewable_gens = dict(md.elements(element_type='generator', generator_type='renewable'))
+    renewable_gens = dict(md.elements(element_type='generator', generator_type=('renewable','virtual')))
     buses = dict(md.elements(element_type='bus'))
     shunts = dict(md.elements(element_type='shunt'))
     branches = dict(md.elements(element_type='branch'))
@@ -109,7 +109,7 @@ def load_params(model, model_data):
     dc_branches = dict(md.elements(element_type='dc_branch'))
 
     thermal_gen_attrs = md.attributes(element_type='generator', generator_type='thermal')
-    renewable_gen_attrs = md.attributes(element_type='generator', generator_type='renewable')
+    renewable_gen_attrs = md.attributes(element_type='generator', generator_type=('renewable','virtual'))
     bus_attrs = md.attributes(element_type='bus')
     branch_attrs = md.attributes(element_type='branch')
     load_attrs = md.attributes(element_type='load')
@@ -340,14 +340,17 @@ def init_quick_start_generators(m):
     for lname, load in loads.items():
         load_time = TimeMapper(load['p_load'])
         bus = load['bus']
+        # priced loads will be handled separately
+        if 'p_price' in load and load['p_price'] is not None:
+            continue
         if isinstance(bus, dict):
             assert bus['data_type'] == 'load_distribution_factor'
             for bn, multi in bus['values'].items():
-                for t in model.TimePeriods:
-                    bus_loads[bn, t] += multi*load_time[t]
+                for t, load in load_time.items():
+                    bus_loads[bn, t] += multi*load
         else:
-            for t in model.TimePeriods:
-                bus_loads[bus, t] += load_time[t]
+            for t, load in load_time.items():
+                bus_loads[bus, t] += load
     model.Demand = Param(model.Buses, model.TimePeriods, initialize=bus_loads, mutable=True)
 
     def calculate_total_demand(m, t):
@@ -362,6 +365,34 @@ def warn_about_negative_demand_rule(m, b, t):
 
     if warn_neg_load:
         model.WarnAboutNegativeDemand = BuildAction(model.Buses, model.TimePeriods, rule=warn_about_negative_demand_rule)
+
+    _price_responsive_load_by_bus = {}
+    _price_responsive_load_attrs = {'names': [], 'p_price': {}, 'p_load': {}}
+    for ln, load in loads.items():
+        if 'p_price' in load and load['p_price'] is not None:
+            bus = load['bus']
+            if bus in _price_responsive_load_by_bus:
+                _price_responsive_load_by_bus[bus].append(ln)
+            else:
+                _price_responsive_load_by_bus[bus]= [ln]
+            _price_responsive_load_attrs['names'].append(ln)
+            _price_responsive_load_attrs['p_price'][ln] = load['p_price']
+            _price_responsive_load_attrs['p_load'][ln] = load['p_load']
+
+    model.PriceResponsiveLoadAtBus = Set(model.Buses,
+            initialize=lambda m,b : _price_responsive_load_by_bus[b] if b in _price_responsive_load_by_bus else ())
+
+    model.PriceResponsiveLoad = Set(initialize=_price_responsive_load_attrs['names'])
+
+    model.PriceResponsiveLoadPrice = Param(model.PriceResponsiveLoad,
+                                           model.TimePeriods,
+                                           within=Reals,
+                                           initialize=TimeMapper(_price_responsive_load_attrs['p_price']))
+
+    model.PriceResponsiveLoadDemand = Param(model.PriceResponsiveLoad,
+                                            model.TimePeriods,
+                                            within=NonNegativeReals,
+                                            initialize=TimeMapper(_price_responsive_load_attrs['p_load']))
 
     ##################################################################
     # the global system reserve, for each time period. units are MW. #
@@ -450,6 +481,13 @@ def maximum_nd_output_validator(m, v, g, t):
                                             validate=maximum_nd_output_validator,
                                             initialize=TimeMapper(renewable_gen_attrs.get('p_max', dict())))
 
+    model.NondispatchableMarginalCost = Param(model.AllNondispatchableGenerators,
+                                            model.TimePeriods,
+                                            within=Reals, # more permissive; e.g. CSP
+                                            default=0.0,
+                                            mutable=True,
+                                            initialize=TimeMapper(renewable_gen_attrs.get('p_cost', dict())))
+
     #################################################
     # generator ramp up/down rates. units are MW/h. #
     # IMPORTANT: Generator ramp limits can exceed   #

egret/model_library/unit_commitment/power_balance.py

@@ -212,10 +212,8 @@ def _ptdf_dcopf_network_model(block,tm):
         lpu.add_monitored_flow_tracker(block)
 
         ### add initial branches to monitored set
-        lpu.add_initial_monitored_branches(block, branches, branches_in_service, ptdf_options, PTDF)
-
-        ### add initial interfaces to monitored set
-        lpu.add_initial_monitored_interfaces(block, interfaces, ptdf_options, PTDF)
+        lpu.add_initial_monitored_constraints(block, m.model_data, branches_in_service,
+                ptdf_options, PTDF, time=tm)
 
     else: ### add all the dense constraints
         if contingencies:
@@ -578,6 +576,7 @@ def pg_expr_rule(block,b):
                 + sum(m.NondispatchablePowerUsed[g, t] for g in m.NondispatchableGeneratorsAtBus[b]) \
                 + sum(m.HVDCLinePower[k,t] for k in m.HVDCLinesTo[b]) \
                 - sum(m.HVDCLinePower[k,t] for k in m.HVDCLinesFrom[b]) \
+                - sum(m.PriceResponsiveLoadServed[l,t] for l in m.PriceResponsiveLoadAtBus[b]) \
                 + m.LoadGenerateMismatch[b,t]
     return pg_expr_rule
 
@@ -784,6 +783,7 @@ def power_balance(m, b, t):
                 - sum(m.LinePower[l,t] for l in m.LinesFrom[b]) \
                 + sum(m.HVDCLinePower[k,t] for k in m.HVDCLinesTo[b]) \
                 - sum(m.HVDCLinePower[k,t] for k in m.HVDCLinesFrom[b]) \
+                - sum(m.PriceResponsiveLoadServed[l,t] for l in m.PriceResponsiveLoadAtBus[b]) \
                 + m.LoadGenerateMismatch[b,t] \
                 == m.Demand[b, t] 
     model.PowerBalance = Constraint(model.Buses, model.TimePeriods, rule=power_balance)

egret/model_library/unit_commitment/services.py

@@ -164,6 +164,23 @@ def storage_cost_rule(m, s, t):
     return
 ## end storage_services
 
+@add_model_attr('load_service', requires = {'data_loader': None})
+def load_services(model):
+    '''
+    For price-responsive load
+    '''
+
+    model.PriceResponsiveLoadServed = Var(model.PriceResponsiveLoad,
+                                          model.TimePeriods,
+                                          within=NonNegativeReals,
+                                          bounds=lambda m,l,t:(0, m.PriceResponsiveLoadDemand[l,t]))
+
+    def _price_responsive_load_cost_rule(m,l,t):
+        return -1*m.TimePeriodLengthHours*m.PriceResponsiveLoadPrice[l,t]*m.PriceResponsiveLoadServed[l,t]
+    model.PriceResponsiveLoadCost = Expression(model.PriceResponsiveLoad, model.TimePeriods,
+                                               rule=_price_responsive_load_cost_rule)
+
+    return
 
 ## NOTE: when moving to a Real-Time market, ramping limits need to be also considered
 ##       It seems MISO did not in its DA as of 2009 [1], but definitely does in RT as of 2016 [2].

egret/model_library/unit_commitment/uc_model_generator.py

@@ -125,6 +125,7 @@ def _generate_model( model_data,
     getattr(uptime_downtime, _uptime_downtime)(model)
     getattr(startup_costs, _startup_costs)(model)
     services.storage_services(model)
+    services.load_services(model)
     services.ancillary_services(model)
     getattr(power_balance, _power_balance)(model)
     getattr(reserve_requirement, _reserve_requirement)(model)

egret/models/dcopf.py

@@ -327,7 +327,7 @@ def create_ptdf_dcopf_model(model_data, include_feasibility_slack=False, base_po
         lpu.add_monitored_flow_tracker(model)
 
         ### add initial branches to monitored set
-        lpu.add_initial_monitored_branches(model, branches, branches_idx, ptdf_options, PTDF)
+        lpu.add_initial_monitored_constraints(model, md, branches_idx, ptdf_options, PTDF)
 
     else:
         p_max = {k: branches[k]['rating_long_term'] for k in branches.keys()}

egret/models/scopf.py

@@ -166,7 +166,7 @@ def create_scopf_model(model_data, include_feasibility_slack=False,
     lpu.add_monitored_flow_tracker(model)
 
     ### add initial branches to monitored set
-    lpu.add_initial_monitored_branches(model, branches, branches_idx, ptdf_options, PTDF)
+    lpu.add_initial_monitored_constraints(model, md, branches_idx, ptdf_options, PTDF)
 
     ### declare the generator cost objective
     p_costs = gen_attrs['p_cost']

egret/models/tests/test_unit_commitment.py

@@ -169,9 +169,15 @@ def test_uc_runner():
 def test_uc_transmission_models():
 
     ## the network tests can optionally specify some kwargs so we can pass them into solve_unit_commitment
-    tc_networks = {'btheta_power_flow': [dict()], 'ptdf_power_flow':[{'ptdf_options': {'lazy':False}}, dict()], 'power_balance_constraints':[dict()],}
+    tc_networks = {'btheta_power_flow': [dict()],
+                   'ptdf_power_flow':[
+                       {'ptdf_options': {'lazy':False}},
+                       dict(),
+                      ],
+                   'power_balance_constraints':[dict()],
+                   }
     no_network = 'copperplate_power_flow'
-    test_names = ['tiny_uc_tc'] + ['tiny_uc_tc_{}'.format(i) for i in range(2,7+1)]
+    test_names = ['tiny_uc_tc'] + ['tiny_uc_tc_{}'.format(i) for i in range(2,11+1)]
     ## based on tiny_uc, tiny_uc_tc_2 has an interface, tiny_uc_tc_3 has a relaxed interface, tiny_uc_tc_4 has a relaxed flow limit; tiny_uc_tc_7 has a HVDC line
 
     for test_name in test_names:
@@ -261,3 +267,4 @@ def test_scuc():
     md_baseline = ModelData.read(os.path.join(
                             current_dir, 'uc_test_instances', 'test_scuc_sparse_enforce_sol.json'))
     assert math.isclose( md_results.data['system']['total_cost'], md_baseline.data['system']['total_cost'], rel_tol=rel_tol)
+