# -*- coding: utf-8 -*-
"""
Created on Wed Dec  2 09:53:35 2020

@author: Andrea Bellè
"""

import numpy as np
import collections
import os
import pandapower as pp
import pandapower.networks
import networkx as nx
import time
import math
import sys
#import matplotlib as mpl
#import matplotlib.pyplot as plt
#import seaborn as sb
#import multiprocessing  
#import pandapower.plotting as pplt
#from mpl_toolkits.axes_grid1 import make_axes_locatable
#from sklearn.linear_model import LinearRegression
#from sklearn.metrics import mean_squared_error, r2_score
#from matplotlib.patches import Patch
#from matplotlib.lines import Line2D
#from glob import glob
#from pylab import *
#from matplotlib import cm
#from mpl_toolkits.mplot3d import Axes3D
#from mpl_toolkits.basemap import Basemap
#from mpi4py import MPI
#import random
#import itertools
'''
import conda
conda_file_dir = conda.__file__
conda_dir = conda_file_dir.split('lib')[0]
proj_lib = os.path.join(os.path.join(conda_dir, 'share'), 'proj')
os.environ["PROJ_LIB"] = proj_lib
'''

##############################################################################
#time_start = time.process_time()
##############################################################################

# INITIALIZATION OF THE POWER NETWORK

# Call the Great Britain reduced power network from Pandapower
net = pandapower.networks.GBreducednetwork()

# Replace static generators with standard generators
pp.replace_sgen_by_gen(net)
list_ext_grid  = list(net.ext_grid.index)
initial_slack = pp.replace_ext_grid_by_gen(net, list_ext_grid, slack=True) 
    
# Set the maximum loading of each line to 100%   
net.line.loc[:,'max_loading_percent'] = 100
net.trafo.loc[:,'max_loading_percent'] = 100

# Set the min/max power values of each load
net.load.loc[:,'max_p_mw'] = net.load.p_mw
net.load.loc[:,'min_p_mw'] = 0

# Set the min/max power values of each generator
net.gen.loc[:,'max_p_mw'] = net.gen.max_p_mw/1.438
net.gen.loc[:,'p_mw'] = net.gen.max_p_mw

# Set generators and loads as controllable (for DC OPF)
net.gen.loc[:,'controllable'] = True
net.load.loc[:,'controllable'] = True

# Set cost function for each generator
net.poly_cost.loc[:,'cp0_eur'] = 0
net.poly_cost.loc[:,'cp1_eur_per_mw'] = 1
net.poly_cost.loc[:,'cp2_eur_per_mw2'] = 0

# Set cost function for each load
list_load =list(net.load.index)
for load in list_load:
    pp.create_poly_cost(net, load, 'load',cp0_eur=0,
                        cp1_eur_per_mw = -100,cp2_eur_per_mw2=0)

# Compute the initial power of each load to determine the load shedding
initial_state_load = dict(net.load.loc[:,'p_mw'])
order_dict_load = collections.OrderedDict(sorted(initial_state_load.items()))
initial_load_power = list(order_dict_load.values())


