import time import pandas as pd import pandapower as pp import numpy as np import networkx as nx import multiprocessing as mp import os # 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()) def create_case3(): net = pp.create_empty_network() net.sn_mva = 100 b0 = pp.create_bus(net, vn_kv=345., name='bus 0') b1 = pp.create_bus(net, vn_kv=345., name='bus 1') b2 = pp.create_bus(net, vn_kv=345., name='bus 2') pp.create_ext_grid(net, bus=b0, vm_pu=1.02, name="Grid Connection") pp.create_load(net, bus=b2, p_mw=10.3, q_mvar=3, name="Load") # pp.create_gen(net, bus=b1, p_mw=0.5, vm_pu=1.03, name="Gen", max_p_mw=1) pp.create_line(net, from_bus=b0, to_bus=b1, length_km=10, name='line 01', std_type='NAYY 4x50 SE') pp.create_line(net, from_bus=b1, to_bus=b2, length_km=5, name='line 01', std_type='NAYY 4x50 SE') pp.create_line(net, from_bus=b2, to_bus=b0, length_km=20, name='line 01', std_type='NAYY 4x50 SE') net.line['c_nf_per_km'] = pd.Series(0., index=net.line['c_nf_per_km'].index, name=net.line['c_nf_per_km'].name) return net def remove_c_nf(net): net.line['c_nf_per_km'] = pd.Series(0., index=net.line['c_nf_per_km'].index, name=net.line['c_nf_per_km'].name) def unify_vn(net): for node_id in range(net.bus['vn_kv'].shape[0]): net.bus['vn_kv'][node_id] = max(net.bus['vn_kv']) def get_trafo_z_pu(net): # if random: # replace_line_r_ohm = max(0., 0.1*np.random.normal(replace_line_r_ohm, np.abs(replace_line_r_ohm)*0.05)) # replace_line_x_ohm = max(0., 0.1*np.random.normal(replace_line_x_ohm, np.abs(replace_line_x_ohm)*0.05)) # for trafo_id in net.trafo.index: # from_bus, to_bus = net.trafo.iloc[trafo_id]['hv_bus'], net.trafo.iloc[trafo_id]['lv_bus'] # pp.create_line_from_parameters(net, from_bus=from_bus, to_bus=to_bus, length_km=1., # r_ohm_per_km=replace_line_r_ohm, x_ohm_per_km=replace_line_r_ohm, # c_nf_per_km=0., max_i_ka=3.e4, # max_loading_percent=100., # type='ol') # pp.drop_trafos(net, net.trafo.index, table='trafo') for trafo_id in net.trafo.index: net.trafo['i0_percent'][trafo_id] = 0. net.trafo['pfe_kw'][trafo_id] = 0. z_pu = net.trafo['vk_percent'].values / 100. * 1000. / net.sn_mva r_pu = net.trafo['vkr_percent'].values / 100. * 1000. / net.sn_mva x_pu = np.sqrt(z_pu**2 - r_pu**2) return x_pu, r_pu # raise NotImplementedError def get_line_z_pu(net): r = net.line['r_ohm_per_km'].values * net.line['length_km'] x = net.line['x_ohm_per_km'].values * net.line['length_km'] from_bus = net.line['from_bus'] to_bus = net.line['to_bus'] vn_kv_to = net.bus['vn_kv'][to_bus].to_numpy() vn_kv_to = pd.Series(vn_kv_to) zn = vn_kv_to**2 / net.sn_mva r_pu = r/zn x_pu = x/zn return r_pu, x_pu number_of_samples = 30000 number_of_processes = 1 test_case = 'case14v3' base_net_create = pp.networks.case14 # base_net_create = create_case3 base_net = base_net_create() base_net.bus['name'] = base_net.bus.index print(base_net.bus) print(base_net.line) # 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 multi_graph = pp.topology.create_nxgraph(base_net) A = nx.adjacency_matrix(multi_graph).todense() def generate_data(sublist_size): edge_features_list = [] node_features_x_list = [] node_features_y_list = [] # graph_feature_list = [] while True: # net = base_net net = base_net_create() # remove_c_nf(net) # unify_vn(net) trafo_r_pu, trafo_x_pu = get_trafo_z_pu(net) net.bus['name'] = base_net.bus.index 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 Pg = net.gen['p_mw'].values # Pmin = 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]) # c = np.random.uniform(0.8*c, 1.2*c, c.shape[0]) le = np.random.uniform(0.8*le, 1.2*le, le.shape[0]) # tau = np.random.uniform(0.8*tau, 1.2*tau, case['branch'].shape[0]) # angle = np.random.uniform(-0.2, 0.2, case['branch'].shape[0]) Vg = np.random.uniform(1.00, 1.05, net.gen['vm_pu'].shape[0]) Pg = np.random.normal(Pg, 0.1*np.abs(Pg), net.gen['p_mw'].shape[0]) # Pd = np.random.uniform(0.5*Pd, 1.5*Pd, net.load['p_mw'].shape[0]) Pd = np.random.normal(Pd, 0.1*np.abs(Pd), net.load['p_mw'].shape[0]) # Qd = np.random.uniform(0.5*Qd, 1.5*Qd, net.load['q_mvar'].shape[0]) Qd = np.random.normal(Qd, 0.1*np.abs(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: net['converged'] = False pp.runpp(net, algorithm='nr', init="results", numba=False) except: if not net['converged']: # print(f"net['converged'] = {net['converged']}") print(f'Failed to converge, current sample number: {len(edge_features_list)}') import pandapower as pp continue # Graph feature # baseMVA = x[0]['baseMVA'] # === DEBUG === ybus_matrix = net._ppc["internal"]["Ybus"] pandapower_bus_idx = 3 ppc_index = net._pd2ppc_lookups["bus"][pandapower_bus_idx] # === DEBUG === # 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], edge_features[:, 3] = get_line_z_pu(net) edge_features[:, 4] = 0 edge_features[:, 5] = 0 edge_features[:, 6] = 0 trafo_edge_features = np.zeros((net.trafo.shape[0], 7)) trafo_edge_features[:, 0] = net.trafo['hv_bus'].values + 1 trafo_edge_features[:, 1] = net.trafo['lv_bus'].values + 1 trafo_edge_features[:, 2], trafo_edge_features[:, 3] = get_trafo_z_pu(net) trafo_edge_features[:, 4] = 0 trafo_edge_features[:, 5] = 0 trafo_edge_features[:, 6] = 0 edge_features = np.concatenate((edge_features, trafo_edge_features), axis=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[:, 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] / net.sn_mva # Pg / pu 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] / net.sn_mva # Pd / pu node_features_x[index, 5] = Qd[j] / net.sn_mva # Qd / pu # 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 # if net.res_bus['vm_pu'].shape[0] == 0: # pass 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'] / net.sn_mva # P / pu node_features_y[:, 5] = net.res_bus['q_mvar'] / net.sn_mva # Q / pu # 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) == sublist_size: break elif len(edge_features_list) % 10 == 0: print(f'[Process {os.getpid()}] Current sample number: {len(edge_features_list)}') return edge_features_list, node_features_x_list, node_features_y_list def generate_data_parallel(num_samples, num_processes): sublist_size = num_samples // num_processes pool = mp.Pool(processes=num_processes) results = pool.map(generate_data, [sublist_size]*num_processes) pool.close() pool.join() edge_features_list = [] node_features_x_list = [] node_features_y_list = [] for sub_res in results: edge_features_list += sub_res[0] node_features_x_list += sub_res[1] node_features_y_list += sub_res[2] return edge_features_list, node_features_x_list, node_features_y_list if __name__ == '__main__': # Generate data edge_features_list, node_features_x_list, node_features_y_list = generate_data_parallel(number_of_samples, number_of_processes) # 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) # 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])}') # print(f"A. {A}") # print(f"edge_features. {edge_features}") # print(f"node_features_x. {node_features_x}") # print(f"node_features_y. {node_features_y}") # save the features os.makedirs("./data/raw", exist_ok=True) 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) # 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) exit() # Computation time experimental comparison beginning (will be moved to other file later on) # calculate power flow for every algorithm and calculate time algorithms = ["nr", "iwamoto_nr", "gs", "fdbx", "fdxb"] times = [] for a in algorithms: t0 = time.time() # pp.runpp(net, algorithm=a) pp.runpp(net, algorithm=a, init="results", numba=False) t1 = time.time() times.append(t1 - t0) for a in algorithms: print(f"{a}: {times[algorithms.index(a)]}") # print(net.res_bus.vm_pu) # print(net.res_line.loading_percent) # calculate power flow for every algorithm and calculate time 1000 times # algorithms = ["nr", "iwamoto_nr", "gs", "fdbx", "fdxb"] algorithms = ["nr", "iwamoto_nr", "fdbx", "fdxb"] times = [] for a in algorithms: print(a) t0 = time.time() for i in range(1000): pp.runpp(net, algorithm=a, init="auto", numba=False) t1 = time.time() times.append(t1 - t0) for a in algorithms: print(f"{a}: {times[algorithms.index(a)]/1000}")