diff --git a/egret/benchmark/benchmark_dcopf.py b/egret/benchmark/benchmark_dcopf.py new file mode 100644 index 00000000..9041b931 --- /dev/null +++ b/egret/benchmark/benchmark_dcopf.py @@ -0,0 +1,32 @@ +import sys + +from egret.models.dcopf import (solve_dcopf, + create_btheta_dcopf_model, + create_ptdf_dcopf_model, + ) +from egret.data.model_data import ModelData + +solver = "xpress_persistent" +if len(sys.argv) < 3: + print("Useage: python benchmark_ptdf.py matpower_file dcopf_type") + print("""dcopf_type is one of "btheta", "fullptdf", "virtualptdf" """) + sys.exit(0) + +matpower_file = sys.argv[1] +dcopf_type = sys.argv[2].lower().replace("-","").replace("_","") + +md = ModelData.read(matpower_file) + +if dcopf_type == "btheta": + mdo = solve_dcopf(md, solver, dcopf_model_generator=create_btheta_dcopf_model) +elif "ptdf" in dcopf_type: + kwargs = {} + if dcopf_type == "fullptdf": + kwargs["virtual_ptdf"] = False + elif dcopf_type == "virtualptdf": + kwargs["virtual_ptdf"] = True + else: + raise RuntimeError(f"Unrecognized dcopf_type {dcopf_type}") + mdo = solve_dcopf(md, solver, dcopf_model_generator=create_ptdf_dcopf_model, **kwargs) +else: + raise RuntimeError(f"Unrecognized dcopf_type {dcopf_type}") diff --git a/egret/common/lazy_ptdf_utils.py b/egret/common/lazy_ptdf_utils.py index 7dfe4d75..3ef25483 100644 --- a/egret/common/lazy_ptdf_utils.py +++ b/egret/common/lazy_ptdf_utils.py @@ -151,6 +151,8 @@ def __init__(self, max_viol_add, md, prepend_str, time=None): self.violations_store = {} self.total_violations = 0 self.monitored_violations = 0 + self.min_flow_screen = True + self.min_violation = 0. def get_violations_named(self, name): for key in self.violations_store: @@ -158,8 +160,8 @@ def get_violations_named(self, name): yield key[1] def min_flow_violation(self): - if self.violations_store: - return min(self.violations_store.values()) + if len(self.violations_store) == self.max_viol_add: + return self.min_violation else: return 0. @@ -167,6 +169,9 @@ def _min_violation_key(self): d = self.violations_store return min(d, key=d.get) + def update_min_flow_violation(self): + self.min_violation = min(self.violations_store.values()) + def _add_violation(self, name, other_name, index, val): if other_name: key = ( name, (other_name, index) ) @@ -184,12 +189,13 @@ def _add_violation(self, name, other_name, index, val): f"violations_store: {self.violations_store}") del self.violations_store[min_key] + self.update_min_flow_violation() def _add_violations( self, name, other_name, viol_array, viol_indices): while viol_indices: idx = np.argmax(viol_array[viol_indices]) val = viol_array[viol_indices[idx]] - if val < self.min_flow_violation() and len(self.violations_store) >= self.max_viol_add: + if len(self.violations_store) >= self.max_viol_add and val < self.min_flow_violation(): break # If this violation is close in value to # one already in the set, it is likely @@ -201,17 +207,20 @@ def _add_violations( self, name, other_name, viol_array, viol_indices): # 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) - continue + # NOTE: the above assumes that this line will not be violated + # Since they often are (e.g. contingencies), we leave this + # logic out + #close_to_existing = any( math.isclose( val, existing ) for existing in self.violations_store.values() ) + #if close_to_existing: + # viol_indices.pop(idx) + # continue self._add_violation( name, other_name, viol_indices[idx], val ) viol_indices.pop(idx) def check_and_add_violations(self, name, flow_array, flow_variable, upper_lazy_limits, upper_enforced_limits, lower_lazy_limits, lower_enforced_limits, - monitored_indices, index_names, outer_name=None, PFV=None): + monitored_indices, index_names, outer_name=None, PFV=None, get_counts=False): if outer_name: # contingencies are named by cn, branch_idx, reduce to @@ -224,21 +233,28 @@ def check_and_add_violations(self, name, flow_array, flow_variable, ## get the indices of the violation ## here filter by least violation in violations_store ## in the limit, this will become 0 eventually -- - upper_viol_lazy_idx = np.nonzero(upper_viol_lazy_array > self.min_flow_violation())[0] + upper_viol_lazy_idx = np.nonzero(upper_viol_lazy_array > 0.0)[0] + + if len(upper_viol_lazy_idx) > 0: + upper_viol_array = flow_array[upper_viol_lazy_idx] - upper_enforced_limits[upper_viol_lazy_idx] + upper_viol_idx = self._calculate_total_and_monitored_violations(upper_viol_array, upper_viol_lazy_idx, monitored_indices, + flow_variable, flow_array, index_names, upper_enforced_limits, + name, outer_name, PFV) - upper_viol_array = flow_array[upper_viol_lazy_idx] - upper_enforced_limits[upper_viol_lazy_idx] - self._calculate_total_and_monitored_violations(upper_viol_array, upper_viol_lazy_idx, monitored_indices, - flow_variable, flow_array, index_names, upper_enforced_limits, - name, outer_name, PFV) + ## viol_lazy_idx will hold the lines we're adding + ## this iteration -- don't want to add lines + ## that are already in the monitored set - ## viol_lazy_idx will hold the lines we're adding - ## this iteration -- don't want to add lines - ## that are already in the monitored set + # eliminate lines in the monitored set + upper_viol_lazy_idx = set(upper_viol_lazy_idx).difference(monitored_indices) + upper_viol_lazy_idx = np.fromiter(upper_viol_lazy_idx, dtype='int64', count=len(upper_viol_lazy_idx)) - # eliminate lines in the monitored set - upper_viol_lazy_idx = list(set(upper_viol_lazy_idx).difference(monitored_indices)) + upper_viol_lazy_idx_idx = np.nonzero(upper_viol_lazy_array[upper_viol_lazy_idx] > self.min_flow_violation())[0] + upper_viol_lazy_idx = list(upper_viol_lazy_idx[upper_viol_lazy_idx_idx]) - self._add_violations( name, outer_name, upper_viol_lazy_array, upper_viol_lazy_idx ) + self._add_violations( name, outer_name, upper_viol_lazy_array, upper_viol_lazy_idx ) + else: + upper_viol_idx = [] ## check lower bound lower_viol_lazy_array = lower_lazy_limits - flow_array @@ -246,21 +262,34 @@ def check_and_add_violations(self, name, flow_array, flow_variable, ## get the indices of the violation ## here filter by least violation in violations_store ## in the limit, this will become 0 eventually -- - lower_viol_lazy_idx = np.nonzero(lower_viol_lazy_array > self.min_flow_violation())[0] + lower_viol_lazy_idx = np.nonzero(lower_viol_lazy_array > 0.0)[0] - lower_viol_array = lower_enforced_limits[lower_viol_lazy_idx] - flow_array[lower_viol_lazy_idx] - self._calculate_total_and_monitored_violations(lower_viol_array, lower_viol_lazy_idx, monitored_indices, - flow_variable, flow_array, index_names, lower_enforced_limits, - name, outer_name, PFV) + if len(lower_viol_lazy_idx) > 0: + lower_viol_array = lower_enforced_limits[lower_viol_lazy_idx] - flow_array[lower_viol_lazy_idx] + lower_viol_idx = self._calculate_total_and_monitored_violations(lower_viol_array, lower_viol_lazy_idx, monitored_indices, + flow_variable, flow_array, index_names, lower_enforced_limits, + name, outer_name, PFV) - ## viol_lazy_idx will hold the lines we're adding - ## this iteration -- don't want to add lines - ## that are already in the monitored set + ## viol_lazy_idx will hold the lines we're adding + ## this iteration -- don't want to add lines + ## that are already in the monitored set - # eliminate lines in the monitored set - lower_viol_lazy_idx = list(set(lower_viol_lazy_idx).difference(monitored_indices)) + # eliminate lines in the monitored set + lower_viol_lazy_idx = set(lower_viol_lazy_idx).difference(monitored_indices) + lower_viol_lazy_idx = np.fromiter(lower_viol_lazy_idx, dtype='int64', count=len(lower_viol_lazy_idx)) + + lower_viol_lazy_idx_idx = np.nonzero(lower_viol_lazy_array[lower_viol_lazy_idx] > self.min_flow_violation())[0] + lower_viol_lazy_idx = list(lower_viol_lazy_idx[lower_viol_lazy_idx_idx]) + + self._add_violations( name, outer_name, lower_viol_lazy_array, lower_viol_lazy_idx ) + else: + lower_viol_idx = [] - self._add_violations( name, outer_name, lower_viol_lazy_array, lower_viol_lazy_idx ) + if get_counts: + viol_idx = set(upper_viol_idx) + viol_idx.update(lower_viol_idx) + return len(viol_idx.difference(monitored_indices)) + return def _calculate_total_and_monitored_violations(self, viol_array, viol_lazy_idx, monitored_indices, @@ -313,6 +342,8 @@ def _calculate_total_and_monitored_violations(self, viol_array, viol_lazy_idx, m print(f'model: {pyo.value(flow_variable[element_name])}') print('') + return viol_idx + ## to hold the indicies of the violations ## in the model or block @@ -359,7 +390,7 @@ def check_violations(mb, md, PTDF, max_viol_add, time=None, prepend_str=""): mb._interfaces_monitored, PTDF.interface_keys) if PTDF.contingencies and \ - violations_store.total_violations == 0: + violations_store.total_violations < violations_store.monitored_violations + violations_store.max_viol_add: ## NOTE: checking contingency constraints in general could be very expensive ## we probably want to delay doing so until we have a nearly transmission feasible ## solution @@ -378,18 +409,40 @@ def check_violations(mb, md, PTDF, max_viol_add, time=None, prepend_str=""): ## In this way, we avoid (number of contingenies) adds PFV+PFV_delta_c - logger.debug("Checking contingency flows...") + logger.info("Checking contingency flows...") lazy_contingency_limits_upper = PTDF.lazy_contingency_limits - PFV lazy_contingency_limits_lower = -PTDF.lazy_contingency_limits - PFV enforced_contingency_limits_upper = PTDF.enforced_contingency_limits - PFV enforced_contingency_limits_lower = -PTDF.enforced_contingency_limits - PFV - for cn in PTDF.contingency_compensators: + if time in PTDF.contingency_compensators._order: + cn_iterator = PTDF.contingency_compensators._order[time] + else: + cn_iterator = PTDF.contingency_compensators + total = 0 + for c_checked, cn in enumerate(cn_iterator): PFV_delta = PTDF.calculate_masked_PFV_delta(cn, PFV, VA) - violations_store.check_and_add_violations('contingency', PFV_delta, mb.pfc, - lazy_contingency_limits_upper, enforced_contingency_limits_upper, - lazy_contingency_limits_lower, enforced_contingency_limits_lower, - mb._contingencies_monitored, PTDF.branches_keys_masked, - outer_name = cn, PFV = PFV) + PTDF.contingency_compensators._count[time,cn] = violations_store.check_and_add_violations('contingency', PFV_delta, mb.pfc, + lazy_contingency_limits_upper, enforced_contingency_limits_upper, + lazy_contingency_limits_lower, enforced_contingency_limits_lower, + mb._contingencies_monitored, PTDF.branches_keys_masked, + outer_name = cn, PFV = PFV, get_counts=True) + + total += PTDF.contingency_compensators._count[time,cn] + if time in PTDF.contingency_compensators._order: + # give us some choice in which violations we add + if total >= violations_store.max_viol_add: + c_checked += 1 + break + if (c_checked+1)%1000 == 0: + logger.info(f"Checked {c_checked+1} contingency flows...") + + PTDF.contingency_compensators._order[time] = list(sorted(PTDF.contingency_compensators, key=lambda cn: PTDF.contingency_compensators._count[time,cn], reverse=True)) + if len(PTDF.contingency_compensators) == c_checked+1: + logger.info(f"Checked **ALL** {c_checked} contingencies") + else: + logger.info(f"Checked {c_checked} contingencies") + else: + logger.info("Checked no contingencies") logger.debug(f"branches_monitored: {mb._idx_monitored}\n" f"interfaces_monitored: {mb._interfaces_monitored}\n" @@ -404,6 +457,7 @@ def check_violations(mb, md, PTDF, max_viol_add, time=None, prepend_str=""): return flows, violations_store.total_violations, violations_store.monitored_violations, viol_lazy + def _generate_flow_monitor_remove_message(flow_type, bn, slack, baseMVA, time): ret_str = "removing {0} {1} from monitored set".format(flow_type, bn) if time is not None: diff --git a/egret/common/solver_interface.py b/egret/common/solver_interface.py index ac83f5a7..a138b60d 100644 --- a/egret/common/solver_interface.py +++ b/egret/common/solver_interface.py @@ -82,7 +82,8 @@ def _solve_model(model, solve_method_options = None, return_solver = False, vars_to_load = None, - set_instance = True): + set_instance = True, + ): ''' Create and solve an Egret power system optimization model @@ -118,7 +119,6 @@ def _solve_model(model, ------- ''' - results = None ## termination conditions which are acceptable @@ -169,6 +169,7 @@ def _solve_model(model, else: model.solutions.load_from(results) + if return_solver: return model, results, solver return model, results diff --git a/egret/data/ptdf_utils.py b/egret/data/ptdf_utils.py index 2684bdc4..e9a8cb9f 100644 --- a/egret/data/ptdf_utils.py +++ b/egret/data/ptdf_utils.py @@ -51,6 +51,10 @@ def get_interface_ptdf_iterator(self, interface_name): def calculate_PFV(self,mb): pass + def _insert_reference_bus(self, bus_array, val): + return np.insert(bus_array, self._busname_to_index_map[self._reference_bus], val) + + class VirtualPTDFMatrix(_PTDFManagerBase): ''' Helper class which *looks like* a PTDF matrix, but doesn't @@ -142,6 +146,9 @@ def __init__(self, branches, buses, reference_bus, base_point, # dense array write buffer self._bus_sensi_buffer = np.empty((1,len(self.buses_keys_no_ref)), dtype=np.float64) + # helper + self._first_time_contingency = True + def _calculate_ptdf_factorization(self): logger.info("Calculating PTDF Matrix Factorization") MLU, B_dA, ref_bus_mask, contingency_compensators, B_dA_I, I = \ @@ -408,9 +415,6 @@ def get_interface_ptdf_iterator(self, interface_name): ''' yield from zip(self.buses_keys_no_ref, self._get_interface_row(interface_name)) - def _insert_reference_bus(self, bus_array, val): - return np.insert(bus_array, self._busname_to_index_map[self._reference_bus], val) - def _calculate_PFV_delta(self, cn, PFV, VA, masked): comp = self.contingency_compensators[cn] @@ -666,6 +670,8 @@ def __init__(self, branches, buses, reference_bus, base_point, if interfaces is None: interfaces = dict() + # no contingency analysis for Full PTDF matrix + self.contingencies = dict() self._calculate_ptdf_interface(interfaces) self._set_lazy_limits(ptdf_options) @@ -874,7 +880,7 @@ def get_interface_ptdf_iterator(self, interface_name): def bus_iterator(self): yield from self.buses_keys - def calculate_PFV(self,mb): + def calculate_masked_PFV(self,mb): NWV = np.fromiter((value(mb.p_nw[b]) for b in self.buses_keys), float, count=len(self.buses_keys)) NWV += self.phi_adjust_array @@ -884,7 +890,40 @@ def calculate_PFV(self,mb): PFV_I = self.PTDFM_I.dot(NWV) PFV_I += self.PTDFM_I_const - return PFV, PFV_I + return PFV, PFV_I, None + + def calculate_PFV(self,mb): + NWV = np.fromiter((value(mb.p_nw[b]) for b in self.buses_keys), float, count=len(self.buses_keys)) + NWV += self.phi_adjust_array + + PFV = self.PTDFM.dot(NWV) + PFV += self.phase_shift_array + + PFV_I = self.PTDFM_I.dot(NWV) + PFV_I += self.PTDFM_I_const + + return PFV, PFV_I, None + + def calculate_LMP(self, mb, dual, bus_balance_constr): + ## NOTE: unmonitored lines cannot contribute to LMPC + PFD = np.fromiter( ( value(dual[mb.ineq_pf_branch_thermal_bounds[bn]]) + if bn in mb.ineq_pf_branch_thermal_bounds else + 0. for i,bn in enumerate(self.branches_keys_masked) ), + float, count=len(self.branches_keys_masked)) + + ## interface constributes to LMP + PFID = np.fromiter( ( value(dual[mb.ineq_pf_interface_bounds[i_n]]) + if i_n in mb.ineq_pf_interface_bounds else + 0. for i,i_n in enumerate(self.interface_keys) ), + float, count=len(self.interface_keys)) + LMPC = -self.PTDFM.T.dot(PFD) + LMPI = -self.PTDFM_I.T.dot(PFID) + + LMPE = value(dual[bus_balance_constr]) + + LMP = LMPE + LMPC + LMPI + + return self._insert_reference_bus(LMP, LMPE) class PTDFLossesMatrix(PTDFMatrix): diff --git a/egret/model_library/transmission/tx_calc.py b/egret/model_library/transmission/tx_calc.py index 88418a1e..8155206a 100644 --- a/egret/model_library/transmission/tx_calc.py +++ b/egret/model_library/transmission/tx_calc.py @@ -766,6 +766,8 @@ def __init__(self, compensators, L, U, Pr, Pc): self._U = U self._Pr = Pr self._Pc = Pc + self._order = {} + self._count = {} def __getitem__(self, key): return self._compensators[key] diff --git a/egret/models/dcopf.py b/egret/models/dcopf.py index 4dd6036c..dbea5f2a 100644 --- a/egret/models/dcopf.py +++ b/egret/models/dcopf.py @@ -12,6 +12,7 @@ #TODO: document this with examples """ +from pyomo.common.timing import HierarchicalTimer import pyomo.environ as pe import numpy as np import egret.model_library.transmission.tx_utils as tx_utils @@ -50,14 +51,15 @@ def _include_feasibility_slack(model, bus_names, bus_p_loads, gens_by_bus, gen_a return p_rhs_kwargs, penalty_expr def create_btheta_dcopf_model(model_data, include_angle_diff_limits=False, include_feasibility_slack=False, pw_cost_model='delta', - keep_vars_for_out_of_service_elements=False): + keep_vars_for_out_of_service_elements=False, timer=None): + timer.start("data preprocessing") if keep_vars_for_out_of_service_elements: out_of_service_gens = tx_utils._get_out_of_service_gens(model_data) out_of_service_branches = tx_utils._get_out_of_service_branches(model_data) else: out_of_service_gens = list() out_of_service_branches = list() - + md = model_data.clone_in_service() tx_utils.scale_ModelData_to_pu(md, inplace = True) @@ -76,7 +78,9 @@ def create_btheta_dcopf_model(model_data, include_angle_diff_limits=False, inclu inlet_branches_by_bus, outlet_branches_by_bus = \ tx_utils.inlet_outlet_branches_by_bus(branches, buses) gens_by_bus = tx_utils.gens_by_bus(buses, gens) + timer.stop("data preprocessing") + timer.start("Pyomo model construction") model = pe.ConcreteModel() ### declare (and fix) the loads at the buses @@ -254,10 +258,13 @@ def create_btheta_dcopf_model(model_data, include_angle_diff_limits=False, inclu model_data.data['elements']['branch'][branch_name]['in_service'] = False md.data['elements']['branch'][branch_name]['in_service'] = False + timer.stop("Pyomo model construction") + return model, md -def create_ptdf_dcopf_model(model_data, include_feasibility_slack=False, base_point=BasePointType.FLATSTART, ptdf_options=None, pw_cost_model='delta'): +def create_ptdf_dcopf_model(model_data, include_feasibility_slack=False, base_point=BasePointType.FLATSTART, ptdf_options=None, pw_cost_model='delta', virtual_ptdf=True, timer=None): + timer.start("data preprocessing") ptdf_options = lpu.populate_default_ptdf_options(ptdf_options) baseMVA = model_data.data['system']['baseMVA'] @@ -283,7 +290,9 @@ def create_ptdf_dcopf_model(model_data, include_feasibility_slack=False, base_po inlet_branches_by_bus, outlet_branches_by_bus = \ tx_utils.inlet_outlet_branches_by_bus(branches, buses) gens_by_bus = tx_utils.gens_by_bus(buses, gens) + timer.stop("data preprocessing") + timer.start("Pyomo model construction") model = pe.ConcreteModel() ### declare (and fix) the loads at the buses @@ -353,7 +362,12 @@ def create_ptdf_dcopf_model(model_data, include_feasibility_slack=False, base_po ## Do and store PTDF calculation reference_bus = md.data['system']['reference_bus'] - PTDF = ptdf_utils.VirtualPTDFMatrix(branches, buses, reference_bus, base_point, ptdf_options, branches_keys=branches_idx, buses_keys=buses_idx) + timer.start("PTDF Pre-computation") + if virtual_ptdf: + PTDF = ptdf_utils.VirtualPTDFMatrix(branches, buses, reference_bus, base_point, ptdf_options, branches_keys=branches_idx, buses_keys=buses_idx) + else: + PTDF = ptdf_utils.PTDFMatrix(branches, buses, reference_bus, base_point, ptdf_options, branches_keys=branches_idx, buses_keys=buses_idx) + timer.stop("PTDF Pre-computation") model._PTDF = PTDF model._ptdf_options = ptdf_options @@ -410,11 +424,11 @@ def create_ptdf_dcopf_model(model_data, include_feasibility_slack=False, base_po obj_expr += penalty_expr model.obj = pe.Objective(expr=obj_expr) - + timer.stop("Pyomo model construction") return model, md -def _lazy_ptdf_dcopf_model_solve_loop(m, md, solver, solver_tee=True, symbolic_solver_labels=False, iteration_limit=100000): +def _lazy_ptdf_dcopf_model_solve_loop(m, md, solver, solver_tee=True, symbolic_solver_labels=False, iteration_limit=100000, timer=None): ''' The lazy PTDF DCOPF solver loop. This function iteratively adds violated transmission constraints until either the result is @@ -443,10 +457,12 @@ def _lazy_ptdf_dcopf_model_solve_loop(m, md, solver, solver_tee=True, symbolic_s int : The number of iterations before termination ''' + timer.start("Lazy PTDF Loop") from pyomo.solvers.plugins.solvers.persistent_solver import PersistentSolver from egret.common.solver_interface import _solve_model PTDF = m._PTDF + PTDF.timer = timer ptdf_options = m._ptdf_options @@ -454,8 +470,10 @@ def _lazy_ptdf_dcopf_model_solve_loop(m, md, solver, solver_tee=True, symbolic_s for i in range(iteration_limit): + timer.start("Violation Check") flows, viol_num, mon_viol_num, viol_lazy \ = lpu.check_violations(m, md, PTDF, ptdf_options['max_violations_per_iteration']) + timer.stop("Violation Check") iter_status_str = "iteration {0}, found {1} violation(s)".format(i,viol_num) if mon_viol_num: @@ -468,28 +486,33 @@ def _lazy_ptdf_dcopf_model_solve_loop(m, md, solver, solver_tee=True, symbolic_s ## load the duals now too, if we're using a persistent solver if persistent_solver: solver.load_duals() + timer.stop("Lazy PTDF Loop") return lpu.LazyPTDFTerminationCondition.NORMAL elif viol_num == mon_viol_num: logger.warning('WARNING: Terminating with monitored violations! Result is not transmission feasible.') if persistent_solver: solver.load_duals() + timer.stop("Lazy PTDF Loop") return lpu.LazyPTDFTerminationCondition.FLOW_VIOLATION + timer.start("Add Constraints") lpu.add_violations(viol_lazy, flows, m, md, solver, ptdf_options, PTDF) total_flow_constr_added = len(viol_lazy) logger.info( "iteration {0}, added {1} flow constraint(s)".format(i,total_flow_constr_added)) + timer.stop("Add Constraints") if persistent_solver: - m, results, solver = _solve_model(m, solver, solver_tee=solver_tee, return_solver=True, vars_to_load=[], set_instance=False) + m, results, solver = _solve_model(m, solver, solver_tee=solver_tee, return_solver=True, vars_to_load=[], set_instance=False, timer=timer) solver.load_vars() else: - m, results, solver = _solve_model(m, solver, solver_tee=solver_tee, return_solver=True) + m, results, solver = _solve_model(m, solver, solver_tee=solver_tee, return_solver=True, timer=timer) else: # we hit the iteration limit logger.warning('WARNING: Exiting on maximum iterations for lazy PTDF model. Result is not transmission feasible.') if persistent_solver: solver.load_duals() + timer.stop("Lazy PTDF Loop") return lpu.LazyPTDFTerminationCondition.ITERATION_LIMIT @@ -531,25 +554,31 @@ def solve_dcopf(model_data, kwargs : dictionary (optional) Additional arguments for building model ''' + timer = HierarchicalTimer() + timer.start("Solve DCOPF") + kwargs['timer'] = timer - import pyomo.environ as pe - import pyomo.opt as po from pyomo.environ import value from egret.common.solver_interface import _solve_model from egret.model_library.transmission.tx_utils import \ scale_ModelData_to_pu, unscale_ModelData_to_pu + timer.start("Model Generation") m, md = dcopf_model_generator(model_data, **kwargs) + timer.stop("Model Generation") m.dual = pe.Suffix(direction=pe.Suffix.IMPORT) + timer.start("Model Solution") m, results, solver = _solve_model(m,solver,timelimit=timelimit,solver_tee=solver_tee, - symbolic_solver_labels=symbolic_solver_labels,solver_options=options, return_solver=True) + symbolic_solver_labels=symbolic_solver_labels,solver_options=options, return_solver=True, timer=timer) if dcopf_model_generator == create_ptdf_dcopf_model and m._ptdf_options['lazy']: iter_limit = m._ptdf_options['iteration_limit'] - term_cond = _lazy_ptdf_dcopf_model_solve_loop(m, md, solver, solver_tee=solver_tee, symbolic_solver_labels=symbolic_solver_labels,iteration_limit=iter_limit) + term_cond = _lazy_ptdf_dcopf_model_solve_loop(m, md, solver, solver_tee=solver_tee, symbolic_solver_labels=symbolic_solver_labels,iteration_limit=iter_limit, timer=timer) + timer.stop("Model Solution") + timer.start("Model Post-Processing") # save results data to ModelData object gens = dict(md.elements(element_type='generator')) buses = dict(md.elements(element_type='bus')) @@ -564,6 +593,7 @@ def solve_dcopf(model_data, ## calculate the power flows from our PTDF matrix for maximum precision ## calculate the LMPC (LMP congestion) using numpy + timer.start("Flow calculation") if dcopf_model_generator == create_ptdf_dcopf_model: PTDF = m._PTDF @@ -575,7 +605,9 @@ def solve_dcopf(model_data, else: for k, k_dict in branches.items(): k_dict['pf'] = value(m.pf[k]) + timer.stop("Flow calculation") + timer.start("LMP calculation") if dcopf_model_generator == create_ptdf_dcopf_model: if hasattr(m, 'p_load_shed'): md.data['system']['p_balance_violation'] = value(m.p_load_shed) - value(m.p_over_generation) @@ -585,7 +617,9 @@ def solve_dcopf(model_data, b_dict = buses[b] b_dict['lmp'] = LMP[i] b_dict['pl'] = value(m.pl[b]) - b_dict['va'] = degrees(VA[i]) + if VA is not None: + for i,b in enumerate(buses_idx): + b_dict['va'] = degrees(VA[i]) else: for b,b_dict in buses.items(): if hasattr(m, 'p_load_shed'): @@ -596,19 +630,25 @@ def solve_dcopf(model_data, b_dict['va'] = degrees(value(m.va[b])) else: raise Exception("Unrecognized dcopf_model_generator {}".format(dcopf_model_generator)) + timer.stop("LMP calculation") for k, k_dict in dc_branches.items(): k_dict['pf'] = value(m.dcpf[k]) + timer.start("Unscaling PU") unscale_ModelData_to_pu(md, inplace=True) + timer.stop("Unscaling PU") + + timer.stop("Model Post-Processing") + timer.stop("Solve DCOPF") if return_model and return_results: - return md, m, results + return md, m, results, timer elif return_model: - return md, m + return md, m, timer elif return_results: - return md, results - return md + return md, results, timer + return md, timer def check_instance_feasibility(instance, tolerance, active_only=True): infeasibilities = list() diff --git a/egret/models/scopf.py b/egret/models/scopf.py index 0d282a8d..00437a76 100644 --- a/egret/models/scopf.py +++ b/egret/models/scopf.py @@ -12,6 +12,7 @@ #TODO: document this with examples """ +from pyomo.common.timing import HierarchicalTimer import pyomo.environ as pyo import numpy as np import egret.model_library.transmission.tx_utils as tx_utils @@ -32,8 +33,9 @@ def create_scopf_model(model_data, include_feasibility_slack=False, - base_point=BasePointType.FLATSTART, ptdf_options=None, pw_cost_model='delta'): + base_point=BasePointType.FLATSTART, ptdf_options=None, pw_cost_model='delta', timer=None): + timer.start("data preprocessing") ptdf_options = lpu.populate_default_ptdf_options(ptdf_options) baseMVA = model_data.data['system']['baseMVA'] @@ -60,6 +62,9 @@ def create_scopf_model(model_data, include_feasibility_slack=False, inlet_branches_by_bus, outlet_branches_by_bus = \ tx_utils.inlet_outlet_branches_by_bus(branches, buses) gens_by_bus = tx_utils.gens_by_bus(buses, gens) + timer.stop("data preprocessing") + + timer.start("Pyomo model construction") model = pyo.ConcreteModel() @@ -134,12 +139,16 @@ def create_scopf_model(model_data, include_feasibility_slack=False, ### as we find violations model._contingency_set = pyo.Set(within=model._contingencies*model._branches) model.pfc = pyo.Expression(model._contingency_set) + model.pfc_slack_pos = pyo.Var(model._contingency_set, domain=pyo.NonNegativeReals, dense=False) + model.pfc_slack_neg = pyo.Var(model._contingency_set, domain=pyo.NonNegativeReals, dense=False) ## Do and store PTDF calculation reference_bus = md.data['system']['reference_bus'] + timer.start("PTDF Pre-computation") PTDF = ptdf_utils.VirtualPTDFMatrix(branches, buses, reference_bus, base_point, ptdf_options,\ contingencies=contingencies, branches_keys=branches_idx, buses_keys=buses_idx) + timer.stop("PTDF Pre-computation") model._PTDF = PTDF model._ptdf_options = ptdf_options @@ -184,8 +193,13 @@ def create_scopf_model(model_data, include_feasibility_slack=False, if include_feasibility_slack: obj_expr += penalty_expr - model.obj = pyo.Objective(expr=obj_expr) + model.ContingencyViolationCost = pyo.Expression(expr=0.) + model.TimePeriodLengthHours = pyo.Param(initialize=1) + model.SystemContingencyLimitPenalty = pyo.Param(initialize=500.0*baseMVA) + obj_expr += model.ContingencyViolationCost + model.obj = pyo.Objective(expr=obj_expr) + timer.stop("Pyomo model construction") return model, md @@ -228,6 +242,9 @@ def solve_scopf(model_data, kwargs : dictionary (optional) Additional arguments for building model ''' + timer = HierarchicalTimer() + timer.start("Solve SCOPF") + kwargs['timer'] = timer import pyomo.environ as pe import pyomo.opt as po @@ -236,17 +253,22 @@ def solve_scopf(model_data, from egret.model_library.transmission.tx_utils import \ scale_ModelData_to_pu, unscale_ModelData_to_pu + timer.start("Model Generation") m, md = scopf_model_generator(model_data, **kwargs) + timer.stop("Model Generation") m.dual = pyo.Suffix(direction=pyo.Suffix.IMPORT) + timer.start("Model Solution") m, results, solver = _solve_model(m,solver,timelimit=timelimit,solver_tee=solver_tee, - symbolic_solver_labels=symbolic_solver_labels,solver_options=options, return_solver=True) + symbolic_solver_labels=symbolic_solver_labels,solver_options=options, return_solver=True, timer=timer) if m._ptdf_options['lazy']: iter_limit = m._ptdf_options['iteration_limit'] - term_cond = _lazy_ptdf_dcopf_model_solve_loop(m, md, solver, solver_tee=solver_tee, symbolic_solver_labels=symbolic_solver_labels,iteration_limit=iter_limit) + term_cond = _lazy_ptdf_dcopf_model_solve_loop(m, md, solver, solver_tee=solver_tee, symbolic_solver_labels=symbolic_solver_labels,iteration_limit=iter_limit, timer=timer) + timer.stop("Model Solution") + timer.start("Model Post-Processing") # save results data to ModelData object gens = dict(md.elements(element_type='generator')) buses = dict(md.elements(element_type='bus')) @@ -261,6 +283,7 @@ def solve_scopf(model_data, ## calculate the power flows from our PTDF matrix for maximum precision ## calculate the LMPC (LMP congestion) using numpy + timer.start("Flow calculation") PTDF = m._PTDF PFV, _, VA = PTDF.calculate_PFV(m) @@ -268,20 +291,24 @@ def solve_scopf(model_data, branches_idx = PTDF.branches_keys for i,bn in enumerate(branches_idx): branches[bn]['pf'] = PFV[i] + timer.stop("Flow calculation") if hasattr(m, 'p_load_shed'): md.data['system']['p_balance_violation'] = value(m.p_load_shed) - value(m.p_over_generation) buses_idx = PTDF.buses_keys + timer.start("LMP calculation") LMP = PTDF.calculate_LMP(m, m.dual, m.eq_p_balance) for i,b in enumerate(buses_idx): b_dict = buses[b] b_dict['lmp'] = LMP[i] b_dict['pl'] = value(m.pl[b]) b_dict['va'] = degrees(VA[i]) + timer.stop("LMP calculation") for k, k_dict in dc_branches.items(): k_dict['pf'] = value(m.dcpf[k]) + timer.start("Monitored Contingency Flow Computation") contingencies = dict(md.elements(element_type='contingency')) contingency_flows = PTDF.calculate_monitored_contingency_flows(m) for (cn,bn), flow in contingency_flows.items(): @@ -289,8 +316,14 @@ def solve_scopf(model_data, if 'monitored_branches' not in c_dict: c_dict['monitored_branches'] = {} c_dict['monitored_branches'][bn] = {'pf': flow} + timer.stop("Monitored Contingency Flow Computation") + timer.start("Unscaling PU") unscale_ModelData_to_pu(md, inplace=True) + timer.stop("Unscaling PU") + + timer.stop("Model Post-Processing") + timer.stop("Solve SCOPF") if return_model and return_results: return md, m, results @@ -298,4 +331,4 @@ def solve_scopf(model_data, return md, m elif return_results: return md, results - return md + return md, timer diff --git a/egret/parsers/prescient_dat_parser.py b/egret/parsers/prescient_dat_parser.py index 299b8604..f9c97d54 100644 --- a/egret/parsers/prescient_dat_parser.py +++ b/egret/parsers/prescient_dat_parser.py @@ -94,8 +94,8 @@ def create_model_data_dict_params(params, keep_names=False): 'to_bus' : params.BusTo[l], 'reactance' : params.Impedence[l], 'rating_long_term' : params.ThermalLimit[l], - 'rating_short_term' : params.ThermalLimit[l], - 'rating_emergency' : params.ThermalLimit[l], + 'rating_short_term' : 1.15*params.ThermalLimit[l], + 'rating_emergency' : 1.3*params.ThermalLimit[l], 'in_service' : True, 'branch_type' : 'line', 'angle_diff_min': -90,