diff --git a/py/__init__.py b/py/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 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 Binary files /dev/null and b/py/__pycache__/power_flow.cpython-35.pyc differ diff --git a/py/percolation_threshold_simulation.py b/py/percolation_threshold_simulation.py index cb9b0c2b68ed752930a74e240b1b9fcb9f679ff0..2ef8492f5f444a28300ba9175231c3e8d6b3e36a 100644 --- a/py/percolation_threshold_simulation.py +++ b/py/percolation_threshold_simulation.py @@ -116,6 +116,7 @@ def step_montecarlo(graph, q, n=2, flow_str="flow"): remove_edges_percolation(G, q, flow_str=flow_str) return size_n_largest_component(G, n) + def simulate_montecarlo(n_sim, graph, q_l, n=2, flow_str="flow", print_time=True): """ Returns the matrix of all Monte Carlo simuations, each column corresponding to one @@ -126,7 +127,7 @@ def simulate_montecarlo(n_sim, graph, q_l, n=2, flow_str="flow", print_time=True Parameters ---------- n_sim: int - Number of Monte Carlo simulation + Number of Monte Carlo simulation graph: Networkx.Graph() object Graph to which an edge is removed. q_l: list @@ -146,15 +147,15 @@ def simulate_montecarlo(n_sim, graph, q_l, n=2, flow_str="flow", print_time=True Returns --------- size_matrix: ndarray - Matrix containing all the Monte Carlo simulation (1 per line) for every q value - in q_l (columns) + Matrix containing all the Monte Carlo simulation (1 per line) for every q value + in q_l (columns) mean_flat_matrix: 1darray - 1 * N matrix containing the average value given by the Monte Carlo simulation for - each value of q given by q_l + 1 * N matrix containing the average value given by the Monte Carlo simulation for + each value of q given by q_l max_q_l: list - List of all q (given through q_l) that maximize the size of the nth component - (see n argument). Usually only one value of q is returned. If only one value - of q is passed, then q is returned. + List of all q (given through q_l) that maximize the size of the nth component + (see n argument). Usually only one value of q is returned. If only one value + of q is passed, then q is returned. """ start = time.clock() size_matrix = [[step_montecarlo(graph, q, n=n, flow_str=flow_str) for q in q_l] for j in range(n_sim)] @@ -167,4 +168,4 @@ def simulate_montecarlo(n_sim, graph, q_l, n=2, flow_str="flow", print_time=True max_q_l = [q_l[i] for i in range(len(q_l)) if mean_flat_matrix[i] == max(mean_flat_matrix)] return size_matrix, mean_flat_matrix, max_q_l else: - return size_matrix, mean_flat_matrix, q_l[0] \ No newline at end of file + return size_matrix, mean_flat_matrix, q_l[0] \ No newline at end of file diff --git a/py/power_flow.py b/py/power_flow.py index b761b20c6d14671e064bcffb829cdff9072e7387..b1582251cf7c75ee3733b23800e5faf075ea1406 100644 --- a/py/power_flow.py +++ b/py/power_flow.py @@ -1,9 +1,9 @@ import networkx as nx -from networkx.classes.function import get_edge_attributes, get_node_attributes, set_edge_attributes, set_node_attributes +from networkx.classes.function import get_node_attributes, set_edge_attributes import numpy as np -def power_flow(graph, size, PCC_node, load_str = 'load', gen_str = 'gen', distance_str = 'weight'): +def power_flow(graph, pcc_node, load_str='load', gen_str='gen', distance_str='weight'): """ Find the power flow for each edge on a graph. Return the absolute flow (float) and graph (Networkx.Graph() object) with the following attribute: @@ -15,7 +15,7 @@ def power_flow(graph, size, PCC_node, load_str = 'load', gen_str = 'gen', distan Graph for which to calculate DC power flow. size: int Number of nodes in the graph - PCC_node: str + pcc_node: str Name of the PCC node. load_str: str (optional) Name of the attribute that stores the power load of each node of the graph. @@ -24,71 +24,53 @@ def power_flow(graph, size, PCC_node, load_str = 'load', gen_str = 'gen', distan distance_str: str (optional) Name of the attribute that stores the length of each edge of the graph. """ - PDC = np.zeros((size, 1)) + size = len(graph) + + pdc = np.zeros((size, 1)) load = get_node_attributes(graph, load_str) for node in load: - PDC[node] = graph.node[node][gen_str] - graph.node[node][load_str] - - grid_compensation = find_grid_compensation(graph, load_str, gen_str) + pdc[node] = graph.node[node][gen_str] - graph.node[node][load_str] - if grid_compensation > 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 5b0e5df7f992a40acba97580b6ec3ce3ae049bb2..c781e2f01e4b83cd3ae462696248465df2df573a 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))