main.py

import numpy as np
import pandapower as pp
import pandapower.networks as networks
import pandapower.topology as top
import networkx as nx
from choquet_main import choquet_integral
from plotting_functions import plotting_function
from operator import itemgetter
import random
import matplotlib.pyplot as plt
from G_to_net import G_to_net
#import numba
import matplotlib as mpl
mpl.use('Qt5Agg')
random.seed(1)


# Function to create and initialize a network graph from the pandapower network data
def create_network_graph(net):
    G = nx.Graph()
    # Add nodes with voltage level and set initial generator/load flags
    for bus in net.bus.itertuples():
        G.add_node(bus.Index, vn_kv=bus.vn_kv, is_gen_bus=False, is_load_bus=False)

    # Mark nodes as generator buses
    for element in [net.gen, net.sgen, net.ext_grid]:
        for item in element.itertuples():
            if hasattr(item, 'bus'):
                G.nodes[item.bus]['is_gen_bus'] = True

    # Mark nodes as load buses
    for load in net.load.itertuples():
        G.nodes[load.bus]['is_load_bus'] = True

    # Add edges with calculated impedance and maximum power flow
    for line in net.line.itertuples():
        if line.in_service:
            from_bus_vn_kv = G.nodes[line.from_bus]['vn_kv']
            to_bus_vn_kv = G.nodes[line.to_bus]['vn_kv']
            average_vn_kv = (from_bus_vn_kv + to_bus_vn_kv) / 2
            impedance = np.sqrt(line.r_ohm_per_km ** 2 + line.x_ohm_per_km ** 2) * line.length_km
            max_power_flow = line.max_i_ka * np.sqrt(3) * average_vn_kv
            G.add_edge(line.from_bus, line.to_bus, impedance=impedance, max_power_flow=max_power_flow)
    return G


# Calculate the admittance matrix Ybus
def calculate_Ybus(net):
    Ybus = net._ppc['internal']['Ybus']
    return Ybus.toarray()


# Ybus = calculate_Ybus(net)

# Calculate the Zbus matrix by inverting Ybus
def calculate_Zbus(Ybus):
    Zbus = np.linalg.inv(Ybus)
    return Zbus


# Zbus = calculate_Zbus(Ybus)

# Now calculate Zdg using Zii - 2Zij + Zjj
def calculate_Zdg(Zbus):
    num_buses = Zbus.shape[0]
    Zdg = np.zeros((num_buses, num_buses), dtype=complex)
    for i in range(num_buses):
        for j in range(num_buses):
            if i != j:
                Zdg[i, j] = abs(Zbus[i, i] - 2 * Zbus[i, j] + Zbus[j, j])
    return Zdg


# Calculate Zdg matrix
# Zdg = calculate_Zdg(Zbus)

def calculate_B_prime_inv(Ybus):
    B = Ybus.imag
    B_prime = B[0:, 0:]  # make it 0
    B_prime_inv = np.linalg.inv(B_prime)
    return [B, B_prime_inv]


# [B, B_prime_inv] = calculate_B_prime_inv(Ybus)

def count_network_elements(net):
    num_lines = len(net.line.index)
    num_buses = len(net.bus.index)
    return [num_lines, num_buses]


# [num_lines,num_buses] = count_network_elements(net)

def calculate_Hprime(B, net):
    num_buses = len(net.bus.index)
    num_lines = len(net.line.index)
    H = np.zeros((num_lines, num_buses))  # Rows = lines, Columns = buses

    # Populate H matrix based on the connectivity and B matrix values
    for idx, line in enumerate(net.line.itertuples()):
        if line.in_service:
            from_bus = line.from_bus
            to_bus = line.to_bus
            # Hij equals to B[i, j] for the from_bus to to_bus
            # Hji equals to -B[i, j] for the to_bus to from_bus
            try:
                H[idx, from_bus] = B[from_bus, to_bus]
                H[idx, to_bus] = -B[from_bus, to_bus]
            except:
                H[idx, from_bus] = 0
                H[idx, to_bus] = 0
    # return H

    H_prime = H[:, 0:]  # This slices out the first column # make 0

    return H_prime


# H_prime = calculate_Hprime (B,net)