# CASCADING FAILURES FUNCTION
'''
INPUT: 
      - net: initial power network configuration in Pandapower format
      - f_bus: indexes of initial failed buses 
      - f_gen: indexes of initial failed generators
      - f_load: indexes of initial failed loads
      - f_line: indexes of initial failed lines
      - f_trafo: indexes of initial failed transformers 
OUTPUT: 
      - final power production at each generator
      - final power demand at each load 
'''
def OPA_model(net, f_bus, f_gen, f_load, f_line, f_trafo):  
    
    # This variable allows to stop the function when a total blackout situation is reached
    total_blackout = 0
    
    # Set failed elements as out of service.
    # In the case of failed buses, all the elements connected to the bus are also removed
    f_bus = list(f_bus)
    pp.set_element_status(net, f_bus, in_service=False)
    f_gen = list(f_gen)
    net.gen.loc[f_gen,'in_service'] = False
    f_load = list(f_load)
    net.load.loc[f_load,'in_service'] = False
    f_line = list(f_line)
    net.line.loc[f_line,'in_service'] = False
    f_trafo = list(f_trafo)
    net.trafo.loc[f_trafo,'in_service'] = False
    
    # Variable for entering in the while loop
    cascade = True
    
    # Store the initial power level of each generator
    
    intermediate_gen_power = dict(net.gen.loc[:,'p_mw'])
    order_dict_gen = collections.OrderedDict(sorted(intermediate_gen_power.items()))
    intermediate_gen_power = list(order_dict_gen.values())
    
    intermediate_load_power = dict(net.load.loc[:,'p_mw'])
    order_dict_load = collections.OrderedDict(sorted(intermediate_load_power.items()))
    intermediate_load_power = list(order_dict_load.values())
    
    while cascade == True:
        
        # Create a networkx graph of the power network
        G=pp.topology.create_nxgraph(net)
        # Create a list of sets containing the islands of the network
        list_islands = list(nx.connected_components(G))  
        
        # Create a list of buses with generators in service
        list_bus_with_gen_in_service = list(net.gen.bus[net.gen.in_service==True])
        list_bus_with_gen_in_service = list(set(list_bus_with_gen_in_service))
        list_bus_with_gen_in_service.sort()
        list_bus_with_slack_gen_in_service = list(net.gen.bus[net.gen.in_service==True][net.gen.slack==True])

        # Check the configuration of each island in the power network
        for island in list_islands:
            #print(island)
            island = list(island)
            
            # Check if the island has an external grid. If yes, it is ok and you can proceed to the next island.
            # If no, check if there is an availabl generator in the island. If yes, turn it into an external grid.
            # If no, go to the next island.
            if set(list_bus_with_slack_gen_in_service).isdisjoint(island):
                if set(list_bus_with_gen_in_service).isdisjoint(island):
                    pass
                else:
                    set_island = set(island)
                    bus_with_gen_in_island = list(set_island.intersection(list_bus_with_gen_in_service))
                    slack_gen = list(net.gen.index[net.gen.bus==bus_with_gen_in_island[0]])
                    net.gen.loc[slack_gen[0],'slack'] = True
            else:
                pass
            
        # Set the isolated islands (without ext. grid/slack gen.) out of service
        pp.set_isolated_areas_out_of_service(net)   
        
        # Try to run a DC OPF
        try:
            pp.rundcopp(net, verbose=False)
            # If it converges, save a the power level at each generators and ext grids
            intermediate_gen_power = dict(net.res_gen.loc[:,'p_mw'])
            order_dict_gen = collections.OrderedDict(sorted(intermediate_gen_power.items()))
            intermediate_gen_power = list(order_dict_gen.values())
            
            intermediate_load_power = dict(net.res_load.loc[:,'p_mw'])
            order_dict_load = collections.OrderedDict(sorted(intermediate_load_power.items()))
            intermediate_load_power = list(order_dict_load.values())
            
        # If the DC OPF does not converge, increase the loads costs 
        # A while loop which decreases the loads costs until the DC OPF converges or the costs reach a value of 0
        except pp.OPFNotConverged:
            #print(f_line, f_bus, f_trafo)
            print('CONVERGENCE PROBLEMS')
            cost = 1
            load_cost = -100
            while cost == 1:
                load_cost += 2
                #print(load_cost)
                net.poly_cost.cp1_eur_per_mw[net.poly_cost.et == 'load'] = load_cost
                try:
                    pp.rundcopp(net)
                    #pp.rundcopp(net, SCPDIPM_RED_IT=100, PDIPM_COSTTOL = 1e-3, PDIPM_GRADTOL = 1e-3)
                    intermediate_gen_power = dict(net.res_gen.loc[:,'p_mw'])
                    order_dict_gen = collections.OrderedDict(sorted(intermediate_gen_power.items()))
                    intermediate_gen_power = list(order_dict_gen.values())

                    intermediate_load_power = dict(net.res_load.loc[:,'p_mw'])
                    order_dict_load = collections.OrderedDict(sorted(intermediate_load_power.items()))
                    intermediate_load_power = list(order_dict_load.values())

                    cost = 0
                except pp.OPFNotConverged:
                    pass
                # If the loads costs are set to 0, run a DC power flow with the last available generators power levels
                if load_cost == 0:
                    net.gen.loc[:,'p_mw'] = intermediate_gen_power
                    net.load.loc[:,'p_mw'] = intermediate_load_power
                    pp.rundcpp(net)
                    break
                
        # This is in case to generators is available (everything is failed basically)
        except UserWarning:
            total_blackout = 1
            break
        # If there is a total blackout situation, break the loop
        if total_blackout == 1:
            break
        
        # Create a list of lines still in service
        list_line_in_service = list(net.line.loc[net.line.in_service== True].index)
        list_line_in_service.sort()
        # This variable is to check if there are new failures due to overloads
        new_failure = 0
        for line in list_line_in_service:
            level_loading = net.res_line.loc[line, 'loading_percent']
            # If a line a loading >= the 99% of its maximum capacity, it is considered overloaded
            if level_loading >= 99:
                print('OVERLOADS')
                new_failure = 1
                net.line.loc[line, 'in_service'] = False
        
        # Same for transformer lines
        list_trafo_in_service = list(net.trafo.loc[net.trafo.in_service== True].index)
        list_trafo_in_service.sort()
        for trafo in list_trafo_in_service:
            level_loading = net.res_trafo.loc[trafo, 'loading_percent']
            if level_loading >= 99:
                print('OVERLOADS')
                new_failure = 1
                net.trafo.loc[trafo, 'in_service'] = False
        
        # If there are no overloads, break the while loop
        if new_failure == 0:
            break
    
    # Save the final power levels of each generator and load
    if total_blackout == 0:        
        final_state_load = dict(net.res_load.loc[:,'p_mw'])
        order_dict_load = collections.OrderedDict(sorted(final_state_load.items()))
        list_load_power = list(order_dict_load.values())
        
        final_state_gen = dict(net.res_gen.loc[:,'p_mw'])
        order_dict_gen = collections.OrderedDict(sorted(final_state_gen.items()))
        list_gen_power = list(order_dict_gen.values())
        
        list_gen_power = np.array(list_gen_power)
        list_load_power = np.array(list_load_power)       
        
        for i,j in enumerate(list_gen_power):
            if math.isnan(j):
                list_gen_power[i] = 0
        for i,j in enumerate(list_load_power):
            if math.isnan(j):
                list_load_power[i] = 0
    elif total_blackout == 1:
        list_gen_power = np.zeros(len(net.gen.index))
        list_load_power = np.zeros(len(net.load.index))
        
    return(list_gen_power,list_load_power)



'''
g,l = OPA_model(net, [8,9], [], [], [], [])   
print(np.around(l/initial_load_power, 2))
'''


6f6884bc5e664f3e4cab05cae6abb06198daf998
<mask token>


