From 3d2d6e252f504611e44bfb9127537b72732884a9 Mon Sep 17 00:00:00 2001 From: Thomas Date: Mon, 8 Aug 2016 18:25:08 -0400 Subject: [PATCH] update python files power_flow, simulation_annealing, and percolation --- py/__init__.py | 0 py/__pycache__/power_flow.cpython-35.pyc | Bin 0 -> 2702 bytes py/percolation_threshold_simulation.py | 19 ++--- py/power_flow.py | 66 ++++++---------- py/simulated_annealing_algo.py | 96 +++++++++++++++++------ 5 files changed, 108 insertions(+), 73 deletions(-) create mode 100644 py/__init__.py create mode 100644 py/__pycache__/power_flow.cpython-35.pyc diff --git a/py/__init__.py b/py/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/py/__pycache__/power_flow.cpython-35.pyc b/py/__pycache__/power_flow.cpython-35.pyc new file mode 100644 index 0000000000000000000000000000000000000000..d05ede7a1f74aa3cd25e88910574038ae668f0bd GIT binary patch literal 2702 zcmb_d&2Jk;6o0e3{)l6@B_E}sMQkaB7^#T_LPa40)I{afs3Jv?tX9^>yS6u7?>e*N z#;u$S?TJ5v8;71b@LzD@HaEnD9=URX-+LQ7X#sKRCNrLW``*0o-@Djo)PDA!?|ok* z`i;&!C7i#+$_=&<|B4c#mXfnWi9-XIS}qMdYI&qiI7-wiQ|MCYQCOl7=L$UoUqxQ4 z6wXmtqh}7aaK`^Ug>`zSs0Dh1!gKTtI%=$Jw|V{rG}BJgv5URPw9_n%+JP}T-XEDL zw~IN5QP_Kl=p|VYTDKRawtN)Dy}s$NBm6Y5%2-=i`5)Ne*o5e*qT>=7MF$Q|V6vjF zBg8Tpml-B_b?CT4#-mAvDFaf z2G379iH_?u!9TH6yG<|n)gh|+cy7mO=$n4z53{34``siv^1GS#qoC9GIjVk^`hnlm z!LYyP`(H<9q*DVWSk zHhsCKnr5QhL?(})L>qpbn%NCIqrpBxm34hC#T?(n5%E}4u*h zeJTdqS%`yucWX<83l&_eZLkWzznTqAoTWk1{8!OglyFgCpdr(+;WUGM8LN>5vVl{b zVqc*UZSmUl;jTyrZ1MkT4ddJdX(xIOs*@<~p?;=r7Ij*L_Kce^6uR4bVPFF7;n=KO z7n0UVhqm-2(pj$224qV^jr`h@q|WG#0XcEnD7j;&I+tZLiFqbxB2 zz-`zM^2j!(T_?((lfDerKX=!^$PYHPot8Q%iWmpQGQ#x8BF{ z{PIP|WmY_MSQRXR+MjSZg<5wCwZK+)3bh`Mj|6Q2qYB^|3nY6Cv<$!{0pK!?cNprn zAKYd%W&|z(wFGFbU@0tzmM#M+A;PrH>9GVzou0uuf79JOyhD43YrtM0e}xW!u>e@i z0yzM(1u)$m&3?T23r}-SvtJs6lJq<_Thc+=i%?c$*1cw_fCPY7@^9>+hDk?|<$ z`GJmQI36I-NDt-|1S&NCt{8paN)QppGClN2-{XlUGY~&z9(*!E`N&j^&XcsgOjgJq ziUsrVRr(5yxCt{kq5t@H6P6#|Lux;vJvL5i43lridfcX(p+R#Dnib9ETIBQ{KKZ;| zK%2}(abl?6&uIKN)IuLA&R?AwIi1|-J+nlTv}CjjKTF4z4p>+C$}HP`6x)@ zC({3SAeJ)(lG07=&3Ri*PbRj6!LCeS*E&g%=TW}a9i<({4uMVy9e{H*7(Ol&zK`cdj%}Ex;e 0: - PDC[PCC_node] += grid_compensation + pdc[pcc_node] -= sum(pdc) - else: - PDC[PCC_node] -= grid_compensation - - B = find_B_matrix(graph, size, distance_str) - voltage_phase = np.linalg.solve(B, PDC) - absolute_flow = 0 - graph, absolute_flow = calculate_power_flow(graph, size, voltage_phase) + B, _ = find_b_matrix(graph, distance_str) + voltage_phase = np.linalg.solve(B, pdc) + graph, absolute_flow = calculate_power_flow(graph, voltage_phase) return absolute_flow, graph -def find_grid_compensation(graph, load_str = 'load', gen_str = 'gen'): - total_load = 0 - total_gen = 0 - grid_compensation = 0 - load = [] - gen = [] - load = get_node_attributes(graph, load_str) - gen = get_node_attributes(graph, gen_str) - print load - print gen - for node in load: - total_load += graph.node[node][load_str] - total_gen += graph.node[node][gen_str] - - grid_compensation = total_load - total_gen - - return grid_compensation - -def find_B_matrix(graph, size, distance_str = 'weight', reactance_per_thousand_feet = 0.3): +def find_b_matrix(graph, distance_str='weight', reactance_per_thousand_feet=0.3): + size = len(graph) B = np.zeros((size, size)) for i in range(size): - neighbors = nx.all_neighbors(graph, i) for node in neighbors: - B[i][node] = -1/(graph.edge[i][node][distance_str] * reactance_per_thousand_feet) + B[i][node] = -1 / (graph.edge[i][node][distance_str] * reactance_per_thousand_feet) B[i][i] += -B[i][node] + det = np.linalg.det(B) + + return B, det - return B -def calculate_power_flow(graph, size, voltage_phase, reactance_per_thousand_feet = 0.3): +def calculate_power_flow(graph, voltage_phase, reactance_per_thousand_feet=0.3): set_edge_attributes(graph, 'flow', 0) absolute_flow = 0 + total_flow = 0 + size = len(graph) for i in range(size): neighbors = nx.all_neighbors(graph, i) - for node in neighbors: - if node > i: - graph.edge[i][node]['flow'] = (voltage_phase[node] - voltage_phase[i])/(graph.edge[i][node]['weight'] * reactance_per_thousand_feet) + graph.edge[i][node]['flow'] = (voltage_phase[node] - voltage_phase[i]) / \ + (graph.edge[i][node]['weight'] * reactance_per_thousand_feet) absolute_flow += abs(graph.edge[i][node]['flow']) + total_flow += graph.edge[i][node]['flow'] return graph, absolute_flow + def power_flow_visualization(graph, size, flow_str): return size diff --git a/py/simulated_annealing_algo.py b/py/simulated_annealing_algo.py index 5b0e5df..c781e2f 100644 --- a/py/simulated_annealing_algo.py +++ b/py/simulated_annealing_algo.py @@ -2,20 +2,17 @@ Simulated Annealing """ -import sys -path = "/Users/Thomas/Documents/BlocPower/Microgrids/Simulation/py/" -sys.path.append(path) - import numpy as np import networkx as nx import matplotlib.pyplot as plt - +from power_flow import power_flow, find_b_matrix # TODO: Power flow calculation module to include as the hamiltonian for the simulation # TODO: Gather information about existing simulated annealing algo in python # TODO: Remove PCC definition after each neighbor generated? The PCC should be taken # TODO: into account in the power flow calculation + class MyException(Exception): pass @@ -264,46 +261,101 @@ def neighbor_generator(graph, weighted_matrix, weight_matrix0, flattened_weight_ i*n and i*(n+1), the j indices varying in between from 0 to n-1. This list is used if the graph is connected. distance_matrix_: ndarray Attribute matrix for each edge i,j. - pcc_str: str - Name of the attribute in which to store the PCC (Point of Common Coupling) position. - pos_str: str - Name of the attribute in which the position of each node is stored. """ remove_edge(graph, weight_matrix0) add_edge(graph, weighted_matrix, flattened_weight_matrix, distance_matrix_) +def annealing_algorithm(G, X_coord, Y_coord, Temp_init=1, n_sim=1000, temp_factor=0.95, temp_threshold=0.01, + pcc_node_str="pcc_node"): + """ + Annealing Algorithm: + A probabilistic technique for approximating the global optimum of a given function. + ---------- + G: graph + X_coord, Y_coord: coordinates + Temp: list of temperature + iter_: number of iterations + Returns + ------- + G_new: graph + """ + Distance_Matrix_tmp = distance_matrix(X_coord, Y_coord) + weight_matrix_tmp, weight_matrix0_tmp, flattened_weight_matrix_tmp = weight_matrix(X_coord, Y_coord) + add_random_edge_budget(G, flattened_weight_matrix_tmp, Distance_Matrix_tmp, budget=2) + + pcc_node = G.graph[pcc_node_str] + i = 0 + best = G + pf_best = power_flow(G, pcc_node)[0] + temp = Temp_init + nb_error = 0 + nb_sim_tot = 0 + while temp > temp_threshold: + n = 0 + while n < n_sim: + nb_sim_tot += 1 + G0 = G + neighbor_generator(G0, weight_matrix_tmp, weight_matrix0_tmp, flattened_weight_matrix_tmp, + Distance_Matrix_tmp) + + try: + pf_new = power_flow(G0, pcc_node)[0] + pf_old = power_flow(G, pcc_node)[0] + if pf_new < pf_old: + G = G0 + if pf_new < pf_best: + best = G0 + pf_best = pf_new + else: + p = np.exp(- (pf_new - pf_old) / pf_old / temp) + if np.random.rand() < p: + G = G0 + except np.linalg.LinAlgError: + print("""Error when solving the linear system. The determinant might be 0 and the matrix singular. + Therefore the neighbor is not selected and the next neighbor is tested""") + nb_error += 1 + pass + + n += 1 + i += 1 + temp *= temp_factor ** (i + 1) + print("Proportion of networks for which no solution has been found: %f" % (nb_error/nb_sim_tot)) + + return best, pf_best + + if __name__ == "__main__": import time + from centrality_pcc import set_highest_betweenness_node + from networkx.classes.function import get_node_attributes ######## Neighbor generator testing ####### - size = 500 + size = 20 np.random.seed(2959) X_coord = np.random.rand(size) Y_coord = np.random.rand(size) # Matrix construction Distance_Matrix = distance_matrix(X_coord, Y_coord) - weight_matrix, weight_matrix0, flattened_weight_matrix = weight_matrix(X_coord, Y_coord) - # Graph construction + graph = nx.Graph() G = nx.Graph() for i in range(size): - G.add_node(i, pos=(X_coord[i], Y_coord[i])) + gen = round(np.random.normal(3, .5), 5) + load = round(np.random.normal(3, .5), 5) + G.add_node(i, pos=(X_coord[i], Y_coord[i]), load=load, gen=gen) pos = nx.get_node_attributes(G, 'pos') for i in range(size): for j in range(i + 1, size): G.add_edge(i, j, weight=Distance_Matrix[i][j]) - T = nx.minimum_spanning_tree(G) - # Add 2 edge to the original Graph - add_random_edge_budget(T, flattened_weight_matrix, Distance_Matrix) - #add_random_edge(T) + T = nx.minimum_spanning_tree(G) + set_highest_betweenness_node(T, "PCC_pos", "PCC_node") # Looped generation, removing then adding an edge, and accepting the new state each time - graph = T.copy() - start_time = time.time() - for i in range(1000): - neighbor_generator(graph, weight_matrix, weight_matrix0, flattened_weight_matrix, Distance_Matrix) + T_new = T.copy() + start_time = time.clock() + b, pf_b = annealing_algorithm(T_new, X_coord, Y_coord, n_sim=1000) - print("--- %s seconds ---" % (time.time() - start_time), len(graph.edges())) + print("--- %f minutes ---" % ((time.clock() - start_time)/60)) -- GitLab