# Compute PTDF
def calculate_PTDF(H_prime, B_prime_inv):

    try:
        F = np.dot(H_prime, B_prime_inv)
    except:
        temp = np.zeros((300,300)) # making size to 300 to get product of matrix
        temp[:len(B_prime_inv),:len(B_prime_inv)] = B_prime_inv
        F = np.dot(H_prime, temp)

    return F


# F = calculate_PTDF (H_prime, B_prime_inv)


def calculate_Fldg(F, d, g):
    # Ensure d and g are within the bounds of the matrix columns
    if d >= F.shape[1] or g >= F.shape[1]:
        raise ValueError("Column indices d or g are out of bounds.")

    # Calculate Fldg for all rows l
    Fldg = F[:, d] - F[:, g]
    return Fldg  # vector


def get_individual_buses(G):
    # Retrieve all generator buses
    G_tilda = [bus for bus, data in G.nodes(data=True) if data.get('is_gen_bus', True)]
    # Retrieve all load buses
    L_tilda = [bus for bus, data in G.nodes(data=True) if data.get('is_load_bus', True)]

    return [G_tilda, L_tilda]


# [G_tilda, L_tilda] = get_individual_buses(G)


def calculate_Cdg(G, fldg):
    # Initialize an empty list to store the Cdg values for each line
    Cdg_values = []

    # Iterate over each line in the graph
    for idx, (u, v, data) in enumerate(G.edges(data=True)):
        max_power_flow = data['max_power_flow']
        # fldg = fldg[idx] if idx < len(fldg) else 0  # Safeguard against index out of range
        fldg_l = fldg[idx]

        # Prevent division by zero
        if fldg_l == 0:  # if 0, append infinity if  zero, use formula and append the value
            Cdg_line = float('inf')
            Cdg_values.append(Cdg_line)
        else:

            # Calculate Cdg for this line
            Cdg_line = abs(max_power_flow / fldg_l)
            Cdg_values.append(Cdg_line)

    # Calculate the minimum Cdg value from the list
    Cdg = min(Cdg_values) if Cdg_values else float('inf')  # Handle empty list scenario

    return Cdg


def compute_net_ability(NG, NL, F, Zdg, G_tilda, L_tilda,G):
    A = 0
    # Loop over all generator-load pairs to calculate their contribution
    for i in range(len(G_tilda)):  # G_tilde
        g = G_tilda[i]
        for j in range(len(L_tilda)):  # L_tilde
            d = L_tilda[j]
            try:
                Z_dg = abs(Zdg[g, g] - 2 * Zdg[g, d] + Zdg[d, d])
            except:
                temp = np.ones((300, 300))*np.max(Zdg)  # making size to 300 to get product of matrix and replacing missing entries by max
                temp[:len(Zdg), :len(Zdg)] = Zdg
                Zdg = temp
                Z_dg = abs(Zdg[g, g] - 2 * Zdg[g, d] + Zdg[d, d])

            fl_dg = calculate_Fldg(F=F, d=d, g=g)
            C_dg = calculate_Cdg(G=G, fldg=fl_dg)
            if Z_dg == 0 or C_dg == float('inf'):
                A += 0
            else:
                A += C_dg / Z_dg  # Ensure to take only the real part
            #A /= (NG * NL)
    return A/(NG * NL)


