dataset_generator_pandapower.py

import time
import pandapower as pp
import numpy as np
import pickle
# write file documentation here




# dict_keys(['bus', 'load', 'sgen', 'motor', 'asymmetric_load', 'asymmetric_sgen', 'storage', 'gen', 'switch', 'shunt', 'svc', 'ext_grid', 'line', 'trafo', 'trafo3w', 'impedance', 'tcsc', 'dcline', 'ward', 'xward', 'measurement', 'pwl_cost', 'poly_cost', 'characteristic', 'controller', 'group', 'line_geodata', 'bus_geodata', '_empty_res_bus', '_empty_res_ext_grid', '_empty_res_line', '_empty_res_trafo', '_empty_res_load', '_empty_res_asymmetric_load', '_empty_res_asymmetric_sgen', '_empty_res_motor', '_empty_res_sgen', '_empty_res_shunt', '_empty_res_svc', '_empty_res_switch', '_empty_res_impedance', '_empty_res_tcsc', '_empty_res_dcline', '_empty_res_ward', '_empty_res_xward', '_empty_res_trafo_3ph', '_empty_res_trafo3w', '_empty_res_bus_3ph', '_empty_res_ext_grid_3ph', '_empty_res_line_3ph', '_empty_res_asymmetric_load_3ph', '_empty_res_asymmetric_sgen_3ph', '_empty_res_storage', '_empty_res_storage_3ph', '_empty_res_gen',
        #   '_ppc', '_ppc0', '_ppc1', '_ppc2', '_is_elements', '_pd2ppc_lookups', 'version', 'format_version', 'converged', 'OPF_converged', 'name', 'f_hz', 'sn_mva', '_empty_res_load_3ph', '_empty_res_sgen_3ph', 'std_types', 'res_bus', 'res_line', 'res_trafo', 'res_trafo3w', 'res_impedance', 'res_ext_grid', 'res_load', 'res_motor', 'res_sgen', 'res_storage', 'res_shunt', 'res_gen', 'res_ward', 'res_xward', 'res_dcline', 'res_asymmetric_load', 'res_asymmetric_sgen', 'res_switch', 'res_tcsc', 'res_svc', 'res_bus_est', 'res_line_est', 'res_trafo_est', 'res_trafo3w_est', 'res_impedance_est', 'res_switch_est', 'res_bus_sc', 'res_line_sc', 'res_trafo_sc', 'res_trafo3w_sc', 'res_ext_grid_sc', 'res_gen_sc', 'res_sgen_sc', 'res_switch_sc', 'res_bus_3ph', 'res_line_3ph', 'res_trafo_3ph', 'res_ext_grid_3ph', 'res_shunt_3ph', 'res_load_3ph', 'res_sgen_3ph', 'res_storage_3ph', 'res_asymmetric_load_3ph', 'res_asymmetric_sgen_3ph', 'user_pf_options'])

# net = pp.networks.GBnetwork()

# algorithm (str, “nr”) - algorithm that is used to solve the power flow problem.

# The following algorithms are available:

# “nr” Newton-Raphson (pypower implementation with numba accelerations)

# “iwamoto_nr” Newton-Raphson with Iwamoto multiplier (maybe slower than NR but more robust)

# “bfsw” backward/forward sweep (specially suited for radial and weakly-meshed networks)

# “gs” gauss-seidel (pypower implementation)

# “fdbx” fast-decoupled (pypower implementation)

# “fdxb” fast-decoupled (pypower implementation)

# print(net)
# print(net.keys())

number_of_samples = 100000

test_case = 'case118'
# base_net = pp.networks.case6470rte()
base_net = pp.networks.case118()
base_net.bus['name'] = base_net.bus.index
print(base_net.bus)
print(base_net.line)
print(base_net.gen)

# Get Adjacency Matrix
bus_names = base_net.bus['name'].values.tolist()
n = base_net.bus.values.shape[0]
A = np.zeros((n, n))
for edge1,edge2 in base_net.line[['from_bus', 'to_bus']].values:
    
    edge_1 = bus_names.index(edge1)
    edge_2 = bus_names.index(edge2)

    A[edge_1, edge_2] = 1
    A[edge_2, edge_1] = 1

edge_features_list = []
node_features_x_list = []
node_features_y_list = []
graph_feature_list = []
reconstruction_case_list = []