def OPA_model(net, f_bus, f_gen, f_load, f_line, f_trafo):
    total_blackout = 0
    f_bus = list(f_bus)
    pp.set_element_status(net, f_bus, in_service=False)
    f_gen = list(f_gen)
    net.gen.loc[f_gen, 'in_service'] = False
    f_load = list(f_load)
    net.load.loc[f_load, 'in_service'] = False
    f_line = list(f_line)
    net.line.loc[f_line, 'in_service'] = False
    f_trafo = list(f_trafo)
    net.trafo.loc[f_trafo, 'in_service'] = False
    cascade = True
    intermediate_gen_power = dict(net.gen.loc[:, 'p_mw'])
    order_dict_gen = collections.OrderedDict(sorted(intermediate_gen_power.
        items()))
    intermediate_gen_power = list(order_dict_gen.values())
    intermediate_load_power = dict(net.load.loc[:, 'p_mw'])
    order_dict_load = collections.OrderedDict(sorted(
        intermediate_load_power.items()))
    intermediate_load_power = list(order_dict_load.values())
    while cascade == True:
        G = pp.topology.create_nxgraph(net)
        list_islands = list(nx.connected_components(G))
        list_bus_with_gen_in_service = list(net.gen.bus[net.gen.in_service ==
            True])
        list_bus_with_gen_in_service = list(set(list_bus_with_gen_in_service))
        list_bus_with_gen_in_service.sort()
        list_bus_with_slack_gen_in_service = list(net.gen.bus[net.gen.
            in_service == True][net.gen.slack == True])
        for island in list_islands:
            island = list(island)
            if set(list_bus_with_slack_gen_in_service).isdisjoint(island):
                if set(list_bus_with_gen_in_service).isdisjoint(island):
                    pass
                else:
                    set_island = set(island)
                    bus_with_gen_in_island = list(set_island.intersection(
                        list_bus_with_gen_in_service))
                    slack_gen = list(net.gen.index[net.gen.bus ==
                        bus_with_gen_in_island[0]])
                    net.gen.loc[slack_gen[0], 'slack'] = True
            else:
                pass
        pp.set_isolated_areas_out_of_service(net)
        try:
            pp.rundcopp(net, verbose=False)
            intermediate_gen_power = dict(net.res_gen.loc[:, 'p_mw'])
            order_dict_gen = collections.OrderedDict(sorted(
                intermediate_gen_power.items()))
            intermediate_gen_power = list(order_dict_gen.values())
            intermediate_load_power = dict(net.res_load.loc[:, 'p_mw'])
            order_dict_load = collections.OrderedDict(sorted(
                intermediate_load_power.items()))
            intermediate_load_power = list(order_dict_load.values())
        except pp.OPFNotConverged:
            print('CONVERGENCE PROBLEMS')
            cost = 1
            load_cost = -100
            while cost == 1:
                load_cost += 2
                net.poly_cost.cp1_eur_per_mw[net.poly_cost.et == 'load'
                    ] = load_cost
                try:
                    pp.rundcopp(net)
                    intermediate_gen_power = dict(net.res_gen.loc[:, 'p_mw'])
                    order_dict_gen = collections.OrderedDict(sorted(
                        intermediate_gen_power.items()))
                    intermediate_gen_power = list(order_dict_gen.values())
                    intermediate_load_power = dict(net.res_load.loc[:, 'p_mw'])
                    order_dict_load = collections.OrderedDict(sorted(
                        intermediate_load_power.items()))
                    intermediate_load_power = list(order_dict_load.values())
                    cost = 0
                except pp.OPFNotConverged:
                    pass
                if load_cost == 0:
                    net.gen.loc[:, 'p_mw'] = intermediate_gen_power
                    net.load.loc[:, 'p_mw'] = intermediate_load_power
                    pp.rundcpp(net)
                    break
        except UserWarning:
            total_blackout = 1
            break
        if total_blackout == 1:
            break
        list_line_in_service = list(net.line.loc[net.line.in_service == 
            True].index)
        list_line_in_service.sort()
        new_failure = 0
        for line in list_line_in_service:
            level_loading = net.res_line.loc[line, 'loading_percent']
            if level_loading >= 99:
                print('OVERLOADS')
                new_failure = 1
                net.line.loc[line, 'in_service'] = False
        list_trafo_in_service = list(net.trafo.loc[net.trafo.in_service == 
            True].index)
        list_trafo_in_service.sort()
        for trafo in list_trafo_in_service:
            level_loading = net.res_trafo.loc[trafo, 'loading_percent']
            if level_loading >= 99:
                print('OVERLOADS')
                new_failure = 1
                net.trafo.loc[trafo, 'in_service'] = False
        if new_failure == 0:
            break
    if total_blackout == 0:
        final_state_load = dict(net.res_load.loc[:, 'p_mw'])
        order_dict_load = collections.OrderedDict(sorted(final_state_load.
            items()))
        list_load_power = list(order_dict_load.values())
        final_state_gen = dict(net.res_gen.loc[:, 'p_mw'])
        order_dict_gen = collections.OrderedDict(sorted(final_state_gen.
            items()))
        list_gen_power = list(order_dict_gen.values())
        list_gen_power = np.array(list_gen_power)
        list_load_power = np.array(list_load_power)
        for i, j in enumerate(list_gen_power):
            if math.isnan(j):
                list_gen_power[i] = 0
        for i, j in enumerate(list_load_power):
            if math.isnan(j):
                list_load_power[i] = 0
    elif total_blackout == 1:
        list_gen_power = np.zeros(len(net.gen.index))
        list_load_power = np.zeros(len(net.load.index))
    return list_gen_power, list_load_power


<mask token>

<mask token>
pp.replace_sgen_by_gen(net)
<mask token>
for load in list_load:
    pp.create_poly_cost(net, load, 'load', cp0_eur=0, cp1_eur_per_mw=-100,
        cp2_eur_per_mw2=0)
<mask token>