def calculate_electrical_extended_betweenness(G, Zdg, num_lines, F, G_tilda, L_tilda):  # major correction needed here
    EEB = {bus: 0 for bus in G.nodes()}  # Initialize EEB for all buses

    # Get the list of edges with an assumed order or index
    edges_with_indices = {edge: idx for idx, edge in enumerate(G.edges())}

    # for loop for v belong s to G.nodes

    '''for k in range (len(G.nodes())):
        v = list(G.nodes())[k]
        EEB_temp = 0
    # Loop over all generator-load pairs to calculate their contribution
        for i in range(len(G_tilda)): #G_tilde
            g = G_tilda[i]
            for j in range(len(L_tilda)): # L_tilde
                d = L_tilda[j]
                #Z_dg= Zdg[g,g] - 2*Zdg[g,d] + Zdg[d,d]
                fl_dg = calculate_Fldg (F=F, d=d, g=g)
                C_dg = calculate_Cdg (G=G, fldg=fl_dg)
                # For each bus V, iterate over all lines connected to V
                Lv = list(G.edges(v))  # need index of line ??
                temp = 0
                for (u, w) in Lv:
                    if (u == g and w == d) or (u == d and w == g):  # updated
                        # Find the index for the line (u, w)
                        line_index = edges_with_indices[(u, w)] if (u, w) in edges_with_indices else edges_with_indices[(w, u)]
                        # Use the line index to fetch the appropriate values from F matrix
                        temp += abs(F[line_index, g] - F[line_index, d])
                #temp += fl_dg[L_tilda]  ## find f_ldg  # This assumes C_dg is scalar and the same for all lines, adjust if needed
                EEB_temp += 0.5 * C_dg * temp
        EEB[v] = EEB_temp'''

    # Calculate the average extended betweenness
    #average_EEB = sum(EEB.values()) / num_buses if EEB else 0

    # FIND index of line which is most critical using topological approach for given graph X,

    # Find the bus with the maximum EEB value
    # max_EEB_bus = max(EEB, key=EEB.get) if EEB else None
    edge_betweenness = nx.edge_betweenness_centrality(G, weight=None)  # weight=None to consider unweighted
    average_EEB = sum(edge_betweenness.values())/len(edge_betweenness.values())
    # Find the edge with the highest betweenness centrality
    edge_to_remove = max(edge_betweenness, key=edge_betweenness.get)

    return [average_EEB, edge_to_remove[0], edge_to_remove[1]]


def calculate_entropic_degree(G, num_buses):
    entropic_degrees = {}
    for i in G.nodes():
        wij_list = [G[i][j]['max_power_flow'] for j in G.neighbors(i)]  # i. and ii.
        sum_wij = sum(wij_list)
        pij_list = [wij / sum_wij for wij in wij_list]  # iii.

        # iv. Calculate entropic degree for node i
        Si = (1 - sum(pij * np.log(pij) if pij > 0 else 0 for pij in pij_list)) * sum_wij
        entropic_degrees[i] = Si

    # Calculate average entropic degree
    average_entropic_degree = sum(entropic_degrees.values()) / num_buses  # divide it by number of buses in base case
    return average_entropic_degree


# Calculate and print the average entropic degree
# average_entropic_degree = calculate_entropic_degree(G)
# print(f'Average Entropic Degree: {average_entropic_degree}')


