From b15175138ab8502b493fcaa395ca61579c6e515a Mon Sep 17 00:00:00 2001 From: RyanVaughan Date: Thu, 4 May 2023 18:57:28 -0500 Subject: [PATCH 1/2] adding glodf --- andes/routines/glodf.py | 244 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 244 insertions(+) create mode 100644 andes/routines/glodf.py diff --git a/andes/routines/glodf.py b/andes/routines/glodf.py new file mode 100644 index 000000000..17f8fd81d --- /dev/null +++ b/andes/routines/glodf.py @@ -0,0 +1,244 @@ +# -*- coding: utf-8 -*- +""" +Created on Mon Apr 24 19:59:15 2023 + +@author: rvaug +""" + +import andes +import pandas as pd +import numpy as np + +class GLODF: + """ + Computes the Global Line Outage Distribution Factors (GLODF) for a power system. + + Parameters + system : The power system model. + """ + def __init__(self, system): + self.system = system + self.Bs = self.get_branch_susceptance() + self.Ar = self.get_reduced_incidence() + self.Bn = self.get_reduced_nodal_susceptance() + + def get_branch_susceptance(self): + """ + Returns the branch susceptance matrix. + + Returns + ------- + numpy.ndarray : Bs + The L x L diagonal matrix of branch susceptances. + """ + return np.diag(1/self.system.Line.x.v) + + def get_reduced_incidence(self): + """ + Returns the reduced incidence matrix. + + Returns + ------- + numpy.ndarray : Ar + The L x (N-s) reduced incidence matrix where 's' is the number of + slack buses. + """ + num_line = self.system.Line.n + num_bus = self.system.Bus.n + + A = np.zeros((num_line, num_bus)) + + for l in range(num_line): + A[l,self.system.Line.bus1.v[l]-1] = -1 + A[l,self.system.Line.bus2.v[l]-1] = 1 + + # delete all slack rows as required + Ar = np.delete(A, np.asarray(self.system.Slack.bus.v)-1, axis=1) + return Ar + + def get_reduced_nodal_susceptance(self): + """ + Returns the reduced nodal susceptance matrix. + + Returns + ------- + numpy.ndarray : Bn + The (N-s) x (N-s) reduced nodal susceptance matrix where 's' is the + number of slack buses. + """ + return -1.0 * np.dot(np.dot(self.Ar.T, self.Bs), self.Ar) + + def get_isf(self): + """ + Returns the injection shift factor (ISF) matrix. + + Returns + ------- + numpy.ndarray : isf + The L x (N-s) injection shift factor matrix. + """ + psi = np.dot(np.dot(self.Bs, self.Ar), np.linalg.inv(self.Bn)) + + return psi + + def get_ptdf(self, change_lines): + """ + Returns the power transfer distribution factors for the given lines. + + Returns + ------- + numpy.ndarray : ptdf + The L x (lines) matrix. + """ + change_lines = np.atleast_1d(change_lines) + psi = self.get_isf() + + slack = np.array(self.system.Slack.bus.v) + non_slack = np.delete(np.arange(self.system.Line.n), slack - 1) + + # minus ones to make zero indexed + bfrom = np.asarray(self.system.Line.bus1.v)[change_lines - 1] + bto = np.asarray(self.system.Line.bus2.v)[change_lines - 1] + + bfrom_idx = np.zeros_like(bfrom) + bto_idx = np.zeros_like(bto) + for i in range(np.size(change_lines)): + bfrom_idx[i] = np.argwhere(non_slack == bfrom[i] - 1) + bto_idx[i] = np.argwhere(non_slack == bto[i] - 1) + + phi = psi[:,bfrom_idx] - psi[:,bto_idx] + return phi + + def lodf(self, change_line): + """ + Returns the line outage distribution factors (LODFs) matrix for single + line outage + + Returns + ------- + numpy.ndarray : glodf + The 1 x L injection shift factor matrix, where 'o' is the number of + outages. + """ + phi = self.get_ptdf(change_line) + + sigma = phi / (1 - phi[change_line - 1]) + sigma[change_line - 1] = 0 + + return sigma #[np.arange(sigma.size) != change_line - 1] + + def flow_after_lodf(self, change_line): + """ + Returns the line flows after the line outage + + Returns + ------- + numpy.ndarray : flow_after + The length L vector of line flows after the outage + """ + sigma = np.squeeze(self.lodf(change_line)) + + flow_before = self.system.Line.a1.e + + flow_after = flow_before + flow_before[change_line-1] * sigma + flow_after[change_line-1] = 0 + + return flow_after + + def glodf(self, change_lines): + """ + Returns the generalized line outage distribution factors (GLODFs) + matrix for the line outages in the change_lines list. + + Returns + ------- + numpy.ndarray : glodf + The o x L injection shift factor matrix, where 'o' is the number of + outages. + """ + change_lines = np.atleast_1d(change_lines) + + phi = self.get_ptdf(change_lines) + + # right side of equation is all lines + right_side = phi.T + + # left side is identity - Phi of change lines + Phi = right_side[:,change_lines-1] + left_side = (np.eye(np.shape(Phi)[0]) - Phi) + + xi = np.linalg.solve(left_side, right_side) + + return xi + + def flow_after_glodf(self, change_lines): + """ + Returns the line flows after the line outages given in change_lines + + Returns + ------- + numpy.ndarray : flow_after + The length L vector of line flows after the outages + """ + change_lines = np.atleast_1d(change_lines) + + flow_before = self.system.Line.a1.e + + xi = self.glodf(change_lines) + + #GLODFs times flow before + delta_flow = xi.T @ flow_before[change_lines-1] + flow_after = flow_before + delta_flow + + # ensure lines that are out have no flow + flow_after[change_lines-1] = 0 + + return flow_after + +if __name__ == "__main__": + """ + Example code to test the GLODF class + """ + print("GLODF Test") + # load system + ss = andes.load(andes.get_case("ieee14/ieee14_linetrip.xlsx")) + + # solve system + ss.PFlow.run() + + # create GLODF object + g = GLODF(ss) + + # lines to be taken out + change_lines = [5, 6, 12] + + lines_before = np.copy(ss.Line.a1.e) + lines_after = g.flow_after_glodf(change_lines) + + np.set_printoptions(precision=5) + # print("flow Before:\n" + str(ss.Line.a1.e)) + # print("flow after GLODF:\n" + str(lines_after)) + + # turn off lines and resolve + for i in range(np.size(change_lines)): + ss.Line.u.v[change_lines[i]-1] = 0 + ss.PFlow.run() + # print("flow Re-solved:\n" + str(ss.Line.a1.e)) + + lineflows = { + "bus1": ss.Line.bus1.v, + "bus2": ss.Line.bus2.v, + "P1 before": lines_before, + "P1 GLODF": lines_after, + "P1 re-solved": ss.Line.a1.e, + "error": np.abs(lines_after - ss.Line.a1.e) + } + + df_lineflow = pd.DataFrame(lineflows, index=ss.Line.idx.v) + + print(df_lineflow) + + mask = ss.Line.a1.e != 0 + mape = np.mean(np.abs((ss.Line.a1.e[mask] - lines_after[mask]) / ss.Line.a1.e[mask])) * 100 + print("mean absolute percent error: {:.3f}%".format(mape)) + \ No newline at end of file From a88975ac6f1b9734e5fa07736639cea26a7a9e01 Mon Sep 17 00:00:00 2001 From: RyanVaughan Date: Wed, 10 May 2023 16:40:15 -0500 Subject: [PATCH 2/2] converting GLODF code to use line indeces and adding test code --- andes/routines/__init__.py | 1 + andes/routines/glodf.py | 136 ++++++++++++++++++------------------- tests/test_glodf.py | 60 ++++++++++++++++ 3 files changed, 128 insertions(+), 69 deletions(-) create mode 100644 tests/test_glodf.py diff --git a/andes/routines/__init__.py b/andes/routines/__init__.py index 1cb121554..664ef412d 100644 --- a/andes/routines/__init__.py +++ b/andes/routines/__init__.py @@ -10,6 +10,7 @@ all_routines = OrderedDict([('pflow', ['PFlow']), ('tds', ['TDS']), ('eig', ['EIG']), + ('glodf', ['GLODF']), ]) class_names = list_flatten(list(all_routines.values())) diff --git a/andes/routines/glodf.py b/andes/routines/glodf.py index 17f8fd81d..2f4a5e7fb 100644 --- a/andes/routines/glodf.py +++ b/andes/routines/glodf.py @@ -1,8 +1,13 @@ -# -*- coding: utf-8 -*- """ -Created on Mon Apr 24 19:59:15 2023 +Generalized Line Outage Distribution Factors (GLODF) -@author: rvaug +The GLODFs are useful for analyzing the impact of the failure of multiple +transmission lines on the power flow of the entire system. The GLODFs are a set +of coefficients that quantify the proportional change in power flows on all +other lines due to a change in power flow on a specific line. They are +computed using the power transfer distribution factors (PTDFs), which describe +the impact of a change in power flow on a specific generator or load on the +power flows of all other generators or loads. """ import andes @@ -53,7 +58,7 @@ def get_reduced_incidence(self): A[l,self.system.Line.bus2.v[l]-1] = 1 # delete all slack rows as required - Ar = np.delete(A, np.asarray(self.system.Slack.bus.v)-1, axis=1) + Ar = np.delete(A, np.asarray(self.system.Bus.idx2uid(self.system.Slack.bus.v)), axis=1) return Ar def get_reduced_nodal_susceptance(self): @@ -85,6 +90,11 @@ def get_ptdf(self, change_lines): """ Returns the power transfer distribution factors for the given lines. + Parameters + ---------- + change_lines : line indeces + List of line indices for which the GLODFs are to be computed. + Returns ------- numpy.ndarray : ptdf @@ -93,20 +103,30 @@ def get_ptdf(self, change_lines): change_lines = np.atleast_1d(change_lines) psi = self.get_isf() - slack = np.array(self.system.Slack.bus.v) - non_slack = np.delete(np.arange(self.system.Line.n), slack - 1) + slack = np.array(self.system.Bus.idx2uid(self.system.Slack.bus.v)) + non_slack = np.delete(np.arange(self.system.Line.n), slack) - # minus ones to make zero indexed - bfrom = np.asarray(self.system.Line.bus1.v)[change_lines - 1] - bto = np.asarray(self.system.Line.bus2.v)[change_lines - 1] + bfrom = self.system.Bus.idx2uid(self.system.Line.get('bus1', change_lines)) + bto = self.system.Bus.idx2uid(self.system.Line.get('bus2', change_lines)) bfrom_idx = np.zeros_like(bfrom) bto_idx = np.zeros_like(bto) + for i in range(np.size(change_lines)): - bfrom_idx[i] = np.argwhere(non_slack == bfrom[i] - 1) - bto_idx[i] = np.argwhere(non_slack == bto[i] - 1) + if (bfrom[i] in slack): + bfrom_idx[i] = -1 + else: + bfrom_idx[i] = np.argwhere(non_slack == bfrom[i]) - phi = psi[:,bfrom_idx] - psi[:,bto_idx] + if (bto[i] in slack): + bto_idx[i] = -1 + else: + bto_idx[i] = np.argwhere(non_slack == bto[i]) + + # zeros row is needed because power injection at slack bus should yield psi = 0 + zeros_row = np.zeros((psi.shape[0], 1)) + psi_zero = np.hstack((psi, zeros_row)) + phi = psi_zero[:,bfrom_idx] - psi_zero[:,bto_idx] return phi def lodf(self, change_line): @@ -114,6 +134,11 @@ def lodf(self, change_line): Returns the line outage distribution factors (LODFs) matrix for single line outage + Parameters + ---------- + change_line : line index + line indix for which the LODFs are to be computed. + Returns ------- numpy.ndarray : glodf @@ -122,8 +147,10 @@ def lodf(self, change_line): """ phi = self.get_ptdf(change_line) - sigma = phi / (1 - phi[change_line - 1]) - sigma[change_line - 1] = 0 + uid_lines = self.system.Line.idx2uid(change_line) + + sigma = phi / (1 - phi[uid_lines]) + sigma[uid_lines] = 0 return sigma #[np.arange(sigma.size) != change_line - 1] @@ -131,6 +158,11 @@ def flow_after_lodf(self, change_line): """ Returns the line flows after the line outage + Parameters + ---------- + change_line : line index + line indix for which the LODFs are to be computed. + Returns ------- numpy.ndarray : flow_after @@ -138,10 +170,12 @@ def flow_after_lodf(self, change_line): """ sigma = np.squeeze(self.lodf(change_line)) + uid_lines = self.system.Line.idx2uid(change_line) + flow_before = self.system.Line.a1.e - flow_after = flow_before + flow_before[change_line-1] * sigma - flow_after[change_line-1] = 0 + flow_after = flow_before + flow_before[uid_lines] * sigma + flow_after[uid_lines] = 0 return flow_after @@ -150,13 +184,18 @@ def glodf(self, change_lines): Returns the generalized line outage distribution factors (GLODFs) matrix for the line outages in the change_lines list. + Parameters + ---------- + change_lines : line indices + List of line indices for which the GLODFs are to be computed. + Returns ------- numpy.ndarray : glodf The o x L injection shift factor matrix, where 'o' is the number of outages. """ - change_lines = np.atleast_1d(change_lines) + uid_lines = np.atleast_1d(self.system.Line.idx2uid(change_lines)) phi = self.get_ptdf(change_lines) @@ -164,7 +203,7 @@ def glodf(self, change_lines): right_side = phi.T # left side is identity - Phi of change lines - Phi = right_side[:,change_lines-1] + Phi = right_side[:,uid_lines] left_side = (np.eye(np.shape(Phi)[0]) - Phi) xi = np.linalg.solve(left_side, right_side) @@ -175,70 +214,29 @@ def flow_after_glodf(self, change_lines): """ Returns the line flows after the line outages given in change_lines + Parameters + ---------- + change_lines : line indices + List of line indices for which the GLODFs are to be computed. + Returns ------- numpy.ndarray : flow_after The length L vector of line flows after the outages """ - change_lines = np.atleast_1d(change_lines) + xi = self.glodf(change_lines) + + uid_lines = np.atleast_1d(self.system.Line.idx2uid(change_lines)) flow_before = self.system.Line.a1.e - xi = self.glodf(change_lines) #GLODFs times flow before - delta_flow = xi.T @ flow_before[change_lines-1] + delta_flow = xi.T @ flow_before[uid_lines] flow_after = flow_before + delta_flow # ensure lines that are out have no flow - flow_after[change_lines-1] = 0 + flow_after[uid_lines] = 0 return flow_after - -if __name__ == "__main__": - """ - Example code to test the GLODF class - """ - print("GLODF Test") - # load system - ss = andes.load(andes.get_case("ieee14/ieee14_linetrip.xlsx")) - - # solve system - ss.PFlow.run() - - # create GLODF object - g = GLODF(ss) - - # lines to be taken out - change_lines = [5, 6, 12] - - lines_before = np.copy(ss.Line.a1.e) - lines_after = g.flow_after_glodf(change_lines) - - np.set_printoptions(precision=5) - # print("flow Before:\n" + str(ss.Line.a1.e)) - # print("flow after GLODF:\n" + str(lines_after)) - - # turn off lines and resolve - for i in range(np.size(change_lines)): - ss.Line.u.v[change_lines[i]-1] = 0 - ss.PFlow.run() - # print("flow Re-solved:\n" + str(ss.Line.a1.e)) - - lineflows = { - "bus1": ss.Line.bus1.v, - "bus2": ss.Line.bus2.v, - "P1 before": lines_before, - "P1 GLODF": lines_after, - "P1 re-solved": ss.Line.a1.e, - "error": np.abs(lines_after - ss.Line.a1.e) - } - - df_lineflow = pd.DataFrame(lineflows, index=ss.Line.idx.v) - - print(df_lineflow) - - mask = ss.Line.a1.e != 0 - mape = np.mean(np.abs((ss.Line.a1.e[mask] - lines_after[mask]) / ss.Line.a1.e[mask])) * 100 - print("mean absolute percent error: {:.3f}%".format(mape)) \ No newline at end of file diff --git a/tests/test_glodf.py b/tests/test_glodf.py new file mode 100644 index 000000000..7d1363232 --- /dev/null +++ b/tests/test_glodf.py @@ -0,0 +1,60 @@ +""" +GLODF Test Code + +Test to ensure IEEE 14 bus system has less than 5% error with the following +line outages: + ['Line_1', 'Line_3', 'Line_5', 'Line_7', 'Line_9',] +""" + +from andes.routines.glodf import GLODF +import andes +import numpy as np +import pandas as pd + +def test_flow_after_glodf(): + """ + Example code to test the GLODF class + """ + print("GLODF Test") + # load system + ss = andes.load(andes.get_case("ieee14/ieee14_linetrip.xlsx")) + + # solve system + ss.PFlow.run() + + # create GLODF object + g = GLODF(ss) + + # lines to be taken out + change_lines = ['Line_1', 'Line_3', 'Line_5', 'Line_7', 'Line_9',] + + lines_before = np.copy(ss.Line.a1.e) + lines_after = g.flow_after_glodf(change_lines) + + np.set_printoptions(precision=5) + + uid_lines = np.atleast_1d(ss.Line.idx2uid(change_lines)) + # turn off lines and resolve1 + for i in range(np.size(change_lines)): + ss.Line.u.v[uid_lines] = 0 + ss.PFlow.run() + + lineflows = {"bus1": ss.Line.bus1.v, + "bus2": ss.Line.bus2.v, + "P1 before": lines_before, + "P1 GLODF": lines_after, + "P1 re-solved": ss.Line.a1.e, + "error": np.abs(lines_after - ss.Line.a1.e) } + + df_lineflow = pd.DataFrame(lineflows, index=ss.Line.idx.v) + + print(df_lineflow) + + percent_error = np.mean(np.abs((ss.Line.a1.e - lines_after))) / np.mean(np.abs(ss.Line.a1.e)) * 100 + print("normalized absolute percent error: {:.3f}%".format(percent_error)) + + # Check if lines_after is close to ss.Line.a1.e + assert np.allclose(percent_error, 0, atol=5) # error should be less than 5% + +if __name__ == "__main__": + test_flow_after_glodf()