def OPA_model(net, f_bus, f_gen, f_load, f_line, f_trafo):
    total_blackout = 0
    f_bus = list(f_bus)
    pp.set_element_status(net, f_bus, in_service=False)
    f_gen = list(f_gen)
    net.gen.loc[f_gen, 'in_service'] = False
    f_load = list(f_load)
    net.load.loc[f_load, 'in_service'] = False
    f_line = list(f_line)
    net.line.loc[f_line, 'in_service'] = False
    f_trafo = list(f_trafo)
    net.trafo.loc[f_trafo, 'in_service'] = False
    cascade = True
    intermediate_gen_power = dict(net.gen.loc[:, 'p_mw'])
    order_dict_gen = collections.OrderedDict(sorted(intermediate_gen_power.
        items()))
    intermediate_gen_power = list(order_dict_gen.values())
    intermediate_load_power = dict(net.load.loc[:, 'p_mw'])
    order_dict_load = collections.OrderedDict(sorted(
        intermediate_load_power.items()))
    intermediate_load_power = list(order_dict_load.values())
    while cascade == True:
        G = pp.topology.create_nxgraph(net)
        list_islands = list(nx.connected_components(G))
        list_bus_with_gen_in_service = list(net.gen.bus[net.gen.in_service ==
            True])
        list_bus_with_gen_in_service = list(set(list_bus_with_gen_in_service))
        list_bus_with_gen_in_service.sort()
        list_bus_with_slack_gen_in_service = list(net.gen.bus[net.gen.
            in_service == True][net.gen.slack == True])
        for island in list_islands:
            island = list(island)
            if set(list_bus_with_slack_gen_in_service).isdisjoint(island):
                if set(list_bus_with_gen_in_service).isdisjoint(island):
                    pass
                else:
                    set_island = set(island)
                    bus_with_gen_in_island = list(set_island.intersection(
                        list_bus_with_gen_in_service))
                    slack_gen = list(net.gen.index[net.gen.bus ==
                        bus_with_gen_in_island[0]])
                    net.gen.loc[slack_gen[0], 'slack'] = True
            else:
                pass
        pp.set_isolated_areas_out_of_service(net)
        try:
            pp.rundcopp(net, verbose=False)
            intermediate_gen_power = dict(net.res_gen.loc[:, 'p_mw'])
            order_dict_gen = collections.OrderedDict(sorted(
                intermediate_gen_power.items()))
            intermediate_gen_power = list(order_dict_gen.values())
            intermediate_load_power = dict(net.res_load.loc[:, 'p_mw'])
            order_dict_load = collections.OrderedDict(sorted(
                intermediate_load_power.items()))
            intermediate_load_power = list(order_dict_load.values())
        except pp.OPFNotConverged:
            print('CONVERGENCE PROBLEMS')
            cost = 1
            load_cost = -100
            while cost == 1:
                load_cost += 2
                net.poly_cost.cp1_eur_per_mw[net.poly_cost.et == 'load'
                    ] = load_cost
                try:
                    pp.rundcopp(net)
                    intermediate_gen_power = dict(net.res_gen.loc[:, 'p_mw'])
                    order_dict_gen = collections.OrderedDict(sorted(
                        intermediate_gen_power.items()))
                    intermediate_gen_power = list(order_dict_gen.values())
                    intermediate_load_power = dict(net.res_load.loc[:, 'p_mw'])
                    order_dict_load = collections.OrderedDict(sorted(
                        intermediate_load_power.items()))
                    intermediate_load_power = list(order_dict_load.values())
                    cost = 0
                except pp.OPFNotConverged:
                    pass
                if load_cost == 0:
                    net.gen.loc[:, 'p_mw'] = intermediate_gen_power
                    net.load.loc[:, 'p_mw'] = intermediate_load_power
                    pp.rundcpp(net)
                    break
        except UserWarning:
            total_blackout = 1
            break
        if total_blackout == 1:
            break
        list_line_in_service = list(net.line.loc[net.line.in_service == 
            True].index)
        list_line_in_service.sort()
        new_failure = 0
        for line in list_line_in_service:
            level_loading = net.res_line.loc[line, 'loading_percent']
            if level_loading >= 99:
                print('OVERLOADS')
                new_failure = 1
                net.line.loc[line, 'in_service'] = False
        list_trafo_in_service = list(net.trafo.loc[net.trafo.in_service == 
            True].index)
        list_trafo_in_service.sort()
        for trafo in list_trafo_in_service:
            level_loading = net.res_trafo.loc[trafo, 'loading_percent']
            if level_loading >= 99:
                print('OVERLOADS')
                new_failure = 1
                net.trafo.loc[trafo, 'in_service'] = False
        if new_failure == 0:
            break
    if total_blackout == 0:
        final_state_load = dict(net.res_load.loc[:, 'p_mw'])
        order_dict_load = collections.OrderedDict(sorted(final_state_load.
            items()))
        list_load_power = list(order_dict_load.values())
        final_state_gen = dict(net.res_gen.loc[:, 'p_mw'])
        order_dict_gen = collections.OrderedDict(sorted(final_state_gen.
            items()))
        list_gen_power = list(order_dict_gen.values())
        list_gen_power = np.array(list_gen_power)
        list_load_power = np.array(list_load_power)
        for i, j in enumerate(list_gen_power):
            if math.isnan(j):
                list_gen_power[i] = 0
        for i, j in enumerate(list_load_power):
            if math.isnan(j):
                list_load_power[i] = 0
    elif total_blackout == 1:
        list_gen_power = np.zeros(len(net.gen.index))
        list_load_power = np.zeros(len(net.load.index))
    return list_gen_power, list_load_power


<mask token>

<mask token>
net = pandapower.networks.GBreducednetwork()
pp.replace_sgen_by_gen(net)
list_ext_grid = list(net.ext_grid.index)
initial_slack = pp.replace_ext_grid_by_gen(net, list_ext_grid, slack=True)
net.line.loc[:, 'max_loading_percent'] = 100
net.trafo.loc[:, 'max_loading_percent'] = 100
net.load.loc[:, 'max_p_mw'] = net.load.p_mw
net.load.loc[:, 'min_p_mw'] = 0
net.gen.loc[:, 'max_p_mw'] = net.gen.max_p_mw / 1.438
net.gen.loc[:, 'p_mw'] = net.gen.max_p_mw
net.gen.loc[:, 'controllable'] = True
net.load.loc[:, 'controllable'] = True
net.poly_cost.loc[:, 'cp0_eur'] = 0
net.poly_cost.loc[:, 'cp1_eur_per_mw'] = 1
net.poly_cost.loc[:, 'cp2_eur_per_mw2'] = 0
list_load = list(net.load.index)
for load in list_load:
    pp.create_poly_cost(net, load, 'load', cp0_eur=0, cp1_eur_per_mw=-100,
        cp2_eur_per_mw2=0)
initial_state_load = dict(net.load.loc[:, 'p_mw'])
order_dict_load = collections.OrderedDict(sorted(initial_state_load.items()))
initial_load_power = list(order_dict_load.values())
<mask token>