def perform_attacks(NG, NL, num_buses, num_lines, num_lines_to_fail=30):
    # removed_lines_list_index = []
    target_attack_results = {"Entropic": [], "net_ability": [], "electrical_betweenness": []}

    for attack_round in range(0, num_lines_to_fail):
        print(f"target_attack iteration:{attack_round}")
        if attack_round == 0:

            net = networks.case300()
            net.line.x_ohm_per_km = abs(net.line.x_ohm_per_km)
            pp.runpp(net,tolerance_mva=1e100)
            G = create_network_graph(net)
            Ybus = calculate_Ybus(net)
            Zbus = calculate_Zbus(Ybus)
            Zdg = calculate_Zdg(Zbus)
            [B, B_prime_inv] = calculate_B_prime_inv(Ybus)
            H_prime = calculate_Hprime(B, net)
            F = calculate_PTDF(H_prime, B_prime_inv)
            [G_tilda, L_tilda] = get_individual_buses(G)

            average_entropic_degree = calculate_entropic_degree(G=G, num_buses=num_buses)
            [EEB_avg, u, v] = calculate_electrical_extended_betweenness(G=top.create_nxgraph(net), Zdg=Zdg, F=F, num_lines = num_lines, G_tilda=G_tilda,
                                                                        L_tilda=L_tilda)

            average_netability = compute_net_ability(NG=NG, NL=NL, F=F, Zdg=Zdg, G_tilda=G_tilda, L_tilda=L_tilda,G = G)

            target_attack_results["Entropic"].append(average_entropic_degree)
            target_attack_results["net_ability"].append(average_netability)
            target_attack_results["electrical_betweenness"].append(EEB_avg)

            # removing critical line
            #G.remove_edge(u, v)  # u and v is index of bus to which line is connected
            try:
                net.line = net.line.drop(net.line[((net.line.from_bus == u) & (net.line.to_bus == v))].index)
            except:
                net.line = net.line.drop(net.line[((net.line.from_bus == v) & (net.line.to_bus == u))].index)
        else:

            #net = G_to_net(G)
            pp.runpp(net,tolerance_mva=1e100)
            G = create_network_graph(net)
            Ybus = calculate_Ybus(net)
            Zbus = calculate_Zbus(Ybus)
            Zdg = calculate_Zdg(Zbus)
            [B, B_prime_inv] = calculate_B_prime_inv(Ybus=Ybus)
            H_prime = calculate_Hprime(B=B, net=net)
            F = calculate_PTDF(H_prime=H_prime, B_prime_inv=B_prime_inv)
            [G_tilda, L_tilda] = get_individual_buses(G=G)

            average_entropic_degree = calculate_entropic_degree(G=G, num_buses=num_buses)
            [EEB_avg, u, v] = calculate_electrical_extended_betweenness(G=top.create_nxgraph(net), Zdg=Zdg, F=F, num_lines = num_lines, G_tilda=G_tilda,
                                                                        L_tilda=L_tilda)
            average_netability = compute_net_ability(NG=NG, NL=NL, F=F, Zdg=Zdg, G_tilda=G_tilda, L_tilda=L_tilda,G = G)

            # removing critical line
            #G.remove_edge(u, v)  # u and v is index of bus to which line is connected

            try:
                net.line = net.line.drop(net.line[((net.line.from_bus == u) & (net.line.to_bus == v))].index)
            except:
                net.line = net.line.drop(net.line[((net.line.from_bus == v) & (net.line.to_bus == u))].index)

            target_attack_results["Entropic"].append(average_entropic_degree / target_attack_results["Entropic"][0])
            target_attack_results["net_ability"].append(average_netability / target_attack_results["net_ability"][0])
            target_attack_results["electrical_betweenness"].append(target_attack_results["electrical_betweenness"][0]/EEB_avg)
    return target_attack_results
    # Output the results for each attack round
    # print(f'After failing {attack_round} lines:')
    # print(f'Average Entropic Degree: {average_entropic_degree}')
    # print(f'Average Electrical Betweenness: {EEB_avg}')
    # print(f'Net Ability: {average_netability}')
    # print('-----------------------')

### for random attack

def random_attacks(NG, NL, num_buses, num_lines, num_lines_to_fail):
    # removed_lines_list_index = []
    random_attack_results = {"Entropic": [], "net_ability": [], "electrical_betweenness": []}

    for attack_round in range(0, num_lines_to_fail):
        print(f"random_attack iteration:{attack_round}")
        if attack_round == 0:

            net = networks.case300()
            net.line.x_ohm_per_km = abs(net.line.x_ohm_per_km)
            pp.runpp(net, tolerance_mva=1e100)
            G = create_network_graph(net)
            Ybus = calculate_Ybus(net)
            Zbus = calculate_Zbus(Ybus)
            Zdg = calculate_Zdg(Zbus)
            [B, B_prime_inv] = calculate_B_prime_inv(Ybus)
            H_prime = calculate_Hprime(B, net)
            F = calculate_PTDF(H_prime, B_prime_inv)
            [G_tilda, L_tilda] = get_individual_buses(G)

            average_entropic_degree = calculate_entropic_degree(G=G, num_buses=num_buses)
            [EEB_avg, u, v] = calculate_electrical_extended_betweenness(G=top.create_nxgraph(net), Zdg=Zdg, F=F, G_tilda=G_tilda,
                                                                        L_tilda=L_tilda,num_lines = num_lines)

            average_netability = compute_net_ability(NG=NG, NL=NL, F=F, Zdg=Zdg, G_tilda=G_tilda, L_tilda=L_tilda, G=G)

            random_attack_results["Entropic"].append(average_entropic_degree)
            random_attack_results["net_ability"].append(average_netability)
            random_attack_results["electrical_betweenness"].append(EEB_avg)

            # removing critical line
            # G.remove_edge(u, v)  # u and v is index of bus to which line is connected
            net.line = net.line.drop(random.choice(list(net.line.index)))
        else:

            # net = G_to_net(G)
            pp.runpp(net, tolerance_mva=1e100)
            G = create_network_graph(net)
            Ybus = calculate_Ybus(net)
            Zbus = calculate_Zbus(Ybus)
            Zdg = calculate_Zdg(Zbus)
            [B, B_prime_inv] = calculate_B_prime_inv(Ybus=Ybus)
            H_prime = calculate_Hprime(B=B, net=net)
            F = calculate_PTDF(H_prime=H_prime, B_prime_inv=B_prime_inv)
            [G_tilda, L_tilda] = get_individual_buses(G=G)

            average_entropic_degree = calculate_entropic_degree(G=G, num_buses=num_buses)
            [EEB_avg, u, v] = calculate_electrical_extended_betweenness(G=top.create_nxgraph(net), Zdg=Zdg, F=F, num_lines = num_lines, G_tilda=G_tilda,
                                                                        L_tilda=L_tilda)
            average_netability = compute_net_ability(NG=NG, NL=NL, F=F, Zdg=Zdg, G_tilda=G_tilda, L_tilda=L_tilda, G=G)

            # removing critical line
            # G.remove_edge(u, v)  # u and v is index of bus to which line is connected


            random_attack_results["Entropic"].append(average_entropic_degree / random_attack_results["Entropic"][0])
            random_attack_results["net_ability"].append(average_netability / random_attack_results["net_ability"][0])
            random_attack_results["electrical_betweenness"].append(random_attack_results["electrical_betweenness"][0]/EEB_avg)
            # choose line to damage randomly

            net.line = net.line.drop(random.choice(list(net.line.index)))

            #random_edge = random.choice(list(G.edges()))
            #G.remove_edge(*random_edge)

    return random_attack_results


