From 9e4dcc5696603ea955c841cc3d494e4862b28249 Mon Sep 17 00:00:00 2001 From: Bernard Knueven Date: Wed, 27 Jul 2022 12:59:44 -0600 Subject: [PATCH 01/11] setting up dcopf benchmark --- egret/benchmark/benchmark_dcopf.py | 32 +++++++++++++++ egret/common/solver_interface.py | 15 ++++++- egret/data/ptdf_utils.py | 46 ++++++++++++++++++--- egret/models/dcopf.py | 66 ++++++++++++++++++++++++------ 4 files changed, 139 insertions(+), 20 deletions(-) create mode 100644 egret/benchmark/benchmark_dcopf.py 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/solver_interface.py b/egret/common/solver_interface.py index ac83f5a7..f247a40a 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, + timer = None): ''' Create and solve an Egret power system optimization model @@ -118,7 +119,7 @@ def _solve_model(model, ------- ''' - + timer.start("model solve") results = None ## termination conditions which are acceptable @@ -149,16 +150,23 @@ def _solve_model(model, if isinstance(solver, PersistentSolver): if set_instance: + timer.start("set instance") solver.set_instance(model, symbolic_solver_labels=symbolic_solver_labels) + timer.stop("set instance") + timer.start("solve call") results = solver.solve(model, tee=solver_tee, load_solutions=False, save_results=False, **solve_method_options) + timer.stop("solve call") else: + timer.start("solve call") results = solver.solve(model, tee=solver_tee, \ symbolic_solver_labels=symbolic_solver_labels, load_solutions=False, **solve_method_options) + timer.stop("solve call") if results.solver.termination_condition not in safe_termination_conditions: raise Exception('Problem encountered during solve, termination_condition {}'.format(results.solver.termination_condition)) + timer.start("load solution") if isinstance(solver, PersistentSolver): solver.load_vars(vars_to_load) if vars_to_load is None: @@ -168,6 +176,9 @@ def _solve_model(model, solver.load_slacks() else: model.solutions.load_from(results) + timer.stop("load solution") + + timer.stop("model solve") if return_solver: return model, results, solver diff --git a/egret/data/ptdf_utils.py b/egret/data/ptdf_utils.py index 2684bdc4..2cf22bda 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 @@ -408,9 +412,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 +667,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 +877,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 +887,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/models/dcopf.py b/egret/models/dcopf.py index 4dd6036c..1a18b1e4 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,6 +457,7 @@ 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 @@ -454,8 +469,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 +485,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 +553,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 +592,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 +604,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 +616,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,11 +629,18 @@ 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") + print(timer) if return_model and return_results: return md, m, results From 981122bedc1cf115d2bb59f7fd8ffd5f314cd53d Mon Sep 17 00:00:00 2001 From: Bernard Knueven Date: Tue, 27 Sep 2022 09:41:28 -0600 Subject: [PATCH 02/11] returning timer --- egret/models/dcopf.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/egret/models/dcopf.py b/egret/models/dcopf.py index 1a18b1e4..a322472b 100644 --- a/egret/models/dcopf.py +++ b/egret/models/dcopf.py @@ -640,15 +640,14 @@ def solve_dcopf(model_data, timer.stop("Model Post-Processing") timer.stop("Solve DCOPF") - print(timer) 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() From 9ea1c02bdf42608e19f1ab603f67bd79c01e29a3 Mon Sep 17 00:00:00 2001 From: Bernard Knueven Date: Wed, 28 Sep 2022 12:01:09 -0600 Subject: [PATCH 03/11] tracking row time --- egret/data/ptdf_utils.py | 4 ++++ egret/models/dcopf.py | 1 + 2 files changed, 5 insertions(+) diff --git a/egret/data/ptdf_utils.py b/egret/data/ptdf_utils.py index 2cf22bda..e2b6e14d 100644 --- a/egret/data/ptdf_utils.py +++ b/egret/data/ptdf_utils.py @@ -302,6 +302,7 @@ def _set_lazy_limits(self, ptdf_options): self.enforced_interface_min_limits) def _get_ptdf_row(self, branch_name): + self.timer.start("Generate PTDF row") if branch_name in self._ptdf_rows: PTDF_row = self._ptdf_rows[branch_name] else: @@ -309,6 +310,7 @@ def _get_ptdf_row(self, branch_name): branch_idx = self._branchname_to_index_map[branch_name] PTDF_row = self.MLU.solve(self.B_dA[branch_idx].toarray(out=self._bus_sensi_buffer)[0], trans='T') self._ptdf_rows[branch_name] = PTDF_row + self.timer.stop("Generate PTDF row") return PTDF_row def _get_contingency_row(self, contingency_name, branch_name): @@ -834,7 +836,9 @@ def _set_lazy_limits(self, ptdf_options): def get_branch_ptdf_iterator(self, branch_name): row_idx = self._branchname_to_index_map[branch_name] ## get the row slice + self.timer.start("Generate PTDF row") PTDF_row = self.PTDFM[row_idx] + self.timer.stop("Generate PTDF row") yield from zip(self.buses_keys, PTDF_row) def get_branch_ptdf_abs_max(self, branch_name): diff --git a/egret/models/dcopf.py b/egret/models/dcopf.py index a322472b..dbea5f2a 100644 --- a/egret/models/dcopf.py +++ b/egret/models/dcopf.py @@ -462,6 +462,7 @@ def _lazy_ptdf_dcopf_model_solve_loop(m, md, solver, solver_tee=True, symbolic_s from egret.common.solver_interface import _solve_model PTDF = m._PTDF + PTDF.timer = timer ptdf_options = m._ptdf_options From 3d11c3dc27d2843226697d4bbcd0a0ed57130fca Mon Sep 17 00:00:00 2001 From: Bernard Knueven Date: Mon, 10 Oct 2022 12:07:42 -0600 Subject: [PATCH 04/11] benchmark SCOPF --- egret/models/scopf.py | 43 ++++++++++++++++++++++++++++++++++++++----- 1 file changed, 38 insertions(+), 5 deletions(-) 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 From 20df54dd9d629cc66b1a66617f06dcb3c0c03d82 Mon Sep 17 00:00:00 2001 From: Bernard Knueven Date: Mon, 10 Oct 2022 13:43:37 -0600 Subject: [PATCH 05/11] trying new violation checker --- egret/common/lazy_ptdf_utils.py | 21 ++++++++++++++++++++- egret/data/ptdf_utils.py | 3 +++ 2 files changed, 23 insertions(+), 1 deletion(-) diff --git a/egret/common/lazy_ptdf_utils.py b/egret/common/lazy_ptdf_utils.py index 7dfe4d75..487d4f3b 100644 --- a/egret/common/lazy_ptdf_utils.py +++ b/egret/common/lazy_ptdf_utils.py @@ -189,7 +189,7 @@ 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 @@ -383,6 +383,9 @@ def check_violations(mb, md, PTDF, max_viol_add, time=None, prepend_str=""): 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 + if PTDF._first_time_contingency: + _old_max_viol_add = violations_store.max_viol_add + violations_store.max_viol_add = 1e100 for cn in PTDF.contingency_compensators: PFV_delta = PTDF.calculate_masked_PFV_delta(cn, PFV, VA) violations_store.check_and_add_violations('contingency', PFV_delta, mb.pfc, @@ -390,6 +393,21 @@ def check_violations(mb, md, PTDF, max_viol_add, time=None, prepend_str=""): lazy_contingency_limits_lower, enforced_contingency_limits_lower, mb._contingencies_monitored, PTDF.branches_keys_masked, outer_name = cn, PFV = PFV) + if not PTDF._first_time_contingency: + if len(violations_store.violations_store) == violations_store.max_viol_add: + break + if PTDF._first_time_contingency: + # sort PTDF.contingency_compensators + sorter = {cn : 0. for cn in PTDF.contingency_compensators} + for key, val in violations_store.violations_store.items(): + if isinstance(key[1], tuple): + sorter[key[1][0]] += val + PTDF.contingency_compensators._compensators = dict(sorted(PTDF.contingency_compensators._compensators.items(), key=lambda x: sorter[x[0]], reverse=True)) + # take violations_store back down to max_viol_add + violations_store.violations_store = dict(sorted(violations_store.violations_store.items(), key=lambda x: x[1], reverse=True)[:_old_max_viol_add]) + violations_store.max_viol_add = _old_max_viol_add + + PTDF._first_time_contingency = False logger.debug(f"branches_monitored: {mb._idx_monitored}\n" f"interfaces_monitored: {mb._interfaces_monitored}\n" @@ -404,6 +422,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/data/ptdf_utils.py b/egret/data/ptdf_utils.py index e2b6e14d..f12769ea 100644 --- a/egret/data/ptdf_utils.py +++ b/egret/data/ptdf_utils.py @@ -146,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 = \ From 5b5be1fa90348fa9070ff7fcaf539416dc45a07d Mon Sep 17 00:00:00 2001 From: Bernard Knueven Date: Mon, 10 Oct 2022 15:28:13 -0600 Subject: [PATCH 06/11] updates to contingency violation checker --- egret/common/lazy_ptdf_utils.py | 48 ++++++++++----------- egret/model_library/transmission/tx_calc.py | 2 + 2 files changed, 25 insertions(+), 25 deletions(-) diff --git a/egret/common/lazy_ptdf_utils.py b/egret/common/lazy_ptdf_utils.py index 487d4f3b..a6dfae4a 100644 --- a/egret/common/lazy_ptdf_utils.py +++ b/egret/common/lazy_ptdf_utils.py @@ -151,6 +151,7 @@ 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 def get_violations_named(self, name): for key in self.violations_store: @@ -158,7 +159,7 @@ def get_violations_named(self, name): yield key[1] def min_flow_violation(self): - if self.violations_store: + if len(self.violations_store) == self.max_viol_add: return min(self.violations_store.values()) else: return 0. @@ -211,7 +212,7 @@ def _add_violations( self, name, other_name, viol_array, viol_indices): 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 @@ -259,9 +260,16 @@ def check_and_add_violations(self, name, flow_array, flow_variable, # eliminate lines in the monitored set lower_viol_lazy_idx = list(set(lower_viol_lazy_idx).difference(monitored_indices)) + count_lower = len(lower_viol_lazy_idx) self._add_violations( name, outer_name, lower_viol_lazy_array, lower_viol_lazy_idx ) + if get_counts: + viol_indices = set(np.nonzero(lower_viol_lazy_array > 0)[0]) + viol_indices.update(np.nonzero(upper_viol_lazy_array > 0)[0]) + return len(viol_indices.difference(monitored_indices)) + return + def _calculate_total_and_monitored_violations(self, viol_array, viol_lazy_idx, monitored_indices, flow_variable, flow_array, index_names, limits, @@ -359,7 +367,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: ## 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 @@ -383,31 +391,21 @@ def check_violations(mb, md, PTDF, max_viol_add, time=None, prepend_str=""): 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 - if PTDF._first_time_contingency: - _old_max_viol_add = violations_store.max_viol_add - violations_store.max_viol_add = 1e100 - 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 + for cn in 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) - if not PTDF._first_time_contingency: + 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) + if time in PTDF.contingency_compensators._order: if len(violations_store.violations_store) == violations_store.max_viol_add: break - if PTDF._first_time_contingency: - # sort PTDF.contingency_compensators - sorter = {cn : 0. for cn in PTDF.contingency_compensators} - for key, val in violations_store.violations_store.items(): - if isinstance(key[1], tuple): - sorter[key[1][0]] += val - PTDF.contingency_compensators._compensators = dict(sorted(PTDF.contingency_compensators._compensators.items(), key=lambda x: sorter[x[0]], reverse=True)) - # take violations_store back down to max_viol_add - violations_store.violations_store = dict(sorted(violations_store.violations_store.items(), key=lambda x: x[1], reverse=True)[:_old_max_viol_add]) - violations_store.max_viol_add = _old_max_viol_add - - PTDF._first_time_contingency = False + PTDF.contingency_compensators._order[time] = list(sorted(PTDF.contingency_compensators, key=lambda cn: PTDF.contingency_compensators._count[time,cn], reverse=True)) logger.debug(f"branches_monitored: {mb._idx_monitored}\n" f"interfaces_monitored: {mb._interfaces_monitored}\n" 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] From 3ae6ea540b6f65c4b64d7b94fad1f7b7bc83219e Mon Sep 17 00:00:00 2001 From: Bernard Knueven Date: Mon, 10 Oct 2022 15:54:22 -0600 Subject: [PATCH 07/11] debug --- egret/common/lazy_ptdf_utils.py | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/egret/common/lazy_ptdf_utils.py b/egret/common/lazy_ptdf_utils.py index a6dfae4a..9adfe273 100644 --- a/egret/common/lazy_ptdf_utils.py +++ b/egret/common/lazy_ptdf_utils.py @@ -386,11 +386,12 @@ 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...") + print("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 + all_contingencies = set(PTDF.contingency_compensators) if time in PTDF.contingency_compensators._order: cn_iterator = PTDF.contingency_compensators._order[time] else: @@ -402,10 +403,20 @@ def check_violations(mb, md, PTDF, max_viol_add, time=None, prepend_str=""): lazy_contingency_limits_lower, enforced_contingency_limits_lower, mb._contingencies_monitored, PTDF.branches_keys_masked, outer_name = cn, PFV = PFV, get_counts=True) + all_contingencies.remove(cn) if time in PTDF.contingency_compensators._order: if len(violations_store.violations_store) == violations_store.max_viol_add: + print(violations_store.violations_store) + print(violations_store.max_viol_add) + print("BREAKING") break PTDF.contingency_compensators._order[time] = list(sorted(PTDF.contingency_compensators, key=lambda cn: PTDF.contingency_compensators._count[time,cn], reverse=True)) + if len(all_contingencies) > 0: + print(f"ONLY checked a subset of contingencies! Number checked: {len(PTDF.contingency_compensators)-len(all_contingencies)}") + else: + print(f"CHECKED ALL CONTINGENCIES, Number checked: {len(PTDF.contingency_compensators) - len(all_contingencies)}") + else: + print("Skipping contingency flow check") logger.debug(f"branches_monitored: {mb._idx_monitored}\n" f"interfaces_monitored: {mb._interfaces_monitored}\n" From 8281ec032d0393cc9fb686108cadd838204491e9 Mon Sep 17 00:00:00 2001 From: Bernard Knueven Date: Mon, 10 Oct 2022 16:07:28 -0600 Subject: [PATCH 08/11] fixing TODO --- egret/common/lazy_ptdf_utils.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/egret/common/lazy_ptdf_utils.py b/egret/common/lazy_ptdf_utils.py index 9adfe273..1b87695b 100644 --- a/egret/common/lazy_ptdf_utils.py +++ b/egret/common/lazy_ptdf_utils.py @@ -265,8 +265,8 @@ def check_and_add_violations(self, name, flow_array, flow_variable, self._add_violations( name, outer_name, lower_viol_lazy_array, lower_viol_lazy_idx ) if get_counts: - viol_indices = set(np.nonzero(lower_viol_lazy_array > 0)[0]) - viol_indices.update(np.nonzero(upper_viol_lazy_array > 0)[0]) + viol_indices = set(np.nonzero((lower_enforced_limits - flow_array) > 0)[0]) + viol_indices.update(np.nonzero((flow_array - upper_enforced_limits) > 0)[0]) return len(viol_indices.difference(monitored_indices)) return @@ -396,6 +396,7 @@ def check_violations(mb, md, PTDF, max_viol_add, time=None, prepend_str=""): cn_iterator = PTDF.contingency_compensators._order[time] else: cn_iterator = PTDF.contingency_compensators + total = 0 for cn in cn_iterator: PFV_delta = PTDF.calculate_masked_PFV_delta(cn, PFV, VA) PTDF.contingency_compensators._count[time,cn] = violations_store.check_and_add_violations('contingency', PFV_delta, mb.pfc, @@ -403,9 +404,11 @@ def check_violations(mb, md, PTDF, max_viol_add, time=None, prepend_str=""): 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] all_contingencies.remove(cn) if time in PTDF.contingency_compensators._order: - if len(violations_store.violations_store) == violations_store.max_viol_add: + if total >= violations_store.max_viol_add: print(violations_store.violations_store) print(violations_store.max_viol_add) print("BREAKING") From bad1a41d95bf4de23e6fee9cd45cb99a23403c35 Mon Sep 17 00:00:00 2001 From: Bernard Knueven Date: Tue, 11 Oct 2022 16:50:35 -0600 Subject: [PATCH 09/11] cleanup --- egret/common/lazy_ptdf_utils.py | 109 +++++++++++++++++++------------- 1 file changed, 66 insertions(+), 43 deletions(-) diff --git a/egret/common/lazy_ptdf_utils.py b/egret/common/lazy_ptdf_utils.py index 1b87695b..3ef25483 100644 --- a/egret/common/lazy_ptdf_utils.py +++ b/egret/common/lazy_ptdf_utils.py @@ -152,6 +152,7 @@ def __init__(self, max_viol_add, md, prepend_str, time=None): 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: @@ -160,7 +161,7 @@ def get_violations_named(self, name): def min_flow_violation(self): if len(self.violations_store) == self.max_viol_add: - return min(self.violations_store.values()) + return self.min_violation else: return 0. @@ -168,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) ) @@ -185,6 +189,7 @@ 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: @@ -202,10 +207,13 @@ 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) @@ -225,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 @@ -247,27 +262,33 @@ 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)) - count_lower = len(lower_viol_lazy_idx) + # 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)) - self._add_violations( name, outer_name, lower_viol_lazy_array, 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 = [] if get_counts: - viol_indices = set(np.nonzero((lower_enforced_limits - flow_array) > 0)[0]) - viol_indices.update(np.nonzero((flow_array - upper_enforced_limits) > 0)[0]) - return len(viol_indices.difference(monitored_indices)) + viol_idx = set(upper_viol_idx) + viol_idx.update(lower_viol_idx) + return len(viol_idx.difference(monitored_indices)) return @@ -321,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 @@ -367,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 == violations_store.monitored_violations: + 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 @@ -386,18 +409,17 @@ 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 - print("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 - all_contingencies = set(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 cn in cn_iterator: + for c_checked, cn in enumerate(cn_iterator): PFV_delta = PTDF.calculate_masked_PFV_delta(cn, PFV, VA) 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, @@ -406,20 +428,21 @@ def check_violations(mb, md, PTDF, max_viol_add, time=None, prepend_str=""): outer_name = cn, PFV = PFV, get_counts=True) total += PTDF.contingency_compensators._count[time,cn] - all_contingencies.remove(cn) if time in PTDF.contingency_compensators._order: + # give us some choice in which violations we add if total >= violations_store.max_viol_add: - print(violations_store.violations_store) - print(violations_store.max_viol_add) - print("BREAKING") + 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(all_contingencies) > 0: - print(f"ONLY checked a subset of contingencies! Number checked: {len(PTDF.contingency_compensators)-len(all_contingencies)}") + if len(PTDF.contingency_compensators) == c_checked+1: + logger.info(f"Checked **ALL** {c_checked} contingencies") else: - print(f"CHECKED ALL CONTINGENCIES, Number checked: {len(PTDF.contingency_compensators) - len(all_contingencies)}") + logger.info(f"Checked {c_checked} contingencies") else: - print("Skipping contingency flow check") + logger.info("Checked no contingencies") logger.debug(f"branches_monitored: {mb._idx_monitored}\n" f"interfaces_monitored: {mb._interfaces_monitored}\n" From 4b49dfa00c61336b56d28d480e588d61876f49ad Mon Sep 17 00:00:00 2001 From: Bernard Knueven Date: Mon, 22 May 2023 12:41:56 -0600 Subject: [PATCH 10/11] remove timers from elsewhere --- egret/common/solver_interface.py | 12 +----------- egret/data/ptdf_utils.py | 4 ---- 2 files changed, 1 insertion(+), 15 deletions(-) diff --git a/egret/common/solver_interface.py b/egret/common/solver_interface.py index f247a40a..a138b60d 100644 --- a/egret/common/solver_interface.py +++ b/egret/common/solver_interface.py @@ -83,7 +83,7 @@ def _solve_model(model, return_solver = False, vars_to_load = None, set_instance = True, - timer = None): + ): ''' Create and solve an Egret power system optimization model @@ -119,7 +119,6 @@ def _solve_model(model, ------- ''' - timer.start("model solve") results = None ## termination conditions which are acceptable @@ -150,23 +149,16 @@ def _solve_model(model, if isinstance(solver, PersistentSolver): if set_instance: - timer.start("set instance") solver.set_instance(model, symbolic_solver_labels=symbolic_solver_labels) - timer.stop("set instance") - timer.start("solve call") results = solver.solve(model, tee=solver_tee, load_solutions=False, save_results=False, **solve_method_options) - timer.stop("solve call") else: - timer.start("solve call") results = solver.solve(model, tee=solver_tee, \ symbolic_solver_labels=symbolic_solver_labels, load_solutions=False, **solve_method_options) - timer.stop("solve call") if results.solver.termination_condition not in safe_termination_conditions: raise Exception('Problem encountered during solve, termination_condition {}'.format(results.solver.termination_condition)) - timer.start("load solution") if isinstance(solver, PersistentSolver): solver.load_vars(vars_to_load) if vars_to_load is None: @@ -176,9 +168,7 @@ def _solve_model(model, solver.load_slacks() else: model.solutions.load_from(results) - timer.stop("load solution") - timer.stop("model solve") if return_solver: return model, results, solver diff --git a/egret/data/ptdf_utils.py b/egret/data/ptdf_utils.py index f12769ea..e9a8cb9f 100644 --- a/egret/data/ptdf_utils.py +++ b/egret/data/ptdf_utils.py @@ -305,7 +305,6 @@ def _set_lazy_limits(self, ptdf_options): self.enforced_interface_min_limits) def _get_ptdf_row(self, branch_name): - self.timer.start("Generate PTDF row") if branch_name in self._ptdf_rows: PTDF_row = self._ptdf_rows[branch_name] else: @@ -313,7 +312,6 @@ def _get_ptdf_row(self, branch_name): branch_idx = self._branchname_to_index_map[branch_name] PTDF_row = self.MLU.solve(self.B_dA[branch_idx].toarray(out=self._bus_sensi_buffer)[0], trans='T') self._ptdf_rows[branch_name] = PTDF_row - self.timer.stop("Generate PTDF row") return PTDF_row def _get_contingency_row(self, contingency_name, branch_name): @@ -839,9 +837,7 @@ def _set_lazy_limits(self, ptdf_options): def get_branch_ptdf_iterator(self, branch_name): row_idx = self._branchname_to_index_map[branch_name] ## get the row slice - self.timer.start("Generate PTDF row") PTDF_row = self.PTDFM[row_idx] - self.timer.stop("Generate PTDF row") yield from zip(self.buses_keys, PTDF_row) def get_branch_ptdf_abs_max(self, branch_name): From f5a197bec2dbed44c498795ae723e15f6782407e Mon Sep 17 00:00:00 2001 From: Bernard Knueven Date: Mon, 22 May 2023 15:32:34 -0600 Subject: [PATCH 11/11] setting more resonable emergency ratings --- egret/parsers/prescient_dat_parser.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) 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,