def OPA_model(net, f_bus, f_gen, f_load, f_line, f_trafo):
    total_blackout = 0
    f_bus = list(f_bus)
    pp.set_element_status(net, f_bus, in_service=False)
    f_gen = list(f_gen)
    net.gen.loc[f_gen, 'in_service'] = False
    f_load = list(f_load)
    net.load.loc[f_load, 'in_service'] = False
    f_line = list(f_line)
    net.line.loc[f_line, 'in_service'] = False
    f_trafo = list(f_trafo)
    net.trafo.loc[f_trafo, 'in_service'] = False
    cascade = True
    intermediate_gen_power = dict(net.gen.loc[:, 'p_mw'])
    order_dict_gen = collections.OrderedDict(sorted(intermediate_gen_power.
        items()))
    intermediate_gen_power = list(order_dict_gen.values())
    intermediate_load_power = dict(net.load.loc[:, 'p_mw'])
    order_dict_load = collections.OrderedDict(sorted(
        intermediate_load_power.items()))
    intermediate_load_power = list(order_dict_load.values())
    while cascade == True:
        G = pp.topology.create_nxgraph(net)
        list_islands = list(nx.connected_components(G))
        list_bus_with_gen_in_service = list(net.gen.bus[net.gen.in_service ==
            True])
        list_bus_with_gen_in_service = list(set(list_bus_with_gen_in_service))
        list_bus_with_gen_in_service.sort()
        list_bus_with_slack_gen_in_service = list(net.gen.bus[net.gen.
            in_service == True][net.gen.slack == True])
        for island in list_islands:
            island = list(island)
            if set(list_bus_with_slack_gen_in_service).isdisjoint(island):
                if set(list_bus_with_gen_in_service).isdisjoint(island):
                    pass
                else:
                    set_island = set(island)
                    bus_with_gen_in_island = list(set_island.intersection(
                        list_bus_with_gen_in_service))
                    slack_gen = list(net.gen.index[net.gen.bus ==
                        bus_with_gen_in_island[0]])
                    net.gen.loc[slack_gen[0], 'slack'] = True
            else:
                pass
        pp.set_isolated_areas_out_of_service(net)
        try:
            pp.rundcopp(net, verbose=False)
            intermediate_gen_power = dict(net.res_gen.loc[:, 'p_mw'])
            order_dict_gen = collections.OrderedDict(sorted(
                intermediate_gen_power.items()))
            intermediate_gen_power = list(order_dict_gen.values())
            intermediate_load_power = dict(net.res_load.loc[:, 'p_mw'])
            order_dict_load = collections.OrderedDict(sorted(
                intermediate_load_power.items()))
            intermediate_load_power = list(order_dict_load.values())
        except pp.OPFNotConverged:
            print('CONVERGENCE PROBLEMS')
            cost = 1
            load_cost = -100
            while cost == 1:
                load_cost += 2
                net.poly_cost.cp1_eur_per_mw[net.poly_cost.et == 'load'
                    ] = load_cost
                try:
                    pp.rundcopp(net)
                    intermediate_gen_power = dict(net.res_gen.loc[:, 'p_mw'])
                    order_dict_gen = collections.OrderedDict(sorted(
                        intermediate_gen_power.items()))
                    intermediate_gen_power = list(order_dict_gen.values())
                    intermediate_load_power = dict(net.res_load.loc[:, 'p_mw'])
                    order_dict_load = collections.OrderedDict(sorted(
                        intermediate_load_power.items()))
                    intermediate_load_power = list(order_dict_load.values())
                    cost = 0
                except pp.OPFNotConverged:
                    pass
                if load_cost == 0:
                    net.gen.loc[:, 'p_mw'] = intermediate_gen_power
                    net.load.loc[:, 'p_mw'] = intermediate_load_power
                    pp.rundcpp(net)
                    break
        except UserWarning:
            total_blackout = 1
            break
        if total_blackout == 1:
            break
        list_line_in_service = list(net.line.loc[net.line.in_service == 
            True].index)
        list_line_in_service.sort()
        new_failure = 0
        for line in list_line_in_service:
            level_loading = net.res_line.loc[line, 'loading_percent']
            if level_loading >= 99:
                print('OVERLOADS')
                new_failure = 1
                net.line.loc[line, 'in_service'] = False
        list_trafo_in_service = list(net.trafo.loc[net.trafo.in_service == 
            True].index)
        list_trafo_in_service.sort()
        for trafo in list_trafo_in_service:
            level_loading = net.res_trafo.loc[trafo, 'loading_percent']
            if level_loading >= 99:
                print('OVERLOADS')
                new_failure = 1
                net.trafo.loc[trafo, 'in_service'] = False
        if new_failure == 0:
            break
    if total_blackout == 0:
        final_state_load = dict(net.res_load.loc[:, 'p_mw'])
        order_dict_load = collections.OrderedDict(sorted(final_state_load.
            items()))
        list_load_power = list(order_dict_load.values())
        final_state_gen = dict(net.res_gen.loc[:, 'p_mw'])
        order_dict_gen = collections.OrderedDict(sorted(final_state_gen.
            items()))
        list_gen_power = list(order_dict_gen.values())
        list_gen_power = np.array(list_gen_power)
        list_load_power = np.array(list_load_power)
        for i, j in enumerate(list_gen_power):
            if math.isnan(j):
                list_gen_power[i] = 0
        for i, j in enumerate(list_load_power):
            if math.isnan(j):
                list_load_power[i] = 0
    elif total_blackout == 1:
        list_gen_power = np.zeros(len(net.gen.index))
        list_load_power = np.zeros(len(net.load.index))
    return list_gen_power, list_load_power


<mask token>