# Load the network and run power flow
net = networks.case300()
pp.runpp(net,tolerance_mva=1e100)

# Initialize the network graph
G = create_network_graph(net)

## global main
num_lines = len(net.line.index)
num_buses = len(net.bus.index)
NG = sum(1 for _, node_data in G.nodes(data=True) if node_data['is_gen_bus'])
NL = sum(1 for _, node_data in G.nodes(data=True) if node_data['is_load_bus'])
num_iterations = 10
# targeted attack
target_attack_results = perform_attacks(NG=NG, NL=NL, num_buses=num_buses, num_lines = num_lines, num_lines_to_fail=num_iterations)

# Perform the random attacks
random_attack_results = random_attacks(NG=NG, NL=NL, num_lines = num_lines, num_buses=num_buses, num_lines_to_fail=num_iterations)

# making base value to 0 as it is normalized wrt base case
target_attack_results["Entropic"][0] = 1
target_attack_results["net_ability"][0] = 1
target_attack_results["electrical_betweenness"][0] = 1
random_attack_results["Entropic"][0] = 1
random_attack_results["net_ability"][0] = 1
random_attack_results["electrical_betweenness"][0] = 1

## Getting aggregated value
target_attack_aggregated = []
random_attack_aggregated = []
for i in range(num_iterations):
    target_attack_values = [target_attack_results["net_ability"][i], target_attack_results["Entropic"][i], target_attack_results["electrical_betweenness"][i]]
    target_attack_aggregated.append(choquet_integral(target_attack_values))

    random_attack_values = [random_attack_results["net_ability"][i], random_attack_results["Entropic"][i],
                            random_attack_results["electrical_betweenness"][i]]
    random_attack_aggregated.append(choquet_integral(random_attack_values))


## analysis of result
plotting_function(plot_title="Net-Ability",target_data = list(target_attack_results["net_ability"]),random_data = list(random_attack_results["net_ability"]),figure_number=1)
plotting_function(plot_title="Entropic Degree",target_data = list(target_attack_results["Entropic"]),random_data = list(random_attack_results["Entropic"]),figure_number=2)
plotting_function(plot_title="Electrical Betweenness",target_data = list(target_attack_results["electrical_betweenness"]),random_data = list(random_attack_results["electrical_betweenness"]),figure_number=3)
plotting_function(plot_title="Aggregated Metric",target_data = target_attack_aggregated, random_data = random_attack_aggregated,figure_number=4)
plt.show()

1+1