Select Git revision
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
RBO_Service_Restoration.py 51.54 KiB
# Copyright © 2018, RWTH Aachen University
# Authors: Alberto Dognini, Abhinav Sadu
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# A copy of the GNU General Public License is in the LICENSE file
# in the top level directory of this source tree.
# coding=utf8
import sys
import os
import time
import xlsxwriter
import pandas as pd
import networkx as net
import math
import cmath
from numpy import linalg
import numpy as np
import psycopg2
from datetime import datetime,timedelta
from dsse_functions_FLISR_1 import connDB, DSSEVRI1
#CODE AT THE END OF THE FILE; AFTER FUNCTIONS AND CLASSES
#FLISR algorithm function, iterate the research of restoration path and check multiple faults/modifications of the grid:
def flisr_algorithm(prin_v, switches_excel, lines_excel, loads_excel, nodes_excel, ltype_excel, general):
######################### run the FLISR ##################################
subgrid = creation_grid(switches_excel, lines_excel)
[set_SS, alive_loads, grid_lines, grid_switches, grid_loads, grid_nodes, tripped_nodes, all_possible_loads, alive_grid] = flisr_loads(subgrid, nodes_excel, loads_excel, lines_excel, switches_excel, ltype_excel)
if all_possible_loads == 0:
[to_be_closed, final_path, load, subst, all_possible_loads] = flisr_dijkstra(set_SS, alive_loads, grid_lines, grid_switches, grid_loads, grid_nodes, ltype_excel, loads_excel, nodes_excel, alive_grid, general, prin_v)
else:
to_be_closed = 0
return to_be_closed, subgrid, set_SS, switches_excel, lines_excel, loads_excel, nodes_excel, all_possible_loads
print('No more faults happened during FLISR')
print('FLISR completed successfully')
#########################################################################################
##################### Final Updation of FLISR completion results of FLISR #########
if not final_path == []:
print('the connection between load', load, 'and substation', subst, ' is:', final_path)
print('close the switches', to_be_closed, 'immediately!!!!')
return to_be_closed, subgrid, set_SS, switches_excel, lines_excel, loads_excel, nodes_excel, all_possible_loads
else:
to_be_closed = 0
all_possible_loads == 1
print('The FLISR is concluded, no more loads can be reconnected with a suitable path')
return to_be_closed, subgrid, set_SS, switches_excel, lines_excel, loads_excel, nodes_excel, all_possible_loads
###############################################################################
### FUNCTION TO CREATE THE GRID OF CLOSED OR NO SWITCHES; TO BE USED FOR STATE ESTIMATION ###
def creation_grid(switches_excel, lines_excel):
# building the graph and verifying that the line is operating
for i, elrow in lines_excel.iterrows():
if elrow[5] == 1:
grid_lines.add_edge(int(elrow[3]), int(elrow[4]), attr_dict=elrow[0:].to_dict())
#####additional graph for switches (to check which tripped)#####
# creating the graph of switches
for i, elrow in switches_excel.iterrows():
grid_switches.add_edge(int(elrow[2]), int(elrow[3]), attr_dict=elrow[0:].to_dict())
subgrid = net.Graph((source, target, attr) for source, target, attr in grid_lines.edges.data(True)) # copy of total graph
for (p, x, d) in grid_switches.edges(data=True):
if (d['attr_dict']['status'] == 0): # removal of open switch lines
subgrid.remove_edge(p, x)
print('edges', subgrid.edges())
print('FLISR : Step 1 : Creating the real time grid topology')
return subgrid;
###############################################################################
############## FUNCTION TO IDENTIFY THE LOADS TO BE CONNECTED #################
def flisr_loads(subgrid, nodes_excel, loads_excel, lines_excel, switches_excel, ltype_excel):
# building the graph and verifying that the line is operating
for i, elrow in lines_excel.iterrows():
if elrow[5] == 1:
grid_lines.add_edge(int(elrow[3]), int(elrow[4]), attr_dict=elrow[0:].to_dict())
# add attribute column containing the line impedance and complex impedance
net.set_edge_attributes(grid_lines, 'impedance line', 1) # at the beginning, adding a new column attribute with a random value
net.set_edge_attributes(grid_lines, 'complex impedance', 1)
net.set_edge_attributes(grid_lines, 'admittance', 1)
net.set_edge_attributes(grid_lines, 'complex admittance', 1)
for (a, b, c) in grid_lines.edges(data=True):
for o, elrow in ltype_excel.iterrows():
if elrow[1] == c['attr_dict']['l_code']:
length = float(c['attr_dict']['length'])
R_direct = float(elrow[4]) * length
X_direct = float(elrow[5]) * length
B_direct = float(elrow[6]) * length
G_direct = float(elrow[7]) * length
complex_impedance = R_direct + 1j * X_direct
impedance_mag = abs(complex_impedance)
c['attr_dict']['impedance line'] = impedance_mag # then, change the value with the true one
c['attr_dict']['complex impedance'] = complex_impedance
complex_admittance = G_direct + 1j * B_direct
c['attr_dict']['complex admittance'] = complex_admittance
######additional graph for nodes (to check the reference)#####
# creating the graph of nodes
for i, elrow in nodes_excel.iterrows():
grid_nodes.add_node(int(elrow[0]), attr_dict=elrow[1:].to_dict())
n_key = []
ref_ss_val=[]
for (n, n_attr) in grid_nodes.nodes(data=True):
n_key.append(n)
ref_ss_val.append(n_attr['attr_dict']['reference_ss_node'])
ref_node= dict(zip(n_key, ref_ss_val))
# determine the number of substation present in the grid
SS = []
for (p, d) in grid_nodes.nodes(data=True):
SS.append(ref_node[p]) # list with all the reference nodes (substations) with duplicates
set_SS = set(SS) # set in order to take each element only once
num_SS = len(set_SS) # number of substations in the grid
#####additional graph for switches (to check which tripped)#####
# creating the graph of switches
for i, elrow in switches_excel.iterrows():
grid_switches.add_edge(int(elrow[2]), int(elrow[3]), attr_dict=elrow[1:].to_dict())
######additional graph for loads (to check the importance)#####
# creating the graph of loads
for i, elrow in loads_excel.iterrows():
grid_loads.add_node(int(elrow[4]), attr_dict=elrow[1:].to_dict()) # changed the dict from 3 to 1
####################################################################
######### IDENTIFY THE BRANCHES AFFECTED BY THE FAULT ##############
# finding the nodes involved, from the switches that tripped --> TO BE USED FOR THE EXTERNAL CHECK IN "FLISR Algorithm"
tripped_nodes = []
for (p, x, d) in grid_switches.edges(data=True):
if d['attr_dict']['switch_tripped'] == 0:
tripped_nodes.extend([p, x])
##### IMPORTANT!!!! in case of bus-tie, it takes also nodes on healthy branch. But successively, they are excluded since there is a path between nodes and SS (feeder switch closed)
# determine the lost loads, checking connectivity with all substations
lost_loads = []
for (p, d) in grid_nodes.nodes(data=True):
count = 0
for i in set_SS:
if not (net.has_path(subgrid, p, i)):
count += 1
if count == num_SS:
lost_loads.append(p)
##############################################################################
####### eliminate from the graph the switches tripped due to the fault #######
for (p, x, d) in grid_switches.edges(data=True):
if d['attr_dict']['switch_tripped'] == 0:
grid_lines.remove_edge(p, x)
#### remove the loads that are permanent lost (in tripped zone) ####
multiple_loads = []
for i in lost_loads: # iterates the load(s) to be restored
for x in set_SS: # iterates the substations
if net.has_path(grid_lines, x, i, ):
multiple_loads.append(i) # loads that can be connected
alive_loads = set(multiple_loads)
### define the graph to be used for the algorithm #####
# remove the tie-unit that connects two alive portions of the network
for (p, x, d) in grid_switches.edges(data=True):
if d['attr_dict']['status'] == 0: # identify the tie unit
if p not in lost_loads and x not in lost_loads: # verify that both the sides are energized
grid_lines.remove_edge(p, x)
#### remove the lost loads from grid_lines ######
permanent_lost = [val for val in lost_loads if val not in alive_loads]
alive_grid_iter = grid_lines.copy()
alive_grid = grid_lines.copy()
print('alive grid',alive_grid.edges(data=True))
for (a, b) in alive_grid_iter.nodes(data=True):
for c in permanent_lost:
if c in alive_grid.nodes() and a == c:
alive_grid.remove_node(c)
if alive_loads: # list not empty
print(
'FLISR : Step 2 : Creating the complete grid with the nodes , lines, and loads of the disconnected grid and identifying the non faulty section of the fauty feeder')
print(grid_lines.edges())
all_possible_loads = 0
else:
all_possible_loads = 1
print('all loads have been reconnected')
return set_SS, alive_loads, grid_lines, grid_switches, grid_loads, grid_nodes, tripped_nodes, all_possible_loads, alive_grid;
############## SELECT THE LOAD TO BE RESTORED #######################
def selecting_loads(grid_loads, alive_loads):
# creating a graph with all the lost loads
graph_lost_loads = grid_loads.subgraph(alive_loads)
# create a list with nodes having the highest cost grade
cost_loads = []
# increasing grade from 1 until finding the corresponding of lost loads
i = 1
while not cost_loads:
for (p, d) in graph_lost_loads.nodes(data=True):
if d['attr_dict']['cost'] == i:
cost_loads.append(int(p))
i = i + 1
# creating a graph with the max cost grade loads
graph_cost_loads = graph_lost_loads.subgraph(cost_loads)
restored_loads = [] # initialize a new list, definitive loads
# verify if there are more loads with the same grade, so to select the max power
if len(cost_loads) > 1:
i = 2e6
while not restored_loads:
for (p, d) in graph_cost_loads.nodes(data=True):
if d['attr_dict']['active_p'] == i:
restored_loads.append(p)
i = i - 1
else: # if there is only one load, doesn't check the power
restored_loads = cost_loads
print('Loads to be restored:', restored_loads)
return restored_loads
#############################################################################
def flisr_dijkstra(set_SS, alive_loads, grid_lines, grid_switches, grid_loads, grid_nodes, line_type, loads_updated, nodes_updated, alive_grid, general, prin_v):
to_be_closed = []
solution = 0
vbase = nodes_updated.loc[1, 'vbase']
pbase = float(general['pbase']) # define base power, to compute p.u. losses
###### call the function to determine the loads to be restored############
restored_loads = selecting_loads(grid_loads, alive_loads)
while solution == 0:
dec_mat = np.empty((0, 4), int)
container = []
###### CARRY OUT THE DIJKSTRA ALGORITHM TO FIND THE MINIMUM PATH #############
for x in set_SS: # iterates the substations
for i in restored_loads: # iterates the load(s) to be restored
if net.has_path(grid_lines, x,i): # necessary to verify the existence of path, Dijkstra cannot be empty
shortest_path = net.dijkstra_path(grid_lines, x, i,'impedance line') # Dijkstra compute the path according to the specific line impedance
shortest_graph = grid_lines.subgraph(shortest_path) # creating a subgraph with only lines involved in the dijkstra path
print('proposed shortest path:', shortest_graph.edges())
######### to verify that the closing path doesn't include two switches (two tie units): network no more radial ###########
count = 0 #
for (p, m, d) in shortest_graph.edges(data=True):
for (k, n, v) in grid_switches.edges(data=True):
if (p == k and m == n) or (p == n and m == k):
if (v['attr_dict']['status'] == 0):
count += 1
if count > 1:
print('No radial path between substation', x, 'and load', i) #
continue # stop the current for loop, look for a new dijkstra path
##########################################################################################################################
# remove from the graph "alive_grid" the open switches and add the proposed tie-switch
alive_grid_iter_1 = alive_grid.copy()
for (a, b, c) in grid_switches.edges(data=True):
for (d, e, f) in alive_grid_iter_1.edges(data=True):
if (a == d and b == e) and c['attr_dict']['status'] == 0: # and ((a,b) in alive_grid.nodes()):
alive_grid.remove_edge(a, b)
if (a == e and b == d) and c['attr_dict']['status'] == 0: # and ((a,b) in alive_grid.nodes()):
alive_grid.remove_edge(b, a)
alive_grid = (net.compose(alive_grid,shortest_graph)).copy() # merge witht the new shortest graph, in order to include the proposed closed tie-switch
##### INITIALIZATION OF STATE ESTIMATION MEASUREMENTS ####
error = 0
power_loss = 0
pload = 0
min_line_relat = [9999, 1, 1, 9999, 1, 1, 9999, 1, 1]
#### SUB-DIVISION OF THE NETWORK; CREATING ONE GRAPH FOR EACH SUBSTATION #####
for z in set_SS:
SS_alive_grid_iter = alive_grid.copy()
SS_alive_grid = alive_grid.copy()
for (y, t) in SS_alive_grid_iter.nodes(data=True):
if (not net.has_path(SS_alive_grid, z, y)) and y in SS_alive_grid.nodes():
SS_alive_grid.remove_node(y)
print('grid', z, 'composed by:', SS_alive_grid.edges())
### CALL THE FUNCTION TO DEFINE THE INPUT FOR THE STATE ESTIMATOR
graph_to_SE, mapping, reverse_mapping = graph_mapping_for_SE(SS_alive_grid, grid_nodes)
griddata = convert_graph_to_SE_linedata(graph_to_SE, line_type, grid_nodes, pbase, vbase)
measdata_to_SE, pload = Measurements_to_SE(mapping, grid_loads, SS_alive_grid, pbase, pload)
### STATE ESTIMATION FUNCTION - for the portion of the grid
dssedata1 = DSSEVRI1(griddata, measdata_to_SE)
temp1 = np.zeros((dssedata1.nstate, 1))
temp2 = np.zeros((dssedata1.nstate, 1))
for h in range(dssedata1.nstate):
for key in reverse_mapping.keys():
if int(dssedata1.dsseinstance[h].itemid1) == key:
temp1[h, 0] = reverse_mapping[key]
# dssedata2.dsseinstance[h].itemid1=reverse_mapping[key]
if int(dssedata1.dsseinstance[h].itemid2) == key:
temp2[h, 0] = reverse_mapping[key]
# dssedata2.dsseinstance[h].itemid2=reverse_mapping[key]
for h in range(dssedata1.nstate):
dssedata1.dsseinstance[h].itemid1 = temp1[h, 0]
dssedata1.dsseinstance[h].itemid2 = temp2[h, 0]
error, power_loss, min_line_relat = check_dijkstra(dssedata1, grid_nodes, nodes_updated,line_type, SS_alive_grid, error, power_loss, min_line_relat, pbase)
if error == 0:
### GET STATE ESTIMATION RESULTS!!!!!
# carry out computation of path losses, and add the weighting factor u1
graph_path = grid_switches.subgraph(shortest_path) # create a subgraph with edges having switches
graph_path_unfrozen = net.Graph()
graph_path_unfrozen = graph_path.copy()
for (a, b, c) in graph_path.edges(data=True):
if c['attr_dict']['status'] == 1: # keep only the open switches (to be closed)
graph_path_unfrozen.remove_edge(a, b)
graph_path = graph_path_unfrozen.copy()
to_be_closed = list(net.get_edge_attributes(graph_path, 'name').values())
subst = x
load = i
final_path = shortest_path
container.append(to_be_closed) # save the specific info of this solution
container.append(shortest_path)
container.append(i)
container.append(x)
temp_mat = np.array([power_loss, min_line_relat[0], min_line_relat[3], min_line_relat[6]]) # add row to decision matrix as new possible solution
print('most consumed line bewteen nodes', min_line_relat[1], min_line_relat[2],'is equal to', min_line_relat[0])
print('most consumed line bewteen nodes', min_line_relat[4], min_line_relat[5],'is equal to', min_line_relat[3])
print('most consumed line bewteen nodes', min_line_relat[7], min_line_relat[8],'is equal to', min_line_relat[6])
dec_mat = np.vstack((dec_mat, temp_mat))
solution = 1 # solution found
### REPEAT THE CHECK WITH HIGHER AMPACITY LIMIT
if solution == 0:
print('increase the ampacity limit by 20%')
# increase the amax limit of 20% to stay in temporary overload condition
for row, elrow in line_type.iterrows():
old_amax = elrow[2]
line_type.loc[row, 'amax'] = old_amax * 1.2
###### CARRY OUT THE DIJKSTRA ALGORITHM TO FIND THE MINIMUM PATH #############
for x in set_SS: # iterates the substations
for i in restored_loads: # iterates the load(s) to be restored
if net.has_path(grid_lines, x,i): # necessary to verify the existence of path, Dijkstra cannot be empty
shortest_path = net.dijkstra_path(grid_lines, x, i,'impedance line') # Dijkstra compute the path according to the specific line impedance
shortest_graph = grid_lines.subgraph(shortest_path) # creating a subgraph with only lines involved in the dijkstra path
print('proposed shortest path:', shortest_graph.edges())
######### to verify that the closing path doesn't include two switches (two tie units): network no more radial ###########
count = 0 #
for (p, m, d) in shortest_graph.edges(data=True):
for (k, n, v) in grid_switches.edges(data=True):
if (p == k and m == n) or (p == n and m == k):
if (v['attr_dict']['status'] == 0):
count += 1
if count > 1:
print('No radial path between substation', x, 'and load', i) #
continue # stop the current for loop, look for a new dijkstra path
##########################################################################################################################
# remove from the graph "alive_grid" the open switches and add the proposed tie-switch
alive_grid_iter = alive_grid.copy()
for (a, b, c) in grid_switches.edges(data=True):
for (d, e, f) in alive_grid_iter.edges(data=True):
if (a == d and b == e) and c['attr_dict']['status'] == 0: # and ((a,b) in alive_grid.nodes()):
alive_grid.remove_edge(a, b)
if (a == e and b == d) and c['attr_dict']['status'] == 0: # and ((a,b) in alive_grid.nodes()):
alive_grid.remove_edge(b, a)
alive_grid = (net.compose(alive_grid,shortest_graph)).copy() # merge witht the new shortest graph, in order to include the proposed closed tie-switch
##### INITIALIZATION OF STATE ESTIMATION MEASUREMENTS ####
error = 0
power_loss = 0
pload = 0
min_line_relat = [9999, 1, 1, 9999, 1, 1, 9999, 1, 1]
#### SUB-DIVISION OF THE NETWORK; CREATING ONE GRAPH FOR EACH SUBSTATION #####
for z in set_SS:
SS_alive_grid = alive_grid.copy()
SS_alive_grid_iter = alive_grid.copy()
for (y, t) in SS_alive_grid_iter.nodes(data=True):
if (not net.has_path(SS_alive_grid, z, y)) and y in SS_alive_grid.nodes():
SS_alive_grid.remove_node(y)
print('grid', z, 'composed by:', SS_alive_grid.edges())
### CALL THE FUNCTION TO DEFINE THE INPUT FOR THE STATE ESTIMATOR
graph_to_SE, mapping, reverse_mapping = graph_mapping_for_SE(SS_alive_grid, grid_nodes)
griddata = convert_graph_to_SE_linedata(graph_to_SE, line_type, grid_nodes, pbase,
vbase)
measdata_to_SE, pload = Measurements_to_SE(mapping, grid_loads, SS_alive_grid, pbase,
pload)
### STATE ESTIMATION FUNCTION - for the portion of the grid
dssedata1 = DSSEVRI1(griddata, measdata_to_SE)
temp1 = np.zeros((dssedata1.nstate, 1))
temp2 = np.zeros((dssedata1.nstate, 1))
for h in range(dssedata1.nstate):
for key in reverse_mapping.keys():
if int(dssedata1.dsseinstance[h].itemid1) == key:
temp1[h, 0] = reverse_mapping[key]
# dssedata2.dsseinstance[h].itemid1=reverse_mapping[key]
if int(dssedata1.dsseinstance[h].itemid2) == key:
temp2[h, 0] = reverse_mapping[key]
# dssedata2.dsseinstance[h].itemid2=reverse_mapping[key]
for h in range(dssedata1.nstate):
dssedata1.dsseinstance[h].itemid1 = temp1[h, 0]
dssedata1.dsseinstance[h].itemid2 = temp2[h, 0]
error, power_loss, min_line_relat = check_dijkstra(dssedata1, grid_nodes, nodes_updated,line_type, SS_alive_grid, error, power_loss, min_line_relat, pbase)
if error == 0:
### GET STATE ESTIMATION RESULTS!!!!!
# carry out computation of path losses, and add the weighting factor u1
graph_path = grid_switches.subgraph(
shortest_path) # create a subgraph with edges having switches
for (a, b, c) in graph_path.edges(data=True):
if c['attr_dict']['status'] == 1: # keep only the open switches (to be closed)
graph_path.remove_edge(a, b)
to_be_closed = list(net.get_edge_attributes(graph_path, 'name').values())
subst = x
load = i
final_path = shortest_path
container.append(to_be_closed) # save the specific info of this solution
container.append(shortest_path)
container.append(i)
container.append(x)
temp_mat = np.array([power_loss, min_line_relat[0], min_line_relat[3], min_line_relat[
6]]) # add row to decision matrix as new possible solution
dec_mat = np.vstack((dec_mat, temp_mat))
solution = 1 # solution found
if solution == 1:
print('ATTENTION!!! This load can be restored only with overload condition of the network, this is an emergency temporary configuration!')
for row, elrow in line_type.iterrows(): # restoring original nominal amax value
new_amax = elrow[2]
line_type.loc[row, 'amax'] = new_amax / 1.2
if solution == 1: # compute the best solution according to TOPSIS method
print('dec_mat', dec_mat)
row_sums = dec_mat.sum(axis=0) # normalization of decision matrix according to columns
ndec_mat = dec_mat / row_sums[np.newaxis, :]
print('ndec_mat', ndec_mat)
weig_mat = np.diag(prin_v) # creation of diagonal matrix with weights from principal eigenvector
wdec_mat = np.matmul(ndec_mat, weig_mat) # calculation of weighted normalized decision matrix
print('wdec_mat', wdec_mat)
pos_neg = [0, 1, 1,1] # 0 means negative criteria (to be minimized), LOSSES, 1 means positive criteria (to be maximized), LINE UTILIZATION FACTOR
pis = []
nis = []
for j in range(len(wdec_mat[0])): # iterate the columns (criteria)
if pos_neg[j] == 1:
pis.append(np.amax(wdec_mat[:, j])) # compute the ideal positive and negative solutions
nis.append(np.amin(wdec_mat[:, j]))
else:
pis.append(np.amin(wdec_mat[:, j]))
nis.append(np.amax(wdec_mat[:, j]))
pis_mat = np.asarray(pis) # convertion from list to array
nis_mat = np.asarray(nis)
pis_mat_1 = pis_mat # duplication for following operations
nis_mat_1 = nis_mat
print('pis_mat', pis_mat)
print('nis_mat', nis_mat)
for i in range(len(wdec_mat) - 1): # add new identical array
pis_mat = np.vstack(
(pis_mat, pis_mat_1)) # create matrix with ideal solution as number of possibilities
nis_mat = np.vstack((nis_mat, nis_mat_1))
pos_dif = wdec_mat - pis_mat
neg_dif = wdec_mat - nis_mat
sqr_pos_dif = np.square(pos_dif) # square of the matrix
sqr_neg_dif = np.square(neg_dif)
s_pos = np.sqrt(sqr_pos_dif.sum(axis=1)) # root square of the sum along the row
s_neg = np.sqrt(sqr_neg_dif.sum(axis=1))
sol = []
for i in range(len(s_pos)): # iterate the columns (possibilities)
closeness = s_neg[i] / (s_pos[i] + s_neg[i]) # calculation of closeness to ideal solution
sol.append(closeness)
print('Closeness to ideal solution:', sol)
idx_sol = np.argmax(sol) # index of solution having the highest closeness to ideal solution
to_be_closed = container[
4 * idx_sol] # retrieve the information related to the identified best solution
final_path = container[4 * idx_sol + 1]
load = container[4 * idx_sol + 2]
subst = container[4 * idx_sol + 3]
if solution == 0: # this loads cannot be connected, look for other loads
for w in restored_loads: # remove the inspected loads from the selected list
alive_loads.remove(w)
if alive_loads: # if the list is not empty, look for new loads and repeat the dijkstra
restored_loads = selecting_loads(grid_loads, alive_loads)
else: # otherwise stop the process
solution = 1
if not to_be_closed:
print('no more load can be reconnected with a suitable path')
all_possible_loads = 1
final_path = []
load = []
subst = []
else:
print(
'FLISR : Step 3 : identifying the Optimal reconfigurable path for restoring the non faulty section of faulty feeder')
all_possible_loads = 0
return to_be_closed, final_path, load, subst, all_possible_loads;
###############################################################################
###### FUNCTION TO CHECK THE RESPECT OF NETWORK RADIALITY AND LINE CAPABILITY##########
def check_dijkstra(dssedata, grid_nodes, nodes_excel,ltype_excel, SS_alive_grid, error, power_loss, min_line_relat, pbase):
#########################################################################
# HERE THE VOLTAGE LEVEL CHECK TAKES PLACE, of the whole grid
# import outcomes from STATE ESTIMATION
for o, elrow in nodes_excel.iterrows(): # iterate the nodes, all the grid must be checked
for q in range(dssedata.nstate):
if int(dssedata.dsseinstance[q].itemid1) == elrow[0] and int(dssedata.dsseinstance[q].mtype) == 1: # it should be the LN voltage, TO BE CHECKED!!!
if (float(dssedata.dsseinstance[q].estimationvalue) + 3 * float(dssedata.dsseinstance[q].estimationaccuracy)) > 1.1 or (float(dssedata.dsseinstance[q].estimationvalue) + 3 * float(dssedata.dsseinstance[q].estimationaccuracy)) < 0.9: # removed the conversion to absolute value: * float(elrow[4]): #Check voltage level in range 0.9 -> 1.1
print('Voltage level exceeds the limits at node', elrow[1], ' being ', (float(dssedata.dsseinstance[q].estimationvalue) + 3 * float(dssedata.dsseinstance[q].estimationaccuracy)),'% so this network configuration is not acceptable')
error = 1
return error, power_loss, min_line_relat
# print('In this network configuration, voltage limits are respected for all the nodes')
####### load the nodes voltage
no_key = []
vbase_val=[]
for (n, n_attr) in grid_nodes.nodes(data=True):
no_key.append(n)
vbase_val.append(n_attr['attr_dict']['vbase'])
voltages_complete = dict(zip(no_key, vbase_val))
# HERE THE CURRENT CAPACITY CHECK TAKES PLACE, of the whole grid
for (a, b, c) in SS_alive_grid.edges(data=True):
Ibase = pbase / (math.sqrt(3) * float(voltages_complete[a]))
for o, elrow in ltype_excel.iterrows():
if elrow[1] == c['attr_dict']['l_code']:
amax = elrow[2]
for r in range(dssedata.nstate):
# acquire voltage magnitude and phase (with uncertainties), for node a and b
if int(dssedata.dsseinstance[r].itemid1) == a and int(dssedata.dsseinstance[r].mtype) == 1: # and int(dssedata.dsseinstance[r].phase)=='AN': #it should be the LN voltage, TO BE CHECKED!!! phase A taken only as reference (everything symmetrical)
V1_mag = float(dssedata.dsseinstance[r].estimationvalue) * float(voltages_complete[a]) # get the p.u. voltage and multiply by nominal voltage
u_V1_mag = float(dssedata.dsseinstance[r].estimationaccuracy) * float(voltages_complete[a])
if int(dssedata.dsseinstance[r].itemid1) == a and int(dssedata.dsseinstance[r].mtype) == 3: # and int(dssedata.dsseinstance[r].phase)==A:
V1_phase = float(dssedata.dsseinstance[r].estimationvalue)
u_V1_phase = float(dssedata.dsseinstance[r].estimationaccuracy)
if int(dssedata.dsseinstance[r].itemid1) == b and int(dssedata.dsseinstance[r].mtype) == 1: # and int(dssedata.dsseinstance[r].phase)==A:
V2_mag = float(dssedata.dsseinstance[r].estimationvalue) * float(voltages_complete[a])
u_V2_mag = float(dssedata.dsseinstance[r].estimationaccuracy) * float(voltages_complete[a])
if int(dssedata.dsseinstance[r].itemid1) == b and int(dssedata.dsseinstance[r].mtype) == 3: # and int(dssedata.dsseinstance[r].phase)==A:
V2_phase = float(dssedata.dsseinstance[r].estimationvalue)
u_V2_phase = float(dssedata.dsseinstance[r].estimationaccuracy)
# acquire current magnitude and phase (with uncertainties) by comparign nodes a and b with itemid1 and itemid2 (considering both the combinations)
if ((int(dssedata.dsseinstance[r].itemid1) == a and int(dssedata.dsseinstance[r].itemid2) == b) or (int(dssedata.dsseinstance[r].itemid1) == b and int(dssedata.dsseinstance[r].itemid2) == a)) and int(dssedata.dsseinstance[r].mtype) == 5:
current_mag = float(dssedata.dsseinstance[r].estimationvalue) * Ibase # get the p.u. voltage and multiply by nominal voltage
u_current_mag = float(dssedata.dsseinstance[r].estimationaccuracy) * Ibase
# print('curr mag + acc:',current_mag, u_current_mag, 'nodes:',a, b)
if ((int(dssedata.dsseinstance[r].itemid1) == a and int(dssedata.dsseinstance[r].itemid2) == b) or (int(dssedata.dsseinstance[r].itemid1) == b and int(dssedata.dsseinstance[r].itemid2) == a)) and int(dssedata.dsseinstance[r].mtype) == 7:
current_phase = float(dssedata.dsseinstance[r].estimationvalue)
u_current_phase = float(dssedata.dsseinstance[r].estimationaccuracy)
# transform the voltage measurement from polar to rectangular form
V1_rect = cmath.rect(V1_mag, V1_phase)
V2_rect = cmath.rect(V2_mag, V2_phase)
###complex impedance and admittance
Z = c['attr_dict']['complex impedance']
Y = c['attr_dict']['complex admittance']
if (current_mag + 3 * u_current_mag) > amax: # for each edge, check that the current limit is respected
print('too much current:', current_mag + 3 * u_current_mag, 'in percentage:',
(current_mag + 3 * u_current_mag) / amax)
print('Line current exceeds the thermal limits at nodes', a, 'and', b,
':this network configuration is not acceptable')
error = 1
return error, power_loss, min_line_relat
# compute the power losses in each line
power_loss += 3 * np.real(V1_rect*(np.conjugate(V1_rect*Y/2)) + (V1_rect - V2_rect)*(np.conjugate((V1_rect - V2_rect)/Z)) + V2_rect*(np.conjugate(V2_rect*Y/2))) # sum the power losses of each edge
# compute line exploitation by current reference
line_relat = (amax - current_mag) / amax
if line_relat < min_line_relat[0]: # if supergraph.has_edge(a,b):
min_line_relat[6] = min_line_relat[3]
min_line_relat[7] = min_line_relat[4]
min_line_relat[8] = min_line_relat[5]
min_line_relat[3] = min_line_relat[0]
min_line_relat[4] = min_line_relat[1]
min_line_relat[5] = min_line_relat[2]
min_line_relat[0] = line_relat
min_line_relat[1] = a
min_line_relat[2] = b
elif line_relat < min_line_relat[3]:
min_line_relat[6] = min_line_relat[3]
min_line_relat[7] = min_line_relat[4]
min_line_relat[8] = min_line_relat[5]
min_line_relat[3] = line_relat
min_line_relat[4] = a
min_line_relat[5] = b
elif line_relat < min_line_relat[6]:
min_line_relat[6] = line_relat
min_line_relat[7] = a
min_line_relat[8] = b
return error, power_loss, min_line_relat
#functions to map the single feeders accordingly to state estimation function
def graph_mapping_for_SE(SS_alive_grid, grid_nodes):
n_key = []
ref_ss_val=[]
for (n, n_attr) in grid_nodes.nodes(data=True):
n_key.append(n)
ref_ss_val.append(n_attr['attr_dict']['reference_ss_node'])
ref_node= dict(zip(n_key, ref_ss_val))
mapping = {}
keys = SS_alive_grid.nodes()
count = 0
for bus in keys:
if (bus == ref_node[bus]):
mapping[bus] = 0
else:
count = count + 1
mapping[bus] = count
graph_to_SE = net.relabel_nodes(SS_alive_grid, mapping)
reverse_keys = list(mapping.values())
reverse_values = list(mapping.keys())
reverse_mapping = dict(zip(reverse_keys, reverse_values))
return graph_to_SE, mapping, reverse_mapping
############################ function for converting graph to topology data ########################################
def convert_graph_to_SE_linedata(graph_to_SE, ltype_excel, grid_nodes, pbasegen, vbase):
####### load the nodes voltage
n_key = []
voltages_complete_val=[]
for (n, n_attr) in grid_nodes.nodes(data=True):
n_key.append(n)
voltages_complete_val.append(n_attr['attr_dict']['vbase'])
voltages_complete = dict(zip(n_key, voltages_complete_val))
i = 0
griddata1 = griddata() # pbase is 100 kW and gridtype
griddata1.pbase = pbasegen;
griddata1.vbasenom = vbase;
griddata1.gtype = 1;
griddata1.nline = 0; # initialization to zero
griddata1.nnode = 0; # initialization to zero #iteration to populate the griddata1.nodes class
for (a, b) in graph_to_SE.nodes(data=True):
griddata1.addnode(a, voltages_complete[a])
griddata1.nnode = griddata1.nnode + 1 # iteration to populate the griddata1.edges class
for (a, b, c) in graph_to_SE.edges(data=True):
for o, elrow in ltype_excel.iterrows():
if elrow[1] == c['attr_dict']['l_code']:
griddata1.addline(c['attr_dict']['id'], 1, 1, 1, a, b, elrow[2]) # added id_line and max current value
griddata1.line[i].addphaseP11(float(elrow[4] * c['attr_dict']['length']), float(elrow[5] * c['attr_dict']['length']),float(elrow[6] * c['attr_dict']['length']), float(elrow[7] * c['attr_dict']['length']))
griddata1.line[i].addphaseP12(float(elrow[4] * 0.28 * c['attr_dict']['length']), float(elrow[5] * c['attr_dict']['length']), float(elrow[6] * c['attr_dict']['length']), float(elrow[7] * c['attr_dict']['length']))
griddata1.nline = griddata1.nline + 1
i += 1
# positive sequence in Ohm (later will be converted by the estimator in pu)
# assuming r_pos_m1 is already multiplied by the length
return griddata1
####################################### function to build the measurement object to be given to the SE ##################################
def Measurements_to_SE(mapping, grid_loads, SS_alive_grid, pbase, pload):
Nmeas = 0;
measdata_to_SE = measdata()
############# retrieve the active_power and reactive_power from the graph of loadst mimicing the power injection values (i real time)####
n_key = []
active_power_complete_val=[]
reactive_power_complete_val = []
for (n, n_attr) in grid_loads.nodes(data=True):
n_key.append(n)
active_power_complete_val.append(n_attr['attr_dict']['active_p'])
reactive_power_complete_val.append(n_attr['attr_dict']['active_p'])
active_power_complete = dict(zip(n_key, active_power_complete_val))
reactive_power_complete = dict(zip(n_key, reactive_power_complete_val))
#active_power_complete = net.get_node_attributes(grid_loads, 'active_p')
#reactive_power_complete = net.get_node_attributes(grid_loads, 'reactive_q')
# Setting the measurement class object of SE with the voltage mag and anng at substation buses and injection values at all other buses
for (node, b) in SS_alive_grid.nodes(data=True):
if (mapping[node] == 0): # include the voltage measurement for the SS node or nodes with 0 injections
measdata_to_SE.addmeasurementinstance(Nmeas, 1, 1, mapping[node], 0.01,1) # append with a voltage magnitude instance
Nmeas = Nmeas + 1;
measdata_to_SE.addmeasurementinstance(Nmeas, 3, 1, mapping[node], 0.01,0) # append with a voltage angle instance
Nmeas = Nmeas + 1;
elif (active_power_complete[node] != 0 and reactive_power_complete[node] != 0):
measdata_to_SE.addmeasurementinstance(Nmeas, 9, 1, mapping[node], 0.01414, active_power_complete[node] / pbase) # append with an active power injection instance
Nmeas = Nmeas + 1;
pload += active_power_complete[node]
measdata_to_SE.addmeasurementinstance(Nmeas, 11, 1, mapping[node], 0.01414, reactive_power_complete[node] / pbase) # append with a reactive power injection instance
Nmeas = Nmeas + 1;
measdata_to_SE.nmeas = Nmeas
return measdata_to_SE, pload
#detect fault function:
def detect_fault(switches_excel):
is_fault_detected = 0 #variable initiated
for i, elrow in switches_excel.iterrows():
if elrow[2] == 0: # if there is at least a CB that tripped, the FLISR is started
is_fault_detected = 1
return is_fault_detected
########### DEFINE CLASSES IN ORDER TO INITIALIZE STATE ESTIMATOR MEASUREMENTS #########
class dssedata(object):
def __init__(self, nstate=0, numnr=20, time=1):
self.nstate = nstate
self.numnr = numnr
self.time = time
self.dsseinstance = []
def adddsseinstance(self, estimationid, mtype, phase, itemid1, itemid2, estimationaccuracy, estimationvalue):
self.dsseinstance.append(
dsseinstance(estimationid, mtype, phase, itemid1, itemid2, estimationaccuracy, estimationvalue))
class dsseinstance(object):
def __init__(self, estimationid, mtype, phase, itemid1, itemid2, estimationaccuracy, estimationvalue):
self.estimationid = estimationid
self.mtype = mtype
self.phase = phase
self.itemid1 = itemid1
self.itemid2 = itemid2
self.estimationaccuracy = estimationaccuracy
self.estimationvalue = estimationvalue
################## definition of class objects for generating the measurement data for SE ##########################
class measdata(object):
def __init__(self, nmeas=0, time=0, switchphmeas=0):
self.nmeas = nmeas
self.time = time
self.switchphmeas = switchphmeas
self.measurementinstance = []
def addmeasurementinstance(self, measurementid, mtype, phase, itemid, measurementaccuracy, measurementvalue):
self.measurementinstance.append(
measurementinstance(measurementid, mtype, phase, itemid, measurementaccuracy, measurementvalue))
class measurementinstance(object):
def __init__(self, measurementid, mtype, phase, itemid, measurementaccuracy, measurementvalue):
self.measurementid = measurementid
self.mtype = mtype
self.phase = phase
self.itemid = itemid
self.measurementaccuracy = measurementaccuracy
self.measurementvalue = measurementvalue
class griddata(object):
def __init__(self, pbase=1000000, gtype=1, nnode=0, nline=0, nswitch=0, vbasenom=8660):
self.pbase = pbase
self.vbasenom = vbasenom
self.nline = nline
self.nnode = nnode
self.nswitch = nswitch
self.gtype = gtype
self.node = []
self.line = []
def addnode(self, nodeid, vbase):
self.node.append(node(nodeid, vbase))
def addline(self, lineid, phasea, phaseb, phasec, linefrom, lineto, amax):
self.line.append(line(lineid, phasea, phaseb, phasec, linefrom, lineto, amax))
def addswitch(self, switchid, nodefrom, nodeto, switchtripped, status, amax):
self.node.append(switch(switchid, nodefrom, nodeto, switchtripped, status, amax))
class node(object):
def __init__(self, nodeid, vbase=8660):
self.nodeid = nodeid
self.vbase = vbase
class switch(object):
def __init__(self, switchid, nodefrom, nodeto, switchtripped, status, amax):
self.switchid = switchid
self.nodefrom = nodefrom
self.nodeto = nodeto
self.switchtripped = switchtripped
self.status = status
self.amax = amax
class line(object):
def __init__(self, lineid, phasea, phaseb, phasec, linefrom, lineto, amax):
self.lineid = lineid
self.phasea = phasea
self.phaseb = phaseb
self.phasec = phasec
self.linefrom = linefrom
self.lineto = lineto
self.amax = amax
self.phaseP11 = []
self.phaseP12 = []
self.phaseP13 = []
self.phaseP22 = []
self.phaseP23 = []
self.phaseP33 = []
def addphaseP11(self, resistance, reactance, susceptance, conductance):
self.phaseP11.append(phaseP11(resistance, reactance, susceptance, conductance))
def addphaseP12(self, resistance, reactance, susceptance, conductance):
self.phaseP12.append(phaseP12(resistance, reactance, susceptance, conductance))
def addphaseP13(self, resistance, reactance, susceptance, conductance):
self.phaseP13.append(phaseP13(resistance, reactance, susceptance, conductance))
def addphaseP22(self, resistance, reactance, susceptance, conductance):
self.phaseP22.append(phaseP22(resistance, reactance, susceptance, conductance))
def addphaseP23(self, resistance, reactance, susceptance, conductance):
self.phaseP23.append(phaseP23(resistance, reactance, susceptance, conductance))
def addphaseP33(self, resistance, reactance, susceptance, conductance):
self.phaseP33.append(phaseP33(resistance, reactance, susceptance, conductance))
class phaseP11(object):
def __init__(self, resistance, reactance, susceptance, conductance):
self.resistance = resistance
self.reactance = reactance
self.susceptance = susceptance
self.conductance = conductance
class phaseP12(object):
def __init__(self, resistance, reactance, susceptance, conductance):
self.resistance = resistance
self.reactance = reactance
self.susceptance = susceptance
self.conductance = conductance
class phaseP13(object):
def __init__(self, resistance, reactance, susceptance, conductance):
self.resistance = resistance
self.reactance = reactance
self.susceptance = susceptance
self.conductance = conductance
class phaseP22(object):
def __init__(self, resistance, reactance, susceptance, conductance):
self.resistance = resistance
self.reactance = reactance
self.susceptance = susceptance
self.conductance = conductance
class phaseP23(object):
def __init__(self, resistance, reactance, susceptance, conductance):
self.resistance = resistance
self.reactance = reactance
self.susceptance = susceptance
self.conductance = conductance
class phaseP33(object):
def __init__(self, resistance, reactance, susceptance, conductance):
self.resistance = resistance
self.reactance = reactance
self.susceptance = susceptance
self.conductance = conductance
class databaseconn(object):
def __init__(self,DBname,user,password,ipaddress,port):
self.DBname = DBname
self.user = user
self.password = password
self.ipaddress = ipaddress
self.port = port
directory=os.path.dirname(os.path.abspath(__file__))
file_name='/networkdata_FLISR.xlsx'
real_status_excel=directory+file_name
############################### update the database ##############################
switches_excel = pd.read_excel(real_status_excel,sheet_name='switches')
lines_excel = pd.read_excel(real_status_excel,sheet_name='lines')
loads_excel = pd.read_excel(real_status_excel,sheet_name='loads')
nodes_excel = pd.read_excel(real_status_excel,sheet_name='nodes')
ltype_excel = pd.read_excel(real_status_excel,sheet_name='ltype')
general = pd.read_excel(real_status_excel,sheet_name='general')
###################################### Create empty graph#####################################################
grid_lines=net.Graph()
grid_switches=net.Graph()
grid_loads=net.Graph()
grid_nodes=net.Graph()
#### COMPARISON FACTORS TO BE DEFINED!!!!!
# In the scale from 1 (equally important) to 9 (extremely more important) or to 1/9 (extremely less important)
uL1=7 #Power losses VS 1st line utilization
uL2=8 #Power losses VS 2nd line utilization
uL3=9 #Power losses VS 3rd line utilization
u12=2 #1st line utilization VS 2nd line utilization (>1)
u13=2 #1st line utilization VS 3rd line utilization (>1)
u23=2 #2nd line utilization VS 3rd line utilization (>1)
# random consistency index
rand_idx=0.9 ##ATTENTION: this value is 0.9 ONLY in case of 4 attributes to be compared
#### ANALYTICAL HIERARCHY PROCESS
# Creation of comparison matrix
# Criteria: Power Losses, 1st utilized line, 2nd utilized line, 3rd utilized line
comp_mat= np.array([[1, uL1, uL2, uL3],[1/uL1, 1, u12, u13],[1/uL2, 1/u12, 1, u23],[1/uL3, 1/u13, 1/u23, 1]])
dim=len(comp_mat)
w,v = linalg.eig(comp_mat) #function to extract eigenvalues and eigenvectors
prin_w=np.amax(w) #calculation of principal eigenvalue
idx=np.argmax(w) #index of max eigenvalue
prin_v=v[:,idx] #calculation of principal eigenvector
prin_v=prin_v/prin_v.sum(axis=0,keepdims=1) #normalization of principal eigenvector
print('this is the principal, normalized AHP eigenvector:', prin_v)
const_ratio=(prin_w - dim)/(rand_idx*(dim-1)) #check of consistency ratio
if const_ratio >= 0.1:
print('Consistency Ratio:',const_ratio)
print('The comparison factors are not consistent, chenge them!')
print('Service Restoration aborted')
sys.exit()
############################
################################# detect in real time if there was a fault ##########################################
is_fault_detected=detect_fault(switches_excel)
################## if fault detected update flag and start the FLISR ##########################################
if is_fault_detected==1:
all_possible_loads=0
while all_possible_loads==0:
is_FLISR_in_operation =1
to_be_closed, subgrid, set_SS, switches_excel, lines_excel, loads_excel, nodes_excel, all_possible_loads = flisr_algorithm(prin_v, switches_excel, lines_excel, loads_excel, nodes_excel, ltype_excel, general)
if all_possible_loads==0:
for i, elrow in switches_excel.iterrows():
################ update the switches status (close the corresponding tie-unit) in the DB ###################
for j in range(len(to_be_closed)):
switches_excel.loc[switches_excel['name'] == to_be_closed[j], 'status'] = 1
############### update the SS reference for each node (if a new connection exists) in the DB #################
grid_nodes=net.Graph()
for i, elrow in nodes_excel.iterrows():
grid_nodes.add_node(int(elrow[0]), attr_dict=elrow[1:].to_dict())
for i in set_SS:
for (p, d) in grid_nodes.nodes(data=True):
if (net.has_path(subgrid, p, i)):
nodes_excel.loc[nodes_excel['id'] == p, 'reference_ss_node'] = i
else:
print('no fault identified: enjoy the healthy grid :-)')
sys.exit()