<mask token>
import numpy as np
import collections
import os
import pandapower as pp
import pandapower.networks
import networkx as nx
import time
import math
import sys
<mask token>
net = pandapower.networks.GBreducednetwork()
pp.replace_sgen_by_gen(net)
list_ext_grid = list(net.ext_grid.index)
initial_slack = pp.replace_ext_grid_by_gen(net, list_ext_grid, slack=True)
net.line.loc[:, 'max_loading_percent'] = 100
net.trafo.loc[:, 'max_loading_percent'] = 100
net.load.loc[:, 'max_p_mw'] = net.load.p_mw
net.load.loc[:, 'min_p_mw'] = 0
net.gen.loc[:, 'max_p_mw'] = net.gen.max_p_mw / 1.438
net.gen.loc[:, 'p_mw'] = net.gen.max_p_mw
net.gen.loc[:, 'controllable'] = True
net.load.loc[:, 'controllable'] = True
net.poly_cost.loc[:, 'cp0_eur'] = 0
net.poly_cost.loc[:, 'cp1_eur_per_mw'] = 1
net.poly_cost.loc[:, 'cp2_eur_per_mw2'] = 0
list_load = list(net.load.index)
for load in list_load:
    pp.create_poly_cost(net, load, 'load', cp0_eur=0, cp1_eur_per_mw=-100,
        cp2_eur_per_mw2=0)
initial_state_load = dict(net.load.loc[:, 'p_mw'])
order_dict_load = collections.OrderedDict(sorted(initial_state_load.items()))
initial_load_power = list(order_dict_load.values())
<mask token>


def OPA_model(net, f_bus, f_gen, f_load, f_line, f_trafo):
    total_blackout = 0
    f_bus = list(f_bus)
    pp.set_element_status(net, f_bus, in_service=False)
    f_gen = list(f_gen)
    net.gen.loc[f_gen, 'in_service'] = False
    f_load = list(f_load)
    net.load.loc[f_load, 'in_service'] = False
    f_line = list(f_line)
    net.line.loc[f_line, 'in_service'] = False
    f_trafo = list(f_trafo)
    net.trafo.loc[f_trafo, 'in_service'] = False
    cascade = True
    intermediate_gen_power = dict(net.gen.loc[:, 'p_mw'])
    order_dict_gen = collections.OrderedDict(sorted(intermediate_gen_power.
        items()))
    intermediate_gen_power = list(order_dict_gen.values())
    intermediate_load_power = dict(net.load.loc[:, 'p_mw'])
    order_dict_load = collections.OrderedDict(sorted(
        intermediate_load_power.items()))
    intermediate_load_power = list(order_dict_load.values())
    while cascade == True:
        G = pp.topology.create_nxgraph(net)
        list_islands = list(nx.connected_components(G))
        list_bus_with_gen_in_service = list(net.gen.bus[net.gen.in_service ==
            True])
        list_bus_with_gen_in_service = list(set(list_bus_with_gen_in_service))
        list_bus_with_gen_in_service.sort()
        list_bus_with_slack_gen_in_service = list(net.gen.bus[net.gen.
            in_service == True][net.gen.slack == True])
        for island in list_islands:
            island = list(island)
            if set(list_bus_with_slack_gen_in_service).isdisjoint(island):
                if set(list_bus_with_gen_in_service).isdisjoint(island):
                    pass
                else:
                    set_island = set(island)
                    bus_with_gen_in_island = list(set_island.intersection(
                        list_bus_with_gen_in_service))
                    slack_gen = list(net.gen.index[net.gen.bus ==
                        bus_with_gen_in_island[0]])
                    net.gen.loc[slack_gen[0], 'slack'] = True
            else:
                pass
        pp.set_isolated_areas_out_of_service(net)
        try:
            pp.rundcopp(net, verbose=False)
            intermediate_gen_power = dict(net.res_gen.loc[:, 'p_mw'])
            order_dict_gen = collections.OrderedDict(sorted(
                intermediate_gen_power.items()))
            intermediate_gen_power = list(order_dict_gen.values())
            intermediate_load_power = dict(net.res_load.loc[:, 'p_mw'])
            order_dict_load = collections.OrderedDict(sorted(
                intermediate_load_power.items()))
            intermediate_load_power = list(order_dict_load.values())
        except pp.OPFNotConverged:
            print('CONVERGENCE PROBLEMS')
            cost = 1
            load_cost = -100
            while cost == 1:
                load_cost += 2
                net.poly_cost.cp1_eur_per_mw[net.poly_cost.et == 'load'
                    ] = load_cost
                try:
                    pp.rundcopp(net)
                    intermediate_gen_power = dict(net.res_gen.loc[:, 'p_mw'])
                    order_dict_gen = collections.OrderedDict(sorted(
                        intermediate_gen_power.items()))
                    intermediate_gen_power = list(order_dict_gen.values())
                    intermediate_load_power = dict(net.res_load.loc[:, 'p_mw'])
                    order_dict_load = collections.OrderedDict(sorted(
                        intermediate_load_power.items()))
                    intermediate_load_power = list(order_dict_load.values())
                    cost = 0
                except pp.OPFNotConverged:
                    pass
                if load_cost == 0:
                    net.gen.loc[:, 'p_mw'] = intermediate_gen_power
                    net.load.loc[:, 'p_mw'] = intermediate_load_power
                    pp.rundcpp(net)
                    break
        except UserWarning:
            total_blackout = 1
            break
        if total_blackout == 1:
            break
        list_line_in_service = list(net.line.loc[net.line.in_service == 
            True].index)
        list_line_in_service.sort()
        new_failure = 0
        for line in list_line_in_service:
            level_loading = net.res_line.loc[line, 'loading_percent']
            if level_loading >= 99:
                print('OVERLOADS')
                new_failure = 1
                net.line.loc[line, 'in_service'] = False
        list_trafo_in_service = list(net.trafo.loc[net.trafo.in_service == 
            True].index)
        list_trafo_in_service.sort()
        for trafo in list_trafo_in_service:
            level_loading = net.res_trafo.loc[trafo, 'loading_percent']
            if level_loading >= 99:
                print('OVERLOADS')
                new_failure = 1
                net.trafo.loc[trafo, 'in_service'] = False
        if new_failure == 0:
            break
    if total_blackout == 0:
        final_state_load = dict(net.res_load.loc[:, 'p_mw'])
        order_dict_load = collections.OrderedDict(sorted(final_state_load.
            items()))
        list_load_power = list(order_dict_load.values())
        final_state_gen = dict(net.res_gen.loc[:, 'p_mw'])
        order_dict_gen = collections.OrderedDict(sorted(final_state_gen.
            items()))
        list_gen_power = list(order_dict_gen.values())
        list_gen_power = np.array(list_gen_power)
        list_load_power = np.array(list_load_power)
        for i, j in enumerate(list_gen_power):
            if math.isnan(j):
                list_gen_power[i] = 0
        for i, j in enumerate(list_load_power):
            if math.isnan(j):
                list_load_power[i] = 0
    elif total_blackout == 1:
        list_gen_power = np.zeros(len(net.gen.index))
        list_load_power = np.zeros(len(net.load.index))
    return list_gen_power, list_load_power