ref_bus = base_net.ext_grid['bus'].values[0]
print(f'ref_bus: {ref_bus}')

r_original = base_net.line['r_ohm_per_km'].values
x_original = base_net.line['x_ohm_per_km'].values

vm_original = base_net.gen['vm_pu'].values
vg_original = base_net.gen['p_mw'].values

pd_original = base_net.load['p_mw'].values
qd_original = base_net.load['q_mvar'].values
net = base_net

counter = 0
while True:

    counter += 1
    if counter % 100 == 0:
        base_net = pp.networks.case118()    
    else:
        net = base_net
    # net = pp.networks.case6470rte()
    # net = pp.networks.case14()
    net.bus['name'] = base_net.bus.index

    # net.line['r_ohm_per_km'] = r_original
    # net.line['x_ohm_per_km'] = x_original

    # net.gen['vm_pu'] = vm_original
    # net.gen['p_mw'] = vg_original
    
    # net.load['p_mw'] = pd_original
    # net.load['q_mvar'] = qd_original



    r = net.line['r_ohm_per_km'].values    
    x = net.line['x_ohm_per_km'].values
    # c = net.line['c_nf_per_km'].values
    le = net.line['length_km'].values
    # x = case['branch'][:, 3]
    # b = case['branch'][:, 4]
    # tau = case['branch'][:, 8]  # ratio

    Pmax = net.gen['max_p_mw'].values
    Pmin = net.gen['min_p_mw'].values
    Pmw = net.gen['p_mw'].values

    Pd = net.load['p_mw'].values
    Qd = net.load['q_mvar'].values

    r = np.random.uniform(0.8*r, 1.2*r, r.shape[0])    
    x = np.random.uniform(0.8*x, 1.2*x, x.shape[0])
    # le = np.random.uniform(0.8*le, 1.2*le, le.shape[0])
   
    vg = np.random.uniform(0.95, 1.05, net.gen['vm_pu'].shape[0])
    # Pg = np.random.uniform(0.25*Pmax, 0.75*Pmax, net.gen['p_mw'].shape[0])
    Pg = np.random.uniform(0.25*Pmw, 1.75*Pmw, net.gen['p_mw'].shape[0])
    
    
    Pd = np.random.uniform(0.5*Pd, 1.5*Pd, net.load['p_mw'].shape[0])
    Qd = np.random.uniform(0.5*Qd, 1.5*Qd, net.load['q_mvar'].shape[0])
    
    net.line['r_ohm_per_km'] = r
    net.line['x_ohm_per_km'] = x

    net.gen['vm_pu'] = vg
    net.gen['p_mw'] = Pg

    net.load['p_mw'] = Pd
    net.load['q_mvar'] = Qd

    

    try:    
        pp.runpp(net, algorithm='nr', init="results", numba=False)
    except:
        print(f'Failed to converge, current sample number: {len(edge_features_list)}')        
        continue        
    
    reconstruction_case_list.append((r,x,vg,Pg,Pd,Qd))

    # Graph feature
    # baseMVA = x[0]['baseMVA']

    # Create a vector od branch features including start and end nodes,r,x,b,tau,angle
    edge_features = np.zeros((net.line.shape[0], 7))
    edge_features[:, 0] = net.line['from_bus'].values + 1
    edge_features[:, 1] = net.line['to_bus'].values + 1
    edge_features[:, 2] = net.line['r_ohm_per_km'].values * net.line['length_km'].values
    edge_features[:, 3] = net.line['x_ohm_per_km'].values * net.line['length_km'].values
    edge_features[:, 4] = 0
    edge_features[:, 5] = 0
    edge_features[:, 6] = 0

    # Create a vector of node features including index, type, Vm, Va, Pd, Qd, Gs, Bs, Pg
    # case['bus'] = x[0]['bus']

    node_features_x = np.zeros((n, 9))
    node_features_x[:, 0] = net.bus['name'].values + 1# index
    # Va ----This changes for every bus excecpt slack bus
    node_features_x[:, 3] = np.zeros((n, )) #Va

    node_features_x[ref_bus, 3] = 0 #Va for reference bus
    
    # node_features_x[:, 6] = np.zeros((n,1)) # Gs
    # node_features_x[:, 7] = np.zeros((n,1)) # Bs
    # Vm is 1 if type is not "generator" else it is case['gen'][:,j]
    vm = np.ones(n)
    types = np.ones(n)*2
    for j in range(net.gen.shape[0]):    
        # find index of case['gen'][j,0] in case['bus'][:,0]
        index = np.where(net.gen['bus'].values[j] == net.bus['name'])[0][0]        
        vm[index] = net.gen['vm_pu'].values[j]  # Vm = Vg
        types[index] = 1  # type = generator
        node_features_x[index, 8] = net.gen['p_mw'].values[j]  # Pg
    
    types[ref_bus] = 3 # 
    node_features_x[:, 2] = vm  # Vm
    node_features_x[:, 1] = types  # type
    
    for j in range(net.load.shape[0]):    
        # find index of case['gen'][j,0] in case['bus'][:,0]
        index = np.where(net.load['bus'].values[j] == net.bus['name'])[0][0]        
        node_features_x[index, 4] = Pd[j]  # Pd
        node_features_x[index, 5] = Qd[j] # Qd

    # Create a vector of node features including index, type, Vm, Va, Pd, Qd, Gs, Bs    
    node_features_y = np.zeros((n, 8))
    node_features_y[:, 0] = net.bus['name'].values + 1 # index
    node_features_y[:, 1] = types  # type
    # Vm ----This changes for Load Buses
    node_features_y[:, 2] = net.res_bus['vm_pu']  # Vm
    # Va ----This changes for every bus excecpt slack bus
    node_features_y[:, 3] = net.res_bus['va_degree']  # Va
    node_features_y[:, 4] = net.res_bus['p_mw']  # P
    node_features_y[:, 5] = net.res_bus['q_mvar']  # Q
    # node_features_y[:, 6] = case['bus'][:, 4]  # Gs
    # node_features_y[:, 7] = case['bus'][:, 5]  # Bs

    edge_features_list.append(edge_features)
    node_features_x_list.append(node_features_x)
    node_features_y_list.append(node_features_y)
    # graph_feature_list.append(baseMVA)

    if len(edge_features_list) == number_of_samples:
        break
    elif len(edge_features_list) % 100 == 0:
        print(f'Current sample number: {len(edge_features_list)}')

    if len(edge_features_list) % 1000 == 0:
        # Turn the lists into numpy arrays
        edge_features = np.array(edge_features_list)
        node_features_x = np.array(node_features_x_list)
        node_features_y = np.array(node_features_y_list)
        # graph_features = np.array(graph_feature_list)
        # reconstruction_case = np.array(reconstruction_case_list)

        with open("./data/raw/"+test_case+"_edge_features.npy", 'wb') as f:
            np.save(f, edge_features)

        with open("./data/raw/"+test_case+"_node_features_x.npy", 'wb') as f:
            np.save(f, node_features_x)

        with open("./data/raw/"+test_case+"_node_features_y.npy", 'wb') as f:
            np.save(f, node_features_y)

        #pickle the reconstructed cases
        with open("./data/raw/"+test_case+"_reconstruction_case.pkl", 'wb') as f:
            pickle.dump(reconstruction_case_list, f)
        

        # with open("./data/"+test_case+"_graph_features.npy", 'wb') as f:
        #     np.save(f, graph_features)

        with open("./data/raw/"+test_case+"_adjacency_matrix.npy", 'wb') as f:
            np.save(f, A)        

edge_features = np.array(edge_features_list)
node_features_x = np.array(node_features_x_list)
node_features_y = np.array(node_features_y_list)
# graph_features = np.array(graph_feature_list)

# Print the shapes
print(f'Adjacency matrix shape: {A.shape}')
print(f'edge_features shape: {edge_features.shape}')
print(f'node_features_x shape: {node_features_x.shape}')
print(f'node_features_y shape: {node_features_y.shape}')
# print(f'graph_features shape: {graph_features.shape}')

print(f'range of edge_features "from": {np.min(edge_features[:,:,0])} - {np.max(edge_features[:,:,0])}')
print(f'range of edge_features "to": {np.min(edge_features[:,:,1])} - {np.max(edge_features[:,:,1])}')

print(f'range of node_features_x "index": {np.min(node_features_x[:,:,0])} - {np.max(node_features_x[:,:,0])}')

print(f'range of node_features_y "index": {np.min(node_features_y[:,:,0])} - {np.max(node_features_y[:,:,0])}')