added all generation case

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

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")

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()