<mask token>

# -*- coding: utf-8 -*-
"""
Created on Wed Dec  2 09:53:35 2020

@author: Andrea Bellè
"""

import numpy as np
import collections
import os
import pandapower as pp
import pandapower.networks
import networkx as nx
import time
import math
import sys
#import matplotlib as mpl
#import matplotlib.pyplot as plt
#import seaborn as sb
#import multiprocessing  
#import pandapower.plotting as pplt
#from mpl_toolkits.axes_grid1 import make_axes_locatable
#from sklearn.linear_model import LinearRegression
#from sklearn.metrics import mean_squared_error, r2_score
#from matplotlib.patches import Patch
#from matplotlib.lines import Line2D
#from glob import glob
#from pylab import *
#from matplotlib import cm
#from mpl_toolkits.mplot3d import Axes3D
#from mpl_toolkits.basemap import Basemap
#from mpi4py import MPI
#import random
#import itertools
'''
import conda
conda_file_dir = conda.__file__
conda_dir = conda_file_dir.split('lib')[0]
proj_lib = os.path.join(os.path.join(conda_dir, 'share'), 'proj')
os.environ["PROJ_LIB"] = proj_lib
'''

##############################################################################
#time_start = time.process_time()
##############################################################################

# INITIALIZATION OF THE POWER NETWORK

# Call the Great Britain reduced power network from Pandapower
net = pandapower.networks.GBreducednetwork()

# Replace static generators with standard generators
pp.replace_sgen_by_gen(net)
list_ext_grid  = list(net.ext_grid.index)
initial_slack = pp.replace_ext_grid_by_gen(net, list_ext_grid, slack=True) 
    
# Set the maximum loading of each line to 100%   
net.line.loc[:,'max_loading_percent'] = 100
net.trafo.loc[:,'max_loading_percent'] = 100

# Set the min/max power values of each load
net.load.loc[:,'max_p_mw'] = net.load.p_mw
net.load.loc[:,'min_p_mw'] = 0

# Set the min/max power values of each generator
net.gen.loc[:,'max_p_mw'] = net.gen.max_p_mw/1.438
net.gen.loc[:,'p_mw'] = net.gen.max_p_mw

# Set generators and loads as controllable (for DC OPF)
net.gen.loc[:,'controllable'] = True
net.load.loc[:,'controllable'] = True

# Set cost function for each generator
net.poly_cost.loc[:,'cp0_eur'] = 0
net.poly_cost.loc[:,'cp1_eur_per_mw'] = 1
net.poly_cost.loc[:,'cp2_eur_per_mw2'] = 0

# Set cost function for each load
list_load =list(net.load.index)
for load in list_load:
    pp.create_poly_cost(net, load, 'load',cp0_eur=0,
                        cp1_eur_per_mw = -100,cp2_eur_per_mw2=0)

# Compute the initial power of each load to determine the load shedding
initial_state_load = dict(net.load.loc[:,'p_mw'])
order_dict_load = collections.OrderedDict(sorted(initial_state_load.items()))
initial_load_power = list(order_dict_load.values())


