diff --git a/.idea/workspace.xml b/.idea/workspace.xml new file mode 100644 index 0000000..6ef02cb --- /dev/null +++ b/.idea/workspace.xml @@ -0,0 +1,85 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 1745565553589 + + + + + + + + file://$PROJECT_DIR$/main.py + 72 + + + + + + + + + + + \ No newline at end of file diff --git a/__pycache__/create_graphs.cpython-310.pyc b/__pycache__/create_graphs.cpython-310.pyc new file mode 100644 index 0000000..31cb710 Binary files /dev/null and b/__pycache__/create_graphs.cpython-310.pyc differ diff --git a/__pycache__/cutcalculator.cpython-310.pyc b/__pycache__/cutcalculator.cpython-310.pyc new file mode 100644 index 0000000..f12d9d1 Binary files /dev/null and b/__pycache__/cutcalculator.cpython-310.pyc differ diff --git a/__pycache__/loadflowcalculator.cpython-310.pyc b/__pycache__/loadflowcalculator.cpython-310.pyc new file mode 100644 index 0000000..13eb2c2 Binary files /dev/null and b/__pycache__/loadflowcalculator.cpython-310.pyc differ diff --git a/__pycache__/radialgraph.cpython-310.pyc b/__pycache__/radialgraph.cpython-310.pyc new file mode 100644 index 0000000..9602748 Binary files /dev/null and b/__pycache__/radialgraph.cpython-310.pyc differ diff --git a/create_graphs.py b/create_graphs.py index 0d7db40..9826c95 100644 --- a/create_graphs.py +++ b/create_graphs.py @@ -93,7 +93,7 @@ def create_linear_graph() -> RadialGraph: return graph -def create_IEEE33bus() -> RadialGraph: +def create_IEEE33bus_load() -> RadialGraph: # Create a radial graph graph = RadialGraph() # Create nodes @@ -202,4 +202,115 @@ def create_IEEE33bus() -> RadialGraph: #graph.add_edge(node8, node14, 2.0000, 2.0000) #graph.add_edge(node11, node21, 2.0000, 2.0000) #graph.add_edge(node17, node32, 0.5000, 0.5000) + return graph + +def create_IEEE33bus_gen() -> RadialGraph: + # Create a radial graph + graph = RadialGraph() + # Create nodes + root = Root("Root") + node1 = Node("Node 1", -100, -60) # p, q in kW, 3 phases (balanced load is assumed) + node2 = Node("Node 2", -90, -40) + node3 = Node("Node 3", -120, -80) + node4 = Node("Node 4", -60, -30) # p, q in kW, 3 phases (balanced load is assumed) + node5 = Node("Node 5", -60, -20) + node6 = Node("Node 6", -200, -100) + node7 = Node("Node 7", -200, -100) # p, q in kW, 3 phases (balanced load is assumed) + node8 = Node("Node 8", -100, -50) + node9 = Node("Node 9", -60, -20) + node10 = Node("Node 10", -45, -30) # p, q in kW, 3 phases (balanced load is assumed) + node11 = Node("Node 11", -60, -35) + node12 = Node("Node 12", -60, -35) + node13 = Node("Node 13", -120, -80) # p, q in kW, 3 phases (balanced load is assumed) + node14 = Node("Node 14", -70, -30) + node15 = Node("Node 15", -60, -20) + node16 = Node("Node 16", -90, -30) # p, q in kW, 3 phases (balanced load is assumed) + node17 = Node("Node 17", -100, -50) + node18 = Node("Node 18", -90, -40) + node19 = Node("Node 19", -90, -40) # p, q in kW, 3 phases (balanced load is assumed) + node20 = Node("Node 20", -90, -40) + node21 = Node("Node 21", -90, -40) + node22 = Node("Node 22", -90, -50) # p, q in kW, 3 phases (balanced load is assumed) + node23 = Node("Node 23", -420, -200) + node24 = Node("Node 24", -420, -200) + node25 = Node("Node 25", -60, -25) + node26 = Node("Node 26", -60, -25) # p, q in kW, 3 phases (balanced load is assumed) + node27 = Node("Node 27", -60, -20) + node28 = Node("Node 28", -120, -70) + node29 = Node("Node 29", -200, -600) # p, q in kW, 3 phases (balanced load is assumed) + node30 = Node("Node 30", -150, -70) + node31 = Node("Node 31", -210, -100) + node32 = Node("Node 32", -60, -40) # p, q in kW, 3 phases (balanced load is assumed) + # Add nodes to the graph + graph.add_node(root) + graph.add_node(node1) + graph.add_node(node2) + graph.add_node(node3) + graph.add_node(node4) + graph.add_node(node5) + graph.add_node(node6) + graph.add_node(node7) + graph.add_node(node8) + graph.add_node(node9) + graph.add_node(node10) + graph.add_node(node11) + graph.add_node(node12) + graph.add_node(node13) + graph.add_node(node14) + graph.add_node(node15) + graph.add_node(node16) + graph.add_node(node17) + graph.add_node(node18) + graph.add_node(node19) + graph.add_node(node20) + graph.add_node(node21) + graph.add_node(node22) + graph.add_node(node23) + graph.add_node(node24) + graph.add_node(node25) + graph.add_node(node26) + graph.add_node(node27) + graph.add_node(node28) + graph.add_node(node29) + graph.add_node(node30) + graph.add_node(node31) + graph.add_node(node32) + # Add edges to the graph + graph.add_edge(root, node1, 0.0922, 0.0470) + graph.add_edge(node1, node2, 0.4930, 0.2511) + graph.add_edge(node2, node3, 0.3660, 0.1864) + graph.add_edge(node3, node4, 0.3811, 0.1941) + graph.add_edge(node4, node5, 0.8190, 0.7070) + graph.add_edge(node5, node6, 0.1872, 0.6188) + graph.add_edge(node6, node7, 0.7114, 0.2351) + graph.add_edge(node7, node8, 1.0300, 0.7400) + graph.add_edge(node8, node9, 1.0440, 0.7400) + graph.add_edge(node9, node10, 0.1966, 0.0650) + graph.add_edge(node10, node11, 0.3744, 0.1238) + graph.add_edge(node11, node12, 1.4680, 1.1550) + graph.add_edge(node12, node13, 0.5416, 0.7129) + graph.add_edge(node13, node14, 0.5910, 0.5260) + graph.add_edge(node14, node15, 0.7463, 0.5450) + graph.add_edge(node15, node16, 1.2890, 1.7210) + graph.add_edge(node16, node17, 0.7320, 0.5740) + graph.add_edge(node1, node18, 0.1640, 0.1565) + graph.add_edge(node18, node19, 1.5042, 1.3554) + graph.add_edge(node19, node20, 0.4095, 0.4784) + graph.add_edge(node20, node21, 0.7089, 0.9373) + graph.add_edge(node2, node22, 0.4512, 0.3083) + graph.add_edge(node22, node23, 0.8980, 0.7091) + #graph.add_edge(node24, node28, 0.8960, 0.7011) + graph.add_edge(node23, node24, 0.8960, 0.7011) + graph.add_edge(node5, node25, 0.2030, 0.1034) + graph.add_edge(node25, node26, 0.2842, 0.1447) + graph.add_edge(node26, node27, 1.0590, 0.9337) + graph.add_edge(node27, node28, 0.8042, 0.7006) + graph.add_edge(node28, node29, 0.5075, 0.2585) + graph.add_edge(node29, node30, 0.9744, 0.9630) + graph.add_edge(node30, node31, 0.3105, 0.3619) + graph.add_edge(node31, node32, 0.3410, 0.5302) + #graph.add_edge(node20, node7, 2.0000, 2.0000) + #graph.add_edge(node8, node14, 2.0000, 2.0000) + #graph.add_edge(node11, node21, 2.0000, 2.0000) + #graph.add_edge(node17, node32, 0.5000, 0.5000) return graph \ No newline at end of file diff --git a/cutcalculator.py b/cutcalculator.py index 8b199d2..f132978 100644 --- a/cutcalculator.py +++ b/cutcalculator.py @@ -39,14 +39,21 @@ def get_cut_slopes_intercepts(self, v_properties: dict, include_losses: bool, lo loss_parameter = 1: losses are calculated based on the branch flows at the end of the iteration loss_parameter = 0.5: branch flows are the weighted sum of the flows at the beginning/end of the iteration """ - # First calculate the voltage term - v_term_intercept = self._calculate_voltage_intercept_term(v_properties) - # Calculate the bucket filling information for p and q total_p_load = sum([n.p_max for n in self.radial_graph.nodes]) total_q_load = sum([n.q_max for n in self.radial_graph.nodes]) active_info = self.radial_graph.get_active_bucket_filling_info(total_p_load) reactive_info = self.radial_graph.get_reactive_bucket_filling_info(total_q_load) + + # First calculate the voltage term + if total_p_load > 0.0 and total_q_load > 0.0: + loads = True + elif total_p_load < 0.0 and total_p_load < 0.0: + loads = False + else: + raise ValueError(f'Total p and q should either be positive or negative') + + v_term_intercept = self._calculate_voltage_intercept_term(v_properties, loads) print("Bucket fill information:") print("Active power: ", active_info) @@ -88,8 +95,11 @@ def get_cut_slopes_intercepts(self, v_properties: dict, include_losses: bool, lo return slopes, intercepts @staticmethod - def _calculate_voltage_intercept_term(v_properties: dict) -> float: - return (v_properties['v_base'] ** 2 - v_properties['v_min'] ** 2) / 2 + def _calculate_voltage_intercept_term(v_properties: dict, loads: bool) -> float: + if loads: + return (v_properties['v_base'] ** 2 - v_properties['v_min'] ** 2) / 2 + else: + return (v_properties['v_base'] ** 2 - v_properties['v_max'] ** 2) / 2 def _calculate_active_power_terms(self, active_bucket_filling_information: dict) -> Tuple[List[float]]: print("Calculate P terms") diff --git a/main.py b/main.py index aeb61fc..d99ff3e 100644 --- a/main.py +++ b/main.py @@ -1,14 +1,25 @@ -from create_graphs import create_linear_graph, create_simple_branch_graph, create_graph, create_IEEE33bus +from create_graphs import create_linear_graph, create_simple_branch_graph, create_graph, create_IEEE33bus_gen, create_IEEE33bus_load from loadflowcalculator import LoadFlowCalculator from cutcalculator import CutCalculator import numpy as np import warnings import matplotlib.pyplot as plt +def get_P_from_lf_results(lf_results: dict) -> float: + return max(list(lf_results['P_edges']), key=abs) + +def get_Q_from_lf_results(lf_results: dict) -> float: + return max(list(lf_results['Q_edges']), key=abs) if __name__ == "__main__": # Create a radial graph - graph = create_IEEE33bus() + + all_generation = True # Either all loads or all generation! + # graph = create_IEEE33bus_load() + if all_generation: + graph = create_IEEE33bus_gen() + else: + graph = create_IEEE33bus_load() # Create loadflow calculator load_flow_calculator = LoadFlowCalculator(graph) @@ -19,7 +30,8 @@ # Set voltage parameters v_base = 12.66 * 1.0e3 / np.sqrt(3) # 12.66 * 1.0e3 v_min = 0.95 * v_base - v_properties = {'v_base': v_base, 'v_min': v_min} + v_max = 1.05 * v_base + v_properties = {'v_base': v_base, 'v_min': v_min, 'v_max': v_max} ''' Calculate the sampled FOR @@ -30,8 +42,11 @@ max_active_power = sum([n.p_max for n in graph.nodes]) max_reactive_power = sum([n.q_max for n in graph.nodes]) - p_totals = np.linspace(0, max_active_power, n_points) - q_totals = np.linspace(0, max_reactive_power, n_points) + p_totals = np.linspace(min(max_active_power, 0.0), max(0.0, max_active_power), n_points) + q_totals = np.linspace(min(max_reactive_power, 0.0), max(0.0, max_reactive_power), n_points) + + # p_totals = np.linspace(0.0, max(0.0, max_active_power), n_points) + # q_totals = np.linspace(min(max_reactive_power, 0.0), max(0.0, max_active_power), n_points) # ''' # Allocate result lists @@ -42,8 +57,8 @@ Qs_transformer_no_loss = [] count = 0 # keep track of how many load flow calculations we have done - for i, p_total in enumerate(p_totals): - for q_total in q_totals: + for i, p_total in enumerate(p_totals[::-1]): + for q_total in q_totals[::-1]: # Calculate individual loads by means of bucket filling loads_p_dict = graph.apply_active_bucket_filling(p_total) loads_p = [load for node_name, load in loads_p_dict.items() if node_name != 'Root'] @@ -60,15 +75,18 @@ # Check if voltage constraint is violated and safe P and Q through the interconnection # Pandapower loadflow: v_nodes = lf_pp_results['v_nodes'] - if min(v_nodes) > v_min: - Ps_transformer.append(max(lf_pp_results['P_edges'])) - Qs_transformer.append(max(lf_pp_results['Q_edges'])) + if min(v_nodes) > v_min and max(v_nodes) < v_max: + Ps_transformer.append(get_P_from_lf_results(lf_pp_results)) + Qs_transformer.append(get_Q_from_lf_results(lf_pp_results)) # Lossless loadflow: v_nodes = lf_lossless_results['v_nodes'] - if min(v_nodes) > v_min: - Ps_transformer_no_loss.append(max(lf_lossless_results['P_edges'])) - Qs_transformer_no_loss.append(max(lf_lossless_results['Q_edges'])) + if min(v_nodes) > v_min and max(v_nodes) < v_max: + print(min(v_nodes), max(v_nodes)) + print(v_min, v_max) + Ps_transformer_no_loss.append(get_P_from_lf_results(lf_lossless_results)) + Qs_transformer_no_loss.append(get_Q_from_lf_results(lf_lossless_results)) + except: # If the pandapower load flow did not converge we skip this point warnings.warn("Pandapower loadflow did not converge") @@ -94,16 +112,12 @@ slopes, intercepts = cut_calculator.get_cut_slopes_intercepts(v_properties, include_losses=False) slopes_losses, intercepts_losses = cut_calculator.get_cut_slopes_intercepts(v_properties, include_losses=True, - loss_parameter=0.0) + loss_parameter=1.0 * all_generation) # ''' ''' Plotting ''' - # Samples - plt.scatter(Ps_transformer, Qs_transformer, color='b', label='Full Distflow') - plt.scatter(Ps_transformer_no_loss, Qs_transformer_no_loss, color='r', label='Linear Distflow') - # Cuts # ''' cut_ids = slopes.keys() @@ -121,10 +135,16 @@ q_values = intercept + slope * 1000 * p_totals # p_totals is in kW q_values_losses = intercept_loss + slope * 1000 * p_totals - # plt.plot(p_totals, q_values/1000.0, 'r', label=str(cut_id)) - plt.plot(p_totals, q_values_losses / 1000.0, 'b', label=str(cut_id)) + # plt.plot(p_totals, q_values/1000.0, 'r', label=str(cut_id), alpha=0.5) + plt.plot(p_totals, q_values_losses / 1000.0, 'b', label=str(cut_id), alpha=0.1) + + # Samples + plt.scatter(Ps_transformer, Qs_transformer, color='b', label='Full Distflow') + plt.scatter(Ps_transformer_no_loss, Qs_transformer_no_loss, color='r', label='Linear Distflow') + + - # plt.legend() + plt.legend() # ''' - plt.ylim(0, 1.1 * max_reactive_power) + plt.ylim(min(0, 1.1 * max_reactive_power), max(0, 1.1 * max_reactive_power)) plt.show() diff --git a/radialgraph.py b/radialgraph.py index 0ff3c2f..e0da73f 100644 --- a/radialgraph.py +++ b/radialgraph.py @@ -8,8 +8,8 @@ class Node: def __init__(self, name: str, p_max: float, q_max: float): # Apply checks - assert p_max >= 0, "Please don't use only positive values (loads) for p_max" - assert q_max >= 0, "Please don't use only positive values (loads) for q_max" + # assert p_max >= 0, "Please don't use only positive values (loads) for p_max" + # assert q_max >= 0, "Please don't use only positive values (loads) for q_max" self.name = name self.p_max = p_max @@ -188,7 +188,7 @@ def _get_bucket_filling_info(self, load_type: str, edge_prop_type: str, load_to_ max_load_type = load_type + '_max' # either p_max or q_max # Do some input tests - assert load_to_fill <= sum([getattr(n, max_load_type) for n in self.nodes]), "provided load_to_fill is too big" + assert load_to_fill <= abs(sum([getattr(n, max_load_type) for n in self.nodes])), "provided load_to_fill is too big" # Allocate output: how much load to put at every node loads = {n.name: 0.0 for n in self.nodes} @@ -208,7 +208,7 @@ def _get_bucket_filling_info(self, load_type: str, edge_prop_type: str, load_to_ iteration_information['loads_end'] = loads.copy() information_dict['iteration_' + str(iteration)] = iteration_information - while round(load_to_fill, 5) > 0.0: # round is for preventing small residual load due to floating point errors + while abs(round(load_to_fill, 5)) > 0.0: # round is for preventing small residual load due to floating point errors # print(f"load_to_fill: {load_to_fill}") # print(f"nodes_to_fill: {[n.name for n in nodes_to_fill]}") # Set up dict to store outputs @@ -221,13 +221,14 @@ def _get_bucket_filling_info(self, load_type: str, edge_prop_type: str, load_to_ fill_factors = self._calculate_fill_factors(edge_prop_paths) # Find the loads that are first the be full given the fill factors - total_load_full = [remaining_bucket_sizes[n.name] / fill_factors[i] for i, n in enumerate(nodes_to_fill)] + total_load_full = [abs(remaining_bucket_sizes[n.name]) / fill_factors[i] for i, n in enumerate(nodes_to_fill)] # print(f'Total load full: {total_load_full}') min_total_load_full = min(total_load_full) nodes_max_fill = [n for i, n in enumerate(nodes_to_fill) if total_load_full[i] == min_total_load_full] # Constrain the total amount of filling by the amount that is left to fill - total_fill = min(min_total_load_full, load_to_fill) + # take care of sign. for loads, we can remove the abs and the sign + total_fill = np.sign(load_to_fill) * min(abs(min_total_load_full), abs(load_to_fill)) fills = [factor * total_fill for factor in fill_factors] # print(f'Total fill: {total_fill}') # print(f'fills: {fills}') @@ -246,7 +247,7 @@ def _get_bucket_filling_info(self, load_type: str, edge_prop_type: str, load_to_ # print(f'New remaining_bucket_sizes: {remaining_bucket_sizes}') # Check whether we have leftover load to fill. If so, update nodes_to_fill - if load_to_fill > total_fill: + if abs(load_to_fill) > abs(total_fill): # print("NOT FULL THIS ROUND", load_to_fill, total_fill) # print(f'Old nodes to fill: {[n.name for n in nodes_to_fill]}') # print(f'Nodes full: {[n.name for n in nodes_max_fill]}')