# CASCADING FAILURES FUNCTION
'''
INPUT: 
      - net: initial power network configuration in Pandapower format
      - f_bus: indexes of initial failed buses 
      - f_gen: indexes of initial failed generators
      - f_load: indexes of initial failed loads
      - f_line: indexes of initial failed lines
      - f_trafo: indexes of initial failed transformers 
OUTPUT: 
      - final power production at each generator
      - final power demand at each load 
'''
def OPA_model(net, f_bus, f_gen, f_load, f_line, f_trafo):  
    
    # This variable allows to stop the function when a total blackout situation is reached
    total_blackout = 0
    
    # Set failed elements as out of service.
    # In the case of failed buses, all the elements connected to the bus are also removed
    f_bus = list(f_bus)
    pp.set_element_status(net, f_bus, in_service=False)
    f_gen = list(f_gen)
    net.gen.loc[f_gen,'in_service'] = False
    f_load = list(f_load)
    net.load.loc[f_load,'in_service'] = False
    f_line = list(f_line)
    net.line.loc[f_line,'in_service'] = False
    f_trafo = list(f_trafo)
    net.trafo.loc[f_trafo,'in_service'] = False
    
    # Variable for entering in the while loop
    cascade = True
    
    # Store the initial power level of each generator
    
    intermediate_gen_power = dict(net.gen.loc[:,'p_mw'])
    order_dict_gen = collections.OrderedDict(sorted(intermediate_gen_power.items()))
    intermediate_gen_power = list(order_dict_gen.values())
    
    intermediate_load_power = dict(net.load.loc[:,'p_mw'])
    order_dict_load = collections.OrderedDict(sorted(intermediate_load_power.items()))
    intermediate_load_power = list(order_dict_load.values())
    
    while cascade == True:
        
        # Create a networkx graph of the power network
        G=pp.topology.create_nxgraph(net)
        # Create a list of sets containing the islands of the network
        list_islands = list(nx.connected_components(G))  
        
        # Create a list of buses with generators in service
        list_bus_with_gen_in_service = list(net.gen.bus[net.gen.in_service==True])
        list_bus_with_gen_in_service = list(set(list_bus_with_gen_in_service))
        list_bus_with_gen_in_service.sort()
        list_bus_with_slack_gen_in_service = list(net.gen.bus[net.gen.in_service==True][net.gen.slack==True])

        # Check the configuration of each island in the power network
        for island in list_islands:
            #print(island)
            island = list(island)
            
            # Check if the island has an external grid. If yes, it is ok and you can proceed to the next island.
            # If no, check if there is an availabl generator in the island. If yes, turn it into an external grid.
            # If no, go to the next island.
            if set(list_bus_with_slack_gen_in_service).isdisjoint(island):
                if set(list_bus_with_gen_in_service).isdisjoint(island):
                    pass
                else:
                    set_island = set(island)
                    bus_with_gen_in_island = list(set_island.intersection(list_bus_with_gen_in_service))
                    slack_gen = list(net.gen.index[net.gen.bus==bus_with_gen_in_island[0]])
                    net.gen.loc[slack_gen[0],'slack'] = True
            else:
                pass
            
        # Set the isolated islands (without ext. grid/slack gen.) out of service
        pp.set_isolated_areas_out_of_service(net)   
        
        # Try to run a DC OPF
        try:
            pp.rundcopp(net, verbose=False)
            # If it converges, save a the power level at each generators and ext grids
            intermediate_gen_power = dict(net.res_gen.loc[:,'p_mw'])
            order_dict_gen = collections.OrderedDict(sorted(intermediate_gen_power.items()))
            intermediate_gen_power = list(order_dict_gen.values())
            
            intermediate_load_power = dict(net.res_load.loc[:,'p_mw'])
            order_dict_load = collections.OrderedDict(sorted(intermediate_load_power.items()))
            intermediate_load_power = list(order_dict_load.values())
            
        # If the DC OPF does not converge, increase the loads costs 
        # A while loop which decreases the loads costs until the DC OPF converges or the costs reach a value of 0
        except pp.OPFNotConverged:
            #print(f_line, f_bus, f_trafo)
            print('CONVERGENCE PROBLEMS')
            cost = 1
            load_cost = -100
            while cost == 1:
                load_cost += 2
                #print(load_cost)
                net.poly_cost.cp1_eur_per_mw[net.poly_cost.et == 'load'] = load_cost
                try:
                    pp.rundcopp(net)
                    #pp.rundcopp(net, SCPDIPM_RED_IT=100, PDIPM_COSTTOL = 1e-3, PDIPM_GRADTOL = 1e-3)
                    intermediate_gen_power = dict(net.res_gen.loc[:,'p_mw'])
                    order_dict_gen = collections.OrderedDict(sorted(intermediate_gen_power.items()))
                    intermediate_gen_power = list(order_dict_gen.values())

                    intermediate_load_power = dict(net.res_load.loc[:,'p_mw'])
                    order_dict_load = collections.OrderedDict(sorted(intermediate_load_power.items()))
                    intermediate_load_power = list(order_dict_load.values())

                    cost = 0
                except pp.OPFNotConverged:
                    pass
                # If the loads costs are set to 0, run a DC power flow with the last available generators power levels
                if load_cost == 0:
                    net.gen.loc[:,'p_mw'] = intermediate_gen_power
                    net.load.loc[:,'p_mw'] = intermediate_load_power
                    pp.rundcpp(net)
                    break
                
        # This is in case to generators is available (everything is failed basically)
        except UserWarning:
            total_blackout = 1
            break
        # If there is a total blackout situation, break the loop
        if total_blackout == 1:
            break
        
        # Create a list of lines still in service
        list_line_in_service = list(net.line.loc[net.line.in_service== True].index)
        list_line_in_service.sort()
        # This variable is to check if there are new failures due to overloads
        new_failure = 0
        for line in list_line_in_service:
            level_loading = net.res_line.loc[line, 'loading_percent']
            # If a line a loading >= the 99% of its maximum capacity, it is considered overloaded
            if level_loading >= 99:
                print('OVERLOADS')
                new_failure = 1
                net.line.loc[line, 'in_service'] = False
        
        # Same for transformer lines
        list_trafo_in_service = list(net.trafo.loc[net.trafo.in_service== True].index)
        list_trafo_in_service.sort()
        for trafo in list_trafo_in_service:
            level_loading = net.res_trafo.loc[trafo, 'loading_percent']
            if level_loading >= 99:
                print('OVERLOADS')
                new_failure = 1
                net.trafo.loc[trafo, 'in_service'] = False
        
        # If there are no overloads, break the while loop
        if new_failure == 0:
            break
    
    # Save the final power levels of each generator and load
    if total_blackout == 0:        
        final_state_load = dict(net.res_load.loc[:,'p_mw'])
        order_dict_load = collections.OrderedDict(sorted(final_state_load.items()))
        list_load_power = list(order_dict_load.values())
        
        final_state_gen = dict(net.res_gen.loc[:,'p_mw'])
        order_dict_gen = collections.OrderedDict(sorted(final_state_gen.items()))
        list_gen_power = list(order_dict_gen.values())
        
        list_gen_power = np.array(list_gen_power)
        list_load_power = np.array(list_load_power)       
        
        for i,j in enumerate(list_gen_power):
            if math.isnan(j):
                list_gen_power[i] = 0
        for i,j in enumerate(list_load_power):
            if math.isnan(j):
                list_load_power[i] = 0
    elif total_blackout == 1:
        list_gen_power = np.zeros(len(net.gen.index))
        list_load_power = np.zeros(len(net.load.index))
        
    return(list_gen_power,list_load_power)



'''
g,l = OPA_model(net, [8,9], [], [], [], [])   
print(np.around(l/initial_load_power, 2))
'''

